Skip to content

Commit

Permalink
Degno full for SI
Browse files Browse the repository at this point in the history
  • Loading branch information
@brianz98
brianz98 committed Nov 15, 2024
1 parent 426cc7b commit bb89c72
Show file tree
Hide file tree
Showing 17 changed files with 278 additions and 295 deletions.
1 change: 0 additions & 1 deletion forte/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -217,7 +217,6 @@ mrdsrg-spin-adapted/sadsrg.cc
mrdsrg-spin-adapted/sadsrg_amps_analysis.cc
mrdsrg-spin-adapted/sadsrg_block_labels.cc
mrdsrg-spin-adapted/sadsrg_comm.cc
mrdsrg-spin-adapted/eom_dsrg.cc
mrdsrg-spin-integrated/dsrg_mrpt2.cc
mrdsrg-spin-integrated/dsrg_mrpt2_grad/dsrg_mrpt2_deriv_multipliers.cc
mrdsrg-spin-integrated/dsrg_mrpt2_grad/dsrg_mrpt2_deriv_setb.cc
Expand Down
135 changes: 57 additions & 78 deletions forte/api/forte_python_module.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,71 +181,6 @@ PYBIND11_MODULE(_forte, m) {
},
"Return the cumulants of the RDMs in a spinorbital basis. Spinorbitals follow the ordering "
"abab...");
m.def(
"Heff_dict",
[](std::vector<ambit::BlockedTensor> Heff) {
py::dict pyints;
for (const auto& int_ : Heff) {
auto labels = int_.block_labels();
for (const auto& label : labels) {
pyints[py::str(label)] = ambit_to_np(int_.block(label));
}
}
return pyints;
},
"Return the full Heff in a dictionary");
// m.def(
// "Mbar_dict",
// [](std::vector<ambit::BlockedTensor> Mbar) {
// std::vector<py::dict> pymbar;
// py::dict pyints;
// for (size_t i = 0; i < Mbar.size(); i++) {
// pyints = {};
// for (const auto& int_ : Mbar[i]) {
// auto labels = int_.block_labels();
// for (const auto& label : labels) {
// pyints[py::str(label)] = ambit_to_np(int_.block(label));
// }
// }
// pymbar.push_back(pyints);
// }
// return pymbar;
// },
// "Return the full Meff in a list of dictionary");
m.def(
"blocktensor_to_dict",
[](ambit::BlockedTensor rdm) {
py::dict pyrdm;
auto labels = rdm.block_labels();
for (const auto& label : labels) {
pyrdm[py::str(label)] = ambit_to_np(rdm.block(label));
}
return pyrdm;
},
"Return the BlockTensor in a dictionary");
m.def(
"L3_dict",
[](std::vector<ambit::Tensor> rdm) {
py::dict pyrdm;
pyrdm[py::str("aaaaaa")] = ambit_to_np(rdm[0]);
pyrdm[py::str("aaAaaA")] = ambit_to_np(rdm[1]);
pyrdm[py::str("aAAaAA")] = ambit_to_np(rdm[2]);
pyrdm[py::str("AAAAAA")] = ambit_to_np(rdm[3]);
return pyrdm;
},
"Return the L3 in a dictionary");
m.def(
"L4_dict", // L4aaaa_, L4aaab_, L4aabb_, L4abbb_, L4bbbb_
[](std::vector<ambit::Tensor> rdm) {
py::dict pyrdm;
pyrdm[py::str("aaaaaaaa")] = ambit_to_np(rdm[0]);
pyrdm[py::str("aaaAaaaA")] = ambit_to_np(rdm[1]);
pyrdm[py::str("aaAAaaAA")] = ambit_to_np(rdm[2]);
pyrdm[py::str("aAAAaAAA")] = ambit_to_np(rdm[3]);
pyrdm[py::str("AAAAAAAA")] = ambit_to_np(rdm[4]);
return pyrdm;
},
"Return the L3 in a dictionary");
m.def("get_gas_occupation", &get_gas_occupation);
m.def("get_ci_occupation_patterns", &get_ci_occupation_patterns);

