5
$\begingroup$

I am trying to calculate the UCCSD energy and UCCSD(T) correction for the S atom using PySCF. Using sto-3g basis set, I found that the UCCSD calculation would be completed successfully, but the calculation of UCCSD(T) correction failed with the error message reproduced below.

sto-3g input file

#!/usr/bin/python

import pyscf    as pyscf
import numpy    as np
import sys      as sys
from pyscf      import gto, scf, ao2mo, fci, ci, cc
from pyscf.lib  import logger

def calculate_energy() -> float:
 
    neutralMolecule     = pyscf.M(atom = [['S', (0,0,0)]], basis = 'sto-3g', verbose = 9, charge = 0, spin = 2, symmetry = True)
    neutralMoleculeHF   = scf.UHF(neutralMolecule).run()

    logObject               =   logger.new_logger(neutralMoleculeHF)                # Creates a new logger object
    orbitals,_,stability,_  =   neutralMoleculeHF.stability(return_status = True)   # Returns four variables: Internal MO, External MO, Internal Stability, External Stability
    optimizationCycles      =   0
    while (not stability and optimizationCycles < 200): 
        logObject.note("Currently at {}-th step optimizing the wavefunction".format(optimizationCycles))
    densityMatrix           =   neutralMoleculeHF.make_rdm1(orbitals, neutralMoleculeHF.mo_occ)
    neutralMoleculeHF       =   neutralMoleculeHF.run(densityMatrix) 
    orbitals,_,stability,_  =   neutralMoleculeHF.stability(return_status = True)
    optimizationCycles      =   optimizationCycles + 1
    if not stability:
        print("Unable to achieve stable wavefunction after {} steps. Exiting".format(optimizationCycles))
        sys.exit(-1)
    
    # Now start Post-HF calculations
    neutralMoleculeCC           =   cc.UCCSD(neutralMoleculeHF).set(conv_tol = 1e-7, frozen = 5, direct = True)
    neutralMoleculeCC.run()
    neutralMoleculeCC_results   =   neutralMoleculeCC.kernel()
    neutralMoleculeCC_CCSDT     =   neutralMoleculeCC.ccsd_t()   
    return neutralMoleculeCC.e_tot + neutralMoleculeCC_CCSDT

energy  =   calculate_energy()

sto-3g stack trace

DIIS for step 0
cycle = 1  E_corr(CCSD) = 0  dE = 0  norm(t1,t2) = 0
    CPU time for CCSD iter      0.07 sec, wall time      0.02 sec
    CPU time for CCSD      0.07 sec, wall time      0.02 sec
UCCSD converged
E(UCCSD) = -393.130219398113  E_corr = 0
    CPU time for transpose 0:2      0.06 sec, wall time      0.02 sec
    CPU time for transpose 0:2      0.01 sec, wall time      0.00 sec
    CPU time for UCCSD(T) sort_eri      0.08 sec, wall time      0.03 sec
max_memory 3874 MB (125 MB in use)
    CPU time for contract_aaa      0.00 sec, wall time      0.00 sec
max_memory 3874 MB (125 MB in use)
    CPU time for contract 0:2,0:2      0.04 sec, wall time      0.01 sec
    CPU time for contract_bbb      0.04 sec, wall time      0.01 sec
Traceback (most recent call last):
  File "/home/hemanth/Desktop/hpqc/cs/cs_pes_scan.py", line 34, in <module>
    energy  =   calculate_energy(scanDistance=0)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/hemanth/Desktop/hpqc/cs/cs_pes_scan.py", line 31, in calculate_energy
    neutralMoleculeCC_CCSDT     =   neutralMoleculeCC.ccsd_t()   
                                    ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/hemanth/miniconda3/envs/hpqc/lib/python3.11/site-packages/pyscf/cc/uccsd.py", line 635, in ccsd_t
    return uccsd_t.kernel(self, eris, t1, t2, self.verbose)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/hemanth/miniconda3/envs/hpqc/lib/python3.11/site-packages/pyscf/cc/uccsd_t.py", line 139, in kernel
    for a0, a1 in lib.prange(0, nvirb, int(bufsize/nvira+1)):
                                           ~~~~~~~^~~~~~
ZeroDivisionError: division by zero

However, if I change the basis set to 6-31g, the UCCSD and UCCSD(T) calculations would be completed without any issues.

6-31g input file

#!/usr/bin/python

import pyscf    as pyscf
import numpy    as np
import sys      as sys
from pyscf      import gto, scf, ao2mo, fci, ci, cc
from pyscf.lib  import logger

