Spin adaptation of DETCI (#336)

* Theory meeting, adding spin adaptation to sparce_ci_solver

* Enable SpinAdapter

* Ksa (#340)

* SA for sparce_ci

* Debugging ncmo

* fixed sa for sparceCI

* Remove swp files

* Fix sa-fci-1 input

* Remove ref output files

* Add test cases detci-5,6,7

* Remove unnecessary files

* Fix bug and remove s2_labels

* Fix test case

* Delete setup.cfg

* Improvements to the DL class and remove obsolete option

* Fix test cases

* Fix python test case and document function

* Fix a bug and testcases. Remove unnecessary call to dls.get_result()

* Code cleanup

* Loosen residual condition

* Factor code for initial guess

* Introduce the DL_GUESS_PER_ROOT option

* Fix failing test cases

* Make sure all codes push options

* Pretty printing

* Rename folder

* Use same code for detci and fci CSF guess

* Minor fixes

* Add missing files

* Squash warnings and implement suggested change

* Add interface to new functions

---------

Co-authored-by: marink2 <kmarin3@emory.edu>
Co-authored-by: Kevin Marin <45517485+marink2@users.noreply.github.com>

docs/notebooks/methods/nb_01_fci.ipynb

@@ -201,6 +201,52 @@
     "```"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "b78a6658",
+   "metadata": {},
+   "source": [
+    "    \n",
+    "## Initial guess and parameters of the Davidson-Liu algorithm\n",
+    "\n",
+    "The Davidson-Liu procedure requires an initial set of guess vectors. In Forte, the initial guesses are determined by diagonalizing the Hamiltonian in a small space of determinants $\\{ |\\Phi'_I\\rangle \\}$. This space is spin complete. The procedure starts with the diagonalization of the matrix representation of the $\\hat{S}^2$ operator\n",
+    "$$\n",
+    "(\\mathbf{S}^2)_{IJ} = \\langle \\Phi'_I | \\hat{S}^2 |\\Phi'_J\\rangle\n",
+    "$$\n",
+    "which allows to group the solution according to the multiplicity.\n",
+    "After this step, we transform the Hamiltonian to the basis of eigenstates of $\\hat{S}^2$ and diagonalize each block separately.\n",
+    "The guess vectors are taken from the lowest energy solutions with the correct multiplicity.\n",
+    "If a guess vector with wrong multiplicity happens to fall below the energy of guess vector with the correct multiplicity, we project the incorrect guess during the diagonalization procedure.\n",
+    "\n",
+    "The number of guess vectors and the size of the determinant space used to generate the initial guess is determined by options controlled by the user.\n",
+    "Once the user specifies the number of target states (roots) and the multiplicity, then the option `DL_GUESS_PER_ROOT` (default = 1) controls the number of guess states requested, that is:\n",
+    "\n",
+    "Number of guess states =  Number of roots $\\times$  `DL_GUESS_PER_ROOT`\n",
+    "\n",
+    "The subspace of determinants selected for the initial guess has size controlled by the option `DL_DETS_PER_GUESS` (default = 50) , which determines the number of determinants included in the following way\n",
+    "\n",
+    "Number of determinants used to build guesses = Number of guess states $\\times$ `DL_DETS_PER_GUESS`\n",
+    "\n",
+    "Note that the number of determinants used in the initial guess procedure could be larger than the one determined by the equation above. This happens if the determinants selected do not form a spin-complete set. In this case, additional determinants are added to ensure that the determinant space is spin-complete.\n",
+    "\n",
+    "Two other options control the Davidson-Liu procedure:\n",
+    "- `DL_SUBSPACE_PER_ROOT`: this option controls the maximum number of subspace vectors stored by the DL algorithm:\n",
+    "\n",
+    "    Maximum number of subspace vectors =  Number of roots $\\times$ `DL_SUBSPACE_PER_ROOT`\n",
+    "\n",
+    "- `DL_COLLAPSE_PER_ROOT`: once the maximum number of subspace vectors is reached, the DL procedure forms the best solution vectors and collapses (resets) the subspace size.\n",
+    "This option controls the number of vectors kept after collapes:\n",
+    "\n",
+    "    Number of subspace vectors after collapse =  Number of roots $\\times$ `DL_COLLAPSE_PER_ROOT`\n",
+    "\n",
+    "**When should you modify these options?**\n",
+    "`DL_GUESS_PER_ROOT` and `DL_DETS_PER_GUESS` should be modified **if the DL procedure has trouble finding the correct initial guess**. For example, if the determinant space is too small, the initial guess procedure may produce a list of states with incorrect energetic ordering. See the test case `forte/tests/methods/detci-5` for an example.\n",
+    "\n",
+    "`DL_COLLAPSE_PER_ROOT` and `DL_SUBSPACE_PER_ROOT` should be modified **if the DL procedure has trouble converging**.\n",
+    "\n",
+    "Note that these options need to be changed carefully. If you ask for more guess vectors than those available in a given space the code will fail. Additionally, `DL_COLLAPSE_PER_ROOT` should be greater or equal to `DL_GUESS_PER_ROOT`, and `DL_SUBSPACE_PER_ROOT` should be greater than `DL_COLLAPSE_PER_ROOT`.\n"
+   ]
+  },
   {
    "attachments": {},
    "cell_type": "markdown",
@@ -308,6 +354,12 @@
     "    --------------------------------------------------------\n",
     "```\n"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e5a79d8e",
+   "metadata": {},
+   "source": []
   }
  ],
  "metadata": {

docs/source/methods/nb_01_fci.rst

@@ -238,6 +238,83 @@ quadrupole operators and the occupation numer of natural orbitals
            1B3u    0.143579      1B2u    0.143579      2B1u    0.024482
            ...
 
+Initial guess and parameters of the Davidson-Liu algorithm
+----------------------------------------------------------
+
+The Davidson-Liu procedure requires an initial set of guess vectors. In
+Forte, the initial guesses are determined by diagonalizing the
+Hamiltonian in a small space of determinants
+:math:`\{ |\Phi'_I\rangle \}`. This space is spin complete. The
+procedure starts with the diagonalization of the matrix representation
+of the :math:`\hat{S}^2` operator
+
+.. math::
+
+
+   (\mathbf{S}^2)_{IJ} = \langle \Phi'_I | \hat{S}^2 |\Phi'_J\rangle
+
+which allows to group the solution according to the multiplicity. After
+this step, we transform the Hamiltonian to the basis of eigenstates of
+:math:`\hat{S}^2` and diagonalize each block separately. The guess
+vectors are taken from the lowest energy solutions with the correct
+multiplicity. If a guess vector with wrong multiplicity happens to fall
+below the energy of guess vector with the correct multiplicity, we
+project the incorrect guess during the diagonalization procedure.
+
+The number of guess vectors and the size of the determinant space used
+to generate the initial guess is determined by options controlled by the
+user. Once the user specifies the number of target states (roots) and
+the multiplicity, then the option ``DL_GUESS_PER_ROOT`` (default = 1)
+controls the number of guess states requested, that is:
+
+Number of guess states = Number of roots :math:`\times`
+``DL_GUESS_PER_ROOT``
+
+The subspace of determinants selected for the initial guess has size
+controlled by the option ``DL_DETS_PER_GUESS`` (default = 50) , which
+determines the number of determinants included in the following way
+
+Number of determinants used to build guesses = Number of guess states
+:math:`\times` ``DL_DETS_PER_GUESS``
+
+Note that the number of determinants used in the initial guess procedure
+could be larger than the one determined by the equation above. This
+happens if the determinants selected do not form a spin-complete set. In
+this case, additional determinants are added to ensure that the
+determinant space is spin-complete.
+
+Two other options control the Davidson-Liu procedure: -
+``DL_SUBSPACE_PER_ROOT``: this option controls the maximum number of
+subspace vectors stored by the DL algorithm:
+
+::
+
+   Maximum number of subspace vectors =  Number of roots $\times$ `DL_SUBSPACE_PER_ROOT`
+
+-  ``DL_COLLAPSE_PER_ROOT``: once the maximum number of subspace vectors
+   is reached, the DL procedure forms the best solution vectors and
+   collapses (resets) the subspace size. This option controls the number
+   of vectors kept after collapes:
+
+   Number of subspace vectors after collapse = Number of roots
+   :math:`\times` ``DL_COLLAPSE_PER_ROOT``
+
+**When should you modify these options?** ``DL_GUESS_PER_ROOT`` and
+``DL_DETS_PER_GUESS`` should be modified **if the DL procedure has
+trouble finding the correct initial guess**. For example, if the
+determinant space is too small, the initial guess procedure may produce
+a list of states with incorrect energetic ordering. See the test case
+``forte/tests/methods/detci-5`` for an example.
+
+``DL_COLLAPSE_PER_ROOT`` and ``DL_SUBSPACE_PER_ROOT`` should be modified
+**if the DL procedure has trouble converging**.
+
+Note that these options need to be changed carefully. If you ask for
+more guess vectors than those available in a given space the code will
+fail. Additionally, ``DL_COLLAPSE_PER_ROOT`` should be greater or equal
+to ``DL_GUESS_PER_ROOT``, and ``DL_SUBSPACE_PER_ROOT`` should be greater
+than ``DL_COLLAPSE_PER_ROOT``.
+
 Spin-adapted FCI
 ----------------
 
@@ -365,3 +442,5 @@ composition and in the final energy summary that reports the value of
        --------------------------------------------------------
           5  (  0)    Ag     0      -12.596862494551   6.000000
        --------------------------------------------------------
+
+

forte/CMakeLists.txt

@@ -123,6 +123,7 @@ helpers/cube_file.cc
 helpers/disk_io.cc
 helpers/helpers.cc
 helpers/symmetry.cc
+helpers/determinant_helpers.cc
 helpers/iterative_solvers.cc
 helpers/lbfgs/lbfgs.cc
 helpers/lbfgs/lbfgs_param.cc
@@ -223,6 +224,7 @@ sparse_ci/sigma_vector_full.cc
 sparse_ci/sigma_vector_sparse_list.cc
 sparse_ci/sorted_string_list.cc
 sparse_ci/sparse_ci_solver.cc
+sparse_ci/sparse_initial_guess.cc
 sparse_ci/sparse_operator.cc
 sparse_ci/sparse_exp.cc
 sparse_ci/sparse_fact_exp.cc

forte/api/sci_api.cc

@@ -32,6 +32,8 @@
 #include "psi4/libmints/vector.h"
 #include "psi4/libmints/matrix.h"
 
+#include "helpers/determinant_helpers.h"
+
 #include "sparse_ci/sigma_vector.h"
 #include "sparse_ci/sparse_ci_solver.h"
 #include "integrals/active_space_integrals.h"
@@ -336,6 +338,11 @@ void export_Determinant(py::module& m) {
     m.def("get_projection", &get_projection);
     m.def("overlap", &overlap);
     m.def("spin2", &spin2<Determinant::nbits>);
+    m.def("make_hamiltonian_matrix", &make_hamiltonian_matrix, "dets"_a, "as_ints"_a,
+          "Make a Hamiltonian matrix (psi::Matrix) from a list of determinants and an "
+          "ActiveSpaceIntegrals object");
+    m.def("make_s2_matrix", &make_s2_matrix, "dets"_a,
+          "Make a matrix (psi::Matrix) of the S^2 operator from a list of determinants");
 }
 
 void export_SigmaVector(py::module& m) {

forte/attic/fci_mo.cc

@@ -913,7 +913,7 @@ void FCI_MO::Diagonalize_H(const vecdet& p_space, const int& multi, const int& n
     sparse_solver.set_r_convergence(rconv_);
     sparse_solver.set_spin_project(options_->get_bool("SCI_PROJECT_OUT_SPIN_CONTAMINANTS"));
     sparse_solver.set_maxiter_davidson(options_->get_int("DL_MAXITER"));
-    sparse_solver.set_guess_dimension(options_->get_int("DL_GUESS_SIZE"));
+    sparse_solver.set_ndets_per_guess_state(options_->get_int("DL_DETS_PER_GUESS"));
     if (!projected_roots_.empty()) {
         sparse_solver.set_root_project(true);
         sparse_solver.add_bad_states(projected_roots_);

forte/ci_ex_states/excited_state_solver.cc

@@ -88,11 +88,14 @@ void ExcitedStateSolver::set_options(std::shared_ptr<ForteOptions> options) {
     sparse_solver_->set_parallel(true);
     sparse_solver_->set_force_diag(options->get_bool("FORCE_DIAG_METHOD"));
     sparse_solver_->set_e_convergence(options->get_double("E_CONVERGENCE"));
+    sparse_solver_->set_r_convergence(options->get_double("R_CONVERGENCE"));
+    sparse_solver_->set_guess_per_root(options->get_int("DL_GUESS_PER_ROOT"));
+    sparse_solver_->set_ndets_per_guess_state(options->get_int("DL_DETS_PER_GUESS"));
+    sparse_solver_->set_collapse_per_root(options->get_int("DL_COLLAPSE_PER_ROOT"));
+    sparse_solver_->set_subspace_per_root(options->get_int("DL_SUBSPACE_PER_ROOT"));
     sparse_solver_->set_maxiter_davidson(options->get_int("DL_MAXITER"));
     sparse_solver_->set_spin_project(options->get_bool("SCI_PROJECT_OUT_SPIN_CONTAMINANTS"));
     sparse_solver_->set_spin_project_full(options->get_bool("SCI_PROJECT_OUT_SPIN_CONTAMINANTS"));
-    sparse_solver_->set_guess_dimension(options->get_int("DL_GUESS_SIZE"));
-    sparse_solver_->set_num_vecs(options->get_int("N_GUESS_VEC"));
     sci_->set_options(options);
 }
 
@@ -502,8 +505,6 @@ void ExcitedStateSolver::print_final(DeterminantHashVec& dets,
 void ExcitedStateSolver::print_wfn(DeterminantHashVec& space, std::shared_ptr<psi::Matrix> evecs,
                                    int nroot) {
     std::string state_label;
-    std::vector<std::string> s2_labels({"singlet", "doublet", "triplet", "quartet", "quintet",
-                                        "sextet", "septet", "octet", "nonet", "decatet"});
 
     //    std::vector<std::pair<double, double>> spins = compute_spin(space, op, evecs, nroot);
 
@@ -523,7 +524,7 @@ void ExcitedStateSolver::print_wfn(DeterminantHashVec& space, std::shared_ptr<ps
                                  tmp_evecs[I] * tmp_evecs[I], space.get_idx(tmp.get_det(I)),
                                  str(tmp.get_det(I), nact_).c_str());
         }
-        //        state_label = s2_labels[std::round(spins[n].first * 2.0)];
+        //        state_label = s2_label(std::round(spins[n].first * 2.0));
         //        psi::outfile->Printf("\n\n  Spin state for root %zu: S^2 = %5.6f, S = %5.3f,
         //        %s", n,
         //                             spins[n].first, spins[n].second, state_label.c_str());

forte/fci/fci_solver.cc

@@ -65,6 +65,10 @@ FCISolver::FCISolver(StateInfo state, size_t nroot, std::shared_ptr<MOSpaceInfo>
 
 void FCISolver::set_fci_iterations(int value) { fci_iterations_ = value; }
 
+void FCISolver::set_ndets_per_guess_state(size_t value) { ndets_per_guess_ = value; }
+
+void FCISolver::set_guess_per_root(int value) { guess_per_root_ = value; }
+
 void FCISolver::set_collapse_per_root(int value) { collapse_per_root_ = value; }
 
 void FCISolver::set_subspace_per_root(int value) { subspace_per_root_ = value; }
@@ -122,8 +126,10 @@ void FCISolver::set_options(std::shared_ptr<ForteOptions> options) {
     set_root(options->get_int("ROOT"));
     set_test_rdms(options->get_bool("FCI_TEST_RDMS"));
     set_fci_iterations(options->get_int("FCI_MAXITER"));
+    set_guess_per_root(options->get_int("DL_GUESS_PER_ROOT"));
     set_collapse_per_root(options->get_int("DL_COLLAPSE_PER_ROOT"));
     set_subspace_per_root(options->get_int("DL_SUBSPACE_PER_ROOT"));
+    set_ndets_per_guess_state(options->get_int("DL_DETS_PER_GUESS"));
     set_print(options->get_int("PRINT"));
     set_e_convergence(options->get_double("E_CONVERGENCE"));
     set_r_convergence(options->get_double("R_CONVERGENCE"));
@@ -171,18 +177,18 @@ double FCISolver::compute_energy() {
     dls.set_subspace_per_root(subspace_per_root_);
 
     // determine the number of guess vectors
-    const size_t guess_size = std::min(collapse_per_root_ * nroot_, basis_size);
+    const size_t num_guess_states = std::min(guess_per_root_ * nroot_, basis_size);
 
     // Form the diagonal of the Hamiltonian and the initial guess
     if (spin_adapt_) {
         auto Hdiag_vec = form_Hdiag_csf(as_ints_, spin_adapter_);
         dls.startup(Hdiag_vec);
-        initial_guess_csf(Hdiag_vec, guess_size, dls, sigma_basis);
+        initial_guess_csf(Hdiag_vec, num_guess_states, dls);
     } else {
         Hdiag.form_H_diagonal(as_ints_);
         Hdiag.copy_to(sigma);
         dls.startup(sigma);
-        initial_guess_det(Hdiag, guess_size, as_ints_, dls, sigma);
+        initial_guess_det(Hdiag, num_guess_states, as_ints_, dls);
     }
 
     // Set a variable to track the convergence of the solver
@@ -260,9 +266,6 @@ double FCISolver::compute_energy() {
         throw std::runtime_error("FCI did not converge. Try increasing FCI_MAXITER.");
     }
 
-    // Compute final eigenvectors
-    dls.get_results();
-
     // Copy eigenvalues and eigenvectors from the Davidson-Liu solver
     evals_ = dls.eigenvalues();
     energies_ = std::vector<double>(nroot_, 0.0);
@@ -282,7 +285,7 @@ double FCISolver::compute_energy() {
 
     // Print determinants
     if (print_) {
-        print_solutions(guess_size, b, b_basis, dls);
+        print_solutions(num_guess_states, b, b_basis, dls);
     }
 
     // Optionally, test the RDMs