Expand Down Expand Up @@ -314,19 +249,57 @@ PYBIND11_MODULE(_forte, m) {
.def("compute_gradient", &MASTER_DSRG::compute_gradient, "Compute the DSRG gradient")
.def("compute_Heff_actv", &MASTER_DSRG::compute_Heff_actv,
"Return the DSRG dressed ActiveSpaceIntegrals")
.def("compute_Heff_full", &MASTER_DSRG::compute_Heff_full,
"Return full transformed Hamiltonian")
.def("compute_Heff_full", [](MASTER_DSRG& self) {
const auto Heff = self.compute_Heff_full();
return py::make_tuple(blockedtensor_to_np(Heff.at(0)), blockedtensor_to_np(Heff.at(1)));
})
.def("compute_Heff_full_degno", [](MASTER_DSRG& self) {
const auto Heff = self.compute_Heff_full_degno();
return py::make_tuple(blockedtensor_to_np(Heff.at(0)), blockedtensor_to_np(Heff.at(1)));
})
.def("compute_Mbar0_full", &MASTER_DSRG::compute_Mbar0_full,
"Return full transformed zero-body dipole integrals")
.def("compute_Mbar1_full", &MASTER_DSRG::compute_Mbar1_full,
"Return full transformed one-body dipole integrals")
.def("compute_Mbar2_full", &MASTER_DSRG::compute_Mbar2_full,
"Return full transformed two-body dipole integrals")
.def("get_gamma1", &MASTER_DSRG::get_gamma1, "Return the gamma1 tensor")
.def("get_eta1", &MASTER_DSRG::get_eta1, "Return the eta1 tensor")
.def("get_lambda2", &MASTER_DSRG::get_lambda2, "Return the lambda2 tensor")
.def("get_lambda3", &MASTER_DSRG::get_lambda3, "Return the lambda3 tensor")
.def("get_lambda4", &MASTER_DSRG::get_lambda4, "Return the lambda4 tensor")
.def("compute_Mbar1_full", [](MASTER_DSRG& self) {
const auto Mbar1 = self.compute_Mbar1_full();
return py::make_tuple(blockedtensor_to_np(Mbar1.at(0)), blockedtensor_to_np(Mbar1.at(1)),
blockedtensor_to_np(Mbar1.at(2)));
})
.def("compute_Mbar2_full", [](MASTER_DSRG& self) {
const auto Mbar2 = self.compute_Mbar2_full();
return py::make_tuple(blockedtensor_to_np(Mbar2.at(0)), blockedtensor_to_np(Mbar2.at(1)),
blockedtensor_to_np(Mbar2.at(2)));
})
.def("get_gamma1", [](MASTER_DSRG& self) {
const auto gamma1 = self.get_gamma1();
return blockedtensor_to_np(gamma1);
})
.def("get_eta1", [](MASTER_DSRG& self) {
const auto eta1 = self.get_eta1();
return blockedtensor_to_np(eta1);
})
.def("get_lambda2", [](MASTER_DSRG& self) {
const auto lambda2 = self.get_lambda2();
return blockedtensor_to_np(lambda2);
})
.def("get_lambda3", [](MASTER_DSRG& self) {
const auto lambda3 = self.get_lambda3();
py::dict pyrdm;
pyrdm[py::str("aaaaaa")] = ambit_to_np(lambda3.at(0));
pyrdm[py::str("aaAaaA")] = ambit_to_np(lambda3.at(1));
pyrdm[py::str("aAAaAA")] = ambit_to_np(lambda3.at(2));
pyrdm[py::str("AAAAAA")] = ambit_to_np(lambda3.at(3));
return pyrdm;
})
.def("get_lambda4", [](MASTER_DSRG& self) {
const auto lambda4 = self.get_lambda4();
py::dict pyrdm;
pyrdm[py::str("aaaaaaaa")] = ambit_to_np(lambda4.at(0));
pyrdm[py::str("aaaAaaaA")] = ambit_to_np(lambda4.at(1));
pyrdm[py::str("aaAAaaAA")] = ambit_to_np(lambda4.at(2));
pyrdm[py::str("aAAAaAAA")] = ambit_to_np(lambda4.at(3));
pyrdm[py::str("AAAAAAAA")] = ambit_to_np(lambda4.at(4));
return pyrdm;
})
.def("deGNO_DMbar_actv", &MASTER_DSRG::deGNO_DMbar_actv,
"Return the DSRG dressed dipole integrals")
.def("nuclear_dipole", &MASTER_DSRG::nuclear_dipole,
Expand All @@ -352,8 +325,14 @@ PYBIND11_MODULE(_forte, m) {
.def("compute_energy", &SADSRG::compute_energy, "Compute the DSRG energy")
.def("compute_Heff_actv", &SADSRG::compute_Heff_actv,
"Return the DSRG dressed ActiveSpaceIntegrals")
.def("compute_Heff_full", &SADSRG::compute_Heff_full,
"Return full transformed Hamiltonian")
.def("compute_Heff_full", [](SADSRG& self) {
const auto Heff = self.compute_Heff_full();
return py::make_tuple(blockedtensor_to_np(Heff.at(0)), blockedtensor_to_np(Heff.at(1)));
})
.def("compute_Heff_full_degno", [](SADSRG& self) {
const auto Heff = self.compute_Heff_full_degno();
return py::make_tuple(blockedtensor_to_np(Heff.at(0)), blockedtensor_to_np(Heff.at(1)));
})
.def("compute_mp_eff_actv", &SADSRG::compute_mp_eff_actv,
"Return the DSRG dressed ActiveMultipoleIntegrals")
.def("set_Uactv", &SADSRG::set_Uactv, "Ua"_a,
Expand Down
9 changes: 9 additions & 0 deletions forte/helpers/helpers.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,15 @@ py::array_t<double> vector_to_np(const std::vector<double>& v, const std::vector
return py::array_t<double>(dims, &(v.data()[0]));
}

py::dict blockedtensor_to_np(ambit::BlockedTensor bt) {
py::dict pybt;
auto labels = bt.block_labels();
for (const auto& label : labels) {
pybt[py::str(label)] = ambit_to_np(bt.block(label));
}
return pybt;
}

std::shared_ptr<psi::Matrix> tensor_to_matrix(ambit::Tensor t) {
size_t size1 = t.dim(0);
size_t size2 = t.dim(1);
Expand Down
8 changes: 8 additions & 0 deletions forte/helpers/helpers.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,14 @@ enum class Spin3 { aaa, aab, abb, bbb };
*/
py::array_t<double> ambit_to_np(ambit::Tensor t);

/**
* @brief Convert an ambit BlockedTensor to a dictionary of numpy ndarray's.
* The returned tensor is stored according to the C storage convention.
* @param t The input BlockedTensor
* @return A dictionary of numpy arrays
*/
py::dict blockedtensor_to_np(ambit::BlockedTensor bt);

/**
* @brief Convert a std::vector<double> to a numpy ndarray.
* The returned tensor is stored according to the C storage convention.
Expand Down
111 changes: 16 additions & 95 deletions forte/mrdsrg-so/mrdsrg_so.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,13 @@ void MRDSRG_SO::startup() {

source_ = foptions_->get_str("SOURCE");

dsrg_trans_type_ = foptions_->get_str("DSRG_TRANS_TYPE");
if (dsrg_trans_type_ == "CC" && foptions_->get_str("CORR_LEVEL") == "QDSRG2"){
outfile->Printf("\n Warning: DSRG_TRANS_TYPE option CC is not supported with CORR_LEVEL QDSRG2.");
outfile->Printf("\n Changed DSRG_TRANS_TYPE option to UNITARY");
dsrg_trans_type_ = "UNITARY";
}

ntamp_ = foptions_->get_int("NTAMP");
intruder_tamp_ = foptions_->get_double("INTRUDER_TAMP");

Expand Down Expand Up @@ -645,99 +652,11 @@ double MRDSRG_SO::compute_energy() {
outfile->Printf("\n\n\n MR-DSRG(2) correlation energy = %25.15f", Etotal - Eref);
outfile->Printf("\n * MR-DSRG(2) total energy = %25.15f\n", Etotal);

// compute_eom();

psi::Process::environment.globals["CURRENT ENERGY"] = Etotal;

return Etotal;
}

// double MRDSRG_SO::compute_eom() {
// // IP singles
// size_t ncore = nc_ / 2;
// size_t nactv = na_ / 2;
// size_t nmo = nso_ / 2;
// size_t nhole = ncore + nactv;

// auto EOM_Hbar = BTF_->build(tensor_type_, "EOM-Hbar", {"hh"});
// std::shared_ptr<psi::Matrix> EOM_Hbar_mat_ =
// std::make_shared<psi::Matrix>("EOM-Hbar-Matrx", nh_, nh_);

// EOM_Hbar["nm"] = -Hbar1["nm"];
// EOM_Hbar["mv"] = -Hbar1["mu"] * Gamma1["uv"];
// EOM_Hbar["mv"] += 0.5 * Hbar2["mwux"] * Lambda2["uxwv"];
// EOM_Hbar["vm"] = EOM_Hbar["mv"];
// EOM_Hbar["wx"] = -Hbar1["vu"] * Gamma1["ux"] * Gamma1["wv"];
// EOM_Hbar["wx"] += Hbar1["vu"] * Lambda2["uwvx"];
// EOM_Hbar["wx"] += 0.5 * Hbar2["yzuv"] * Gamma1["wy"] * Lambda2["uvzx"];
// EOM_Hbar["wx"] -= 0.5 * Hbar2["yzuv"] * Gamma1["vx"] * Lambda2["uwyz"];
// EOM_Hbar["wx"] += 0.25 * Hbar2["yzuv"] * Lambda3["uvwyzx"];

// EOM_Hbar.iterate(
// [&](const std::vector<size_t>& i, const std::vector<SpinType>&, double& value) {
// if (i[0] < nmo && i[1] < nmo) {
// EOM_Hbar_mat_->set(i[0], i[1], value);
// }
// if (i[0] >= nmo && i[1] >= nmo) {
// EOM_Hbar_mat_->set(i[0] - nmo + nhole, i[1] - nmo + nhole, value);
// }
// });

// EOM_Hbar_mat_->print();

// /// Overlap matrix
// std::shared_ptr<psi::Matrix> S = std::make_shared<psi::Matrix>("EOM-S-sub", nh_, nh_);
// std::shared_ptr<psi::Matrix> Sevec = std::make_shared<psi::Matrix>("S-evec", nh_, nh_);
// std::shared_ptr<psi::Vector> Seval = std::make_shared<psi::Vector>("S-eval", nh_);
// // std::shared_ptr<psi::Matrix> Eevec = std::make_shared<psi::Matrix>("E-evec", nh_, nh_);
// // std::shared_ptr<psi::Vector> Eeval = std::make_shared<psi::Vector>("E-eval", nh_);
// S->identity();
// (rdms_->g1a()).citerate([&](const std::vector<size_t>& i, const double& value) {
// S->set(i[0] + ncore, i[1] + ncore, value);
// });
// (rdms_->g1b()).citerate([&](const std::vector<size_t>& i, const double& value) {
// std::cout << value << std::endl;
// S->set(i[0] + 2 * ncore + nactv, i[1] + 2 * ncore + nactv, value);
// });
// S->print();
// S->diagonalize(Sevec, Seval);
// // auto Sevec = S->partial_cholesky_factorize(1e-10);
// // auto Seval = std::make_shared<psi::Vector>("S-eval", Sevec->rowdim());

// std::vector<size_t> s_list;
// for (size_t i = 0; i < Sevec->coldim(); i++) {
// double value = abs(Seval->get(i));
// if (value > 1e-8) {
// s_list.push_back(i);
// }
// }

// size_t num_s_list = s_list.size();

// std::shared_ptr<psi::Matrix> S_new = std::make_shared<psi::Matrix>("S-new", nh_, num_s_list);
// for (size_t i = 0; i < num_s_list; i++) {
// psi::SharedVector temp = Sevec->get_column(0, s_list[i]);
// temp->scale(1.0 / Seval->get(s_list[i]));
// S_new->set_column(0, i, temp);
// }
// std::shared_ptr<psi::Matrix> Eevec = std::make_shared<psi::Matrix>("E-evec", nh_,
// num_s_list); std::shared_ptr<psi::Vector> Eeval = std::make_shared<psi::Vector>("E-eval",
// num_s_list);

// EOM_Hbar_mat_ = psi::linalg::triplet(S_new, EOM_Hbar_mat_, S_new, true, false, false);

// EOM_Hbar_mat_->diagonalize(Eevec, Eeval);

// outfile->Printf("\n "
// "----------------------------------------------------------"
// "----------------------------------------");
// outfile->Printf("\n\n\n EOM-DSRG ");

// Eeval->print();

// return Eeval->get(0);
// }

void MRDSRG_SO::compute_hbar() {

// outfile->Printf("\n\n Computing the similarity-transformed
Expand Down Expand Up @@ -799,13 +718,15 @@ void MRDSRG_SO::compute_hbar() {
// outfile->Printf("\n |H2| = %20.12f", C2.norm(1));
// outfile->Printf("\n --------------------------------");

// [H, A] = [H, T] + [H, T]^dagger
C0 *= 2.0;
O1["pq"] = C1["pq"];
C1["pq"] += O1["qp"];
O2["pqrs"] = C2["pqrs"];
C2["pqrs"] += O2["rspq"];

if (dsrg_trans_type_ == "UNITARY"){
// [H, A] = [H, T] + [H, T]^dagger
C0 *= 2.0;
O1["pq"] = C1["pq"];
C1["pq"] += O1["qp"];
O2["pqrs"] = C2["pqrs"];
C2["pqrs"] += O2["rspq"];
}

// Hbar += C
Hbar0 += C0;
Hbar1["pq"] += C1["pq"];
Expand Down
3 changes: 3 additions & 0 deletions forte/mrdsrg-so/mrdsrg_so.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,6 +120,9 @@ class MRDSRG_SO : public DynamicCorrelationSolver {
/// Order of the Taylor expansion of f(z) = (1-exp(-z^2))/z
int taylor_order_;

/// DSRG transformation type
std::string dsrg_trans_type_;

std::shared_ptr<BlockedTensorFactory> BTF_;
TensorType tensor_type_;

Expand Down
Loading

0 comments on commit bb89c72

Please sign in to comment.