diff --git a/src/madness/chem/CC2.cc b/src/madness/chem/CC2.cc index 262b53faea8..e7c3e0bb7bc 100644 --- a/src/madness/chem/CC2.cc +++ b/src/madness/chem/CC2.cc @@ -701,9 +701,7 @@ CC2::solve_cc2(CC_vecfunction& singles, Pairs& doubles, Info& info) cons MADNESS_ASSERT(singles.type == PARTICLE); CCTimer time(world, "CC2 Ground State"); - double omega = CCPotentials::compute_cc2_correlation_energy(world, singles, doubles, info); - if (world.rank() == 0) - std::cout << std::fixed << std::setprecision(10) << "Current Correlation Energy = " << omega << "\n"; + 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"); @@ -719,8 +717,8 @@ 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); + // 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++) { CCTimer time_miter(world, "Macroiteration " + std::to_string(int(iter)) + " of CC2"); @@ -771,16 +769,20 @@ CC2::solve_cc2(CC_vecfunction& singles, Pairs& doubles, Info& info) cons // save latest iteration if (world.rank()==0) print("saving latest iteration to file"); for (const auto& pair : pair_vec) { + pair.constant_part.reconstruct(); save(pair.constant_part, pair.name() + "_const"); + pair.function().reconstruct(); save(pair.function(), pair.name()); singles.save_restartdata(world,madness::name(singles.type)); } + timer1.tag("saving cc2 singles and doubles"); auto [rmsrnorm,maxrnorm]=CCPotentials::residual_stats(residual); bool doubles_converged=rmsrnorm& pairs, const CCState ftype, const CalcType 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); - if (ctype==CT_CC2) tmp=CCPotentials::make_pair_cc2(utmp, tau, i, j, info); + 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); tmp.constant_part = const_part; pairs.insert(i, j, tmp); @@ -1315,7 +1317,7 @@ void CC2::update_reg_residues_gs(World& world, const CC_vecfunction& singles, Pa 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); + 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); diff --git a/src/madness/chem/CCPotentials.cc b/src/madness/chem/CCPotentials.cc index 146d5666e54..9a0d3588f7d 100644 --- a/src/madness/chem/CCPotentials.cc +++ b/src/madness/chem/CCPotentials.cc @@ -372,24 +372,37 @@ CCPotentials::compute_pair_correlation_energy(World& world, const CCPair& u, } double -CCPotentials::compute_cc2_correlation_energy(World& world, const CC_vecfunction& singles, const Pairs& doubles, const Info& info) +CCPotentials::compute_cc2_correlation_energy(World& world, const CC_vecfunction& singles, const Pairs& doubles, const Info& info, const std::string msg) { - MADNESS_ASSERT(singles.type == PARTICLE); - CCTimer time(world, "Computing CC2 Correlation Energy"); - // output.section("Computing CC2 Correlation Energy"); - double result = 0.0; - for (const auto& tmp : doubles.allpairs) { - const size_t i = tmp.second.i; - const size_t j = tmp.second.j; - const double omega = compute_pair_correlation_energy(world, tmp.second, singles, info); - result += omega; - if (world.rank() == 0) - std::cout << std::fixed << "omega " << i << j << " =" << std::setprecision(10) << omega << "\n"; - } - if (world.rank() == 0) std::cout << std::fixed << "sum " << " =" << std::setprecision(10) << result << "\n"; - - time.info(); - return result; + auto triangular_map=PairVectorMap::triangular_map(info.parameters.freeze(),info.mo_ket.size()); + std::vector pair_vec=Pairs::pairs2vector(doubles,triangular_map); + MacroTaskComputeCorrelationEnergy t; + MacroTask task1(world, t); + auto pair_energies=task1(pair_vec, singles, info); + // pair_energies is now scattered over the universe + + double total_energy=0.0; + for ( auto& pair_energy : pair_energies) total_energy += pair_energy.get(); + // pair_energy.get() invokes a broadcast from rank 0 to all other ranks + + if (not msg.empty() and world.rank()==0) printf("%s %12.8f\n", msg.c_str(), total_energy); + return total_energy; +// MADNESS_ASSERT(singles.type == PARTICLE); +// CCTimer time(world, "Computing CC2 Correlation Energy"); +// // output.section("Computing CC2 Correlation Energy"); +// double result = 0.0; +// for (const auto& tmp : doubles.allpairs) { +// const size_t i = tmp.second.i; +// const size_t j = tmp.second.j; +// const double omega = compute_pair_correlation_energy(world, tmp.second, singles, info); +// result += omega; +// if (world.rank() == 0) +// std::cout << std::fixed << "omega " << i << j << " =" << std::setprecision(10) << omega << "\n"; +// } +// if (world.rank() == 0) std::cout << std::fixed << "sum " << " =" << std::setprecision(10) << result << "\n"; +// +// time.info(); +// return result; } double diff --git a/src/madness/chem/CCPotentials.h b/src/madness/chem/CCPotentials.h index 37dea590591..01c6efbfd5f 100644 --- a/src/madness/chem/CCPotentials.h +++ b/src/madness/chem/CCPotentials.h @@ -212,11 +212,11 @@ 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 = true); + 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=true); + 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, @@ -279,7 +279,7 @@ class CCPotentials { /// @param info static double compute_cc2_correlation_energy(World& world, const CC_vecfunction& singles, const Pairs& doubles, - const Info& info); + const Info& info, const std::string msg=""); static double diff --git a/src/madness/chem/CCStructures.cc b/src/madness/chem/CCStructures.cc index a125ff4e055..f8a9f9cc178 100644 --- a/src/madness/chem/CCStructures.cc +++ b/src/madness/chem/CCStructures.cc @@ -597,7 +597,8 @@ MacroTaskSinglesPotentialEx::operator()(const std::vector& result_index, for (auto& x : doubles_gs1.allpairs) { auto& tau=x.second; MADNESS_CHECK_THROW(tau.functions.size()==1,"doubles in MacroTaskSinglesPotentialsEx should only contain one function"); - tau=CCPotentials::make_pair_cc2(tau.function(),singles_gs,tau.i,tau.j,info); + 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); } // the doubles currently only contain the full 6d function -> complete it with the Q12 f12 |ti tj> part for (auto& x : doubles_ex1.allpairs) { @@ -658,10 +659,18 @@ MacroTaskComputeCorrelationEnergy::operator()(const std::vector& pairs, const Info& info) const { World &world = pairs[0].function().world(); auto result=scalar_result_vector(world,pairs.size()); + CalcType ctype=pairs[0].ctype; 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); - result[i]=CCPotentials::compute_pair_correlation_energy(world,pair,singles_gs,info); + if (ctype==CT_MP2) { + // when serialized the Qf12 |ij> 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); + 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); + result[i]=CCPotentials::compute_pair_correlation_energy(world,pair,singles_gs,info); + } else { + MADNESS_EXCEPTION("MacroTaskComputeCorrelationEnergy: unknown ctype",1); + } } return result; }