Molecules in high pressure environments

python/mrchem/helpers.py

@@ -126,6 +126,16 @@ def write_scf_fock(user_dict, wf_dict, origin):
             "r_O": origin,
         }
 
+    # Confinement Potential
+    if len(user_dict["ConfinementPotential"]) > 0:
+        fock_dict["confinement_operator"] = {
+            "r_0": user_dict["ConfinementPotential"]["radius"],
+            "N": user_dict["ConfinementPotential"]["stiffness"],
+            "slope": user_dict["ConfinementPotential"]["slope"],
+            "cavity_radii": user_dict["ConfinementPotential"]["cavity_radii"],
+            #"centers": user_dict["ConfinementPotential"]["centers"],
+        }
+
     return fock_dict
 
 
@@ -184,6 +194,7 @@ def write_scf_guess(user_dict, wf_dict):
         "relativity": wf_dict["relativity_name"],
         "environment": wf_dict["environment_name"],
         "external_field": wf_dict["external_name"],
+        "confinement_potential": wf_dict["confinement_name"],
         "screen": scf_dict["guess_screen"],
         "localize": scf_dict["localize"],
         "restricted": user_dict["WaveFunction"]["restricted"],
@@ -217,6 +228,7 @@ def write_scf_solver(user_dict, wf_dict):
         "relativity": wf_dict["relativity_name"],
         "environment": wf_dict["environment_name"],
         "external_field": wf_dict["external_name"],
+        "confinement_potential": wf_dict["confinement_name"],
         "kain": scf_dict["kain"],
         "max_iter": scf_dict["max_iter"],
         "rotation": scf_dict["rotation"],
@@ -519,10 +531,16 @@ def parse_wf_method(user_dict):
         # Labels to aggregate
         external_name = f"Electric field ({x}, {y}, {z})"
 
+    # Determine confinement potential name
+    confinement_name = "None"
+    if len(user_dict["ConfinementPotential"]) > 0:
+        confinement_name = "Confinement"
+
     wf_dict = {
         "relativity_name": relativity_name,
         "environment_name": environment_name,
         "external_name": external_name,
+        "confinement_name": confinement_name,
         "method_name": method_name,
         "method_type": method_type,
         "dft_funcs": dft_funcs,

python/mrchem/input_parser/api.py

@@ -365,6 +365,20 @@ def stencil() -> JSONDict:
                                                               '== 3'],
                                             'type': 'List[float]'}],
                         'name': 'ExternalFields'},
