Added option to ask for multiple eigenstates in selected_ci_diag. Cal…

…ls LOBPCG with random guess vectors

include/macis/asci/grow.hpp

@@ -128,7 +128,7 @@ auto asci_grow(ASCISettings asci_settings, MCSCFSettings mcscf_settings,
       selected_ci_diag(
           wfn.begin(), wfn.end(), ham_gen, mcscf_settings.ci_matel_tol,
           mcscf_settings.ci_max_subspace, mcscf_settings.ci_res_tol,
-          X_local MACIS_MPI_CODE(, comm));
+          X_local MACIS_MPI_CODE(, comm), false, mcscf_settings.ci_nstates );
 
       if(world_size > 1) {
 #ifdef MACIS_ENABLE_MPI

include/macis/asci/iteration.hpp

@@ -34,7 +34,7 @@ auto asci_iter(ASCISettings asci_settings, MCSCFSettings mcscf_settings,
   auto E = selected_ci_diag<N, index_t>(
       wfn.begin(), wfn.end(), ham_gen, mcscf_settings.ci_matel_tol,
       mcscf_settings.ci_max_subspace, mcscf_settings.ci_res_tol,
-      X_local MACIS_MPI_CODE(, comm));
+      X_local MACIS_MPI_CODE(, comm), false, mcscf_settings.ci_nstates);
 
 #ifdef MACIS_ENABLE_MPI
   auto world_size = comm_size(comm);

include/macis/solvers/lobpcg.hpp

@@ -0,0 +1,114 @@
+/*
+ * MACIS Copyright (c) 2023, The Regents of the University of California,
+ * through Lawrence Berkeley National Laboratory (subject to receipt of
+ * any required approvals from the U.S. Dept. of Energy). All rights reserved.
+ *
+ * See LICENSE.txt for details
+ */
+
+/**
+ * @file lobpcg.h++
+ *
+ * @brief Implements simple call to LOBPCG to get the lowest few
+ *        eigenstates of a given matrix.
+ *
+ * @author Carlos Mejuto Zaera
+ * @date 25/10/2024
+ */
+#pragma once
+#include <assert.h>
+#include <map>
+#include <complex>
+#include <iomanip>
+#include <limits>
+#include <fstream>
+#include <utility>
+#include <sys/stat.h>
+#include <lobpcgxx/lobpcg.hpp>
+#include <sparsexx/matrix_types/csr_matrix.hpp>
+#include <sparsexx/spblas/spmbv.hpp>
+#include <sparsexx/spblas/pspmbv.hpp>
+
+namespace macis {
+
+/**
+ * @brief Perform a band Lanczos calculation on the Hamiltonian operator H, starting from vectors qs, for at most nLanIts
+ *        iterations. Returns the first len(qs) eigenvalues, converged to some accuracy. Note that this implementation
+ *        does not account for deflations (i.e., pruning the span of the qs for linear dependencies in higher powers of H).
+ *
+ * @param[in] typename Functor &H: Hamiltonian oprator. Just needs to implement a matrix vector product.
+ * @param[in] std::vector<std::vector<double> > &qs: Initial set of vetors to perform the band Lanczos on. Deleted on exit.
+ * @param[out] std::vector<double> &evals: Lowest len(qs) eigenvalues.
+ * @param[out] std::vector<std::vector<double> > &evecs: Lowest len(qs) eigenvectors, in the Krylov basis.
+ * @param[in] int &nLanIts: Number of Lanczos iterations to perform.
+ * @param[in] double tol: Target tolerance for the eigenvalue convergence.
+ * @param[in] double thres: Threshold determining when to ignore beta's for being too small.
+ * @param[in] bool print: If true, write intermediate results to file.
+ *
+ * @author Carlos Mejuto Zaera
+ * @date 25/10/2024
+ */ 
+template <typename Functor>
+void LobpcgGS(
+  size_t dimH,
+  size_t nstates,
+  const Functor& H, 
+  std::vector<double> &evals,
+  std::vector<double> &X,
+  int maxIts, 
+  double tol  = 1.E-8,
+  bool print = false )
+{
+  // Run LOBPCG to get the first few eigenstates of a given
+  // Hamiltonian
+  evals.clear(); evals.resize( nstates, 0. );
+  X.clear(); X.resize( dimH * nstates, 0. );
+  // Hop
+  lobpcgxx::operator_action_type<double> Hop = 
+    [&]( int64_t n , int64_t k , const double* x , int64_t ldx ,
+         double* y , int64_t ldy ) -> void {
+
+      for( int ik = 0; ik < k; ik++)
+        H.operator_action(1, 1., x+ik*n, n, 0., y+ik*n, n );
+    };
+
+  // Random vectors 
+  std::default_random_engine gen;
+  std::normal_distribution<> dist(0., 1.);
+  auto rand_gen = [&](){ return dist(gen); };
+  std::generate( X.begin(), X.end(), rand_gen );
+  lobpcgxx::cholqr( dimH, nstates, X.data(), dimH ); // Orthogonalize
+
+  lobpcgxx::lobpcg_settings settings;
+  settings.conv_tol = tol;
+  settings.maxiter  = maxIts;
+  lobpcgxx::lobpcg_operator<double> lob_op( Hop );
+
+
+
+  std::vector<double> res(nstates);
+
+  try
+  {
+    lobpcgxx::lobpcg( settings, dimH , nstates , nstates, 
+                      lob_op, evals.data(), X.data(), dimH, 
+      		  res.data() ); 
+  }
+  catch( std::runtime_error e )
+  {
+    std::cout << "Runtime error during lobpcg: " << e.what() << std::endl;
+  }
+
+  if( print )
+  {
+    std::cout << std::scientific << std::setprecision(10) << std::endl ;
+    for ( int64_t i = 0; i < nstates; ++i ) {
+      std::cout << "  evals[" << i << "] = " << evals[i]
+                << ",   res[" << i << "] = " << res[i]
+                << std::endl ;
+    }
+  }
+
+}
+
+}// namespace macis

include/macis/solvers/selected_ci_diag.hpp

@@ -11,6 +11,7 @@
 #include <macis/csr_hamiltonian.hpp>
 #include <macis/hamiltonian_generator.hpp>
 #include <macis/solvers/davidson.hpp>
+#include <macis/solvers/lobpcg.hpp>
 #include <macis/types.hpp>
 #include <macis/util/mpi.hpp>
 #include <sparsexx/matrix_types/dense_conversions.hpp>
@@ -22,7 +23,8 @@ namespace macis {
 template <typename SpMatType>
 double parallel_selected_ci_diag(const SpMatType& H, size_t davidson_max_m,
                                  double davidson_res_tol,
-                                 std::vector<double>& C_local, MPI_Comm comm) {
+                                 std::vector<double>& C_local, MPI_Comm comm,
+				 const size_t nstates = 1) {
   auto logger = spdlog::get("ci_solver");
   if(!logger) {
     logger = spdlog::stdout_color_mt("ci_solver");
@@ -56,10 +58,28 @@ double parallel_selected_ci_diag(const SpMatType& H, size_t davidson_max_m,
   // Solve EVP
   MPI_Barrier(comm);
   auto dav_st = clock_type::now();
-
-  auto [niter, E] =
-      p_davidson(H.local_row_extent(), davidson_max_m, op, D_local.data(),
-                 davidson_res_tol, C_local.data() MACIS_MPI_CODE(, H.comm()));
+  double E = 0.;
+  size_t niter = 0;
+
+  if( nstates == 1 )
+  {
+    auto [niter0, E0] =
+        p_davidson(H.local_row_extent(), davidson_max_m, op, D_local.data(),
+                   davidson_res_tol, C_local.data() MACIS_MPI_CODE(, H.comm()));
+    niter = niter0; E = E0;
+  }
+  else
+  {
+    std::vector<double> evals;
+    std::vector<double> evecs;
+    LobpcgGS( H.m(), nstates, op, evals, evecs, 10000, davidson_res_tol );
+    logger->info("  Hamiltonian Eigenvalues:");
+    for( int ii = 0; ii < nstates; ii++ )
+      logger->info("    --> Eval #{:4}: {:.6e}", ii, evals[ii]);
+    E = evals[0];
+    for( int ii = 0; ii < H.m(); ii++ )
+      C_local[ii] = evecs[ii];
+  }
 
   MPI_Barrier(comm);
   auto dav_en = clock_type::now();
@@ -75,7 +95,8 @@ double parallel_selected_ci_diag(const SpMatType& H, size_t davidson_max_m,
 template <typename SpMatType>
 double serial_selected_ci_diag(const SpMatType& H, size_t davidson_max_m,
                                double davidson_res_tol,
-                               std::vector<double>& C) {
+                               std::vector<double>& C,
+			       const size_t nstates = 1) {
   auto logger = spdlog::get("ci_solver");
   if(!logger) {
     logger = spdlog::stdout_color_mt("ci_solver");
@@ -108,9 +129,27 @@ double serial_selected_ci_diag(const SpMatType& H, size_t davidson_max_m,
 
   // Solve EVP
   auto dav_st = clock_type::now();
-
-  auto [niter, E] =
-      davidson(H.m(), davidson_max_m, op, D.data(), davidson_res_tol, C.data());
+  double E = 0.;
+  size_t niter = 0;
+
+  if( nstates == 1 )
+  {
+    auto [niter0, E0] =
+        davidson(H.m(), davidson_max_m, op, D.data(), davidson_res_tol, C.data());
+    niter = niter0; E = E0;
+  }
+  else
+  {
+    std::vector<double> evals;
+    std::vector<double> evecs;
+    LobpcgGS( H.m(), nstates, op, evals, evecs, 10000, davidson_res_tol );
+    logger->info("  Hamiltonian Eigenvalues:");
+    for( int ii = 0; ii < nstates; ii++ )
+      logger->info("    --> Eval #{:4}: {:.6e}", ii, evals[ii]);
+    E = evals[0];
+    for( int ii = 0; ii < H.m(); ii++ )
+      C[ii] = evecs[ii];
+  }
 
   auto dav_en = clock_type::now();
 
@@ -128,7 +167,8 @@ double selected_ci_diag(wavefunction_iterator_t<N> dets_begin,
                         size_t davidson_max_m, double davidson_res_tol,
                         std::vector<double>& C_local,
                         MACIS_MPI_CODE(MPI_Comm comm, )
-                            const bool quiet = false) {
+                            const bool quiet = false,
+			    const size_t nstates = 1) {
   auto logger = spdlog::get("ci_solver");
   if(!logger) {
     logger = spdlog::stdout_color_mt("ci_solver");
@@ -197,10 +237,10 @@ double selected_ci_diag(wavefunction_iterator_t<N> dets_begin,
   // Solve EVP
 #ifdef MACIS_ENABLE_MPI
   auto E = parallel_selected_ci_diag(H, davidson_max_m, davidson_res_tol,
-                                     C_local, comm);
+                                     C_local, comm, nstates);
 #else
   auto E =
-      serial_selected_ci_diag(H, davidson_max_m, davidson_res_tol, C_local);
+      serial_selected_ci_diag(H, davidson_max_m, davidson_res_tol, C_local, nstates);
 #endif
 
   return E;

include/macis/util/cas.hpp

@@ -53,7 +53,7 @@ double compute_casci_rdms(
   double E0 =
       selected_ci_diag(dets.begin(), dets.end(), ham_gen, settings.ci_matel_tol,
                        settings.ci_max_subspace, settings.ci_res_tol, C,
-                       MACIS_MPI_CODE(comm, ) true);
+                       MACIS_MPI_CODE(comm, ) true, settings.ci_nstates);
 
   // Compute RDMs
   ham_gen.form_rdms(dets.begin(), dets.end(), dets.begin(), dets.end(),

include/macis/util/mcscf.hpp

@@ -27,6 +27,8 @@ struct MCSCFSettings {
   double ci_res_tol = 1e-8;
   size_t ci_max_subspace = 20;
   double ci_matel_tol = std::numeric_limits<double>::epsilon();
+
+  size_t ci_nstates = 1;
 };
 
 double casscf_diis(MCSCFSettings settings, NumElectron nalpha,

tests/sd_build.cxx

@@ -104,7 +104,7 @@ TEST_CASE("Single Double Build") {
     double E0 = macis::selected_ci_diag(
         dets.begin(), dets.end(), ham_gen, mcscf_settings.ci_matel_tol,
         mcscf_settings.ci_max_subspace, mcscf_settings.ci_res_tol, C_local,
-        MACIS_MPI_CODE(MPI_COMM_WORLD, ) true);
+        MACIS_MPI_CODE(MPI_COMM_WORLD, ) true, mcscf_settings.ci_nstates);
     E0 += E_inactive + E_core;
     // Compute RDMs
     ham_gen.form_rdms(

tests/standalone_driver.cxx

@@ -167,6 +167,8 @@ int main(int argc, char** argv) {
     OPT_KEYWORD("MCSCF.CI_MAX_SUB", mcscf_settings.ci_max_subspace, size_t);
     OPT_KEYWORD("MCSCF.CI_MATEL_TOL", mcscf_settings.ci_matel_tol, double);
 
+    OPT_KEYWORD("MCSCF.CI_NSTATES", mcscf_settings.ci_nstates, size_t);
+
     // ASCI Settings
     macis::ASCISettings asci_settings;
     std::string asci_wfn_fname, asci_wfn_out_fname;