Skip to content

Commit

Permalink
gs cc2 works with macrotasks and the water dimer, lrcc2 to go
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Oct 30, 2024
1 parent 734cedb commit a480733
Show file tree
Hide file tree
Showing 4 changed files with 57 additions and 33 deletions.
20 changes: 11 additions & 9 deletions src/madness/chem/CC2.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -701,9 +701,7 @@ CC2::solve_cc2(CC_vecfunction& singles, Pairs<CCPair>& 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");
Expand All @@ -719,8 +717,8 @@ CC2::solve_cc2(CC_vecfunction& singles, Pairs<CCPair>& 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");
Expand Down Expand Up @@ -771,16 +769,20 @@ CC2::solve_cc2(CC_vecfunction& singles, Pairs<CCPair>& 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<parameters.dconv_6D();

// check if singles converged
const bool singles_converged = iterate_cc2_singles(world, singles, doubles, info);
timer1.tag("computing cc2 singles");

// check if energy converged
const double omega_new = CCPotentials::compute_cc2_correlation_energy(world, singles, doubles, info);
Expand Down Expand Up @@ -1188,7 +1190,7 @@ bool CC2::iterate_singles(World& world, CC_vecfunction& singles, const CC_vecfun
// else if(ctype==CT_LRCC2) update_reg_residues_ex(world, singles2,singles,ex_doubles, info);

if (world.rank()==0) CCPotentials::print_convergence(singles.name(0),rmsresidual,
rmsresidual,omega-old_omega,iter);
maxresidual,omega-old_omega,iter);
converged = (R2vector_error < info.parameters.dconv_3D());

// time.info();
Expand Down Expand Up @@ -1274,8 +1276,8 @@ CC2::initialize_pairs(Pairs<CCPair>& 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);

Expand Down Expand Up @@ -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);
Expand Down
47 changes: 30 additions & 17 deletions src/madness/chem/CCPotentials.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<CCPair>& doubles, const Info& info)
CCPotentials::compute_cc2_correlation_energy(World& world, const CC_vecfunction& singles, const Pairs<CCPair>& 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<CCPair> pair_vec=Pairs<CCPair>::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
Expand Down
6 changes: 3 additions & 3 deletions src/madness/chem/CCPotentials.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -279,7 +279,7 @@ class CCPotentials {
/// @param info
static double
compute_cc2_correlation_energy(World& world, const CC_vecfunction& singles, const Pairs<CCPair>& doubles,
const Info& info);
const Info& info, const std::string msg="");


static double
Expand Down
17 changes: 13 additions & 4 deletions src/madness/chem/CCStructures.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -597,7 +597,8 @@ MacroTaskSinglesPotentialEx::operator()(const std::vector<int>& 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) {
Expand Down Expand Up @@ -658,10 +659,18 @@ MacroTaskComputeCorrelationEnergy::operator()(const std::vector<CCPair>& pairs,
const Info& info) const {
World &world = pairs[0].function().world();
auto result=scalar_result_vector<double>(world,pairs.size());
CalcType ctype=pairs[0].ctype;
for (int i=0; i<pairs.size(); ++i) {
// 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);
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;
}
Expand Down

0 comments on commit a480733

Please sign in to comment.