CCSD Quadratic Response Function

.github/workflows/CI.yaml

@@ -21,8 +21,8 @@ jobs:
     strategy:
       fail-fast: true
       matrix:
-        os: [macOS-latest, ubuntu-latest]
-        python-version: ['3.11']
+        os: [macOS-latest]
+        python-version: ['3.8','3.11']
 
     steps:
     - uses: actions/checkout@v2
@@ -49,6 +49,7 @@ jobs:
         hash -r
         env
         conda install psi4 python=${{ matrix.python-version }} -c conda-forge
+        conda install pytorch
         ls -l $CONDA
 
     - name: Test Psi4 Python Loading
@@ -67,7 +68,7 @@ jobs:
       run: |
         export KMP_DUPLICATE_LIB_OK=TRUE
         pip install coverage
-        coverage run -m pytest
+        coverage run -m pytest -s pycc/tests/test_038_sim_quadresp.py::test_PNOpp_ccsd_OR
         coverage xml
         ls -l
 

pycc/ccwfn.py

@@ -87,13 +87,14 @@ def __init__(self, scf_wfn, **kwargs):
 
         self.make_t3_density = kwargs.pop('make_t3_density', False)
 
-        valid_local_models = [None, 'PNO', 'PAO','PNO++']
+        valid_local_models = [None, 'PNO', 'PAO','CPNO++','PNO++']
         local = kwargs.pop('local', None)
         # TODO: case-protect this kwarg
         if local not in valid_local_models:
             raise Exception("%s is not an allowed local-CC model." % (local))
         self.local = local
         self.local_cutoff = kwargs.pop('local_cutoff', 1e-5)
+        self.ed_omega = kwargs.pop('omega', 0) 
 
         valid_local_MOs = ['PIPEK_MEZEY', 'BOYS']
         local_MOs = kwargs.pop('local_mos', 'PIPEK_MEZEY')
@@ -151,7 +152,7 @@ def __init__(self, scf_wfn, **kwargs):
         self.H = Hamiltonian(self.ref, self.C, self.C, self.C, self.C)
 
         if local is not None:
-            self.Local = Local(local, self.C, self.nfzc, self.no, self.nv, self.H, self.local_cutoff,self.it2_opt)
+            self.Local = Local(local, self.C, self.nfzc, self.no, self.nv, self.H, self.local_cutoff,self.it2_opt, self.ed_omega)
             if filter is not True:
                 self.Local.trans_integrals(self.o, self.v)
                 self.Local.overlaps(self.Local.QL)

pycc/data/molecules.py

@@ -26,6 +26,16 @@
         symmetry c1
         """
 
+# CO
+co = """
+C  0.000000  0.000000  0.000000
+O  0.000000  0.000000  -2.132000
+noreorient
+nocom
+units au
+symmetry c1
+"""
+
 # Håkon's H2O test case
 h2o_hek="""
         units au
@@ -43,15 +53,24 @@
 symmetry c1
 units bohr
 """
-### Water molecule
 
+### Water molecule
 h2o = """
 O
 H 1 1.1
 H 1 1.1 2 104
 symmetry c1
 """
 
+### Water molecule from Dalton
+h2o_dalton = """
+O 0.000000000000     0.000000000000     0.000000000000
+H 0.000000000000    -0.756689920000     0.585891940000
+H 0.000000000000     0.756689920000     0.585891940000
+units Angstrom
+symmetry c1
+"""
+
 ### Water cluster
 ## Number of water molecules
 ## 2
@@ -269,10 +288,12 @@
 moldict["He"] = he
 moldict["Be"] = be
 moldict["LiH"] = lih
+moldict["CO"] = co
 moldict["H2"] = h2
 moldict["H2O_HEK"] = h2o_hek
 moldict["H2O_Teach"] = h2o_tutorial
 moldict["H2O"] = h2o
+moldict["H2O_D"] = h2o_dalton
 moldict["(H2O)_2"] = h2o_2
 moldict["(H2O)_3"] = h2o_3
 moldict["(H2O)_4"] = h2o_4

pycc/local.py

@@ -49,6 +49,7 @@ class Local(object):
     _build_PAO(): build PAO orbital rotation tensors
     _build_PNO(): build PNO orbital rotation tensors   
     _build_PNOpp(): build PNO++ orbital rotation tensors
+    _build_cPNOpp(): build cPNO++ orbital rotation tensors 
     trans_integrals(): transform Fock matrix and ERI from the MO basis to a local basis 
 
     Notes
@@ -58,7 +59,7 @@ class Local(object):
      to run local MP2, uncomment the necessary lines within the _build_"local" functions which are at the end 
     """
 
