Skip to content

Commit

Permalink
making mp2 work again
Browse files Browse the repository at this point in the history
  • Loading branch information
@fbischoff
fbischoff committed Oct 27, 2023
1 parent 9ff2820 commit ed21aa1
Show file tree
Hide file tree
Showing 5 changed files with 53 additions and 110 deletions.
1 change: 1 addition & 0 deletions src/apps/mp2/mp2.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,7 @@ int main(int argc, char** argv) {
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) {
Expand Down
1 change: 0 additions & 1 deletion src/madness/chem/CCPotentials.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3420,7 +3420,6 @@ void CCPotentials::test() {
CCState test5 = GROUND_STATE;
PotentialType test6 = POT_F3D_;
assign_name(test2);
assign_name(test3);
assign_name(test4);
assign_name(test5);
assign_name(test6);
Expand Down
101 changes: 8 additions & 93 deletions src/madness/chem/CCStructures.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,7 +304,7 @@ CCConvolutionOperator::operator()(const CCFunction& bra, const CCFunction& ket,
real_function_3d result;
if (not use_im) {
if (world.rank() == 0)
std::cout << "Recalculating <" << bra.name() << "|" << assign_name(operator_type) << "|" << ket.name()
std::cout << "Recalculating <" << bra.name() << "|" << name() << "|" << ket.name()
<< ">\n";
result = ((*op)(bra.function * ket.function)).truncate();
} else if (bra.type == HOLE and ket.type == HOLE and not imH.allpairs.empty()) result = imH(bra.i, ket.i);
Expand All @@ -321,18 +321,16 @@ 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 == OpType::OT_G12);
MADNESS_ASSERT(op);
MADNESS_CHECK(particle == 1 or particle == 2);
MADNESS_CHECK(op);
op->particle() = particle;
return (*op)(u);
}

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 == OpType::OT_G12);
MADNESS_ASSERT(op);
MADNESS_CHECK(particle == 1 or particle == 2);
MADNESS_CHECK(op);
const real_function_6d tmp = multiply(copy(u), copy(bra.function), particle);
op->particle() = particle;
const real_function_6d g_tmp = (*op)(tmp);
Expand All @@ -348,7 +346,7 @@ void CCConvolutionOperator::update_elements(const CC_vecfunction& bra, const CC_
<< std::endl;
if (bra.type != HOLE)
error("Can not create intermediate of type " + operation_name + " , bra-element has to be of type HOLE");
op.reset(init_op(operator_type, parameters));
op.reset(init_op(type(), parameters));
intermediateT xim;
for (auto tmpk : bra.functions) {
const CCFunction& k = tmpk.second;
Expand Down Expand Up @@ -408,60 +406,8 @@ 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 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 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 OpType::OT_F212 : {
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 SlaterF12sqOperatorPtr(world, parameters.gamma, parameters.lo, parameters.thresh_op);
}
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 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 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;
}
default : {
error("Unknown operatorype " + assign_name(type));
MADNESS_EXCEPTION("error", 1);
}
}

print("init_op: creating",type,"with thresh, lo, gamma",parameters.thresh_op,parameters.gamma,parameters.lo);
return new SeparatedConvolution<double,3>(world,OperatorInfo(parameters.gamma,parameters.lo,parameters.thresh_op,type));
}

/// Assigns strings to enums for formated output
Expand All @@ -481,37 +427,6 @@ assign_name(const CCState& input) {
return "unknown pairtype";
}

/// Assigns strings to enums for formated output
std::string
assign_name(const OpType& input) {
switch (input) {
case OpType::OT_G12:
return "g12";
case OpType::OT_F12:
return "f12";
case OpType::OT_SLATER:
return "slater";
case OpType::OT_FG12:
return "fg12";
case OpType::OT_BSH:
return "bsh";
case OpType::OT_ONE:
return "identity";
case OpType::OT_GAUSS: /// exp(-r2)
return "gauss";
case OpType::OT_F212: /// (1-exp(-r))^2
return "f12^2";
case OpType::OT_F2G12: /// (1-exp(-r))^2/r = 1/r + exp(-2r)/r - 2 exp(-r)/r
return "f12^2g";
default: {
MADNESS_EXCEPTION("Unvalid enum assignement!", 1);
return "undefined";
}
}
MADNESS_EXCEPTION("assign_name:optype, should not end up here", 1);
return "unknown operatortype";
}

/// Assigns enum to string
CalcType
assign_calctype(const std::string name) {
Expand Down
25 changes: 11 additions & 14 deletions src/madness/chem/CCStructures.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,10 +54,6 @@ enum PotentialType {
std::string
assign_name(const CCState& input);

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

/// Assigns enum to string
CalcType
assign_calctype(const std::string name);
Expand Down Expand Up @@ -722,8 +718,7 @@ struct CCConvolutionOperator {
/// @param[in] optype: the operatortype (can be g12_ or f12_)
/// @param[in] param: the parameters of the current CC-Calculation (including function and operator thresholds and the exponent for f12)
CCConvolutionOperator(World& world, const OpType type, Parameters param) : parameters(param), world(world),
operator_type(type),
op(init_op(operator_type, parameters)) {
op(init_op(type, parameters)) {
}

CCConvolutionOperator(const CCConvolutionOperator& other) = default;
Expand Down Expand Up @@ -806,7 +801,11 @@ struct CCConvolutionOperator {
void update_elements(const CC_vecfunction& bra, const CC_vecfunction& ket);

/// @param[out] prints the name of the operator (convenience) which is g12 or f12 or maybe other things like gf in the future
std::string name() const { return assign_name(operator_type); }
std::string name() const {
std::stringstream ss;
ss << type();
return ss.str();
}

/// @param[in] the type of which intermediates will be deleted
/// e.g if(type==HOLE) then all intermediates of type <mo_bra_k|op|HOLE> will be deleted
Expand All @@ -822,15 +821,15 @@ struct CCConvolutionOperator {
void print_intermediate(const FuncType type) const {
if (type == HOLE)
for (const auto& tmp:imH.allpairs)
tmp.second.print_size("<H" + std::to_string(tmp.first.first) + "|" + assign_name(operator_type) + "|H" +
tmp.second.print_size("<H" + std::to_string(tmp.first.first) + "|" + name() + "|H" +
std::to_string(tmp.first.second) + "> intermediate");
else if (type == PARTICLE)
for (const auto& tmp:imP.allpairs)
tmp.second.print_size("<H" + std::to_string(tmp.first.first) + "|" + assign_name(operator_type) + "|P" +
tmp.second.print_size("<H" + std::to_string(tmp.first.first) + "|" + name() + "|P" +
std::to_string(tmp.first.second) + "> intermediate");
else if (type == RESPONSE)
for (const auto& tmp:imR.allpairs)
tmp.second.print_size("<H" + std::to_string(tmp.first.first) + "|" + assign_name(operator_type) + "|R" +
tmp.second.print_size("<H" + std::to_string(tmp.first.first) + "|" + name() + "|R" +
std::to_string(tmp.first.second) + "> intermediate");
}

Expand All @@ -844,7 +843,7 @@ struct CCConvolutionOperator {
return TwoElectronFactory(world);
}

OpType type() const { return operator_type; }
OpType type() const { return get_op()->info.type; }

const Parameters parameters;

Expand All @@ -853,8 +852,6 @@ struct CCConvolutionOperator {
private:
/// the world
World& world;
/// the operatortype, currently this can be g12_ or f12_
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 @@ -870,7 +867,7 @@ struct CCConvolutionOperator {
/// the function will throw an MADNESS_EXCEPTION
void error(const std::string& msg) const {
if (world.rank() == 0)
std::cout << "\n\n!!!!ERROR in CCConvolutionOperator " << assign_name(operator_type) << ": " << msg
std::cout << "\n\n!!!!ERROR in CCConvolutionOperator " << name() << ": " << msg
<< "!!!!!\n\n" << std::endl;
MADNESS_EXCEPTION(msg.c_str(), 1);
}
Expand Down
35 changes: 33 additions & 2 deletions src/madness/mra/operator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,10 +139,41 @@ namespace madness {
OT_BSH /// exp(-r)/r
};

/// operator type to string
template<std::size_t N=1> // dummy template argument to avoid ambiguity with the other operator<<
std::ostream& operator<<(std::ostream& os, const OpType type) {
auto name = [](OpType type) {
switch (type) {
case OpType::OT_UNDEFINED:
return "undefined";
case OpType::OT_ONE:
return "identity";
case OpType::OT_G12:
return "g12";
case OpType::OT_SLATER:
return "slater";
case OpType::OT_GAUSS:
return "gauss";
case OpType::OT_F12:
return "f12";
case OpType::OT_FG12:
return "fg12";
case OpType::OT_F212:
return "f12^2";
case OpType::OT_F2G12:
return "f12^2g";
case OpType::OT_BSH:
return "bsh";
default:
return "undefined";
}
};
os << name(type);
return os;
}

struct OperatorInfo {
OperatorInfo() = default;
// template<std::size_t NDIM>
// OperatorInfo(double mu, double lo=1.e-5, double thresh=FunctionDefaults<NDIM>::get_thresh()) : mu(mu), thresh(thresh), lo(lo) {}
OperatorInfo(double mu, double lo, double thresh, OpType type) : mu(mu), thresh(thresh), lo(lo), type(type) {}
OpType type=OT_UNDEFINED; ///< introspection
double mu=0.0; ///< some introspection
Expand Down

0 comments on commit ed21aa1

Please sign in to comment.