diff --git a/exachem/scf/scf_hartree_fock.cpp b/exachem/scf/scf_hartree_fock.cpp index df51913..9b52ebc 100644 --- a/exachem/scf/scf_hartree_fock.cpp +++ b/exachem/scf/scf_hartree_fock.cpp @@ -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; @@ -301,27 +300,6 @@ void SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) { /*** =========================== ***/ scf_compute.compute_hamiltonian(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::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}; @@ -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(ec, chem_env, scf_vars, ttensors); - } - TensorType gauxc_exc = 0.; #if defined(USE_GAUXC) if(chem_env.sys_data.is_ks) { @@ -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(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) @@ -845,13 +812,6 @@ void SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) { 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 @@ -930,18 +890,6 @@ void SCFHartreeFock::scf_hf(ExecutionContext& exc, ChemEnv& chem_env) { 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); - 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 diff --git a/exachem/scf/scf_hartree_fock.hpp b/exachem/scf/scf_hartree_fock.hpp index 7f49983..069fb4c 100644 --- a/exachem/scf/scf_hartree_fock.hpp +++ b/exachem/scf/scf_hartree_fock.hpp @@ -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 diff --git a/exachem/scf/scf_outputs.cpp b/exachem/scf/scf_outputs.cpp index aff57bf..3a83f10 100644 --- a/exachem/scf/scf_outputs.cpp +++ b/exachem/scf/scf_outputs.cpp @@ -23,16 +23,13 @@ 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); @@ -40,54 +37,7 @@ void SCFIO::print_energies(ExecutionContext& ec, const ChemEnv& chem_env, TAMMTe 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 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}; - 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::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); @@ -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; - } } } diff --git a/exachem/scf/scf_tamm_tensors.hpp b/exachem/scf/scf_tamm_tensors.hpp index 18f2f95..837cbf3 100644 --- a/exachem/scf/scf_tamm_tensors.hpp +++ b/exachem/scf/scf_tamm_tensors.hpp @@ -20,22 +20,10 @@ 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 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 H1; // core hamiltonian + Tensor S1; // overlap ints + Tensor T1; // kinetic ints + Tensor V1; // nuclear ints Tensor X_alpha; Tensor F_alpha; // H1+F_alpha_tmp