Fix issues with apps and mcscf

Author: fevangelista

forte/apps/__init__.py

@@ -1,2 +1,2 @@
-from .mcscf import run_mcscf
-from .hf import run_hf
+from .mcscf import mcscf
+from .hf import hf

forte/apps/hf.py

@@ -6,16 +6,21 @@
 from forte.modules.hf import HF
 
 
-def run_hf(geom, basis, state):
+def hf(geom, basis, state, e_convergence=1e-8, d_convergence=1e-8, options=None):
     charge, multiplicity, sym = parse_state(state)
     hf_workflow = Workflow(
         [
-            OptionsFactory(),
+            # Pass the options to the ForteOptions object
+            OptionsFactory(options),
+            # Generate the molecule
             MoleculeFactory(geom),
+            # Generate the state
             StateFactory(charge=charge, multiplicity=multiplicity, sym=sym),
-            HF(basis=basis),
+            # Run the HF calculation
+            HF(basis=basis, e_convergence=e_convergence, d_convergence=d_convergence),
         ],
         name="HF Workflow",
     )
     data = hf_workflow.run()
+    data.options.set_str("basis", basis)
     return data

forte/apps/mcscf.py

@@ -1,13 +1,20 @@
 from .apps_helpers import parse_active_space
-from .hf import run_hf
+from .hf import hf
 from forte.modules.workflow import Workflow
 from forte.modules.objects_factory_psi4 import ObjectsFromPsi4
 from forte.modules.active_space_selector import ActiveSpaceSelector
 from forte.modules.mcscf import MCSCF
+from forte.modules.options_factory import OptionsFactory
 
 
-def run_mcscf(geom, basis, state, active_space, solver_type="FCI"):
-    data = run_hf(geom, basis, state)
+def mcscf(active_space, solver_type="FCI", options=None, data=None, geom=None, basis=None, state=None):
+    if data is None and geom is not None and basis is not None:
+        # If data is None, create a new data object
+        data = hf(geom, basis, state, options=options)
+    elif data is None:
+        raise ValueError("Either data or geom and basis must be provided.")
+
+    # Setup a workflow to run MCSCF
     mcscf_workflow = Workflow(
         [
             # Parse the active space before defining the MOSpaceInfo object
@@ -19,5 +26,8 @@ def run_mcscf(geom, basis, state, active_space, solver_type="FCI"):
         ],
         name="MCSCF Workflow",
     )
+
+    # Run the workflow
     data = mcscf_workflow.run(data)
+
     return data

forte/mcscf/mcscf_2step.cc

@@ -425,16 +425,7 @@ double MCSCF_2STEP::compute_energy() {
         diis_manager.delete_diis_file();
     } // end of MCSCF macro iterations
 
