From 6a5a6fc1902aae186854b31d1c6c72e954449956 Mon Sep 17 00:00:00 2001 From: Daniel Mejia-Rodriguez Date: Sun, 27 Oct 2024 14:52:34 -0700 Subject: [PATCH] QED HF --- exachem/scf/scf.cmake | 2 + exachem/scf/scf_hartree_fock.cpp | 55 ++++++ exachem/scf/scf_hartree_fock.hpp | 1 + exachem/scf/scf_iter.cpp | 1 + exachem/scf/scf_outputs.cpp | 82 +++++++- exachem/scf/scf_qed.cpp | 318 +++++++++++++++++++++++++++++++ exachem/scf/scf_qed.hpp | 33 ++++ exachem/scf/scf_tamm_tensors.hpp | 20 +- 8 files changed, 500 insertions(+), 12 deletions(-) create mode 100644 exachem/scf/scf_qed.cpp create mode 100644 exachem/scf/scf_qed.hpp diff --git a/exachem/scf/scf.cmake b/exachem/scf/scf.cmake index a01809d..83e86c5 100644 --- a/exachem/scf/scf.cmake +++ b/exachem/scf/scf.cmake @@ -4,6 +4,7 @@ include(TargetMacros) set(SCF_SRCDIR scf) set(SCF_INCLUDES + ${SCF_SRCDIR}/scf_qed.hpp ${SCF_SRCDIR}/scf_vars.hpp ${SCF_SRCDIR}/scf_iter.hpp ${SCF_SRCDIR}/scf_main.hpp @@ -21,6 +22,7 @@ set(SCF_INCLUDES ) set(SCF_SRCS + ${SCF_SRCDIR}/scf_qed.cpp ${SCF_SRCDIR}/scf_iter.cpp ${SCF_SRCDIR}/scf_main.cpp ${SCF_SRCDIR}/scf_gauxc.cpp diff --git a/exachem/scf/scf_hartree_fock.cpp b/exachem/scf/scf_hartree_fock.cpp index f94f812..88b7414 100644 --- a/exachem/scf/scf_hartree_fock.cpp +++ b/exachem/scf/scf_hartree_fock.cpp @@ -18,6 +18,7 @@ void exachem::scf::SCFHartreeFock::initialize(ExecutionContext& exc, ChemEnv& ch void exachem::scf::SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) { SCFCompute scf_compute; + SCFQed scf_qed; SCFIter scf_iter; SCFGuess scf_guess; SCFRestart scf_restart; @@ -336,6 +337,27 @@ void exachem::scf::SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_e /*** compute 1-e integrals ***/ /*** =========================== ***/ scf_compute.compute_hamiltonian(ec, scf_vars, chem_env, 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::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(ec, chem_env, scf_vars, ttensors); + if(chem_env.sys_data.do_qed) + scf_qed.compute_QED_1body(ec, chem_env, scf_vars, ttensors); + } if(chem_env.sys_data.has_ecp) { Tensor ECP{tAO, tAO}; @@ -653,6 +675,11 @@ void exachem::scf::SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_e 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(ec, chem_env, scf_vars, ttensors); + } + #if defined(USE_GAUXC) // Add snK contribution if(chem_env.sys_data.do_snK) { @@ -707,6 +734,7 @@ void exachem::scf::SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_e auto debug = chem_env.ioptions.scf_options.debug; if(rank == 0 && debug) std::cout << std::fixed << std::setprecision(2) << "xcf: " << xcf_time << "s, "; + if(sys_data.is_qed && !sys_data.do_qed) { scf_vars.eqed = gauxc_exc; } } ehf += gauxc_exc; @@ -770,6 +798,11 @@ void exachem::scf::SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_e 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(ec, chem_env, scf_vars, ttensors); + } + std::tie(ehf, rmsd) = scf_iter.scf_iter_body(ec, chem_env, scalapack_info, iter, scf_vars, ttensors, etensors #if defined(USE_GAUXC) @@ -920,6 +953,13 @@ void exachem::scf::SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_e scf_iter.compute_2bf(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(ec, chem_env, scf_vars, ttensors); + } } else if(scf_vars.lshift > 0) { // Remove level shift from Fock matrix @@ -1001,6 +1041,20 @@ void exachem::scf::SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_e Tensor::deallocate(ttensors.VXC_alpha); if(chem_env.sys_data.is_unrestricted) Tensor::deallocate(ttensors.VXC_beta); } + if(chem_env.sys_data.is_qed) { + scf_output.rw_mat_disk(ttensors.QED_Dx, files_prefix + ".QED_Dx", + chem_env.ioptions.scf_options.debug); + scf_output.rw_mat_disk(ttensors.QED_Dy, files_prefix + ".QED_Dy", + chem_env.ioptions.scf_options.debug); + scf_output.rw_mat_disk(ttensors.QED_Dz, files_prefix + ".QED_Dz", + chem_env.ioptions.scf_options.debug); + scf_output.rw_mat_disk(ttensors.QED_Qxx, files_prefix + ".QED_Qxx", + chem_env.ioptions.scf_options.debug); + Tensor::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 @@ -1064,6 +1118,7 @@ void exachem::scf::SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_e chem_env.sys_data.update(); if(rank == 0 && chem_env.ioptions.scf_options.debug) chem_env.sys_data.print(); + chem_env.sys_data.scf_energy = ehf; // iter not broadcasted, but fine since only rank 0 writes to json if(rank == 0) { chem_env.sys_data.results["output"]["SCF"]["final_energy"] = ehf; diff --git a/exachem/scf/scf_hartree_fock.hpp b/exachem/scf/scf_hartree_fock.hpp index e763ef3..1c00284 100644 --- a/exachem/scf/scf_hartree_fock.hpp +++ b/exachem/scf/scf_hartree_fock.hpp @@ -14,6 +14,7 @@ #include "exachem/scf/scf_compute.hpp" #include "exachem/scf/scf_iter.hpp" #include "exachem/scf/scf_outputs.hpp" +#include "exachem/scf/scf_qed.hpp" #include "exachem/scf/scf_restart.hpp" #include "exachem/scf/scf_taskmap.hpp" #include diff --git a/exachem/scf/scf_iter.cpp b/exachem/scf/scf_iter.cpp index 9c90cec..8f403b0 100644 --- a/exachem/scf/scf_iter.cpp +++ b/exachem/scf/scf_iter.cpp @@ -109,6 +109,7 @@ std::tuple exachem::scf::SCFIter::scf_iter_body( std::chrono::duration_cast>((xcf_stop - xcf_start)).count(); if(rank == 0 && debug) std::cout << std::fixed << std::setprecision(2) << "xcf: " << xcf_time << "s, "; + if(sys_data.is_qed && !sys_data.do_qed) { scf_vars.eqed = gauxc_exc; } } ehf += gauxc_exc; diff --git a/exachem/scf/scf_outputs.cpp b/exachem/scf/scf_outputs.cpp index bb6548e..b54ca13 100644 --- a/exachem/scf/scf_outputs.cpp +++ b/exachem/scf/scf_outputs.cpp @@ -26,18 +26,22 @@ double exachem::scf::SCFIO::tt_trace(ExecutionContext& ec, Tensor& T void exachem::scf::SCFIO::print_energies(ExecutionContext& ec, ChemEnv& chem_env, TAMMTensors& ttensors, EigenTensors& etensors, SCFVars& scf_vars, ScalapackInfo& scalapack_info) { - const SystemData& sys_data = chem_env.sys_data; + const SystemData& sys_data = chem_env.sys_data; + const SCFOptions& scf_options = chem_env.ioptions.scf_options; 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 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; if(is_rhf) { nelectrons = tt_trace(ec, ttensors.D_last_alpha, ttensors.S1); kinetic_1e = tt_trace(ec, ttensors.D_last_alpha, ttensors.T1); @@ -45,7 +49,62 @@ void exachem::scf::SCFIO::print_energies(ExecutionContext& ec, ChemEnv& chem_env 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) { energy_2e += scf_vars.exc; } + 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 X_a; + +#if defined(USE_SCALAPACK) + X_a = {scf_vars.tAO, scf_vars.tAO_ortho}; + Tensor::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 polvec = {0.0, 0.0, 0.0}; + Scheduler sch{ec}; + Tensor ehf_tmp{scf_vars.tAO, scf_vars.tAO}; + Tensor QED_Qxx{scf_vars.tAO, scf_vars.tAO}; + Tensor QED_Qyy{scf_vars.tAO, scf_vars.tAO}; + sch.allocate(ehf_tmp, QED_Qxx, QED_Qyy).execute(); + + for(int i = 0; i < sys_data.qed_nmodes; i++) { + polvec = scf_options.qed_polvecs[i]; + + // clang-format off + sch + (ehf_tmp("i", "j") = X_a("i", "k") * X_a("j", "k")) + (ehf_tmp("i", "j") -= 0.5 * ttensors.D_last_alpha("i", "j")) + (QED_Qxx("i", "j") = polvec[0] * ttensors.QED_Dx("i", "j")) + (QED_Qxx("i", "j") += polvec[1] * ttensors.QED_Dy("i", "j")) + (QED_Qxx("i", "j") += polvec[2] * ttensors.QED_Dz("i", "j")) + (QED_Qyy("i", "j") = QED_Qxx("i", "k") * ehf_tmp("k", "j")) + (ehf_tmp("i", "j") = QED_Qyy("i", "k") * 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); + } + sch.deallocate(ehf_tmp, QED_Qxx, QED_Qyy).execute(); + +#if defined(USE_SCALAPACK) + Tensor::deallocate(X_a); +#endif + } } if(is_uhf) { nelectrons = tt_trace(ec, ttensors.D_last_alpha, ttensors.S1); @@ -72,6 +131,13 @@ void exachem::scf::SCFIO::print_energies(ExecutionContext& ec, ChemEnv& chem_env chem_env.sys_data.results["output"]["SCF"]["kinetic_1e"] = kinetic_1e; chem_env.sys_data.results["output"]["SCF"]["energy_1e"] = energy_1e; chem_env.sys_data.results["output"]["SCF"]["energy_2e"] = energy_2e; + + if(is_qed) { + std::cout << "QED energy = " << energy_qed << std::endl; + std::cout << "QED eT energy = " << energy_qed_et << std::endl; + chem_env.sys_data.results["output"]["SCF"]["energy_qed"] = energy_qed; + chem_env.sys_data.results["output"]["SCF"]["energy_qed_et"] = energy_qed_et; + } } } diff --git a/exachem/scf/scf_qed.cpp b/exachem/scf/scf_qed.cpp new file mode 100644 index 0000000..de63b0a --- /dev/null +++ b/exachem/scf/scf_qed.cpp @@ -0,0 +1,318 @@ +/* + * ExaChem: Open Source Exascale Computational Chemistry Software. + * + * Copyright 2023-2024 Pacific Northwest National Laboratory, Battelle Memorial Institute. + * + * See LICENSE.txt for details + */ + +#include "exachem/scf/scf_qed.hpp" + +void exachem::scf::SCFQed::qed_functionals_setup(std::vector& params, ChemEnv& chem_env) { + SCFOptions& scf_options = chem_env.ioptions.scf_options; + + params[0] = 1.0; + // SCFOptions& scf_options = scf_options; + + std::copy(scf_options.qed_lambdas.begin(), scf_options.qed_lambdas.end(), params.begin() + 1); + if(!scf_options.qed_omegas.empty()) { + std::copy(scf_options.qed_omegas.begin(), scf_options.qed_omegas.end(), params.begin() + 1025); + } + std::ofstream libxc_params("./libxc.params"); + libxc_params << std::fixed << std::setprecision(16); + for(auto const& val: params) libxc_params << val << std::endl; + libxc_params.close(); +} + +template +void exachem::scf::SCFQed::compute_QED_1body(ExecutionContext& ec, ChemEnv& chem_env, + const SCFVars& scf_vars, TAMMTensors& ttensors) { + auto [mu, nu] = scf_vars.tAO.labels<2>("all"); + Scheduler sch{ec}; + + SystemData& sys_data = chem_env.sys_data; + SCFOptions& scf_options = chem_env.ioptions.scf_options; + std::vector& atoms = chem_env.atoms; + + const int nmodes = sys_data.qed_nmodes; + const std::vector lambdas = scf_options.qed_lambdas; + const std::vector> polvecs = scf_options.qed_polvecs; + + // compute nuclear dipole operator + std::vector mu_nuc = {0.0, 0.0, 0.0}; + for(size_t i = 0; i < atoms.size(); i++) { + mu_nuc[0] += atoms[i].x * atoms[i].atomic_number; + mu_nuc[1] += atoms[i].y * atoms[i].atomic_number; + mu_nuc[2] += atoms[i].z * atoms[i].atomic_number; + } + + double s0_scaling = 0.0; + for(int i = 0; i < nmodes; i++) { + s0_scaling += + lambdas[i] * + (mu_nuc[0] * polvecs[i][0] + mu_nuc[1] * polvecs[i][1] + mu_nuc[2] * polvecs[i][2]) / + sys_data.nelectrons; + } + s0_scaling *= 0.5 * s0_scaling; + + // clang-format off + sch + (ttensors.QED_1body(mu,nu) = s0_scaling*ttensors.S1(mu,nu)) + (ttensors.QED_1body(mu,nu) += ttensors.QED_Qxx(mu,nu)); + // clang-format on + + for(int i = 0; i < nmodes; i++) { + double s1_scaling = + (mu_nuc[0] * polvecs[i][0] + mu_nuc[1] * polvecs[i][1] + mu_nuc[2] * polvecs[i][2]) * + lambdas[i] * lambdas[i] / sys_data.nelectrons; + + // clang-format off + sch + (ttensors.QED_2body(mu,nu) = polvecs[i][0] * ttensors.QED_Dx(mu,nu)) + (ttensors.QED_2body(mu,nu) += polvecs[i][1] * ttensors.QED_Dy(mu,nu)) + (ttensors.QED_2body(mu,nu) += polvecs[i][2] * ttensors.QED_Dz(mu,nu)) + (ttensors.QED_1body(mu,nu) -= s1_scaling * ttensors.QED_2body(mu,nu)); + // clang-format on + } + + // clang-format off + sch + (ttensors.H1(mu,nu) += ttensors.QED_1body(mu,nu)) + .execute(); + // clang-format on +} + +template +void exachem::scf::SCFQed::compute_QED_2body(ExecutionContext& ec, ChemEnv& chem_env, + const SCFVars& scf_vars, TAMMTensors& ttensors) { + auto [mu, nu, ku] = scf_vars.tAO.labels<3>("all"); + Tensor tensor{ttensors.QED_1body.tiled_index_spaces()}; //{tAO, tAO}; + Scheduler sch{ec}; + + auto& atoms = chem_env.atoms; + SystemData& sys_data = chem_env.sys_data; + SCFOptions& scf_options = chem_env.ioptions.scf_options; + + const int nmodes = sys_data.qed_nmodes; + const std::vector lambdas = scf_options.qed_lambdas; + const std::vector> polvecs = scf_options.qed_polvecs; + + // compute nuclear dipole operator + std::vector mu_nuc = {0.0, 0.0, 0.0}; + for(size_t i = 0; i < atoms.size(); i++) { + mu_nuc[0] += atoms[i].x * atoms[i].atomic_number; + mu_nuc[1] += atoms[i].y * atoms[i].atomic_number; + mu_nuc[2] += atoms[i].z * atoms[i].atomic_number; + } + + // clang-format off + sch + .allocate(tensor) + (ttensors.QED_2body(mu,nu) = 0.0); + // clang-format on + + for(int i = 0; i < nmodes; i++) { + double mu_nuc_ope = 0.0; + mu_nuc_ope = + (mu_nuc[0] * polvecs[i][0] + mu_nuc[1] * polvecs[i][1] + mu_nuc[2] * polvecs[i][2]) * + lambdas[i] / sys_data.nelectrons; + + double t1 = -0.5 * pow(lambdas[i], 2); + double t2 = 0.5 * mu_nuc_ope * lambdas[i]; + double t3 = -0.5 * pow(mu_nuc_ope, 2); + + // clang-format off + sch + (tensor() = polvecs[i][0]*ttensors.QED_Dx()) + (tensor() += polvecs[i][1]*ttensors.QED_Dy()) + (tensor() += polvecs[i][2]*ttensors.QED_Dz()) + (ttensors.ehf_tmp(mu,nu) = tensor(mu,ku)*ttensors.D_last_alpha(ku,nu)) + (ttensors.QED_2body(mu,nu) += t1*ttensors.ehf_tmp(mu,ku)*tensor(ku,nu)) + (ttensors.QED_2body(mu,nu) += t2*ttensors.ehf_tmp(mu,ku)*ttensors.S1(ku,nu)) + (ttensors.QED_2body(mu,nu) += t2*ttensors.S1(mu,ku)*ttensors.ehf_tmp(nu,ku)) + (ttensors.ehf_tmp(mu,nu) = t3*ttensors.S1(mu,ku)*ttensors.D_alpha(ku,nu)) + (ttensors.QED_2body(mu,nu) += ttensors.ehf_tmp(mu,ku)*ttensors.S1(ku,nu)); + // clang-format on + } + + // clang-format off + sch + (ttensors.F_alpha(mu,nu) += ttensors.QED_2body(mu,nu)) + .deallocate(tensor) + .execute(); + // clang-format on +} + +template +void exachem::scf::SCFQed::compute_qed_emult_ints(ExecutionContext& ec, ChemEnv& chem_env, + const SCFVars& spvars, TAMMTensors& ttensors) { + using libint2::Atom; + using libint2::BasisSet; + using libint2::Engine; + using libint2::Operator; + using libint2::Shell; + + // auto& atoms = chem_env.atoms; + SystemData& sys_data = chem_env.sys_data; + SCFOptions& scf_options = chem_env.ioptions.scf_options; + const libint2::BasisSet& shells = chem_env.shells; + + const std::vector& AO_tiles = spvars.AO_tiles; + const std::vector& shell_tile_map = spvars.shell_tile_map; + + const int nmodes = sys_data.qed_nmodes; + const std::vector lambdas = scf_options.qed_lambdas; + const std::vector> polvecs = scf_options.qed_polvecs; + + Engine engine(Operator::emultipole2, max_nprim(shells), max_l(shells), 0); + + // engine.set(otype); + + // auto& buf = (engine.results()); + + auto compute_qed_emult_ints_lambda = [&](const IndexVector& blockid) { + auto bi0 = blockid[0]; + auto bi1 = blockid[1]; + + const TAMM_SIZE size = ttensors.QED_Dx.block_size(blockid); + auto block_dims = ttensors.QED_Dx.block_dims(blockid); + std::vector dbufx(size); + std::vector dbufy(size); + std::vector dbufz(size); + std::vector dbufQ(size); + + auto bd1 = block_dims[1]; + + // cout << "blockid: [" << blockid[0] <<"," << blockid[1] << "], dims(0,1) = " << + // block_dims[0] << ", " << block_dims[1] << endl; + + // auto s1 = blockid[0]; + auto s1range_end = shell_tile_map[bi0]; + decltype(s1range_end) s1range_start = 0l; + if(bi0 > 0) s1range_start = shell_tile_map[bi0 - 1] + 1; + + // cout << "s1-start,end = " << s1range_start << ", " << s1range_end << endl; + for(auto s1 = s1range_start; s1 <= s1range_end; ++s1) { + // auto bf1 = shell2bf[s1]; //shell2bf[s1]; // first basis function in + // this shell + auto n1 = shells[s1].size(); + + auto s2range_end = shell_tile_map[bi1]; + decltype(s2range_end) s2range_start = 0l; + if(bi1 > 0) s2range_start = shell_tile_map[bi1 - 1] + 1; + + // cout << "s2-start,end = " << s2range_start << ", " << s2range_end << endl; + + // cout << "screend shell pair list = " << s2spl << endl; + for(auto s2 = s2range_start; s2 <= s2range_end; ++s2) { + // for (auto s2: spvars.obs_shellpair_list.at(s1)) { + // auto s2 = blockid[1]; + // if (s2>s1) continue; + + if(s2 > s1) { + auto s2spl = spvars.obs_shellpair_list.at(s2); + if(std::find(s2spl.begin(), s2spl.end(), s1) == s2spl.end()) continue; + } + else { + auto s2spl = spvars.obs_shellpair_list.at(s1); + if(std::find(s2spl.begin(), s2spl.end(), s2) == s2spl.end()) continue; + } + + // auto bf2 = shell2bf[s2]; + auto n2 = shells[s2].size(); + + std::vector tbufX(n1 * n2); + std::vector tbufY(n1 * n2); + std::vector tbufZ(n1 * n2); + std::vector tbufXX(n1 * n2); + std::vector tbufXY(n1 * n2); + std::vector tbufXZ(n1 * n2); + std::vector tbufYY(n1 * n2); + std::vector tbufYZ(n1 * n2); + std::vector tbufZZ(n1 * n2); + + // compute shell pair; return is the pointer to the buffer + const auto& buf = engine.compute(shells[s1], shells[s2]); + EXPECTS(buf.size() >= 10); + if(buf[0] == nullptr) continue; + // "map" buffer to a const Eigen Matrix, and copy it to the + // corresponding blocks of the result + + // cout << buf[1].size() << endl; + // cout << buf[2].size() << endl; + // cout << buf[3].size() << endl; + Eigen::Map buf_mat_X(buf[1], n1, n2); + Eigen::Map(&tbufX[0], n1, n2) = buf_mat_X; + Eigen::Map buf_mat_Y(buf[2], n1, n2); + Eigen::Map(&tbufY[0], n1, n2) = buf_mat_Y; + Eigen::Map buf_mat_Z(buf[3], n1, n2); + Eigen::Map(&tbufZ[0], n1, n2) = buf_mat_Z; + Eigen::Map buf_mat_XX(buf[4], n1, n2); + Eigen::Map(&tbufXX[0], n1, n2) = buf_mat_XX; + Eigen::Map buf_mat_XY(buf[5], n1, n2); + Eigen::Map(&tbufXY[0], n1, n2) = buf_mat_XY; + Eigen::Map buf_mat_XZ(buf[6], n1, n2); + Eigen::Map(&tbufXZ[0], n1, n2) = buf_mat_XZ; + Eigen::Map buf_mat_YY(buf[7], n1, n2); + Eigen::Map(&tbufYY[0], n1, n2) = buf_mat_YY; + Eigen::Map buf_mat_YZ(buf[8], n1, n2); + Eigen::Map(&tbufYZ[0], n1, n2) = buf_mat_YZ; + Eigen::Map buf_mat_ZZ(buf[9], n1, n2); + Eigen::Map(&tbufZZ[0], n1, n2) = buf_mat_ZZ; + + // cout << "buf_mat_X :" << buf_mat_X << endl; + // cout << "buf_mat_Y :" << buf_mat_Y << endl; + // cout << "buf_mat_Z :" << buf_mat_Z << endl; + + auto curshelloffset_i = 0U; + auto curshelloffset_j = 0U; + for(auto x = s1range_start; x < s1; x++) curshelloffset_i += AO_tiles[x]; + for(auto x = s2range_start; x < s2; x++) curshelloffset_j += AO_tiles[x]; + + size_t c = 0; + auto dimi = curshelloffset_i + AO_tiles[s1]; + auto dimj = curshelloffset_j + AO_tiles[s2]; + + for(size_t i = curshelloffset_i; i < dimi; i++) { + for(size_t j = curshelloffset_j; j < dimj; j++, c++) { + dbufx[i * bd1 + j] = tbufX[c]; + dbufy[i * bd1 + j] = tbufY[c]; + dbufz[i * bd1 + j] = tbufZ[c]; + dbufQ[i * bd1 + j] = 0.0; + for(int k = 0; k < nmodes; k++) { + dbufQ[i * bd1 + j] += + 0.5 * pow(lambdas[k], 2) * + (tbufXX[c] * pow(polvecs[k][0], 2) + tbufYY[c] * pow(polvecs[k][1], 2) + + tbufZZ[c] * pow(polvecs[k][2], 2) + + 2.0 * tbufXY[c] * polvecs[k][0] * polvecs[k][1] + + 2.0 * tbufXZ[c] * polvecs[k][0] * polvecs[k][2] + + 2.0 * tbufYZ[c] * polvecs[k][1] * polvecs[k][2]); + } + } + } + } // s2 + } // s1 + + ttensors.QED_Dx.put(blockid, dbufx); + ttensors.QED_Dy.put(blockid, dbufy); + ttensors.QED_Dz.put(blockid, dbufz); + ttensors.QED_Qxx.put(blockid, dbufQ); + }; + + block_for(ec, ttensors.QED_Dx(), compute_qed_emult_ints_lambda); +} + +template void exachem::scf::SCFQed::compute_QED_1body(ExecutionContext& ec, + ChemEnv& chem_env, + const SCFVars& scf_vars, + TAMMTensors& ttensors); + +template void exachem::scf::SCFQed::compute_QED_2body(ExecutionContext& ec, + ChemEnv& chem_env, + const SCFVars& scf_vars, + TAMMTensors& ttensors); + +template void exachem::scf::SCFQed::compute_qed_emult_ints(ExecutionContext& ec, + ChemEnv& chem_env, + const SCFVars& spvars, + TAMMTensors& ttensors); diff --git a/exachem/scf/scf_qed.hpp b/exachem/scf/scf_qed.hpp new file mode 100644 index 0000000..8a6f910 --- /dev/null +++ b/exachem/scf/scf_qed.hpp @@ -0,0 +1,33 @@ +/* + * ExaChem: Open Source Exascale Computational Chemistry Software. + * + * Copyright 2023-2024 Pacific Northwest National Laboratory, Battelle Memorial Institute. + * + * See LICENSE.txt for details + */ + +#pragma once + +#include + +#include "exachem/scf/scf_common.hpp" + +namespace exachem::scf { + +class SCFQed { +public: + void qed_functionals_setup(std::vector& params, ChemEnv& chem_env); + + template + void compute_QED_1body(ExecutionContext& ec, ChemEnv& chem_env, const SCFVars& scf_vars, + TAMMTensors& ttensors); + + template + void compute_QED_2body(ExecutionContext& ec, ChemEnv& chem_env, const SCFVars& scf_vars, + TAMMTensors& ttensors); + + template + void compute_qed_emult_ints(ExecutionContext& ec, ChemEnv& chem_env, const SCFVars& spvars, + TAMMTensors& ttensors); +}; +} // namespace exachem::scf diff --git a/exachem/scf/scf_tamm_tensors.hpp b/exachem/scf/scf_tamm_tensors.hpp index bacbeb2..5da58dd 100644 --- a/exachem/scf/scf_tamm_tensors.hpp +++ b/exachem/scf/scf_tamm_tensors.hpp @@ -30,10 +30,22 @@ class TAMMTensors { Tensor ehf_tmp; Tensor ehf_beta_tmp; - Tensor H1; // core hamiltonian - Tensor S1; // overlap ints - Tensor T1; // kinetic ints - Tensor V1; // nuclear ints + Tensor H1; // core hamiltonian + Tensor S1; // overlap ints + Tensor T1; // kinetic ints + Tensor V1; // nuclear ints + Tensor QED_Dx; // dipole ints + Tensor QED_Dy; // dipole ints + Tensor QED_Dz; // dipole ints + Tensor QED_Qxx; // quadrupole ints + Tensor QED_Qxy; // quadrupole ints + Tensor QED_Qxz; // quadrupole ints + Tensor QED_Qyy; // quadrupole ints + Tensor QED_Qyz; // quadrupole ints + Tensor QED_Qzz; // quadrupole ints + Tensor QED_1body; + Tensor QED_2body; + Tensor QED_energy; Tensor X_alpha; Tensor F_alpha; // H1+F_alpha_tmp