Skip to content

Commit

Permalink
ref_relax
Browse files Browse the repository at this point in the history
  • Loading branch information
@shuhangli98
shuhangli98 committed Nov 12, 2024
1 parent 65e2b98 commit 11eb4dd
Show file tree
Hide file tree
Showing 4 changed files with 176 additions and 83 deletions.
4 changes: 4 additions & 0 deletions forte/api/orbital_api.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,8 +78,12 @@ void export_SemiCanonical(py::module& m) {
"Semicanonicalize the orbitals and transform the integrals and reference")
.def("Ua", &SemiCanonical::Ua, "Return the alpha rotation matrix")
.def("Ub", &SemiCanonical::Ub, "Return the alpha rotation matrix")
.def("Ua_c", &SemiCanonical::Ua_c, "Return the alpha rotation matrix in the core space")
.def("Ub_c", &SemiCanonical::Ub_c, "Return the beta rotation matrix in the core space")
.def("Ua_t", &SemiCanonical::Ua_t, "Return the alpha rotation matrix in the active space")
.def("Ub_t", &SemiCanonical::Ub_t, "Return the beta rotation matrix in the active space")
.def("Ua_v", &SemiCanonical::Ua_v, "Return the alpha rotation matrix in the virtual space")
.def("Ub_v", &SemiCanonical::Ub_v, "Return the beta rotation matrix in the virtual space")
.def("fix_orbital_success", &SemiCanonical::fix_orbital_success,
"Return if the orbital ordering and phases are fixed successfully");
}
Expand Down
74 changes: 74 additions & 0 deletions forte/orbital-helpers/semi_canonicalize.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,16 +74,26 @@ void SemiCanonical::read_options(const std::shared_ptr<ForteOptions>& options) {
void SemiCanonical::startup() {
nirrep_ = mo_space_info_->nirrep();
nmopi_ = mo_space_info_->dimension("ALL");
ncore_ = mo_space_info_->size("CORE");
nact_ = mo_space_info_->size("ACTIVE");
nvirt_ = mo_space_info_->size("VIRTUAL");

// Prepare orbital rotation matrix, which transforms all MOs
Ua_ = std::make_shared<psi::Matrix>("Ua", nmopi_, nmopi_);
Ub_ = std::make_shared<psi::Matrix>("Ub", nmopi_, nmopi_);

// Prepare orbital rotation matrix, which transforms only core MOs
Ua_c_ = ambit::Tensor::build(ambit::CoreTensor, "Ua", {ncore_, ncore_});
Ub_c_ = ambit::Tensor::build(ambit::CoreTensor, "Ub", {ncore_, ncore_});

// Prepare orbital rotation matrix, which transforms only active MOs
Ua_t_ = ambit::Tensor::build(ambit::CoreTensor, "Ua", {nact_, nact_});
Ub_t_ = ambit::Tensor::build(ambit::CoreTensor, "Ub", {nact_, nact_});

// Prepare orbital rotation matrix, which transforms only virtual MOs
Ua_v_ = ambit::Tensor::build(ambit::CoreTensor, "Ua", {nvirt_, nvirt_});
Ub_v_ = ambit::Tensor::build(ambit::CoreTensor, "Ub", {nvirt_, nvirt_});

// Initialize U to identity
set_U_to_identity();

Expand Down Expand Up @@ -120,11 +130,23 @@ void SemiCanonical::set_U_to_identity() {
Ua_->identity();
Ub_->identity();

Ua_c_.iterate(
[&](const std::vector<size_t>& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; });

Ub_c_.iterate(
[&](const std::vector<size_t>& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; });

Ua_t_.iterate(
[&](const std::vector<size_t>& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; });

Ub_t_.iterate(
[&](const std::vector<size_t>& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; });

Ua_v_.iterate(
[&](const std::vector<size_t>& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; });

Ub_v_.iterate(
[&](const std::vector<size_t>& i, double& value) { value = (i[0] == i[1]) ? 1.0 : 0.0; });
}

void SemiCanonical::semicanonicalize(std::shared_ptr<RDMs> rdms, const bool& build_fock,
Expand Down Expand Up @@ -274,12 +296,64 @@ void SemiCanonical::build_transformation_matrices(const bool& semi) {
// keep phase and order unchanged
fix_orbital_success_ = ints_->fix_orbital_phases(Ua_, true);

// fill in Ua_c_
fill_Ucore(Ua_, Ua_c_);
// fill in Ua_t_
fill_Uactv(Ua_, Ua_t_);
// fill in Ua_v_
fill_Uvirt(Ua_, Ua_v_);

// pass to Ub
Ub_->copy(Ua_);
Ub_c_.copy(Ua_c_);
Ub_t_.copy(Ua_t_);
Ub_v_.copy(Ua_v_);
}

void SemiCanonical::fill_Ucore(const std::shared_ptr<psi::Matrix>& U, ambit::Tensor& Ut) {
auto core_names = mo_space_info_->composite_space_names()["CORE"];
auto& Ut_data = Ut.data();

for (const std::string& name : core_names) {
auto size = mo_space_info_->size(name);
if (size == 0)
continue;

auto pos = mo_space_info_->pos_in_space(name, "CORE");
auto relative_mos = mo_space_info_->relative_mo(name);
for (size_t p = 0; p < size; ++p) {
const auto& [hp, np] = relative_mos[p];
for (size_t q = 0; q < size; ++q) {
const auto& [hq, nq] = relative_mos[q];
if (hp != hq)
continue;
Ut_data[pos[p] * ncore_ + pos[q]] = U->get(hp, np, nq);
}
}
}
}

void SemiCanonical::fill_Uvirt(const std::shared_ptr<psi::Matrix>& U, ambit::Tensor& Ut) {
auto virtual_names = mo_space_info_->composite_space_names()["VIRTUAL"];
auto& Ut_data = Ut.data();

for (const std::string& name : virtual_names) {
auto size = mo_space_info_->size(name);
if (size == 0)
continue;

auto pos = mo_space_info_->pos_in_space(name, "VIRTUAL");
auto relative_mos = mo_space_info_->relative_mo(name);
for (size_t p = 0; p < size; ++p) {
const auto& [hp, np] = relative_mos[p];
for (size_t q = 0; q < size; ++q) {
const auto& [hq, nq] = relative_mos[q];
if (hp != hq)
continue;
Ut_data[pos[p] * nvirt_ + pos[q]] = U->get(hp, np, nq);
}
}
}
}

void SemiCanonical::fill_Uactv(const std::shared_ptr<psi::Matrix>& U, ambit::Tensor& Ut) {
Expand Down
49 changes: 41 additions & 8 deletions forte/orbital-helpers/semi_canonicalize.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,16 +42,16 @@ namespace forte {
class ForteOptions;

/// @brief The SemiCanonical class
/// This class computes semi-canonical orbitals from the 1RDM and optionally transforms the integrals and RDMs
/// Semi-canonical orbitals are obtained by diagonalizing the Fock matrix in each orbital space separately
/// The class can also produce natural orbitals. These differ by the semi-canonical orbital only in the active space
/// where they are defined to be eigenvectors of the 1RDM
///
/// The final orbitals are ordered by increasing energy within each irrep and space. Natural orbitals are ordered
/// by decreasing occupation number
/// This class computes semi-canonical orbitals from the 1RDM and optionally transforms the
/// integrals and RDMs Semi-canonical orbitals are obtained by diagonalizing the Fock matrix in each
/// orbital space separately The class can also produce natural orbitals. These differ by the
/// semi-canonical orbital only in the active space where they are defined to be eigenvectors of the
/// 1RDM
///
/// The final orbitals are ordered by increasing energy within each irrep and space. Natural
/// orbitals are ordered by decreasing occupation number
class SemiCanonical {
public:

/// @brief SemiCanonical Constructor
/// @param mo_space_info The MOSpaceInfo object
/// @param ints The ForteIntegrals object
Expand All @@ -78,12 +78,22 @@ class SemiCanonical {
/// @return the beta rotation matrix
std::shared_ptr<psi::Matrix> Ub() { return Ub_; }

/// @return the alpha rotation matrix in the core space
ambit::Tensor Ua_c() const { return Ua_c_.clone(); }
/// @return the beta rotation matrix in the core space
ambit::Tensor Ub_c() const { return Ub_c_.clone(); }

/// @return the alpha rotation matrix in the active space
ambit::Tensor Ua_t() const { return Ua_t_.clone(); }

/// @return the beta rotation matrix in the active space
ambit::Tensor Ub_t() const { return Ub_t_.clone(); }

/// @return the alpha rotation matrix in the virtual space
ambit::Tensor Ua_v() const { return Ua_v_.clone(); }
/// @return the beta rotation matrix in the virtual space
ambit::Tensor Ub_v() const { return Ub_v_.clone(); }

/// @return if the orbital ordering and phases are fixed successfully
bool fix_orbital_success() const { return fix_orbital_success_; }

Expand Down Expand Up @@ -117,21 +127,38 @@ class SemiCanonical {
/// Offset of GAS orbitals within ACTIVE
std::map<std::string, psi::Dimension> actv_offsets_;

/// Number of core MOs
size_t ncore_;

/// Number of active MOs
size_t nact_;

/// Number of virtual MOs
size_t nvirt_;

/// Number of irreps
size_t nirrep_;

/// Unitary matrix for alpha orbital rotation
std::shared_ptr<psi::Matrix> Ua_;
/// Unitary matrix for beta orbital rotation
std::shared_ptr<psi::Matrix> Ub_;

/// Unitary matrix for alpha orbital rotation in the core space
ambit::Tensor Ua_c_;
/// Unitary matrix for beta orbital rotation in the core space
ambit::Tensor Ub_c_;

/// Unitary matrix for alpha orbital rotation in the active space
ambit::Tensor Ua_t_;
/// Unitary matrix for beta orbital rotation in the active space
ambit::Tensor Ub_t_;

/// Unitary matrix for alpha orbital rotation in the virtual space
ambit::Tensor Ua_v_;
/// Unitary matrix for beta orbital rotation in the virtual space
ambit::Tensor Ub_v_;

/// Set Ua_, Ub_, Ua_t_, and Ub_t_ to identity
void set_U_to_identity();

Expand All @@ -153,9 +180,15 @@ class SemiCanonical {
/// Builds unitary matrices used to diagonalize diagonal blocks of Fock
void build_transformation_matrices(const bool& semi);

/// Fill ambit::Tensor Ua_c_ (Ub_c_) using std::shared_ptr<psi::Matrix> Ua_ (Ub_)
void fill_Ucore(const std::shared_ptr<psi::Matrix>& U, ambit::Tensor& Ut);

/// Fill ambit::Tensor Ua_t_ (Ub_t_) using std::shared_ptr<psi::Matrix> Ua_ (Ub_)
void fill_Uactv(const std::shared_ptr<psi::Matrix>& U, ambit::Tensor& Ut);

/// Fill ambit::Tensor Ua_v_ (Ub_v_) using std::shared_ptr<psi::Matrix> Ua_ (Ub_)
void fill_Uvirt(const std::shared_ptr<psi::Matrix>& U, ambit::Tensor& Ut);

/// Successfully fix the orbital ordering and phases
bool fix_orbital_success_;
};
Expand Down
132 changes: 57 additions & 75 deletions forte/proc/dsrg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -246,41 +246,11 @@ def compute_energy(self):
if not self.Heff_implemented:
self.relax_maxiter = 0

# if self.options.get_bool('FULL_HBAR'):
# self.relax_maxiter = 0

if self.options.get_bool("FULL_HBAR") and self.relax_maxiter == 0:
psi4.core.print_out("\n =>** Saving Full Hbar **<=\n")
psi4.core.print_out("\n =>** Saving Full Hbar (unrelax) **<=\n")
Heff = self.dsrg_solver.compute_Heff_full()
Heff_dict = forte.Heff_dict(Heff)
np.savez("save_Hbar", **Heff_dict)
psi4.core.print_out("\n =>** Getting gamma1 **<=\n")
gamma1 = self.dsrg_solver.get_gamma1()
psi4.core.print_out("\n =>** Getting eta1 **<=\n")
eta1 = self.dsrg_solver.get_eta1()
psi4.core.print_out("\n =>** Getting lambda2 **<=\n")
lambda2 = self.dsrg_solver.get_lambda2()
psi4.core.print_out("\n =>** Getting lambda3 **<=\n")
lambda3 = self.dsrg_solver.get_lambda3()
gamma1_dict = forte.blocktensor_to_dict(gamma1)
eta_1_dict = forte.blocktensor_to_dict(eta1)
lambda2_dict = forte.blocktensor_to_dict(lambda2)
if self.solver_type in ["MRDSRG_SO", "MRDSRG-SO"]:
lambda3_dict = forte.blocktensor_to_dict(lambda3)
else:
lambda3_dict = forte.L3_dict(lambda3)

np.savez("save_gamma1", **gamma1_dict)
np.savez("save_eta1", **eta_1_dict)
np.savez("save_lambda2", **lambda2_dict)
np.savez("save_lambda3", **lambda3_dict)
del gamma1, eta1, lambda2, lambda3, gamma1_dict, eta_1_dict, lambda2_dict, lambda3_dict
if self.options.get_str("FOURPDC") != "ZERO":
psi4.core.print_out("\n =>** Getting lambda4 **<=\n")
lambda4 = self.dsrg_solver.get_lambda4()
lambda4_dict = forte.L4_dict(lambda4)
np.savez("save_lambda4", **lambda4_dict)
del lambda4, lambda4_dict

if self.solver_type not in ["MRDSRG_SO", "MRDSRG-SO"]:
psi4.core.print_out("\n =>** Getting dipole integral **<=\n")
Expand All @@ -305,6 +275,27 @@ def compute_energy(self):
# so that the CI vectors are comparable before and after DSRG dressing.
# However, the ForteIntegrals object and the dipole integrals always refer to the current semi-canonical basis.
# so to compute the dipole moment correctly, we need to make the RDMs and orbital basis consistent

if self.options.get_bool("FULL_HBAR") and n == self.relax_maxiter - 1:
psi4.core.print_out("\n =>** Saving Full Hbar (relax) **<=\n")
Heff = self.dsrg_solver.compute_Heff_full()
Heff_dict = forte.Heff_dict(Heff)
np.savez("save_Hbar", **Heff_dict)

if self.solver_type not in ["MRDSRG_SO", "MRDSRG-SO"]:
psi4.core.print_out("\n =>** Getting dipole integral **<=\n")
Mbar0 = self.dsrg_solver.compute_Mbar0_full()
print(Mbar0)
np.save("Mbar0", Mbar0)
Mbar1 = self.dsrg_solver.compute_Mbar1_full()
Mbar2 = self.dsrg_solver.compute_Mbar2_full()

for i in range(3):
np.savez(f"Mbar1_{i}", **forte.blocktensor_to_dict(Mbar1[i]))
np.savez(f"Mbar2_{i}", **forte.blocktensor_to_dict(Mbar2[i]))

del Mbar0, Mbar1, Mbar2

ints_dressed = self.dsrg_solver.compute_Heff_actv()
if self.fno_pt2_Heff_shift is not None:
ints_dressed.add(self.fno_pt2_Heff_shift, 1.0)
Expand Down Expand Up @@ -434,52 +425,43 @@ def compute_energy(self):
psi4.core.print_out(f"\n\n DSRG-MRPT2 FNO energy correction: {self.fno_pt2_energy_shift:20.15f}")
psi4.core.print_out(f"\n DSRG-MRPT2 FNO corrected energy: {e_dsrg:20.15f}")

if self.options.get_bool("FULL_HBAR"):
self.rdms.rotate(self.Ua, self.Ub)
psi4.core.print_out("\n =>** Saving Full Hbar **<=\n")
Heff = self.dsrg_solver.compute_Heff_full()
Heff_dict = forte.Heff_dict(Heff)
np.savez("save_Hbar", **Heff_dict)
psi4.core.print_out("\n =>** Getting gamma1 **<=\n")
gamma1 = self.dsrg_solver.get_gamma1()
psi4.core.print_out("\n =>** Getting eta1 **<=\n")
eta1 = self.dsrg_solver.get_eta1()
psi4.core.print_out("\n =>** Getting lambda2 **<=\n")
lambda2 = self.dsrg_solver.get_lambda2()
psi4.core.print_out("\n =>** Getting lambda3 **<=\n")
lambda3 = self.dsrg_solver.get_lambda3()
gamma1_dict = forte.blocktensor_to_dict(gamma1)
eta_1_dict = forte.blocktensor_to_dict(eta1)
lambda2_dict = forte.blocktensor_to_dict(lambda2)
if self.solver_type in ["MRDSRG_SO", "MRDSRG-SO"]:
lambda3_dict = forte.blocktensor_to_dict(lambda3)
else:
lambda3_dict = forte.L3_dict(lambda3)

np.savez("save_gamma1", **gamma1_dict)
np.savez("save_eta1", **eta_1_dict)
np.savez("save_lambda2", **lambda2_dict)
np.savez("save_lambda3", **lambda3_dict)
del gamma1, eta1, lambda2, lambda3, gamma1_dict, eta_1_dict, lambda2_dict, lambda3_dict
if self.options.get_str("FOURPDC") != "ZERO":
psi4.core.print_out("\n =>** Getting lambda4 **<=\n")
lambda4 = self.dsrg_solver.get_lambda4()
lambda4_dict = forte.L4_dict(lambda4)
np.savez("save_lambda4", **lambda4_dict)
del lambda4, lambda4_dict
if self.solver_type not in ["MRDSRG_SO", "MRDSRG-SO"]:
psi4.core.print_out("\n =>** Getting dipole integral **<=\n")
Mbar0 = self.dsrg_solver.compute_Mbar0_full()
print(Mbar0)
np.save("Mbar0", Mbar0)
Mbar1 = self.dsrg_solver.compute_Mbar1_full()
Mbar2 = self.dsrg_solver.compute_Mbar2_full()
if self.options.get_bool("FULL_HBAR"):
if self.relax_maxiter != 0:
self.rdms = self.active_space_solver.compute_average_rdms(
self.state_weights_map, self.max_rdm_level, self.rdm_type
)
self.rdms.rotate(self.Ua, self.Ub) # To previous semi-canonical basis

for i in range(3):
np.savez(f"Mbar1_{i}", **forte.blocktensor_to_dict(Mbar1[i]))
np.savez(f"Mbar2_{i}", **forte.blocktensor_to_dict(Mbar2[i]))
self.make_dsrg_solver()

del Mbar0, Mbar1, Mbar2
psi4.core.print_out("\n =>** Getting gamma1 **<=\n")
gamma1 = self.dsrg_solver.get_gamma1()
psi4.core.print_out("\n =>** Getting eta1 **<=\n")
eta1 = self.dsrg_solver.get_eta1()
psi4.core.print_out("\n =>** Getting lambda2 **<=\n")
lambda2 = self.dsrg_solver.get_lambda2()
psi4.core.print_out("\n =>** Getting lambda3 **<=\n")
lambda3 = self.dsrg_solver.get_lambda3()

gamma1_dict = forte.blocktensor_to_dict(gamma1)
eta1_dict = forte.blocktensor_to_dict(eta1)
lambda2_dict = forte.blocktensor_to_dict(lambda2)
if self.solver_type in ["MRDSRG_SO", "MRDSRG-SO"]:
lambda3_dict = forte.blocktensor_to_dict(lambda3)
else:
lambda3_dict = forte.L3_dict(lambda3)

np.savez("save_gamma1", **gamma1_dict)
np.savez("save_eta1", **eta1_dict)
np.savez("save_lambda2", **lambda2_dict)
np.savez("save_lambda3", **lambda3_dict)
del gamma1, eta1, lambda2, lambda3, gamma1_dict, eta1_dict, lambda2_dict, lambda3_dict
if self.options.get_str("FOURPDC") != "ZERO":
psi4.core.print_out("\n =>** Getting lambda4 **<=\n")
lambda4 = self.dsrg_solver.get_lambda4()
lambda4_dict = forte.L4_dict(lambda4)
np.savez("save_lambda4", **lambda4_dict)
del lambda4, lambda4_dict

self.dsrg_cleanup()

Expand Down

0 comments on commit 11eb4dd

Please sign in to comment.