-    // perform final active space computation using converged orbitals
-    if (print_ >= PrintLevel::Default)
-        psi::outfile->Printf("\n\n  Performing final CI Calculation using converged orbitals");
-
-    as_solver_->set_maxiter(as_maxiter);
-    as_solver_->set_die_if_not_converged(true);
-    energy_ =
-        diagonalize_hamiltonian(cas_grad.active_space_ints(),
-                                {print_, e_conv_, r_conv, options_->get_bool("DUMP_ACTIVE_WFN")});
-
+    // Transform the orbitals to canonical or natural basis
     if (ints_->integral_type() != Custom) {
         // fix orbitals for redundant pairs
         rdms = as_solver_->compute_average_rdms(state_weights_map_, 1, RDMsType::spin_free);
@@ -456,6 +447,9 @@ double MCSCF_2STEP::compute_energy() {
         semi.semicanonicalize(rdms, false, actv_orb_type, false);
         cas_grad.canonicalize_final(semi.Ua());
 
+        // after this step, the integrals, CI coefficients, and rdms are no longer valid.
+        // Below we rediagonalize the Hamiltonian to get the new CI coefficients.
+
         // pass the MO coefficients to the wave function
         auto Ca = cas_grad.Ca();
         ints_->wfn()->Ca()->copy(Ca);
@@ -518,6 +512,17 @@ double MCSCF_2STEP::compute_energy() {
             throw_convergence_error();
     }
 
+    // Perform final active space computation using converged transformed orbitals
+    if (der_type_ != "FIRST" or options_->get_str("CORRELATION_SOLVER") != "NONE") {
+        if (print_ >= PrintLevel::Default)
+            psi::outfile->Printf("\n\n  Performing final CI Calculation using converged orbitals");
+        as_solver_->set_maxiter(as_maxiter);
+        as_solver_->set_die_if_not_converged(true);
+        energy_ = diagonalize_hamiltonian(
+            cas_grad.active_space_ints(),
+            {print_, e_conv_, r_conv, options_->get_bool("DUMP_ACTIVE_WFN")});
+    }
+
     return energy_;
 }
 

forte/modules/active_space_selector.py

@@ -7,7 +7,13 @@
 
 class ActiveSpaceSelector(Module):
     """
-    A module to prepare an ActiveSpaceIntegral
+    A module to prepare the active space for the MCSCF calculation.
+    This code will use the active space information from the input file to select the active space.
+
+    The active space can be specified in several ways:
+    - active_space = {"nel": 6, "norb": 6}  # 6 electrons in 6 orbitals
+    - active_space = {"nel": 6, "active_orbitals": ["1 A1", "2 A1", "3 A1"]}  # 6 electrons in the specified orbitals
+    - active_space = {"active": [2, 0, 1, 1]}  # select 2 x A1, 1 x B1, and 1 x B2 active orbitals
     """
 
     # accept as input a dictionary of string,int list pairs or a list of integers
@@ -19,7 +25,6 @@ def __init__(self, active_space: Union[Dict[str, List[int]], List[int]] = None):
             The active space to be used for the calculation
         """
         super().__init__()
-        # assert isinstance(active_space, dict) or isinstance(active_space, list)
 
         self.active_space = active_space
 
@@ -29,7 +34,7 @@ def _run(self, data: ForteData) -> ForteData:
         Selects the active space based on the input
         """
 
-        # in the case where we have no active space selection, we just assemble this object from the options
+        # if no active space is provided, use the options object
         self.active_space = self.active_space or self.make_dict_from_data(data)
 
         # select the active space using the information in the active_space dictionary

forte/modules/hf.py

@@ -55,23 +55,24 @@ def _run(self, data: ForteData) -> ForteData:
         molecule.set_multiplicity(state.multiplicity())
 
         # # prepare options for psi4
-        # scf_type_dict = {
-        #     "CONVENTIONAL": "PK",
-        #     "STD": "PK",
-        # }
-        # # convert to psi4 terminology
-        # int_type = self.model.int_type.upper()
-        # if int_type in scf_type_dict:
-        #     scf_type = scf_type_dict[int_type]
-        # else:
-        #     scf_type = int_type
-        scf_type = "PK"
+        scf_type_dict = {
+            "CONVENTIONAL": "PK",
+            "STD": "PK",
+        }
+        # convert to psi4 terminology
+        int_type = data.options.get_str("int_type")
+        if int_type in scf_type_dict:
+            scf_type = scf_type_dict[int_type]
+        else:
+            scf_type = int_type
 
         if self.restricted:
             ref = "RHF" if state.multiplicity() == 1 else "ROHF"
         else:
             ref = "UHF"
 
+        print(f"HF: using {scf_type} for the SCF calculation")
+
         options = {
             "BASIS": self.basis,
             "REFERENCE": ref,

forte/modules/objects_factory_psi4.py

@@ -368,5 +368,6 @@ def prepare_forte_objects_from_wavefunction(self, data):
         # Build Forte SCFInfo object
         data.scf_info = SCFInfo(data.psi_wfn)
 
-        # Build a map from Forte StateInfo to the weights
-        data.state_weights_map = make_state_weights_map(data.options, data.mo_space_info)
+        # Build a map of StateInfo and weights from the Options and MOSpaceInfo objects if not already done elsewhere
+        if data.state_weights_map is None:
+            data.state_weights_map = make_state_weights_map(data.options, data.mo_space_info)

tests/pytest/apps/test_hf_app.py

@@ -7,7 +7,7 @@ def test_rhf_app():
     H 0.0 0.0 0.0
     H 0.0 0.0 1.0
     """
-    data = forte.apps.run_hf(geom=xyz, basis="cc-pVDZ", state={"charge": 0, "multiplicity": 1, "sym": "ag"})
+    data = forte.apps.hf(geom=xyz, basis="cc-pVDZ", state={"charge": 0, "multiplicity": 1, "sym": "ag"})
     assert data.results.value("hf energy") == pytest.approx(-1.10015376479352, 1.0e-10)
 
 

tests/pytest/apps/test_mcscf_app.py

@@ -7,7 +7,7 @@ def test_mcscf_app1():
     H 0.0 0.0 0.0
     H 0.0 0.0 1.5
     """
-    data = forte.apps.run_mcscf(
+    data = forte.apps.mcscf(
         geom=xyz,
         basis="cc-pVDZ",
         state={"charge": 0, "multiplicity": 1, "sym": "ag"},
@@ -22,7 +22,7 @@ def test_mcscf_app2():
     H 0.0 0.0 0.0
     H 0.0 0.0 1.5
     """
-    data = forte.apps.run_mcscf(
+    data = forte.apps.mcscf(
         geom=xyz,
         basis="cc-pVDZ",
         state={"charge": 0, "multiplicity": 1, "sym": "ag"},
@@ -37,7 +37,7 @@ def test_mcscf_app3():
     H 0.0 0.0 0.0
     H 0.0 0.0 1.5
     """
-    data = forte.apps.run_mcscf(
+    data = forte.apps.mcscf(
         geom=xyz,
         basis="cc-pVDZ",
         state={"charge": 0, "multiplicity": 1, "sym": "ag"},
@@ -48,6 +48,22 @@ def test_mcscf_app3():
     data.scf_info.reorder_orbitals([[0, 1, 2], [], [0], [0], [], [1, 0, 2], [0], [0]])
 
 
+# def test_mcscf_app4():
+#     xyz = """
+#     H 0.0 0.0 0.0
+#     H 0.0 0.0 1.5
+#     """
+#     data = forte.apps.run_mcscf(
+#         geom=xyz,
+#         basis="cc-pVDZ",
+#         states=[{"charge": 0, "multiplicity": 1, "sym": "ag", "weights": [0.5, 0.5]}],
+#         active_space={"active_orbitals": ["1 Ag", "1 B1u"], "nel": 2},
+#         solver_type="GENCI",
+#     )
+#     assert data.results.value("mcscf energy") == pytest.approx(-1.0561253825629822, 1.0e-10)
+#     data.scf_info.reorder_orbitals([[0, 1, 2], [], [0], [0], [], [1, 0, 2], [0], [0]])
+
+
 if __name__ == "__main__":
     test_mcscf_app1()
     test_mcscf_app2()