Orbitals should be multiplied by sqrt(occ) not occ

psi4/src/psi4/libscf_solver/sad.cc

@@ -523,7 +523,8 @@ void SADGuess::get_uhf_atomic_density(std::shared_ptr<BasisSet> bas, std::shared
     int iteration = 0;
 
     // Setup DIIS
-    DIISManager diis_manager(6, "SAD DIIS", DIISManager::RemovalPolicy::LargestError, DIISManager::StoragePolicy::InCore);
+    DIISManager diis_manager(6, "SAD DIIS", DIISManager::RemovalPolicy::LargestError,
+                             DIISManager::StoragePolicy::InCore);
     diis_manager.set_error_vector_size(gradient_a.get(), gradient_b.get());
     diis_manager.set_vector_size(Fa.get(), Fb.get());
 
@@ -711,9 +712,9 @@ void SADGuess::form_C_and_D(SharedMatrix X, SharedMatrix F, SharedMatrix C, Shar
     for (int i = 0; i < nbf; i++) {
         C_DCOPY(nocc, Cp[i], 1, Coccp[i], 1);
     }
-    // Scale by occ
+    // Scale by sqrt(occ)
     for (int i = 0; i < nocc; i++) {
-        C_DSCAL(nbf, occ->get(i), &Coccp[0][i], nocc);
+        C_DSCAL(nbf, std::sqrt(occ->get(i)), &Coccp[0][i], nocc);
     }
     // Form D = Cocc*Cocc'
     D->gemm(false, true, 1.0, Cocc, Cocc, 0.0);

tests/cc10/input.dat

@@ -13,7 +13,8 @@ set {
   basis       cc-pVDZ
   docc        [4, 0, 1, 1]
   socc        [1, 0, 0, 0]
-  freeze_core = true
+  freeze_core true
+  d_convergence 10
 }
 
 energy('ccsd')

tests/dftd3/nbody-cp-gradient/input.dat

@@ -117,8 +117,8 @@ for result in [  #TEST
     compare_values(ref_2cp_ene + ref_disp_2cp_ene, result, 8, 'CP-Corrected Energy')  #TEST
 
 for result in [  #TEST
-    wfn.variable("CP-CORRECTED 2-BODY INTERACTION ENERGY"),  #TEST
-    variable("CP-CORRECTED 2-BODY INTERACTION ENERGY"),  #TEST
+    wfn.variable("CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY"),  #TEST
+    variable("CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY"),  #TEST
 ]:  #TEST
     compare_values(ref_2cp_ie + ref_disp_2cp_ie, result, 8, 'CP-Corrected Interaction Energy')  #TEST
 
@@ -128,7 +128,7 @@ for result in [  #TEST
     wfn.variable("CURRENT GRADIENT"),  #TEST
     variable("CURRENT GRADIENT"),  #TEST
 ]:  #TEST
-    compare_values(ref_2cp_grad + ref_disp_2cp_grad, result, 8, 'CP-Corrected Gradient')  #TEST
+    compare_values(ref_2cp_grad + ref_disp_2cp_grad, result, 7, 'CP-Corrected Gradient')  #TEST
 
 # compute pieces
 cp_scheme = {  #TEST
@@ -191,7 +191,7 @@ for result in [  #TEST
     compare_values(ref_2cp_ene + ref_disp_2cp_ene, result, 8, 'CP-Corrected Energy')  #TEST
 
 for result in [  #TEST
-    wfn.variable("CP-CORRECTED 2-BODY INTERACTION ENERGY"),  #TEST
+    wfn.variable("CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY"),  #TEST
 ]:  #TEST
     compare_values(ref_2cp_ie + ref_disp_2cp_ie, result, 8, 'CP-Corrected Interaction Energy')  #TEST
 
@@ -200,7 +200,7 @@ for result in [  #TEST
     wfn.variable("CURRENT GRADIENT"),  #TEST
     grad,  #TEST
 ]:  #TEST
-    compare_values(ref_2cp_grad + ref_disp_2cp_grad, result, 8, 'CP-Corrected Gradient')  #TEST
+    compare_values(ref_2cp_grad + ref_disp_2cp_grad, result, 7, 'CP-Corrected Gradient')  #TEST
 
 
 # compute nbody again with return_total_data=False
@@ -217,8 +217,8 @@ for result in [  #TEST
     compare_values(ref_2cp_ene + ref_disp_2cp_ene, result, 8, 'rtd=F CP-Corrected Energy')  #TEST
 
 for result in [  #TEST
-    wfn.variable("CP-CORRECTED 2-BODY INTERACTION ENERGY"),  #TEST
-    variable("CP-CORRECTED 2-BODY INTERACTION ENERGY"),  #TEST
+    wfn.variable("CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY"),  #TEST
+    variable("CP-CORRECTED INTERACTION ENERGY THROUGH 2-BODY"),  #TEST
     wfn.variable("CURRENT ENERGY"),  #TEST
     variable("CURRENT ENERGY"),  #TEST
 ]:  #TEST
@@ -227,11 +227,11 @@ for result in [  #TEST
 for result in [  #TEST
     wfn.variable('GRADIENT 2'),  #TEST
 ]:  #TEST
-    compare_values(ref_2cp_grad + ref_disp_2cp_grad, result, 8, 'rtd=F CP-Corrected Gradient')  #TEST
+    compare_values(ref_2cp_grad + ref_disp_2cp_grad, result, 7, 'rtd=F CP-Corrected Gradient')  #TEST
 
 for result in [  #TEST
     g,  #TEST
     wfn.variable("CURRENT GRADIENT"),  #TEST
     variable("CURRENT GRADIENT"),  #TEST
 ]:  #TEST
-    compare_values(ref_2cp_grad_ie + ref_disp_2cp_grad_ie, result, 8, 'rtd=F CP-Corrected Interaction Gradient')  #TEST
+    compare_values(ref_2cp_grad_ie + ref_disp_2cp_grad_ie, result, 7, 'rtd=F CP-Corrected Interaction Gradient')  #TEST

tests/linK-3/input.dat

@@ -36,9 +36,11 @@ set {
 
 uhf_energy = energy('scf')
 psi4.compare_values(uhf_energy, uhf_ref, 8, "UHF Energy (using LinK algo)")
-compare(1, variable("SCF ITERATIONS") < 19.0, "UHF LinK Incfock Efficient")
+print(f'{variable("SCF ITERATIONS")=}')
+compare(1, variable("SCF ITERATIONS") <= 19.0, "UHF LinK Incfock Efficient")
 
 set reference rohf
 rohf_energy = energy('scf')
 psi4.compare_values(rohf_energy, rohf_ref, 8, "ROHF Energy (using LinK algo)")
-compare(1, variable("SCF ITERATIONS") < 19.0, "ROHF LinK Incfock Efficient")
+print(f'{variable("SCF ITERATIONS")=}')
+compare(1, variable("SCF ITERATIONS") <= 19.0, "ROHF LinK Incfock Efficient")

tests/nbody-cp-gradient/input.dat

@@ -103,7 +103,7 @@ for result in [  #TEST
     wfn.variable("CURRENT GRADIENT"),  #TEST
     variable("CURRENT GRADIENT"),  #TEST
 ]:  #TEST
-    compare_values(ref_2cp_grad, result, 8, 'CP-Corrected Gradient')  #TEST
+    compare_values(ref_2cp_grad, result, 7, 'CP-Corrected Gradient')  #TEST
 
 # compute pieces
 cp_scheme = {  #TEST

tests/pytests/test_ddx.py

@@ -310,7 +310,7 @@ def test_ddx_rhf_consistency(inp):
     for key in inp["ddx"].keys():
         psi4.set_options({"ddx_" + key: inp["ddx"][key]})
     scf_e, wfn = psi4.energy('SCF', return_wfn=True, molecule=mol)
-    assert compare_values(inp["ref"], scf_e, 9, "Total SCF energy with DDX versus reference data")
+    assert compare_values(inp["ref"], scf_e, 8, "Total SCF energy with DDX versus reference data")
 
 
 @uusing("ddx")

tests/pytests/test_erisieve.py

@@ -230,7 +230,7 @@ def test_schwarz_vs_density_quartets():
     density_computed_shells = density_wfn.jk().computed_shells_per_iter("Quartets")
 
     schwarz_computed_shells_expected = [20290, 20290, 20290, 20290, 20290, 20290, 20290, 20290, 20290]
-    density_computed_shells_expected = [13171, 19618, 19665, 19657, 19661, 19661, 19663, 19663, 19663]
+    density_computed_shells_expected = [13187, 19683, 19644, 19663, 19661, 19661, 19663, 19663, 19663]
 
     # compare iteration counts of runs with computed shell quartet array lengths
     # iteration_+1 is used to account for computed_shells arrays including SAD guess results
@@ -290,7 +290,7 @@ def test_rhf_vs_uhf_screening():
     rhf_computed_shells = rhf_wfn.jk().computed_shells_per_iter("Quartets")
     uhf_computed_shells = uhf_wfn.jk().computed_shells_per_iter("Quartets")
 
-    computed_shells_expected = [13171, 19618, 19665, 19657, 19661, 19661, 19663, 19663, 19663]
+    computed_shells_expected = [13187, 19683, 19644, 19663, 19661, 19661, 19663, 19663, 19663]
 
     # compare iteration counts of runs with computed shell quartet array lengths
     # iteration_+1 is used to account for computed_shells arrays including SAD guess results

tests/pytests/test_large_atoms.py

@@ -12,29 +12,35 @@
         'Ba2': {
             'grad': np.array([[0, 0,  2.08544743e-02],
                 [ 0,  0, -2.08544743e-02]]),
-            'energy': -50.12131931846057},
+            'energy': -50.12131931846057,
+            'docc' : [ 3, 0, 1, 1, 0, 3, 1, 1 ]},
         'In2': {
             'grad': np.array([[0, 0,  -1.46964113e-02],
                 [ 0,  0, 1.46964113e-02]]),
-            'energy': -378.3986291158715},
+            'energy': -378.3986291158715,
+            'docc' : [ 6, 1, 2, 2, 1, 5, 2, 2 ]},
         'Sr2': {
             'grad': np.array([[0, 0,  9.74576806e-03],
                 [ 0,  0, -9.74576806e-03]]),
-            'energy': -60.66907844700816},
+            'energy': -60.66907844700816,
+            'docc': [ 3, 0, 1, 1, 0, 3, 1, 1 ]}
      },
      'mp2': {
         'Ba2': {
             'grad': np.array([[0, 0, 1.6421392746403734e-02],
                      [0, 0, -1.6421392746403734e-02]]),
-            'energy': -50.34353063117147},
+            'energy': -50.34353063117147,
+            'docc' : [ 3, 0, 1, 1, 0, 3, 1, 1 ]},
         'In2': {
             'grad': np.array([[0, 0, -1.4464867319381867e-02],
                      [0, 0, 1.4464867319381867e-02]]),
-            'energy': -378.5493603721844},
+            'energy': -378.5493603721844,
+            'docc' : [ 6, 1, 2, 2, 1, 5, 2, 2 ]},
         'Sr2': {
             'grad': np.array([[0, 0, 6.469022590311186e-03],
                      [0, 0, -6.469022590311186e-03]]),
-            'energy': -60.83705034568311},
+            'energy': -60.83705034568311,
+            'docc': [ 3, 0, 1, 1, 0, 3, 1, 1 ]}
      }
 }
 
@@ -79,6 +85,11 @@ def test_large_atoms(inp, mol, mols):
 
     psi4.set_options(inp['options'])
     psi4.set_options({'points': 5})
+    # The atoms are so far apart SCF convergence can land on different
+    # solutions, and this is why it is important to set the wanted
+    # occupations
+    psi4.set_options({'docc' : inp['ref'][mol]['docc']})
+
     method = inp['method']
     analytic_grad, wfn = psi4.gradient(method, molecule=mols[mol], dertype=1, return_wfn=True)
     findif_grad = psi4.gradient(method, molecule=mols[mol], dertype=0)