Skip to content

Commit

Permalink
LRCC2 seems to work for water dimer, but singles are very slow
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Nov 5, 2024
1 parent df7d831 commit 84c23ef
Show file tree
Hide file tree
Showing 7 changed files with 124 additions and 78 deletions.
88 changes: 37 additions & 51 deletions src/madness/chem/CC2.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,7 @@ CC2::solve() {
}

if (need_cc2) {
if (world.rank()==0) print_header2("computing the CC2 correlation energy");
// check if singles or/and doubles to restart are there
cc2singles=initialize_singles(world,CT_CC2, PARTICLE);

Expand All @@ -112,11 +113,11 @@ CC2::solve() {
output_calc_info_schema("cc2",cc2_energy);
}

output.section(assign_name(CT_CC2) + " Calculation Ended !");
if (world.rank() == 0) {
printf_msg_energy_time("CC2 correlation energy",cc2_energy,wall_time());
std::cout << std::fixed << std::setprecision(10) << " CC2 Correlation Energy =" << cc2_energy << "\n";
}
if (world.rank()==0) print_header2("end computing the CC2 correlation energy");
}

if (ctype == CT_LRCCS) {
Expand Down Expand Up @@ -264,16 +265,17 @@ CC2::solve() {
} else if (ctype == CT_LRCC2) {
CCTimer time(world, "Whole LRCC2 Calculation");

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");
// iterate_ccs_singles(vccs[0], info);
// 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");
// print_header1("skipping reiteration of CC2 singles for the potential");
iterate_cc2_singles(world, cc2singles, cc2pairs, info);
if (world.rank()==0) print_header3("end reiterating CC2 singles");
}
Expand Down Expand Up @@ -365,8 +367,8 @@ std::vector<CC_vecfunction> CC2::solve_ccs() const {
1);
result.push_back(excitations[x]);
}
if (world.rank()==0) print_header3("Solution of the CCS equations");
tdhf->analyze(result);
// if (world.rank()==0) print_header3("Solution of the CCS equations");
// tdhf->analyze(result);
return result;
}

