A/E DIIS (psi4#2320)

* A/EDIIS

.azure-pipelines/azure-pipelines-linux.yml

@@ -115,7 +115,8 @@ jobs:
         conda-forge::qcengine \
         conda-forge::pymdi \
         pandas!=1.3.0 \
-        adcc::adcc
+        adcc::adcc \
+        scipy
       source activate p4env
       which python
       pip install git+https://github.com/i-pi/i-pi.git@master-py3

.azure-pipelines/azure-pipelines-windows.yml

@@ -135,7 +135,8 @@ jobs:
                         conda-forge::mpmath ^
                         conda-forge::msgpack-python ^
                         conda-forge::qcelemental ^
-                        conda-forge::qcengine
+                        conda-forge::qcengine ^
+                        scipy
           pip install git+https://github.com/i-pi/i-pi.git@master-py3
           conda install --only-deps anaconda-client
           pip install git+https://github.com/loriab/anaconda-client.git@upload-catch-silent

doc/sphinxman/source/scf.rst

@@ -579,17 +579,23 @@ that |PSIfour| expects the numpy file on disk to have the ``.npy`` extension, no
 Convergence Stabilization
 ~~~~~~~~~~~~~~~~~~~~~~~~~
 
-With regard to convergence stabilization, Pulay's Direct Inversion of the
-Iterative Subspace (DIIS) extrapolation,  Gill's Maximum Overlap Method (MOM),
-and damping are all implemented. A summary of each is presented below,
+A summary of Psi's supported convergence stabilization techniques is presented below:
 
 DIIS [On by Default]
-    DIIS uses previous iterates of the Fock Matrix together
+    DIIS uses previous iterates of the Fock matrix together
     with an error criterion based on the orbital gradient to produce an informed
     estimate of the next Fock Matrix. DIIS is almost always necessary to converge
     the SCF procedure and is therefore turned on by default. In rare cases, the
     DIIS algorithm may need to be modified or turned off altogether, which may be
     accomplished via :term:`options <DIIS (SCF)>`.
+ADIIS [On by Default]
+    ADIIS uses previous iterates of the Fock and density matrices to produce an
+    informed estimate of the next Fock matrix. ADIIS estimates are based on minimizing
+    an energy estimate rather than zeroing the residual, so this performs best in the early
+    iterations. By default, Psi will start using ADIIS before blending the ADIIS step with
+    the DIIS step, eventually using the pure DIIS step. The closely-related EDIIS procedure
+    may be used instead by setting |scf__scf_initial_accelerator|. This is formally identical
+    to ADIIS for HF, but the methods will differ for more general DFT.
 MOM [Off by Default]
     MOM was developed to combat a particular class of convergence failure:
     occupation flipping. In some cases, midway though the SCF procedure, a partially

psi4/driver/procrouting/scf_proc/scf_iterator.py

@@ -1,4 +1,4 @@
-#
+
 # @BEGIN LICENSE
 #
 # Psi4: an open-source quantum chemistry software package
@@ -254,7 +254,7 @@ def scf_iterate(self, e_conv=None, d_conv=None):
     reference = core.get_option('SCF', "REFERENCE")
 
     # self.member_data_ signals are non-local, used internally by c-side fns
-    self.diis_enabled_ = _validate_diis()
+    self.diis_enabled_ = self.validate_diis()
     self.MOM_excited_ = _validate_MOM()
     self.diis_start_ = core.get_option('SCF', 'DIIS_START')
     damping_enabled = _validate_damping()
@@ -400,10 +400,8 @@ def scf_iterate(self, e_conv=None, d_conv=None):
                 Dnorm = self.compute_orbital_gradient(add_to_diis_subspace, core.get_option('SCF', 'DIIS_MAX_VECS'))
 
                 if add_to_diis_subspace:
-                    diis_performed = self.diis()
-
-                if diis_performed:
-                    status.append("DIIS")
+                    for engine_used in self.diis(Dnorm):
+                        status.append(engine_used)
 
                 core.timer_off("HF: DIIS")
 
@@ -823,7 +821,7 @@ def _validate_damping():
     return enabled
 
 
-def _validate_diis():
+def _validate_diis(self):
     """Sanity-checks DIIS control options
 
     Raises
@@ -834,10 +832,24 @@ def _validate_diis():
     Returns
     -------
     bool
-        Whether DIIS is enabled during scf.
+        Whether some form of DIIS is enabled during SCF.
 
     """
-    enabled = bool(core.get_option('SCF', 'DIIS'))
+
+    restricted_open = self.same_a_b_orbs() and not self.same_a_b_dens()
+    aediis_active = core.get_option('SCF', 'SCF_INITIAL_ACCELERATOR') != "NONE" and not restricted_open
+
+    if aediis_active:
+        start = core.get_option('SCF', 'SCF_INITIAL_START_DIIS_TRANSITION')
+        stop = core.get_option('SCF', 'SCF_INITIAL_FINISH_DIIS_TRANSITION')
+        if start < stop:
+            raise ValidationError('SCF_INITIAL_START_DIIS_TRANSITION error magnitude cannot be less than SCF_INITIAL_FINISH_DIIS_TRANSITION.')
+        elif start < 0:
+            raise ValidationError('SCF_INITIAL_START_DIIS_TRANSITION cannot be negative.')
+        elif stop < 0:
+            raise ValidationError('SCF_INITIAL_FINISH_DIIS_TRANSITION cannot be negative.')
+
+    enabled = bool(core.get_option('SCF', 'DIIS')) or aediis_active
     if enabled:
         start = core.get_option('SCF', 'DIIS_START')
         if start < 1:
@@ -927,6 +939,7 @@ def _validate_soscf():
 
     return enabled
 
+core.HF.validate_diis = _validate_diis
 
 def efp_field_fn(xyz):
     """Callback function for PylibEFP to compute electric field from electrons

psi4/driver/procrouting/scf_proc/subclass_methods.py

@@ -3,18 +3,35 @@
 from psi4 import core
 from ..diis import DIIS, StoragePolicy, RemovalPolicy
 
+def diis_engine_helper(self):
+    engines = set()
+    if core.get_option('SCF', 'DIIS'):
+        engines.add('diis')
+    restricted_open = self.same_a_b_orbs() and not self.same_a_b_dens()
+    if not restricted_open:
+        aediis = core.get_option('SCF', 'SCF_INITIAL_ACCELERATOR')
+        if aediis != "NONE":
+            engines.add(aediis.lower())
+    return engines
+
 def _RHF_orbital_gradient(self, save_fock: bool, max_diis_vectors: int) -> float:
     gradient = self.form_FDSmSDF(self.Fa(), self.Da())
 
     if save_fock:
         if not self.initialized_diis_manager_:
             storage_policy = StoragePolicy.InCore if self.scf_type() == "DIRECT" else StoragePolicy.OnDisk
-            self.diis_manager_ = DIIS(max_diis_vectors, "HF DIIS vector", RemovalPolicy.LargestError, storage_policy)
-            self.diis_manager_.set_error_vector_size(gradient)
-            self.diis_manager_.set_vector_size(self.Fa())
+            self.diis_manager_ = DIIS(max_diis_vectors, "HF DIIS vector", RemovalPolicy.LargestError, storage_policy, engines=diis_engine_helper(self))
             self.initialized_diis_manager_ = True
-
-        self.diis_manager_.add_entry(gradient, self.Fa())
+
+        entry = {"target": [self.Fa()]}
+        if core.get_option('SCF', 'DIIS'):
+            entry["error"] = [gradient]
+        aediis = core.get_option('SCF', 'SCF_INITIAL_ACCELERATOR')
+        if aediis != "NONE":
+            entry["densities"] = [self.Da()]
+            if aediis == "EDIIS":
+                entry["energy"] = [self.compute_E()]
+        self.diis_manager_.add_entry(entry)
 
     if self.options().get_bool("DIIS_RMS_ERROR"):
         return gradient.rms()
@@ -28,12 +45,18 @@ def _UHF_orbital_gradient(self, save_fock: bool, max_diis_vectors: int) -> float
     if save_fock:
         if not self.initialized_diis_manager_:
             self.diis_manager_ = DIIS(max_diis_vectors, "HF DIIS vector", RemovalPolicy.LargestError,
-                                                          StoragePolicy.OnDisk)
-            self.diis_manager_.set_error_vector_size(gradient_a, gradient_b)
-            self.diis_manager_.set_vector_size(self.Fa(), self.Fb())
+                                                          StoragePolicy.OnDisk, False, engines=diis_engine_helper(self))
             self.initialized_diis_manager_ = True
 
-        self.diis_manager_.add_entry(gradient_a, gradient_b, self.Fa(), self.Fb())
+        entry = {"target": [self.Fa(), self.Fb()]}
+        if core.get_option('SCF', 'DIIS'):
+            entry["error"] = [gradient_a, gradient_b]
+        aediis = core.get_option('SCF', 'SCF_INITIAL_ACCELERATOR')
+        if aediis != "NONE":
+            entry["densities"] = [self.Da(), self.Db()]
+            if aediis == "EDIIS":
+                entry["energy"] = [self.compute_E()]
+        self.diis_manager_.add_entry(entry)
 
     if self.options().get_bool("DIIS_RMS_ERROR"):
         return math.sqrt(0.5 * (gradient_a.rms() ** 2 + gradient_b.rms() ** 2))
@@ -68,12 +91,12 @@ def _ROHF_orbital_gradient(self, save_fock: bool, max_diis_vectors: int) -> floa
 
     if save_fock:
         if not self.initialized_diis_manager_:
-            self.diis_manager_ = DIIS(max_diis_vectors, "HF DIIS vector", RemovalPolicy.LargestError, StoragePolicy.OnDisk)
+            self.diis_manager_ = DIIS(max_diis_vectors, "HF DIIS vector", RemovalPolicy.LargestError, StoragePolicy.OnDisk, engines=diis_engine_helper(self))
             self.diis_manager_.set_error_vector_size(gradient)
             self.diis_manager_.set_vector_size(self.soFeff())
             self.initialized_diis_manager_ = True
 
-        self.diis_manager_.add_entry(gradient, self.soFeff())
+        self.diis_manager_.add_entry({"error": [gradient], "target": [self.soFeff()]})
 
     if self.options().get_bool("DIIS_RMS_ERROR"):
         return gradient.rms()
@@ -83,3 +106,16 @@ def _ROHF_orbital_gradient(self, save_fock: bool, max_diis_vectors: int) -> floa
 core.RHF.compute_orbital_gradient = _RHF_orbital_gradient
 core.UHF.compute_orbital_gradient = core.CUHF.compute_orbital_gradient = _UHF_orbital_gradient
 core.ROHF.compute_orbital_gradient = _ROHF_orbital_gradient
+
+def _RHF_diis(self, Dnorm):
+    return self.diis_manager_.extrapolate(self.Fa(), Dnorm=Dnorm)
+
+def _UHF_diis(self, Dnorm):
+    return self.diis_manager_.extrapolate(self.Fa(), self.Fb(), Dnorm=Dnorm)
+
+def _ROHF_diis(self, Dnorm):
+    return self.diis_manager_.extrapolate(self.soFeff(), Dnorm=Dnorm)
+
+core.RHF.diis = _RHF_diis
+core.UHF.diis = core.CUHF.diis = _UHF_diis
+core.ROHF.diis = _ROHF_diis

psi4/src/psi4/libdiis/diismanager.h

@@ -80,8 +80,7 @@ class PSI_API DIISManager {
     };
     template <typename... types>
     bool extrapolate(types... arrays) {
-        auto success = pydiis.attr("extrapolate")(arrays...);
-        return success.template cast<bool>();
+        return py::len(pydiis.attr("extrapolate")(arrays...));
     };
     template <typename ... types>
     bool add_entry(types... arrays) {

psi4/src/psi4/libscf_solver/cuhf.cc

@@ -388,8 +388,6 @@ double CUHF::compute_E() {
     return Etotal;
 }
 
-bool CUHF::diis() { return diis_manager_.attr("extrapolate")(Fa_.get(), Fb_.get()).cast<bool>(); }
-
 bool CUHF::stability_analysis() {
     throw PSIEXCEPTION("CUHF stability analysis has not been implemented yet.  Sorry :(");
     return false;

psi4/src/psi4/libscf_solver/cuhf.h

@@ -92,7 +92,6 @@ class CUHF final : public HF {
          std::shared_ptr<PSIO> psio);
     ~CUHF() override;
 
-    bool diis() override;
     void save_density_and_energy() override;
 
     void form_C(double shift = 0.0) override;

psi4/src/psi4/libscf_solver/hf.h

@@ -320,7 +320,7 @@ class HF : public Wavefunction {
     void find_occupation();
 
     /** Performs DIIS extrapolation */
-    virtual bool diis() { return false; }
+    virtual bool diis(double dnorm) { return false; }
 
     /** Compute the orbital gradient */
     virtual double compute_orbital_gradient(bool save_diis, int max_diis_vectors) { return 0.0; }

psi4/src/psi4/libscf_solver/rhf.cc

@@ -233,8 +233,6 @@ void RHF::form_G() {
     }
 }
 
-bool RHF::diis() { return diis_manager_.attr("extrapolate")(Fa_.get()).cast<bool>(); }
-
 void RHF::form_F() {
     Fa_->copy(H_);
     Fa_->add(G_);

psi4/src/psi4/libscf_solver/rhf.h

@@ -60,7 +60,6 @@ class RHF : public HF {
     virtual bool same_a_b_orbs() const { return true; }
     virtual bool same_a_b_dens() const { return true; }
 
-    bool diis() override;
     void save_density_and_energy() override;
 
     void form_C(double shift = 0.0) override;

psi4/src/psi4/libscf_solver/rohf.cc

@@ -290,8 +290,6 @@ void ROHF::save_density_and_energy() {
     Dt_old_->copy(Dt_);
 }
 
-bool ROHF::diis() { return diis_manager_.attr("extrapolate")(soFeff_.get()).cast<bool>(); }
-
 void ROHF::form_initial_F() {
     // Form the initial Fock matrix, closed and open variants
     Fa_->copy(H_);

psi4/src/psi4/libscf_solver/rohf.h

@@ -77,7 +77,6 @@ class ROHF : public HF {
     SharedMatrix moFb() const { return moFb_; }
     SharedMatrix Ct() const {return Ct_; }
 
-    bool diis() override;
     void save_density_and_energy() override;
 
     void form_C(double shift = 0.0) override;

psi4/src/psi4/libscf_solver/uhf.cc

@@ -1089,8 +1089,6 @@ int UHF::soscf_update(double soscf_conv, int soscf_min_iter, int soscf_max_iter,
     return cphf_nfock_builds_;
 }
 
-bool UHF::diis() { return diis_manager_.attr("extrapolate")(Fa_.get(), Fb_.get()).cast<bool>(); }
-
 bool UHF::stability_analysis() {
     if (functional_->needs_xc()) {
         throw PSIEXCEPTION("Stability analysis not yet supported for XC functionals.");

psi4/src/psi4/libscf_solver/uhf.h

@@ -72,7 +72,6 @@ class UHF : public HF {
     virtual bool same_a_b_orbs() const { return false; }
     virtual bool same_a_b_dens() const { return false; }
 
-    bool diis() override;
     void save_density_and_energy() override;
 
     void form_C(double shift = 0.0) override;

psi4/src/read_options.cc

@@ -1394,7 +1394,7 @@ int read_options(const std::string &name, Options &options, bool suppress_printi
         options.add_str("DF_BASIS_GUESS", "FALSE", "");
         /*- Use RMS error instead of the more robust absolute error? -*/
         options.add_bool("DIIS_RMS_ERROR", true);
-        /*- The minimum iteration to start storing DIIS vectors -*/
+        /*- The minimum iteration to start storing DIIS vectors and performing ADIIS/EDIIS. -*/
         options.add_int("DIIS_START", 1);
         /*- Minimum number of error vectors stored for DIIS extrapolation. Will be removed in v1.7. -*/
         options.add_int("DIIS_MIN_VECS", 2);
@@ -1439,6 +1439,15 @@ int read_options(const std::string &name, Options &options, bool suppress_printi
         /*- When using |scf__stability_analysis| ``FOLLOW``, maximum number of orbital optimization attempts
             to make the wavefunction stable. !expert -*/
         options.add_int("MAX_ATTEMPTS", 1);
+        /*- Use a method to accelerate initial SCF convergence? Use ``NONE`` for DIIS alone (if enabled) and ``EDIIS`` or ``ADIIS``
+            to have both the chosen accelerator and DIIS (if enabled). For restricted-open references, ``EDIIS`` and ``ADIIS`` have no effect. -*/
+        options.add_str("SCF_INITIAL_ACCELERATOR", "ADIIS", "NONE EDIIS ADIIS");
+        /*- SCF error at which to start the linear interpolation between DIIS steps and steps of the initial SCF accelerator.
+            Value taken from Garza and Scuseria, DOI: 10.1063/1.4740249 -*/
+        options.add_double("SCF_INITIAL_START_DIIS_TRANSITION", 1.0E-1);
+        /*- SCF error at which to complete the linear interpolation between DIIS steps and steps of the initial SCF accelerator
+            Value taken from Garza and Scuseria, DOI: 10.1063/1.4740249 -*/
+        options.add_double("SCF_INITIAL_FINISH_DIIS_TRANSITION", 1.0E-4);
         /*- Do Perform Incremental Fock Build? -*/
         options.add_bool("INCFOCK", false);
         /*- Frequency with which to compute the full Fock matrix if using |scf__incfock| . 

tests/CMakeLists.txt

@@ -33,7 +33,8 @@ find_package(Perl QUIET)
 #set(py3_fail_list pywrap-freq-e-sowreap
 #                  pywrap-freq-g-sowreap pywrap-opt-sowreap
 #                  pywrap-db2)
-foreach(test_name adc1 adc2 casscf-fzc-sp casscf-semi casscf-sa-sp ao-casscf-sp casscf-sp castup1
+foreach(test_name adc1 adc2 aediis-1 aediis-2
+                  casscf-fzc-sp casscf-semi casscf-sa-sp ao-casscf-sp casscf-sp castup1
                   castup2 castup3 cbs-delta-energy cbs-parser cbs-xtpl-alpha cbs-xtpl-energy
                   cbs-xtpl-freq cbs-xtpl-gradient cbs-xtpl-opt cbs-xtpl-func cbs-xtpl-nbody
                   cbs-xtpl-wrapper cbs-xtpl-dict cc1 cc10 cc11 cc12 cc13 cc13a cc13b cc13c

tests/aediis-1/CMakeLists.txt

@@ -0,0 +1,3 @@
+include(TestingMacros)
+
+add_regression_test(aediis-1 "psi;scf;aediis")

tests/aediis-1/input.dat

@@ -0,0 +1,24 @@
+#! ADIIS test case, from 10.1063/1.3304922 
+
+molecule {
+2 1
+Cd       0.0000000000      0.0000000000      0.0000000000    
+N        0.0000000000      0.0000000000     -2.2600010635    
+N       -0.6854437203      0.0000000000     -4.3480354859    
+C        0.6760533303      0.0000000000     -4.3850694419    
+C        1.0852404753      0.0000000000     -3.0912317094    
+C       -1.0447523879      0.0000000000     -3.0602202054    
+H        1.2315299626      0.0000000000     -5.3007589304    
+H        2.0886410569      0.0000000000     -2.7110772355    
+H       -2.0687499916      0.0000000000     -2.7265149874    
+H       -1.3131700977      0.0000000000     -5.1747180497    
+}
+
+set guess core    # With this guess, un-accelerated DIIS will fail.
+set scf_type pk
+set basis 3-21G
+set scf_initial_accelerator adiis
+set maxiter 15
+energy = energy('scf')
+
+compare_values(-5663.1433914266744978, energy, 6, "SCF Energy")