+                    {   'keywords': [   {   'default': 10.0,
+                                            'name': 'radius',
+                                            'type': 'float',},
+                                        {   'default': 2,
+                                            'name': 'stiffness',
+                                            'type': 'int'},
+                                        {   'default': 0.2,
+                                            'name': 'slope',
+                                            'type': 'float'},
+                                        {   'name': 'cavity_radii',
+                                            'type': 'List[float]'}],
+#                                        {   'name': 'centers',
+#                                            'type': 'List[List[float]]'}],
+                        'name': 'ConfinementPotential'},
                     {   'keywords': [   {   'default': [0.0],
                                             'name': 'frequency',
                                             'type': 'List[float]'}],

src/analyticfunctions/CMakeLists.txt

@@ -3,4 +3,5 @@ target_sources(mrchem PRIVATE
     ${CMAKE_CURRENT_SOURCE_DIR}/NuclearFunction.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/NuclearGradientFunction.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/CUBEfunction.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/ConfinementFunction.cpp
     )

src/analyticfunctions/ConfinementFunction.cpp

@@ -0,0 +1,42 @@
+#include "ConfinementFunction.h"
+#include "environment/Cavity.h"
+#include "utils/math_utils.h"
+
+namespace mrchem {
+
+  ConfinementFunction::ConfinementFunction(double r_0, const int N, double s, std::vector<double> R, std::vector<mrcpp::Coord<3>> centers)
+    : radius(r_0)
+    , slope(s)
+    , stiffness(N)
+    , cavity_radii(R)
+    , centers(centers) {}
+
+  double ConfinementFunction::evalf(const mrcpp::Coord<3> &r) const {
+    double f = 0.0;
+
+    const int N = this->stiffness;
+    auto r_0 = this->radius;
+    auto s = this->slope;
+    auto R = this->cavity_radii;
+    auto coords = this->centers;
+
+    Cavity sphere(coords, R, s);
+
+    //double min_dist = 1000.0;
+    //double &rad = min_dist;
+    //for (auto& center : coords) {
+    //  double distance = math_utils::calc_distance(r, center);
+    // if (distance < min_dist) {
+    //rad = distance;
+    //  }
+    //}
+
+    //f += std::pow(min_dist / r_0, N);
+    mrcpp::Coord<3> origin = {0.0, 0.0, 0.0};
+    f += std::pow(math_utils::calc_distance(origin, r) / r_0, N);
+    f *= (1 - sphere.evalf(r));
+
+    return f;
+  }
+
+} // namespace mrchem

src/analyticfunctions/ConfinementFunction.h

@@ -0,0 +1,29 @@
+#pragma once
+
+#include <MRCPP/MWFunctions>
+
+namespace mrchem {
+
+  class ConfinementFunction : public mrcpp::RepresentableFunction<3> {
+  public:
+    ConfinementFunction(double r_0, const int N, double slope, std::vector<double> R, std::vector<mrcpp::Coord<3>> centers);
+
+    auto getRadius() { return this->radius; }
+    auto getStiffness() { return this->stiffness; }
+    auto getSlope() { return this->slope; }
+    auto getCavityradii() { return this->cavity_radii; }
+    auto getCenters() { return this->centers; }
+
+  protected:
+
+    double radius; //
+    const int stiffness; // stffness parameter
+    double slope; // cavity boundary slope
+    std::vector<double> cavity_radii;
+    std::vector<mrcpp::Coord<3>> centers;
+
+    double evalf(const mrcpp::Coord<3> &r) const override;
+
+  };
+
+} // namespace mrchem

src/driver.cpp

@@ -54,6 +54,7 @@
 #include "qmoperators/one_electron/NuclearGradientOperator.h"
 #include "qmoperators/one_electron/NuclearOperator.h"
 #include "qmoperators/one_electron/ZoraOperator.h"
+#include "qmoperators/one_electron/ConfinementOperator.h"
 
 #include "qmoperators/one_electron/H_BB_dia.h"
 #include "qmoperators/one_electron/H_BM_dia.h"
@@ -268,6 +269,7 @@ json driver::scf::run(const json &json_scf, Molecule &mol) {
         auto relativity = json_scf["scf_solver"]["relativity"];
         auto environment = json_scf["scf_solver"]["environment"];
         auto external_field = json_scf["scf_solver"]["external_field"];
+	auto confinement_potential = json_scf["scf_solver"]["confinement_potential"];
         auto max_iter = json_scf["scf_solver"]["max_iter"];
         auto rotation = json_scf["scf_solver"]["rotation"];
         auto localize = json_scf["scf_solver"]["localize"];
@@ -287,6 +289,7 @@ json driver::scf::run(const json &json_scf, Molecule &mol) {
         solver.setRelativityName(relativity);
         solver.setEnvironmentName(environment);
         solver.setExternalFieldName(external_field);
+	solver.setConfinementPotentialName(confinement_potential);
         solver.setCheckpoint(checkpoint);
         solver.setCheckpointFile(file_chk);
         solver.setMaxIterations(max_iter);
@@ -401,16 +404,18 @@ bool driver::scf::guess_energy(const json &json_guess, Molecule &mol, FockBuilde
     auto relativity = json_guess["relativity"];
     auto environment = json_guess["environment"];
     auto external_field = json_guess["external_field"];
+    auto confinement_potential = json_guess["confinement_potential"];
     auto localize = json_guess["localize"];
 
     mrcpp::print::separator(0, '~');
-    print_utils::text(0, "Calculation    ", "Compute initial energy");
-    print_utils::text(0, "Method         ", method);
-    print_utils::text(0, "Relativity     ", relativity);
-    print_utils::text(0, "Environment    ", environment);
-    print_utils::text(0, "External fields", external_field);
-    print_utils::text(0, "Precision      ", print_utils::dbl_to_str(prec, 5, true));
-    print_utils::text(0, "Localization   ", (localize) ? "On" : "Off");
+    print_utils::text(0, "Calculation          ", "Compute initial energy");
+    print_utils::text(0, "Method               ", method);
+    print_utils::text(0, "Relativity           ", relativity);
+    print_utils::text(0, "Environment          ", environment);
+    print_utils::text(0, "External fields      ", external_field);
+    print_utils::text(0, "Confinement potential", confinement_potential);
+    print_utils::text(0, "Precision            ", print_utils::dbl_to_str(prec, 5, true));
+    print_utils::text(0, "Localization         ", (localize) ? "On" : "Off");
     mrcpp::print::separator(0, '~', 2);
 
     Timer t_scf;
@@ -1123,6 +1128,23 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
         auto V_ext = std::make_shared<ElectricFieldOperator>(field, r_O);
         F.getExtOperator() = V_ext;
     }
+    ///////////////////////////////////////////////////////////
+    ///////////////   Confinement Potential   /////////////////
+    ///////////////////////////////////////////////////////////
+    if (json_fock.contains("confinement_operator")) {
+        auto r_0 = json_fock["confinement_operator"]["r_0"];
+	auto N = json_fock["confinement_operator"]["N"];
+	auto s = json_fock["confinement_operator"]["slope"];
+	auto R = json_fock["confinement_operator"]["cavity_radii"];
+	std::vector<mrcpp::Coord<3>> centers;
+	for (auto i = 0; i < mol.getNNuclei(); i++) {
+	    const auto &nuc = mol.getNuclei()[i];
+	    const auto &r_i = nuc.getCoord();
+	    centers.push_back(r_i);
+       }
+	auto V_N = std::make_shared<ConfinementOperator>(r_0, N, s, R, centers);
+	F.getConfinementOperator() = V_N;
+    }
     F.build(exx);
 }
 

src/properties/SCFEnergy.h

@@ -47,11 +47,12 @@ class SCFEnergy final {
                        double en = 0.0, double ee = 0.0,
                        double x = 0.0, double xc = 0.0,
                        double next = 0.0, double eext = 0.0,
-                       double rt = 0.0, double rn = 0.0, double re = 0.0) :
+                       double rt = 0.0, double rn = 0.0,
+		       double re = 0.0, double cp = 0.0) :
         E_kin(kin), E_nn(nn), E_en(en), E_ee(ee),
-          E_x(x), E_xc(xc), E_next(next), E_eext(eext), Er_tot(rt), Er_nuc(rn), Er_el(re) {
+	E_x(x), E_xc(xc), E_next(next), E_eext(eext), Er_tot(rt), Er_nuc(rn), Er_el(re), E_cp(cp) {
             E_nuc = E_nn + E_next + Er_nuc;
-            E_el = E_kin + E_en + E_ee + E_xc + E_x + E_eext + Er_el;
+            E_el = E_kin + E_en + E_ee + E_xc + E_x + E_eext + Er_el, E_cp;
         }
 
     double getTotalEnergy() const { return this->E_nuc + this->E_el; }
@@ -69,6 +70,7 @@ class SCFEnergy final {
     double getReactionEnergy() const { return this->Er_tot; }
     double getElectronReactionEnergy() const { return this->Er_el; }
     double getNuclearReactionEnergy() const { return this->Er_nuc; }
+    double getConfinementPotentialEnergy() const { return this-> E_cp; }
 
     void print(const std::string &id) const {
         auto E_au = E_nuc + E_el;
@@ -78,6 +80,7 @@ class SCFEnergy final {
 
         bool has_ext = (std::abs(E_eext) > mrcpp::MachineZero) || (std::abs(E_next) > mrcpp::MachineZero);
         bool has_react = (std::abs(Er_el) > mrcpp::MachineZero) || (std::abs(Er_nuc) > mrcpp::MachineZero);
+	bool has_conf = (std::abs(E_cp) > mrcpp::MachineZero);
 
         auto pprec = 2 * mrcpp::Printer::getPrecision();
         mrcpp::print::header(0, "Molecular Energy (" + id + ")");
@@ -99,6 +102,10 @@ class SCFEnergy final {
             print_utils::scalar(0, "Reaction energy (nuc) ", Er_nuc,  "(au)", pprec, false);
             print_utils::scalar(0, "Reaction energy (tot) ", Er_tot,  "(au)", pprec, false);
         }
+	if (has_conf) {
+	    mrcpp::print::separator(0, '-');
+	    print_utils::scalar(0, "Confinement potential ", E_cp, "(au)", pprec, false);
+	}
         mrcpp::print::separator(0, '-');
         print_utils::scalar(0, "Electronic energy", E_el,   "(au)", pprec, false);
         print_utils::scalar(0, "Nuclear energy   ", E_nuc,  "(au)", pprec, false);
@@ -124,6 +131,7 @@ class SCFEnergy final {
             {"Er_tot", Er_tot},
             {"Er_el", Er_el},
             {"Er_nuc", Er_nuc},
+	    {"E_cp", E_cp},
             {"E_nuc", E_nuc},
             {"E_tot", E_nuc + E_el}
         };
@@ -144,6 +152,7 @@ class SCFEnergy final {
     double Er_tot{0.0};
     double Er_nuc{0.0};
     double Er_el{0.0};
+    double E_cp{0.0};
 };
 // clang-format on
 

src/qmoperators/one_electron/CMakeLists.txt

@@ -1,4 +1,5 @@
 target_sources(mrchem PRIVATE
     ${CMAKE_CURRENT_SOURCE_DIR}/NuclearOperator.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/ZoraOperator.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/ConfinementPotential.cpp
     )

src/qmoperators/one_electron/ConfinementOperator.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include "tensor/RankZeroOperator.h"
+#include "ConfinementPotential.h"
+#include "chemistry/Molecule.h"
+
+/** @class ConfinementOperator
+ * 
+ * @brief Operator containing a ConfinementPotential
+ *
+ * A TensorOperator realization of @class ConfinementPotential.
+ *
+ */
+
+namespace mrchem {
+
+  class ConfinementOperator final : public RankZeroOperator {
+  public:
+    ConfinementOperator(double r_0, const int N, double s, std::vector<double> R, std::vector<mrcpp::Coord<3>> centers) {
+      potential = std::make_shared<ConfinementPotential>(r_0, N, s, R, centers);
+
+      // Invoke operator= to assign *this operator
+      RankZeroOperator &Co = (*this);
+      Co = potential;
+      Co.name() = "Co";
+    }
+    ~ConfinementOperator() override = default;
+
+    auto getRadius() { return this->potential->getRadius(); }
+    auto getStiffness() { return this->potential->getStiffness(); }
+
+  private:
+    std::shared_ptr<ConfinementPotential> potential{nullptr};
+  };
+
+} // namespace mrchem

src/qmoperators/one_electron/ConfinementPotential.cpp

@@ -0,0 +1,41 @@
+#include "ConfinementPotential.h"
+#include "analyticfunctions/ConfinementFunction.h"
+
+#include "utils/math_utils.h"
+#include "chemistry/Molecule.h"
+
+namespace mrchem {
+
+  ConfinementPotential::ConfinementPotential(double r_0, const int N, double s, std::vector<double> R, std::vector<mrcpp::Coord<3>> coords)
+    : QMPotential(1, false)
+    , radius(r_0)
+    , stiffness(N)
+    , slope(s)
+    , cavity_radii(R)
+    , centers(coords) {}
+
+
+  void ConfinementPotential::setup(double prec) {
+    // Initalize the function representing the confinement
+    ConfinementFunction *f_loc = nullptr;
+
+    f_loc = new ConfinementFunction(this->getRadius(), this->getStiffness(), this->getSlope(), this->getCavityradii(), this->getCenters());
+
+    // Project the potential onto the function representation
+    mrcpp::ComplexFunction V_loc(false);
+    mrcpp::cplxfunc::project(V_loc, *f_loc, NUMBER::Real, prec);
+
+    mrcpp::ComplexFunction &V_c = (*this);
+    setApplyPrec(prec);
+    V_c = V_loc;
+    delete f_loc;
+
+  }
+
+  void ConfinementPotential::clear() {
+    mrcpp::ComplexFunction::free(NUMBER::Total);
+    clearApplyPrec();
+  }
+
+
+} // namespace mrchem