-    def __init__(self, local, C, nfzc, no, nv, H, cutoff, it2_opt,
+    def __init__(self, local, C, nfzc, no, nv, H, cutoff, it2_opt, omega, 
             core_cut=5E-2,
             lindep_cut=1E-6,
             e_conv=1e-12,
@@ -72,6 +73,7 @@ def __init__(self, local, C, nfzc, no, nv, H, cutoff, it2_opt,
         self.C = C.to_array()
         self.local = local
         self.it2_opt = it2_opt
+        self.omega = omega
         self.core_cut = core_cut
         self.lindep_cut = lindep_cut
         self.e_conv = e_conv
@@ -86,6 +88,8 @@ def _build(self):
             self._build_PAO()
         elif self.local.upper() in ["PNO++"]:
             self._build_PNOpp()
+        elif self.local.upper() in ["CPNO++"]:
+            self._build_cPNOpp()
         else:
             raise Exception("Not a valid local type!")
 
@@ -326,14 +330,14 @@ def _build_PNO(self):
         # MP2 loop (optional)
         if self.it2_opt:
             self._MP2_loop(t2,self.H.F,self.H.ERI,self.H.L,Dijab)
-         
+
         print("Computing PNOs.  Canonical VMO dim: %d" % (self.nv))
 
         #Constructing the PNO density
         D = self._pairdensity(t2) 
-
+        
         # Now obtain the Q and L
-        Q, L, eps, dim = self.QL_tensors(v,self.local,t2,D)
+        Q, L, eps, dim = self.QL_tensors(v,t2,D,local ='PNO')
 
         self.Q = Q  # transform between canonical VMO and local spaces
         self.L = L  # transform between local and semicanonical local spaces      
@@ -375,17 +379,18 @@ def _build_PNOpp(self):
         Dijab = eps_occ.reshape(-1,1,1,1) + eps_occ.reshape(-1,1,1) - eps_vir.reshape(-1,1) - eps_vir
 
         # initial guess amplitudes
-        t2 = self.H.ERI[o,o,v,v]/Dijab
+        t2 = self.H.ERI[o,o,v,v]/(Dijab + self.omega)
+
 
         # MP2 loop (optional) 
         if self.it2_opt:
-            self._MP2_loop(t2,self.H.F,self.H.ERI,self.H.L,Dijab)
+            self._MP2_loop(t2,self.H.F,self.H.ERI,self.H.L,Dijab + self.omega)
 
         # Construct the perturbed pair density, Eqn. 10  
         D = self._pert_pairdensity(t2)
 
         # Now obtain Q and L 
-        Q, L, eps, dim = self.QL_tensors(v,self.local,t2,D)       
+        Q, L, eps, dim = self.QL_tensors(v,t2,D,local ='PNO++')       
 
         self.Q = Q  # transform between canonical VMO and local spaces
         self.L = L  # transform between local and semicanonical local spaces 
@@ -401,6 +406,66 @@ def _build_PNOpp(self):
                 self.Q[ji] = self.Q[ij]
                 self.L[ji] = self.L[ij]
 
+    def _build_cPNOpp(self):  
+        """
+        Perform MP2 loop in non-canonical MO basis, then construct pair density based on t2 amplitudes
+        then build MP2-level PNO
+
+        Attributes
+        ----------
+        Q: transform between canonical VMO and PNO++ spaces
+        L: transform between LPNO and semicanonical PNO++ spaces
+        dim: dimension of cPNO++ space
+        eps: semicananonical cPNO++ energies
+
+        Notes
+        -----
+        Equations from D'Cunha & Crawford 2021 [10.1021/acs.jctc.0c01086]
+        """ 
+        v = slice(self.no, self.no+self.nv)
+
+        self._build_PNO()
+        Q_PNO = self.Q
+
+        self._build_PNOpp()
+        Q_PNOpp = self.Q 
+
+        self.Q = []  # truncated PNO list
+        self.dim = np.zeros((self.no*self.no), dtype=int)  # dimension of local space for each pair
+        self.L = [] # semicanonical PNO list
+        self.eps = [] # approximated virtual orbital energies
+        T2_local = 0
+        for ij in range(self.no*self.no):
+            i = ij // self.no
+            j = ij % self.no
+
+            Q_comb = np.hstack((Q_PNO[ij], Q_PNOpp[ij]))
+            Q_ortho, trash = np.linalg.qr(Q_comb)
+            self.Q.append(Q_ortho)
+            # Compute semicanonical virtual space
+            F = Q_ortho.T @ self.H.F[v,v] @ Q_ortho  # Fock matrix in local basis
+            eval, evec = np.linalg.eigh(F)
+            self.eps.append(eval)
+            self.L.append(evec)
+            self.dim[ij] = Q_ortho.shape[1] 
+            T2_local += self.dim[ij] * self.dim[ij]
+            print(self.local + " dimension of pair %d = %d" % (ij, self.dim[ij]))
+
+        print("Average " + self.local + " dimension: %2.3f" % (np.average(self.dim)))
+        print("T2 " +  self.local + ": %d" % (T2_local))
+        T2_full = (self.no*self.no)*(self.nv*self.nv)
+        print("T2 full: %d" % (T2_full))
+        print("T2 Ratio: %3.12f" % (T2_local/T2_full))
+
+        #temporary way to generate make sure the phase factor of Q_ij and L_ij matches with Q_ji and L_ji
+        for i in range(self.no):
+            for j in range(0,i):
+                ij = i*self.no + j
+                ji = j*self.no + i
+
+                self.Q[ji] = self.Q[ij]
+                self.L[ji] = self.L[ij]
+
     def _pert_pairdensity(self,t2):
         '''
          Constructing the approximated perturbed pair density
@@ -483,7 +548,7 @@ def _pairdensity(self, t_ijab):
             D[ij] *= 0.5
         return D
 
-    def QL_tensors(self,v,local,t2,D):
+    def QL_tensors(self,v,t2,D,local):
         # Create list for Q, L and eps
         Q_full = np.zeros_like(t2.copy().reshape((self.no*self.no,self.nv,self.nv)))
         Q = []  # truncated PNO list
@@ -496,7 +561,7 @@ def QL_tensors(self,v,local,t2,D):
         for ij in range(self.no*self.no):
             i = ij // self.no
             j = ij % self.no
-  
+
             # Compute local and truncate
             occ[ij], Q_full[ij] = np.linalg.eigh(D[ij])
             if (occ[ij] < 0).any(): # Check for negative occupation numbers
@@ -505,7 +570,7 @@ def QL_tensors(self,v,local,t2,D):
                       Using absolute values - please check if your input is correct.".format(neg))
             dim[ij] = (np.abs(occ[ij]) > self.cutoff).sum()
             Q.append(Q_full[ij, :, (self.nv-dim[ij]):])
-
+            
             # Compute semicanonical virtual space
             F = Q[ij].T @ self.H.F[v,v] @ Q[ij]  # Fock matrix in local basis
             eval, evec = np.linalg.eigh(F)
@@ -780,6 +845,41 @@ def filter_amps(self, r1, r2):
 
         return t1, t2
 
+    def filter_pertamps(self, r1, r2, eps_occ, eps_vir, omega):
+        no = self.no
+        nv = self.nv
+        dim = self.dim
+
+        t1 = np.zeros((no,nv))
+        for i in range(no):
+            ii = i * no + i
+
+            X = self.Q[ii].T @ r1[i]
+            Y = self.L[ii].T @ X
+
+            for a in range(dim[ii]):
+                Y[a] = Y[a]/(eps_occ[i] - eps_vir[ii][a] + omega)
+
+            X = self.L[ii] @ Y
+            t1[i] = self.Q[ii] @ X
+
+        t2 = np.zeros((no,no,nv,nv))
+        for ij in range(no*no):
+            i = ij // no
+            j = ij % no
+
+            X = self.Q[ij].T @ r2[i,j] @ self.Q[ij]
+            Y = self.L[ij].T @ X @ self.L[ij]
+
+            for a in range(dim[ij]):
+                for b in range(dim[ij]):
+                    Y[a,b] = Y[a,b]/(eps_occ[i] + eps_occ[j] - eps_vir[ij][a] - eps_vir[ij][b] + omega)
+
+            X = self.L[ij] @ Y @ self.L[ij].T
+            t2[i,j] = self.Q[ij] @ X @ self.Q[ij].T
+
+        return t1, t2
+
     def filter_res(self, r1, r2):
         no = self.no
         nv = self.nv

pycc/tests/test_036_cpnoppcc.py

@@ -0,0 +1,80 @@
+"""
+Test basic CPNO++-CCSD energy
+"""
+
+# Import package, test suite, and other packages as needed
+import psi4
+import pycc
+import pytest
+from ..data.molecules import *
+
+
+def test_cpnopp_ccsd():
+    """H2O CPNO++-CCSD Test"""
+    # Psi4 Setup
+    psi4.set_memory('2 GB')
+    psi4.core.set_output_file('output.dat', False)
+    psi4.set_options({'basis': 'cc-pVDZ',
+                      'scf_type': 'pk',
+                      'mp2_type': 'conv',
+                      'freeze_core': 'false',
+                      'e_convergence': 1e-13,
+                      'd_convergence': 1e-13,
+                      'r_convergence': 1e-13,
+                      'diis': 8})
+    mol = psi4.geometry(moldict["H2O"])
+    rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
+
+    maxiter = 75
+    e_conv = 1e-12
+    r_conv = 1e-12
+    max_diis = 8
+
+    ccsd = pycc.ccwfn(rhf_wfn, local='CPNO++', local_cutoff=1e-7, it2_opt=False)
+    eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter)
+
+    hbar = pycc.cchbar(ccsd)
+    cclambda = pycc.cclambda(ccsd, hbar)
+    lccsd = cclambda.solve_lambda(e_conv, r_conv, maxiter, max_diis)
+
+    #Ruhee's ccsd_lpno code 
+    esim = -0.22303320613504354
+    lsim = -0.21890326836263854
+
+    assert (abs(esim - eccsd) < 1e-7)
+    assert (abs(lsim - lccsd) < 1e-7)   
+
+def test_pnopp_ccsd_opt():
+    """H2O CPNO++-CCSD with Optimized Initial T2 Amplitudes"""
+    # Psi4 Setup
+    psi4.set_memory('2 GB')
+    psi4.core.set_output_file('output.dat', False)
+    psi4.set_options({'basis': 'cc-pVDZ',
+                      'scf_type': 'pk',
+                      'mp2_type': 'conv',
+                      'freeze_core': 'false',
+                      'e_convergence': 1e-13,
+                      'd_convergence': 1e-13,
+                      'r_convergence': 1e-13,
+                      'diis': 8})
+    mol = psi4.geometry(moldict["H2O"])
+    rhf_e, rhf_wfn = psi4.energy('SCF', return_wfn=True)
+
+    maxiter = 75
+    e_conv = 1e-12
+    r_conv = 1e-12
+    max_diis = 8
+
+    ccsd = pycc.ccwfn(rhf_wfn, local='CPNO++', local_cutoff=1e-7)
+    eccsd = ccsd.solve_cc(e_conv, r_conv, maxiter)
+
+    hbar = pycc.cchbar(ccsd)
+    cclambda = pycc.cclambda(ccsd, hbar)
+    lccsd = cclambda.solve_lambda(e_conv, r_conv, maxiter, max_diis)
+
+    #Comparing against simulation code 
+    esim = -0.223866820104919 
+    lsim = -0.21966259490352782 
+
+    assert (abs(esim - eccsd) < 1e-7)
+    assert (abs(lsim - lccsd) < 1e-7)