Skip to content

Commit

Permalink
some cleanup for LRCC2
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Nov 14, 2024
1 parent 84c23ef commit 81497a4
Show file tree
Hide file tree
Showing 6 changed files with 463 additions and 206 deletions.
78 changes: 48 additions & 30 deletions src/madness/chem/CC2.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ CC2::solve() {
// check if all pair functions have been loaded or computed
auto all_pairs_exist = [](const Pairs<CCPair>& pairs) {
for (const auto& p : pairs.allpairs) {
if (p.second.functions.size()==0) return false;
if (not p.second.function().is_initialized()) return false;
if (not p.second.constant_part.is_initialized()) return false;
}
Expand Down Expand Up @@ -386,7 +387,7 @@ double CC2::solve_mp2_coupled(Pairs<CCPair>& doubles, Info& info) {
for (auto& c : pair_vec) {
MADNESS_CHECK_THROW(c.constant_part.is_initialized(), "could not find constant part");
// constant part is zero-order guess for pair.function
if (not c.function().is_initialized()) c.update_u(c.constant_part);
if (not c.function_exists()) c.update_u(c.constant_part);
}

} else {
Expand All @@ -406,15 +407,14 @@ double CC2::solve_mp2_coupled(Pairs<CCPair>& doubles, Info& info) {

// transform vector back to Pairs structure
for (size_t i = 0; i < pair_vec.size(); i++) {
auto pair=pair_vec[i];
auto& pair=pair_vec[i];
pair.constant_part = result_vec[i];
pair.constant_part.truncate().reduce_rank();
pair.constant_part.print_size("constant_part");
pair.function().truncate().reduce_rank();
save(pair.constant_part, pair.name() + "_const");

// initialize pair function to constant part
if (not pair.function().is_initialized()) pair.update_u(pair.constant_part);
pair.update_u(pair.constant_part);
save(pair.constant_part, pair.name() + "_const");

}
}
Expand Down Expand Up @@ -654,6 +654,8 @@ CC2::iterate_lrcc2_pairs(World& world, const CC_vecfunction& cc2_s,
cc2_s.reconstruct();
lrcc2_s.reconstruct();

if (1) {

// make new constant part
MacroTaskConstantPart tc;
MacroTask task(world, tc);
Expand All @@ -664,9 +666,17 @@ CC2::iterate_lrcc2_pairs(World& world, const CC_vecfunction& cc2_s,
pair_vec[i].constant_part=cp[i];
save(pair_vec[i].constant_part, pair_vec[i].name() + "_const");
}
} else {
print("reading LRCC2 constant part from file");
for (int i=0; i<pair_vec.size(); ++i) {
real_function_6d tmp=real_factory_6d(world);
load(tmp, pair_vec[i].name() + "_const");
pair_vec[i].constant_part=tmp;
}
}

// if no function has been computed so far use the constant part (first iteration)
for (auto& pair : pair_vec) if (not pair.function().is_initialized()) pair.update_u(pair.constant_part);
for (auto& pair : pair_vec) if (not pair.function_exists()) pair.update_u(pair.constant_part);

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");
Expand Down Expand Up @@ -862,13 +872,13 @@ CC2::solve_lrcc2(Pairs<CCPair>& gs_doubles, const CC_vecfunction& gs_singles, co
Pairs<CCPair> ex_doubles;
bool found_lrcc2d = initialize_pairs(ex_doubles, EXCITED_STATE, CT_LRCC2, gs_singles, ex_singles, excitation, info);

// if ((not found_singles) and found_lrcc2d) {
if ((not found_singles) and found_lrcc2d) {
iterate_lrcc2_singles(world, gs_singles, gs_doubles, ex_singles, ex_doubles, info);

std::string filename1=singles_name(CT_LRCC2,ex_singles.type,excitation);
if (world.rank()==0) print("saving singles to disk",filename1);
ex_singles.save_restartdata(world,filename1);
// }
}

if (not info.parameters.no_compute_lrcc2()) {
for (size_t iter = 0; iter < parameters.iter_max(); iter++) {
Expand Down Expand Up @@ -1273,37 +1283,45 @@ CC2::initialize_singles_to_zero(World& world, const CalcType& ctype, const FuncT


bool
CC2::initialize_pairs(Pairs<CCPair>& pairs, const CCState ftype, const CalcType ctype, const CC_vecfunction& tau,
const CC_vecfunction& x, const size_t excitation, const Info& info) const {
MADNESS_ASSERT(tau.type == PARTICLE);
MADNESS_ASSERT(x.type == RESPONSE);
CC2::initialize_pairs(Pairs<CCPair>& pairs, const CCState ftype, const CalcType ctype,
const CC_vecfunction& gs_singles, const CC_vecfunction& ex_singles,
const size_t excitation, const Info& info) const {

MADNESS_ASSERT(gs_singles.type == PARTICLE);
MADNESS_ASSERT(ex_singles.type == RESPONSE);
MADNESS_ASSERT(pairs.empty());

auto builder=CCPairBuilder(world,info);
builder.set_gs_singles(gs_singles).set_ex_singles(ex_singles).set_ctype(ctype);

std::string name1 = CCPair(0, 0, ftype, ctype).name();
if (world.rank()==0) print("initializing doubles",assign_name(ftype), " --- reading from file(s)",name1);

bool restarted = false;

for (size_t i = parameters.freeze(); i < CCOPS.mo_ket().size(); i++) {
for (size_t j = i; j < CCOPS.mo_ket().size(); j++) {

std::string name = CCPair(i, j, ftype, ctype).name();
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(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);
CCPair pair=builder.make_bare_pair_from_file(i,j);

// std::string name = CCPair(i, j, ftype, ctype).name();
// 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(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, pair);
}
}
// it's a restarted calculation if all pairs are found on disk
bool restarted=std::all_of(pairs.allpairs.begin(),pairs.allpairs.end(),
[](const auto& datum){return datum.second.functions.size()>0 && datum.second.function().is_initialized();});

return restarted;
}

Expand Down
45 changes: 38 additions & 7 deletions src/madness/chem/CCPotentials.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,6 +100,7 @@ CCPair CCPotentials::make_pair_cc2(World& world, const real_function_6d& u, cons

// construct Q12 f12 |ij>
if (compute_Q12_f12) {
timer t1(world);
auto phi=info.mo_ket;
auto phi_bra=info.mo_bra;
auto t=make_full_t_intermediate(gs_singles,info).get_vecfunction();
Expand All @@ -110,6 +111,7 @@ CCPair CCPotentials::make_pair_cc2(World& world, const real_function_6d& u, cons
CCPairFunction<double,6> fij(f12, t[i], t[j]);
std::vector<CCPairFunction<double,6>> tmp=Q12(std::vector<CCPairFunction<double,6>>(1,fij));
functions+=tmp;
t1.tag("make low-rank parts in make_pair_cc2");
} else {
if (world.rank()==0) print("skipping the computation of Q12 f12 |ij> in make_pair_cc2");
}
Expand All @@ -133,6 +135,7 @@ CCPair CCPotentials::make_pair_lrcc2(World& world, const CalcType& ctype, const
auto functions=std::vector<cpT>(1,cpT(u));

if (compute_Q12_f12) {
timer t1(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()));
Expand Down Expand Up @@ -160,7 +163,9 @@ CCPair CCPotentials::make_pair_lrcc2(World& world, const CalcType& ctype, const
auto f_tt=std::vector<cpT>(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
t1.tag("make low-rank parts in make_pair_lrcc2");
functions=consolidate(functions);
t1.tag("consolidate low-rank parts in make_pair_lrcc2");
}

CCPair pair(i, j, EXCITED_STATE, ctype, functions);
Expand Down Expand Up @@ -2413,6 +2418,12 @@ CCPotentials::get_CC2_singles_potential_ex(World& world, const CC_vecfunction& g
MacroTaskSinglesPotentialEx taskex;
MacroTask<MacroTaskSinglesPotentialEx> mtaskex(world,taskex,taskq);

info.reconstruct();
for (auto& p : gs_doubles.allpairs) p.second.reconstruct();
for (auto& p : ex_doubles.allpairs) p.second.reconstruct();
gs_singles.reconstruct();
ex_singles.reconstruct();

std::vector<int> result_index(ex_singles.size());
for (int i=0; i<result_index.size(); ++i) result_index[i]=i+info.parameters.freeze();
auto [fock_residue, dum1] = mtaskex(result_index, gs_singles, gs_doubles_vec, ex_singles, ex_doubles_vec, int(POT_F3D_), info);
Expand Down Expand Up @@ -2551,6 +2562,8 @@ CCPotentials::potential_singles_gs(World& world, const std::vector<int>& 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);

if (name == POT_F3D_) {
result = fock_residue_closed_shell(world, singles_external, info);
} else if (name == POT_ccs_) { // or S3c
Expand All @@ -2574,9 +2587,9 @@ CCPotentials::potential_singles_gs(World& world, const std::vector<int>& result_
// << get_size(world,result_s4a) << " (GB), " << timer.current_time().first << "s (wall), " << timer.current_time().second << "s (cpu)\n";
// result = add(world,result_s2b,result_s4a);
// returns the s2b potential (unprojected)
std::tie(result,intermediate) = s2b(world, result_index, singles, doubles, info);
std::tie(result,intermediate) = s2b(world, result_index, singles, doubles, builder, info);
} else if (name == POT_s2c_) {
std::tie(result,intermediate) = s2c(world, result_index, singles, doubles, info);
std::tie(result,intermediate) = s2c(world, result_index, singles, doubles, builder, info);
} else if (name == POT_s4a_) {
error("potential_singles: Demanded s4a potential -> this is calculated along with the s2b potential");
} else if (name == POT_s4b_) {
Expand Down Expand Up @@ -2662,6 +2675,8 @@ CCPotentials::potential_singles_ex(World& world, const std::vector<int> 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);

vector_real_function_3d result, intermediate;

// collect all singles that correspond to external lines
Expand Down Expand Up @@ -2694,9 +2709,9 @@ CCPotentials::potential_singles_ex(World& world, const std::vector<int> 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, info);
std::tie(result, intermediate) = s2b(world, result_index, singles_ex, doubles_ex, builder, info);
} else if (name == POT_s2c_) {
std::tie(result, intermediate) = s2c(world, result_index, singles_ex, doubles_ex, info);
std::tie(result, intermediate) = s2c(world, result_index, singles_ex, doubles_ex, builder, info);
} else if (name == POT_s4a_) {
error("potential_singles: Demanded s4a potential -> this is calculated from the s2b potential");
} else if (name == POT_s4b_) {
Expand Down Expand Up @@ -3193,9 +3208,10 @@ CCPotentials::x_s4c(const CC_vecfunction& x, const CC_vecfunction& t, const Pair

std::tuple<madness::vector_real_function_3d, madness::vector_real_function_3d>
CCPotentials::s2b(World& world, const std::vector<int> external_indices, const CC_vecfunction& singles,
const Pairs<CCPair>& doubles, const Info& info)
const Pairs<CCPair>& doubles, const CCPairBuilder& builder, const Info& info)
{
vector_real_function_3d result;
// auto builder=PairBuilder(world,info).set_gs_singles(singles);

// 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)
Expand All @@ -3210,7 +3226,12 @@ CCPotentials::s2b(World& world, const std::vector<int> external_indices, const C
real_function_3d resulti_r = real_factory_3d(world);
for (const auto& ktmp : singles.functions) {
const size_t k = ktmp.first;
std::vector<CCPairFunction<double,6>> uik = get_pair_function(doubles, i, k);
// next line might swap particles 1 and 2
std::vector<CCPairFunction<double,6>> uik1 = get_pair_function(doubles, i, k);
// make a ful pair function using the pure 6D part uik and complete with the rest
CCPair pair_ik=builder.complete_pair_with_low_rank_parts(builder.make_pair(i,k,uik1));
std::vector<CCPairFunction<double,6>> uik = pair_ik.functions;

// check if the first function in the vector is really the pure 6D part
MADNESS_ASSERT(uik[0].is_pure());
if (recalc_u_part) {
Expand All @@ -3235,15 +3256,25 @@ CCPotentials::s2b(World& world, const std::vector<int> external_indices, const C
}

std::tuple<madness::vector_real_function_3d, madness::vector_real_function_3d>
CCPotentials::s2c(World& world, const std::vector<int> external_index, const CC_vecfunction& singles, const Pairs<CCPair>& doubles, const Info& info) {
CCPotentials::s2c(World& world, const std::vector<int> external_index, const CC_vecfunction& singles,
const Pairs<CCPair>& doubles1, const CCPairBuilder& builder, const Info& info) {
vector_real_function_3d result;
// auto builder=PairBuilder(world,info).set_gs_singles(singles);

// 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_,false);
vector_real_function_3d result_u_return;
bool recalc_u_part = false;
if (result_u.empty()) recalc_u_part = true;
auto g12=CCConvolutionOperator<double,3>(world,OT_G12,info.parameters);

// add the regularization terms Q12 f12 |xy> to the pair
Pairs<CCPair> doubles=doubles1;
for (auto& itmp : doubles.allpairs) {
auto& pair=itmp.second;
pair=builder.complete_pair_with_low_rank_parts(pair);
}

// for (const auto& itmp : singles.functions) {
for (const int i : external_index) {
// const size_t i = itmp.first;
Expand Down
8 changes: 6 additions & 2 deletions src/madness/chem/CCPotentials.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -892,12 +892,14 @@ class CCPotentials {
///@param external_indices
///@param[in] singles:CC_vecfunction fof type response or particle (depending on this the correct intermediates will be used) the functions themselves are not needed
///@param[in] doubles:Pairs of CC_Pairs (GS or Response)
///@param builder
///@param info
///@param[out] \f$ \sum_k( 2<k|g|uik>_2 - <k|g|uik>_1 ) \f$
/// Q-Projector is not applied, sign is correct
/// if the s2b potential has already been calculated it will be loaded from the intermediate_potentials structure
static std::tuple<madness::vector_real_function_3d, madness::vector_real_function_3d>
s2b(World& world, std::vector<int> external_indices, const CC_vecfunction& singles, const Pairs<CCPair>& doubles, const Info& info);
s2b(World& world, std::vector<int> external_indices, const CC_vecfunction& singles, const Pairs<CCPair>& doubles, const
CCPairBuilder& builder, const Info& info);

/// result: -\sum_k( <l|kgi|ukl>_2 - <l|kgi|ukl>_1)

Expand All @@ -908,11 +910,13 @@ class CCPotentials {
///@param external_index
///@param[in] singles:CC_vecfunction fof type response or particle (depending on this the correct intermediates will be used) the functions themselves are not needed
///@param[in] doubles:Pairs of CC_Pairs (GS or Response)
///@param builder
///@param info
///@param[out] \f$ -\sum_k( <l|kgi|ukl>_2 - <l|kgi|ukl>_1) \f$
/// Q-Projector is not applied, sign is correct
static std::tuple<vector<Function<double, 3>>, vector<Function<double, 3>>>
s2c(World& world, std::vector<int> external_index, const CC_vecfunction& singles, const Pairs<CCPair>& doubles, const Info& info);
s2c(World& world, std::vector<int> external_index, const CC_vecfunction& singles, const Pairs<CCPair>& doubles, const
CCPairBuilder& builder, const Info& info);

/// the S4a potential can be calcualted from the S2b potential
/// result is \f$ s4a_i = - <l|s2b_i>*|tau_l> \f$
Expand Down
Loading

0 comments on commit 81497a4

Please sign in to comment.