Expand Down Expand Up @@ -647,6 +649,11 @@ CC2::iterate_lrcc2_pairs(World& world, const CC_vecfunction& cc2_s,
auto triangular_map=PairVectorMap::triangular_map(info.parameters.freeze(),info.mo_ket.size());
auto pair_vec=Pairs<CCPair>::pairs2vector(lrcc2_d,triangular_map);

for (auto& p: pair_vec) p.reconstruct();
info.reconstruct();
cc2_s.reconstruct();
lrcc2_s.reconstruct();

// make new constant part
MacroTaskConstantPart tc;
MacroTask task(world, tc);
Expand All @@ -656,7 +663,6 @@ CC2::iterate_lrcc2_pairs(World& world, const CC_vecfunction& cc2_s,
for (int i=0; i<pair_vec.size(); ++i) {
pair_vec[i].constant_part=cp[i];
save(pair_vec[i].constant_part, pair_vec[i].name() + "_const");

}

// if no function has been computed so far use the constant part (first iteration)
Expand All @@ -670,10 +676,7 @@ CC2::iterate_lrcc2_pairs(World& world, const CC_vecfunction& cc2_s,
Pairs<real_function_6d> coupling=compute_local_coupling(pair_vec, info);
auto coupling_vec=Pairs<real_function_6d>::pairs2vector(coupling,triangular_map);
reconstruct(world,coupling_vec);
for (auto& p : pair_vec) {
p.constant_part.reconstruct();
p.function().reconstruct();
}
for (auto& p : pair_vec) p.reconstruct();

if (info.parameters.debug()) print_size(world, coupling_vec, "couplingvector");

Expand Down Expand Up @@ -846,36 +849,34 @@ CC2::solve_lrcc2(Pairs<CCPair>& gs_doubles, const CC_vecfunction& gs_singles, co
/// read LRCC2 singles from file or use the CIS vectors as guess
auto ex_singles=initialize_singles(world,CT_LRCC2,RESPONSE, excitation);
bool found_singles=(not (ex_singles.get_vecfunction().empty()));
if (not found_singles) {
// if (not found_singles) {
if (world.rank()==0) print("using CIS vectors as guess for the LRCC2 singles");
ex_singles = copy(cis);
iterate_ccs_singles(ex_singles, info);

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

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 filename=singles_name(CT_LRCC2,ex_singles.type,excitation);
if (world.rank()==0) print("saving singles to disk",filename);
ex_singles.save_restartdata(world,filename);
}
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++) {
if (world.rank()==0) print_header2("Macroiteration " + std::to_string(int(iter)) + " of LRCC2 for excitation energy "+std::to_string(ex_singles.omega));
// update_reg_residues_ex(world, gs_singles, ex_singles, ex_doubles, info);

bool dconv = iterate_lrcc2_pairs(world, gs_singles, ex_singles, ex_doubles, info);
bool sconv = iterate_lrcc2_singles(world, gs_singles, gs_doubles, ex_singles, ex_doubles, info);

// update_reg_residues_ex(world, gs_singles, ex_singles, ex_doubles, info);

std::string filename=singles_name(CT_LRCC2,ex_singles.type,excitation);
if (world.rank()==0) print("saving singles to disk",filename);
ex_singles.save_restartdata(world,filename);
Expand Down Expand Up @@ -1287,35 +1288,20 @@ CC2::initialize_pairs(Pairs<CCPair>& 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;
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);

// } 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");
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);
}
}
return restarted;
Expand Down
19 changes: 14 additions & 5 deletions src/madness/chem/CCPotentials.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1380,8 +1380,8 @@ std::vector<CCPairFunction<double,6>>
// calculate the regularized potential
real_function_6d V=real_factory_6d(world);
std::vector<CCPairFunction<double,6>> 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);
Expand Down Expand Up @@ -2291,6 +2291,8 @@ CCPotentials::get_CC2_singles_potential_gs(World& world, const CC_vecfunction& s
//
// s2b and s2c return the potential and an intermediate which needs to be stored

for (auto& p : doubles.allpairs) p.second.reconstruct();

// set up taskq
MacroTaskSinglesPotentialGs taskgs;
taskgs.partitioner->min_batch_size=2;
Expand All @@ -2303,7 +2305,7 @@ CCPotentials::get_CC2_singles_potential_gs(World& world, const CC_vecfunction& s
// partition macrotasks over active orbitals
std::vector<int> result_index(singles.size());
for (int i=0; i<result_index.size(); ++i) result_index[i]=i+info.parameters.freeze();
print_header1("not clearing intermediate potentials");
// print_header1("not clearing intermediate potentials");
// info.intermediate_potentials.clear_all();

auto [fock_residue, dum1] = mtaskgs(result_index, singles, doubles_vec, int(POT_F3D_), info);
Expand Down Expand Up @@ -2412,6 +2414,7 @@ CCPotentials::get_CC2_singles_potential_ex(World& world, const CC_vecfunction& g
MacroTask<MacroTaskSinglesPotentialEx> mtaskex(world,taskex,taskq);

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);
auto [Vccs, dum2] = mtaskex(result_index, gs_singles, gs_doubles_vec, ex_singles, ex_doubles_vec, int(POT_ccs_), info);
auto [Vs2b, imed_s2b] = mtaskex(result_index, gs_singles, gs_doubles_vec, ex_singles, ex_doubles_vec, int(POT_s2b_), info);
Expand Down Expand Up @@ -2831,12 +2834,18 @@ CCPotentials::apply_K_macrotask(World& world, const std::vector<real_function_3d
real_convolution_3d g12 = CoulombOperator(world, parameters.lo(), parameters.thresh_poisson());
g12.particle() = particle;
for (size_t k = 0; k < mo_ket.size(); k++) {
// print("k",k);
real_function_6d copyu = copy(u);
// copyu.print_size("copyu");
real_function_6d X = (multiply(copyu, copy(mo_bra[k]), particle)).truncate();
// real_function_6d Y=(*poisson)(X);
// X.print_size("X");
real_function_6d Y = g12(X); // overwrite X to save space
result += (multiply(copy(Y), copy(mo_ket[k]),
particle)).truncate(); // this will destroy X, but I d not intend to use it again so I choose here to save this copy
// Y.print_size("Y");
// Y.print_tree();
auto tmp=(multiply(copy(Y), copy(mo_ket[k]),particle)).truncate(); // this will destroy X, but I d not intend to use it again so I choose here to save this copy
// tmp.print_size("tmp");
result += tmp;
}
return result.truncate(parameters.tight_thresh_3D()*3.0).reduce_rank(parameters.tight_thresh_6D()*3.0);
}
Expand Down
3 changes: 2 additions & 1 deletion src/madness/chem/CCPotentials.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,9 +77,10 @@ class CCPotentials {
ar & f;
if (do_print) f.print_size(name);
if (f.is_compressed()) {
if (do_print) print("function is compressed -- reconstructing");
if (world.rank()==0 and do_print) print("function is compressed -- reconstructing");
f.change_tree_state(reconstructed);
if (do_print) f.print_size(name+" reconstructed");
save(f, name);
}
f.set_thresh(FunctionDefaults<NDIM>::get_thresh());
f.truncate();
Expand Down
31 changes: 31 additions & 0 deletions src/madness/chem/CCStructures.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -592,6 +592,9 @@ struct CC_vecfunction : public archive::ParallelSerializableObject {
return tmp;
}

void reconstruct() const {
for (auto& x : functions) x.second.function.reconstruct();
}

//madness::CC_vecfunction
//CC_vecfunction::copy() const {
Expand Down Expand Up @@ -1053,6 +1056,14 @@ class CCPair : public archive::ParallelSerializableObject {
if (cexist) ar & constant_part;
}

/// reconstruct constant part and all functions
void reconstruct() const {
constant_part.reconstruct();
for (auto& f : functions) {
if (f.is_assigned() and f.is_pure()) f.get_function().reconstruct();
}
}

bool load_pair(World& world) {
std::string fname=this->name();
bool exists = archive::ParallelInputArchive<archive::BinaryFstreamInputArchive>::exists(world, fname.c_str());
Expand Down Expand Up @@ -1136,6 +1147,16 @@ struct CCIntermediatePotentials {
Function<double,3>
operator()(const CCFunction<double,3>& f, const PotentialType& type, const bool throw_if_empty) const;

void reconstruct() const {
madness::reconstruct(current_s2b_potential_ex_);
madness::reconstruct(current_s2b_potential_gs_);
madness::reconstruct(current_s2c_potential_ex_);
madness::reconstruct(current_s2c_potential_gs_);
madness::reconstruct(current_singles_potential_ex_);
madness::reconstruct(current_singles_potential_gs_);
madness::reconstruct(unprojected_cc2_projector_response_);
}

/// deltes all stored potentials
void clear_all() {
current_singles_potential_gs_.clear();
Expand Down Expand Up @@ -1252,6 +1273,16 @@ struct Info {
return result;
}

void reconstruct() const {
madness::reconstruct(mo_bra);
madness::reconstruct(mo_ket);
R_square.reconstruct();
madness::reconstruct(U1);
U2.reconstruct();
intermediate_potentials.reconstruct();

}

/// customized function to store this to the cloud

/// functions and constant_part can be very large and we want to split them and store them in different records
Expand Down
26 changes: 14 additions & 12 deletions src/madness/chem/TDHF.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,14 +174,14 @@ void TDHF::prepare_calculation() {
if (get_calcparam().do_localize()) {
Fock<double,3> F(world,get_reference().get());
F_occ = F(get_active_mo_bra(), get_active_mo_ket());
for (size_t i = 0; i < get_active_mo_ket().size(); ++i) {
msg << std::scientific << std::setprecision(10);
msg << "F(" << i << "," << i << ")=" << F_occ(i, i) << "\n";
if (std::fabs(get_orbital_energy(i + parameters.freeze()) - F_occ(i, i)) > 1.e-5) {
msg << "eps(" << i << ")=" << get_orbital_energy(i) << " | diff="
<< get_orbital_energy(i + parameters.freeze()) - F_occ(i, i) << "\n";
}
}
// for (size_t i = 0; i < get_active_mo_ket().size(); ++i) {
// msg << std::scientific << std::setprecision(10);
// msg << "F(" << i << "," << i << ")=" << F_occ(i, i) << "\n";
// if (std::fabs(get_orbital_energy(i + parameters.freeze()) - F_occ(i, i)) > 1.e-5) {
// msg << "eps(" << i << ")=" << get_orbital_energy(i) << " | diff="
// << get_orbital_energy(i + parameters.freeze()) - F_occ(i, i) << "\n";
// }
// }
} else {
F_occ = Tensor<double>(get_active_mo_bra().size(), get_active_mo_ket().size());
F_occ *= 0.0;
Expand Down Expand Up @@ -343,10 +343,12 @@ std::vector<CC_vecfunction> TDHF::solve_cis() const {
bool skip_solve=(parameters.restart()=="no_compute") or (converged_roots.size()>=parameters.nexcitations());

if (skip_solve) {
print("skipping the solution of the CIS equations");
if (parameters.restart()=="no_compute") print(" -> no_compute is set");
if (converged_roots.size()>=parameters.nexcitations())
print(" -> number of converged excitations from disk is sufficient:", converged_roots.size());
if (world.rank()==0) {
print("skipping the solution of the CIS equations");
if (parameters.restart()=="no_compute") print(" -> no_compute is set");
if (converged_roots.size()>=parameters.nexcitations())
print(" -> number of converged excitations from disk is sufficient:", converged_roots.size());
}
} else {
for (size_t macrocycle = 0; macrocycle < 1; ++macrocycle) {
//msg.section("CIS Macroiteration " + std::to_string(macrocycle));
Expand Down
29 changes: 20 additions & 9 deletions src/madness/mra/mraimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3291,15 +3291,26 @@ template <typename T, std::size_t NDIM>
else {
// special case: tree has only root node: keep sum coeffs and make zero diff coeffs
if (key.level()==0) {
coeffT result(node.coeff());
coeffT sdcoeff(cdata.v2k,this->get_tensor_type());
sdcoeff(cdata.s0)+=node.coeff();
node.coeff()=sdcoeff;
double snorm=node.coeff().normf();
node.set_dnorm(0.0);
node.set_snorm(snorm);
node.set_norm_tree(snorm);
return Future< std::pair<GenTensor<T>,double> >(std::make_pair(result,node.coeff().normf()));
if (redundant1) {
// with only the root node existing redundant and reconstructed are the same
coeffT result(node.coeff());
double snorm=node.coeff().normf();
node.set_dnorm(0.0);
node.set_snorm(snorm);
node.set_norm_tree(snorm);
return Future< std::pair<GenTensor<T>,double> >(std::make_pair(result,snorm));
} else {
// compress
coeffT result(node.coeff());
coeffT sdcoeff(cdata.v2k,this->get_tensor_type());
sdcoeff(cdata.s0)+=node.coeff();
node.coeff()=sdcoeff;
double snorm=node.coeff().normf();
node.set_dnorm(0.0);
node.set_snorm(snorm);
node.set_norm_tree(snorm);
return Future< std::pair<GenTensor<T>,double> >(std::make_pair(result,node.coeff().normf()));
}

} else { // this is a leaf node
Future<coeffT > result(node.coeff());
Expand Down
Loading

0 comments on commit 84c23ef

Please sign in to comment.