From 072796c2725878d73acc1ece743979fe6d96ff6b Mon Sep 17 00:00:00 2001 From: fbischoff Date: Tue, 13 Aug 2024 11:06:52 -0400 Subject: [PATCH] fixing expectation value of using the singles --- src/madness/chem/BSHApply.h | 9 +++-- src/madness/chem/CC2.h | 57 +++++++++++++++++--------------- src/madness/chem/CCPotentials.cc | 14 ++++---- src/madness/chem/CCPotentials.h | 6 ++-- 4 files changed, 48 insertions(+), 38 deletions(-) diff --git a/src/madness/chem/BSHApply.h b/src/madness/chem/BSHApply.h index f2038877db6..f3f5dc727d7 100644 --- a/src/madness/chem/BSHApply.h +++ b/src/madness/chem/BSHApply.h @@ -29,7 +29,7 @@ class BSHApply { double levelshift=0.0; double lo=1.e-6; double bshtol=1.e-5; - bool printme=false; + bool printme=true; bool destroy_Vpsi=false; Function metric; return_value ret_value=residual; // return the new orbitals/functions or the residuals @@ -64,6 +64,7 @@ class BSHApply { std::vector < std::shared_ptr > > ops(psi.size()); for (int i=0; i >( BSHOperatorPtr(world, sqrt(-2.0*eps_in_green(e_i)), lo, bshtol)); ops[i]->destructive()=true; @@ -130,7 +131,7 @@ class BSHApply { const Tensor fock1) const { // check dimensions - bool consistent=(psi.size()==size_t(fock1.dim(0))); + bool consistent=(psi.size()==size_t(fock1.dim(0))); if ((fock1.ndim()==2) and not (psi.size()==size_t(fock1.dim(1)))) consistent=false; if (not consistent) { @@ -150,6 +151,10 @@ class BSHApply { for (int i=0; i bsh_apply(world); -// bsh_apply.ret_value=BSHApply::update; // return the new singles functions, not the residual -// bsh_apply.metric=info.R_square; -// -// // coupling between singles involves the active fock matrix shifted by the excitation energy -// auto nfreeze=info.parameters.freeze(); -// Tensor fock=info.fock(Slice(nfreeze,-1),Slice(nfreeze,-1)); -// for (int i=0; i > > G(singles.size()); for (size_t i = 0; i < G.size(); i++) { const double bsh_eps = info.orbital_energies[i + info.parameters.freeze()] + omega; @@ -309,14 +301,11 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface { BSHOperatorPtr3D(world, sqrt(-2.0 * bsh_eps), info.parameters.lo(), info.parameters.thresh_bsh_3D())); } world.gop.fence(); - time_makebsh.info(); -// -// // apply bsh operators + + // apply bsh operators CCTimer time_applyG(world, "Apply G-Operators"); - // auto [GV, energy_update] = bsh_apply(singles.get_vecfunction(), fock, V); vector_real_function_3d GV = apply, double, 3>(world, G, V); -// world.gop.fence(); - // auto GV=res-singles.get_vecfunction(); + world.gop.fence(); time_applyG.info(); // apply Q-Projector to result @@ -513,9 +502,23 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface { for (auto& tmp_pair : vpairs) pairs.insert(tmp_pair.i, tmp_pair.j, tmp_pair); auto ccpair2function = [](const CCPair& a) {return a.function();}; return compute_local_coupling(pairs.convert(pairs,ccpair2function), info); - }; + /// compute the coupling of singles function if orbitals are localized + + /// @return the coupling terms c_i = -\sum_(j\neq i) f_ij |\phi_j> (for whatever phi is) + static std::vector compute_local_coupling(const std::vector& singles, + const Info& info) { + + MADNESS_CHECK_THROW(singles.size()>0,"compute_local_coupling: singles vector is empty"); + World& world=singles.front().world(); + auto active=Slice(info.parameters.freeze(),-1); + Tensor Fact=info.fock(active,active); + for (int i=0; i + \sum_{l\neq j} f_lj |u_il> diff --git a/src/madness/chem/CCPotentials.cc b/src/madness/chem/CCPotentials.cc index 92d5500ba7c..e524b0ed4bb 100644 --- a/src/madness/chem/CCPotentials.cc +++ b/src/madness/chem/CCPotentials.cc @@ -486,10 +486,16 @@ CCPotentials::compute_kinetic_energy(World& world, const vector_real_function_3d return kinetic; } + double CCPotentials::compute_cis_expectation_value(World& world, const CC_vecfunction& x, const vector_real_function_3d& V, const bool print, const Info& info) { + // following eq. (34) of the CIS paper Kottmann et al, PCCP, 17, 31453, (2015) + // doi: https://doi.org/10.1039/C5CP00345H + // the expectation value of the CIS wave function is computed by projecting the + // CIS wave function onto eq. (22) + // the potential V must contain the coupling term when using localized orbitals const vector_real_function_3d xbra = info.R_square*(x.get_vecfunction()); const vector_real_function_3d xket = x.get_vecfunction(); const double kinetic = compute_kinetic_energy(world, xbra, xket); @@ -2879,7 +2885,7 @@ CCPotentials::get_CC2_singles_potential_gs(World& world, const CC_vecfunction& s } madness::vector_real_function_3d -CCPotentials::get_CCS_potential_ex(World& world, CC_vecfunction& x, const bool print, Info& info) { +CCPotentials::get_CCS_potential_ex(World& world, const CC_vecfunction& x, const bool print, Info& info) { if (x.type != RESPONSE) error("get_CCS_response_potential: Wrong type of input singles"); Pairs empty_doubles; @@ -2895,14 +2901,12 @@ CCPotentials::get_CCS_potential_ex(World& world, CC_vecfunction& x, const bool p info.intermediate_potentials.insert(copy(world, potential), x, POT_singles_); vector_real_function_3d result = add(world, fock_residue, potential); truncate(world, result); - const double omega = compute_cis_expectation_value(world, x, result, print, info); - x.omega = omega; return result; } madness::vector_real_function_3d CCPotentials::get_CC2_singles_potential_ex(World& world, const CC_vecfunction& gs_singles, - const Pairs& gs_doubles, CC_vecfunction& ex_singles, + const Pairs& gs_doubles, const CC_vecfunction& ex_singles, const Pairs& response_doubles, Info& info) { MADNESS_ASSERT(gs_singles.type == PARTICLE); @@ -2949,8 +2953,6 @@ CCPotentials::get_CC2_singles_potential_ex(World& world, const CC_vecfunction& g info.intermediate_potentials.insert(copy(world, potential), ex_singles, POT_singles_); vector_real_function_3d result = add(world, fock_residue, potential); truncate(world, result); - const double omega = compute_cis_expectation_value(world, ex_singles, result, true, info); - ex_singles.omega = omega; return result; } diff --git a/src/madness/chem/CCPotentials.h b/src/madness/chem/CCPotentials.h index e2f3ccd6b84..4fa36bcb859 100644 --- a/src/madness/chem/CCPotentials.h +++ b/src/madness/chem/CCPotentials.h @@ -276,7 +276,7 @@ class CCPotentials { static double compute_kinetic_energy(World& world, const vector_real_function_3d& xbra, const vector_real_function_3d& xket); - /// returns \f$ + \f$ + /// compute the expectation value excitation energy using the CIS/CCS/CC2 singles static double compute_cis_expectation_value(World& world, const CC_vecfunction& x, const vector_real_function_3d& V, const bool print, const Info& info); @@ -681,13 +681,13 @@ class CCPotentials { /// the V part is stored in the intermediate_potentials structure /// the expectation value is calculated and updated static vector_real_function_3d - get_CCS_potential_ex(World& world, CC_vecfunction& x, const bool print, Info& info); + get_CCS_potential_ex(World& world, const CC_vecfunction& x, const bool print, Info& info); /// Calculates the CC2 singles potential for the Excited state: result = Fock_residue + V /// the V part is stored in the intermediate_potentials structure static vector_real_function_3d get_CC2_singles_potential_ex(World& world, const CC_vecfunction& gs_singles, - const Pairs& gs_doubles, CC_vecfunction& ex_singles, + const Pairs& gs_doubles, const CC_vecfunction& ex_singles, const Pairs& response_doubles, Info& info); /// Calculates the CC2 singles potential for the Excited state: result = Fock_residue + V