Skip to content

Commit

Permalink
refactor
Browse files Browse the repository at this point in the history
  • Loading branch information
@ajaypanyala
ajaypanyala committed Apr 14, 2024
1 parent 6e7c151 commit e7faa91
Show file tree
Hide file tree
Showing 4 changed files with 11 additions and 130 deletions.
52 changes: 0 additions & 52 deletions exachem/scf/scf_hartree_fock.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@ void SCFHartreeFock::initialize(ExecutionContext& exc, ChemEnv& chem_env) { scf_

void SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) {
SCFCompute scf_compute;
ISCFQed scf_qed;
ISCFIter scf_iter;
ISCFGuess scf_guess;
ISCFRestart scf_restart;
Expand Down Expand Up @@ -301,27 +300,6 @@ void SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) {
/*** =========================== ***/
scf_compute.compute_hamiltonian<TensorType>(ec, scf_vars, chem_env.atoms, chem_env.shells,
ttensors, etensors);
if(chem_env.sys_data.is_qed) {
ttensors.QED_Dx = {tAO, tAO};
ttensors.QED_Dy = {tAO, tAO};
ttensors.QED_Dz = {tAO, tAO};
ttensors.QED_Qxx = {tAO, tAO};
ttensors.QED_Qxy = {tAO, tAO};
ttensors.QED_Qxz = {tAO, tAO};
ttensors.QED_Qyy = {tAO, tAO};
ttensors.QED_Qyz = {tAO, tAO};
ttensors.QED_Qzz = {tAO, tAO};
ttensors.QED_1body = {tAO, tAO};
ttensors.QED_2body = {tAO, tAO};
Tensor<TensorType>::allocate(&ec, ttensors.QED_Dx, ttensors.QED_Dy, ttensors.QED_Dz,
ttensors.QED_Qxx, ttensors.QED_Qxy, ttensors.QED_Qxz,
ttensors.QED_Qyy, ttensors.QED_Qyz, ttensors.QED_Qzz,
ttensors.QED_1body, ttensors.QED_2body);

scf_qed.compute_qed_emult_ints<TensorType>(ec, chem_env, scf_vars, ttensors);
if(chem_env.sys_data.do_qed)
scf_qed.compute_QED_1body<TensorType>(ec, chem_env, scf_vars, ttensors);
}

if(chem_env.sys_data.has_ecp) {
Tensor<TensorType> ECP{tAO, tAO};
Expand Down Expand Up @@ -659,12 +637,6 @@ void SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) {
shell2bf, SchwarzK, max_nprim4, ttensors, etensors,
is_3c_init, do_density_fitting, xHF);

// Add QED contribution

if(chem_env.sys_data.do_qed) {
scf_qed.compute_QED_2body<TensorType>(ec, chem_env, scf_vars, ttensors);
}

TensorType gauxc_exc = 0.;
#if defined(USE_GAUXC)
if(chem_env.sys_data.is_ks) {
Expand Down Expand Up @@ -738,11 +710,6 @@ void SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) {
shell2bf, SchwarzK, max_nprim4, ttensors, etensors,
is_3c_init, do_density_fitting, xHF);

// Add QED contribution
if(chem_env.sys_data.do_qed) {
scf_qed.compute_QED_2body<TensorType>(ec, chem_env, scf_vars, ttensors);
}

std::tie(ehf, rmsd) = scf_iter.scf_iter_body<TensorType>(ec, chem_env, scalapack_info, iter,
scf_vars, ttensors, etensors
#if defined(USE_GAUXC)
Expand Down Expand Up @@ -845,13 +812,6 @@ void SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) {
scf_iter.compute_2bf<TensorType>(ec, chem_env, scalapack_info, scf_vars, do_schwarz_screen,
shell2bf, SchwarzK, max_nprim4, ttensors, etensors,
is_3c_init, do_density_fitting, xHF_adjust);

// Add QED contribution;
// CHECK

if(chem_env.sys_data.do_qed) {
scf_qed.compute_QED_2body<TensorType>(ec, chem_env, scf_vars, ttensors);
}
}
else if(scf_vars.lshift > 0) {
// Remove level shift from Fock matrix
Expand Down Expand Up @@ -930,18 +890,6 @@ void SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) {
Tensor<TensorType>::deallocate(ttensors.VXC_alpha);
if(chem_env.sys_data.is_unrestricted) Tensor<TensorType>::deallocate(ttensors.VXC_beta);
}
if(chem_env.sys_data.is_qed) {
scf_output.rw_mat_disk<TensorType>(ttensors.QED_Dx, files_prefix + ".QED_Dx",
chem_env.ioptions.scf_options.debug);
scf_output.rw_mat_disk<TensorType>(ttensors.QED_Dy, files_prefix + ".QED_Dy",
chem_env.ioptions.scf_options.debug);
scf_output.rw_mat_disk<TensorType>(ttensors.QED_Dz, files_prefix + ".QED_Dz",
chem_env.ioptions.scf_options.debug);
Tensor<TensorType>::deallocate(ttensors.QED_Dx, ttensors.QED_Dy, ttensors.QED_Dz,
ttensors.QED_Qxx, ttensors.QED_Qxy, ttensors.QED_Qxz,
ttensors.QED_Qyy, ttensors.QED_Qyz, ttensors.QED_Qzz,
ttensors.QED_1body, ttensors.QED_2body);
}
#if SCF_THROTTLE_RESOURCES
ec.flush_and_sync();
#endif
Expand Down
1 change: 0 additions & 1 deletion exachem/scf/scf_hartree_fock.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@
#include "scf_compute.hpp"
#include "scf_iter.hpp"
#include "scf_outputs.hpp"
#include "scf_qed.hpp"
#include "scf_restart.hpp"
#include "scf_taskmap.hpp"
#include <variant>
Expand Down
68 changes: 7 additions & 61 deletions exachem/scf/scf_outputs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,71 +23,21 @@ void SCFIO::print_energies(ExecutionContext& ec, const ChemEnv& chem_env, TAMMTe
const bool is_uhf = sys_data.is_unrestricted;
const bool is_rhf = sys_data.is_restricted;
const bool is_ks = sys_data.is_ks;
const bool is_qed = sys_data.is_qed;
const bool do_qed = sys_data.do_qed;

double nelectrons = 0.0;
double kinetic_1e = 0.0;
double NE_1e = 0.0;
double energy_1e = 0.0;
double energy_2e = 0.0;
double energy_qed = 0.0;
double energy_qed_et = 0.0;
double nelectrons = 0.0;
double kinetic_1e = 0.0;
double NE_1e = 0.0;
double energy_1e = 0.0;
double energy_2e = 0.0;

if(is_rhf) {
nelectrons = tt_trace(ec, ttensors.D_last_alpha, ttensors.S1);
kinetic_1e = tt_trace(ec, ttensors.D_last_alpha, ttensors.T1);
NE_1e = tt_trace(ec, ttensors.D_last_alpha, ttensors.V1);
energy_1e = tt_trace(ec, ttensors.D_last_alpha, ttensors.H1);
energy_2e = 0.5 * tt_trace(ec, ttensors.D_last_alpha, ttensors.F_alpha_tmp);

if(is_ks) {
if((!is_qed) || (is_qed && do_qed)) { energy_2e += scf_vars.exc; }
}

if(is_qed) {
if(do_qed) {
energy_qed = tt_trace(ec, ttensors.D_last_alpha, ttensors.QED_1body);
energy_qed += 0.5 * tt_trace(ec, ttensors.D_last_alpha, ttensors.QED_2body);
}
else { energy_qed = scf_vars.eqed; }

Tensor<double> X_a;

#if defined(USE_SCALAPACK)
X_a = {scf_vars.tAO, scf_vars.tAO_ortho};
Tensor<double>::allocate(&ec, X_a);
if(scalapack_info.comm != MPI_COMM_NULL) {
tamm::from_block_cyclic_tensor(ttensors.X_alpha, X_a);
}
#else
X_a = ttensors.X_alpha;
#endif

energy_qed_et = 0.0;
std::vector<double> polvec = {0.0, 0.0, 0.0};
for(int i = 0; i < sys_data.qed_nmodes; i++) {
polvec = scf_options.qed_polvecs[i];

// clang-format off
Scheduler{ec}
(ttensors.ehf_tmp("i", "j") = X_a("i", "k") * X_a("j", "k"))
(ttensors.ehf_tmp("i", "j") -= 0.5 * ttensors.D_last_alpha("i", "j"))
(ttensors.QED_Qxx("i", "j") = polvec[0] * ttensors.QED_Dx("i", "j"))
(ttensors.QED_Qxx("i", "j") += polvec[1] * ttensors.QED_Dy("i", "j"))
(ttensors.QED_Qxx("i", "j") += polvec[2] * ttensors.QED_Dz("i", "j"))
(ttensors.QED_Qyy("i", "j") = ttensors.QED_Qxx("i", "k") * ttensors.ehf_tmp("k", "j"))
(ttensors.ehf_tmp("i", "j") = ttensors.QED_Qyy("i", "k") * ttensors.QED_Qxx("k", "j"))
.execute();
// clang-format on

const double coupl_strength = pow(scf_options.qed_lambdas[i], 2);
energy_qed_et +=
0.5 * coupl_strength * tt_trace(ec, ttensors.D_last_alpha, ttensors.ehf_tmp);
}
#if defined(USE_SCALAPACK)
Tensor<double>::deallocate(X_a);
#endif
}
if(is_ks) { energy_2e += scf_vars.exc; }
}
if(is_uhf) {
nelectrons = tt_trace(ec, ttensors.D_last_alpha, ttensors.S1);
Expand All @@ -109,10 +59,6 @@ void SCFIO::print_energies(ExecutionContext& ec, const ChemEnv& chem_env, TAMMTe
std::cout << "1e energy N-e = " << NE_1e << endl;
std::cout << "1e energy = " << energy_1e << endl;
std::cout << "2e energy = " << energy_2e << std::endl;
if(is_qed) {
std::cout << "QED energy = " << energy_qed << std::endl;
std::cout << "QED eT energy = " << energy_qed_et << std::endl;
}
}
}

Expand Down
20 changes: 4 additions & 16 deletions exachem/scf/scf_tamm_tensors.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,22 +20,10 @@ class TAMMTensors {
Tensor<TensorType> ehf_tmp;
Tensor<TensorType> ehf_beta_tmp;

Tensor<TensorType> H1; // core hamiltonian
Tensor<TensorType> S1; // overlap ints
Tensor<TensorType> T1; // kinetic ints
Tensor<TensorType> V1; // nuclear ints
Tensor<TensorType> QED_Dx; // dipole ints
Tensor<TensorType> QED_Dy; // dipole ints
Tensor<TensorType> QED_Dz; // dipole ints
Tensor<TensorType> QED_Qxx; // quadrupole ints
Tensor<TensorType> QED_Qxy; // quadrupole ints
Tensor<TensorType> QED_Qxz; // quadrupole ints
Tensor<TensorType> QED_Qyy; // quadrupole ints
Tensor<TensorType> QED_Qyz; // quadrupole ints
Tensor<TensorType> QED_Qzz; // quadrupole ints
Tensor<TensorType> QED_1body;
Tensor<TensorType> QED_2body;
Tensor<TensorType> QED_energy;
Tensor<TensorType> H1; // core hamiltonian
Tensor<TensorType> S1; // overlap ints
Tensor<TensorType> T1; // kinetic ints
Tensor<TensorType> V1; // nuclear ints

Tensor<TensorType> X_alpha;
Tensor<TensorType> F_alpha; // H1+F_alpha_tmp
Expand Down

0 comments on commit e7faa91

Please sign in to comment.