From df7d831ff5861bc583e34ac48493308fb7eb515d Mon Sep 17 00:00:00 2001 From: fbischoff Date: Fri, 1 Nov 2024 13:04:02 +0100 Subject: [PATCH] LRCC2 works for He, singles also for water dimer, but a tensor exception occurs at iterating doubles --- src/apps/cc2/cc2.cc | 6 +- src/madness/chem/CC2.cc | 178 +++++++++++++------------------ src/madness/chem/CC2.h | 13 +-- src/madness/chem/CCPotentials.cc | 133 +++++++++++++---------- src/madness/chem/CCPotentials.h | 15 +-- src/madness/chem/CCStructures.cc | 74 ++++++------- src/madness/chem/CCStructures.h | 25 ++++- 7 files changed, 228 insertions(+), 216 deletions(-) diff --git a/src/apps/cc2/cc2.cc b/src/apps/cc2/cc2.cc index 8fe169735ae..5eca090a897 100644 --- a/src/apps/cc2/cc2.cc +++ b/src/apps/cc2/cc2.cc @@ -97,7 +97,11 @@ int main(int argc, char **argv) { nemo->molecule().print(); } double hf_energy = nemo->value(); - cc2.solve(); + try { + cc2.solve(); + } catch (std::exception& e) { + print("Caught exception: ", e.what()); + } if (world.rank() == 0) printf("\nfinished at time %.1fs\n\n", wall_time()); world.gop.fence(); diff --git a/src/madness/chem/CC2.cc b/src/madness/chem/CC2.cc index e7c3e0bb7bc..88bee4e8f7c 100644 --- a/src/madness/chem/CC2.cc +++ b/src/madness/chem/CC2.cc @@ -38,6 +38,15 @@ CC2::solve() { info=CCOPS.update_info(parameters,nemo); info.intermediate_potentials=CCIntermediatePotentials(parameters); + // check if all pair functions have been loaded or computed + auto all_pairs_exist = [](const Pairs& pairs) { + for (const auto& p : pairs.allpairs) { + if (not p.second.function().is_initialized()) return false; + if (not p.second.constant_part.is_initialized()) return false; + } + return true; + }; + // doubles for ground state Pairs mp2pairs, cc2pairs; // singles for ground state @@ -56,9 +65,8 @@ CC2::solve() { // check for restart data for CC2, otherwise use MP2 as guess if (need_cc2) { - Pairs dummypairs; - bool found_cc2d = initialize_pairs(dummypairs, GROUND_STATE, CT_CC2, cc2singles, CC_vecfunction(RESPONSE), 0, info); - if (not found_cc2d) need_mp2=true; + initialize_pairs(cc2pairs, GROUND_STATE, CT_CC2, cc2singles, CC_vecfunction(RESPONSE), 0, info); + if (not all_pairs_exist(cc2pairs)) need_mp2=true; } if (need_tdhf) { @@ -85,12 +93,11 @@ CC2::solve() { if (need_cc2) { // check if singles or/and doubles to restart are there - cc2singles=initialize_singles(CT_CC2,PARTICLE, true); - const bool load_doubles = initialize_pairs(cc2pairs, GROUND_STATE, CT_CC2, cc2singles, CC_vecfunction(RESPONSE), 0, info); + cc2singles=initialize_singles(world,CT_CC2, PARTICLE); - // nothing to restart -> make MP2 - if (not load_doubles) { - // use mp2 as cc2 guess + // use mp2 as cc2 guess + if (not all_pairs_exist(cc2pairs)) { + if (world.rank()==0) print("could not find CC2 pairs -- using MP2 as guess"); for (auto& tmp:mp2pairs.allpairs) { const size_t i = tmp.second.i; const size_t j = tmp.second.j; @@ -98,8 +105,12 @@ CC2::solve() { } } - cc2_energy = solve_cc2(cc2singles, cc2pairs, info); - output_calc_info_schema("cc2",cc2_energy); + if (parameters.no_compute_cc2()) { + if (world.rank()==0) print("found no_compute_cc2 key -- no recomputation of singles or doubles or energy"); + } else { + cc2_energy = solve_cc2(cc2singles, cc2pairs, info); + output_calc_info_schema("cc2",cc2_energy); + } output.section(assign_name(CT_CC2) + " Calculation Ended !"); if (world.rank() == 0) { @@ -253,16 +264,26 @@ CC2::solve() { } else if (ctype == CT_LRCC2) { CCTimer time(world, "Whole LRCC2 Calculation"); - auto vccs=solve_ccs(); - - if (world.rank()==0) print_header3("reiterating CCS"); + if (world.rank() == 0) print_header3("reiterating CCS"); + auto vccs = solve_ccs(); iterate_ccs_singles(vccs[0], info); - if (world.rank()==0) print_header3("end reiterating CCS"); + if (world.rank() == 0) print_header3("end reiterating CCS"); + + // for solving the LRCC2 equations we need converged CC2 singles and the applied potential + bool need_reiterate_cc2=(cc2singles.size()==0) + or (not info.intermediate_potentials.potential_exists(cc2singles, POT_singles_)); + if (need_reiterate_cc2) { + if (world.rank()==0) print_header3("reiterating CC2 singles for intermediate potentials"); + iterate_cc2_singles(world, cc2singles, cc2pairs, info); + if (world.rank()==0) print_header3("end reiterating CC2 singles"); + } - for (size_t iexcitation = 0; iexcitation < vccs.size(); iexcitation++) { - if (world.rank()==0) print_header1("Solving LRCC2 for excitation " + std::to_string(iexcitation) - + " with omega "+std::to_string(vccs[iexcitation].omega)); - solve_lrcc2(cc2pairs,cc2singles,vccs[iexcitation],iexcitation,info); + + for (size_t iexcitation = 0; iexcitation < vccs.size(); iexcitation++) { + if (world.rank() == 0) + print_header1("Solving LRCC2 for excitation " + std::to_string(iexcitation) + + " with omega " + std::to_string(vccs[iexcitation].omega)); + solve_lrcc2(cc2pairs,cc2singles,vccs[iexcitation],iexcitation,info); } } else MADNESS_EXCEPTION(("Unknown Calculation Type: " + assign_name(ctype)).c_str(), 1); @@ -697,18 +718,14 @@ CC2::solve_cc2(CC_vecfunction& singles, Pairs& doubles, Info& info) cons { output.section("Solving CC2 Ground State"); + if (singles.size()==0) singles=initialize_singles_to_zero(world,CT_CC2, PARTICLE, info); MADNESS_ASSERT(singles.type == PARTICLE); CCTimer time(world, "CC2 Ground State"); + if (world.rank()==0) print("computing CC2 energy"); double omega = CCPotentials::compute_cc2_correlation_energy(world, singles, doubles, info, "initial CC2 energy"); - if (parameters.no_compute_cc2()) { - if (world.rank()==0) print("found no_compute_cc2 key -- recompute singles for the singles-potentials"); - iterate_cc2_singles(world, singles, doubles, info); - return omega; - } - CC_vecfunction ex_singles_dummy; // first singles iteration @@ -716,8 +733,6 @@ CC2::solve_cc2(CC_vecfunction& singles, Pairs& doubles, Info& info) cons // given the doubles, we can solve the singles equations iterate_cc2_singles(world, singles, doubles, info); - // the doubles ansatz depends on the singles and must be updated: |\tau_ij> = |u_ij> + Q12 f12 |t_i t_j> - // update_reg_residues_gs(world, singles, doubles, info); omega = CCPotentials::compute_cc2_correlation_energy(world, singles, doubles, info, "CC2 energy with converged singles"); for (size_t iter = 0; iter < parameters.iter_max(); iter++) { @@ -773,7 +788,7 @@ CC2::solve_cc2(CC_vecfunction& singles, Pairs& doubles, Info& info) cons save(pair.constant_part, pair.name() + "_const"); pair.function().reconstruct(); save(pair.function(), pair.name()); - singles.save_restartdata(world,madness::name(singles.type)); + singles.save_restartdata(world,CC2::singles_name(CT_CC2,singles.type)); } timer1.tag("saving cc2 singles and doubles"); @@ -829,7 +844,7 @@ CC2::solve_lrcc2(Pairs& gs_doubles, const CC_vecfunction& gs_singles, co CCTimer time(world, "Whole LRCC2 Calculation"); /// read LRCC2 singles from file or use the CIS vectors as guess - auto ex_singles=initialize_singles(CT_LRCC2,RESPONSE,false, excitation); + auto ex_singles=initialize_singles(world,CT_LRCC2,RESPONSE, excitation); bool found_singles=(not (ex_singles.get_vecfunction().empty())); if (not found_singles) { if (world.rank()==0) print("using CIS vectors as guess for the LRCC2 singles"); @@ -1025,6 +1040,7 @@ bool CC2::iterate_singles(World& world, CC_vecfunction& singles, const CC_vecfun CCMessenger output(world); if (world.rank()==0) print_header2("Iterating Singles for "+assign_name(ctype)); CCTimer time_all(world, "Overall Iteration of " + assign_name(ctype) + "-Singles"); + if (singles.size()==0) singles=initialize_singles_to_zero(world,ctype, singles.type, info); // consistency checks switch (ctype) { @@ -1227,27 +1243,29 @@ bool CC2::iterate_singles(World& world, CC_vecfunction& singles, const CC_vecfun } CC_vecfunction -CC2::initialize_singles(const CalcType& ctype, const FuncType type, const bool default_to_zero, const int ex) const { +CC2::initialize_singles(World& world, const CalcType& ctype, const FuncType type, const int ex) { MADNESS_CHECK_THROW((type==PARTICLE) or (ex>=0), "Invalid type/excitation combination in initialize_singles"); std::string fname=singles_name(ctype,type,ex); - if (world.rank()==0) print("initializing singles",fname); + if (world.rank()==0) print("initializing singles from file:",fname); CC_vecfunction singles(type); try { singles=CC_vecfunction::load_restartdata(world,fname); if (world.rank()==0) print(" .. singles found on file"); - return singles; } catch (...) { if (world.rank()==0) print(" .. singles not found on file"); } - if (default_to_zero) { - if (world.rank()==0) print(" .. initializing singles to zero functions"); - for (size_t i = parameters.freeze(); i < CCOPS.mo_ket().size(); i++) { - real_function_3d tmpi = real_factory_3d(world); - CCFunction single_i(tmpi, i, type); - singles.insert(i,single_i); - } + return singles; +} + +CC_vecfunction +CC2::initialize_singles_to_zero(World& world, const CalcType& ctype, const FuncType type, const Info& info) { + CC_vecfunction singles(type); + for (size_t i = info.parameters.freeze(); i < info.mo_ket.size(); i++) { + real_function_3d tmpi = real_factory_3d(world); + CCFunction single_i(tmpi, i, type); + singles.insert(i,single_i); } return singles; } @@ -1269,80 +1287,38 @@ CC2::initialize_pairs(Pairs& pairs, const CCState ftype, const CalcType for (size_t j = i; j < CCOPS.mo_ket().size(); j++) { std::string name = CCPair(i, j, ftype, ctype).name(); - if (ftype == GROUND_STATE) { - real_function_6d utmp = real_factory_6d(world); + // if (ftype == GROUND_STATE) { + real_function_6d utmp; const bool found = CCOPS.load_function(utmp, name, info.parameters.debug()); if (found) restarted = true; // if a single pair was found then the calculation is not from scratch real_function_6d const_part; CCOPS.load_function(const_part, name + "_const", info.parameters.debug()); CCPair tmp; - if (ctype==CT_MP2) tmp=CCPotentials::make_pair_mp2(utmp, i, j, info, false); - if (ctype==CT_CC2) tmp=CCPotentials::make_pair_cc2(utmp, tau, i, j, info, false); + if (ctype==CT_MP2) + tmp=CCPotentials::make_pair_mp2(world, utmp, i, j, info, false); + if (ctype==CT_CC2) + tmp=CCPotentials::make_pair_cc2(world, utmp, tau, i, j, info, false); + if (ctype==CT_LRCC2 or ctype==CT_ADC2 or ctype==CT_CISPD) + tmp=CCPotentials::make_pair_lrcc2(world, ctype, utmp, tau, x, i, j, info, false); tmp.constant_part = const_part; pairs.insert(i, j, tmp); - } else if (ftype == EXCITED_STATE) { - // name = std::to_string(int(excitation)) + "_" + name; - real_function_6d utmp = real_factory_6d(world); - const bool found = CCOPS.load_function(utmp, name, info.parameters.debug()); - if (found) restarted = true; - real_function_6d const_part; - CCOPS.load_function(const_part, name + "_const", info.parameters.debug()); - CCPair tmp = CCOPS.make_pair_ex(utmp, tau, x, i, j, ctype); - - { - CCPair tmp2=CCPotentials::make_pair_lrcc2(ctype, utmp, tau, x, i, j, info); - std::swap(tmp,tmp2); - print("going on with Florian's pair"); - // print("going on with Jakob's pair"); - } - - tmp.constant_part = const_part; - pairs.insert(i, j, tmp); - // CCPotentials::compute_excited_pair_energy(world, pairs(i, j), x, info); - } else error("Unknown pairtype"); +// } else if (ftype == EXCITED_STATE) { +// // name = std::to_string(int(excitation)) + "_" + name; +// real_function_6d utmp = real_factory_6d(world); +// const bool found = CCOPS.load_function(utmp, name, info.parameters.debug()); +// if (found) restarted = true; +// real_function_6d const_part; +// CCOPS.load_function(const_part, name + "_const", info.parameters.debug()); +// // CCPair tmp = CCOPS.make_pair_ex(utmp, tau, x, i, j, ctype); +// CCPair tmp=CCPotentials::make_pair_lrcc2(world, ctype, utmp, tau, x, i, j, info); +// print("going on with Florian's pair"); +// tmp.constant_part = const_part; +// pairs.insert(i, j, tmp); +// } else error("Unknown pairtype"); } } return restarted; } -void CC2::update_reg_residues_gs(World& world, const CC_vecfunction& singles, Pairs& doubles, const Info& info) -{ - CCTimer time(world, "Updated Regularization Residues of the Ground State"); - MADNESS_ASSERT(singles.type == PARTICLE); - Pairs updated_pairs; - for (auto& tmp:doubles.allpairs) { - MADNESS_ASSERT(tmp.second.type == GROUND_STATE); - CCPair& pair = tmp.second; - const size_t i = pair.i; - const size_t j = pair.j; - // const CCPair updated_pair = CCOPS.make_pair_gs(pair.function(), singles, i, j); - const CCPair updated_pair = CCPotentials::make_pair_cc2(pair.function(), singles, i, j, info, false); - updated_pairs.insert(i, j, updated_pair); - } - doubles.swap(updated_pairs); - time.info(); -} - -void CC2::update_reg_residues_ex(World& world, const CC_vecfunction& singles, - const CC_vecfunction& response, Pairs& doubles, const Info& info) -{ - CCTimer time(world, "Updated Regularization Residues of the Excited State"); - MADNESS_ASSERT(singles.type == PARTICLE); - MADNESS_ASSERT(response.type == RESPONSE); - CalcType ctype = doubles.allpairs.begin()->second.ctype; - Pairs updated_pairs; - for (auto& tmp:doubles.allpairs) { - MADNESS_ASSERT(tmp.second.type == EXCITED_STATE); - CCPair& pair = tmp.second; - // CCPair updated_pair = CCPotentials::make_pair_ex(pair.function(), singles, response, i, j, pair.ctype); - CCPair updated_pair = - CCPotentials::make_pair_lrcc2(ctype, pair.function(), singles, response, pair.i, pair.j, info); - updated_pairs.insert(pair.i, pair.j, updated_pair); - } - doubles.swap(updated_pairs); - time.info(); -} - - } /* namespace madness */ diff --git a/src/madness/chem/CC2.h b/src/madness/chem/CC2.h index 1deb91862c1..5be4096d0d8 100644 --- a/src/madness/chem/CC2.h +++ b/src/madness/chem/CC2.h @@ -226,19 +226,16 @@ class CC2 : public OptimizationTargetInterface, public QCPropertyInterface { /// type: PARTICLE (cc2) or RESPONSE (lrcc2) /// default_to_zero: if true, initialize with zero functions, otherwise return empty vector /// ex: if type is RESPONSE, the excitation number - CC_vecfunction initialize_singles(const CalcType& ctype, const FuncType type, const bool default_to_zero, int ex=-1) const; + static CC_vecfunction + initialize_singles(World&, const CalcType& ctype, const FuncType type, int ex=-1); + + static CC_vecfunction + initialize_singles_to_zero(World& world, const CalcType& ctype, const FuncType type, const Info& info); /// read pairs from file or initialize new ones bool initialize_pairs(Pairs& pairs, const CCState ftype, const CalcType ctype, const CC_vecfunction& tau, const CC_vecfunction& x, const size_t extitation, const Info& info) const; - static void - update_reg_residues_gs(World& world, const CC_vecfunction& singles, Pairs& doubles, const Info& info); - - static void - update_reg_residues_ex(World& world, const CC_vecfunction& singles, const CC_vecfunction& response, Pairs& doubles, - const Info& info); - /// Iterates a pair of the CC2 doubles equations bool iterate_pair(CCPair& pair, const CC_vecfunction& singles = CC_vecfunction(UNDEFINED)) const; diff --git a/src/madness/chem/CCPotentials.cc b/src/madness/chem/CCPotentials.cc index 9a0d3588f7d..f58cab7eec3 100644 --- a/src/madness/chem/CCPotentials.cc +++ b/src/madness/chem/CCPotentials.cc @@ -63,9 +63,8 @@ CCPotentials::init_orbital_energies(const Nemo& nemo) const { return eps; } -CCPair CCPotentials::make_pair_mp2(const real_function_6d& u, const size_t i, const size_t j, const Info& info, - const bool compute_Q12_f12_ij) { - World& world=u.world(); +CCPair CCPotentials::make_pair_mp2(World& world, const real_function_6d& u, const size_t i, const size_t j, + const Info& info, const bool compute_Q12_f12_ij) { timer t(world, true); // first term is the 6d function u, then follows Q12 f12 |ij> @@ -92,9 +91,8 @@ CCPair CCPotentials::make_pair_mp2(const real_function_6d& u, const size_t i, co return pair; } -CCPair CCPotentials::make_pair_cc2(const real_function_6d& u, const CC_vecfunction& gs_singles, const size_t i, const size_t j, - const Info& info, bool compute_Q12_f12) { - World& world=u.world(); +CCPair CCPotentials::make_pair_cc2(World& world, const real_function_6d& u, const CC_vecfunction& gs_singles, const size_t i, + const size_t j, const Info& info, bool compute_Q12_f12) { // first term is the 6d function u, then follows Q12 f12 |t_i t_j> std::vector> functions; @@ -122,45 +120,48 @@ CCPair CCPotentials::make_pair_cc2(const real_function_6d& u, const CC_vecfuncti } /// follow eq. (23) of Kottmann, JCTC 13, 5956 (2017) -CCPair CCPotentials::make_pair_lrcc2(const CalcType& ctype, const real_function_6d& u, - const CC_vecfunction& gs_singles, const CC_vecfunction& ex_singles, const size_t i, const size_t j, const Info& info) { +CCPair CCPotentials::make_pair_lrcc2(World& world, const CalcType& ctype, const real_function_6d& u, + const CC_vecfunction& gs_singles, const CC_vecfunction& ex_singles, const size_t i, const size_t j, const Info& info, + const bool compute_Q12_f12) { MADNESS_ASSERT(gs_singles.type == PARTICLE || gs_singles.type == HOLE); MADNESS_ASSERT(ex_singles.type == RESPONSE); MADNESS_ASSERT(ctype == CT_CISPD || ctype == CT_LRCC2 || ctype == CT_ADC2); MADNESS_ASSERT(!(i < info.parameters.freeze())); MADNESS_ASSERT(!(j < info.parameters.freeze())); - World& world=u.world(); - // compute the t intermediates for active orbitals only -- they go into the ansatz - const auto t = CC_vecfunction(info.get_active_mo_ket()+gs_singles.get_vecfunction(),MIXED,info.parameters.freeze()); - MADNESS_ASSERT(t.size() == (info.mo_ket.size()-info.parameters.freeze())); + typedef CCPairFunction cpT; + auto functions=std::vector(1,cpT(u)); - // compute the t intermediates for all orbitals -- they go into the projector - const CC_vecfunction pt = copy(make_full_t_intermediate(gs_singles,info)); - MADNESS_ASSERT(pt.size() == info.mo_ket.size()); + if (compute_Q12_f12) { + // compute the t intermediates for active orbitals only -- they go into the ansatz + const auto t = CC_vecfunction(info.get_active_mo_ket()+gs_singles.get_vecfunction(),MIXED,info.parameters.freeze()); + MADNESS_ASSERT(t.size() == (info.mo_ket.size()-info.parameters.freeze())); - auto f12=CCConvolutionOperatorPtr(world,OT_F12,info.parameters); + // compute the t intermediates for all orbitals -- they go into the projector + const CC_vecfunction pt = copy(make_full_t_intermediate(gs_singles,info)); + MADNESS_ASSERT(pt.size() == info.mo_ket.size()); - // set up projectors -- they project out the occupied space from the response pair function + auto f12=CCConvolutionOperatorPtr(world,OT_F12,info.parameters); - // dQ12t = -(Qt(1) Ox(2) + Ox(1) Qt(2)) eq. (22) of the excited state paper - QProjector Qt(info.mo_bra,pt.get_vecfunction()); - Projector Ox(info.get_active_mo_bra(),ex_singles.get_vecfunction()); // this works on active orbitals only - auto dQt_1 = outer(Qt,Ox); - auto dQt_2 = outer(Ox,Qt); + // set up projectors -- they project out the occupied space from the response pair function - StrongOrthogonalityProjector Q12t(world); // eq. (21) of the ground state paper - Q12t.set_spaces(info.mo_bra,pt.get_vecfunction(),info.mo_bra,pt.get_vecfunction()); + // dQ12t = -(Qt(1) Ox(2) + Ox(1) Qt(2)) eq. (22) of the excited state paper + QProjector Qt(info.mo_bra,pt.get_vecfunction()); + Projector Ox(info.get_active_mo_bra(),ex_singles.get_vecfunction()); // this works on active orbitals only + auto dQt_1 = outer(Qt,Ox); + auto dQt_2 = outer(Ox,Qt); - typedef CCPairFunction cpT; - auto functions=std::vector(1,cpT(u)); + StrongOrthogonalityProjector Q12t(world); // eq. (21) of the ground state paper + Q12t.set_spaces(info.mo_bra,pt.get_vecfunction(),info.mo_bra,pt.get_vecfunction()); - auto f_xt=std::vector(1,cpT(f12, ex_singles(i), t(j))); - auto f_tx=std::vector(1,cpT(f12, t(i), ex_singles(j))); - auto f_tt=std::vector(1,cpT(f12, t(i), t(j))); - functions+=(Q12t(f_xt) + Q12t(f_tx) - dQt_1(f_tt) -dQt_2(f_tt)); // note the sign change in the last two terms - functions=consolidate(functions); + auto f_xt=std::vector(1,cpT(f12, ex_singles(i), t(j))); + auto f_tx=std::vector(1,cpT(f12, t(i), ex_singles(j))); + auto f_tt=std::vector(1,cpT(f12, t(i), t(j))); + + functions+=(Q12t(f_xt) + Q12t(f_tx) - dQt_1(f_tt) -dQt_2(f_tt)); // note the sign change in the last two terms + functions=consolidate(functions); + } CCPair pair(i, j, EXCITED_STATE, ctype, functions); MADNESS_ASSERT(ex_singles.omega != 0.0); @@ -1273,7 +1274,7 @@ CCPotentials::make_constant_part_cispd_Qt(const CCPair& u, const CC_vecfunction& CCTimer time_cpr(world, "Commutator-Projector Response"); { time_CPR.start(); - const vector_real_function_3d Vxtmp = sub(world, get_potentials(x, POT_singles_), + const vector_real_function_3d Vxtmp = sub(world, get_potentials(x, POT_singles_,true), x.omega * x.get_vecfunction()); const CC_vecfunction Vx(Vxtmp, UNDEFINED, parameters.freeze()); CCPairFunction ftt(f12, moi, moj); @@ -1379,8 +1380,8 @@ std::vector> // calculate the regularized potential real_function_6d V=real_factory_6d(world); std::vector> V_lowrank; - if (exists("Ue")) V += apply_Ue(world,ti,tj,info,&Gscreen); - if (exists("KffK")) V -= apply_KffK(world,ti,tj,info,&Gscreen); + // if (exists("Ue")) V += apply_Ue(world,ti,tj,info,&Gscreen); + // if (exists("KffK")) V -= apply_KffK(world,ti,tj,info,&Gscreen); if (exists("reduced_Fock")) V += apply_reduced_F(world,ti,tj,info,&Gscreen); if (exists("comm_F_Qt_f12")) { V_lowrank += apply_commutator_F_Qt_f12(world,ti,tj,gs_singles,ex_singles,info,&Gscreen); @@ -1465,8 +1466,8 @@ CCPotentials::apply_reduced_F1(const CCFunction& ti, const CCFunction< //CC_Timer time(world,"(F-eij)|"+ti.name()+tj.name()+">"); // get singles potential const bool symmetric = (ti.type == tj.type && ti.i == tj.i); - const real_function_3d Vti = get_potentials(ti, POT_singles_); - const real_function_3d Vtj = get_potentials(tj, POT_singles_); + const real_function_3d Vti = get_potentials(ti, POT_singles_,true); + const real_function_3d Vtj = get_potentials(tj, POT_singles_,true); const real_function_6d Vt = make_f_xy(Vti, tj, Gscreen); real_function_6d tV; if (symmetric) tV = swap_particles(Vt); @@ -1488,8 +1489,8 @@ CCPotentials::apply_reduced_F(World& world, const CCFunction& ti, cons //CC_Timer time(world,"(F-eij)|"+ti.name()+tj.name()+">"); // get singles potential const bool symmetric = (ti == tj); - const real_function_3d Vti = info.intermediate_potentials(ti, POT_singles_); - const real_function_3d Vtj = info.intermediate_potentials(tj, POT_singles_); + const real_function_3d Vti = info.intermediate_potentials(ti, POT_singles_,true); + const real_function_3d Vtj = info.intermediate_potentials(tj, POT_singles_,true); const real_function_6d Vt = make_f_xy(world, Vti, tj, info, Gscreen); real_function_6d tV; if (symmetric) tV = madness::swap_particles(Vt); @@ -1791,7 +1792,7 @@ CCPotentials::apply_commutator_F_Qt_f12(World& world, const CCFunction auto f12=CCConvolutionOperatorPtr(world,OT_F12,parameters); auto ftt=std::vector>({CCPairFunction(f12, phi_i.function, phi_j.function)}); - const vector_real_function_3d Vtau=info.intermediate_potentials(gs_singles, POT_singles_); + const vector_real_function_3d Vtau=info.intermediate_potentials(gs_singles, POT_singles_,true); Projector OVtau(info.get_active_mo_bra(),Vtau); QProjector Qt(info.get_active_mo_bra(),gs_singles.get_vecfunction()); @@ -1821,8 +1822,8 @@ CCPotentials::apply_commutator_F_dQt_f12(World& world, const CCFunction>({CCPairFunction(f12, phi_i.function, phi_j.function)}); auto t=CCPotentials::make_active_t_intermediate(gs_singles,info); - const vector_real_function_3d Vtau=info.intermediate_potentials(gs_singles, POT_singles_); - const vector_real_function_3d Vx=info.intermediate_potentials(ex_singles, POT_singles_); + const vector_real_function_3d Vtau=info.intermediate_potentials(gs_singles, POT_singles_,true); + const vector_real_function_3d Vx=info.intermediate_potentials(ex_singles, POT_singles_,true); auto bra=info.get_active_mo_bra(); Projector OVtau(bra,Vtau); @@ -2292,6 +2293,8 @@ CCPotentials::get_CC2_singles_potential_gs(World& world, const CC_vecfunction& s // set up taskq MacroTaskSinglesPotentialGs taskgs; + taskgs.partitioner->min_batch_size=2; + taskgs.partitioner->max_batch_size=2; int printlevel=(info.parameters.debug()) ? 3 : 0; auto taskq=std::shared_ptr(new MacroTaskQ(world,world.size(),printlevel)); MacroTask mtaskgs(world,taskgs,taskq); @@ -2365,6 +2368,7 @@ CCPotentials::get_CCS_potential_ex(World& world, const CC_vecfunction& x, const auto taskq=std::shared_ptr(new MacroTaskQ(world,world.size(),printlevel)); MacroTask mtask(world,task,taskq); std::vector result_index(x.size()); + for (int i=0; i& result_ if (name == POT_F3D_) { result = fock_residue_closed_shell(world, singles_external, info); - } else if (name == POT_ccs_) { + } else if (name == POT_ccs_) { // or S3c + // this term corresponds to the diagrams S3c, S5b, S5c and S6, see eq. (28) of the GS paper const CC_vecfunction t = make_active_t_intermediate(singles,info); QProjector Qt(info.get_active_mo_bra(),t.get_vecfunction()); result = Qt(ccs_unprojected(world, t_external, singles, info)); @@ -2654,24 +2659,37 @@ CCPotentials::potential_singles_ex(World& world, const std::vector result_i MADNESS_ASSERT(singles_gs.type == PARTICLE); MADNESS_ASSERT(singles_ex.type == RESPONSE); + vector_real_function_3d result, intermediate; + + // collect all singles that correspond to external lines + auto get_batch_of_external_singles = [](const CC_vecfunction& singles, const std::vector& result_index) { + std::vector> singles_tmp; + for (int i : result_index) singles_tmp.push_back(singles(i)); + return CC_vecfunction(singles_tmp); + }; + + const CC_vecfunction singles_ex_external=get_batch_of_external_singles(singles_ex,result_index); + Projector Ox(info.get_active_mo_bra(),singles_ex.get_vecfunction()); - vector_real_function_3d result, intermediate; CCTimer timer(world, "timer-ex-potential"); if (name == POT_F3D_) { - result = fock_residue_closed_shell(world, singles_ex, info); + result = fock_residue_closed_shell(world, singles_ex_external, info); } else if (name == POT_ccs_) { - // const CC_vecfunction t = make_t_intermediate(singles_gs,info.parameters); + const CC_vecfunction singles_gs_external=get_batch_of_external_singles(singles_gs,result_index); + const CC_vecfunction t_gs_external = make_active_t_intermediate(singles_gs_external,info); + // this is the functional derivative of the GS terms S3c + S5b + S5c + S6 = Qt(S3c^(t_i)) + // see Tab. 3 line 2 and Eq. (19) of the excited state paper const CC_vecfunction t = make_active_t_intermediate(singles_gs,info); QProjector Qt(info.get_active_mo_bra(),t.get_vecfunction()); - // vector_real_function_3d part1 = apply_Qt(ccs_unprojected(world, t, singles_ex, info), t); - // vector_real_function_3d part2 = apply_Qt(ccs_unprojected(world, singles_ex, singles_gs, info), t); - vector_real_function_3d part1 = Qt(ccs_unprojected(world, t, singles_ex, info)); - vector_real_function_3d part2 = Qt(ccs_unprojected(world, singles_ex, singles_gs, info)); - // vector_real_function_3d part3 = apply_projector(ccs_unprojected(world, t, singles_gs, info), singles_ex); - vector_real_function_3d part3 = Ox(ccs_unprojected(world, t, singles_gs, info)); - vector_real_function_3d tmp = add(world, part1, part2); - result = sub(world, tmp, part3); + // part1 + part2 correspond to eq. (19) third term: Qt(dS5b(t_i)) (missing superscript t_i in the equation seems a typo) + // ccs_unprojected corresponds to Tab 3 line 2 first diagram + vector_real_function_3d part1 = Qt(ccs_unprojected(world, t_gs_external, singles_ex, info)); + // ccs_unprojected corresponds to Tab 3 line 2 second diagram + vector_real_function_3d part2 = Qt(ccs_unprojected(world, singles_ex_external, singles_gs, info)); + // corresponds to eq. (19) second but last term: dQt S5b(t_i) = -Ox(S5b(ti)) + vector_real_function_3d part3 = Ox(ccs_unprojected(world, t_gs_external, singles_gs, info)); + result=part1 + part2 - part3; } else if (name == POT_s2b_) { std::tie(result, intermediate) = s2b(world, result_index, singles_ex, doubles_ex, info); } else if (name == POT_s2c_) { @@ -2687,10 +2705,11 @@ CCPotentials::potential_singles_ex(World& world, const std::vector result_i vector_real_function_3d s4c_part2 = s4c(world, result_index, singles_ex, doubles_gs, info); result = add(world, s4c_part1, s4c_part2); } else if (name == POT_cis_) { - result = ccs_unprojected(world, CC_vecfunction(info.get_active_mo_ket(), HOLE, info.parameters.freeze()), singles_ex, info); + CC_vecfunction orbitals(info.mo_ket,HOLE); + const CC_vecfunction orbitals_external=get_batch_of_external_singles(orbitals,result_index); + result = ccs_unprojected(world, orbitals_external, singles_ex, info); } else MADNESS_EXCEPTION(("potential_singles: Unknown potential " + assign_name(name)).c_str(), 1) - ; const double size = get_size(world, result); const double norm = norm2(world, result); const std::pair time = timer.current_time(); @@ -3171,7 +3190,7 @@ CCPotentials::s2b(World& world, const std::vector external_indices, const C // which "external_indices" refers to a subset of the full external indices (active orbitals) only (through Macrotasks), // the intermediates result_u is stored as the full vector [nfreeze,nocc) - vector_real_function_3d result_u = info.intermediate_potentials(singles, POT_s2b_); + vector_real_function_3d result_u = info.intermediate_potentials(singles, POT_s2b_,false); vector_real_function_3d result_u_return; bool recalc_u_part = false; if (result_u.empty()) recalc_u_part = true; @@ -3210,7 +3229,7 @@ std::tuple CCPotentials::s2c(World& world, const std::vector external_index, const CC_vecfunction& singles, const Pairs& doubles, const Info& info) { vector_real_function_3d result; // see if we can skip the recalculation of the pure 6D part since this does not change during the singles iteration - vector_real_function_3d result_u = info.intermediate_potentials(singles, POT_s2c_); + vector_real_function_3d result_u = info.intermediate_potentials(singles, POT_s2c_,false); vector_real_function_3d result_u_return; bool recalc_u_part = false; if (result_u.empty()) recalc_u_part = true; diff --git a/src/madness/chem/CCPotentials.h b/src/madness/chem/CCPotentials.h index 01c6efbfd5f..b21060709f9 100644 --- a/src/madness/chem/CCPotentials.h +++ b/src/madness/chem/CCPotentials.h @@ -211,17 +211,18 @@ class CCPotentials { public: /// return the regularized MP2 ansatz: |\tau_ij> = |u_ij> + Q12 f12 |ij> - static CCPair make_pair_mp2(const real_function_6d& u, const size_t i, const size_t j, const Info& info, - bool compute_Q12_f12_ij); + static CCPair make_pair_mp2(World& world, const real_function_6d& u, const size_t i, const size_t j, + const Info& info, bool compute_Q12_f12_ij); /// return the regularized CC2 ansatz: |\tau_ij> = |u_ij> + Q12t f12 |t_i t_j> - static CCPair make_pair_cc2(const real_function_6d& u, const CC_vecfunction& gs_singles, - const size_t i, const size_t j, const Info& info, const bool compute_Q12_f12_ij); + static CCPair make_pair_cc2(World& world, const real_function_6d& u, + const CC_vecfunction& gs_singles, const size_t i, const size_t j, const Info& info, const bool compute_Q12_f12_ij); /// return the regularized CC2 ansatz: |x_ij> = |u_ij> + Q12t f12 |t_i t_j> + ????? - static CCPair make_pair_lrcc2(const CalcType& ctype, const real_function_6d& u, - const CC_vecfunction& gs_singles, const CC_vecfunction& ex_singles, - const size_t i, const size_t j, const Info& info); + static CCPair make_pair_lrcc2(World& world, const CalcType& ctype, + const real_function_6d& u, const CC_vecfunction& gs_singles, + const CC_vecfunction& ex_singles, const size_t i, const size_t j, const Info& info, + const bool compute_Q12_f12); // Pair functions diff --git a/src/madness/chem/CCStructures.cc b/src/madness/chem/CCStructures.cc index f8a9f9cc178..4512ff55ef3 100644 --- a/src/madness/chem/CCStructures.cc +++ b/src/madness/chem/CCStructures.cc @@ -94,50 +94,46 @@ CCPair::info() const { } madness::vector_real_function_3d -CCIntermediatePotentials::operator()(const CC_vecfunction& f, const PotentialType& type) const { - output("Getting " + assign_name(type) + " for " + f.name(0)); +CCIntermediatePotentials::get_potential(const PotentialType& ptype, const FuncType& ftype, const bool throw_if_empty) const { vector_real_function_3d result; - if (type == POT_singles_ and (f.type == PARTICLE or f.type == MIXED)) result= current_singles_potential_gs_; - else if (type == POT_singles_ and f.type == RESPONSE) result= current_singles_potential_ex_; - else if (type == POT_s2b_ and f.type == PARTICLE) result= current_s2b_potential_gs_; - else if (type == POT_s2b_ and f.type == RESPONSE) result= current_s2b_potential_ex_; - else if (type == POT_s2c_ and f.type == PARTICLE) result= current_s2c_potential_gs_; - else if (type == POT_s2c_ and f.type == RESPONSE) result= current_s2c_potential_ex_; - else if (f.type == HOLE) { - output(assign_name(type) + " is zero for HOLE states"); - // result = zero_functions(f.size()); - } else { - output("ERROR: Potential was not supposed to be stored"); - MADNESS_EXCEPTION("Potential was not supposed to be stored", 1); + if (parameters.debug()) print("Getting " , assign_name(ptype) , " for " , madness::name(ftype,0)); + if (ptype == POT_singles_ and (ftype == PARTICLE or ftype == MIXED)) result= current_singles_potential_gs_; + else if (ptype == POT_singles_ and ftype == RESPONSE) result= current_singles_potential_ex_; + else if (ptype == POT_s2b_ and ftype == PARTICLE) result= current_s2b_potential_gs_; + else if (ptype == POT_s2b_ and ftype == RESPONSE) result= current_s2b_potential_ex_; + else if (ptype == POT_s2c_ and ftype == PARTICLE) result= current_s2c_potential_gs_; + else if (ptype == POT_s2c_ and ftype == RESPONSE) result= current_s2c_potential_ex_; + else { + MADNESS_EXCEPTION("unknown potential in CCIntermediatePotentials::get_potential", 1); + } + if (result.empty() and throw_if_empty) { + std::string errmsg="CCIntermediatePotential was not computed/stored "+assign_name(ptype) + " " +assign_name(ftype); + errmsg+="\n --> you might need to iterate the corresponding singles"; + print(errmsg); + MADNESS_EXCEPTION(errmsg.c_str(),1); } - if (result.empty()) { - output("!!!WARNING: Potential is empty!!!"); - } else { + if (not result.empty() and parameters.debug()) { World& world=result.front().world(); if (parameters.debug()) print_size(world,result, "potential"); } - return result; } + +madness::vector_real_function_3d +CCIntermediatePotentials::operator()(const CC_vecfunction& f, const PotentialType& type, + const bool throw_if_empty) const { + return get_potential(type,f.type,throw_if_empty); +} + madness::real_function_3d -CCIntermediatePotentials::operator()(const CCFunction& f, const PotentialType& type) const { - output("Getting " + assign_name(type) + " for " + f.name()); - std::vector result; - if (type == POT_singles_ and (f.type == PARTICLE or f.type == MIXED)) result= current_singles_potential_gs_; - else if (type == POT_singles_ and f.type == RESPONSE) result= current_singles_potential_ex_; - else if (type == POT_s2b_ and f.type == PARTICLE) result= current_s2b_potential_gs_; - else if (type == POT_s2b_ and f.type == RESPONSE) result= current_s2b_potential_ex_; - else if (type == POT_s2c_ and f.type == PARTICLE) result= current_s2c_potential_gs_; - else if (type == POT_s2c_ and f.type == RESPONSE) result= current_s2c_potential_ex_; - else if (f.type == HOLE) output(assign_name(type) + " is zero for HOLE states"); - else MADNESS_EXCEPTION("Potential was not supposed to be stored", 1); - - std::string errmsg="CCIntermediatePotential was not computed/stored "+assign_name(type) + " " +assign_name(f.type); - errmsg+="\n --> you might need to iterate the corresponding singles"; - MADNESS_CHECK_THROW(result.size()>(f.i-parameters.freeze()),errmsg.c_str()); - return result[f.i-parameters.freeze()]; +CCIntermediatePotentials::operator()(const CCFunction& f, const PotentialType& type, + const bool throw_if_empty) const { + vector_real_function_3d result=get_potential(type,f.type,throw_if_empty); + long iact=f.i-parameters.freeze(); // active index + MADNESS_CHECK_THROW(iact& result_index, auto& tau=x.second; MADNESS_CHECK_THROW(tau.functions.size()==1,"doubles in MacroTaskSinglesPotentialsEx should only contain one function"); bool compute_Q12_F12=(PotentialType(name)==POT_s2b_ or PotentialType(name)==POT_s2c_); - tau=CCPotentials::make_pair_cc2(tau.function(),singles_gs,tau.i,tau.j,info, compute_Q12_F12); + x.second=CCPotentials::make_pair_cc2(world,tau.function(),singles_gs,tau.i,tau.j, info, compute_Q12_F12); } // the doubles currently only contain the full 6d function -> complete it with the Q12 f12 |ti tj> part for (auto& x : doubles_ex1.allpairs) { auto& tau=x.second; MADNESS_CHECK_THROW(tau.functions.size()==1,"doubles in MacroTaskSinglesPotentialsEx should only contain one function"); - tau=CCPotentials::make_pair_lrcc2(tau.ctype,tau.function(),singles_gs,singles_ex,tau.i,tau.j,info); + x.second=tau=CCPotentials::make_pair_lrcc2(world,tau.ctype,tau.function(),singles_gs,singles_ex,tau.i,tau.j, info, true); } resultT result=CCPotentials::potential_singles_ex(world, @@ -636,7 +632,7 @@ MacroTaskSinglesPotentialGs::operator()(const std::vector& result_index, auto& tau=x.second; MADNESS_CHECK_THROW(tau.functions.size()==1,"doubles in MacroTaskSinglesPotentialsGS should only contain one function"); bool compute_Q12_F12=(PotentialType(name)==POT_s2b_ or PotentialType(name)==POT_s2c_); - tau=CCPotentials::make_pair_cc2(tau.function(),singles_gs,tau.i,tau.j,info,compute_Q12_F12); + tau=CCPotentials::make_pair_cc2(world,tau.function(),singles_gs,tau.i,tau.j,info, compute_Q12_F12); } resultT result=CCPotentials::potential_singles_gs(world, result_index, @@ -663,10 +659,10 @@ MacroTaskComputeCorrelationEnergy::operator()(const std::vector& pairs, for (int i=0; i part is not stored in the cloud, so recompute it here - auto pair=CCPotentials::make_pair_mp2(pairs[i].function(),pairs[i].i,pairs[i].j,info,true); + auto pair=CCPotentials::make_pair_mp2(world,pairs[i].function(),pairs[i].i,pairs[i].j,info, true); result[i]=CCPotentials::compute_pair_correlation_energy(world,pair,singles_gs,info); } else if (ctype==CT_CC2) { - auto pair=CCPotentials::make_pair_cc2(pairs[i].function(),singles_gs,pairs[i].i,pairs[i].j,info,true); + auto pair=CCPotentials::make_pair_cc2(world,pairs[i].function(),singles_gs,pairs[i].i,pairs[i].j,info, true); result[i]=CCPotentials::compute_pair_correlation_energy(world,pair,singles_gs,info); } else { MADNESS_EXCEPTION("MacroTaskComputeCorrelationEnergy: unknown ctype",1); diff --git a/src/madness/chem/CCStructures.h b/src/madness/chem/CCStructures.h index c7e90fae492..266a758d20e 100644 --- a/src/madness/chem/CCStructures.h +++ b/src/madness/chem/CCStructures.h @@ -1107,17 +1107,34 @@ class CCPair : public archive::ParallelSerializableObject { struct CCIntermediatePotentials { CCIntermediatePotentials() = default; CCIntermediatePotentials(const CCParameters& p) : parameters(p) {}; - CCIntermediatePotentials(const CCIntermediatePotentials& other) = default; CCIntermediatePotentials& operator=(const CCIntermediatePotentials& other) = default; + /// check if the intermediate potential exists + bool potential_exists(const CC_vecfunction& f, const PotentialType& type) const { + return potential_exists(type,f.type); + } + + /// check if the intermediate potential exists + bool potential_exists(const PotentialType& type,const FuncType& ftype) const { + bool exists=get_potential(type,ftype,false).size()>0; + return exists; + } + + /// return a vector of the intermediate potentials + + /// @param[in] ptype: the potential type (POT_SINGLES, POT_S2B, ..) + /// @param[in] ftype: the function type (HOLE, PARTICLE, RESPONSE) + vector_real_function_3d + get_potential(const PotentialType& ptype, const FuncType& ftype, const bool throw_if_empty) const; + /// fetches the correct stored potential or throws an exception vector_real_function_3d - operator()(const CC_vecfunction& f, const PotentialType& type) const; + operator()(const CC_vecfunction& f, const PotentialType& type, const bool throw_if_empty) const; /// fetch the potential for a single function Function - operator()(const CCFunction& f, const PotentialType& type) const; + operator()(const CCFunction& f, const PotentialType& type, const bool throw_if_empty) const; /// deltes all stored potentials void clear_all() { @@ -1474,6 +1491,8 @@ class MacroTaskSinglesPotentialEx : public MacroTaskOperationBase { std::string basename="SinglesPotentialEx"; MacroTaskSinglesPotentialEx() { name="SinglesPotentialEx"; + partitioner->max_batch_size=2; + partitioner->min_batch_size=2; } typedef std::tuple<