FCI and GENCI initial core guess determinants for Davidson-Liu solver

docs/source/options.rst

@@ -1863,6 +1863,14 @@ Type: int
 
 Default value: 10
 
+**CORE_GUESS**
+
+Integer boolean. Include determinants with single- and none-occupation in the first orbital position as the initial guess space for the DL solver
+
+Type: int_list
+
+Default value: [0]
+
 **PRINT_NO**
 
 Print the NO from the rdm of FCI

forte/base_classes/active_space_method.cc

@@ -78,6 +78,8 @@ void ActiveSpaceMethod::set_r_convergence(double value) { r_convergence_ = value
 
 void ActiveSpaceMethod::set_maxiter(size_t value) { maxiter_ = value; }
 
+void ActiveSpaceMethod::set_core_guess(bool core_guess) { core_guess_ = core_guess; }
+
 void ActiveSpaceMethod::set_read_wfn_guess(bool read) { read_wfn_guess_ = read; }
 
 void ActiveSpaceMethod::set_dump_wfn(bool dump) { dump_wfn_ = dump; }

forte/base_classes/active_space_method.h

@@ -312,6 +312,9 @@ class ActiveSpaceMethod {
     /// @param value the maximum number of iterations
     void set_maxiter(size_t value);
 
+    /// Enable/disable core determinants as initial guess for Davidson-Liu
+    void set_core_guess(bool core_guess);
+
     /// Set if we dump the wave function to disk
     void set_read_wfn_guess(bool read);
 
@@ -379,6 +382,9 @@ class ActiveSpaceMethod {
     /// The maximum number of iterations
     size_t maxiter_ = 100;
 
+    /// Use core determinants as initial guess for Davidson-Liu?
+    bool core_guess_ = false;
+
     /// The root used to compute properties (zero based, default = 0)
     int root_ = 0;
 

forte/base_classes/active_space_solver.cc

@@ -142,6 +142,9 @@ const std::map<StateInfo, std::vector<double>>& ActiveSpaceSolver::compute_energ
         method->set_e_convergence(e_convergence_);
         method->set_r_convergence(r_convergence_);
         method->set_maxiter(maxiter_);
+
+        // set boolean for using core determinants in Davidson-Liu algorithm
+        method->set_core_guess(state.core_guess());
 
         if (read_initial_guess_) {
             state_filename_map_[state] = method->wfn_filename();
@@ -423,6 +426,9 @@ make_state_weights_map(std::shared_ptr<ForteOptions> options,
     // check if the user provided a AVG_STATE list
     py::list avg_state = options->get_gen_list("AVG_STATE");
 
+    bool core_guess;
+    auto core_guess_list = options->get_int_list("CORE_GUESS");
+
     std::vector<size_t> gas_min(6, 0);
     std::vector<size_t> gas_max(6);
     for (int i = 0; i < 6; ++i) {
@@ -458,8 +464,21 @@ make_state_weights_map(std::shared_ptr<ForteOptions> options,
                 gas_max[gasn] = gas_space_max[0];
             }
         }
+
+        // if 'CORE_GUESS' list is given convert valid first 'CORE_GUESS' entry to boolean
+        if (core_guess_list.empty()){
+            core_guess = false;
+        } else if (!(core_guess_list[0] != 0 || core_guess_list[0] != 1)) {
+            psi::outfile->Printf("\n  Error: wrong entry value for CORE_GUESS (%d). "
+                                    "Only values of 0 or 1 are acceptable",
+                                        core_guess_list[0]);
+            throw std::runtime_error("Wrong input value for CORE_GUESS.");
+        } else {
+            core_guess = static_cast<bool>(core_guess_list[0]);
+        }
+
         StateInfo state_this(state.na(), state.nb(), state.multiplicity(), state.twice_ms(),
-                             state.irrep(), state.irrep_label(), gas_min, gas_max);
+                             state.irrep(), state.irrep_label(), gas_min, gas_max, core_guess);
         state_weights_map[state_this] = weights;
     } else {
         double sum_of_weights = 0.0;
@@ -560,8 +579,26 @@ make_state_weights_map(std::shared_ptr<ForteOptions> options,
                 }
             }
 
+            // if 'CORE_GUESS' list is given convert valid 'CORE_GUESS' entries to boolean
+            if (core_guess_list.empty()){
+                core_guess = false;
+            } else if (core_guess_list.size() != nentry) {
+                psi::outfile->Printf("\n  Error: mismatched number of entries in AVG_STATE "
+                                        "(%d) and CORE_GUESS (%d).",
+                                        nentry, core_guess_list.size());
+                throw std::runtime_error(
+                    "Mismatched number of entries in AVG_STATE and CORE_GUESS.");
+            } else if (!(core_guess_list[i] == 0 || core_guess_list[i] == 1)) {
+                psi::outfile->Printf("\n  Error: wrong entry value for CORE_GUESS (%d). "
+                                        "Only values of 0 or 1 are acceptable",
+                                          core_guess_list[i]);
+                throw std::runtime_error("Wrong input value for CORE_GUESS.");
+            } else {
+                core_guess = static_cast<bool>(core_guess_list[i]);
+            }
+
             StateInfo state_this(state.na(), state.nb(), multi, state.twice_ms(), irrep,
-                                 irrep_label, gas_min, gas_max);
+                                 irrep_label, gas_min, gas_max, core_guess);
             state_weights_map[state_this] = weights;
         }
 

forte/base_classes/state_info.cc

@@ -39,9 +39,9 @@ namespace forte {
 
 StateInfo::StateInfo(int na, int nb, int multiplicity, int twice_ms, int irrep,
                      const std::string& irrep_label, const std::vector<size_t> gas_min,
-                     const std::vector<size_t> gas_max)
+                     const std::vector<size_t> gas_max, bool core_guess)
     : na_(na), nb_(nb), multiplicity_(multiplicity), twice_ms_(twice_ms), irrep_(irrep),
-      irrep_label_(irrep_label), gas_min_(gas_min), gas_max_(gas_max) {}
+      irrep_label_(irrep_label), gas_min_(gas_min), gas_max_(gas_max), core_guess_(core_guess) {}
 
 const std::vector<std::string> StateInfo::multiplicity_labels{
     "Singlet", "Doublet", "Triplet", "Quartet", "Quintet", "Sextet", "Septet", "Octet",
@@ -58,6 +58,8 @@ int StateInfo::twice_ms() const { return twice_ms_; }
 
 int StateInfo::irrep() const { return irrep_; }
 
+bool StateInfo::core_guess() const { return core_guess_; }
+
 const std::string& StateInfo::irrep_label() const { return irrep_label_; }
 
 const std::string& StateInfo::multiplicity_label() const {
@@ -71,8 +73,8 @@ const std::vector<size_t>& StateInfo::gas_max() const { return gas_max_; }
 bool StateInfo::operator<(const StateInfo& rhs) const {
     // Make sure the roots are in increasing energy order for core-excited state calcualtions
     if ((gas_min_ == rhs.gas_min_) && (gas_max_ == rhs.gas_max_)) {
-        return std::tie(na_, nb_, multiplicity_, twice_ms_, irrep_) <
-               std::tie(rhs.na_, rhs.nb_, rhs.multiplicity_, rhs.twice_ms_, rhs.irrep_);
+        return std::tie(na_, nb_, multiplicity_, twice_ms_, irrep_, core_guess_) <
+               std::tie(rhs.na_, rhs.nb_, rhs.multiplicity_, rhs.twice_ms_, rhs.irrep_, rhs.core_guess_);
     } else if (gas_max_ == rhs.gas_max_) {
         // The state with a smaller gas occupation in the first gas space is 'bigger'.
         // Ground state is smaller than core-excited state under this definition.
@@ -83,15 +85,15 @@ bool StateInfo::operator<(const StateInfo& rhs) const {
 }
 
 bool StateInfo::operator!=(const StateInfo& rhs) const {
-    return std::tie(na_, nb_, multiplicity_, twice_ms_, irrep_, gas_min_, gas_max_) !=
+    return std::tie(na_, nb_, multiplicity_, twice_ms_, irrep_, gas_min_, gas_max_, core_guess_) !=
            std::tie(rhs.na_, rhs.nb_, rhs.multiplicity_, rhs.twice_ms_, rhs.irrep_, rhs.gas_min_,
-                    rhs.gas_max_);
+                    rhs.gas_max_, rhs.core_guess_);
 }
 
 bool StateInfo::operator==(const StateInfo& rhs) const {
-    return std::tie(na_, nb_, multiplicity_, twice_ms_, irrep_, gas_min_, gas_max_) ==
+    return std::tie(na_, nb_, multiplicity_, twice_ms_, irrep_, gas_min_, gas_max_, core_guess_) ==
            std::tie(rhs.na_, rhs.nb_, rhs.multiplicity_, rhs.twice_ms_, rhs.irrep_, rhs.gas_min_,
-                    rhs.gas_max_);
+                    rhs.gas_max_, rhs.core_guess_);
 }
 
 StateInfo make_state_info_from_psi(std::shared_ptr<ForteOptions> options) {
@@ -194,6 +196,13 @@ std::size_t StateInfo::hash() const {
         repr += "_" + std::to_string(i);
     for (size_t i : gas_max_)
         repr += "_" + std::to_string(i);
+
+    if (core_guess_) {
+        repr += "_" + std::to_string(1);
+    } else {
+        repr += "_" + std::to_string(0);
+    }
+
     return std::hash<std::string>{}(repr);
 }
 

forte/base_classes/state_info.h

@@ -45,7 +45,7 @@ class StateInfo {
     StateInfo(int na, int nb, int multiplicity, int twice_ms, int irrep,
               const std::string& irrep_label = "",
               const std::vector<size_t> gas_min = std::vector<size_t>(),
-              const std::vector<size_t> gas_max = std::vector<size_t>());
+              const std::vector<size_t> gas_max = std::vector<size_t>(), bool core_guess = false);
 
     StateInfo() = default;
 
@@ -62,6 +62,8 @@ class StateInfo {
     int twice_ms() const;
     /// return the irrep
     int irrep() const;
+    /// return the bool value for running Davidson-Liu with core determinants
+    bool core_guess() const;
     /// return the multiplicity symbol
     const std::string& multiplicity_label() const;
     /// return the irrep symbol
@@ -98,6 +100,8 @@ class StateInfo {
     int twice_ms_;
     // Irrep
     int irrep_;
+    /// Use core determinants as initial guess for Davidson-Liu?
+    bool core_guess_;
     // Irrep label
     std::string irrep_label_;
     // minimum number of electrons in each gas space

forte/fci/fci_solver_initial_guess.cc

@@ -63,14 +63,26 @@ std::vector<Determinant> FCISolver::initial_guess_generate_dets(std::shared_ptr<
                 vec_e_I.begin(), vec_e_I.end(),
                 [&e](const std::tuple<double, size_t>& t) { return e < std::get<0>(t); });
             vec_e_I.insert(it, std::make_tuple(e, I));
-            emax = std::get<0>(vec_e_I.back());
+            // Do not update the maximum energy threshold if using core determinants as initial guess
+            if (!(core_guess_)){
+                emax = std::get<0>(vec_e_I.back());
+            }
             added++;
         }
     }
 
     std::vector<Determinant> guess_dets;
     for (const auto& [e, I] : vec_e_I) {
-        guess_dets.push_back(lists_->determinant(I, symmetry_));
+        const auto& det = lists_->determinant(I, symmetry_);
+        // If using core determinants as initial guess include those with
+        // single and double holes in first bit position
+        if (core_guess_){
+            if (!(det.get_alfa_bit(0) and det.get_beta_bit(0))){
+                guess_dets.push_back(det);
+            }
+        } else {
+            guess_dets.push_back(det);
+        }
     }
 
     // Make sure that the spin space is complete
@@ -99,7 +111,8 @@ FCISolver::initial_guess_det(std::shared_ptr<psi::Vector> diag, size_t num_guess
     // here we use a standard guess procedure
     return find_initial_guess_det(guess_dets, guess_dets_pos, num_guess_states, fci_ints,
                                   state().multiplicity(), true, print_ >= PrintLevel::Default,
-                                  std::vector<std::vector<std::pair<size_t, double>>>());
+                                  std::vector<std::vector<std::pair<size_t, double>>>(),
+                                  core_guess_);
 }
 
 sparse_mat FCISolver::initial_guess_csf(std::shared_ptr<psi::Vector> diag,

forte/genci/genci_solver_initial_guess.cc

@@ -63,14 +63,26 @@ std::vector<Determinant> GenCISolver::initial_guess_generate_dets(std::shared_pt
                 vec_e_I.begin(), vec_e_I.end(),
                 [&e](const std::tuple<double, size_t>& t) { return e < std::get<0>(t); });
             vec_e_I.insert(it, std::make_tuple(e, I));
-            emax = std::get<0>(vec_e_I.back());
+            // Do not update the maximum energy threshold if using core determinants as initial guess
+            if (!(core_guess_)){
+                emax = std::get<0>(vec_e_I.back());
+            }
             added++;
         }
     }
 
     std::vector<Determinant> guess_dets;
     for (const auto& [e, I] : vec_e_I) {
-        guess_dets.push_back(lists_->determinant(I));
+        const auto& det = lists_->determinant(I);
+        // If using core determinants as initial guess include those with
+        // single and double holes in first bit position
+        if (core_guess_){
+            if (!(det.get_alfa_bit(0) and det.get_beta_bit(0))){
+                guess_dets.push_back(det);
+            }
+        } else {
+            guess_dets.push_back(det);
+        }
     }
 
     // Make sure that the spin space is complete
@@ -99,7 +111,8 @@ GenCISolver::initial_guess_det(std::shared_ptr<psi::Vector> diag, size_t num_gue
     // here we use a standard guess procedure
     return find_initial_guess_det(guess_dets, guess_dets_pos, num_guess_states, fci_ints,
                                   state().multiplicity(), true, print_ >= PrintLevel::Default,
-                                  std::vector<std::vector<std::pair<size_t, double>>>());
+                                  std::vector<std::vector<std::pair<size_t, double>>>(),
+                                  core_guess_);
 }
 
 sparse_mat GenCISolver::initial_guess_csf(std::shared_ptr<psi::Vector> diag,

forte/register_forte_options.py

@@ -560,6 +560,8 @@ def register_aci_options(options):
 
     options.add_str("DIAG_ALGORITHM", "SPARSE", ["DYNAMIC", "FULL", "SPARSE"], "The diagonalization method")
 
+    options.add_int_list("CORE_GUESS", "Add core determinants as initial guess")
+
     options.add_bool("FORCE_DIAG_METHOD", False, "Force the diagonalization procedure?")
 
     options.add_bool("ONE_CYCLE", False, "Doing only one cycle of ACI (FCI) ACI iteration?")

forte/sparse_ci/sparse_initial_guess.cc

@@ -74,12 +74,17 @@ find_initial_guess_det(const std::vector<Determinant>& guess_dets,
                        const std::vector<size_t>& guess_dets_pos, size_t num_guess_states,
                        const std::shared_ptr<ActiveSpaceIntegrals>& as_ints, int multiplicity,
                        bool do_spin_project, bool print,
-                       const std::vector<std::vector<std::pair<size_t, double>>>& user_guess) {
+                       const std::vector<std::vector<std::pair<size_t, double>>>& user_guess, 
+                       bool core_guess) {
     size_t num_guess_dets = guess_dets.size();
 
     if (print) {
         print_h2("Initial Guess");
-        psi::outfile->Printf("\n  Initial guess determinants:         %zu", guess_dets.size());
+        std::string guess = "guess";
+        if (core_guess) {
+            guess = "(core) guess";
+        }
+        psi::outfile->Printf("\n  Initial %s determinants:         %zu", guess.c_str(), guess_dets.size());
     }
 
     auto [HS2full, S2evals, S2evecs] =

forte/sparse_ci/sparse_initial_guess.h

@@ -66,7 +66,8 @@ find_initial_guess_det(const std::vector<Determinant>& guess_dets,
                        const std::vector<size_t>& guess_dets_pos, size_t num_guess_states,
                        const std::shared_ptr<ActiveSpaceIntegrals>& as_ints, int multiplicity,
                        bool do_spin_project, bool print,
-                       const std::vector<std::vector<std::pair<size_t, double>>>& user_guess);
+                       const std::vector<std::vector<std::pair<size_t, double>>>& user_guess,
+                       bool core_guess = false);
 
 /// @brief Generate initial guess vectors for the Davidson-Liu solver starting from a set of guess
 /// configurations

tests/methods/fci-core-1/input.dat

@@ -0,0 +1,41 @@
+# LiH 6-31g basis FCI core excited roots
+import forte
+
+refscf = -7.9791777935853290
+reffci = -5.851395993559
+reffci_avg = -5.5360452551438
+
+molecule {
+0 1
+Li
+H 1 R
+
+R = 3.0
+units bohr 
+}
+
+set {
+  basis 6-31g
+  scf_type pk
+  e_convergence 12
+}
+
+set forte {
+  active_space_solver genci
+  core_guess          1
+}
+
+energy('scf')
+compare_values(refscf, variable("CURRENT ENERGY"),11, "SCF energy") #TEST
+
+energy('forte')
+compare_values(reffci, variable("CURRENT ENERGY"),11, "FCI energy") #TEST
+
+set forte {
+  active_space_solver genci
+  avg_state           [[1,1,3]]
+  core_guess          [1]
+}
+
+energy('forte')
+compare_values(reffci_avg, variable("CURRENT ENERGY"),11, "FCI avg energy") #TEST