def calculate_energy() -> float:
 
    neutralMolecule     = pyscf.M(atom = [['S', (0,0,0)]], basis = '6-31g', verbose = 9, charge = 0, spin = 2, symmetry = True)
    neutralMoleculeHF   = scf.UHF(neutralMolecule).run()

    logObject               =   logger.new_logger(neutralMoleculeHF)                # Creates a new logger object
    orbitals,_,stability,_  =   neutralMoleculeHF.stability(return_status = True)   # Returns four variables: Internal MO, External MO, Internal Stability, External Stability
    optimizationCycles      =   0
    while (not stability and optimizationCycles < 200): 
        logObject.note("Currently at {}-th step optimizing the wavefunction".format(optimizationCycles))
    densityMatrix           =   neutralMoleculeHF.make_rdm1(orbitals, neutralMoleculeHF.mo_occ)
    neutralMoleculeHF       =   neutralMoleculeHF.run(densityMatrix) 
    orbitals,_,stability,_  =   neutralMoleculeHF.stability(return_status = True)
    optimizationCycles      =   optimizationCycles + 1
    if not stability:
        print("Unable to achieve stable wavefunction after {} steps. Exiting".format(optimizationCycles))
        sys.exit(-1)
    
    # Now start Post-HF calculations
    neutralMoleculeCC           =   cc.UCCSD(neutralMoleculeHF).set(conv_tol = 1e-7, frozen = 5, direct = True)
    neutralMoleculeCC.run()
    neutralMoleculeCC_results   =   neutralMoleculeCC.kernel()
    neutralMoleculeCC_CCSDT     =   neutralMoleculeCC.ccsd_t()   
    return neutralMoleculeCC.e_tot + neutralMoleculeCC_CCSDT

energy  =   calculate_energy()

6-31g output file

DIIS for step 7
cycle = 8  E_corr(CCSD) = -0.0229167339083696  dE = 2.82885328e-08  norm(t1,t2) = 7.6597e-07
    CPU time for CCSD iter      0.26 sec, wall time      0.06 sec
    CPU time for CCSD      2.52 sec, wall time      0.65 sec
UCCSD converged
E(UCCSD) = -397.4943311312284  E_corr = -0.02291673390836959
    CPU time for transpose 0:4      0.01 sec, wall time      0.00 sec
    CPU time for transpose 0:6      0.03 sec, wall time      0.01 sec
    CPU time for transpose 0:4      0.03 sec, wall time      0.01 sec
    CPU time for transpose 0:6      0.01 sec, wall time      0.00 sec
    CPU time for UCCSD(T) sort_eri      0.10 sec, wall time      0.03 sec
max_memory 3874 MB (125 MB in use)
    CPU time for contract 0:4,0:4      0.02 sec, wall time      0.00 sec
    CPU time for contract_aaa      0.02 sec, wall time      0.00 sec
max_memory 3874 MB (125 MB in use)
    CPU time for contract 0:6,0:6      0.01 sec, wall time      0.00 sec
    CPU time for contract_bbb      0.01 sec, wall time      0.00 sec
    CPU time for contract 0:6,0:4      0.01 sec, wall time      0.00 sec
    CPU time for contract_baa      0.01 sec, wall time      0.00 sec
    CPU time for contract 0:4,0:6      0.01 sec, wall time      0.00 sec
    CPU time for contract_abb      0.01 sec, wall time      0.00 sec
    CPU time for UCCSD(T)      0.15 sec, wall time      0.04 sec
UCCSD(T) correction = -0.000187450473417524

Any pointers on why this is happening would be much appreciated.

Improve this question
$\endgroup$
1

1 Answer 1

Reset to default
6
$\begingroup$

The stack trace shows that the program tries to divide a number by the number of alpha virtual orbitals (nvira), which is zero in the case of STO-3G triplet sulfur atom. This caused a ZeroDivisionError, leading to a program abort. The number of alpha virtual orbitals can only be zero when the basis set is a minimal basis set.

While one may argue that a decent program should implement a branch that deals with this case (since the CCSD(T)/STO-3G energy is a mathematically valid quantity and has a well-defined value), a CCSD(T)/STO-3G calculation is not physically meaningful, as one can get much better accuracy at much lower cost with a DFT/double zeta basis calculation. Therefore, this ZeroDivisionError does not necessarily count as a bug, because it only happens in calculations whose basis sets are so small that the results are not physically useful.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ +1 but people may be interested in CCSD(T)/STO-3G for benchmarking reasons or for composite methods or other reasons apart from testing. In this case, the OP is learning how to build potential energy curves and eventually we'll use big basis sets with CCSD(T) but I recommended STO-3G to start, in order to get speed and ease. $\endgroup$
    – Nike Dattani
    Jun 12 at 13:39
  • 1
    $\begingroup$ @NikeDattani True. If I were a developer of PySCF and had enough time, I would fix this issue, or at least make the program exit gracefully (with a user-friendly error message) before the zero division happens. $\endgroup$
    – wzkchem5
    Jun 12 at 17:25
  • $\begingroup$ If there are no virtual orbitals, then there is no correlation energy to compute (as is seen in the working case of UCCSD). It would presumably be straightforward to fix (before starting any post-HF calculation, check that there are >0 virtual orbitals), but it is a fairly niche case and arguably the "working" case is worse, as it can be fairly easy to miss that because there are no virtual orbitals, you are just performing a Hartree-Fock calculation. $\endgroup$
    – Tyberius ♦
    Jun 12 at 19:07
  • 1
    $\begingroup$ @Tyberius True, for a calculation with zero virtual orbitals, the program can return zero correlation energy right at the beginning of the post-HF calculation. However, if there are no alpha virtual orbitals but >=3 beta virtual orbitals, the (T) energy is not trivially zero and must be calculated. (If there are 2 beta virtual orbitals, the CCSD correlation energy is non-zero but the (T) energy is zero; if there is 1 beta virtual orbital, the CCSD correlation energy is also zero as long as the HF calculation is fully converged.) $\endgroup$
    – wzkchem5
    Jun 12 at 20:30

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .