diff --git a/src/apps/mp2/mp2.cc b/src/apps/mp2/mp2.cc index b15de2ba34c..154919f1f76 100644 --- a/src/apps/mp2/mp2.cc +++ b/src/apps/mp2/mp2.cc @@ -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) { diff --git a/src/madness/chem/CCPotentials.cc b/src/madness/chem/CCPotentials.cc index b5de6c1b837..0bb0d150d1b 100644 --- a/src/madness/chem/CCPotentials.cc +++ b/src/madness/chem/CCPotentials.cc @@ -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); diff --git a/src/madness/chem/CCStructures.cc b/src/madness/chem/CCStructures.cc index 7181a6b8acf..44cbed09443 100644 --- a/src/madness/chem/CCStructures.cc +++ b/src/madness/chem/CCStructures.cc @@ -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); @@ -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); @@ -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; @@ -408,60 +406,8 @@ SeparatedConvolution * 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(world,OperatorInfo(parameters.gamma,parameters.lo,parameters.thresh_op,type)); } /// Assigns strings to enums for formated output @@ -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) { diff --git a/src/madness/chem/CCStructures.h b/src/madness/chem/CCStructures.h index d214a2f996d..4210e491313 100644 --- a/src/madness/chem/CCStructures.h +++ b/src/madness/chem/CCStructures.h @@ -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); @@ -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; @@ -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 will be deleted @@ -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(" intermediate"); else if (type == PARTICLE) for (const auto& tmp:imP.allpairs) - tmp.second.print_size(" intermediate"); else if (type == RESPONSE) for (const auto& tmp:imR.allpairs) - tmp.second.print_size(" intermediate"); } @@ -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; @@ -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) @@ -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); } diff --git a/src/madness/mra/operator.h b/src/madness/mra/operator.h index b76e12a7f07..91f0573b177 100644 --- a/src/madness/mra/operator.h +++ b/src/madness/mra/operator.h @@ -139,10 +139,41 @@ namespace madness { OT_BSH /// exp(-r)/r }; + /// operator type to string + template // 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 -// OperatorInfo(double mu, double lo=1.e-5, double thresh=FunctionDefaults::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