Merge pull request #120 from BoothGroup/tolerance_fix

Adds occupation tolerance option for DMET cluster DM. Fixes EDMET docs example.

docs/source/quickstart/edmetfinite.py

@@ -26,7 +26,7 @@
 # The KS object needs to be converted to a HF object:
 lda = lda.to_rhf()
 
-emb_lda = vayesta.edmet.EDMET(lda, dmet_threshold=1e-12, solver="EBCCSD",
+emb_lda = vayesta.edmet.EDMET(lda, bath_options=dict(dmet_threshold=1e-12), solver="CCSD-S-1-1",
                               oneshot=True, make_dd_moments=False)
 
 with emb_lda.iao_fragmentation() as f:
@@ -42,7 +42,7 @@
 # The KS opbject needs to be converted to a HF object:
 b3lyp = b3lyp.to_rhf()
 
-emb_b3lyp = vayesta.edmet.EDMET(b3lyp, dmet_threshold=1e-12, solver="EBCCSD", oneshot=True, make_dd_moments=False)
+emb_b3lyp = vayesta.edmet.EDMET(b3lyp, bath_options=dict(dmet_threshold=1e-12), solver="CCSD-S-1-1", oneshot=True, make_dd_moments=False)
 with emb_b3lyp.iao_fragmentation() as f:
     f.add_all_atomic_fragments()
 emb_b3lyp.kernel()
@@ -51,7 +51,7 @@
 hf = pyscf.scf.RHF(mol).density_fit()
 hf.kernel()
 
-emb_hf = vayesta.edmet.EDMET(hf, dmet_threshold=1e-12, solver="EBCCSD",
+emb_hf = vayesta.edmet.EDMET(hf, bath_options=dict(dmet_threshold=1e-12), solver="CCSD-S-1-1",
                              oneshot=True, make_dd_moments=False)
 
 with emb_hf.iao_fragmentation() as f:
@@ -65,7 +65,7 @@
 print("E(LDA)=              %+16.8f Ha" % lda.e_tot)
 print("E(Emb. CCSD @LDA)=   %+16.8f Ha" % emb_lda.e_tot)
 print("E(HF @LDA)=          %+16.8f Ha" % emb_lda.e_mf)
-print("E(CCSD)=             %+16.8f Ha" % cc.)
+print("E(CCSD)=             %+16.8f Ha" % cc.e_tot)
 print("E(B3LYP)=            %+16.8f Ha" % b3lyp.e_tot)
 print("E(HF)=               %+16.8f Ha" % hf.e_tot)
 print("E(HF @B3LYP)=        %+16.8f Ha" % emb_b3lyp.e_mf)

vayesta/core/bath/dmet.py

@@ -47,7 +47,8 @@ def kernel(self):
         # --- Separate occupied and virtual in cluster
         cluster = [self.c_frag, c_dmet]
         self.base._check_orthonormal(*cluster, mo_name='cluster MO')
-        c_cluster_occ, c_cluster_vir = self.fragment.diagonalize_cluster_dm(*cluster, tol=2*self.dmet_threshold)
+        tol = self.base.opts.bath_options['occupation_tolerance']
+        c_cluster_occ, c_cluster_vir = self.fragment.diagonalize_cluster_dm(*cluster, tol=tol)
         # Canonicalize
         c_cluster_occ = self.fragment.canonicalize_mo(c_cluster_occ)[0]
         c_cluster_vir = self.fragment.canonicalize_mo(c_cluster_vir)[0]

vayesta/core/qemb/fragment.py

@@ -517,8 +517,8 @@ def diagonalize_cluster_dm(self, *mo_coeff, dm1=None, norm=2, tol=1e-4):
         sc = np.dot(self.base.get_ovlp(), c_cluster)
         dm = dot(sc.T, dm1, sc)
         e, r = np.linalg.eigh(dm)
-        if tol and not np.allclose(np.fmin(abs(e), abs(e-norm)), 0, atol=tol, rtol=0):
-            raise RuntimeError("Eigenvalues of cluster-DM not all close to 0 or %d:\n%s" % (norm, e))
+        if tol and not np.allclose(np.fmin(abs(e), abs(e-norm)), 0, atol=2*tol, rtol=0):
+            self.log.warn("Eigenvalues of cluster-DM not all close to 0 or %d:\n%s" % (norm, e))
         e, r = e[::-1], r[:,::-1]
         c_cluster = np.dot(c_cluster, r)
         c_cluster = fix_orbital_sign(c_cluster)[0]
@@ -863,7 +863,7 @@ def get_orbitals(occtype):
         def check_occup(mo_coeff, expected):
             occup = self.get_mo_occupation(mo_coeff)
             # RHF
-            atol = self.opts.bath_options['dmet_threshold']
+            atol = self.opts.bath_options['occupation_tolerance']
             if np.ndim(occup[0]) == 0:
                 assert np.allclose(occup, 2*expected, rtol=0, atol=2*atol)
             else:

vayesta/core/qemb/qemb.py

@@ -67,7 +67,7 @@ class Options(OptionsBase):
     # --- Bath options
     bath_options: dict = OptionsBase.dict_with_defaults(
         # General
-        bathtype='dmet', canonicalize=True,
+        bathtype='dmet', canonicalize=True, occupation_tolerance=1e-8,
         # DMET bath
         dmet_threshold=1e-8,
         # R2 bath

vayesta/ewf/fragment.py

@@ -114,7 +114,8 @@ def set_cas(self, iaos=None, c_occ=None, c_vir=None, minao='auto', dmet_threshol
             c_env = frag.get_env_coeff(indices)
             bath = DMET_Bath(self, dmet_threshold=dmet_threshold)
             c_dmet = bath.make_dmet_bath(c_env)[0]
-            c_iao_occ, c_iao_vir = self.diagonalize_cluster_dm(c_iao, c_dmet, tol=2*dmet_threshold)
+            tol = self.opts.bath_options['occupation_tolerance']
+            c_iao_occ, c_iao_vir = self.diagonalize_cluster_dm(c_iao, c_dmet, tol=2*tol)
         else:
             c_iao_occ = c_iao_vir = None