Add initial orbital guess for DMRG CI/SCF

forte/api/orbital_api.cc

@@ -47,24 +47,21 @@ namespace forte {
 void export_OrbitalTransform(py::module& m) {
     py::class_<OrbitalTransform>(m, "OrbitalTransform")
         .def("compute_transformation", &OrbitalTransform::compute_transformation)
+        .def("set_print", &OrbitalTransform::set_print, "Set printing level")
         .def("get_Ua", &OrbitalTransform::get_Ua, "Get Ua rotation")
         .def("get_Ub", &OrbitalTransform::get_Ub, "Get Ub rotation");
 }
 
 void export_Localize(py::module& m) {
-    py::class_<Localize>(m, "Localize")
+    py::class_<Localize, OrbitalTransform>(m, "Localize")
         .def(py::init<std::shared_ptr<ForteOptions>, std::shared_ptr<ForteIntegrals>,
                       std::shared_ptr<MOSpaceInfo>>())
-        .def("compute_transformation", &Localize::compute_transformation,
-             "Compute the transformation")
         .def("set_orbital_space",
              (void(Localize::*)(std::vector<int>&)) & Localize::set_orbital_space,
              "Compute the transformation")
         .def("set_orbital_space",
              (void(Localize::*)(std::vector<std::string>&)) & Localize::set_orbital_space,
-             "Compute the transformation")
-        .def("get_Ua", &Localize::get_Ua, "Get Ua rotation")
-        .def("get_Ub", &Localize::get_Ub, "Get Ub rotation");
+             "Compute the transformation");
 }
 
 void export_SemiCanonical(py::module& m) {

forte/base_classes/active_space_method.h

@@ -228,7 +228,7 @@ class ActiveSpaceMethod {
         throw std::runtime_error(
             "ActiveSpaceMethod::eigenvectors(): Not Implemented for this class!");
     }
-    
+
     /// Dump the wave function to file
     /// @param file name
     virtual void dump_wave_function(const std::string&) {
@@ -292,6 +292,9 @@ class ActiveSpaceMethod {
     /// @param value the maximum number of iterations
     void set_maxiter(size_t value);
 
+    /// Set if throw an error when Davidson-Liu not converged
+    void set_die_if_not_converged(bool value) { die_if_not_converged_ = value; }
+
     /// Set if we dump the wave function to disk
     void set_read_wfn_guess(bool read);
 
@@ -356,6 +359,9 @@ class ActiveSpaceMethod {
     /// The maximum number of iterations
     size_t maxiter_ = 100;
 
+    /// Stop if Davidson-Liu not converged
+    bool die_if_not_converged_ = true;
+
     /// The root used to compute properties (zero based, default = 0)
     int root_ = 0;
 

forte/base_classes/active_space_solver.cc

@@ -133,7 +133,6 @@ const std::map<StateInfo, std::vector<double>>& ActiveSpaceSolver::compute_energ
                 it->second = make_active_space_method(solver_type_, state, nroot, scf_info_,
                                                       mo_space_info_, as_ints_, options_);
             }
-            auto method = it->second;
         } else {
             state_method_map_[state] = make_active_space_method(
                 solver_type_, state, nroot, scf_info_, mo_space_info_, as_ints_, options_);
@@ -144,6 +143,7 @@ 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_);
+        method->set_die_if_not_converged(die_if_not_converged_);
 
         state_filename_map_[state] = method->wfn_filename();
         if (read_initial_guess_) {
@@ -247,7 +247,8 @@ void ActiveSpaceSolver::compute_multipole_moment(std::shared_ptr<ActiveMultipole
         return;
 
     // compute RDMs
-    psi::outfile->Printf("\n  Computing RDMs for %s moments ...", title.c_str());
+    if (print_ > PrintLevel::Brief)
+        psi::outfile->Printf("\n  Computing RDMs for %s moments ...", title.c_str());
     std::map<StateInfo, std::vector<std::shared_ptr<RDMs>>> state_rdms_map;
     for (const auto& state_nroots : state_nroots_map_) {
         const auto& [state, nroots] = state_nroots;
@@ -261,82 +262,112 @@ void ActiveSpaceSolver::compute_multipole_moment(std::shared_ptr<ActiveMultipole
         state_rdms_map[state] = method->rdms(root_list, max_rdm_level, RDMsType::spin_free);
         method->set_print(print_);
     }
-    psi::outfile->Printf(" Done");
+    if (print_ > PrintLevel::Brief)
+        psi::outfile->Printf(" Done");
 
     // compute dipole moments
     auto dipole_nuc = ampints->nuclear_dipole();
-    std::string prefix = ampints->dp_name().empty() ? "" : ampints->dp_name() + " ";
-    print_h2("Summary of " + prefix + "Dipole Moments [e a0] (Nuclear + Electronic)");
-
-    psi::outfile->Printf("\n    %8s %14s %14s %14s %14s", "State", "DM_X", "DM_Y", "DM_Z", "|DM|");
-    std::string dash(68, '-');
-    psi::outfile->Printf("\n    %s", dash.c_str());
-
+    std::map<StateInfo, std::vector<std::shared_ptr<psi::Vector>>> tmp_for_print;
     for (const auto& state_nroot : state_nroots_map_) {
         const auto& [state, nroots] = state_nroot;
         auto irrep_label = state.irrep_label();
         auto multi_label = upper_string(state.multiplicity_label());
-
         for (size_t i = 0; i < nroots; ++i) {
             std::string name = std::to_string(i) + upper_string(irrep_label);
             auto dipole = ampints->compute_electronic_dipole(state_rdms_map[state][i]);
             dipole->add(*dipole_nuc);
-
+            tmp_for_print[state].push_back(dipole);
             auto dx = dipole->get(0);
             auto dy = dipole->get(1);
             auto dz = dipole->get(2);
             auto dm = dipole->norm();
-            psi::outfile->Printf("\n    %8s%15.8f%15.8f%15.8f%15.8f", name.c_str(), dx, dy, dz, dm);
-
             push_to_psi4_env_globals(dx, multi_label + " <" + name + "|DM_X|" + name + ">");
             push_to_psi4_env_globals(dy, multi_label + " <" + name + "|DM_Y|" + name + ">");
             push_to_psi4_env_globals(dz, multi_label + " <" + name + "|DM_Z|" + name + ">");
             push_to_psi4_env_globals(dm, multi_label + " |<" + name + "|DM|" + name + ">|");
         }
+    }
+
+    if (print_ > PrintLevel::Quiet) {
+        std::string prefix = ampints->dp_name().empty() ? "" : ampints->dp_name() + " ";
+        std::string dash(68, '-');
+        print_h2("Summary of " + prefix + "Dipole Moments [e a0] (Nuclear + Electronic)");
+        psi::outfile->Printf("\n    %8s %14s %14s %14s %14s", "State", "DM_X", "DM_Y", "DM_Z",
+                             "|DM|");
+        psi::outfile->Printf("\n    %s", dash.c_str());
+
+        for (const auto& state_nroot : state_nroots_map_) {
+            const auto& [state, nroots] = state_nroot;
+            auto irrep_label = state.irrep_label();
+            for (size_t i = 0; i < nroots; ++i) {
+                std::string name = std::to_string(i) + upper_string(irrep_label);
+                auto dipole = tmp_for_print[state][i];
+                auto dx = dipole->get(0);
+                auto dy = dipole->get(1);
+                auto dz = dipole->get(2);
+                auto dm = dipole->norm();
+                psi::outfile->Printf("\n    %8s%15.8f%15.8f%15.8f%15.8f", name.c_str(), dx, dy, dz,
+                                     dm);
+            }
+            psi::outfile->Printf("\n    %s", dash.c_str());
+        }
+        psi::outfile->Printf("\n    %8s%15.8f%15.8f%15.8f%15.8f", "Nuclear", dipole_nuc->get(0),
+                             dipole_nuc->get(1), dipole_nuc->get(2), dipole_nuc->norm());
         psi::outfile->Printf("\n    %s", dash.c_str());
     }
-    psi::outfile->Printf("\n    %8s%15.8f%15.8f%15.8f%15.8f", "Nuclear", dipole_nuc->get(0),
-                         dipole_nuc->get(1), dipole_nuc->get(2), dipole_nuc->norm());
-    psi::outfile->Printf("\n    %s", dash.c_str());
 
     // compute quadrupole moments
     if (level != 1) {
+        tmp_for_print.clear();
         auto quadrupole_nuc = ampints->nuclear_quadrupole();
-        std::string prefix = ampints->qp_name().empty() ? "" : ampints->qp_name() + " ";
-        print_h2("Summary of " + prefix + "Quadrupole Moments [e a0^2] (Nuclear + Electronic)");
-
         std::vector<std::string> qm_dirs{"QM_XX", "QM_XY", "QM_XZ", "QM_YY", "QM_YZ", "QM_ZZ"};
-        psi::outfile->Printf("\n    %8s", "State");
-        for (int z = 0; z < 6; ++z)
-            psi::outfile->Printf(" %14s", qm_dirs[z].c_str());
-        std::string dash(98, '-');
-        psi::outfile->Printf("\n    %s", dash.c_str());
 
         for (const auto& state_nroot : state_nroots_map_) {
             const auto& [state, nroots] = state_nroot;
             auto irrep_label = state.irrep_label();
             auto multi_label = upper_string(state.multiplicity_label());
-
             for (size_t i = 0; i < nroots; ++i) {
                 auto quadrupole = ampints->compute_electronic_quadrupole(state_rdms_map[state][i]);
                 quadrupole->add(*quadrupole_nuc);
-
+                tmp_for_print[state].push_back(quadrupole);
                 std::string name = std::to_string(i) + upper_string(irrep_label);
-                psi::outfile->Printf("\n    %8s", name.c_str());
-
                 for (int z = 0; z < 6; ++z) {
                     auto value = quadrupole->get(z);
-                    psi::outfile->Printf("%15.8f", value);
                     std::string dir = "|" + qm_dirs[z] + "|";
                     push_to_psi4_env_globals(value, multi_label + " <" + name + dir + name + ">");
                 }
             }
+        }
+
+        if (print_ > PrintLevel::Quiet) {
+            std::string prefix = ampints->qp_name().empty() ? "" : ampints->qp_name() + " ";
+            std::string dash(98, '-');
+            print_h2("Summary of " + prefix + "Quadrupole Moments [e a0^2] (Nuclear + Electronic)");
+            psi::outfile->Printf("\n    %8s", "State");
+            for (int z = 0; z < 6; ++z)
+                psi::outfile->Printf(" %14s", qm_dirs[z].c_str());
+            psi::outfile->Printf("\n    %s", dash.c_str());
+
+            for (const auto& state_nroot : state_nroots_map_) {
+                const auto& [state, nroots] = state_nroot;
+                auto irrep_label = state.irrep_label();
+                auto multi_label = upper_string(state.multiplicity_label());
+                for (size_t i = 0; i < nroots; ++i) {
+                    auto quadrupole = tmp_for_print[state][i];
+                    std::string name = std::to_string(i) + upper_string(irrep_label);
+                    psi::outfile->Printf("\n    %8s", name.c_str());
+                    for (int z = 0; z < 6; ++z) {
+                        auto value = quadrupole->get(z);
+                        psi::outfile->Printf("%15.8f", value);
+                    }
+                }
+                psi::outfile->Printf("\n    %s", dash.c_str());
+            }
+            psi::outfile->Printf("\n    %8s", "Nuclear");
+            for (int z = 0; z < 6; ++z)
+                psi::outfile->Printf("%15.8f", quadrupole_nuc->get(z));
             psi::outfile->Printf("\n    %s", dash.c_str());
         }
-        psi::outfile->Printf("\n    %8s", "Nuclear");
-        for (int z = 0; z < 6; ++z)
-            psi::outfile->Printf("%15.8f", quadrupole_nuc->get(z));
-        psi::outfile->Printf("\n    %s", dash.c_str());
     }
 }
 
@@ -384,7 +415,8 @@ void ActiveSpaceSolver::compute_fosc_same_orbs(std::shared_ptr<ActiveMultipoleIn
             root_lists_map[{state1, state1}] = state_ids;
 
         for (const auto& [state2, nroot2] : state_nroots_map_) {
-            if (state1 == state2 or std::find(_states.begin(), _states.end(), state2) != _states.end())
+            if (state1 == state2 or
+                std::find(_states.begin(), _states.end(), state2) != _states.end())
                 continue;
             // skip for different multiplicity (no spin-orbit coupling)
             if (state1.multiplicity() != state2.multiplicity())

forte/base_classes/active_space_solver.h

@@ -190,6 +190,9 @@ class ActiveSpaceSolver {
     /// Set the maximum number of iterations
     void set_maxiter(int maxiter);
 
+    /// Set if throw an error when Davidson-Liu not converged
+    void set_die_if_not_converged(bool value) { die_if_not_converged_ = value; }
+
     /// Set if read wave function from file as initial guess
     void set_read_initial_guess(bool read_guess) { read_initial_guess_ = read_guess; }
 
@@ -259,6 +262,9 @@ class ActiveSpaceSolver {
     /// The maximum number of iterations
     size_t maxiter_ = 100;
 
+    /// Stop if Davidson-Liu not converged
+    bool die_if_not_converged_ = true;
+
     /// Read wave function from disk as initial guess
     bool read_initial_guess_;
 

forte/base_classes/orbital_transform.cc

@@ -56,9 +56,7 @@ make_orbital_transformation(const std::string& type, std::shared_ptr<SCFInfo> sc
 
     if (type == "LOCAL") {
         orb_t = std::make_unique<Localize>(options, ints, mo_space_info);
-    } else if (type == "CHOLESKY_ACTIVE") {
-        orb_t = std::make_unique<CholeskyLocal>(options, ints, mo_space_info);
-    } else if (type == "MP2NO") {
+    }else if (type == "MP2NO") {
         orb_t = std::make_unique<MP2_NOS>(scf_info, options, ints, mo_space_info);
     } else if (type == "CINO") {
         orb_t = std::make_unique<CINO>(options, ints, mo_space_info);

forte/base_classes/orbital_transform.h

@@ -1,4 +1,6 @@
-#pragma once 
+#pragma once
+
+#include "helpers/printing.h"
 
 namespace psi {
 class Matrix;
@@ -30,6 +32,8 @@ class OrbitalTransform {
 
     std::shared_ptr<psi::Matrix> get_Ub() { return Ub_; };
 
+    void set_print(PrintLevel level) { print_ = level; }
+
   protected:
     // The integrals
     std::shared_ptr<ForteIntegrals> ints_;
@@ -40,6 +44,9 @@ class OrbitalTransform {
     std::shared_ptr<psi::Matrix> Ua_;
     /// @brief Unitary matrix for beta orbital rotations
     std::shared_ptr<psi::Matrix> Ub_;
+
+    /// Printing level
+    PrintLevel print_ = PrintLevel::Brief;
 };
 
 std::unique_ptr<OrbitalTransform>

forte/fci/fci_solver.cc

@@ -64,8 +64,6 @@ FCISolver::FCISolver(StateInfo state, size_t nroot, std::shared_ptr<MOSpaceInfo>
     nb_ = state.nb() - mo_space_info->size("INACTIVE_DOCC");
 }
 
-void FCISolver::set_maxiter_davidson(int value) { maxiter_davidson_ = 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; }
@@ -156,6 +154,7 @@ void FCISolver::startup() {
 }
 
 void FCISolver::set_options(std::shared_ptr<ForteOptions> options) {
+    set_maxiter(options->get_int("DL_MAXITER"));
     set_e_convergence(options->get_double("E_CONVERGENCE"));
     set_r_convergence(options->get_double("R_CONVERGENCE"));
     set_spin_adapt(options->get_bool("CI_SPIN_ADAPT"));
@@ -168,7 +167,6 @@ void FCISolver::set_options(std::shared_ptr<ForteOptions> options) {
     set_ndets_per_guess_state(options->get_int("DL_DETS_PER_GUESS"));
     set_collapse_per_root(options->get_int("DL_COLLAPSE_PER_ROOT"));
     set_subspace_per_root(options->get_int("DL_SUBSPACE_PER_ROOT"));
-    set_maxiter_davidson(options->get_int("DL_MAXITER"));
 
     set_print(int_to_print_level(options->get_int("PRINT")));
 }
@@ -213,7 +211,7 @@ double FCISolver::compute_energy() {
     dl_solver_->set_e_convergence(e_convergence_);
     dl_solver_->set_r_convergence(r_convergence_);
     dl_solver_->set_print_level(print_);
-    dl_solver_->set_maxiter(maxiter_davidson_);
+    dl_solver_->set_maxiter(maxiter_);
 
     // determine the number of guess vectors
     const size_t num_guess_states = std::min(guess_per_root_ * nroot_, basis_size);
@@ -268,7 +266,7 @@ double FCISolver::compute_energy() {
     dl_solver_->add_sigma_builder(sigma_builder);
 
     auto converged = dl_solver_->solve();
-    if (not converged) {
+    if (not converged and die_if_not_converged_) {
         throw std::runtime_error(
             "Davidson-Liu solver did not converge.\nPlease try to increase the number of "
             "Davidson-Liu iterations (DL_MAXITER). You can also try to increase:\n - the "

forte/fci/fci_solver.h

@@ -87,9 +87,6 @@ class FCISolver : public ActiveSpaceMethod {
     /// Set the number of determinants per root to use to form the initial guess
     void set_ndets_per_guess_state(size_t value);
 
-    /// Set the maximum number of DL iterations
-    void set_maxiter_davidson(int value);
-
     /// Set the number of guess vectors to use
     void set_guess_per_root(int value);
 
@@ -173,8 +170,6 @@ class FCISolver : public ActiveSpaceMethod {
     size_t subspace_per_root_ = 4;
     /// The number of determinants selected for each guess vector
     size_t ndets_per_guess_ = 10;
-    /// Iterations for FCI
-    int maxiter_davidson_ = 30;
     /// Test the RDMs?
     bool test_rdms_ = false;
     /// Print the NO from the 1-RDM

forte/helpers/davidson_liu_solver.cc

@@ -311,7 +311,7 @@ void DavidsonLiuSolver::setup_guesses() {
         }
         basis_size_ += added;
     } else if (basis_size_ >= nroot_) {
-        // psi::outfile->Printf("\n\n  Davidson-Liu solver: restarting from previous calculation.");
+        psi::outfile->Printf("\n\n  Davidson-Liu solver: restarting from previous calculation.");
         basis_size_ = nroot_; // first nroot_ vectors from the previous calculation
         sigma_size_ = 0;      // trigger computation of all sigma vectors
     } else {

forte/integrals/psi4_integrals.cc

@@ -251,7 +251,6 @@ void Psi4Integrals::make_psi4_JK() {
         throw psi::PSIEXCEPTION("Unknown Psi4 integral type to initialize JK in Forte");
     }
 
-    JK_->set_cutoff(schwarz_cutoff_);
     jk_initialize();
     JK_->print_header();
 }
@@ -261,6 +260,7 @@ void Psi4Integrals::jk_initialize(double mem_percentage, int print_level) {
         throw std::runtime_error("Invalid mem_percentage: must be 0 < value < 1.");
     }
     JK_->set_print(print_level);
+    JK_->set_cutoff(schwarz_cutoff_);
     JK_->set_memory(psi::Process::environment.get_memory() * mem_percentage / sizeof(double));
     JK_->initialize();
     JK_status_ = JKStatus::initialized;