Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/pr-cc-refactoring' into pr-cc-re…
Browse files Browse the repository at this point in the history 
…factoring
  • Loading branch information
@fbischoff
fbischoff committed Oct 16, 2023
2 parents e4ee892 + 3809700 commit 1c3b19d
Show file tree
Hide file tree
Showing 37 changed files with 4,094 additions and 944 deletions.
18 changes: 13 additions & 5 deletions src/apps/mp2/mp2.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,11 +80,19 @@ int main(int argc, char** argv) {

if (world.rank() == 0) printf("\nstarting at time %.1fs\n", wall_time());

const double hf_energy = mp2.get_hf().value();
const double mp2_energy = mp2.value();
if (world.rank() == 0) {
printf("final hf/mp2/total energy %12.8f %12.8f %12.8f\n",
hf_energy, mp2_energy, hf_energy + mp2_energy);
const double hf_energy=mp2.get_hf().value();
const double mp2_energy=mp2.value();
// const double mp2_energy=0.0;
if(world.rank() == 0) {
printf("final hf/mp2/total energy %12.8f %12.8f %12.8f\n",
hf_energy,mp2_energy,hf_energy+mp2_energy);
}
double mp3_correction=mp2.mp3();
print("mp3 correction",mp3_correction);
double mp3_energy=mp3_correction+mp2_energy;
if(world.rank() == 0) {
printf("final hf/mp3/total energy %12.8f %12.8f %12.8f\n",
hf_energy,mp3_energy,hf_energy+mp3_energy);
}
} catch (std::exception& e) {

Expand Down
154 changes: 77 additions & 77 deletions src/madness/chem/CCPotentials.cc
View file

Large diffs are not rendered by default.

8 changes: 4 additions & 4 deletions src/madness/chem/CCPotentials.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -794,9 +794,9 @@ class CCPotentials {

// update the intermediates
void update_intermediates(const CC_vecfunction& t) {
g12.update_elements(mo_bra_, t);
g12->update_elements(mo_bra_, t);
// g12.sanity();
f12.update_elements(mo_bra_, t);
f12->update_elements(mo_bra_, t);
// f12.sanity();
}

Expand Down Expand Up @@ -829,9 +829,9 @@ class CCPotentials {
std::vector<double> orbital_energies_;
/// the coulomb operator with all intermediates
public:
CCConvolutionOperator g12;
std::shared_ptr<CCConvolutionOperator> g12;
/// the f12 operator with all intermediates
CCConvolutionOperator f12;
std::shared_ptr<CCConvolutionOperator> f12;
/// the correlation factor, holds necessary regularized potentials
CorrelationFactor corrfac;
/// Manager for stored intermediate potentials which are s2c, s2b and the whole singles potentials without fock-residue for GS and EX state
Expand Down
69 changes: 32 additions & 37 deletions src/madness/chem/CCStructures.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -322,7 +322,7 @@ CCConvolutionOperator::operator()(const CCFunction& bra, const CCFunction& ket,

real_function_6d CCConvolutionOperator::operator()(const real_function_6d& u, const size_t particle) const {
MADNESS_ASSERT(particle == 1 or particle == 2);
MADNESS_ASSERT(operator_type == OT_G12);
MADNESS_ASSERT(operator_type == OpType::OT_G12);
MADNESS_ASSERT(op);
op->particle() = particle;
return (*op)(u);
Expand All @@ -331,7 +331,7 @@ real_function_6d CCConvolutionOperator::operator()(const real_function_6d& u, co
real_function_3d
CCConvolutionOperator::operator()(const CCFunction& bra, const real_function_6d& u, const size_t particle) const {
MADNESS_ASSERT(particle == 1 or particle == 2);
MADNESS_ASSERT(operator_type == OT_G12);
MADNESS_ASSERT(operator_type == OpType::OT_G12);
MADNESS_ASSERT(op);
const real_function_6d tmp = multiply(copy(u), copy(bra.function), particle);
op->particle() = particle;
Expand Down Expand Up @@ -406,33 +406,47 @@ size_t CCConvolutionOperator::info() const {

SeparatedConvolution<double, 3> *
CCConvolutionOperator::init_op(const OpType& type, const Parameters& parameters) const {
bool debug=false;
bool printme=(world.rank()==0) and debug;
switch (type) {
case OT_G12 : {
if (world.rank() == 0)
case OpType::OT_G12 : {
if (printme)
std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op
<< " and lo=" << parameters.lo << std::endl;
return CoulombOperatorPtr(world, parameters.lo, parameters.thresh_op);
}
case OT_F12 : {
if (world.rank() == 0)
case OpType::OT_F12 : {
if (printme)
std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op
<< " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl;
return SlaterF12OperatorPtr(world, parameters.gamma, parameters.lo, parameters.thresh_op);
}
case OT_SLATER : {
if (world.rank() == 0)
case OpType::OT_SLATER : {
if (printme)
std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op
<< " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl;
return SlaterOperatorPtr(world, parameters.gamma, parameters.lo, parameters.thresh_op);
}
case OT_BSH : {
if (world.rank() == 0)
case OpType::OT_BSH : {
if (printme)
std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op
<< " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl;
return BSHOperatorPtr3D(world, parameters.gamma, parameters.lo, parameters.thresh_op);
}
case OT_ONE : {
if (world.rank() == 0)
case OpType::OT_FG12: {
if (printme)
std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op
<< " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl;
return FGOperatorPtr(world, parameters.gamma, parameters.lo, parameters.thresh_op);
}
case OpType::OT_F2G12: {
if (printme)
std::cout << "Creating " << assign_name(type) << " Operator with thresh=" << parameters.thresh_op
<< " and lo=" << parameters.lo << " and Gamma=" << parameters.gamma << std::endl;
return F2GOperatorPtr(world, parameters.gamma, parameters.lo, parameters.thresh_op);
}
case OpType::OT_ONE : {
if (printme)
std::cout << "Creating " << assign_name(type) << " Operator " << std::endl;
return nullptr;
}
Expand All @@ -444,25 +458,6 @@ CCConvolutionOperator::init_op(const OpType& type, const Parameters& parameters)

}

/// Assigns strings to enums for formated output
std::string
assign_name(const PairFormat& input) {
switch (input) {
case PT_FULL:
return "full";
case PT_DECOMPOSED:
return "decomposed";
case PT_OP_DECOMPOSED:
return "operator-decomposed";
default: {
MADNESS_EXCEPTION("Unvalid enum assignement!", 1);
return "undefined";
}
}
MADNESS_EXCEPTION("assign_name:pairtype, should not end up here", 1);
return "unknown pairtype";
}

/// Assigns strings to enums for formated output
std::string
assign_name(const CCState& input) {
Expand All @@ -484,17 +479,17 @@ assign_name(const CCState& input) {
std::string
assign_name(const OpType& input) {
switch (input) {
case OT_G12:
case OpType::OT_G12:
return "g12";
case OT_F12:
case OpType::OT_F12:
return "f12";
case OT_SLATER:
case OpType::OT_SLATER:
return "slater";
case OT_FG12:
case OpType::OT_FG12:
return "fg12";
case OT_BSH:
case OpType::OT_BSH:
return "bsh";
case OT_ONE:
case OpType::OT_ONE:
return "identity";
default: {
MADNESS_EXCEPTION("Unvalid enum assignement!", 1);
Expand Down
108 changes: 30 additions & 78 deletions src/madness/chem/CCStructures.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,18 +21,6 @@

namespace madness {

/// Operatortypes used by the CCConvolutionOperator Class
enum OpType {
OT_UNDEFINED,
OT_ONE, /// indicates the identity
OT_G12, /// 1/r
OT_SLATER, /// exp(r)
OT_F12, /// 1-exp(r)
OT_FG12, /// (1-exp(r))/r
OT_F212, /// (1-exp(r))^2
OT_BSH /// exp(r)/r
};

/// Calculation Types used by CC2
enum CalcType {
CT_UNDEFINED, CT_MP2, CT_CC2, CT_LRCCS, CT_LRCC2, CT_CISPD, CT_ADC2, CT_TDHF, CT_TEST
Expand Down Expand Up @@ -62,10 +50,6 @@ enum PotentialType {
POT_singles_
};

/// Assigns strings to enums for formated output
std::string
assign_name(const PairFormat& input);

/// Assigns strings to enums for formated output
std::string
assign_name(const CCState& input);
Expand Down Expand Up @@ -744,66 +728,25 @@ struct CCConvolutionOperator {

CCConvolutionOperator(const CCConvolutionOperator& other) = default;

friend bool can_combine(const CCConvolutionOperator& left, const CCConvolutionOperator& right) {
return (combine_OT(left,right).first!=OT_UNDEFINED);
protected:

friend CCConvolutionOperator combine(const CCConvolutionOperator& a, const CCConvolutionOperator& b) {
auto info= SeparatedConvolution<double,3>::combine_OT((*a.get_op()),(*b.get_op()));
Parameters param;
param.gamma=info.mu;
param.thresh_op=info.thresh;
param.lo=info.lo;
param.freeze=a.parameters.freeze;
return CCConvolutionOperator(a.world, info.type, param);
}

friend std::pair<OpType,Parameters> combine_OT(const CCConvolutionOperator& left, const CCConvolutionOperator& right) {
OpType type=OT_UNDEFINED;
Parameters param=left.parameters;
if ((left.type()==OT_F12) and (right.type()==OT_G12)) {
type=OT_FG12;
}
if ((left.type()==OT_G12) and (right.type()==OT_F12)) {
type=OT_FG12;
param.gamma=right.parameters.gamma;
}
if ((left.type()==OT_F12) and (right.type()==OT_F12)) {
type=OT_F212;
// keep the original gamma
// (f12)^2 = (1- slater12)^2 = 1/(4 gamma) (1 - 2 exp(-gamma) + exp(-2 gamma))
MADNESS_CHECK(right.parameters.gamma == left.parameters.gamma);
}
return std::make_pair(type,param);
}


/// combine 2 convolution operators to one

/// @return a vector of pairs: factor and convolution operator
friend std::vector<std::pair<double,CCConvolutionOperator>> combine(const CCConvolutionOperator& left, const CCConvolutionOperator& right) {
MADNESS_CHECK(can_combine(left,right));
MADNESS_CHECK(left.world.id()==right.world.id());
auto [type,param]=combine_OT(left,right);
std::vector<std::pair<double,CCConvolutionOperator>> result;
if (type==OT_FG12) {
// fg = (1 - exp(-gamma r12)) / r12 = 1/r12 - exp(-gamma r12)/r12 = coulomb - bsh

// coulombfit return 1/r
// we need 1/(2 gamma) 1/r
result.push_back(std::make_pair(1.0/(2.0*param.gamma),CCConvolutionOperator(left.world, OT_G12, param)));

// bshfit returns 1/(4 pi) exp(-gamma r)/r
// we need 1/(2 gamma) exp(-gamma r)/r
const double factor = 4.0 * constants::pi /(2.0*param.gamma);
result.push_back(std::make_pair(-factor,CCConvolutionOperator(left.world, OT_BSH, param)));
} else if (type==OT_F212) {
// we use the slater operator which is S = e^(-y*r12), y=gamma
// the f12 operator is: 1/2y*(1-e^(-y*r12)) = 1/2y*(1-S)
// so the squared f12 operator is: f*f = 1/(4*y*y)(1-2S+S*S), S*S = S(2y) = e(-2y*r12)
// we have then: <xy|f*f|xy> = 1/(4*y*y)*(<xy|xy> - 2*<xy|S|xy> + <xy|SS|xy>)
// we have then: <xy|f*f|xy> =(<xy|f12|xy> - 1/(4*y*y)*2*<xy|S|xy>
MADNESS_CHECK(left.parameters.gamma==right.parameters.gamma);
const double prefactor = 1.0 / (4.0 * param.gamma); // Slater has no 1/(2 gamma) per se.
Parameters param2=param;
param2.gamma*=2.0;
result.push_back(std::make_pair(1.0*prefactor,CCConvolutionOperator(left.world, OT_ONE, param)));
result.push_back(std::make_pair(-2.0*prefactor,CCConvolutionOperator(left.world, OT_SLATER, left.parameters)));
result.push_back(std::make_pair(1.0*prefactor,CCConvolutionOperator(left.world, OT_SLATER, param2)));
}
return result;
friend std::shared_ptr<CCConvolutionOperator> combine(const std::shared_ptr<CCConvolutionOperator>& a,
const std::shared_ptr<CCConvolutionOperator>& b) {
auto bla=combine(*a,*b);
return std::shared_ptr<CCConvolutionOperator>(new CCConvolutionOperator(combine(*a,*b)));
}

public:
/// @param[in] f: a 3D function
/// @param[out] the convolution op(f), no intermediates are used
real_function_3d operator()(const real_function_3d& f) const {
Expand Down Expand Up @@ -834,8 +777,8 @@ struct CCConvolutionOperator {
// @param[in] f: a vector of 3D functions
// @param[out] the convolution of op with each function, no intermeditates are used
vector_real_function_3d operator()(const vector_real_function_3d& f) const {
if (op) return apply<double, double, 3>(world, (*op), f);
return f;
MADNESS_CHECK(op);
return apply<double, double, 3>(world, (*op), f);
}

// @param[in] bra: a 3D CC_function, if nuclear-correlation factors are used they have to be applied before
Expand Down Expand Up @@ -891,9 +834,10 @@ struct CCConvolutionOperator {

/// create a TwoElectronFactory with the operatorkernel
TwoElectronFactory get_kernel() const {
if (type() == OT_G12) return TwoElectronFactory(world).dcut(1.e-7);
else if (type() == OT_F12) return TwoElectronFactory(world).dcut(1.e-7).f12().gamma(parameters.gamma);
else if (type() == OT_FG12) return TwoElectronFactory(world).dcut(1.e-7).BSH().gamma(parameters.gamma);
if (type() == OpType::OT_G12) return TwoElectronFactory(world).dcut(1.e-7);
else if (type() == OpType::OT_F12) return TwoElectronFactory(world).dcut(1.e-7).f12().gamma(parameters.gamma);
else if (type() == OpType::OT_FG12) return TwoElectronFactory(world).dcut(1.e-7).BSH().gamma(parameters.gamma);
else if (type() == OpType::OT_SLATER) return TwoElectronFactory(world).dcut(1.e-7).slater().gamma(parameters.gamma);
else error("no kernel of type " + name() + " implemented");
return TwoElectronFactory(world);
}
Expand All @@ -908,7 +852,7 @@ struct CCConvolutionOperator {
/// the world
World& world;
/// the operatortype, currently this can be g12_ or f12_
const OpType operator_type = OT_UNDEFINED;
const OpType operator_type = OpType::OT_UNDEFINED;

/// @param[in] optype: can be f12_ or g12_ depending on which operator shall be intitialzied
/// @param[in] parameters: parameters (thresholds etc)
Expand All @@ -928,8 +872,16 @@ struct CCConvolutionOperator {
<< "!!!!!\n\n" << std::endl;
MADNESS_EXCEPTION(msg.c_str(), 1);
}
public:
};

template<typename T=double>
std::shared_ptr<CCConvolutionOperator> CCConvolutionOperatorPtr(World& world, const OpType type,
CCConvolutionOperator::Parameters param) {
return std::shared_ptr<CCConvolutionOperator>(new CCConvolutionOperator(world,type,param));
}


class CCPair : public archive::ParallelSerializableObject {
public:
CCPair(){};
Expand Down
5 changes: 3 additions & 2 deletions src/madness/chem/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ set(MADCHEM_HEADERS
gth_pseudopotential.h
GuessFactory.h
localizer.h
lowrankfunction.h
masks_and_boxes.h
molecularbasis.h
molecular_functors.h
Expand Down Expand Up @@ -139,8 +140,8 @@ if(BUILD_TESTING)
SET(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
# The list of unit test source files
set(CHEM_TEST_SOURCES_SHORT test_pointgroupsymmetry.cc test_masks_and_boxes.cc
test_qc.cc test_MolecularOrbitals.cc test_BSHApply.cc)
set(CHEM_TEST_SOURCES_LONG test_localizer.cc test_ccpairfunction.cc)
test_qc.cc test_MolecularOrbitals.cc test_BSHApply.cc)
set(CHEM_TEST_SOURCES_LONG test_localizer.cc test_ccpairfunction.cc test_low_rank_function.cc)
if (LIBXC_FOUND)
list(APPEND CHEM_TEST_SOURCES_SHORT test_dft.cc )
list(APPEND CHEM_TEST_SOURCES_LONG test_SCFOperators.cc)
Expand Down
2 changes: 0 additions & 2 deletions src/madness/chem/PNOF12Potentials.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,8 +81,6 @@ F12Potentials::F12Potentials(World& world,const Nemo& nemo, const BasisFunctions
fop = std::shared_ptr < real_convolution_3d > (SlaterF12OperatorPtr(world, param.gamma(), lo, eps));
slaterop = std::shared_ptr < real_convolution_3d > (SlaterF12OperatorPtr(world, param.gamma(), lo, eps));
slaterop_sq = std::shared_ptr < real_convolution_3d > (SlaterF12OperatorPtr(world, param.gamma() * 2.0, lo, eps));
slaterop->is_slaterf12 = false;
slaterop_sq->is_slaterf12 = false;
fop = std::shared_ptr < real_convolution_3d > (SlaterF12OperatorPtr(world, param.gamma(), lo, eps));

//test
Expand Down
2 changes: 1 addition & 1 deletion src/madness/chem/TDHF.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ void TDHF::initialize() {

msg.section("Initialize TDHF Class");
msg.debug = parameters.debug();
g12=std::make_shared<CCConvolutionOperator>(world, OT_G12, parameters.get_ccc_parameters(get_calcparam().lo()));
g12=std::make_shared<CCConvolutionOperator>(world, OpType::OT_G12, parameters.get_ccc_parameters(get_calcparam().lo()));

const double old_thresh = FunctionDefaults<3>::get_thresh();
if (old_thresh > parameters.thresh() * 0.1 and old_thresh > 1.e-5) {
Expand Down
Loading

0 comments on commit 1c3b19d

Please sign in to comment.