From c216e8ccf3b4812fa2225c629b0e107006fe5fab Mon Sep 17 00:00:00 2001 From: fevangelista Date: Fri, 17 May 2024 14:04:59 -0400 Subject: [PATCH 1/4] Massive changes --- forte/api/forte_python_module.cc | 4 +- forte/api/integrals_api.cc | 4 +- forte/api/orbital_api.cc | 5 +- forte/base_classes/README.md | 1 + forte/base_classes/active_space_solver.cc | 4 +- forte/base_classes/orbital_transform.cc | 26 +++--- forte/base_classes/orbital_transform.h | 15 +-- forte/base_classes/orbitals.cc | 46 +++++++-- forte/base_classes/orbitals.h | 23 ++++- forte/data.py | 14 ++- forte/helpers/determinant_helpers.cc | 3 +- forte/helpers/helpers.cc | 90 +++++++++++++++++- forte/helpers/helpers.h | 22 +++++ forte/integrals/cholesky_integrals.cc | 5 +- forte/integrals/cholesky_integrals.h | 9 +- forte/integrals/conventional_integrals.cc | 4 +- forte/integrals/conventional_integrals.h | 2 +- forte/integrals/custom_integrals.cc | 66 ++++++------- forte/integrals/custom_integrals.h | 12 +-- forte/integrals/df_integrals.cc | 6 +- forte/integrals/df_integrals.h | 10 +- forte/integrals/diskdf_integrals.cc | 15 +-- forte/integrals/diskdf_integrals.h | 10 +- forte/integrals/integrals.cc | 65 ++++++------- forte/integrals/integrals.h | 61 ++++++------ forte/integrals/integrals_psi4_interface.cc | 93 ++++++++++--------- forte/integrals/make_integrals.cc | 24 ++--- forte/integrals/make_integrals.h | 19 ++-- forte/mcscf/mcscf_2step.cc | 13 +-- forte/mcscf/mcscf_2step.h | 8 +- forte/mcscf/mcscf_orb_grad.cc | 54 ++++++----- forte/mcscf/mcscf_orb_grad.h | 7 +- forte/mrdsrg-spin-adapted/sa_dsrgpt.cc | 2 +- forte/mrdsrg-spin-adapted/sa_mrpt2.cc | 38 +++++--- forte/mrdsrg-spin-adapted/sa_mrpt2_oeprop.cc | 76 +++++++++------ forte/mrdsrg-spin-adapted/sadsrg.cc | 15 +-- forte/mrdsrg-spin-adapted/sadsrg.h | 4 +- .../dsrg_mrpt2_deriv_write_rdms.cc | 18 ++-- .../three_dsrg_mrpt2.cc | 22 +++-- forte/orbital-helpers/ci-no/ci-no.cc | 6 +- forte/orbital-helpers/ci-no/ci-no.h | 4 +- forte/orbital-helpers/ci-no/mrci-no.cc | 5 +- forte/orbital-helpers/ci-no/mrci-no.h | 3 +- forte/orbital-helpers/localize.cc | 10 +- forte/orbital-helpers/localize.h | 4 +- forte/orbital-helpers/mp2_nos.cc | 5 +- forte/orbital-helpers/mp2_nos.h | 3 +- forte/orbital-helpers/mrpt2_nos.cc | 8 +- forte/orbital-helpers/mrpt2_nos.h | 6 +- forte/orbital-helpers/semi_canonicalize.cc | 2 +- forte/orbital-helpers/unpaired_density.cc | 27 ++++-- forte/orbital-helpers/unpaired_density.h | 24 ++--- forte/post_process/spin_corr.cc | 22 +++-- forte/post_process/spin_corr.h | 4 +- 54 files changed, 658 insertions(+), 390 deletions(-) diff --git a/forte/api/forte_python_module.cc b/forte/api/forte_python_module.cc index 467710484..14519acb9 100644 --- a/forte/api/forte_python_module.cc +++ b/forte/api/forte_python_module.cc @@ -44,6 +44,7 @@ #include "base_classes/dynamic_correlation_solver.h" #include "base_classes/state_info.h" #include "base_classes/scf_info.h" +#include "base_classes/orbitals.h" #include "integrals/make_integrals.h" @@ -191,7 +192,8 @@ PYBIND11_MODULE(_forte, m) { m.def("make_custom_ints", &make_custom_forte_integrals, "Make a custom Forte integral object from arrays"); m.def("make_ints_from_psi4", &make_forte_integrals_from_psi4, "ref_wfn"_a, "options"_a, - "mo_space_info"_a, "int_type"_a = "", "Make a Forte integral object from psi4"); + "mo_space_info"_a, "orbitals"_a, "int_type"_a = "", + "Make a Forte integral object from psi4"); m.def("make_active_space_method", &make_active_space_method, "Make an active space method"); m.def("make_active_space_solver", &make_active_space_solver, "Make an active space solver", "method"_a, "state_nroots_map"_a, "scf_info"_a, "mo_space_info"_a, "options"_a, diff --git a/forte/api/integrals_api.cc b/forte/api/integrals_api.cc index 9d10b75fa..417ea92cc 100644 --- a/forte/api/integrals_api.cc +++ b/forte/api/integrals_api.cc @@ -31,6 +31,7 @@ #include "helpers/helpers.h" #include "integrals/integrals.h" +#include "base_classes/orbitals.h" namespace py = pybind11; @@ -40,8 +41,7 @@ namespace forte { void export_ForteIntegrals(py::module& m) { py::class_>(m, "ForteIntegrals") .def("rotate_orbitals", &ForteIntegrals::rotate_orbitals, "Rotate MOs during contructor") - .def("Ca", &ForteIntegrals::Ca, "Return the alpha MO coefficients") - .def("Cb", &ForteIntegrals::Cb, "Return the beta MO coefficients") + // .def("orbitals", &ForteIntegrals::orbitals, "Return the orbitals object") .def("nmo", &ForteIntegrals::nmo, "Return the total number of moleuclar orbitals") .def("ncmo", &ForteIntegrals::ncmo, "Return the number of correlated orbitals") .def("frozen_core_energy", &ForteIntegrals::frozen_core_energy, diff --git a/forte/api/orbital_api.cc b/forte/api/orbital_api.cc index 6239333b0..a0a891a9b 100644 --- a/forte/api/orbital_api.cc +++ b/forte/api/orbital_api.cc @@ -34,6 +34,7 @@ #include "base_classes/orbital_transform.h" #include "base_classes/forte_options.h" #include "base_classes/mo_space_info.h" +#include "base_classes/orbitals.h" #include "integrals/integrals.h" #include "orbital-helpers/localize.h" @@ -56,8 +57,8 @@ void export_OrbitalTransform(py::module& m) { /// Export the ForteOptions class void export_Localize(py::module& m) { py::class_(m, "Localize") - .def(py::init, std::shared_ptr, - std::shared_ptr>()) + .def(py::init, std::shared_ptr, + std::shared_ptr, std::shared_ptr>()) .def("compute_transformation", &Localize::compute_transformation, "Compute the transformation") .def("set_orbital_space", diff --git a/forte/base_classes/README.md b/forte/base_classes/README.md index 505c8dc52..2e39e28c8 100644 --- a/forte/base_classes/README.md +++ b/forte/base_classes/README.md @@ -5,6 +5,7 @@ Canonical order of base classes - StateInfo - SCFInfo - MOSpaceInfo + - Orbitals - Integrals/ActiveSpaceIntegrals - ForteOptions diff --git a/forte/base_classes/active_space_solver.cc b/forte/base_classes/active_space_solver.cc index b070f94b2..3e31488d5 100644 --- a/forte/base_classes/active_space_solver.cc +++ b/forte/base_classes/active_space_solver.cc @@ -112,7 +112,7 @@ const std::map>& ActiveSpaceSolver::compute_energ } // initialize multipole integrals - if (as_ints_->ints()->integral_type() != Custom) { + if (as_ints_->ints()->integral_type() != IntegralType::Custom) { if (not as_mp_ints_) { auto mp_ints = std::make_shared(as_ints_->ints(), mo_space_info_); as_mp_ints_ = std::make_shared(mp_ints); @@ -160,7 +160,7 @@ const std::map>& ActiveSpaceSolver::compute_energ } print_energies(); - if (as_ints_->ints()->integral_type() != Custom and + if (as_ints_->ints()->integral_type() != IntegralType::Custom and options_->get_str("ACTIVE_SPACE_SOLVER") != "EXTERNAL") { compute_dipole_moment(as_mp_ints_); compute_quadrupole_moment(as_mp_ints_); diff --git a/forte/base_classes/orbital_transform.cc b/forte/base_classes/orbital_transform.cc index 7247669e8..92189ad18 100644 --- a/forte/base_classes/orbital_transform.cc +++ b/forte/base_classes/orbital_transform.cc @@ -42,26 +42,26 @@ namespace forte { -OrbitalTransform::OrbitalTransform(std::shared_ptr ints, - std::shared_ptr mo_space_info) - : ints_(ints), mo_space_info_(mo_space_info) {} +OrbitalTransform::OrbitalTransform(std::shared_ptr mo_space_info, + std::shared_ptr orbitals, + std::shared_ptr ints) + : mo_space_info_(mo_space_info), orbitals_(orbitals), ints_(ints) {} -std::unique_ptr -make_orbital_transformation(const std::string& type, std::shared_ptr scf_info, - std::shared_ptr options, - std::shared_ptr ints, - std::shared_ptr mo_space_info) { +std::unique_ptr make_orbital_transformation( + const std::string& type, std::shared_ptr scf_info, + std::shared_ptr options, std::shared_ptr ints, + std::shared_ptr mo_space_info, std::shared_ptr orbitals) { std::unique_ptr orb_t; if (type == "LOCAL") { - orb_t = std::make_unique(options, ints, mo_space_info); + orb_t = std::make_unique(options, mo_space_info, orbitals, ints); } else if (type == "MP2NO") { - orb_t = std::make_unique(scf_info, options, ints, mo_space_info); + orb_t = std::make_unique(scf_info, options, mo_space_info, orbitals, ints); } else if (type == "CINO") { - orb_t = std::make_unique(options, ints, mo_space_info); + orb_t = std::make_unique(options, mo_space_info, orbitals, ints); } else if (type == "MRCINO") { - orb_t = std::make_unique(scf_info, options, ints, mo_space_info); + orb_t = std::make_unique(scf_info, options, mo_space_info, orbitals, ints); } else if (type == "MRPT2NO") { // perform reference CI calculation auto as_type = options->get_str("ACTIVE_SPACE_SOLVER"); @@ -78,7 +78,7 @@ make_orbital_transformation(const std::string& type, std::shared_ptr sc as_solver->compute_average_rdms(state_weights_map, rdm_level, RDMsType::spin_free); // initialize - orb_t = std::make_unique(rdms, scf_info, options, ints, mo_space_info); + orb_t = std::make_unique(scf_info, options, mo_space_info, orbitals, ints, rdms); } else { throw std::runtime_error("Orbital type " + type + " is not supported!"); } diff --git a/forte/base_classes/orbital_transform.h b/forte/base_classes/orbital_transform.h index e69dc1992..680c3fba5 100644 --- a/forte/base_classes/orbital_transform.h +++ b/forte/base_classes/orbital_transform.h @@ -1,4 +1,4 @@ -#pragma once +#pragma once namespace psi { class Matrix; @@ -7,16 +7,17 @@ class Matrix; namespace forte { class ForteIntegrals; +class ForteOptions; class MOSpaceInfo; class SCFInfo; -class ForteOptions; +class Orbitals; class OrbitalTransform { public: /// Constructor - OrbitalTransform(std::shared_ptr ints, - std::shared_ptr mo_space_info); + OrbitalTransform(std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints = nullptr); /// Default constructor OrbitalTransform() = default; @@ -31,10 +32,12 @@ class OrbitalTransform { std::shared_ptr get_Ub() { return Ub_; }; protected: - // The integrals - std::shared_ptr ints_; /// The MOSpace info std::shared_ptr mo_space_info_; + /// The orbitals + std::shared_ptr orbitals_; + // The integrals + std::shared_ptr ints_; /// @brief Unitary matrix for alpha orbital rotations std::shared_ptr Ua_; diff --git a/forte/base_classes/orbitals.cc b/forte/base_classes/orbitals.cc index 2ea30a015..aa596ed87 100644 --- a/forte/base_classes/orbitals.cc +++ b/forte/base_classes/orbitals.cc @@ -27,19 +27,51 @@ */ #include "psi4/libmints/wavefunction.h" +#include "helpers/helpers.h" #include "orbitals.h" namespace forte { -Orbitals::Orbitals(const std::shared_ptr& wfn, bool restricted) { - // Grab the MO coefficients from psi and enforce spin restriction if necessary - Ca_ = wfn->Ca()->clone(); - Cb_ = restricted ? wfn->Ca()->clone() : wfn->Cb()->clone(); + +Orbitals::Orbitals(const std::shared_ptr& Ca, const std::shared_ptr& Cb) { + Ca_ = Ca->clone(); + Cb_ = Cb->clone(); +} + +const std::shared_ptr Orbitals::Ca() const { return Ca_; } + +const std::shared_ptr Orbitals::Cb() const { return Cb_; } + +void Orbitals::set(const std::shared_ptr& Ca, const std::shared_ptr& Cb) { + if (not(elementwise_compatible_matrices(Ca_, Ca) and + elementwise_compatible_matrices(Cb_, Cb))) { + throw std::runtime_error( + "Orbital::set: orbital coefficient matrices have different dimensions!"); + } + Ca_ = Ca->clone(); + Cb_ = Cb->clone(); +} + +void Orbitals::rotate(std::shared_ptr Ua, std::shared_ptr Ub) { + // 1. Rotate the orbital coefficients and store them in the ForteIntegral object + auto Ca_rotated = psi::linalg::doublet(Ca_, Ua); + auto Cb_rotated = psi::linalg::doublet(Cb_, Ub); + Ca_->copy(Ca_rotated); + Cb_->copy(Cb_rotated); +} + +void Orbitals::copy(const Orbitals& other) { + Ca_->copy(*other.Ca_); + Cb_->copy(*other.Cb_); +} + +bool Orbitals::are_spin_restricted(double threshold) const { + return matrix_distance(Ca_, Cb_) < threshold; } -std::unique_ptr make_orbitals(const std::shared_ptr& wfn, - bool restricted) { - return std::make_unique(wfn, restricted); +std::unique_ptr +make_orbitals_from_psi(const std::shared_ptr& wfn, bool restricted) { + return std::make_unique(wfn->Ca(), restricted ? wfn->Ca() : wfn->Cb()); } } // namespace forte \ No newline at end of file diff --git a/forte/base_classes/orbitals.h b/forte/base_classes/orbitals.h index 7e3d1f10e..671dc6ef4 100644 --- a/forte/base_classes/orbitals.h +++ b/forte/base_classes/orbitals.h @@ -42,13 +42,26 @@ namespace forte { class Orbitals { public: // ==> Class Constructor <== - Orbitals(const std::shared_ptr& wfn, bool restricted); + Orbitals(const std::shared_ptr& Ca, const std::shared_ptr& Cb); // ==> Class Methods <== /// @return The alpha orbital coefficient matrix - const std::shared_ptr Ca() const { return Ca_; } + const std::shared_ptr Ca() const; /// @return The beta orbital coefficient matrix - const std::shared_ptr Cb() const { return Cb_; } + const std::shared_ptr Cb() const; + + void set(const std::shared_ptr& Ca, const std::shared_ptr& Cb); + + /// @brief Rotate the orbitals using the given transformation matrices + /// The orbital coefficients are rotated according to the following formula: + /// C(new) = C(old) U + /// @param Ua Alpha orbital transformation matrix + /// @param Ub Beta orbital transformation matrix + void rotate(std::shared_ptr Ua, std::shared_ptr Ub); + + void copy(const Orbitals& other); + + bool are_spin_restricted(double threshold = 1.0e-8) const; private: // ==> Class Data <== @@ -58,7 +71,7 @@ class Orbitals { std::shared_ptr Cb_; }; -std::unique_ptr make_orbitals(const std::shared_ptr& wfn, - bool restricted); +std::unique_ptr +make_orbitals_from_psi(const std::shared_ptr& wfn, bool restricted); } // namespace forte diff --git a/forte/data.py b/forte/data.py index 888a30b8d..952c6d67d 100644 --- a/forte/data.py +++ b/forte/data.py @@ -3,14 +3,15 @@ from psi4.core import Wavefunction, Molecule from ._forte import ( - ForteOptions, - SCFInfo, - MOSpaceInfo, - ForteIntegrals, ActiveSpaceIntegrals, - StateInfo, ActiveSpaceSolver, + ForteOptions, + ForteIntegrals, + MOSpaceInfo, + Orbitals, RDMs, + SCFInfo, + StateInfo, Symmetry, ) from forte import Model, Results @@ -34,6 +35,8 @@ class ForteData: The model object. scf_info: SCFInfo The SCF information. + orbitals: Orbitals + The molecular orbitals. mo_space_info: MOSpaceInfo The molecular orbital space information. ints: ForteIntegrals @@ -54,6 +57,7 @@ class ForteData: molecule: Molecule = None model: Model = None scf_info: SCFInfo = None + orbitals: Orbitals = None mo_space_info: MOSpaceInfo = None ints: ForteIntegrals = None as_ints: ActiveSpaceIntegrals = None diff --git a/forte/helpers/determinant_helpers.cc b/forte/helpers/determinant_helpers.cc index 540797776..805abd1d2 100644 --- a/forte/helpers/determinant_helpers.cc +++ b/forte/helpers/determinant_helpers.cc @@ -69,7 +69,8 @@ make_hamiltonian_matrix(const std::vector& dets, auto H = std::make_shared("H", n, n); // If we are running DiskDF then we need to revert to a single thread loop - auto threads = (as_ints->get_integral_type() == DiskDF) ? 1 : omp_get_max_threads(); + auto threads = + (as_ints->get_integral_type() == IntegralType::DiskDF) ? 1 : omp_get_max_threads(); #pragma omp parallel for schedule(dynamic) num_threads(threads) for (size_t I = 0; I < n; I++) { diff --git a/forte/helpers/helpers.cc b/forte/helpers/helpers.cc index da78dcde7..88b903884 100644 --- a/forte/helpers/helpers.cc +++ b/forte/helpers/helpers.cc @@ -122,6 +122,93 @@ std::vector Vector_to_vector_double(const psi::Vector& v) { return v_double; } +ambit::Tensor matrix_to_tensor(const std::shared_ptr& m, const std::string& label) { + auto nirrep = m->nirrep(); + auto rowspi = m->rowspi(); + auto colspi = m->colspi(); + auto nrows = static_cast(rowspi.sum()); + auto ncols = static_cast(colspi.sum()); + auto T = ambit::Tensor::build(ambit::CoreTensor, label, {nrows, ncols}); + auto& T_data = T.data(); + + size_t row_offset = 0; + size_t col_offset = 0; + + for (int h = 0; h < nirrep; ++h) { + for (int p = 0; p < rowspi[h]; ++p) { + for (int q = 0; q < colspi[h]; ++q) { + size_t p_full = p + row_offset; + size_t q_full = q + col_offset; + T_data[p_full * ncols + q_full] = m->get(h, p, q); + } + } + row_offset += rowspi[h]; + col_offset += rowspi[h]; + } + return T; +} + +ambit::Tensor matrix_to_tensor(const std::shared_ptr& m, + const std::string& label) { + auto nirrep = m->nirrep(); + auto rowspi = m->rowspi(); + auto colspi = m->colspi(); + auto nrows = static_cast(rowspi.sum()); + auto ncols = static_cast(colspi.sum()); + auto T = ambit::Tensor::build(ambit::CoreTensor, label, {nrows, ncols}); + auto& T_data = T.data(); + + size_t row_offset = 0; + size_t col_offset = 0; + + for (int h = 0; h < nirrep; ++h) { + for (int p = 0; p < rowspi[h]; ++p) { + for (int q = 0; q < colspi[h]; ++q) { + size_t p_full = p + row_offset; + size_t q_full = q + col_offset; + T_data[p_full * ncols + q_full] = m->get(h, p, q); + } + } + row_offset += rowspi[h]; + col_offset += rowspi[h]; + } + return T; +} + +bool elementwise_compatible_matrices(const std::shared_ptr& A, + const std::shared_ptr& B) { + if (A->nirrep() != B->nirrep()) { + return false; + } + for (int h = 0; h < A->nirrep(); ++h) { + if (A->rowspi()[h] != B->rowspi()[h] or A->colspi()[h] != B->colspi()[h]) { + return false; + } + } + return true; +} + +double matrix_distance(const std::shared_ptr& A, const std::shared_ptr& B, + int p) { + if (not elementwise_compatible_matrices(A, B)) { + throw std::runtime_error("Matrix distance: Matrices have different dimensions!"); + } + + auto A_minus_B = A->clone(); + A_minus_B->subtract(B); + double result{0.0}; + // If p is -1, we calculate the infinity norm + if (p == -1) { + result = A_minus_B->absmax(); + } else if (p == 2) { + result = A_minus_B->rms(); + } else { + std::runtime_error( + "Matrix distance: Invalid norm! Only 2-norm and infinity norm are supported."); + } + return result; +} + std::pair to_xb(size_t nele, size_t type_size) { if (nele == 0) return {0.0, "B"}; @@ -258,7 +345,8 @@ std::vector> find_integer_groups(const std::vect return groups; } -// void view_modified_orbitals(psi::SharedWavefunction wfn, const std::shared_ptr& Ca, +// void view_modified_orbitals(psi::SharedWavefunction wfn, const std::shared_ptr& +// Ca, // const std::shared_ptr& diag_F, // const std::shared_ptr& occupation) { // std::shared_ptr molden(new MoldenWriter(wfn)); diff --git a/forte/helpers/helpers.h b/forte/helpers/helpers.h index fd6d00c8c..ca29aef89 100644 --- a/forte/helpers/helpers.h +++ b/forte/helpers/helpers.h @@ -93,6 +93,28 @@ std::shared_ptr tensor_to_matrix(ambit::Tensor t, psi::Dimension di std::vector Vector_to_vector_double(const psi::Vector& v); +/// @brief Convert a psi::Matrix to an ambit::Tensor +/// @param m The input matrix +/// @param label The label of the tensor (optional) +ambit::Tensor matrix_to_tensor(const std::shared_ptr& m, + const std::string& label = ""); +ambit::Tensor matrix_to_tensor(const std::shared_ptr& m, + const std::string& label = ""); + +/// @brief Check if two matrices are compatible for elementwise operations +/// @param A The first matrix +/// @param B The second matrix +/// @return True if the matrices are compatible +bool elementwise_compatible_matrices(const std::shared_ptr& A, + const std::shared_ptr& B); + +/// @brief Compute the rank-p norm of the difference between two matrices A and B +/// @param A The first matrix +/// @param B The second matrix +/// @param rank The rank of the norm (default is 2) +double matrix_distance(const std::shared_ptr& A, const std::shared_ptr& B, + int p = 2); + // /** // * @brief view_modified_orbitals Write orbitals using molden // * @param Ca The Ca matrix to be viewed with MOLDEN diff --git a/forte/integrals/cholesky_integrals.cc b/forte/integrals/cholesky_integrals.cc index c59b6f0db..6d64dd55c 100644 --- a/forte/integrals/cholesky_integrals.cc +++ b/forte/integrals/cholesky_integrals.cc @@ -55,8 +55,9 @@ namespace forte { CholeskyIntegrals::CholeskyIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, IntegralSpinRestriction restricted) - : Psi4Integrals(options, ref_wfn, mo_space_info, Cholesky, restricted) { + : Psi4Integrals(options, ref_wfn, mo_space_info, orbitals, IntegralType::Cholesky, restricted) { initialize(); } @@ -145,7 +146,7 @@ ambit::Tensor CholeskyIntegrals::three_integral_block(const std::vector& const std::vector& q, ThreeIntsBlockOrder order) { ambit::Tensor ReturnTensor; - if (order == pqQ) { + if (order == ThreeIntsBlockOrder::pqQ) { ReturnTensor = ambit::Tensor::build(tensor_type_, "Return", {p.size(), q.size(), A.size()}); ReturnTensor.iterate([&](const std::vector& i, double& value) { value = three_integral(A[i[2]], p[i[0]], q[i[1]]); diff --git a/forte/integrals/cholesky_integrals.h b/forte/integrals/cholesky_integrals.h index 6a1e3d865..6f4c0ab2f 100644 --- a/forte/integrals/cholesky_integrals.h +++ b/forte/integrals/cholesky_integrals.h @@ -45,7 +45,7 @@ class CholeskyIntegrals : public Psi4Integrals { CholeskyIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, std::shared_ptr mo_space_info, - IntegralSpinRestriction restricted); + std::shared_ptr orbitals, IntegralSpinRestriction restricted); void initialize() override; @@ -71,9 +71,10 @@ class CholeskyIntegrals : public Psi4Integrals { double three_integral(size_t A, size_t p, size_t q) const; double** three_integral_pointer() override; - ambit::Tensor three_integral_block(const std::vector& A, const std::vector& p, - const std::vector& q, - ThreeIntsBlockOrder order = Qpq) override; + ambit::Tensor + three_integral_block(const std::vector& A, const std::vector& p, + const std::vector& q, + ThreeIntsBlockOrder order = ThreeIntsBlockOrder::Qpq) override; ambit::Tensor three_integral_block_two_index(const std::vector&, size_t, const std::vector&) override; /// Do not use this if you are using CD/DF integrals diff --git a/forte/integrals/conventional_integrals.cc b/forte/integrals/conventional_integrals.cc index d6d56c26d..23d892913 100644 --- a/forte/integrals/conventional_integrals.cc +++ b/forte/integrals/conventional_integrals.cc @@ -53,8 +53,10 @@ namespace forte { ConventionalIntegrals::ConventionalIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, IntegralSpinRestriction restricted) - : Psi4Integrals(options, ref_wfn, mo_space_info, Conventional, restricted) { + : Psi4Integrals(options, ref_wfn, mo_space_info, orbitals, IntegralType::Conventional, + restricted) { initialize(); } diff --git a/forte/integrals/conventional_integrals.h b/forte/integrals/conventional_integrals.h index 224e53d20..a8088b753 100644 --- a/forte/integrals/conventional_integrals.h +++ b/forte/integrals/conventional_integrals.h @@ -49,7 +49,7 @@ class ConventionalIntegrals : public Psi4Integrals { ConventionalIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, std::shared_ptr mo_space_info, - IntegralSpinRestriction restricted); + std::shared_ptr orbitals, IntegralSpinRestriction restricted); void initialize() override; diff --git a/forte/integrals/custom_integrals.cc b/forte/integrals/custom_integrals.cc index cf603b4c4..7aef478c5 100644 --- a/forte/integrals/custom_integrals.cc +++ b/forte/integrals/custom_integrals.cc @@ -37,10 +37,12 @@ #include "psi4/psifiles.h" #include "base_classes/mo_space_info.h" +#include "base_classes/orbitals.h" #include "helpers/blockedtensorfactory.h" #include "helpers/string_algorithms.h" #include "helpers/timer.h" #include "helpers/printing.h" +#include "helpers/helpers.h" #include "custom_integrals.h" @@ -54,23 +56,24 @@ namespace forte { CustomIntegrals::CustomIntegrals(std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, IntegralSpinRestriction restricted, double scalar, const std::vector& oei_a, const std::vector& oei_b, const std::vector& tei_aa, const std::vector& tei_ab, const std::vector& tei_bb) - : ForteIntegrals(options, mo_space_info, Custom, restricted), full_aphys_tei_aa_(tei_aa), - full_aphys_tei_ab_(tei_ab), full_aphys_tei_bb_(tei_bb) { + : ForteIntegrals(options, mo_space_info, orbitals, IntegralType::Custom, restricted), + full_aphys_tei_aa_(tei_aa), full_aphys_tei_ab_(tei_ab), full_aphys_tei_bb_(tei_bb) { set_nuclear_repulsion(scalar); set_oei_all(oei_a, oei_b); initialize(); } void CustomIntegrals::initialize() { - Ca_ = std::make_shared(nmopi_, nmopi_); - Cb_ = std::make_shared(nmopi_, nmopi_); - Ca_->identity(); - Cb_->identity(); + // Ca_ = std::make_shared(nmopi_, nmopi_); + // Cb_ = std::make_shared(nmopi_, nmopi_); + // Ca_->identity(); + // Cb_->identity(); nsopi_ = nmopi_; nso_ = nmo_; @@ -471,8 +474,8 @@ void CustomIntegrals::transform_one_electron_integrals() { } // transform the one-electron integrals - Ha->transform(Ca_); - Hb->transform(Cb_); + Ha->transform(*orbitals_->Ca()); + Hb->transform(*orbitals_->Cb()); OneBody_symm_ = Ha; @@ -506,24 +509,27 @@ void CustomIntegrals::transform_two_electron_integrals() { save_original_tei_ = true; } - auto Ca = ambit::Tensor::build(tensor_type_, "Ca", {nmo_, nmo_}); - auto Cb = ambit::Tensor::build(tensor_type_, "Cb", {nmo_, nmo_}); + auto Ca = matrix_to_tensor(orbitals_->Ca(), "Ca"); + auto Cb = matrix_to_tensor(orbitals_->Cb(), "Cb"); - auto& Ca_data = Ca.data(); - auto& Cb_data = Cb.data(); + // auto Ca = ambit::Tensor::build(tensor_type_, "Ca", {nmo_, nmo_}); + // auto Cb = ambit::Tensor::build(tensor_type_, "Cb", {nmo_, nmo_}); - int offset = 0; - for (int h = 0; h < nirrep_; ++h) { - for (int p = 0; p < nmopi_[h]; ++p) { - for (int q = 0; q < nmopi_[h]; ++q) { - int p_full = p + offset; - int q_full = q + offset; - Ca_data[p_full * nmo_ + q_full] = Ca_->get(h, p, q); - Cb_data[p_full * nmo_ + q_full] = Cb_->get(h, p, q); - } - } - offset += nmopi_[h]; - } + // auto& Ca_data = Ca.data(); + // auto& Cb_data = Cb.data(); + + // int offset = 0; + // for (int h = 0; h < nirrep_; ++h) { + // for (int p = 0; p < nmopi_[h]; ++p) { + // for (int q = 0; q < nmopi_[h]; ++q) { + // int p_full = p + offset; + // int q_full = q + offset; + // Ca_data[p_full * nmo_ + q_full] = Ca_->get(h, p, q); + // Cb_data[p_full * nmo_ + q_full] = Cb_->get(h, p, q); + // } + // } + // offset += nmopi_[h]; + // } auto T = ambit::Tensor::build(tensor_type_, "temp", {nmo_, nmo_, nmo_, nmo_}); @@ -537,17 +543,15 @@ void CustomIntegrals::transform_two_electron_integrals() { full_aphys_tei_bb_ = T.data(); } -void CustomIntegrals::update_orbitals(std::shared_ptr Ca, - std::shared_ptr Cb, bool re_transform) { +void CustomIntegrals::update_orbitals(std::shared_ptr orbitals, bool re_transform) { // 1. Copy orbitals and, if necessary, test they meet the spin restriction condition - Ca_->copy(Ca); - Cb_->copy(Cb); + orbitals_->copy(*orbitals); ints_consistent_ = false; if (spin_restriction_ == IntegralSpinRestriction::Restricted) { - if (not test_orbital_spin_restriction(Ca, Cb)) { - Ca->print(); - Cb->print(); + if (not orbitals_->are_spin_restricted()) { + orbitals_->Ca()->print(); + orbitals_->Cb()->print(); auto msg = "CustomIntegrals::update_orbitals was passed two different sets of orbitals" "\n but the integral object assumes restricted orbitals"; throw std::runtime_error(msg); diff --git a/forte/integrals/custom_integrals.h b/forte/integrals/custom_integrals.h index 11a283fdf..ec9294ed7 100644 --- a/forte/integrals/custom_integrals.h +++ b/forte/integrals/custom_integrals.h @@ -49,10 +49,11 @@ class CustomIntegrals : public ForteIntegrals { /// @param tei_ab the alpha-beta two-electron integrals in MO basis /// @param tei_bb the beta-beta two-electron integrals in MO basis CustomIntegrals(std::shared_ptr options, - std::shared_ptr mo_space_info, IntegralSpinRestriction restricted, - double scalar, const std::vector& oei_a, - const std::vector& oei_b, const std::vector& tei_aa, - const std::vector& tei_ab, const std::vector& tei_bb); + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + IntegralSpinRestriction restricted, double scalar, + const std::vector& oei_a, const std::vector& oei_b, + const std::vector& tei_aa, const std::vector& tei_ab, + const std::vector& tei_bb); void initialize() override; @@ -130,8 +131,7 @@ class CustomIntegrals : public ForteIntegrals { aptei_idx_ * r + s; } - void update_orbitals(std::shared_ptr Ca, std::shared_ptr Cb, - bool re_transform = true) override; + void update_orbitals(std::shared_ptr orbitals, bool re_transform = true) override; // ==> Class private virtual functions <== void gather_integrals() override; diff --git a/forte/integrals/df_integrals.cc b/forte/integrals/df_integrals.cc index d2ff76543..819c7d154 100644 --- a/forte/integrals/df_integrals.cc +++ b/forte/integrals/df_integrals.cc @@ -61,8 +61,8 @@ namespace forte { DFIntegrals::DFIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, std::shared_ptr mo_space_info, - IntegralSpinRestriction restricted) - : Psi4Integrals(options, ref_wfn, mo_space_info, DF, restricted) { + std::shared_ptr orbitals, IntegralSpinRestriction restricted) + : Psi4Integrals(options, ref_wfn, mo_space_info, orbitals, IntegralType::DF, restricted) { initialize(); } @@ -160,7 +160,7 @@ ambit::Tensor DFIntegrals::three_integral_block(const std::vector& A, const std::vector& q, ThreeIntsBlockOrder order) { ambit::Tensor ReturnTensor; - if (order == pqQ) { + if (order == ThreeIntsBlockOrder::pqQ) { ReturnTensor = ambit::Tensor::build(tensor_type_, "Return", {p.size(), q.size(), A.size()}); ReturnTensor.iterate([&](const std::vector& i, double& value) { value = three_integral(A[i[2]], p[i[0]], q[i[1]]); diff --git a/forte/integrals/df_integrals.h b/forte/integrals/df_integrals.h index 41793c8c6..f14de26c7 100644 --- a/forte/integrals/df_integrals.h +++ b/forte/integrals/df_integrals.h @@ -42,7 +42,8 @@ class DFIntegrals : public Psi4Integrals { public: /// Contructor of DFIntegrals DFIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, - std::shared_ptr mo_space_info, IntegralSpinRestriction restricted); + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + IntegralSpinRestriction restricted); void initialize() override; @@ -65,9 +66,10 @@ class DFIntegrals : public Psi4Integrals { double three_integral(size_t A, size_t p, size_t q); - ambit::Tensor three_integral_block(const std::vector& A, const std::vector& p, - const std::vector& q, - ThreeIntsBlockOrder order = Qpq) override; + ambit::Tensor + three_integral_block(const std::vector& A, const std::vector& p, + const std::vector& q, + ThreeIntsBlockOrder order = ThreeIntsBlockOrder::Qpq) override; ambit::Tensor three_integral_block_two_index(const std::vector&, size_t, const std::vector&) override; double** three_integral_pointer() override; diff --git a/forte/integrals/diskdf_integrals.cc b/forte/integrals/diskdf_integrals.cc index a80ecc68c..895da7f87 100644 --- a/forte/integrals/diskdf_integrals.cc +++ b/forte/integrals/diskdf_integrals.cc @@ -58,8 +58,9 @@ namespace forte { DISKDFIntegrals::DISKDFIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, IntegralSpinRestriction restricted) - : Psi4Integrals(options, ref_wfn, mo_space_info, DiskDF, restricted) { + : Psi4Integrals(options, ref_wfn, mo_space_info, orbitals, IntegralType::DiskDF, restricted) { initialize(); } @@ -301,7 +302,7 @@ ambit::Tensor DISKDFIntegrals::three_integral_block(const std::vector& Q auto Qqsize = Qsize * qsize; ambit::Tensor out; - if (order == pqQ) { + if (order == ThreeIntsBlockOrder::pqQ) { out = ambit::Tensor::build(tensor_type_, "Return", {psize, qsize, Qsize}); } else { out = ambit::Tensor::build(tensor_type_, "Return", {Qsize, psize, qsize}); @@ -363,12 +364,12 @@ ambit::Tensor DISKDFIntegrals::three_integral_block(const std::vector& Q std::vector q_range{cmotomo[q_vec[0]], cmotomo[q_vec[0]] + qsize}; df_->fill_tensor("B", out_data.data(), Q_range, p_range, q_range); - if (order == pqQ) + if (order == ThreeIntsBlockOrder::pqQ) matrix_transpose_in_place(out_data, Qsize, pqsize); } else if ((not p_contiguous) and q_contiguous) { std::vector q_range{cmotomo[q_vec[0]], cmotomo[q_vec[0]] + qsize}; - if (order == pqQ) { + if (order == ThreeIntsBlockOrder::pqQ) { for (size_t p = 0; p < psize; ++p) { auto np = cmotomo[p_vec[p]]; auto Aq = std::make_shared("Aq", Qsize, qsize); @@ -396,7 +397,7 @@ ambit::Tensor DISKDFIntegrals::three_integral_block(const std::vector& Q } else if (p_contiguous and (not q_contiguous)) { std::vector p_range{cmotomo[p_vec[0]], cmotomo[p_vec[0]] + psize}; - if (order == pqQ) { + if (order == ThreeIntsBlockOrder::pqQ) { for (size_t q = 0; q < qsize; ++q) { auto nq = cmotomo[q_vec[q]]; auto Ap = std::make_shared("Ap", Qsize, psize); @@ -441,7 +442,7 @@ ambit::Tensor DISKDFIntegrals::three_integral_block(const std::vector& Q } if (psize < qsize) { - if (order == pqQ) { + if (order == ThreeIntsBlockOrder::pqQ) { for (size_t i = 0; i < batches[n]; ++i) { for (size_t a = 0; a < Qsize; ++a) { for (size_t q = 0; q < qsize; ++q) { @@ -461,7 +462,7 @@ ambit::Tensor DISKDFIntegrals::three_integral_block(const std::vector& Q } } } else { - if (order == pqQ) { + if (order == ThreeIntsBlockOrder::pqQ) { for (size_t i = 0; i < batches[n]; ++i) { for (size_t a = 0; a < Qsize; ++a) { for (size_t p = 0; p < psize; ++p) { diff --git a/forte/integrals/diskdf_integrals.h b/forte/integrals/diskdf_integrals.h index 45c809bb9..ca362db3b 100644 --- a/forte/integrals/diskdf_integrals.h +++ b/forte/integrals/diskdf_integrals.h @@ -44,7 +44,8 @@ class DISKDFIntegrals : public Psi4Integrals { /// Contructor of DISKDFIntegrals DISKDFIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, - std::shared_ptr mo_space_info, IntegralSpinRestriction restricted); + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + IntegralSpinRestriction restricted); void initialize() override; /// aptei_xy functions are slow. try to use three_integral_block @@ -70,9 +71,10 @@ class DISKDFIntegrals : public Psi4Integrals { double** three_integral_pointer() override; /// Read a block of the DFIntegrals and return an Ambit tensor of size A by p by q - ambit::Tensor three_integral_block(const std::vector& A, const std::vector& p, - const std::vector& q, - ThreeIntsBlockOrder order = Qpq) override; + ambit::Tensor + three_integral_block(const std::vector& A, const std::vector& p, + const std::vector& q, + ThreeIntsBlockOrder order = ThreeIntsBlockOrder::Qpq) override; /// return ambit tensor of size A by q ambit::Tensor three_integral_block_two_index(const std::vector& A, size_t p, const std::vector& q) override; diff --git a/forte/integrals/integrals.cc b/forte/integrals/integrals.cc index b2aa860fb..ffb2c1e05 100644 --- a/forte/integrals/integrals.cc +++ b/forte/integrals/integrals.cc @@ -37,6 +37,7 @@ #include "helpers/blockedtensorfactory.h" #include "base_classes/forte_options.h" #include "base_classes/mo_space_info.h" +#include "base_classes/orbitals.h" #include "helpers/printing.h" #include "helpers/timer.h" #include "integrals.h" @@ -60,25 +61,30 @@ namespace forte { #endif std::map int_type_label{ - {Conventional, "Conventional"}, {DF, "Density fitting"}, - {Cholesky, "Cholesky decomposition"}, {DiskDF, "Disk-based density fitting"}, - {DistDF, "Distributed density fitting"}, {Custom, "Custom"}}; + {IntegralType::Conventional, "Conventional"}, + {IntegralType::DF, "Density fitting"}, + {IntegralType::Cholesky, "Cholesky decomposition"}, + {IntegralType::DiskDF, "Disk-based density fitting"}, + {IntegralType::DistDF, "Distributed density fitting"}, + {IntegralType::Custom, "Custom"}}; ForteIntegrals::ForteIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, std::shared_ptr mo_space_info, - IntegralType integral_type, IntegralSpinRestriction restricted) - : options_(options), mo_space_info_(mo_space_info), wfn_(ref_wfn), - integral_type_(integral_type), spin_restriction_(restricted), frozen_core_energy_(0.0), - scalar_energy_(0.0) { + std::shared_ptr orbitals, IntegralType integral_type, + IntegralSpinRestriction restricted, DFType df_type) + : options_(options), mo_space_info_(mo_space_info), orbitals_(orbitals), wfn_(ref_wfn), + integral_type_(integral_type), spin_restriction_(restricted), df_type_(df_type), + frozen_core_energy_(0.0), scalar_energy_(0.0) { common_initialize(); } ForteIntegrals::ForteIntegrals(std::shared_ptr options, std::shared_ptr mo_space_info, - IntegralType integral_type, IntegralSpinRestriction restricted) - : options_(options), mo_space_info_(mo_space_info), integral_type_(integral_type), - spin_restriction_(restricted) { + std::shared_ptr orbitals, IntegralType integral_type, + IntegralSpinRestriction restricted, DFType df_type) + : options_(options), mo_space_info_(mo_space_info), orbitals_(orbitals), + integral_type_(integral_type), spin_restriction_(restricted), df_type_(df_type) { common_initialize(); } @@ -122,7 +128,7 @@ void ForteIntegrals::read_information() { // skip integral allocation and transformation if doing CASSCF auto job_type = options_->get_str("JOB_TYPE"); - skip_build_ = (job_type == "MCSCF_TWO_STEP") and (integral_type_ != Custom); + skip_build_ = (job_type == "MCSCF_TWO_STEP") and (integral_type_ != IntegralType::Custom); } void ForteIntegrals::allocate() { @@ -134,7 +140,8 @@ void ForteIntegrals::allocate() { one_electron_integrals_a_.assign(ncmo_ * ncmo_, 0.0); one_electron_integrals_b_.assign(ncmo_ * ncmo_, 0.0); - if ((integral_type_ == Conventional) or (integral_type_ == Custom)) { + if ((integral_type_ == IntegralType::Conventional) or + (integral_type_ == IntegralType::Custom)) { // Allocate the memory required to store the two-electron integrals aphys_tei_aa_.assign(num_aptei_, 0.0); aphys_tei_ab_.assign(num_aptei_, 0.0); @@ -144,15 +151,13 @@ void ForteIntegrals::allocate() { } } -std::shared_ptr ForteIntegrals::Ca() const { return Ca_; } - -std::shared_ptr ForteIntegrals::Cb() const { return Cb_; } +std::shared_ptr ForteIntegrals::orbitals() const { return orbitals_; } bool ForteIntegrals::update_ints_if_needed() { if (ints_consistent_) { return false; } - update_orbitals(Ca_, Cb_, true); + update_orbitals(orbitals_, true); return true; } @@ -168,7 +173,7 @@ void ForteIntegrals::jk_finalize() { if (JK_status_ == JKStatus::initialized) { JK_->finalize(); // nothing done in finalize() for PKJK and MemDFJK - if (integral_type_ == DiskDF or integral_type_ == Cholesky) + if (integral_type_ == IntegralType::DiskDF or integral_type_ == IntegralType::Cholesky) JK_status_ = JKStatus::finalized; } } @@ -322,7 +327,7 @@ void ForteIntegrals::set_oei(size_t p, size_t q, double value, bool alpha) { } bool ForteIntegrals::fix_orbital_phases(std::shared_ptr U, bool is_alpha, bool debug) { - if (integral_type_ == Custom) { + if (integral_type_ == IntegralType::Custom) { outfile->Printf("\n Warning: Cannot fix orbital phases (%s) for CustomIntegrals.", is_alpha ? "Ca" : "Cb"); return false; @@ -392,13 +397,6 @@ bool ForteIntegrals::fix_orbital_phases(std::shared_ptr U, bool is_ } } -bool ForteIntegrals::test_orbital_spin_restriction(std::shared_ptr A, - std::shared_ptr B) const { - auto A_minus_B = A->clone(); - A_minus_B->subtract(B); - return A_minus_B->absmax() < 1.0e-7; -} - void ForteIntegrals::freeze_core_orbitals() { local_timer freeze_timer; if (ncmo_ < nmo_) { @@ -491,16 +489,13 @@ void ForteIntegrals::print_ints() { void ForteIntegrals::rotate_orbitals(std::shared_ptr Ua, std::shared_ptr Ub, bool re_transform) { // 1. Rotate the orbital coefficients and store them in the ForteIntegral object - auto Ca_rotated = psi::linalg::doublet(Ca_, Ua); - auto Cb_rotated = psi::linalg::doublet(Cb_, Ub); - - update_orbitals(Ca_rotated, Cb_rotated, re_transform); + orbitals_->rotate(Ua, Ub); + update_orbitals(orbitals_, re_transform); } // The following functions throw an error by default -void ForteIntegrals::update_orbitals(std::shared_ptr, std::shared_ptr, - bool) { +void ForteIntegrals::update_orbitals(std::shared_ptr, bool) { _undefined_function("update_orbitals"); } @@ -547,10 +542,10 @@ std::vector> ForteIntegrals::mo_quadrupole_ints() c } void ForteIntegrals::_undefined_function(const std::string& method) const { - outfile->Printf("\n ForteIntegrals::" + method + "not supported for integral type " + - std::to_string(integral_type())); - throw std::runtime_error("ForteIntegrals::" + method + " not supported for integral type " + - std::to_string(integral_type())); + std::string msg = "ForteIntegrals::" + method + " not supported for integral type " + + int_type_label[integral_type()]; + outfile->Printf("\n %s\n", msg.c_str()); + throw std::runtime_error(msg); } } // namespace forte diff --git a/forte/integrals/integrals.h b/forte/integrals/integrals.h index e613a25fc..99cf34ab3 100644 --- a/forte/integrals/integrals.h +++ b/forte/integrals/integrals.h @@ -50,6 +50,7 @@ namespace forte { class ForteOptions; class MOSpaceInfo; +class Orbitals; /** * @brief The IntegralSpinRestriction enum @@ -63,7 +64,7 @@ enum class IntegralSpinRestriction { Restricted, Unrestricted }; * * This decides the type of integral used in a Forte computation */ -enum IntegralType { Conventional, DF, Cholesky, DiskDF, DistDF, Custom }; +enum class IntegralType { Conventional, DF, Cholesky, DiskDF, DistDF, Custom }; /** * @brief The order of three-index integrals when calling @@ -71,7 +72,14 @@ enum IntegralType { Conventional, DF, Cholesky, DiskDF, DistDF, Custom }; * In Forte, the auxiliary index is by default the first index (i.e., Qpq). * However, pqQ order is more convenient for implementing DF-MRPT2. */ -enum ThreeIntsBlockOrder { Qpq, pqQ }; +enum class ThreeIntsBlockOrder { Qpq, pqQ }; + +/** + * @brief The type of density fitting basis + * + * This is used to distinguish between JK and RI density fitting basis sets + */ +enum class DFType { None, JK, RI }; /** * @brief The ForteIntegrals class is a base class for transforming and storing MO integrals @@ -115,8 +123,9 @@ class ForteIntegrals { */ ForteIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, - std::shared_ptr mo_space_info, IntegralType integral_type, - IntegralSpinRestriction restricted); + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + IntegralType integral_type, IntegralSpinRestriction restricted, + DFType df_type = DFType::JK); /** * @brief Class constructor @@ -125,8 +134,9 @@ class ForteIntegrals { * @param mo_space_info The MOSpaceInfo object */ ForteIntegrals(std::shared_ptr options, - std::shared_ptr mo_space_info, IntegralType integral_type, - IntegralSpinRestriction restricted); + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + IntegralType integral_type, IntegralSpinRestriction restricted, + DFType df_type = DFType::JK); /// Virtual destructor to enable deletion of a Derived* through a Base* virtual ~ForteIntegrals() = default; @@ -149,9 +159,7 @@ class ForteIntegrals { /// Use the function update_ints_if_needed to re-transform the integrals only if they changed. /// /// @return the coefficient matrix for the alpha orbitals used to transform the integrals - std::shared_ptr Ca() const; - /// Return Cb - std::shared_ptr Cb() const; + std::shared_ptr orbitals() const; /// Return nuclear repulsion energy double nuclear_repulsion_energy() const; @@ -263,10 +271,10 @@ class ForteIntegrals { const std::vector& s) = 0; // Three-index integral functions (DF, Cholesky) - virtual ambit::Tensor three_integral_block(const std::vector&, - const std::vector&, - const std::vector&, - ThreeIntsBlockOrder order = Qpq); + virtual ambit::Tensor + three_integral_block(const std::vector&, const std::vector&, + const std::vector&, + ThreeIntsBlockOrder order = ThreeIntsBlockOrder::Qpq); /// This function is only used by DiskDF and it is used to go from a Apq->Aq tensor virtual ambit::Tensor three_integral_block_two_index(const std::vector& A, size_t p, @@ -349,17 +357,15 @@ class ForteIntegrals { /// Rotate the MO coefficients, update psi::Wavefunction, and re-transform integrals /// @param Ua the alpha unitary transformation matrix /// @param Ub the beta unitary transformation matrix - /// @param re_transform re-transform integrals if true + /// @param re_transform if true, re-transform integrals, otherwise just update MOs void rotate_orbitals(std::shared_ptr Ua, std::shared_ptr Ub, bool re_transform = true); - /// Copy these MO coeffs to class variables, update psi::Wavefunction, and re-transform + /// Copy the orbital object to class variables, update psi::Wavefunction, and re-transform /// integrals - /// @param Ca the alpha MO coefficients - /// @param Cb the beta MO coefficients + /// @param orbitals the new orbitals /// @param re_transform re-transform integrals if true - virtual void update_orbitals(std::shared_ptr Ca, std::shared_ptr Cb, - bool re_transform = true); + virtual void update_orbitals(std::shared_ptr orbitals, bool re_transform = true); /// Update the integrals if the MO coefficients have changed but the integrals were not /// re-transformed @@ -418,6 +424,9 @@ class ForteIntegrals { /// The MOSpaceInfo object std::shared_ptr mo_space_info_; + /// The orbitals object + std::shared_ptr orbitals_; + /// The Wavefunction object std::shared_ptr wfn_; @@ -427,11 +436,8 @@ class ForteIntegrals { /// Are we doing a spin-restricted computation? IntegralSpinRestriction spin_restriction_; - // Ca matrix from psi - std::shared_ptr Ca_; - - // Cb matrix from psi - std::shared_ptr Cb_; + /// The type of density fitting basis + DFType df_type_; // AO overlap matrix from psi std::shared_ptr S_; @@ -586,8 +592,8 @@ class ForteIntegrals { class Psi4Integrals : public ForteIntegrals { public: Psi4Integrals(std::shared_ptr options, std::shared_ptr ref_wfn, - std::shared_ptr mo_space_info, IntegralType integral_type, - IntegralSpinRestriction restricted); + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + IntegralType integral_type, IntegralSpinRestriction restricted); /// Make the generalized Fock matrix using Psi4 JK object void make_fock_matrix(ambit::Tensor Da, ambit::Tensor Db) override; @@ -623,8 +629,7 @@ class Psi4Integrals : public ForteIntegrals { void setup_psi4_ints(); void transform_one_electron_integrals(); void compute_frozen_one_body_operator() override; - void update_orbitals(std::shared_ptr Ca, std::shared_ptr Cb, - bool re_transform = true) override; + void update_orbitals(std::shared_ptr orbitals, bool re_transform = true) override; void rotate_mos() override; /// Build AO dipole and quadrupole integrals diff --git a/forte/integrals/integrals_psi4_interface.cc b/forte/integrals/integrals_psi4_interface.cc index 3c47d7cb6..ca46e8ba7 100644 --- a/forte/integrals/integrals_psi4_interface.cc +++ b/forte/integrals/integrals_psi4_interface.cc @@ -30,9 +30,7 @@ #include "psi4/psi4-dec.h" #include "psi4/libpsi4util/PsiOutStream.h" #include "psi4/libmints/vector.h" - #include "psi4/libfock/jk.h" - #include "psi4/libmints/basisset.h" #include "psi4/libmints/integral.h" #include "psi4/libmints/matrix.h" @@ -41,9 +39,12 @@ #include "psi4/libpsi4util/process.h" #include "base_classes/mo_space_info.h" +#include "base_classes/orbitals.h" + +#include "integrals/integrals.h" + #include "helpers/printing.h" #include "helpers/timer.h" -#include "integrals/integrals.h" #ifdef HAVE_GA #include @@ -63,9 +64,10 @@ namespace forte { Psi4Integrals::Psi4Integrals(std::shared_ptr options, std::shared_ptr ref_wfn, - std::shared_ptr mo_space_info, IntegralType integral_type, + std::shared_ptr mo_space_info, + std::shared_ptr orbitals, IntegralType integral_type, IntegralSpinRestriction restricted) - : ForteIntegrals(options, ref_wfn, mo_space_info, integral_type, restricted) { + : ForteIntegrals(options, ref_wfn, mo_space_info, orbitals, integral_type, restricted) { base_initialize_psi4(); } @@ -87,10 +89,12 @@ void Psi4Integrals::setup_psi4_ints() { exit(1); } + // ADDORBITALS // Grab the MO coefficients from psi and enforce spin restriction if necessary - Ca_ = wfn_->Ca()->clone(); - Cb_ = (spin_restriction_ == IntegralSpinRestriction::Restricted ? wfn_->Ca()->clone() - : wfn_->Cb()->clone()); + // Ca_ = wfn_->Ca()->clone(); + // Cb_ = (spin_restriction_ == IntegralSpinRestriction::Restricted ? wfn_->Ca()->clone() + // : wfn_->Cb()->clone()); + S_ = wfn_->S()->clone(); nso_ = wfn_->nso(); @@ -113,8 +117,8 @@ void Psi4Integrals::transform_one_electron_integrals() { auto Ha = wfn_->H()->clone(); auto Hb = wfn_->H()->clone(); - Ha->transform(Ca_); - Hb->transform(Cb_); + Ha->transform(*orbitals_->Ca()); + Hb->transform(*orbitals_->Cb()); OneBody_symm_ = Ha; @@ -152,10 +156,10 @@ void Psi4Integrals::make_psi4_JK() { outfile->Printf("\n\n ==> Primary Basis Set Summary <==\n\n"); basis->print(); - if (integral_type_ == Conventional) { + if (integral_type_ == IntegralType::Conventional) { outfile->Printf("\n JK created using conventional PK integrals\n"); JK_ = JK::build_JK(basis, psi::BasisSet::zero_ao_basis_set(), psi4_options, "PK"); - } else if (integral_type_ == Cholesky) { + } else if (integral_type_ == IntegralType::Cholesky) { if (spin_restriction_ == IntegralSpinRestriction::Unrestricted) { throw psi::PSIEXCEPTION("Unrestricted orbitals not supported for CD integrals"); } @@ -173,7 +177,8 @@ void Psi4Integrals::make_psi4_JK() { } JK_ = JK::build_JK(basis, psi::BasisSet::zero_ao_basis_set(), psi4_options, "CD"); - } else if ((integral_type_ == DF) or (integral_type_ == DiskDF) or (integral_type_ == DistDF)) { + } else if ((integral_type_ == IntegralType::DF) or (integral_type_ == IntegralType::DiskDF) or + (integral_type_ == IntegralType::DistDF)) { if (spin_restriction_ == IntegralSpinRestriction::Unrestricted) { throw psi::PSIEXCEPTION("Unrestricted orbitals not supported for DF integrals"); } @@ -190,7 +195,7 @@ void Psi4Integrals::make_psi4_JK() { if (job_type == "CASSCF" or job_type == "MCSCF_TWO_STEP") basis_aux = wfn_->get_basisset("DF_BASIS_SCF"); - if (integral_type_ == DiskDF) { + if (integral_type_ == IntegralType::DiskDF) { outfile->Printf("\n JK created using DiskDF integrals\n"); JK_ = JK::build_JK(basis, basis_aux, psi4_options, "DISK_DF"); } else { @@ -263,21 +268,20 @@ void Psi4Integrals::compute_frozen_one_body_operator() { } } -void Psi4Integrals::update_orbitals(std::shared_ptr Ca, - std::shared_ptr Cb, bool re_transform) { +void Psi4Integrals::update_orbitals(std::shared_ptr orbitals, bool re_transform) { // 1. Copy orbitals and set the invalid flag - Ca_->copy(Ca); - Cb_->copy(Cb); + orbitals_->copy(*orbitals); ints_consistent_ = false; // if necessary, test they meet the spin restriction condition if (spin_restriction_ == IntegralSpinRestriction::Restricted) { - if (not test_orbital_spin_restriction(Ca, Cb)) { - Ca->print(); - Cb->print(); - auto overlap = psi::linalg::triplet(Ca, S_, Cb, true, false, false); - overlap->set_name("Overlap "); - overlap->print(); + if (not orbitals->are_spin_restricted()) { + orbitals_->Ca()->print(); + orbitals_->Cb()->print(); + auto overlap = + psi::linalg::triplet(*orbitals->Ca(), *S_, *orbitals->Cb(), true, false, false); + overlap.set_name("Overlap "); + overlap.print(); auto msg = "Psi4Integrals::update_orbitals was passed two different sets of orbitals" "\n but the integral object assumes restricted orbitals"; throw std::runtime_error(msg); @@ -285,8 +289,8 @@ void Psi4Integrals::update_orbitals(std::shared_ptr Ca, } // 2. Send a copy to psi::Wavefunction - wfn_->Ca()->copy(Ca_); - wfn_->Cb()->copy(Cb_); + wfn_->Ca()->copy(*orbitals_->Ca()); + wfn_->Cb()->copy(*orbitals_->Cb()); // 3. Re-transform the integrals if (re_transform) { @@ -349,7 +353,7 @@ void Psi4Integrals::rotate_mos() { rotate_mo_group[2]); } - auto C_old = Ca_; + auto C_old = orbitals_->Ca()->clone(); auto C_new(C_old->clone()); const auto& eps_a_old = *wfn_->epsilon_a(); @@ -366,8 +370,7 @@ void Psi4Integrals::rotate_mos() { eps_a_new.set(mo_group[0], mo_group[1], epsilon_mo2); } // Update local copy of the orbitals - Ca_->copy(C_new); - Cb_->copy(C_new); + orbitals_->set(C_new, C_new); // Copy to psi::Wavefunction wfn_->Ca()->copy(C_new); @@ -390,8 +393,9 @@ std::shared_ptr Psi4Integrals::Ca_AO() const { for (int i = 0, nmo_this = nmopi_[h]; i < nmo_this; ++i) { // notes: LDA value is nso (not nao, see libqt/blas_intfc23.cc) - C_DGEMV('N', nao, nso, 1.0, aotoso->pointer(h)[0], nso, &Ca_->pointer(h)[0][i], - nmo_this, 0.0, &Ca_ao->pointer()[0][index++], nmo_); + C_DGEMV('N', nao, nso, 1.0, aotoso->pointer(h)[0], nso, + &orbitals_->Ca()->pointer(h)[0][i], nmo_this, 0.0, + &Ca_ao->pointer()[0][index++], nmo_); } } return Ca_ao; @@ -492,7 +496,7 @@ Psi4Integrals::make_fock_inactive(psi::Dimension dim_start, psi::Dimension dim_e for (int h = 0; h < nirrep_; ++h) { for (int p = 0, offset = dim_start[h]; p < dim[h]; ++p) { - Csub->set_column(h, p, Ca_->get_column(h, p + offset)); + Csub->set_column(h, p, orbitals_->Ca()->get_column(h, p + offset)); } } @@ -513,7 +517,8 @@ Psi4Integrals::make_fock_inactive(psi::Dimension dim_start, psi::Dimension dim_e J->add(wfn_->H()); // transform to MO - auto F_closed = psi::linalg::triplet(Ca_, J, Ca_, true, false, false); + auto F_closed = + psi::linalg::triplet(orbitals_->Ca(), J, orbitals_->Ca(), true, false, false); F_closed->set_name("Fock_closed"); // compute closed-shell energy @@ -534,8 +539,8 @@ Psi4Integrals::make_fock_inactive(psi::Dimension dim_start, psi::Dimension dim_e for (int h = 0; h < nirrep_; ++h) { for (int p = 0, offset = dim_start[h]; p < dim[h]; ++p) { - Ca_sub->set_column(h, p, Ca_->get_column(h, p + offset)); - Cb_sub->set_column(h, p, Cb_->get_column(h, p + offset)); + Ca_sub->set_column(h, p, orbitals_->Ca()->get_column(h, p + offset)); + Cb_sub->set_column(h, p, orbitals_->Cb()->get_column(h, p + offset)); } } @@ -570,9 +575,11 @@ Psi4Integrals::make_fock_inactive(psi::Dimension dim_start, psi::Dimension dim_e K->copy(Fb_closed); // transform to MO basis - Fa_closed = psi::linalg::triplet(Ca_, Fa_closed, Ca_, true, false, false); + Fa_closed = + psi::linalg::triplet(orbitals_->Ca(), Fa_closed, orbitals_->Ca(), true, false, false); Fa_closed->set_name("Fock_closed alpha"); - Fb_closed = psi::linalg::triplet(Cb_, Fb_closed, Cb_, true, false, false); + Fb_closed = + psi::linalg::triplet(orbitals_->Cb(), Fb_closed, orbitals_->Cb(), true, false, false); Fb_closed->set_name("Fock_closed beta"); // compute closed-shell energy using unrestricted equation @@ -678,7 +685,7 @@ Psi4Integrals::make_fock_active_restricted(std::shared_ptr g1) { for (int h = 0; h < nirrep_; ++h) { for (int p = 0, offset = ndoccpi[h]; p < nactvpi[h]; ++p) { - Cactv->set_column(h, p, Ca_->get_column(h, p + offset)); + Cactv->set_column(h, p, orbitals_->Ca()->get_column(h, p + offset)); } } @@ -702,7 +709,7 @@ Psi4Integrals::make_fock_active_restricted(std::shared_ptr g1) { K->add(JK_->J()[0]); // transform to MO - auto F_active = psi::linalg::triplet(Ca_, K, Ca_, true, false, false); + auto F_active = psi::linalg::triplet(orbitals_->Ca(), K, orbitals_->Ca(), true, false, false); F_active->set_name("Fock_active"); // pass AO fock to psi4 Wavefunction @@ -731,8 +738,8 @@ Psi4Integrals::make_fock_active_unrestricted(std::shared_ptr g1a, for (int h = 0; h < nirrep_; ++h) { for (int p = 0, offset = ndoccpi[h]; p < nactvpi[h]; ++p) { - Ca_actv->set_column(h, p, Ca_->get_column(h, p + offset)); - Cb_actv->set_column(h, p, Cb_->get_column(h, p + offset)); + Ca_actv->set_column(h, p, orbitals_->Ca()->get_column(h, p + offset)); + Cb_actv->set_column(h, p, orbitals_->Ca()->get_column(h, p + offset)); } } @@ -766,9 +773,9 @@ Psi4Integrals::make_fock_active_unrestricted(std::shared_ptr g1a, Kb->add(JK_->J()[1]); // transform to MO - auto Fa_active = psi::linalg::triplet(Ca_, Ka, Ca_, true, false, false); + auto Fa_active = psi::linalg::triplet(orbitals_->Ca(), Ka, orbitals_->Ca(), true, false, false); Fa_active->set_name("Fock_active alpha"); - auto Fb_active = psi::linalg::triplet(Cb_, Kb, Cb_, true, false, false); + auto Fb_active = psi::linalg::triplet(orbitals_->Cb(), Kb, orbitals_->Cb(), true, false, false); Fb_active->set_name("Fock_active beta"); // pass AO fock to psi4 Wavefunction diff --git a/forte/integrals/make_integrals.cc b/forte/integrals/make_integrals.cc index 35a288053..dd2e26b34 100644 --- a/forte/integrals/make_integrals.cc +++ b/forte/integrals/make_integrals.cc @@ -46,7 +46,8 @@ namespace forte { std::shared_ptr make_forte_integrals_from_psi4(std::shared_ptr ref_wfn, std::shared_ptr options, - std::shared_ptr mo_space_info, std::string int_type) { + std::shared_ptr mo_space_info, + std::shared_ptr orbitals, std::string int_type) { timer int_timer("Integrals"); std::shared_ptr ints; // passing int_type overrides the value of the option INT_TYPE @@ -56,16 +57,16 @@ make_forte_integrals_from_psi4(std::shared_ptr ref_wfn, int_type = options->get_str("INT_TYPE"); } if (int_type == "CHOLESKY") { - ints = std::make_shared(options, ref_wfn, mo_space_info, + ints = std::make_shared(options, ref_wfn, mo_space_info, orbitals, IntegralSpinRestriction::Restricted); } else if (int_type == "DF") { - ints = std::make_shared(options, ref_wfn, mo_space_info, + ints = std::make_shared(options, ref_wfn, mo_space_info, orbitals, IntegralSpinRestriction::Restricted); } else if (int_type == "DISKDF") { - ints = std::make_shared(options, ref_wfn, mo_space_info, + ints = std::make_shared(options, ref_wfn, mo_space_info, orbitals, IntegralSpinRestriction::Restricted); } else if (int_type == "CONVENTIONAL") { - ints = std::make_shared(options, ref_wfn, mo_space_info, + ints = std::make_shared(options, ref_wfn, mo_space_info, orbitals, IntegralSpinRestriction::Restricted); } else if (int_type == "DISTDF") { #ifdef HAVE_GA @@ -86,13 +87,12 @@ make_forte_integrals_from_psi4(std::shared_ptr ref_wfn, return ints; } -std::shared_ptr -make_custom_forte_integrals(std::shared_ptr options, - std::shared_ptr mo_space_info, double scalar, - const std::vector& oei_a, const std::vector& oei_b, - const std::vector& tei_aa, const std::vector& tei_ab, - const std::vector& tei_bb) { - return std::make_shared(options, mo_space_info, +std::shared_ptr make_custom_forte_integrals( + std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, double scalar, const std::vector& oei_a, + const std::vector& oei_b, const std::vector& tei_aa, + const std::vector& tei_ab, const std::vector& tei_bb) { + return std::make_shared(options, mo_space_info, orbitals, IntegralSpinRestriction::Restricted, scalar, oei_a, oei_b, tei_aa, tei_ab, tei_bb); } diff --git a/forte/integrals/make_integrals.h b/forte/integrals/make_integrals.h index 677599c59..404710677 100644 --- a/forte/integrals/make_integrals.h +++ b/forte/integrals/make_integrals.h @@ -37,18 +37,19 @@ namespace forte { * but if the variable int_type is provided, its value will override * the value read from the options object. */ -std::shared_ptr make_forte_integrals_from_psi4( - std::shared_ptr ref_wfn, std::shared_ptr options, - std::shared_ptr mo_space_info, std::string int_type = ""); +std::shared_ptr +make_forte_integrals_from_psi4(std::shared_ptr ref_wfn, + std::shared_ptr options, + std::shared_ptr mo_space_info, + std::shared_ptr orbitals, std::string int_type = ""); /** * @brief Make a ForteIntegrals object by passing integrals stored in vectors */ -std::shared_ptr -make_custom_forte_integrals(std::shared_ptr options, - std::shared_ptr mo_space_info, double scalar, - const std::vector& oei_a, const std::vector& oei_b, - const std::vector& tei_aa, const std::vector& tei_ab, - const std::vector& tei_bb); +std::shared_ptr make_custom_forte_integrals( + std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, double scalar, const std::vector& oei_a, + const std::vector& oei_b, const std::vector& tei_aa, + const std::vector& tei_ab, const std::vector& tei_bb); } // namespace forte diff --git a/forte/mcscf/mcscf_2step.cc b/forte/mcscf/mcscf_2step.cc index 28c47d6dc..ee06a610b 100644 --- a/forte/mcscf/mcscf_2step.cc +++ b/forte/mcscf/mcscf_2step.cc @@ -57,10 +57,11 @@ MCSCF_2STEP::MCSCF_2STEP(std::shared_ptr as_solver, const std::map>& state_weights_map, std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, std::shared_ptr scf_info, std::shared_ptr ints) : as_solver_(as_solver), state_weights_map_(state_weights_map), options_(options), - mo_space_info_(mo_space_info), scf_info_(scf_info), ints_(ints) { + mo_space_info_(mo_space_info), orbitals_(orbitals), scf_info_(scf_info), ints_(ints) { startup(); } @@ -81,7 +82,7 @@ void MCSCF_2STEP::read_options() { int_type_ = options_->get_str("INT_TYPE"); der_type_ = options_->get_str("DERTYPE"); - if (der_type_ == "FIRST" and ints_->integral_type() == Custom) + if (der_type_ == "FIRST" and ints_->integral_type() == IntegralType::Custom) throw std::runtime_error("MCSCF energy gradient not available for CUSTOM integrals!"); maxiter_ = options_->get_int("MCSCF_MAXITER"); @@ -169,7 +170,7 @@ double MCSCF_2STEP::compute_energy() { // prepare for orbital gradients const bool freeze_core = options_->get_bool("MCSCF_FREEZE_CORE"); - MCSCF_ORB_GRAD cas_grad(options_, mo_space_info_, ints_, freeze_core); + MCSCF_ORB_GRAD cas_grad(options_, mo_space_info_, orbitals_, ints_, freeze_core); auto nrot = cas_grad.nrot(); auto dG = std::make_shared("dG", nrot); @@ -412,7 +413,7 @@ double MCSCF_2STEP::compute_energy() { diagonalize_hamiltonian(as_solver_, cas_grad.active_space_ints(), {print_, e_conv_, r_conv, options_->get_bool("DUMP_ACTIVE_WFN")}); - if (ints_->integral_type() != Custom) { + if (ints_->integral_type() != IntegralType::Custom) { auto final_orbs = options_->get_str("MCSCF_FINAL_ORBITAL"); if (final_orbs != "UNSPECIFIED" or der_type_ == "FIRST") { @@ -594,10 +595,10 @@ std::unique_ptr make_mcscf_two_step(std::shared_ptr as_solver, const std::map>& state_weight_map, std::shared_ptr scf_info, std::shared_ptr options, - std::shared_ptr mo_space_info, + std::shared_ptr mo_space_info, std::shared_ptr orbitals, std::shared_ptr ints) { return std::make_unique(as_solver, state_weight_map, options, mo_space_info, - scf_info, ints); + orbitals, scf_info, ints); } } // namespace forte diff --git a/forte/mcscf/mcscf_2step.h b/forte/mcscf/mcscf_2step.h index ed0a78bd3..936985fa7 100644 --- a/forte/mcscf/mcscf_2step.h +++ b/forte/mcscf/mcscf_2step.h @@ -57,7 +57,8 @@ class MCSCF_2STEP { MCSCF_2STEP(std::shared_ptr as_solver, const std::map>& state_weights_map, std::shared_ptr options, std::shared_ptr mo_space_info, - std::shared_ptr scf_info, std::shared_ptr ints); + std::shared_ptr orbitals, std::shared_ptr scf_info, + std::shared_ptr ints); /// Compute the MCSCF energy double compute_energy(); @@ -75,6 +76,9 @@ class MCSCF_2STEP { /// The MOSpaceInfo object std::shared_ptr mo_space_info_; + /// The orbitals object + std::shared_ptr orbitals_; + /// SCF information std::shared_ptr scf_info_; @@ -177,7 +181,7 @@ std::unique_ptr make_mcscf_two_step(std::shared_ptr as_solver, const std::map>& state_weight_map, std::shared_ptr ref_wfn, std::shared_ptr options, - std::shared_ptr mo_space_info, + std::shared_ptr mo_space_info, std::shared_ptr orbitals, std::shared_ptr ints); } // namespace forte diff --git a/forte/mcscf/mcscf_orb_grad.cc b/forte/mcscf/mcscf_orb_grad.cc index be870a298..c0a7ed98e 100644 --- a/forte/mcscf/mcscf_orb_grad.cc +++ b/forte/mcscf/mcscf_orb_grad.cc @@ -50,6 +50,7 @@ #include "integrals/integrals.h" #include "integrals/active_space_integrals.h" #include "base_classes/rdms.h" +#include "base_classes/orbitals.h" #include "mcscf/mcscf_orb_grad.h" @@ -59,9 +60,11 @@ using namespace ambit; namespace forte { MCSCF_ORB_GRAD::MCSCF_ORB_GRAD(std::shared_ptr options, - std::shared_ptr mo_space_info, - std::shared_ptr ints, bool freeze_core) - : options_(options), mo_space_info_(mo_space_info), ints_(ints), freeze_core_(freeze_core) { + std::shared_ptr mo_space_info, + std::shared_ptr orbitals, + std::shared_ptr ints, bool freeze_core) + : options_(options), mo_space_info_(mo_space_info), orbitals_(orbitals), ints_(ints), + freeze_core_(freeze_core) { startup(); } @@ -406,9 +409,9 @@ void MCSCF_ORB_GRAD::nonredundant_pairs() { void MCSCF_ORB_GRAD::init_tensors() { // save a copy of initial MO - C0_ = ints_->Ca()->clone(); + C0_ = orbitals_->Ca()->clone(); C0_->set_name("MCSCF Initial Orbital Coefficients"); - C_ = ints_->Ca()->clone(); + C_ = orbitals_->Ca()->clone(); C_->set_name("MCSCF Orbital Coefficients"); // Fock matrices @@ -453,7 +456,7 @@ void MCSCF_ORB_GRAD::build_mo_integrals() { // form closed-shell Fock matrix build_fock_inactive(); // form the MO 2e-integrals - if (ints_->integral_type() == Custom) { + if (ints_->integral_type() == IntegralType::Custom) { fill_tei_custom(V_); } else { build_tei_from_ao(); @@ -693,7 +696,7 @@ std::shared_ptr MCSCF_ORB_GRAD::fock(std::shared_ptr rdms) { } double MCSCF_ORB_GRAD::evaluate(std::shared_ptr x, std::shared_ptr g, - bool do_g) { + bool do_g) { // if need to update orbitals and integrals if (update_orbitals(x)) { build_mo_integrals(); @@ -748,15 +751,20 @@ bool MCSCF_ORB_GRAD::update_orbitals(std::shared_ptr x) { U_ = psi::linalg::doublet(U_, dU, false, false); U_->set_name("Orthogonal Transformation"); - // update orbitals + // update local orbitals C_->gemm(false, false, 1.0, C0_, U_, 0.0); - if (ints_->integral_type() == Custom) { - ints_->update_orbitals(C_, C_); - } else { - // here we push the new orbitals to the integrals object but do not update the integrals - // for efficiency - ints_->update_orbitals(C_, C_, false); - } + + // update orbital object + orbitals_->set(C_, C_); + // update integrals object. Except for custom integrals, we update the integrals + // object with the new orbitals but do not update the integrals for efficiency + ints_->update_orbitals(orbitals_, ints_->integral_type() != IntegralType::Custom); + // if (ints_->integral_type() != Custom) { + // ints_->update_orbitals(orbitals_); + // } else { + // // here we push the new orbitals to the integrals object but do not update the integrals + // // for efficiency + // } // printing if (debug_print_) { @@ -812,7 +820,7 @@ std::shared_ptr MCSCF_ORB_GRAD::cayley_trans(const std::shared_ptr< std::vector> MCSCF_ORB_GRAD::test_orbital_rotations(const std::shared_ptr& U, - const std::string& warning_msg) { + const std::string& warning_msg) { // the overlap between new and old orbitals is simply U // O = Cold^T S Cnew = Cold^T S Cold U = U // MOM projection index: Pj = sum_{i}^{actv} O_ij @@ -900,7 +908,7 @@ void MCSCF_ORB_GRAD::compute_orbital_grad() { } void MCSCF_ORB_GRAD::hess_diag(std::shared_ptr, - const std::shared_ptr& h0) { + const std::shared_ptr& h0) { compute_orbital_hess_diag(); h0->copy(*hess_diag_); } @@ -1006,7 +1014,7 @@ void MCSCF_ORB_GRAD::compute_orbital_hess_diag() { } void MCSCF_ORB_GRAD::reshape_rot_ambit(ambit::BlockedTensor bt, - const std::shared_ptr& sv) { + const std::shared_ptr& sv) { size_t vec_size = sv->dimpi().sum(); if (vec_size != nrot_) { throw std::runtime_error( @@ -1084,11 +1092,11 @@ void MCSCF_ORB_GRAD::canonicalize_final(const std::shared_ptr& U) { C_->gemm(false, false, 1.0, C0_, U_, 0.0); - if (ints_->integral_type() == Custom) { - ints_->update_orbitals(C_, C_); - } else { - ints_->update_orbitals(C_, C_, false); - } + // update orbital object + orbitals_->set(C_, C_); + // update integrals object. Except for custom integrals, we update the integrals + // object with the new orbitals but do not update the integrals for efficiency + ints_->update_orbitals(orbitals_, ints_->integral_type() != IntegralType::Custom); build_mo_integrals(); } } // namespace forte diff --git a/forte/mcscf/mcscf_orb_grad.h b/forte/mcscf/mcscf_orb_grad.h index f168c0b3f..94ea787b2 100644 --- a/forte/mcscf/mcscf_orb_grad.h +++ b/forte/mcscf/mcscf_orb_grad.h @@ -59,8 +59,8 @@ class MCSCF_ORB_GRAD { * See J. Chem. Phys. 142, 224103 (2015) and Theor. Chem. Acc. 97, 88-95 (1997) */ MCSCF_ORB_GRAD(std::shared_ptr options, - std::shared_ptr mo_space_info, - std::shared_ptr ints, bool freeze_core); + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints, bool freeze_core); /// Evaluate the energy and orbital gradient double evaluate(std::shared_ptr x, std::shared_ptr g, @@ -103,6 +103,9 @@ class MCSCF_ORB_GRAD { /// The MOSpaceInfo object std::shared_ptr mo_space_info_; + /// The orbitals object + std::shared_ptr orbitals_; + /// The Forte integral std::shared_ptr ints_; diff --git a/forte/mrdsrg-spin-adapted/sa_dsrgpt.cc b/forte/mrdsrg-spin-adapted/sa_dsrgpt.cc index bc035d894..2a5cb1c51 100644 --- a/forte/mrdsrg-spin-adapted/sa_dsrgpt.cc +++ b/forte/mrdsrg-spin-adapted/sa_dsrgpt.cc @@ -57,7 +57,7 @@ void SA_DSRGPT::print_options() { {"Taylor expansion threshold", std::pow(10.0, -double(taylor_threshold_))}, {"Intruder amplitudes threshold", intruder_tamp_}}; - if (ints_->integral_type() == Cholesky) { + if (ints_->integral_type() == IntegralType::Cholesky) { auto cholesky_threshold = foptions_->get_double("CHOLESKY_TOLERANCE"); calculation_info_double.push_back({"Cholesky tolerance", cholesky_threshold}); } diff --git a/forte/mrdsrg-spin-adapted/sa_mrpt2.cc b/forte/mrdsrg-spin-adapted/sa_mrpt2.cc index 73ab94b21..423006c18 100644 --- a/forte/mrdsrg-spin-adapted/sa_mrpt2.cc +++ b/forte/mrdsrg-spin-adapted/sa_mrpt2.cc @@ -98,7 +98,7 @@ void SA_MRPT2::build_ints() { if (eri_df_) { std::vector blocks{"vvaa", "aacc", "avca", "avac", "vaaa", "aaca", "aaaa"}; V_ = BTF_->build(tensor_type_, "V", blocks); - if (ints_->integral_type() != DiskDF) { + if (ints_->integral_type() != IntegralType::DiskDF) { auto B = BTF_->build(tensor_type_, "B 3-idx", {"Lph"}); fill_three_index_ints(B); V_["abij"] = B["gai"] * B["gbj"]; @@ -228,7 +228,7 @@ void SA_MRPT2::check_memory() { auto mem_blocks = dsrg_mem_.compute_memory(blocks); dsrg_mem_.add_entry("2-electron (4-index) integrals", mem_blocks); dsrg_mem_.add_entry("T2 cluster amplitudes", 2 * mem_blocks); - if (ints_->integral_type() != DiskDF) { + if (ints_->integral_type() != IntegralType::DiskDF) { dsrg_mem_.add_entry("3-index auxiliary integrals", {"Lph"}); } } else { @@ -242,7 +242,7 @@ void SA_MRPT2::check_memory() { } // local memory for computing minimal V - if (ints_->integral_type() == DiskDF) { + if (ints_->integral_type() == IntegralType::DiskDF) { dsrg_mem_.add_entry("Local 3-index integrals", {"Lca", "Laa", "Lav"}, 1, false); } @@ -487,9 +487,11 @@ double SA_MRPT2::E_V_T2_CCVV() { auto i_nocc = i_batch_occ_mos.size(); ambit::Tensor Bi; // iaQ if (Bcv_file_exist) - Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, pqQ); + Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, + ThreeIntsBlockOrder::pqQ); else - Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, pqQ); + Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, + ThreeIntsBlockOrder::pqQ); auto& Bi_data = Bi.data(); for (size_t j_batch = i_batch, j_shift = i_shift; j_batch < nbatches; ++j_batch) { @@ -500,9 +502,11 @@ double SA_MRPT2::E_V_T2_CCVV() { Bj = Bi; } else { if (Bcv_file_exist) - Bj = read_Bcanonical("cv", {j_shift, j_shift + j_nocc}, {0, nv}, pqQ); + Bj = read_Bcanonical("cv", {j_shift, j_shift + j_nocc}, {0, nv}, + ThreeIntsBlockOrder::pqQ); else - Bj = ints_->three_integral_block(aux_mos_, j_batch_occ_mos, virt_mos_, pqQ); + Bj = ints_->three_integral_block(aux_mos_, j_batch_occ_mos, virt_mos_, + ThreeIntsBlockOrder::pqQ); } auto& Bj_data = Bj.data(); @@ -660,9 +664,9 @@ void SA_MRPT2::compute_Hbar1V_DF(ambit::Tensor& Hbar1, bool Vr) { // 3-index integrals (P|eu) ambit::Tensor Bu; // euQ if (!semi_v or !semi_a) - Bu = read_Bcanonical("va", {0, nv}, {0, na}, pqQ); + Bu = read_Bcanonical("va", {0, nv}, {0, na}, ThreeIntsBlockOrder::pqQ); else - Bu = ints_->three_integral_block(aux_mos_, virt_mos_, actv_mos_, pqQ); + Bu = ints_->three_integral_block(aux_mos_, virt_mos_, actv_mos_, ThreeIntsBlockOrder::pqQ); auto& Bu_data = Bu.data(); for (size_t i_batch = 0, i_shift = 0; i_batch < nbatches; ++i_batch) { @@ -670,9 +674,11 @@ void SA_MRPT2::compute_Hbar1V_DF(ambit::Tensor& Hbar1, bool Vr) { auto i_nocc = i_batch_occ_mos.size(); ambit::Tensor Bi; // iaQ if (!semi_c or !semi_v) - Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, pqQ); + Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, + ThreeIntsBlockOrder::pqQ); else - Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, pqQ); + Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, + ThreeIntsBlockOrder::pqQ); auto& Bi_data = Bi.data(); // index pairs of i and c @@ -836,9 +842,9 @@ void SA_MRPT2::compute_Hbar1C_DF(ambit::Tensor& Hbar1, bool Vr) { // 3-index integrals (P|mv) ambit::Tensor Bu; // umQ if (!semi_c or !semi_a) - Bu = read_Bcanonical("ac", {0, na}, {0, nc}, pqQ); + Bu = read_Bcanonical("ac", {0, na}, {0, nc}, ThreeIntsBlockOrder::pqQ); else - Bu = ints_->three_integral_block(aux_mos_, actv_mos_, core_mos_, pqQ); + Bu = ints_->three_integral_block(aux_mos_, actv_mos_, core_mos_, ThreeIntsBlockOrder::pqQ); auto& Bu_vec = Bu.data(); for (size_t c_batch = 0, c_shift = 0; c_batch < nbatches; ++c_batch) { @@ -846,9 +852,11 @@ void SA_MRPT2::compute_Hbar1C_DF(ambit::Tensor& Hbar1, bool Vr) { auto c_nvir = c_batch_vir_mos.size(); ambit::Tensor Bc; // aiQ if (!semi_c or !semi_v) - Bc = read_Bcanonical("vc", {c_shift, c_shift + c_nvir}, {0, nc}, pqQ); + Bc = read_Bcanonical("vc", {c_shift, c_shift + c_nvir}, {0, nc}, + ThreeIntsBlockOrder::pqQ); else - Bc = ints_->three_integral_block(aux_mos_, c_batch_vir_mos, core_mos_, pqQ); + Bc = ints_->three_integral_block(aux_mos_, c_batch_vir_mos, core_mos_, + ThreeIntsBlockOrder::pqQ); auto& Bc_vec = Bc.data(); #pragma omp parallel for num_threads(nthreads) diff --git a/forte/mrdsrg-spin-adapted/sa_mrpt2_oeprop.cc b/forte/mrdsrg-spin-adapted/sa_mrpt2_oeprop.cc index 5b9cd7500..7563bd479 100644 --- a/forte/mrdsrg-spin-adapted/sa_mrpt2_oeprop.cc +++ b/forte/mrdsrg-spin-adapted/sa_mrpt2_oeprop.cc @@ -268,9 +268,11 @@ void SA_MRPT2::compute_1rdm_cc_CCVV_DF(ambit::BlockedTensor& D1) { auto c_nvir = c_batch_vir_mos.size(); ambit::Tensor Bc; // ckQ if (Bcv_file_exist) - Bc = read_Bcanonical("vc", {c_shift, c_shift + c_nvir}, {0, nc}, pqQ); + Bc = read_Bcanonical("vc", {c_shift, c_shift + c_nvir}, {0, nc}, + ThreeIntsBlockOrder::pqQ); else - Bc = ints_->three_integral_block(aux_mos_, c_batch_vir_mos, core_mos_, pqQ); + Bc = ints_->three_integral_block(aux_mos_, c_batch_vir_mos, core_mos_, + ThreeIntsBlockOrder::pqQ); auto& Bc_data = Bc.data(); for (size_t d_batch = c_batch, d_shift = c_shift; d_batch < nbatches; ++d_batch) { @@ -282,9 +284,11 @@ void SA_MRPT2::compute_1rdm_cc_CCVV_DF(ambit::BlockedTensor& D1) { } else { ambit::Tensor Bd; // dlQ if (Bcv_file_exist) - Bd = read_Bcanonical("vc", {d_shift, d_shift + d_nvir}, {0, nc}, pqQ); + Bd = read_Bcanonical("vc", {d_shift, d_shift + d_nvir}, {0, nc}, + ThreeIntsBlockOrder::pqQ); else - Bd = ints_->three_integral_block(aux_mos_, d_batch_vir_mos, core_mos_, pqQ); + Bd = ints_->three_integral_block(aux_mos_, d_batch_vir_mos, core_mos_, + ThreeIntsBlockOrder::pqQ); } auto& Bd_data = Bd.data(); @@ -433,9 +437,11 @@ void SA_MRPT2::compute_1rdm_vv_CCVV_DF(ambit::BlockedTensor& D1) { auto i_nocc = i_batch_occ_mos.size(); ambit::Tensor Bi; // iaQ if (Bcv_file_exist) - Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, pqQ); + Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, + ThreeIntsBlockOrder::pqQ); else - Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, pqQ); + Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, + ThreeIntsBlockOrder::pqQ); auto& Bi_data = Bi.data(); for (size_t j_batch = i_batch, j_shift = i_shift; j_batch < nbatches; ++j_batch) { @@ -446,9 +452,11 @@ void SA_MRPT2::compute_1rdm_vv_CCVV_DF(ambit::BlockedTensor& D1) { Bj = Bi; } else { if (Bcv_file_exist) - Bj = read_Bcanonical("cv", {j_shift, j_shift + j_nocc}, {0, nv}, pqQ); + Bj = read_Bcanonical("cv", {j_shift, j_shift + j_nocc}, {0, nv}, + ThreeIntsBlockOrder::pqQ); else - Bj = ints_->three_integral_block(aux_mos_, j_batch_occ_mos, virt_mos_, pqQ); + Bj = ints_->three_integral_block(aux_mos_, j_batch_occ_mos, virt_mos_, + ThreeIntsBlockOrder::pqQ); } auto& Bj_data = Bj.data(); @@ -640,9 +648,9 @@ void SA_MRPT2::compute_1rdm_cc_CCAV_DF(ambit::BlockedTensor& D1, // 3-index integrals (P|mv) ambit::Tensor Bmv; // vmQ if (!semi_c or !semi_a) - Bmv = read_Bcanonical("ac", {0, na}, {0, nc}, pqQ); + Bmv = read_Bcanonical("ac", {0, na}, {0, nc}, ThreeIntsBlockOrder::pqQ); else - Bmv = ints_->three_integral_block(aux_mos_, actv_mos_, core_mos_, pqQ); + Bmv = ints_->three_integral_block(aux_mos_, actv_mos_, core_mos_, ThreeIntsBlockOrder::pqQ); auto& Bu_data = Bmv.data(); for (size_t c_batch = 0, c_shift = 0; c_batch < nbatches; ++c_batch) { @@ -650,9 +658,11 @@ void SA_MRPT2::compute_1rdm_cc_CCAV_DF(ambit::BlockedTensor& D1, auto c_nvir = c_batch_vir_mos.size(); ambit::Tensor Bc; // aiQ if (!semi_c or !semi_v) - Bc = read_Bcanonical("vc", {c_shift, c_shift + c_nvir}, {0, nc}, pqQ); + Bc = read_Bcanonical("vc", {c_shift, c_shift + c_nvir}, {0, nc}, + ThreeIntsBlockOrder::pqQ); else - Bc = ints_->three_integral_block(aux_mos_, c_batch_vir_mos, core_mos_, pqQ); + Bc = ints_->three_integral_block(aux_mos_, c_batch_vir_mos, core_mos_, + ThreeIntsBlockOrder::pqQ); auto& Bc_data = Bc.data(); #pragma omp parallel for num_threads(nthreads) @@ -851,9 +861,9 @@ void SA_MRPT2::compute_1rdm_aa_vv_CCAV_DF(ambit::BlockedTensor& D1, // 3-index integrals (P|mv) ambit::Tensor Bmv; // mvQ if (!semi_c or !semi_a) - Bmv = read_Bcanonical("ca", {0, nc}, {0, na}, pqQ); + Bmv = read_Bcanonical("ca", {0, nc}, {0, na}, ThreeIntsBlockOrder::pqQ); else - Bmv = ints_->three_integral_block(aux_mos_, core_mos_, actv_mos_, pqQ); + Bmv = ints_->three_integral_block(aux_mos_, core_mos_, actv_mos_, ThreeIntsBlockOrder::pqQ); auto& Bu_data = Bmv.data(); for (size_t i_batch = 0, i_shift = 0; i_batch < nbatches; ++i_batch) { @@ -861,9 +871,11 @@ void SA_MRPT2::compute_1rdm_aa_vv_CCAV_DF(ambit::BlockedTensor& D1, auto i_nocc = i_batch_occ_mos.size(); ambit::Tensor Bi; // iaQ if (!semi_c or !semi_v) - Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, pqQ); + Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, + ThreeIntsBlockOrder::pqQ); else - Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, pqQ); + Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, + ThreeIntsBlockOrder::pqQ); auto& Bi_data = Bi.data(); for (size_t j_batch = i_batch, j_shift = i_shift; j_batch < nbatches; ++j_batch) { @@ -874,9 +886,11 @@ void SA_MRPT2::compute_1rdm_aa_vv_CCAV_DF(ambit::BlockedTensor& D1, Bj = Bi; } else { if (!semi_c or !semi_v) - Bj = read_Bcanonical("cv", {j_shift, j_shift + j_nocc}, {0, nv}, pqQ); + Bj = read_Bcanonical("cv", {j_shift, j_shift + j_nocc}, {0, nv}, + ThreeIntsBlockOrder::pqQ); else - Bj = ints_->three_integral_block(aux_mos_, j_batch_occ_mos, virt_mos_, pqQ); + Bj = ints_->three_integral_block(aux_mos_, j_batch_occ_mos, virt_mos_, + ThreeIntsBlockOrder::pqQ); } auto& Bj_data = Bj.data(); @@ -1130,9 +1144,9 @@ void SA_MRPT2::compute_1rdm_cc_aa_CAVV_DF(ambit::BlockedTensor& D1, // 3-index integrals (P|mv) ambit::Tensor Beu; // euQ if (!semi_a or !semi_v) - Beu = read_Bcanonical("va", {0, nv}, {0, na}, pqQ); + Beu = read_Bcanonical("va", {0, nv}, {0, na}, ThreeIntsBlockOrder::pqQ); else - Beu = ints_->three_integral_block(aux_mos_, virt_mos_, actv_mos_, pqQ); + Beu = ints_->three_integral_block(aux_mos_, virt_mos_, actv_mos_, ThreeIntsBlockOrder::pqQ); auto& Bu_data = Beu.data(); for (size_t c_batch = 0, c_shift = 0; c_batch < nbatches; ++c_batch) { @@ -1140,9 +1154,11 @@ void SA_MRPT2::compute_1rdm_cc_aa_CAVV_DF(ambit::BlockedTensor& D1, auto c_nvir = c_batch_vir_mos.size(); ambit::Tensor Bc; // ckQ if (!semi_c or !semi_v) - Bc = read_Bcanonical("vc", {c_shift, c_shift + c_nvir}, {0, nc}, pqQ); + Bc = read_Bcanonical("vc", {c_shift, c_shift + c_nvir}, {0, nc}, + ThreeIntsBlockOrder::pqQ); else - Bc = ints_->three_integral_block(aux_mos_, c_batch_vir_mos, core_mos_, pqQ); + Bc = ints_->three_integral_block(aux_mos_, c_batch_vir_mos, core_mos_, + ThreeIntsBlockOrder::pqQ); auto& Bc_data = Bc.data(); for (size_t d_batch = c_batch, d_shift = c_shift; d_batch < nbatches; ++d_batch) { @@ -1153,9 +1169,11 @@ void SA_MRPT2::compute_1rdm_cc_aa_CAVV_DF(ambit::BlockedTensor& D1, Bd = Bc; } else { if (!semi_c or !semi_v) - Bd = read_Bcanonical("vc", {d_shift, d_shift + d_nvir}, {0, nc}, pqQ); + Bd = read_Bcanonical("vc", {d_shift, d_shift + d_nvir}, {0, nc}, + ThreeIntsBlockOrder::pqQ); else - Bd = ints_->three_integral_block(aux_mos_, d_batch_vir_mos, core_mos_, pqQ); + Bd = ints_->three_integral_block(aux_mos_, d_batch_vir_mos, core_mos_, + ThreeIntsBlockOrder::pqQ); } auto& Bd_data = Bd.data(); @@ -1394,9 +1412,9 @@ void SA_MRPT2::compute_1rdm_vv_CAVV_DF(ambit::BlockedTensor& D1, // 3-index integrals (P|eu) ambit::Tensor Beu; // euQ if (!semi_a or !semi_v) - Beu = read_Bcanonical("va", {0, nv}, {0, na}, pqQ); + Beu = read_Bcanonical("va", {0, nv}, {0, na}, ThreeIntsBlockOrder::pqQ); else - Beu = ints_->three_integral_block(aux_mos_, virt_mos_, actv_mos_, pqQ); + Beu = ints_->three_integral_block(aux_mos_, virt_mos_, actv_mos_, ThreeIntsBlockOrder::pqQ); auto& Bu_data = Beu.data(); for (size_t i_batch = 0, i_shift = 0; i_batch < nbatches; ++i_batch) { @@ -1404,9 +1422,11 @@ void SA_MRPT2::compute_1rdm_vv_CAVV_DF(ambit::BlockedTensor& D1, auto i_nocc = i_batch_occ_mos.size(); ambit::Tensor Bi; // iaQ if (!semi_c or !semi_v) - Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, pqQ); + Bi = read_Bcanonical("cv", {i_shift, i_shift + i_nocc}, {0, nv}, + ThreeIntsBlockOrder::pqQ); else - Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, pqQ); + Bi = ints_->three_integral_block(aux_mos_, i_batch_occ_mos, virt_mos_, + ThreeIntsBlockOrder::pqQ); auto& Bi_data = Bi.data(); // index pairs of i and c diff --git a/forte/mrdsrg-spin-adapted/sadsrg.cc b/forte/mrdsrg-spin-adapted/sadsrg.cc index 84a384289..1450e43c4 100644 --- a/forte/mrdsrg-spin-adapted/sadsrg.cc +++ b/forte/mrdsrg-spin-adapted/sadsrg.cc @@ -279,7 +279,8 @@ void SADSRG::set_ambit_MOSpace() { void SADSRG::check_init_memory() { mem_sys_ = psi::Process::environment.get_memory(); int64_t mem_left = mem_sys_ * 0.9; - if (ints_->integral_type() != DiskDF and ints_->integral_type() != Cholesky) { + if (ints_->integral_type() != IntegralType::DiskDF and + ints_->integral_type() != IntegralType::Cholesky) { mem_left -= ints_->jk()->memory_estimate() * sizeof(double); } @@ -287,7 +288,7 @@ void SADSRG::check_init_memory() { size_t n_ele = 0; auto ng = mo_space_info_->size("CORRELATED"); if (eri_df_) { - if (ints_->integral_type() != DiskDF) { + if (ints_->integral_type() != IntegralType::DiskDF) { auto nQ = aux_mos_.size(); n_ele = nQ * ng * ng; } @@ -1021,9 +1022,9 @@ ambit::Tensor SADSRG::read_Bcanonical(const std::string& block, * depending on the contiguity of the fastest indices. Otherwise, this function tries to * load the transpose block {block[1], block[0]} before it throws an error. * - * If the auxiliary index is desired to be the last (i.e., pqQ), we will still load data - * in the order of Qpq. Then, an in-place matrix transpose is performed to give the correct - * ordering of data storage in memory. + * If the auxiliary index is desired to be the last (i.e., ThreeIntsBlockOrder:pqQ), we will + * still load data in the order of Qpq. Then, an in-place matrix transpose is performed to give + * the correct ordering of data storage in memory. */ if (!eri_df_) throw std::runtime_error("For DF/CD integrals ONLY!"); @@ -1044,7 +1045,7 @@ ambit::Tensor SADSRG::read_Bcanonical(const std::string& block, throw std::runtime_error("Incorrect MO indices! Check mos1_range and mos2_range!"); ambit::Tensor T; - if (order == pqQ) + if (order == ThreeIntsBlockOrder::pqQ) T = ambit::Tensor::build(tensor_type_, "Bcan_" + block, {s1, s2, nQ}); else T = ambit::Tensor::build(tensor_type_, "Bcan_" + block, {nQ, s1, s2}); @@ -1200,7 +1201,7 @@ ambit::Tensor SADSRG::read_Bcanonical(const std::string& block, } // in-place matrix transpose - if (order == pqQ) { + if (order == ThreeIntsBlockOrder::pqQ) { matrix_transpose_in_place(Tdata, nQ, S); } diff --git a/forte/mrdsrg-spin-adapted/sadsrg.h b/forte/mrdsrg-spin-adapted/sadsrg.h index 6f1cf8807..2a399b904 100644 --- a/forte/mrdsrg-spin-adapted/sadsrg.h +++ b/forte/mrdsrg-spin-adapted/sadsrg.h @@ -70,7 +70,7 @@ class SADSRG : public DynamicCorrelationSolver { void set_Uactv(ambit::Tensor& U); /// If the amplitudes are converged or not - bool converged() {return converged_; } + bool converged() { return converged_; } protected: /// Startup function called in constructor @@ -308,7 +308,7 @@ class SADSRG : public DynamicCorrelationSolver { ambit::Tensor read_Bcanonical(const std::string& block, const std::pair& mos1_range, const std::pair& mos2_range, - ThreeIntsBlockOrder order = Qpq); + ThreeIntsBlockOrder order = ThreeIntsBlockOrder::Qpq); /// File names for canonicalized 3-index integrals std::unordered_map Bcan_files_; 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 index e0cac5c59..926b974b0 100644 --- 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 @@ -6,15 +6,19 @@ #include "psi4/libmints/matrix.h" #include "psi4/libmints/wavefunction.h" #include "psi4/psifiles.h" -#include "../dsrg_mrpt2.h" -#include "helpers/timer.h" -#include "gradient_tpdm/backtransform_tpdm.h" #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 "helpers/timer.h" + +#include "base_classes/mo_space_info.h" +#include "base_classes/orbitals.h" + +#include "../dsrg_mrpt2.h" +#include "gradient_tpdm/backtransform_tpdm.h" + using namespace ambit; using namespace psi; @@ -52,7 +56,7 @@ void DSRG_MRPT2::write_lagrangian() { }); } - L->back_transform(ints_->Ca()->clone()); + L->back_transform(ints_->orbitals()->Ca()->clone()); ints_->wfn()->set_lagrangian( std::make_shared("Lagrangian", nirrep, irrep_vec, irrep_vec)); ints_->wfn()->lagrangian()->copy(L); @@ -103,7 +107,7 @@ void DSRG_MRPT2::write_1rdm_spin_dependent() { } }); } - D1->back_transform(ints_->Ca()->clone()); + D1->back_transform(ints_->orbitals()->Ca()->clone()); ints_->wfn()->Da()->copy(D1); ints_->wfn()->Db()->copy(D1); outfile->Printf("Done"); @@ -697,7 +701,7 @@ void DSRG_MRPT2::write_df_rdm() { offset_col += nmopi; } - auto Ca = ints_->Ca(); + auto Ca = ints_->orbitals()->Ca(); // Copy Ca to a matrix without symmetry blocking auto Cat = std::make_shared("Ca temp matrix", nso, nmo); diff --git a/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.cc b/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.cc index d6db5dad5..95abc2fb1 100644 --- a/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.cc +++ b/forte/mrdsrg-spin-integrated/three_dsrg_mrpt2.cc @@ -55,13 +55,17 @@ #include "psi4/libpsio/psio.hpp" #include "psi4/libqt/qt.h" -#include "orbital-helpers/ao_helper.h" +#include "base_classes/orbitals.h" + #include "helpers/blockedtensorfactory.h" #include "helpers/printing.h" #include "helpers/timer.h" +#include "orbital-helpers/ao_helper.h" + #include "fci/fci_solver.h" #include "fci/fci_vector.h" #include "sci/aci.h" + #include "three_dsrg_mrpt2.h" using namespace ambit; @@ -254,7 +258,7 @@ void THREE_DSRG_MRPT2::startup() { T1_ = BTF_->build(tensor_type_, "T1 Amplitudes", spin_cases({"hp"})); // allocate memory for T2 and V - if (integral_type_ != DiskDF) { + if (integral_type_ != IntegralType::DiskDF) { std::vector list_of_pphh_V = BTF_->generate_indices("vac", "pphh"); V_ = BTF_->build(tensor_type_, "V_", BTF_->spin_cases_avoid(list_of_pphh_V, 1)); T2_ = BTF_->build(tensor_type_, "T2 Amplitudes", BTF_->spin_cases_avoid(no_hhpp_, 1)); @@ -361,7 +365,7 @@ double THREE_DSRG_MRPT2::compute_energy() { outfile->Printf("\n Reference Energy = %.15f", Eref_); // compute T2 and renormalize V - if (integral_type_ != DiskDF) { + if (integral_type_ != IntegralType::DiskDF) { compute_t2(); renormalize_V(); } else { @@ -1161,7 +1165,7 @@ double THREE_DSRG_MRPT2::E_VT2_2() { if (my_proc == 0) { outfile->Printf("\n %-40s ...", "Computing <[V, T2]> (C_2)^4 (no ccvv)"); // TODO: Implement these without storing V and/or T2 by using blocking - if (integral_type_ != DiskDF) { + if (integral_type_ != IntegralType::DiskDF) { temp.zero(); temp["vu"] += 0.5 * V_["efmu"] * T2_["mvef"]; temp["vu"] += V_["fEuM"] * T2_["vMfE"]; @@ -1795,7 +1799,7 @@ double THREE_DSRG_MRPT2::E_VT2_2_ambit() { /// variable. Requires the reading from aptei_block which makes code /// general for all, but makes it slow for DiskDF. - if (integral_type_ == DiskDF) { + if (integral_type_ == IntegralType::DiskDF) { std::vector BefVec; std::vector BefJKVec; std::vector RDVec; @@ -2282,7 +2286,7 @@ double THREE_DSRG_MRPT2::E_VT2_2_AO_Slow() { double Ealpha = 0.0; double Emixed = 0.0; double Ebeta = 0.0; - auto Cwfn = ints_->Ca()->clone(); + auto Cwfn = ints_->orbitals()->Ca()->clone(); if (mo_space_info_->nirrep() != 1) throw psi::PSIEXCEPTION("AO-DSRGMPT2 does not work with symmetry"); @@ -3007,7 +3011,7 @@ void THREE_DSRG_MRPT2::form_Hbar() { outfile->Printf("Done. Timing: %10.3f s.", timer2.get()); // Again, the below block assume V[ijab] = V[iJaB] - V[iJbA] - if (integral_type_ == DiskDF) { + if (integral_type_ == IntegralType::DiskDF) { C1.zero(); local_timer timer3; @@ -3409,7 +3413,7 @@ void THREE_DSRG_MRPT2::compute_Heff_2nd_coupling(double& H0, ambit::Tensor& H1a, // reset APTEI because it is renormalized std::vector list_of_pphh_V = BTF_->generate_indices("vac", "pphh"); - if (integral_type_ != DiskDF) { + if (integral_type_ != IntegralType::DiskDF) { V_["abij"] = ThreeIntegral_["gai"] * ThreeIntegral_["gbj"]; V_["abij"] -= ThreeIntegral_["gaj"] * ThreeIntegral_["gbi"]; @@ -3501,7 +3505,7 @@ void THREE_DSRG_MRPT2::compute_Heff_2nd_coupling(double& H0, ambit::Tensor& H1a, H1["VU"] -= 0.5 * V_["AVMN"] * T2_["MNAU"]; H1["VU"] -= V_["aVmN"] * T2_["mNaU"]; - if (integral_type_ == DiskDF) { + if (integral_type_ == IntegralType::DiskDF) { compute_Hbar1C_diskDF(H1, false); compute_Hbar1V_diskDF(H1, false); } diff --git a/forte/orbital-helpers/ci-no/ci-no.cc b/forte/orbital-helpers/ci-no/ci-no.cc index c0bfc382f..28b71e903 100644 --- a/forte/orbital-helpers/ci-no/ci-no.cc +++ b/forte/orbital-helpers/ci-no/ci-no.cc @@ -62,9 +62,9 @@ std::string dimension_to_string(psi::Dimension dim) { return s; } -CINO::CINO(std::shared_ptr options, std::shared_ptr ints, - std::shared_ptr mo_space_info) - : OrbitalTransform(ints, mo_space_info), options_(options) { +CINO::CINO(std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, std::shared_ptr ints) + : OrbitalTransform(mo_space_info, orbitals, ints), options_(options) { fci_ints_ = std::make_shared( ints, mo_space_info_->corr_absolute_mo("ACTIVE"), mo_space_info_->symmetry("ACTIVE"), mo_space_info_->corr_absolute_mo("RESTRICTED_DOCC")); diff --git a/forte/orbital-helpers/ci-no/ci-no.h b/forte/orbital-helpers/ci-no/ci-no.h index 6d77c34d1..4f08aacf1 100644 --- a/forte/orbital-helpers/ci-no/ci-no.h +++ b/forte/orbital-helpers/ci-no/ci-no.h @@ -58,8 +58,8 @@ class CINO : public OrbitalTransform { * @param ints A pointer to an allocated integral object * @param mo_space_info A pointer to the MOSpaceInfo object */ - CINO(std::shared_ptr options, std::shared_ptr ints, - std::shared_ptr mo_space_info); + CINO(std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, std::shared_ptr ints); /// Destructor ~CINO(); diff --git a/forte/orbital-helpers/ci-no/mrci-no.cc b/forte/orbital-helpers/ci-no/mrci-no.cc index 01782f614..b3a885e4c 100644 --- a/forte/orbital-helpers/ci-no/mrci-no.cc +++ b/forte/orbital-helpers/ci-no/mrci-no.cc @@ -71,8 +71,9 @@ namespace forte { //} MRCINO::MRCINO(std::shared_ptr scf_info, std::shared_ptr options, - std::shared_ptr ints, std::shared_ptr mo_space_info) - : OrbitalTransform(ints, mo_space_info), scf_info_(scf_info), options_(options) { + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints) + : OrbitalTransform(mo_space_info, orbitals, ints), scf_info_(scf_info), options_(options) { // Copy the wavefunction information std::vector active_mo(mo_space_info_->size("CORRELATED")); diff --git a/forte/orbital-helpers/ci-no/mrci-no.h b/forte/orbital-helpers/ci-no/mrci-no.h index 92c3bed25..d81d000ca 100644 --- a/forte/orbital-helpers/ci-no/mrci-no.h +++ b/forte/orbital-helpers/ci-no/mrci-no.h @@ -42,7 +42,8 @@ class MRCINO : public OrbitalTransform { // ==> Class Constructor and Destructor <== MRCINO(std::shared_ptr scf_info, std::shared_ptr options, - std::shared_ptr ints, std::shared_ptr mo_space_info); + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints); /// Destructor ~MRCINO(); diff --git a/forte/orbital-helpers/localize.cc b/forte/orbital-helpers/localize.cc index 4371be497..1a0df846f 100644 --- a/forte/orbital-helpers/localize.cc +++ b/forte/orbital-helpers/localize.cc @@ -35,6 +35,7 @@ #include "helpers/printing.h" #include "base_classes/rdms.h" +#include "base_classes/orbitals.h" #include "localize.h" @@ -42,9 +43,10 @@ using namespace psi; namespace forte { -Localize::Localize(std::shared_ptr options, std::shared_ptr ints, - std::shared_ptr mo_space_info) - : OrbitalTransform(ints, mo_space_info) { +Localize::Localize(std::shared_ptr options, + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints) + : OrbitalTransform(mo_space_info, orbitals, ints) { if (ints_->nirrep() > 1) { throw psi::PSIEXCEPTION("\n\n ERROR: Localizer only implemented for C1 symmetry!"); @@ -81,7 +83,7 @@ void Localize::compute_transformation() { exit(1); } - auto Ca = ints_->Ca()->clone(); + auto Ca = orbitals_->Ca()->clone(); Ua_ = std::make_shared("U", Ca->coldim(), Ca->coldim()); Ub_ = std::make_shared("U", Ca->coldim(), Ca->coldim()); diff --git a/forte/orbital-helpers/localize.h b/forte/orbital-helpers/localize.h index 0fe50ca01..2c62e4d81 100644 --- a/forte/orbital-helpers/localize.h +++ b/forte/orbital-helpers/localize.h @@ -45,8 +45,8 @@ namespace forte { class Localize : public OrbitalTransform { public: - Localize(std::shared_ptr options, std::shared_ptr ints, - std::shared_ptr mo_space_info); + Localize(std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, std::shared_ptr ints); // Compute the rotation matrices void compute_transformation() override; diff --git a/forte/orbital-helpers/mp2_nos.cc b/forte/orbital-helpers/mp2_nos.cc index 329629a37..96b91ba33 100644 --- a/forte/orbital-helpers/mp2_nos.cc +++ b/forte/orbital-helpers/mp2_nos.cc @@ -61,8 +61,9 @@ using namespace psi; namespace forte { MP2_NOS::MP2_NOS(std::shared_ptr scf_info, std::shared_ptr options, - std::shared_ptr ints, std::shared_ptr mo_space_info) - : OrbitalTransform(ints, mo_space_info), scf_info_(scf_info), options_(options) {} + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints) + : OrbitalTransform(mo_space_info, orbitals, ints), scf_info_(scf_info), options_(options) {} void MP2_NOS::compute_transformation() { print_method_banner( diff --git a/forte/orbital-helpers/mp2_nos.h b/forte/orbital-helpers/mp2_nos.h index 15d8dc591..ff9ae9d50 100644 --- a/forte/orbital-helpers/mp2_nos.h +++ b/forte/orbital-helpers/mp2_nos.h @@ -45,7 +45,8 @@ class MP2_NOS : public OrbitalTransform { public: // => Constructor <= // MP2_NOS(std::shared_ptr scf_info, std::shared_ptr options, - std::shared_ptr ints, std::shared_ptr mo_space_info); + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints); void compute_transformation() override; diff --git a/forte/orbital-helpers/mrpt2_nos.cc b/forte/orbital-helpers/mrpt2_nos.cc index b6d317671..f842bdadb 100644 --- a/forte/orbital-helpers/mrpt2_nos.cc +++ b/forte/orbital-helpers/mrpt2_nos.cc @@ -46,10 +46,10 @@ using namespace psi; namespace forte { -MRPT2_NOS::MRPT2_NOS(std::shared_ptr rdms, std::shared_ptr scf_info, - std::shared_ptr options, std::shared_ptr ints, - std::shared_ptr mo_space_info) - : OrbitalTransform(ints, mo_space_info), options_(options) { +MRPT2_NOS::MRPT2_NOS(std::shared_ptr scf_info, std::shared_ptr options, + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints, std::shared_ptr rdms) + : OrbitalTransform(mo_space_info, orbitals, ints), options_(options) { mrpt2_ = std::make_shared(rdms, scf_info, options, ints, mo_space_info); } diff --git a/forte/orbital-helpers/mrpt2_nos.h b/forte/orbital-helpers/mrpt2_nos.h index 24581d367..e9990957f 100644 --- a/forte/orbital-helpers/mrpt2_nos.h +++ b/forte/orbital-helpers/mrpt2_nos.h @@ -42,9 +42,9 @@ namespace forte { class MRPT2_NOS : public OrbitalTransform { public: // => Constructor <= // - MRPT2_NOS(std::shared_ptr rdms, std::shared_ptr scf_info, - std::shared_ptr options, std::shared_ptr ints, - std::shared_ptr mo_space_info); + MRPT2_NOS(std::shared_ptr scf_info, std::shared_ptr options, + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints, std::shared_ptr rdms); void compute_transformation(); diff --git a/forte/orbital-helpers/semi_canonicalize.cc b/forte/orbital-helpers/semi_canonicalize.cc index 48c2d5ed6..5ab778603 100644 --- a/forte/orbital-helpers/semi_canonicalize.cc +++ b/forte/orbital-helpers/semi_canonicalize.cc @@ -64,7 +64,7 @@ void SemiCanonical::read_options(const std::shared_ptr& options) { // compute thresholds from options double econv = options->get_double("E_CONVERGENCE"); threshold_tight_ = (econv < 1.0e-12) ? 1.0e-12 : econv; - if (ints_->integral_type() == Cholesky) { + if (ints_->integral_type() == IntegralType::Cholesky) { double cd_tlr = options->get_double("CHOLESKY_TOLERANCE"); threshold_tight_ = (threshold_tight_ < 0.5 * cd_tlr) ? 0.5 * cd_tlr : threshold_tight_; } diff --git a/forte/orbital-helpers/unpaired_density.cc b/forte/orbital-helpers/unpaired_density.cc index b8a26c8c6..3285e070c 100644 --- a/forte/orbital-helpers/unpaired_density.cc +++ b/forte/orbital-helpers/unpaired_density.cc @@ -26,6 +26,14 @@ * @END LICENSE */ +#include "psi4/libmints/wavefunction.h" +#include "psi4/libmints/local.h" + +#include "base_classes/mo_space_info.h" +#include "integrals/integrals.h" +#include "base_classes/rdms.h" +#include "iao_builder.h" + #include "psi4/libmints/molecule.h" #include "psi4/libpsio/psio.h" #include "psi4/libpsio/psio.hpp" @@ -38,6 +46,8 @@ #include "psi4/libmints/vector.h" #include "base_classes/forte_options.h" +#include "base_classes/orbitals.h" + #include "unpaired_density.h" #include "localize.h" @@ -45,11 +55,12 @@ using namespace psi; namespace forte { -UPDensity::UPDensity(std::shared_ptr ints, - std::shared_ptr mo_space_info, - std::shared_ptr options, std::shared_ptr Ua, +UPDensity::UPDensity(std::shared_ptr options, + std::shared_ptr mo_space_info, std::shared_ptr orbitals, + std::shared_ptr ints, std::shared_ptr Ua, std::shared_ptr Ub) - : options_(options), ints_(ints), mo_space_info_(mo_space_info), Uas_(Ua), Ubs_(Ub) {} + : options_(options), mo_space_info_(mo_space_info), orbitals_(orbitals), ints_(ints), Uas_(Ua), + Ubs_(Ub) {} void UPDensity::compute_unpaired_density(std::vector& oprdm_a, std::vector& oprdm_b) { @@ -171,7 +182,7 @@ void UPDensity::compute_unpaired_density(std::vector& oprdm_a, psi::Dimension nactpi = mo_space_info_->dimension("ACTIVE"); psi::Dimension nmopi = mo_space_info_->dimension("ALL"); psi::Dimension ncmopi = mo_space_info_->dimension("CORRELATED"); - size_t nirrep = ints_->nirrep(); + size_t nirrep = mo_space_info_->nirrep(); psi::Dimension rdocc = mo_space_info_->dimension("RESTRICTED_DOCC"); psi::Dimension fdocc = mo_space_info_->dimension("FROZEN_DOCC"); @@ -234,7 +245,7 @@ void UPDensity::compute_unpaired_density(std::vector& oprdm_a, // relocalize to atoms // Grab matrix that takes the transforms from the NO basis to our local basis - auto loc = std::make_shared(options_, ints_, mo_space_info_); + auto loc = std::make_shared(options_, mo_space_info_, orbitals_, ints_); std::vector actmo = mo_space_info_->absolute_mo("ACTIVE"); std::vector loc_mo(2); @@ -270,8 +281,8 @@ void UPDensity::compute_unpaired_density(std::vector& oprdm_a, // Build the density using scaled columns of C - auto Ca = ints_->Ca(); - auto Cb = ints_->Cb(); + auto Ca = orbitals_->Ca(); + auto Cb = orbitals_->Cb(); auto Ca_new = psi::linalg::doublet(Ca->clone(), Ua, false, false); auto Cb_new = psi::linalg::doublet(Cb->clone(), Ub, false, false); diff --git a/forte/orbital-helpers/unpaired_density.h b/forte/orbital-helpers/unpaired_density.h index 8426f109e..f028e3462 100644 --- a/forte/orbital-helpers/unpaired_density.h +++ b/forte/orbital-helpers/unpaired_density.h @@ -28,21 +28,20 @@ #pragma once -#include "psi4/libmints/wavefunction.h" -#include "psi4/libmints/local.h" - -#include "base_classes/mo_space_info.h" -#include "integrals/integrals.h" -#include "base_classes/rdms.h" -#include "iao_builder.h" - +namespace psi { +class Matrix; +} // namespace psi namespace forte { +class MOSpaceInfo; +class ForteOptions; +class Orbitals; + class UPDensity { public: - UPDensity(std::shared_ptr ints, std::shared_ptr mo_space_info, - std::shared_ptr options, std::shared_ptr Ua, - std::shared_ptr Ub); + UPDensity(std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, std::shared_ptr ints, + std::shared_ptr Ua, std::shared_ptr Ub); ~UPDensity(); @@ -50,8 +49,9 @@ class UPDensity { private: std::shared_ptr options_; - std::shared_ptr ints_; std::shared_ptr mo_space_info_; + std::shared_ptr orbitals_; + std::shared_ptr ints_; std::shared_ptr Uas_; std::shared_ptr Ubs_; }; diff --git a/forte/post_process/spin_corr.cc b/forte/post_process/spin_corr.cc index 787c0c4e0..c73153c9a 100644 --- a/forte/post_process/spin_corr.cc +++ b/forte/post_process/spin_corr.cc @@ -32,8 +32,11 @@ #include "psi4/libpsi4util/process.h" #include "base_classes/forte_options.h" +#include "base_classes/orbitals.h" + #include "helpers/printing.h" #include "helpers/helpers.h" + #include "post_process/spin_corr.h" using namespace psi; @@ -41,9 +44,10 @@ using namespace psi; namespace forte { SpinCorr::SpinCorr(std::shared_ptr rdms, std::shared_ptr options, - std::shared_ptr mo_space_info, + std::shared_ptr mo_space_info, std::shared_ptr orbitals, std::shared_ptr as_ints) - : rdms_(rdms), options_(options), mo_space_info_(mo_space_info), as_ints_(as_ints) { + : rdms_(rdms), options_(options), mo_space_info_(mo_space_info), orbitals_(orbitals), + as_ints_(as_ints) { nactpi_ = mo_space_info_->dimension("ACTIVE"); nirrep_ = nactpi_.n(); @@ -191,18 +195,21 @@ void SpinCorr::spin_analysis() { } } - auto CA = as_ints_->ints()->Ca()->clone(); - auto CB = as_ints_->ints()->Cb()->clone(); + auto CA = orbitals_->Ca()->clone(); + auto CB = orbitals_->Cb()->clone(); auto Ca_new = psi::linalg::doublet(CA, Ua_full, false, false); auto Cb_new = psi::linalg::doublet(CB, Ub_full, false, false); - as_ints_->ints()->update_orbitals(Ca_new, Cb_new, false); + orbitals_->set(Ca_new, Cb_new); + + as_ints_->ints()->update_orbitals(orbitals_, false); } else if (options_->get_str("SPIN_BASIS") == "LOCAL") { outfile->Printf("\n Computing spin correlation in local basis \n"); - auto loc = std::make_shared(options_, as_ints_->ints(), mo_space_info_); + auto loc = + std::make_shared(options_, mo_space_info_, orbitals_, as_ints_->ints()); std::vector actmo = mo_space_info_->absolute_mo("ACTIVE"); std::vector loc_mo(2); @@ -363,8 +370,9 @@ void SpinCorr::spin_analysis() { void perform_spin_analysis(std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, std::shared_ptr as_ints) { - SpinCorr spin(rdms, options, mo_space_info, as_ints); + SpinCorr spin(rdms, options, mo_space_info, orbitals, as_ints); spin.spin_analysis(); } diff --git a/forte/post_process/spin_corr.h b/forte/post_process/spin_corr.h index b5d8be573..45f3ed726 100644 --- a/forte/post_process/spin_corr.h +++ b/forte/post_process/spin_corr.h @@ -48,7 +48,7 @@ class SpinCorr { public: SpinCorr(std::shared_ptr rdms, std::shared_ptr options, - std::shared_ptr mo_space_info, + std::shared_ptr mo_space_info, std::shared_ptr orbitals, std::shared_ptr as_ints); std::pair, std::shared_ptr> compute_nos(); @@ -62,6 +62,8 @@ class SpinCorr { std::shared_ptr mo_space_info_; + std::shared_ptr orbitals_; + std::shared_ptr as_ints_; size_t nact_; From 9926b4296b60b353553687207b57a2da36a85016 Mon Sep 17 00:00:00 2001 From: fevangelista Date: Fri, 17 May 2024 23:10:03 -0400 Subject: [PATCH 2/4] Code compiles --- forte/api/forte_python_module.cc | 3 +++ forte/api/orbital_api.cc | 13 +++++++++++++ forte/base_classes/orbital_transform.h | 9 ++++----- forte/base_classes/orbitals.cc | 10 +++++----- forte/base_classes/orbitals.h | 8 ++++---- forte/integrals/integrals_psi4_interface.cc | 10 +++++----- forte/modules/mcscf.py | 2 +- forte/modules/objects_factory_psi4.py | 11 ++++++++--- forte/post_process/spin_corr.h | 1 + 9 files changed, 44 insertions(+), 23 deletions(-) diff --git a/forte/api/forte_python_module.cc b/forte/api/forte_python_module.cc index 14519acb9..f405ad05a 100644 --- a/forte/api/forte_python_module.cc +++ b/forte/api/forte_python_module.cc @@ -85,6 +85,7 @@ void export_StateInfo(py::module& m); void export_SigmaVector(py::module& m); void export_SparseCISolver(py::module& m); void export_ForteCubeFile(py::module& m); +void export_Orbitals(py::module& m); void export_OrbitalTransform(py::module& m); void export_Localize(py::module& m); void export_SemiCanonical(py::module& m); @@ -292,6 +293,8 @@ PYBIND11_MODULE(_forte, m) { export_ForteIntegrals(m); export_Symmetry(m); + export_Orbitals(m); + m.def("make_orbitals_from_psi", &make_orbitals_from_psi); export_OrbitalTransform(m); export_Localize(m); export_SemiCanonical(m); diff --git a/forte/api/orbital_api.cc b/forte/api/orbital_api.cc index a0a891a9b..27d2bbc18 100644 --- a/forte/api/orbital_api.cc +++ b/forte/api/orbital_api.cc @@ -46,6 +46,19 @@ using namespace pybind11::literals; namespace forte { +void export_Orbitals(py::module& m) { + py::class_>(m, "Orbitals") + .def(py::init&, const std::shared_ptr&>()) + .def("Ca", &Orbitals::Ca, "Return the alpha orbital coefficient matrix") + .def("Cb", &Orbitals::Cb, "Return the beta orbital coefficient matrix") + .def("set", &Orbitals::set, "Set the alpha and beta orbital coefficient matrices") + .def("rotate", &Orbitals::rotate, + "Rotate the orbitals using the given transformation matrices") + .def("copy", &Orbitals::copy, "Copy the orbitals") + .def("are_spin_restricted", &Orbitals::are_spin_restricted, + "Check if the orbitals are spin restricted"); +} + /// Export the OrbitalTransform class void export_OrbitalTransform(py::module& m) { py::class_(m, "OrbitalTransform") diff --git a/forte/base_classes/orbital_transform.h b/forte/base_classes/orbital_transform.h index 680c3fba5..6d6584063 100644 --- a/forte/base_classes/orbital_transform.h +++ b/forte/base_classes/orbital_transform.h @@ -45,10 +45,9 @@ class OrbitalTransform { std::shared_ptr Ub_; }; -std::unique_ptr -make_orbital_transformation(const std::string& type, std::shared_ptr scf_info, - std::shared_ptr options, - std::shared_ptr ints, - std::shared_ptr mo_space_info); +std::unique_ptr make_orbital_transformation( + const std::string& type, std::shared_ptr scf_info, + std::shared_ptr options, std::shared_ptr ints, + std::shared_ptr mo_space_info, std::shared_ptr orbitals); } // namespace forte diff --git a/forte/base_classes/orbitals.cc b/forte/base_classes/orbitals.cc index aa596ed87..4a7cb342e 100644 --- a/forte/base_classes/orbitals.cc +++ b/forte/base_classes/orbitals.cc @@ -38,9 +38,9 @@ Orbitals::Orbitals(const std::shared_ptr& Ca, const std::shared_ptr Cb_ = Cb->clone(); } -const std::shared_ptr Orbitals::Ca() const { return Ca_; } +std::shared_ptr Orbitals::Ca() { return Ca_; } -const std::shared_ptr Orbitals::Cb() const { return Cb_; } +std::shared_ptr Orbitals::Cb() { return Cb_; } void Orbitals::set(const std::shared_ptr& Ca, const std::shared_ptr& Cb) { if (not(elementwise_compatible_matrices(Ca_, Ca) and @@ -69,9 +69,9 @@ bool Orbitals::are_spin_restricted(double threshold) const { return matrix_distance(Ca_, Cb_) < threshold; } -std::unique_ptr -make_orbitals_from_psi(const std::shared_ptr& wfn, bool restricted) { - return std::make_unique(wfn->Ca(), restricted ? wfn->Ca() : wfn->Cb()); +std::shared_ptr make_orbitals_from_psi(const std::shared_ptr& wfn, + bool restricted) { + return std::make_shared(wfn->Ca(), restricted ? wfn->Ca() : wfn->Cb()); } } // namespace forte \ No newline at end of file diff --git a/forte/base_classes/orbitals.h b/forte/base_classes/orbitals.h index 671dc6ef4..2e8f17c0e 100644 --- a/forte/base_classes/orbitals.h +++ b/forte/base_classes/orbitals.h @@ -46,9 +46,9 @@ class Orbitals { // ==> Class Methods <== /// @return The alpha orbital coefficient matrix - const std::shared_ptr Ca() const; + std::shared_ptr Ca(); /// @return The beta orbital coefficient matrix - const std::shared_ptr Cb() const; + std::shared_ptr Cb(); void set(const std::shared_ptr& Ca, const std::shared_ptr& Cb); @@ -71,7 +71,7 @@ class Orbitals { std::shared_ptr Cb_; }; -std::unique_ptr -make_orbitals_from_psi(const std::shared_ptr& wfn, bool restricted); +std::shared_ptr make_orbitals_from_psi(const std::shared_ptr& wfn, + bool restricted); } // namespace forte diff --git a/forte/integrals/integrals_psi4_interface.cc b/forte/integrals/integrals_psi4_interface.cc index ca46e8ba7..999ff6086 100644 --- a/forte/integrals/integrals_psi4_interface.cc +++ b/forte/integrals/integrals_psi4_interface.cc @@ -117,8 +117,8 @@ void Psi4Integrals::transform_one_electron_integrals() { auto Ha = wfn_->H()->clone(); auto Hb = wfn_->H()->clone(); - Ha->transform(*orbitals_->Ca()); - Hb->transform(*orbitals_->Cb()); + Ha->transform(orbitals_->Ca()); + Hb->transform(orbitals_->Cb()); OneBody_symm_ = Ha; @@ -279,9 +279,9 @@ void Psi4Integrals::update_orbitals(std::shared_ptr orbitals, bool re_ orbitals_->Ca()->print(); orbitals_->Cb()->print(); auto overlap = - psi::linalg::triplet(*orbitals->Ca(), *S_, *orbitals->Cb(), true, false, false); - overlap.set_name("Overlap "); - overlap.print(); + psi::linalg::triplet(orbitals->Ca(), S_, orbitals->Cb(), true, false, false); + overlap->set_name("Overlap "); + overlap->print(); auto msg = "Psi4Integrals::update_orbitals was passed two different sets of orbitals" "\n but the integral object assumes restricted orbitals"; throw std::runtime_error(msg); diff --git a/forte/modules/mcscf.py b/forte/modules/mcscf.py index 2802722f9..dd3b7a353 100644 --- a/forte/modules/mcscf.py +++ b/forte/modules/mcscf.py @@ -12,7 +12,6 @@ class MCSCF(Module): - """ A module to perform MCSCF calculations. """ @@ -39,6 +38,7 @@ def _run(self, data: ForteData) -> ForteData: data.scf_info, data.options, data.mo_space_info, + data.orbitals, data.ints, ) energy = mcscf.compute_energy() diff --git a/forte/modules/objects_factory_psi4.py b/forte/modules/objects_factory_psi4.py index 04cae1101..f7b3a05e8 100644 --- a/forte/modules/objects_factory_psi4.py +++ b/forte/modules/objects_factory_psi4.py @@ -15,6 +15,7 @@ make_mo_space_info_from_map, make_state_weights_map, SCFInfo, + make_orbitals_from_psi, make_ints_from_psi4, ) @@ -209,7 +210,10 @@ def prepare_forte_objects_from_psi4_wfn(options, wfn, mo_space_info): # Build a map from Forte StateInfo to the weights state_weights_map = make_state_weights_map(options, mo_space_info) - return (state_weights_map, mo_space_info, scf_info) + # Build the Orbitals object + orbitals = make_orbitals_from_psi(wfn, True) # True for restricted orbitals + + return (state_weights_map, mo_space_info, scf_info, orbitals) def prepare_forte_objects(data, name, **kwargs): @@ -226,13 +230,14 @@ 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) forte_objects = prepare_forte_objects_from_psi4_wfn(options, ref_wfn, mo_space_info) - state_weights_map, mo_space_info, scf_info = forte_objects + state_weights_map, mo_space_info, scf_info, orbitals = forte_objects fcidump = None data.mo_space_info = mo_space_info data.scf_info = scf_info data.state_weights_map = state_weights_map data.psi_wfn = ref_wfn + data.orbitals = orbitals return data, fcidump @@ -276,6 +281,6 @@ def _run(self, data: ForteData = None) -> ForteData: else: psi4.core.print_out("\n Forte will use psi4 integrals") # Make an integral object from the psi4 wavefunction object - data.ints = make_ints_from_psi4(data.psi_wfn, data.options, data.mo_space_info) + data.ints = make_ints_from_psi4(data.psi_wfn, data.options, data.mo_space_info, data.orbitals) return data diff --git a/forte/post_process/spin_corr.h b/forte/post_process/spin_corr.h index 45f3ed726..2c67635fc 100644 --- a/forte/post_process/spin_corr.h +++ b/forte/post_process/spin_corr.h @@ -75,6 +75,7 @@ class SpinCorr { void perform_spin_analysis(std::shared_ptr rdms, std::shared_ptr options, std::shared_ptr mo_space_info, + std::shared_ptr orbitals, std::shared_ptr as_ints); } // namespace forte From 555fa1f1839175019b2296cb278a2d3368e5f4aa Mon Sep 17 00:00:00 2001 From: fevangelista Date: Sat, 18 May 2024 15:03:39 -0400 Subject: [PATCH 3/4] Only two test cases remaining --- forte/integrals/custom_integrals.cc | 4 - forte/integrals/integrals_psi4_interface.cc | 3 + forte/mcscf/mcscf_orb_grad.cc | 4 +- forte/modules/objects_factory_fcidump.py | 9 ++ forte/modules/orbital_transform.py | 5 +- forte/pymodule.py | 4 +- tests/methods/external_solver-1/as_ints.json | 8 +- tests/methods/external_solver-2/as_ints.json | 8 +- .../methods/external_solver-2/dsrg_ints.json | 90 +++++++++---------- 9 files changed, 72 insertions(+), 63 deletions(-) diff --git a/forte/integrals/custom_integrals.cc b/forte/integrals/custom_integrals.cc index 7aef478c5..f80aa6673 100644 --- a/forte/integrals/custom_integrals.cc +++ b/forte/integrals/custom_integrals.cc @@ -70,10 +70,6 @@ CustomIntegrals::CustomIntegrals(std::shared_ptr options, } void CustomIntegrals::initialize() { - // Ca_ = std::make_shared(nmopi_, nmopi_); - // Cb_ = std::make_shared(nmopi_, nmopi_); - // Ca_->identity(); - // Cb_->identity(); nsopi_ = nmopi_; nso_ = nmo_; diff --git a/forte/integrals/integrals_psi4_interface.cc b/forte/integrals/integrals_psi4_interface.cc index 999ff6086..a0cb3e4fc 100644 --- a/forte/integrals/integrals_psi4_interface.cc +++ b/forte/integrals/integrals_psi4_interface.cc @@ -113,6 +113,9 @@ void Psi4Integrals::setup_psi4_ints() { } void Psi4Integrals::transform_one_electron_integrals() { + std::cout << " From transform_one_electron_integrals skip_build_: " + << (skip_build_ ? "True" : "False") << std::endl; + // Grab the one-electron integrals from psi4's wave function object auto Ha = wfn_->H()->clone(); auto Hb = wfn_->H()->clone(); diff --git a/forte/mcscf/mcscf_orb_grad.cc b/forte/mcscf/mcscf_orb_grad.cc index c0a7ed98e..6e98eeb3a 100644 --- a/forte/mcscf/mcscf_orb_grad.cc +++ b/forte/mcscf/mcscf_orb_grad.cc @@ -758,7 +758,7 @@ bool MCSCF_ORB_GRAD::update_orbitals(std::shared_ptr x) { orbitals_->set(C_, C_); // update integrals object. Except for custom integrals, we update the integrals // object with the new orbitals but do not update the integrals for efficiency - ints_->update_orbitals(orbitals_, ints_->integral_type() != IntegralType::Custom); + ints_->update_orbitals(orbitals_, ints_->integral_type() == IntegralType::Custom); // if (ints_->integral_type() != Custom) { // ints_->update_orbitals(orbitals_); // } else { @@ -1096,7 +1096,7 @@ void MCSCF_ORB_GRAD::canonicalize_final(const std::shared_ptr& U) { orbitals_->set(C_, C_); // update integrals object. Except for custom integrals, we update the integrals // object with the new orbitals but do not update the integrals for efficiency - ints_->update_orbitals(orbitals_, ints_->integral_type() != IntegralType::Custom); + ints_->update_orbitals(orbitals_, ints_->integral_type() == IntegralType::Custom); build_mo_integrals(); } } // namespace forte diff --git a/forte/modules/objects_factory_fcidump.py b/forte/modules/objects_factory_fcidump.py index 386b85fa5..0d5cfc863 100644 --- a/forte/modules/objects_factory_fcidump.py +++ b/forte/modules/objects_factory_fcidump.py @@ -29,6 +29,7 @@ def _make_ints_from_fcidump(fcidump, data: ForteData): ints = forte.make_custom_ints( data.options, data.mo_space_info, + data.orbitals, fcidump["enuc"], fcidump["hcore"].flatten(), fcidump["hcore"].flatten(), @@ -142,6 +143,14 @@ def _prepare_forte_objects_from_fcidump(data, filename: str = None): data.scf_info = forte.SCFInfo(nmopi, doccpi, soccpi, 0.0, epsilon_a, epsilon_b) + # Create the Orbitals object + # create two sets of orbitals, one for alpha and one for beta of dimnsion + Ca = psi4.core.Matrix("Ca", nmopi, nmopi) + Cb = psi4.core.Matrix("Cb", nmopi, nmopi) + Ca.identity() + Cb.identity() + data.orbitals = forte.Orbitals(Ca, Cb) + state_info = _make_state_info_from_fcidump(fcidump, options) data.state_weights_map = {state_info: [1.0]} data.psi_wfn = None diff --git a/forte/modules/orbital_transform.py b/forte/modules/orbital_transform.py index 0fd341ebe..7b8b0c523 100644 --- a/forte/modules/orbital_transform.py +++ b/forte/modules/orbital_transform.py @@ -6,7 +6,6 @@ class OrbitalTransformation(Module): - """ A module to transform the orbitals and the integrals to a new basis. """ @@ -23,7 +22,9 @@ def __init__(self, orb_type: str, transform_ints: bool = True): self.transform_ints = transform_ints def _run(self, data: ForteData) -> ForteData: - orb_t = make_orbital_transformation(self.orb_type, data.scf_info, data.options, data.ints, data.mo_space_info) + orb_t = make_orbital_transformation( + self.orb_type, data.scf_info, data.options, data.ints, data.mo_space_info, data.orbitals + ) orb_t.compute_transformation() Ua = orb_t.get_Ua() Ub = orb_t.get_Ub() diff --git a/forte/pymodule.py b/forte/pymodule.py index 25b485e84..18e71cc97 100644 --- a/forte/pymodule.py +++ b/forte/pymodule.py @@ -101,7 +101,7 @@ def forte_driver(data: ForteData): if options.get_bool("SPIN_ANALYSIS"): data = ActiveSpaceRDMs(max_rdm_level=2, rdms_type=forte.RDMsType.spin_dependent).run(data) - forte.perform_spin_analysis(data.rdms, options, mo_space_info, as_ints) + forte.perform_spin_analysis(data.rdms, options, mo_space_info, data.orbitals, as_ints) # solver for dynamical correlation from DSRG correlation_solver_type = options.get_str("CORRELATION_SOLVER") @@ -275,7 +275,7 @@ def gradient_forte(name, **kwargs): # Make an integral object time_pre_ints = time.time() - data.ints = forte.make_ints_from_psi4(data.psi_wfn, data.options, data.mo_space_info) + data.ints = forte.make_ints_from_psi4(data.psi_wfn, data.options, data.mo_space_info, data.orbitals) start = time.time() diff --git a/tests/methods/external_solver-1/as_ints.json b/tests/methods/external_solver-1/as_ints.json index 701da0b5c..2b331bdce 100644 --- a/tests/methods/external_solver-1/as_ints.json +++ b/tests/methods/external_solver-1/as_ints.json @@ -16,7 +16,7 @@ [ 0, 0, - -1.2796255434474166 + -1.2780734354814196 ], [ 0, @@ -31,12 +31,12 @@ [ 2, 2, - -0.4816690241888207 + -0.5069808360680924 ], [ 1, 1, - -1.2796255434474166 + -1.2780734354814196 ], [ 1, @@ -51,7 +51,7 @@ [ 3, 3, - -0.4816690241888207 + -0.5069808360680924 ] ], "description": "one-electron integrals as a list of tuples (i,j,)" diff --git a/tests/methods/external_solver-2/as_ints.json b/tests/methods/external_solver-2/as_ints.json index 701da0b5c..2b331bdce 100644 --- a/tests/methods/external_solver-2/as_ints.json +++ b/tests/methods/external_solver-2/as_ints.json @@ -16,7 +16,7 @@ [ 0, 0, - -1.2796255434474166 + -1.2780734354814196 ], [ 0, @@ -31,12 +31,12 @@ [ 2, 2, - -0.4816690241888207 + -0.5069808360680924 ], [ 1, 1, - -1.2796255434474166 + -1.2780734354814196 ], [ 1, @@ -51,7 +51,7 @@ [ 3, 3, - -0.4816690241888207 + -0.5069808360680924 ] ], "description": "one-electron integrals as a list of tuples (i,j,)" diff --git a/tests/methods/external_solver-2/dsrg_ints.json b/tests/methods/external_solver-2/dsrg_ints.json index 693fe0f61..e090d32d3 100644 --- a/tests/methods/external_solver-2/dsrg_ints.json +++ b/tests/methods/external_solver-2/dsrg_ints.json @@ -16,7 +16,7 @@ [ 0, 0, - -1.2941570603785824 + -1.294683588798959 ], [ 0, @@ -31,12 +31,12 @@ [ 2, 2, - -0.44130879303420895 + -0.5076397260343521 ], [ 1, 1, - -1.2941570603785824 + -1.294683588798959 ], [ 1, @@ -51,13 +51,13 @@ [ 3, 3, - -0.44130879303420883 + -0.5076397260343521 ] ], "description": "one-electron integrals as a list of tuples (i,j,)" }, "scalar_energy": { - "data": 0.7519948111691335, + "data": 0.7543669510517765, "description": "scalar energy (sum of nuclear repulsion, frozen core, and scalar contributions" }, "spin": { @@ -96,28 +96,28 @@ 1, 0, 1, - 0.6922779923932212 + 0.6834545927868173 ], [ 0, 1, 1, 0, - -0.6922779923932212 + -0.6834545927868173 ], [ 1, 0, 0, 1, - -0.6922779923932212 + -0.6834545927868173 ], [ 1, 0, 1, 0, - 0.6922779923932212 + 0.6834545927868173 ], [ 1, @@ -222,28 +222,28 @@ 1, 2, 3, - 0.14262265005657665 + 0.05511921767176892 ], [ 0, 1, 3, 2, - -0.14262265005657665 + -0.05511921767176892 ], [ 1, 0, 2, 3, - -0.14262265005657665 + -0.05511921767176892 ], [ 1, 0, 3, 2, - 0.14262265005657665 + 0.05511921767176892 ], [ 1, @@ -299,84 +299,84 @@ 2, 0, 2, - 0.44142483066923077 + 0.32687770302232827 ], [ 0, 3, 0, 3, - 0.5897358195853001 + 0.3760159534345901 ], [ 0, 3, 3, 0, - -0.5897358195853001 + -0.3760159534345901 ], [ 3, 0, 0, 3, - -0.5897358195853001 + -0.3760159534345901 ], [ 3, 0, 3, 0, - 0.5897358195853001 + 0.3760159534345901 ], [ 1, 3, 1, 3, - 0.44142483066923077 + 0.32687770302232827 ], [ 0, 2, 2, 0, - -0.44142483066923077 + -0.32687770302232827 ], [ 0, 3, 2, 1, - 0.1483109889160694 + 0.049138250412261844 ], [ 0, 3, 1, 2, - -0.1483109889160694 + -0.049138250412261844 ], [ 3, 0, 2, 1, - -0.1483109889160694 + -0.049138250412261844 ], [ 3, 0, 1, 2, - 0.1483109889160694 + 0.049138250412261844 ], [ 1, 3, 3, 1, - -0.44142483066923077 + -0.32687770302232827 ], [ 0, @@ -467,84 +467,84 @@ 0, 0, 2, - -0.44142483066923077 + -0.32687770302232827 ], [ 2, 1, 0, 3, - 0.1483109889160694 + 0.049138250412261844 ], [ 2, 1, 3, 0, - -0.1483109889160694 + -0.049138250412261844 ], [ 1, 2, 0, 3, - -0.1483109889160694 + -0.049138250412261844 ], [ 1, 2, 3, 0, - 0.1483109889160694 + 0.049138250412261844 ], [ 3, 1, 1, 3, - -0.44142483066923077 + -0.32687770302232827 ], [ 2, 0, 2, 0, - 0.44142483066923077 + 0.32687770302232827 ], [ 2, 1, 2, 1, - 0.5897358195853002 + 0.3760159534345901 ], [ 2, 1, 1, 2, - -0.5897358195853002 + -0.3760159534345901 ], [ 1, 2, 2, 1, - -0.5897358195853002 + -0.3760159534345901 ], [ 1, 2, 1, 2, - 0.5897358195853002 + 0.3760159534345901 ], [ 3, 1, 3, 1, - 0.44142483066923077 + 0.32687770302232827 ], [ 2, @@ -600,28 +600,28 @@ 3, 0, 1, - 0.14262265005657665 + 0.05511921767176892 ], [ 2, 3, 1, 0, - -0.14262265005657665 + -0.05511921767176892 ], [ 3, 2, 0, 1, - -0.14262265005657665 + -0.05511921767176892 ], [ 3, 2, 1, 0, - 0.14262265005657665 + 0.05511921767176892 ], [ 3, @@ -726,28 +726,28 @@ 3, 2, 3, - 0.651413811309362 + 0.32364110338005236 ], [ 2, 3, 3, 2, - -0.651413811309362 + -0.32364110338005236 ], [ 3, 2, 2, 3, - -0.651413811309362 + -0.32364110338005236 ], [ 3, 2, 3, 2, - 0.651413811309362 + 0.32364110338005236 ], [ 3, From 8b0f3dc6ddae0cbd4b64b649fdc6057ba45a5622 Mon Sep 17 00:00:00 2001 From: fevangelista Date: Sun, 19 May 2024 15:27:25 -0400 Subject: [PATCH 4/4] Adds skeleteon JK class --- forte/CMakeLists.txt | 3 + forte/api/forte_python_module.cc | 2 + forte/integrals/customjk_builder.cc | 40 ++++++++++ forte/integrals/integrals.cc | 10 +-- forte/integrals/integrals.h | 35 +++++++-- forte/integrals/jk_builder.cc | 41 ++++++++++ forte/integrals/jk_builder.h | 90 ++++++++++++++++++++++ forte/integrals/psi4jk_builder.cc | 42 ++++++++++ forte/mcscf/mcscf_2step.cc | 9 ++- forte/mcscf/mcscf_2step.h | 7 +- forte/modules/mcscf.py | 8 ++ forte/proc/external_active_space_solver.py | 2 + 12 files changed, 270 insertions(+), 19 deletions(-) create mode 100644 forte/integrals/customjk_builder.cc create mode 100644 forte/integrals/jk_builder.cc create mode 100644 forte/integrals/jk_builder.h create mode 100644 forte/integrals/psi4jk_builder.cc diff --git a/forte/CMakeLists.txt b/forte/CMakeLists.txt index a06e288b1..0f9953264 100644 --- a/forte/CMakeLists.txt +++ b/forte/CMakeLists.txt @@ -169,6 +169,9 @@ integrals/make_integrals.cc integrals/one_body_integrals.cc integrals/parallel_ccvv_algorithms.cc integrals/paralleldfmo.cc +integrals/jk_builder.cc +integrals/customjk_builder.cc +integrals/psi4jk_builder.cc mrdsrg-helper/dsrg_mem.cc mrdsrg-helper/dsrg_source.cc mrdsrg-helper/dsrg_time.cc diff --git a/forte/api/forte_python_module.cc b/forte/api/forte_python_module.cc index f405ad05a..a21199a73 100644 --- a/forte/api/forte_python_module.cc +++ b/forte/api/forte_python_module.cc @@ -47,6 +47,7 @@ #include "base_classes/orbitals.h" #include "integrals/make_integrals.h" +#include "integrals/jk_builder.h" #include "orbital-helpers/aosubspace.h" #include "orbital-helpers/mp2_nos.h" @@ -190,6 +191,7 @@ PYBIND11_MODULE(_forte, m) { m.def("make_fragment_projector", &make_fragment_projector, "Make a fragment(embedding) projector"); m.def("make_embedding", &make_embedding, "Apply fragment projector to embed"); + m.def("make_psi4jk", &make_psi4jk, "Make a Forte JK object based on psi4"); m.def("make_custom_ints", &make_custom_forte_integrals, "Make a custom Forte integral object from arrays"); m.def("make_ints_from_psi4", &make_forte_integrals_from_psi4, "ref_wfn"_a, "options"_a, diff --git a/forte/integrals/customjk_builder.cc b/forte/integrals/customjk_builder.cc new file mode 100644 index 000000000..d7d4b9542 --- /dev/null +++ b/forte/integrals/customjk_builder.cc @@ -0,0 +1,40 @@ +/* + * @BEGIN LICENSE + * + * Forte: an open-source plugin to Psi4 (https://github.com/psi4/psi4) + * that implements a variety of quantum chemistry methods for strongly + * correlated electrons. + * + * Copyright (c) 2012-2024 by its authors (see COPYING, COPYING.LESSER, AUTHORS). + * + * The copyrights for code used from other parties are included in + * the corresponding files. + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see http://www.gnu.org/licenses/. + * + * @END LICENSE + */ + +#include "jk_builder.h" + +using namespace psi; +using namespace ambit; + +namespace forte { + +CustomJKBuilder::CustomJKBuilder() {} +CustomJKBuilder::~CustomJKBuilder() {} +void CustomJKBuilder::buildJK() {} + +} // namespace forte diff --git a/forte/integrals/integrals.cc b/forte/integrals/integrals.cc index ffb2c1e05..e509bed98 100644 --- a/forte/integrals/integrals.cc +++ b/forte/integrals/integrals.cc @@ -72,25 +72,23 @@ ForteIntegrals::ForteIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, std::shared_ptr mo_space_info, std::shared_ptr orbitals, IntegralType integral_type, - IntegralSpinRestriction restricted, DFType df_type) + IntegralSpinRestriction restricted) : options_(options), mo_space_info_(mo_space_info), orbitals_(orbitals), wfn_(ref_wfn), - integral_type_(integral_type), spin_restriction_(restricted), df_type_(df_type), - frozen_core_energy_(0.0), scalar_energy_(0.0) { + integral_type_(integral_type), spin_restriction_(restricted) { common_initialize(); } ForteIntegrals::ForteIntegrals(std::shared_ptr options, std::shared_ptr mo_space_info, std::shared_ptr orbitals, IntegralType integral_type, - IntegralSpinRestriction restricted, DFType df_type) + IntegralSpinRestriction restricted) : options_(options), mo_space_info_(mo_space_info), orbitals_(orbitals), - integral_type_(integral_type), spin_restriction_(restricted), df_type_(df_type) { + integral_type_(integral_type), spin_restriction_(restricted) { common_initialize(); } void ForteIntegrals::common_initialize() { read_information(); - if (not skip_build_) { allocate(); } diff --git a/forte/integrals/integrals.h b/forte/integrals/integrals.h index 99cf34ab3..b99b0cde8 100644 --- a/forte/integrals/integrals.h +++ b/forte/integrals/integrals.h @@ -124,8 +124,7 @@ class ForteIntegrals { ForteIntegrals(std::shared_ptr options, std::shared_ptr ref_wfn, std::shared_ptr mo_space_info, std::shared_ptr orbitals, - IntegralType integral_type, IntegralSpinRestriction restricted, - DFType df_type = DFType::JK); + IntegralType integral_type, IntegralSpinRestriction restricted); /** * @brief Class constructor @@ -135,8 +134,7 @@ class ForteIntegrals { */ ForteIntegrals(std::shared_ptr options, std::shared_ptr mo_space_info, std::shared_ptr orbitals, - IntegralType integral_type, IntegralSpinRestriction restricted, - DFType df_type = DFType::JK); + IntegralType integral_type, IntegralSpinRestriction restricted); /// Virtual destructor to enable deletion of a Derived* through a Base* virtual ~ForteIntegrals() = default; @@ -436,9 +434,6 @@ class ForteIntegrals { /// Are we doing a spin-restricted computation? IntegralSpinRestriction spin_restriction_; - /// The type of density fitting basis - DFType df_type_; - // AO overlap matrix from psi std::shared_ptr S_; @@ -653,4 +648,30 @@ class Psi4Integrals : public ForteIntegrals { double schwarz_cutoff_; }; +/// ForteJK class +class ForteJK { + public: + ForteJK(std::shared_ptr options, std::shared_ptr wfn, + std::shared_ptr mo_space_info, std::shared_ptr orbitals); + + /// Initialize JK object + void initialize(); + + /// Finalize JK object + void finalize(); + + /// Compute JK + void compute_JK(); + + /// Return the JK object + std::shared_ptr jk(); + + private: + std::shared_ptr JK_; + std::shared_ptr options_; + std::shared_ptr wfn_; + std::shared_ptr mo_space_info_; + std::shared_ptr orbitals_; +}; + } // namespace forte diff --git a/forte/integrals/jk_builder.cc b/forte/integrals/jk_builder.cc new file mode 100644 index 000000000..5f79736c4 --- /dev/null +++ b/forte/integrals/jk_builder.cc @@ -0,0 +1,41 @@ +/* + * @BEGIN LICENSE + * + * Forte: an open-source plugin to Psi4 (https://github.com/psi4/psi4) + * that implements a variety of quantum chemistry methods for strongly + * correlated electrons. + * + * Copyright (c) 2012-2024 by its authors (see COPYING, COPYING.LESSER, AUTHORS). + * + * The copyrights for code used from other parties are included in + * the corresponding files. + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see http://www.gnu.org/licenses/. + * + * @END LICENSE + */ + +#include "jk_builder.h" + +using namespace psi; +using namespace ambit; + +namespace forte { + +std::shared_ptr make_psi4jk(std::shared_ptr options, + std::shared_ptr wfn) { + return std::make_shared(options, wfn); +} + +} // namespace forte diff --git a/forte/integrals/jk_builder.h b/forte/integrals/jk_builder.h new file mode 100644 index 000000000..d3fc5814d --- /dev/null +++ b/forte/integrals/jk_builder.h @@ -0,0 +1,90 @@ +/* + * @BEGIN LICENSE + * + * Forte: an open-source plugin to Psi4 (https://github.com/psi4/psi4) + * that implements a variety of quantum chemistry methods for strongly + * correlated electrons. + * + * Copyright (c) 2012-2024 by its authors (see COPYING, COPYING.LESSER, + * AUTHORS). + * + * The copyrights for code used from other parties are included in + * the corresponding files. + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see http://www.gnu.org/licenses/. + * + * @END LICENSE + */ + +#pragma once + +#include + +#include "psi4/libfock/jk.h" +#include "psi4/libmints/dimension.h" +#include "ambit/blocked_tensor.h" + +class Tensor; + +namespace psi { +class Options; +class Matrix; +class Vector; +class Wavefunction; +class Dimension; +class BasisSet; +} // namespace psi + +namespace forte { + +class ForteOptions; +class MOSpaceInfo; +class Orbitals; + +class JKBuilder { + public: + virtual ~JKBuilder() {} + + // Pure virtual function to build the J and K matrices + virtual void buildJK() = 0; + + // Additional common utility functions can be added here + protected: + // Protected constructor to prevent instantiation of base class + JKBuilder() {} +}; + +class Psi4JKBuilder : public JKBuilder { + public: + Psi4JKBuilder(std::shared_ptr options, std::shared_ptr wfn); + ~Psi4JKBuilder(); + void buildJK() override; + + private: + std::shared_ptr options_; + std::shared_ptr wfn_; + std::shared_ptr JK_; +}; + +class CustomJKBuilder : public JKBuilder { + public: + CustomJKBuilder(); + ~CustomJKBuilder(); + void buildJK() override; +}; + +std::shared_ptr make_psi4jk(std::shared_ptr options, + std::shared_ptr wfn); + +} // namespace forte diff --git a/forte/integrals/psi4jk_builder.cc b/forte/integrals/psi4jk_builder.cc new file mode 100644 index 000000000..30aafeef9 --- /dev/null +++ b/forte/integrals/psi4jk_builder.cc @@ -0,0 +1,42 @@ +/* + * @BEGIN LICENSE + * + * Forte: an open-source plugin to Psi4 (https://github.com/psi4/psi4) + * that implements a variety of quantum chemistry methods for strongly + * correlated electrons. + * + * Copyright (c) 2012-2024 by its authors (see COPYING, COPYING.LESSER, AUTHORS). + * + * The copyrights for code used from other parties are included in + * the corresponding files. + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see http://www.gnu.org/licenses/. + * + * @END LICENSE + */ + +#include "jk_builder.h" + +using namespace psi; +using namespace ambit; + +namespace forte { + +Psi4JKBuilder::Psi4JKBuilder(std::shared_ptr options, + std::shared_ptr wfn) {} + +Psi4JKBuilder::~Psi4JKBuilder() {} +void Psi4JKBuilder::buildJK() {} + +} // namespace forte diff --git a/forte/mcscf/mcscf_2step.cc b/forte/mcscf/mcscf_2step.cc index ee06a610b..7cd22932b 100644 --- a/forte/mcscf/mcscf_2step.cc +++ b/forte/mcscf/mcscf_2step.cc @@ -59,9 +59,10 @@ MCSCF_2STEP::MCSCF_2STEP(std::shared_ptr as_solver, std::shared_ptr mo_space_info, std::shared_ptr orbitals, std::shared_ptr scf_info, - std::shared_ptr ints) + std::shared_ptr ints, std::shared_ptr jk) : as_solver_(as_solver), state_weights_map_(state_weights_map), options_(options), - mo_space_info_(mo_space_info), orbitals_(orbitals), scf_info_(scf_info), ints_(ints) { + mo_space_info_(mo_space_info), orbitals_(orbitals), scf_info_(scf_info), ints_(ints), + jk_(jk) { startup(); } @@ -596,9 +597,9 @@ make_mcscf_two_step(std::shared_ptr as_solver, const std::map>& state_weight_map, std::shared_ptr scf_info, std::shared_ptr options, std::shared_ptr mo_space_info, std::shared_ptr orbitals, - std::shared_ptr ints) { + std::shared_ptr ints, std::shared_ptr jk) { return std::make_unique(as_solver, state_weight_map, options, mo_space_info, - orbitals, scf_info, ints); + orbitals, scf_info, ints, jk); } } // namespace forte diff --git a/forte/mcscf/mcscf_2step.h b/forte/mcscf/mcscf_2step.h index 936985fa7..3afceec82 100644 --- a/forte/mcscf/mcscf_2step.h +++ b/forte/mcscf/mcscf_2step.h @@ -58,7 +58,7 @@ class MCSCF_2STEP { const std::map>& state_weights_map, std::shared_ptr options, std::shared_ptr mo_space_info, std::shared_ptr orbitals, std::shared_ptr scf_info, - std::shared_ptr ints); + std::shared_ptr ints, std::shared_ptr jk); /// Compute the MCSCF energy double compute_energy(); @@ -85,6 +85,9 @@ class MCSCF_2STEP { /// The Forte integral std::shared_ptr ints_; + /// The Forte JK object + std::shared_ptr jk_; + /// Common setup for the class void startup(); @@ -182,6 +185,6 @@ make_mcscf_two_step(std::shared_ptr as_solver, const std::map>& state_weight_map, std::shared_ptr ref_wfn, std::shared_ptr options, std::shared_ptr mo_space_info, std::shared_ptr orbitals, - std::shared_ptr ints); + std::shared_ptr ints, std::shared_ptr jk); } // namespace forte diff --git a/forte/modules/mcscf.py b/forte/modules/mcscf.py index dd3b7a353..e55b7e86f 100644 --- a/forte/modules/mcscf.py +++ b/forte/modules/mcscf.py @@ -8,6 +8,7 @@ make_mcscf_two_step, MOSpaceInfo, make_mo_space_info_from_map, + make_psi4jk, ) @@ -32,6 +33,12 @@ def _run(self, data: ForteData) -> ForteData: data.active_space_solver = make_active_space_solver( self.solver_type, state_map, data.scf_info, data.mo_space_info, data.options ) + + if data.psi_wfn: + data.jk = make_psi4jk(data.options, data.psi_wfn) + else: + raise ValueError("Only psi4 JK is supported for now.") + mcscf = make_mcscf_two_step( data.active_space_solver, data.state_weights_map, @@ -40,6 +47,7 @@ def _run(self, data: ForteData) -> ForteData: data.mo_space_info, data.orbitals, data.ints, + data.jk, ) energy = mcscf.compute_energy() data.results.add("energy", energy, "MCSCF energy", "hartree") diff --git a/forte/proc/external_active_space_solver.py b/forte/proc/external_active_space_solver.py index 230a45c0d..21d990683 100644 --- a/forte/proc/external_active_space_solver.py +++ b/forte/proc/external_active_space_solver.py @@ -39,6 +39,8 @@ def read_wavefunction(data): C_mat = psi4.core.Matrix.from_array([np.asarray(C_list)]) data.psi_wfn.Ca().copy(C_mat) data.psi_wfn.Cb().copy(C_mat) + # set the orbitals in the data object + data.orbitals = forte.Orbitals(C_mat, C_mat) def write_external_active_space_file(as_ints, state_map, mo_space_info, json_file="forte_ints.json"):