Fix #409 : Managing Psi Dependencies and Options (#410)

* Fix test that shouldn't have run * Remove state info environment.molecule dependency * Remove molecule name environment.molecule dependency * Remove nuclear_dipole environment.molecule dependency --------- Co-authored-by: fevangelista <[email protected]>
evangelistalab · Sep 30, 2024 · 4e6d017 · 4e6d017
1 parent 8efa148
commit 4e6d017
Show file tree

Hide file tree

Showing 15 changed files with 61 additions and 56 deletions.
diff --git a/forte/api/forte_python_module.cc b/forte/api/forte_python_module.cc
@@ -116,8 +116,8 @@ PYBIND11_MODULE(_forte, m) {
 
     m.def("make_orbital_transformation", &make_orbital_transformation,
           "Make an orbital transformation");
-    m.def("make_state_info_from_psi", &make_state_info_from_psi,
-          "Make a state info object from a psi4 Wavefunction");
+    m.def("make_state_info_from_options", &make_state_info_from_options,
+          "Make a state info object from ForteOptions");
     m.def("to_state_nroots_map", &to_state_nroots_map,
           "Convert a map of StateInfo to weight lists to a map of StateInfo to number of "
           "states.");

diff --git a/forte/base_classes/active_space_solver.cc b/forte/base_classes/active_space_solver.cc
@@ -35,9 +35,7 @@
 
 #include "psi4/psi4-dec.h"
 #include "psi4/physconst.h"
-#include "psi4/libmints/molecule.h"
 #include "psi4/libmints/vector.h"
-#include "psi4/libpsi4util/process.h"
 
 #include "base_classes/forte_options.h"
 #include "base_classes/rdms.h"
@@ -663,9 +661,9 @@ make_state_weights_map(std::shared_ptr<ForteOptions> options,
                        std::shared_ptr<MOSpaceInfo> mo_space_info) {
     std::map<StateInfo, std::vector<double>> state_weights_map;
 
-    // make a StateInfo object using the information from psi4
+    // make a StateInfo object using the information
     // TODO: need to optimize for spin-free RDMs
-    auto state = make_state_info_from_psi(options); // assumes low-spin
+    auto state = make_state_info_from_options(options, mo_space_info->point_group_label()); // assumes low-spin
 
     // check if the user provided a AVG_STATE list
     py::list avg_state = options->get_gen_list("AVG_STATE");

diff --git a/forte/base_classes/state_info.cc b/forte/base_classes/state_info.cc
@@ -26,9 +26,6 @@
  * @END LICENSE
  */
 
-#include "psi4/libpsi4util/process.h"
-#include "psi4/libmints/molecule.h"
-
 #include "helpers/helpers.h"
 
 #include "base_classes/forte_options.h"
@@ -94,25 +91,10 @@ bool StateInfo::operator==(const StateInfo& rhs) const {
                     rhs.gas_max_);
 }
 
-StateInfo make_state_info_from_psi(std::shared_ptr<ForteOptions> options) {
-    int charge = psi::Process::environment.molecule()->molecular_charge();
-    if (not options->is_none("CHARGE")) {
-        charge = options->get_int("CHARGE");
-    }
+StateInfo make_state_info_from_options(std::shared_ptr<ForteOptions> options, const Symmetry& symmetry) {
+    int nel = options->get_int("NEL");
 
-    int nel = 0;
-    int natom = psi::Process::environment.molecule()->natom();
-    for (int i = 0; i < natom; i++) {
-        nel += static_cast<int>(psi::Process::environment.molecule()->Z(i));
-    }
-    // If the charge has changed, recompute the number of electrons
-    // Or if you cannot find the number of electrons
-    nel -= charge;
-
-    size_t multiplicity = psi::Process::environment.molecule()->multiplicity();
-    if (not options->is_none("MULTIPLICITY")) {
-        multiplicity = options->get_int("MULTIPLICITY");
-    }
+    size_t multiplicity = options->get_int("MULTIPLICITY");
 
     // If the user did not specify ms determine the value from the input or
     // take the lowest value consistent with the value of "MULTIPLICITY"
@@ -126,7 +108,7 @@ StateInfo make_state_info_from_psi(std::shared_ptr<ForteOptions> options) {
     }
 
     if (((nel - twice_ms) % 2) != 0) {
-        throw std::runtime_error("\n\n  make_state_info_from_psi: Wrong value of M_s.\n\n");
+        throw std::runtime_error("\n\n  make_state_info_from_options: Wrong value of M_s.\n\n");
     }
 
     size_t na = (nel + twice_ms) / 2;
@@ -138,7 +120,7 @@ StateInfo make_state_info_from_psi(std::shared_ptr<ForteOptions> options) {
         irrep = options->get_int("ROOT_SYM");
     }
 
-    std::string irrep_label = psi::Process::environment.molecule()->irrep_labels()[irrep];
+    std::string irrep_label = symmetry.irrep_label(irrep);
     return StateInfo(na, nb, multiplicity, twice_ms, irrep, irrep_label);
 }
 

diff --git a/forte/base_classes/state_info.h b/forte/base_classes/state_info.h
@@ -32,6 +32,8 @@
 #include <string>
 #include <vector>
 
+#include "helpers/symmetry.h"
+
 namespace psi {
 class Wavefunction;
 }
@@ -107,10 +109,9 @@ class StateInfo {
 };
 
 /**
- * @brief make_state_info_from_psi Make a StateInfo object by reading variables set in the psi4
- *        environmental variables
+ * @brief make_state_info_from_options Make a StateInfo object by reading ForteOptions
  * @return a StateInfo object
  */
-StateInfo make_state_info_from_psi(std::shared_ptr<ForteOptions> options);
+StateInfo make_state_info_from_options(std::shared_ptr<ForteOptions> options, const Symmetry& symmetry);
 
 } // namespace forte
diff --git a/forte/dmrg/block2_dmrg_solver.cc b/forte/dmrg/block2_dmrg_solver.cc
@@ -233,7 +233,8 @@ Block2DMRGSolver::Block2DMRGSolver(StateInfo state, size_t nroot, std::shared_pt
         options_->get_bool("BLOCK2_SPIN_ADAPTED");
     std::string scratch = psi::PSIOManager::shared_object()->get_default_path() + "forte." +
                           std::to_string(getpid()) + ".block2." +
-                          psi::Process::environment.molecule()->name();
+                          std::to_string(mo_space_info_->size("ACTIVE")) + "." +
+                          state.str_short();
     Block2ScratchManager::manager().scratch_folders.insert(scratch);
     size_t stack_mem =
         static_cast<size_t>(options_->get_double("BLOCK2_STACK_MEM") * 1024 * 1024 * 1024);

diff --git a/forte/integrals/one_body_integrals.h b/forte/integrals/one_body_integrals.h
@@ -45,7 +45,9 @@ namespace forte {
 class MultipoleIntegrals {
   public:
     /**
-     * @brief Construct a new Multipole Integrals object
+     * @brief Construct a new Multipole Integrals object.
+     *        Warning! This is only valid for a ForteIntegrals
+     *        with wfn_ set. It will segfault otherwise.
      *
      * @param ints forte integral object
      * @param mo_space_info The MOSpaceInfo object

diff --git a/forte/modules/objects_factory_psi4.py b/forte/modules/objects_factory_psi4.py
@@ -225,6 +225,25 @@ def prepare_forte_objects(data, name, **kwargs):
 
     psi4.core.print_out("\n\n  Preparing forte objects from a Psi4 Wavefunction object")
     ref_wfn, mo_space_info = prepare_psi4_ref_wfn(options, **kwargs)
+
+    # Copy state information from the wfn where applicable.
+    # This **MUST** be done before the next function prepares states.
+    # Note that we do not require Psi4 and Forte agree. This allows, e.g., using singlet orbitals for a triplet.
+    molecule = ref_wfn.molecule()
+    if (data.options.is_none("CHARGE")):
+        data.options.set_int("CHARGE", molecule.molecular_charge())
+    elif molecule.molecular_charge() != data.options.get_int("CHARGE"):
+        warnings.warn(f"Psi4 and Forte options disagree about the molecular charge ({molecule.molecular_charge()} vs {data.options.get_int('CHARGE')}).", UserWarning)
+    if (data.options.is_none("MULTIPLICITY")):
+        data.options.set_int("MULTIPLICITY", molecule.multiplicity())
+    elif molecule.multiplicity() != data.options.get_int("MULTIPLICITY"):
+        warnings.warn(f"Psi4 and Forte options disagree about the multiplicity ({molecule.multiplicity()}) vs ({data.options.get_int('MULTIPLICITY')}).", UserWarning)
+    nel = int(sum(molecule.Z(i) for i in range(molecule.natom()))) - data.options.get_int("CHARGE")
+    if (data.options.is_none("NEL")):
+        data.options.set_int("NEL", nel)
+    elif nel != data.options.get_int("NEL"):
+        warnings.warn(f"Psi4 and Forte options disagree about the number of electorns ({nel}) vs ({data.options.get_int('NEL')}).", UserWarning)
+
     forte_objects = prepare_forte_objects_from_psi4_wfn(options, ref_wfn, mo_space_info)
     state_weights_map, mo_space_info, scf_info = forte_objects
     fcidump = None

diff --git a/forte/modules/options_factory.py b/forte/modules/options_factory.py
@@ -24,10 +24,10 @@ def __init__(self, options: dict = None):
     def _run(self, data: ForteData = None) -> ForteData:
         if data is None:
             data = ForteData()
-        data.options = forte.forte_options
         # if no options dict is provided then read from psi4
         if self.options is None:
-            # Get the option object
+            # Copy globals into a new object
+            data.options = forte.ForteOptions(forte.forte_options)
             psi4_options = psi4.core.get_options()
             psi4_options.set_current_module("FORTE")
 

diff --git a/forte/modules/tdaci.py b/forte/modules/tdaci.py
@@ -1,7 +1,7 @@
 from typing import List
 
 from forte.data import ForteData
-from forte._forte import make_state_info_from_psi, to_state_nroots_map, make_active_space_method, TDCI
+from forte._forte import to_state_nroots_map, make_active_space_method, TDCI
 
 from .module import Module
 from .active_space_ints import ActiveSpaceInts
@@ -21,8 +21,12 @@ def __init__(self):
         super().__init__()
 
     def _run(self, data: ForteData) -> ForteData:
-        state = make_state_info_from_psi(data.options)
         data = ActiveSpaceInts(active="ACTIVE", core=["RESTRICTED_DOCC"]).run(data)
+        # As long as TDCI takes ActiveSpaceMethod rather than ActiveSpaceSolver, this error message must stay.
+        if len(data.state_weights_map) != 1:
+            raise("TDACI does not support multi-state computations yet. File an issue if you need this!")
+        else:
+            state = next(iter(data.state_weights_map))
         state_map = to_state_nroots_map(data.state_weights_map)
         active_space_method = make_active_space_method(
             "ACI", state, data.options.get_int("NROOT"), data.scf_info, data.mo_space_info, data.as_ints, data.options

diff --git a/forte/mrdsrg-spin-adapted/sa_mrdsrg.cc b/forte/mrdsrg-spin-adapted/sa_mrdsrg.cc
@@ -29,8 +29,6 @@
 #include <unistd.h>
 #include <algorithm>
 
-#include "psi4/libmints/molecule.h"
-#include "psi4/libpsi4util/process.h"
 #include "psi4/libpsi4util/PsiOutStream.h"
 #include "psi4/libpsio/psio.hpp"
 
@@ -112,7 +110,7 @@ void SA_MRDSRG::startup() {
     // determine file names
     chk_filename_prefix_ = psi::PSIOManager::shared_object()->get_default_path() + "forte." +
                            std::to_string(getpid()) + "." +
-                           psi::Process::environment.molecule()->name();
+                           std::to_string(mo_space_info_->size("ACTIVE"));
     t1_file_chk_.clear();
     t2_file_chk_.clear();
     if (restart_amps_ and (relax_ref_ != "NONE")) {

diff --git a/forte/mrdsrg-spin-adapted/sadsrg.cc b/forte/mrdsrg-spin-adapted/sadsrg.cc
@@ -29,7 +29,6 @@
 #include <numeric>
 
 #include "psi4/libmints/matrix.h"
-#include "psi4/libmints/molecule.h"
 #include "psi4/libpsi4util/process.h"
 #include "psi4/libpsi4util/PsiOutStream.h"
 #include "psi4/libmints/vector.h"
@@ -139,7 +138,7 @@ void SADSRG::startup() {
     // setup checkpoint filename prefix
     chk_filename_prefix_ = PSIOManager::shared_object()->get_default_path();
     chk_filename_prefix_ += "forte." + std::to_string(getpid());
-    chk_filename_prefix_ += "." + psi::Process::environment.molecule()->name();
+    chk_filename_prefix_ += "." + std::to_string(mo_space_info_->size("ACTIVE"));
     Bcan_files_.clear();
 }
 

diff --git a/forte/mrdsrg-spin-integrated/dsrg_mrpt2_grad/dsrg_mrpt2_deriv_write_rdms.cc b/forte/mrdsrg-spin-integrated/dsrg_mrpt2_grad/dsrg_mrpt2_deriv_write_rdms.cc
@@ -12,8 +12,8 @@
 #include "psi4/lib3index/3index.h"
 #include "psi4/libmints/mintshelper.h"
 #include "base_classes/mo_space_info.h"
-#include "psi4/libmints/molecule.h"
 #include "psi4/libmints/vector.h"
+#include "integrals/one_body_integrals.h"
 
 using namespace ambit;
 using namespace psi;
@@ -113,8 +113,8 @@ void DSRG_MRPT2::write_1rdm_spin_dependent() {
     std::map<char, std::vector<size_t>> idxmap_abs;
     idxmap_abs = {{'c', core_all_}, {'a', actv_all_}, {'v', virt_all_}};
     std::vector<double> dipole(4, 0.0);
-    Vector3 dm_nuc =
-        psi::Process::environment.molecule()->nuclear_dipole(psi::Vector3(0.0, 0.0, 0.0));
+    MultipoleIntegrals dm_ints(ints_, mo_space_info_);
+    auto dm_nuc = dm_ints.nuclear_dipole(psi::Vector3(0.0, 0.0, 0.0));
 
     for (int i = 0; i < 3; ++i) {
         auto dm_ints = mo_dipole_ints[i];
@@ -127,7 +127,7 @@ void DSRG_MRPT2::write_1rdm_spin_dependent() {
             });
         }
         dipole[i] = 2.0 * D1_temp["pq"] * dipole_ints["pq"];
-        dipole[i] += dm_nuc[i];
+        dipole[i] += (*dm_nuc)[i];
         dipole[3] += dipole[i] * dipole[i];
     }
     dipole[3] = std::sqrt(dipole[3]);

diff --git a/forte/mrdsrg-spin-integrated/master_mrdsrg.cc b/forte/mrdsrg-spin-integrated/master_mrdsrg.cc
@@ -4,12 +4,12 @@
 #include "psi4/libpsi4util/PsiOutStream.h"
 #include "psi4/libpsi4util/process.h"
 
-#include "psi4/libmints/molecule.h"
 #include "psi4/libmints/matrix.h"
 #include "psi4/libmints/vector.h"
 
 #include "base_classes/mo_space_info.h"
 #include "integrals/active_space_integrals.h"
+#include "integrals/one_body_integrals.h"
 
 #include "helpers/printing.h"
 #include "helpers/timer.h"
@@ -448,10 +448,10 @@ double MASTER_DSRG::compute_reference_energy_df(BlockedTensor H, BlockedTensor F
 
 void MASTER_DSRG::init_dm_ints() {
     outfile->Printf("\n    Preparing ambit tensors for dipole moments ...... ");
-    Vector3 dm_nuc =
-        psi::Process::environment.molecule()->nuclear_dipole(psi::Vector3(0.0, 0.0, 0.0));
+    MultipoleIntegrals dm_ints(ints_, mo_space_info_);
+    auto dm_nuc = dm_ints.nuclear_dipole(psi::Vector3(0.0, 0.0, 0.0));
     for (int i = 0; i < 3; ++i) {
-        dm_nuc_[i] = dm_nuc[i];
+        dm_nuc_[i] = (*dm_nuc)[i];
         dm_[i] = BTF_->build(tensor_type_, "Dipole " + dm_dirs_[i], spin_cases({"gg"}));
     }
 

diff --git a/forte/mrdsrg-spin-integrated/mrdsrg.cc b/forte/mrdsrg-spin-integrated/mrdsrg.cc
@@ -31,10 +31,8 @@
 #include <map>
 #include <vector>
 
-#include "psi4/libmints/molecule.h"
-#include "psi4/libpsi4util/process.h"
 #include "psi4/libpsi4util/PsiOutStream.h"
-#include "psi4/libmints/molecule.h"
+#include "psi4/libpsio/psio.hpp"
 
 #define FMT_HEADER_ONLY
 #include "lib/fmt/core.h"
@@ -132,7 +130,7 @@ void MRDSRG::startup() {
     // set up file name prefix
     restart_file_prefix_ = psi::PSIOManager::shared_object()->get_default_path() + "forte." +
                            std::to_string(getpid()) + "." +
-                           psi::Process::environment.molecule()->name();
+                           std::to_string(mo_space_info_->size("ACTIVE"));
     t1_file_chk_.clear();
     t2_file_chk_.clear();
     if (restart_amps_ and (relax_ref_ != "NONE") and

diff --git a/forte/pymodule.py b/forte/pymodule.py
@@ -279,6 +279,10 @@ def gradient_forte(name, **kwargs):
     # if job_type not in {"CASSCF", "MCSCF_TWO_STEP"} and correlation_solver != "DSRG-MRPT2":
     #     raise Exception("Analytic energy gradients are only implemented for" " CASSCF, MCSCF_TWO_STEP, or DSRG-MRPT2.")
 
+    if "FCIDUMP" in int_type:
+        raise Exception("Analytic gradients with FCIDUMP are not theoretically possible.")
+    if int_type == "PYSCF":
+        raise "Analytic gradients with PySCF are not yet implemented."
     # Prepare Forte objects: state_weights_map, mo_space_info, scf_info
     data = ObjectsFromPsi4(**kwargs).run(data)
 
@@ -353,7 +357,6 @@ def mr_dsrg_pt2(job_type, data):
     scf_info = data.scf_info
     ints = data.ints
 
-    state = forte.make_state_info_from_psi(options)
     # generate a list of states with their own weights
     state_map = forte.to_state_nroots_map(state_weights_map)