From 00ea5713d4a419106b784d88151c6a590b912a7c Mon Sep 17 00:00:00 2001 From: fbischoff Date: Mon, 25 Nov 2024 16:51:29 +0100 Subject: [PATCH] LRCC2 works for BH --- src/madness/chem/CC2.cc | 3 +- src/madness/chem/CCPotentials.cc | 49 +++++++++++++++++++++++--------- src/madness/chem/CCPotentials.h | 6 ++-- src/madness/chem/CCStructures.h | 2 +- 4 files changed, 42 insertions(+), 18 deletions(-) diff --git a/src/madness/chem/CC2.cc b/src/madness/chem/CC2.cc index 997e671d92d..c45875042a5 100644 --- a/src/madness/chem/CC2.cc +++ b/src/madness/chem/CC2.cc @@ -676,7 +676,8 @@ CC2::iterate_lrcc2_pairs(World& world, const CC_vecfunction& cc2_s, } // if no function has been computed so far use the constant part (first iteration) - for (auto& pair : pair_vec) if (not pair.function_exists()) pair.update_u(pair.constant_part); + // for (auto& pair : pair_vec) if (not pair.function_exists()) pair.update_u(pair.constant_part); + for (auto& pair : pair_vec) if (not pair.function_exists()) pair.update_u(real_factory_6d(world)); for (const auto& p : pair_vec) p.constant_part.print_size("constant_part before iter"); for (const auto& p : pair_vec) p.function().print_size("u before iter"); diff --git a/src/madness/chem/CCPotentials.cc b/src/madness/chem/CCPotentials.cc index 626aa08d9a4..4524e9bba34 100644 --- a/src/madness/chem/CCPotentials.cc +++ b/src/madness/chem/CCPotentials.cc @@ -2451,6 +2451,10 @@ CCPotentials::get_CC2_singles_potential_ex(World& world, const CC_vecfunction& g // vector_real_function_3d potential = apply_Qt(unprojected, mo_ket_); vector_real_function_3d potential = Q(unprojected); if (info.parameters.debug()) { + CCTimer timer(world, "potential energy of F3D"); + + print("fock residue"); + timer.info(true, norm2(world, fock_residue)); if (world.rank()==0) print_header3("intermediate potentials in CC2 singles potential excited state"); madness::print_size(world,fock_residue,"fock_residue"); @@ -2562,7 +2566,8 @@ CCPotentials::potential_singles_gs(World& world, const std::vector& result_ const CC_vecfunction singles_external(singles_tmp); const CC_vecfunction t_external = make_active_t_intermediate(singles_external,info); - auto builder=CCPairBuilder(world,info).set_gs_singles(singles); + auto ctype=doubles.allpairs.begin()->second.ctype; + auto builder=CCPairBuilder(world,info).set_gs_singles(singles).set_ctype(ctype); if (name == POT_F3D_) { result = fock_residue_closed_shell(world, singles_external, info); @@ -2593,9 +2598,9 @@ CCPotentials::potential_singles_gs(World& world, const std::vector& result_ } else if (name == POT_s4a_) { error("potential_singles: Demanded s4a potential -> this is calculated along with the s2b potential"); } else if (name == POT_s4b_) { - result = s4b(world, result_index, singles, doubles, info); + result = s4b(world, result_index, singles, doubles, builder, info); } else if (name == POT_s4c_) { - result = s4c(world, result_index, singles, doubles, info); + result = s4c(world, result_index, singles, doubles, builder, info); } else MADNESS_EXCEPTION(("potential_singles: Unknown potential " + assign_name(name)).c_str(), 1) madness::print_size(world,result,"result of "+assign_name(name)); @@ -2675,7 +2680,12 @@ CCPotentials::potential_singles_ex(World& world, const std::vector result_i MADNESS_ASSERT(singles_gs.type == PARTICLE); MADNESS_ASSERT(singles_ex.type == RESPONSE); - auto builder=CCPairBuilder(world,info).set_gs_singles(singles_gs).set_ex_singles(singles_ex); + // gs depends on es: CC2 -> LRCC2; MP2 -> ADC(2), etc + auto ctype_ex=doubles_ex.allpairs.begin()->second.ctype; + auto builder_ex=CCPairBuilder(world,info).set_gs_singles(singles_gs).set_ex_singles(singles_ex).set_ctype(ctype_ex); + + auto ctype_gs=doubles_gs.allpairs.begin()->second.ctype; + auto builder_gs=CCPairBuilder(world,info).set_gs_singles(singles_gs).set_ex_singles(singles_ex).set_ctype(ctype_gs); vector_real_function_3d result, intermediate; @@ -2709,18 +2719,18 @@ CCPotentials::potential_singles_ex(World& world, const std::vector result_i 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, builder, info); + std::tie(result, intermediate) = s2b(world, result_index, singles_ex, doubles_ex, builder_ex, info); } else if (name == POT_s2c_) { - std::tie(result, intermediate) = s2c(world, result_index, singles_ex, doubles_ex, builder, info); + std::tie(result, intermediate) = s2c(world, result_index, singles_ex, doubles_ex, builder_ex, info); } else if (name == POT_s4a_) { error("potential_singles: Demanded s4a potential -> this is calculated from the s2b potential"); } else if (name == POT_s4b_) { - vector_real_function_3d s4b_part1 = s4b(world, result_index, singles_gs, doubles_ex, info); - vector_real_function_3d s4b_part2 = s4b(world, result_index, singles_ex, doubles_gs, info); + vector_real_function_3d s4b_part1 = s4b(world, result_index, singles_gs, doubles_ex, builder_ex, info); + vector_real_function_3d s4b_part2 = s4b(world, result_index, singles_ex, doubles_gs, builder_gs, info); result = add(world, s4b_part1, s4b_part2); } else if (name == POT_s4c_) { - vector_real_function_3d s4c_part1 = s4c(world, result_index, singles_gs, doubles_ex, info); - vector_real_function_3d s4c_part2 = s4c(world, result_index, singles_ex, doubles_gs, info); + vector_real_function_3d s4c_part1 = s4c(world, result_index, singles_gs, doubles_ex, builder_ex, info); + vector_real_function_3d s4c_part2 = s4c(world, result_index, singles_ex, doubles_gs, builder_gs, info); result = add(world, s4c_part1, s4c_part2); } else if (name == POT_cis_) { CC_vecfunction orbitals(info.mo_ket,HOLE); @@ -3341,7 +3351,7 @@ CCPotentials::s4a_from_s2b(const vector_real_function_3d& s2b, const CC_vecfunct } madness::vector_real_function_3d -CCPotentials::s4b(World& world, const std::vector external_index, const CC_vecfunction& singles, const Pairs& doubles, const Info& info) +CCPotentials::s4b(World& world, const std::vector external_index, const CC_vecfunction& singles, const Pairs& doubles, const CCPairBuilder& builder, const Info& info) { auto g12=CCConvolutionOperator(world,OT_G12,info.parameters); vector_real_function_3d result; @@ -3357,7 +3367,12 @@ CCPotentials::s4b(World& world, const std::vector external_index, const CC_ truncate(world, l_kgi); for (const auto& ltmp : singles.functions) { const size_t l = ltmp.first; - const std::vector> ukl = get_pair_function(doubles, k, l); + // const std::vector> ukl = get_pair_function(doubles, k, l); + // make a ful pair function using the pure 6D part uik and complete with the rest + std::vector> ukl1 = get_pair_function(doubles, k, l); + CCPair pair_kl=builder.complete_pair_with_low_rank_parts(builder.make_pair(i,l,ukl1)); + std::vector> ukl = pair_kl.functions; + for (size_t mm = 0; mm < ukl.size(); mm++) { resulti += -2.0 * ukl[mm].project_out(l_kgi[l - info.parameters.freeze()], 2); resulti += ukl[mm].project_out(l_kgi[l - info.parameters.freeze()], 1); @@ -3370,7 +3385,7 @@ CCPotentials::s4b(World& world, const std::vector external_index, const CC_ } madness::vector_real_function_3d -CCPotentials::s4c(World& world, const std::vector external_index, const CC_vecfunction& singles, const Pairs& doubles, const Info& info) +CCPotentials::s4c(World& world, const std::vector external_index, const CC_vecfunction& singles, const Pairs& doubles, const CCPairBuilder& builder, const Info& info) { vector_real_function_3d result; auto g12=CCConvolutionOperator(world,OT_G12,info.parameters); @@ -3392,7 +3407,13 @@ CCPotentials::s4c(World& world, const std::vector external_index, const CC_ truncate(world, l_kgtauk); for (const auto& ltmp : singles.functions) { const size_t l = ltmp.first; - const std::vector> uil = get_pair_function(doubles, i, l); +// const std::vector> uil = get_pair_function(doubles, i, l); + + // make a ful pair function using the pure 6D part uik and complete with the rest + std::vector> uil1 = get_pair_function(doubles, i, l); + CCPair pair_il=builder.complete_pair_with_low_rank_parts(builder.make_pair(i,l,uil1)); + std::vector> uil = pair_il.functions; + for (size_t mm = 0; mm < uil.size(); mm++) { part1 += uil[mm].project_out(l_kgtauk[l - info.parameters.freeze()], 2); part2 += uil[mm].project_out(l_kgtauk[l - info.parameters.freeze()], 1); diff --git a/src/madness/chem/CCPotentials.h b/src/madness/chem/CCPotentials.h index f4161bc4ec3..28d7d947e53 100644 --- a/src/madness/chem/CCPotentials.h +++ b/src/madness/chem/CCPotentials.h @@ -932,7 +932,8 @@ class CCPotentials { ///@param[out] \f$ -( _2 - _1) | kgtaui = | taui=singles_i \f$ /// Q-Projector is not applied, sign is correct static vector_real_function_3d - s4b(World& world, std::vector external_index, const CC_vecfunction& singles, const Pairs& doubles, const Info& info); + s4b(World& world, std::vector external_index, const CC_vecfunction& singles, const Pairs& doubles, + const CCPairBuilder& builder, const Info& info); ///@param world @@ -943,7 +944,8 @@ class CCPotentials { ///@param[out] \f$ ( 4_2 - 2_1 - 2_2 + _1 ) \f$ /// Q-Projector is not applied, sign is correct static vector_real_function_3d - s4c(World& world, std::vector external_index, const CC_vecfunction& singles, const Pairs& doubles, const Info& info); + s4c(World& world, std::vector external_index, const CC_vecfunction& singles, const Pairs& doubles, + const CCPairBuilder& builder, const Info& info); // update the intermediates void update_intermediates(const CC_vecfunction& t) { diff --git a/src/madness/chem/CCStructures.h b/src/madness/chem/CCStructures.h index 6838eaed971..9291e455364 100644 --- a/src/madness/chem/CCStructures.h +++ b/src/madness/chem/CCStructures.h @@ -1420,7 +1420,7 @@ class CCPairBuilder { World& world; const Info& info; CC_vecfunction gs_singles, ex_singles; - CalcType ctype=CT_MP2; + CalcType ctype=CT_UNDEFINED; };