diff --git a/cmake/compiler_flags/CXXFlags.cmake b/cmake/compiler_flags/CXXFlags.cmake index a12df3d12..cbcf32898 100644 --- a/cmake/compiler_flags/CXXFlags.cmake +++ b/cmake/compiler_flags/CXXFlags.cmake @@ -13,7 +13,7 @@ option(ENABLE_ARCH_FLAGS "Enable architecture-specific compiler flags" ON) -set(CMAKE_CXX_STANDARD 14) +set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED TRUE) set(CMAKE_CXX_EXTENSIONS FALSE) set(CMAKE_EXPORT_COMPILE_COMMANDS TRUE) diff --git a/doc/users/user_ref.rst b/doc/users/user_ref.rst index 6d6ddbcce..d87e2714d 100644 --- a/doc/users/user_ref.rst +++ b/doc/users/user_ref.rst @@ -905,13 +905,7 @@ User input reference **Default** ``1.0`` - :epsilon_out: Permittivity outside the cavity. This is characteristic of the solvent used. - - **Type** ``float`` - - **Default** ``1.0`` - - :formulation: Formulation of the Permittivity function. Currently only the exponential is used. + :formulation: Formulation of the Permittivity function. Currently only the exponential is available. **Type** ``str`` @@ -920,6 +914,28 @@ User input reference **Predicates** - ``value.lower() in ['exponential']`` + :red:`Sections` + :epsilon_out: Parameters for the continuum solvent outside the cavity. + + :red:`Keywords` + :nonequilibrium: Whether to use the nonequilibrium formulation of response, *i.e.* use the dynamic permittivity for the calculation of the response reaction field. Defaults to false. + + **Type** ``bool`` + + **Default** ``False`` + + :static: Static permittivity outside the cavity. This is characteristic of the solvent used. + + **Type** ``float`` + + **Default** ``1.0`` + + :dynamic: Dynamic permittivity outside the cavity. This is characteristic of the solvent used and relevant only in response calculations. Defaults to the same value as `epsilon_static`. + + **Type** ``float`` + + **Default** ``user['PCM']['Permittivity']['epsilon_out']['static']`` + :Constants: Physical and mathematical constants used by MRChem :red:`Keywords` diff --git a/python/mrchem/helpers.py b/python/mrchem/helpers.py index 506c4f52b..616c7dd3c 100644 --- a/python/mrchem/helpers.py +++ b/python/mrchem/helpers.py @@ -73,16 +73,7 @@ def write_scf_fock(user_dict, wf_dict, origin): # Reaction if user_dict["WaveFunction"]["environment"].lower() == "pcm": - fock_dict["reaction_operator"] = { - "poisson_prec": user_dict["world_prec"], - "kain": user_dict["PCM"]["SCRF"]["kain"], - "max_iter": user_dict["PCM"]["SCRF"]["max_iter"], - "dynamic_thrs": user_dict["PCM"]["SCRF"]["dynamic_thrs"], - "density_type": user_dict["PCM"]["SCRF"]["density_type"], - "epsilon_in": user_dict["PCM"]["Permittivity"]["epsilon_in"], - "epsilon_out": user_dict["PCM"]["Permittivity"]["epsilon_out"], - "formulation": user_dict["PCM"]["Permittivity"]["formulation"], - } + fock_dict["reaction_operator"] = _reaction_operator_handler(user_dict) # Coulomb if wf_dict["method_type"] in ["hartree", "hf", "dft"]: @@ -128,6 +119,32 @@ def write_scf_fock(user_dict, wf_dict, origin): return fock_dict +def _reaction_operator_handler(user_dict, rsp=False): + # convert density_type from string to integer + if user_dict["PCM"]["SCRF"]["density_type"] == "total": + density_type = 0 + elif user_dict["PCM"]["SCRF"]["density_type"] == "electronic": + density_type = 1 + else: + density_type = 2 + + return { + "poisson_prec": user_dict["world_prec"], + "kain": user_dict["PCM"]["SCRF"]["kain"], + "max_iter": user_dict["PCM"]["SCRF"]["max_iter"], + "dynamic_thrs": user_dict["PCM"]["SCRF"]["dynamic_thrs"], + # if doing a response calculation, then density_type is set to 1 (electronic only) + "density_type": 1 if rsp else density_type, + "epsilon_in": user_dict["PCM"]["Permittivity"]["epsilon_in"], + "epsilon_static": user_dict["PCM"]["Permittivity"]["epsilon_out"]["static"], + "epsilon_dynamic": user_dict["PCM"]["Permittivity"]["epsilon_out"]["dynamic"], + "nonequilibrium": user_dict["PCM"]["Permittivity"]["epsilon_out"][ + "nonequilibrium" + ], + "formulation": user_dict["PCM"]["Permittivity"]["formulation"], + } + + def write_scf_guess(user_dict, wf_dict): guess_str = user_dict["SCF"]["guess_type"].lower() guess_type = guess_str.split("_")[0] @@ -305,7 +322,6 @@ def write_rsp_calc(omega, user_dict, origin): rsp_calc["components"] = [] for dir in [0, 1, 2]: - rsp_comp = {} program_guess_type = user_guess_type @@ -405,6 +421,10 @@ def write_rsp_fock(user_dict, wf_dict): }, } + # Reaction + if user_dict["WaveFunction"]["environment"].lower() == "pcm": + fock_dict["reaction_operator"] = _reaction_operator_handler(user_dict, rsp=True) + return fock_dict diff --git a/python/mrchem/input_parser/README.md b/python/mrchem/input_parser/README.md index 4009aaee9..007a1627b 100644 --- a/python/mrchem/input_parser/README.md +++ b/python/mrchem/input_parser/README.md @@ -1,2 +1,2 @@ -This file was automatically generated by parselglossy on 2023-11-15 +This file was automatically generated by parselglossy on 2023-11-22 Editing is *STRONGLY DISCOURAGED* diff --git a/python/mrchem/input_parser/__init__.py b/python/mrchem/input_parser/__init__.py index fd5fdac6b..b2d8b86a3 100644 --- a/python/mrchem/input_parser/__init__.py +++ b/python/mrchem/input_parser/__init__.py @@ -1,4 +1,4 @@ # -*- coding: utf-8 -*- -# This file was automatically generated by parselglossy on 2023-11-15 +# This file was automatically generated by parselglossy on 2023-11-22 # Editing is *STRONGLY DISCOURAGED* diff --git a/python/mrchem/input_parser/api.py b/python/mrchem/input_parser/api.py index 16ca78ceb..742c7ec9e 100644 --- a/python/mrchem/input_parser/api.py +++ b/python/mrchem/input_parser/api.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- -# This file was automatically generated by parselglossy on 2023-11-15 +# This file was automatically generated by parselglossy on 2023-11-22 # Editing is *STRONGLY DISCOURAGED* from copy import deepcopy @@ -607,16 +607,23 @@ def stencil() -> JSONDict: { 'keywords': [ { 'default': 1.0, 'name': 'epsilon_in', 'type': 'float'}, - { 'default': 1.0, - 'name': 'epsilon_out', - 'type': 'float'}, { 'default': 'exponential', 'name': 'formulation', 'predicates': [ 'value.lower() ' 'in ' "['exponential']"], 'type': 'str'}], - 'name': 'Permittivity'}]}, + 'name': 'Permittivity', + 'sections': [ { 'keywords': [ { 'default': False, + 'name': 'nonequilibrium', + 'type': 'bool'}, + { 'default': 1.0, + 'name': 'static', + 'type': 'float'}, + { 'default': "user['PCM']['Permittivity']['epsilon_out']['static']", + 'name': 'dynamic', + 'type': 'float'}], + 'name': 'epsilon_out'}]}]}, { 'keywords': [ { 'default': 78.9451185, 'name': 'hartree2simagnetizability', 'type': 'float'}, diff --git a/python/mrchem/input_parser/cli.py b/python/mrchem/input_parser/cli.py index f32a3fba9..99e5688b6 100644 --- a/python/mrchem/input_parser/cli.py +++ b/python/mrchem/input_parser/cli.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- -# This file was automatically generated by parselglossy on 2023-11-15 +# This file was automatically generated by parselglossy on 2023-11-22 # Editing is *STRONGLY DISCOURAGED* import argparse diff --git a/python/mrchem/input_parser/docs/user_ref.rst b/python/mrchem/input_parser/docs/user_ref.rst index 6d6ddbcce..d87e2714d 100644 --- a/python/mrchem/input_parser/docs/user_ref.rst +++ b/python/mrchem/input_parser/docs/user_ref.rst @@ -905,13 +905,7 @@ User input reference **Default** ``1.0`` - :epsilon_out: Permittivity outside the cavity. This is characteristic of the solvent used. - - **Type** ``float`` - - **Default** ``1.0`` - - :formulation: Formulation of the Permittivity function. Currently only the exponential is used. + :formulation: Formulation of the Permittivity function. Currently only the exponential is available. **Type** ``str`` @@ -920,6 +914,28 @@ User input reference **Predicates** - ``value.lower() in ['exponential']`` + :red:`Sections` + :epsilon_out: Parameters for the continuum solvent outside the cavity. + + :red:`Keywords` + :nonequilibrium: Whether to use the nonequilibrium formulation of response, *i.e.* use the dynamic permittivity for the calculation of the response reaction field. Defaults to false. + + **Type** ``bool`` + + **Default** ``False`` + + :static: Static permittivity outside the cavity. This is characteristic of the solvent used. + + **Type** ``float`` + + **Default** ``1.0`` + + :dynamic: Dynamic permittivity outside the cavity. This is characteristic of the solvent used and relevant only in response calculations. Defaults to the same value as `epsilon_static`. + + **Type** ``float`` + + **Default** ``user['PCM']['Permittivity']['epsilon_out']['static']`` + :Constants: Physical and mathematical constants used by MRChem :red:`Keywords` diff --git a/python/mrchem/input_parser/plumbing/lexer.py b/python/mrchem/input_parser/plumbing/lexer.py index f2903f061..f502c92df 100644 --- a/python/mrchem/input_parser/plumbing/lexer.py +++ b/python/mrchem/input_parser/plumbing/lexer.py @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- -# This file was automatically generated by parselglossy on 2023-11-15 +# This file was automatically generated by parselglossy on 2023-11-22 # Editing is *STRONGLY DISCOURAGED* import json diff --git a/python/template.yml b/python/template.yml index 361767a4c..7242c41b1 100644 --- a/python/template.yml +++ b/python/template.yml @@ -936,12 +936,6 @@ sections: docstring: | Permittivity inside the cavity. 1.0 is the permittivity of free space, anything other than this is undefined behaviour. - - name: epsilon_out - type: float - default: 1.0 - docstring: | - Permittivity outside the cavity. This is characteristic of the - solvent used. - name: formulation type: str default: exponential @@ -949,7 +943,33 @@ sections: - value.lower() in ['exponential'] docstring: | Formulation of the Permittivity function. Currently only the - exponential is used. + exponential is available. + sections: + - name: epsilon_out + docstring: | + Parameters for the continuum solvent outside the cavity. + keywords: + - name: static + type: float + default: 1.0 + docstring: | + Static permittivity outside the cavity. This is characteristic of the + solvent used. + - name: dynamic + type: float + default: user['PCM']['Permittivity']['epsilon_out']['static'] + docstring: | + Dynamic permittivity outside the cavity. This is + characteristic of the solvent used and relevant only in + response calculations. Defaults to the same value as + `epsilon_static`. + - name: nonequilibrium + type: bool + default: false + docstring: | + Whether to use the nonequilibrium formulation of response, + *i.e.* use the dynamic permittivity for the calculation of the + response reaction field. Defaults to false. - name: Constants docstring: Physical and mathematical constants used by MRChem keywords: diff --git a/src/driver.cpp b/src/driver.cpp index 578034711..59c63c79a 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -104,7 +104,7 @@ namespace driver { DerivativeOperator_p get_derivative(const std::string &name); template RankOneOperator get_operator(const std::string &name, const json &json_oper); template RankTwoOperator get_operator(const std::string &name, const json &json_oper); -void build_fock_operator(const json &input, Molecule &mol, FockBuilder &F, int order); +void build_fock_operator(const json &input, Molecule &mol, FockBuilder &F, int order, bool is_dynamic = false); void init_properties(const json &json_prop, Molecule &mol); namespace scf { @@ -746,7 +746,7 @@ json driver::rsp::run(const json &json_rsp, Molecule &mol) { FockBuilder F_1; const auto &json_fock_1 = json_rsp["fock_operator"]; - driver::build_fock_operator(json_fock_1, mol, F_1, 1); + driver::build_fock_operator(json_fock_1, mol, F_1, 1, dynamic); const auto &json_pert = json_rsp["perturbation"]; auto h_1 = driver::get_operator<3>(json_pert["operator"], json_pert); @@ -983,7 +983,7 @@ void driver::rsp::calc_properties(const json &json_prop, Molecule &mol, int dir, * construct all operator which are present in this input. Option to set * perturbation order of the operators. */ -void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuilder &F, int order) { +void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuilder &F, int order, bool is_dynamic) { auto &nuclei = mol.getNuclei(); auto Phi_p = mol.getOrbitals_p(); auto X_p = mol.getOrbitalsX_p(); @@ -1043,26 +1043,66 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild /////////////////////////////////////////////////////////// if (json_fock.contains("reaction_operator")) { + // only perturbaton orders 0 and 1 are implemented + if (order > 1) MSG_ABORT("Invalid perturbation order"); + // preparing Reaction Operator auto poisson_prec = json_fock["reaction_operator"]["poisson_prec"]; auto P_p = std::make_shared(*MRA, poisson_prec); auto D_p = std::make_shared>(*MRA, 0.0, 0.0); + auto cavity_p = mol.getCavity_p(); + cavity_p->printParameters(); auto kain = json_fock["reaction_operator"]["kain"]; auto max_iter = json_fock["reaction_operator"]["max_iter"]; auto dynamic_thrs = json_fock["reaction_operator"]["dynamic_thrs"]; + // the input parser helpers have converted this parameter from string + // to integer, compatibly with the DensityType enum auto density_type = json_fock["reaction_operator"]["density_type"]; + // permittivity inside the cavity auto eps_i = json_fock["reaction_operator"]["epsilon_in"]; - auto eps_o = json_fock["reaction_operator"]["epsilon_out"]; + // static permittivity outside the cavity + auto eps_s = json_fock["reaction_operator"]["epsilon_static"]; + // dynamic permittivity outside the cavity + auto eps_d = json_fock["reaction_operator"]["epsilon_dynamic"]; + // whether nonequilibrium response was requested + auto noneq = json_fock["reaction_operator"]["nonequilibrium"]; auto formulation = json_fock["reaction_operator"]["formulation"]; - Permittivity dielectric_func(mol.getCavity_p(), eps_i, eps_o, formulation); - dielectric_func.printParameters(); + // compute nuclear charge density Density rho_nuc(false); rho_nuc = chemistry::compute_nuclear_density(poisson_prec, nuclei, 100); + // decide which permittivity to use outside of the cavity + auto eps_o = [order, is_dynamic, noneq, eps_s, eps_d] { + if (order >= 1 && noneq && is_dynamic) { + // in response (order >= 1), use dynamic permittivity if: + // a. nonequilibrium was requested (noneq is true), and + // b. the frequency is nonzero (is_dynamic is true) + return eps_d; + } else { + // for the ground state, always use the static permittivity + return eps_s; + } + }(); + + // initialize Permittivity function with static or dynamic epsilon, based on perturbation order + Permittivity dielectric_func(cavity_p, eps_i, eps_o, formulation); + dielectric_func.printParameters(); + + // initialize SCRF object auto scrf_p = std::make_unique(dielectric_func, rho_nuc, P_p, D_p, kain, max_iter, dynamic_thrs, density_type); - auto V_R = std::make_shared(std::move(scrf_p), Phi_p); + + // initialize reaction potential object + auto V_R = [&] { + if (order == 0) { + return std::make_shared(std::move(scrf_p), Phi_p); + } else { + return std::make_shared(std::move(scrf_p), Phi_p, X_p, Y_p); + } + }(); + + // set reaction potential in FockBuilder object F.getReactionOperator() = V_R; } /////////////////////////////////////////////////////////// diff --git a/src/environment/Cavity.cpp b/src/environment/Cavity.cpp index 55d5c9dbf..b7fda58a5 100644 --- a/src/environment/Cavity.cpp +++ b/src/environment/Cavity.cpp @@ -51,7 +51,7 @@ namespace detail { * where the subscript \f$i\f$ is the index related to each * sphere in the cavity, and \f$\operatorname{s}\f$ is the signed normal distance from the surface of each sphere. * @param r The coordinates of a test point in 3D space. - * @param index An integer that defines the variable of differentiation (0->x, 1->z and 2->z). + * @param index An integer that defines the variable of differentiation (0->x, 1->y and 2->z). * @param centers A vector containing the coordinates of the centers of the spheres in the cavity. * @param radii A vector containing the radii of the spheres. * @param width A double value describing the width of the transition at the boundary of the spheres. @@ -112,19 +112,6 @@ Cavity::Cavity(const std::vector> &coords, const std::vector> &coords, const std::vector &R, double sigma) - : Cavity(coords, R, std::vector(R.size(), 1.0), std::vector(R.size(), 0.0), std::vector(R.size(), sigma)) {} - /** @brief Evaluates the value of the cavity at a 3D point \f$\mathbf{r}\f$ * @param r coordinate of 3D point at which the Cavity is to be evaluated at. * @return double value of the Cavity at point \f$\mathbf{r}\f$ @@ -145,6 +132,65 @@ double Cavity::evalf(const mrcpp::Coord<3> &r) const { return C; } +void Cavity::printParameters() const { + // Collect relevant quantities + auto coords = this->centers; + auto radii = this->radii; + auto radii_0 = this->radii_0; + auto alphas = this->alphas; + auto sigmas = this->sigmas; + auto betas = this->betas; + + // Set widths + auto w0 = mrcpp::Printer::getWidth() - 1; + auto w1 = 5; + auto w2 = 9; + auto w3 = 6; + auto w4 = 10; + auto w5 = w0 - w1 - w2 - 3 * w3 - 3 * w4; + + // Build table column headers + std::stringstream o_head; + o_head << std::setw(w1) << "N"; + o_head << std::setw(w2) << "R_0"; + o_head << std::setw(w3 + 1) << "Alpha"; + o_head << std::setw(w3 - 1) << "Beta"; + o_head << std::setw(w3) << "Sigma"; + o_head << std::setw(w5) << "Radius"; + o_head << std::setw(w4) << "x"; + o_head << std::setw(w4) << "y"; + o_head << std::setw(w4) << "z"; + + // Print + mrcpp::print::header(0, "Solvation Cavity"); + println(0, o_head.str()); + mrcpp::print::separator(0, '-'); + for (auto i = 0; i < coords.size(); i++) { + auto coord = coords[i]; + auto x = coord[0]; + auto y = coord[1]; + auto z = coord[2]; + auto r = radii[i]; + auto r_0 = radii_0[i]; + auto alpha = alphas[i]; + auto beta = betas[i]; + auto sigma = sigmas[i]; + + std::stringstream o_coord; + o_coord << std::setw(w1) << i; + o_coord << std::setw(w2) << std::setprecision(4) << std::fixed << r_0; + o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << alpha; + o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << beta; + o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << sigma << " ->"; + o_coord << std::setw(w5 - 4) << std::setprecision(4) << std::fixed << r; + o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << x; + o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << y; + o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << z; + println(0, o_coord.str()); + } + mrcpp::print::separator(0, '=', 2); +} + bool Cavity::isVisibleAtScale(int scale, int nQuadPts) const { auto max_sigma = *std::max_element(this->sigmas.cbegin(), this->sigmas.cend()); diff --git a/src/environment/Cavity.h b/src/environment/Cavity.h index 0e8ca68b2..25e299a35 100644 --- a/src/environment/Cavity.h +++ b/src/environment/Cavity.h @@ -65,15 +65,23 @@ namespace mrchem { class Cavity final : public mrcpp::RepresentableFunction<3> { public: Cavity(const std::vector> &coords, const std::vector &R, const std::vector &alphas, const std::vector &betas, const std::vector &sigmas); - Cavity(const std::vector> &coords, const std::vector &R, double sigma); + /** @brief Initializes the members of the class and constructs the analytical gradient vector of the Cavity. + * + * This CTOR applies a single width factor to the cavity and **does** not modify the radii. That is, in the formula: + * + * \f[ + * R_{i} = \alpha_{i} R_{0,i} + \beta_{i}\sigma_{i} + * \f] + * + * for every atom \f$i\f$, \f$\alpha_{i} = 1.0\f$ and \f$\beta_{i} = 0.0\f$. + */ + Cavity(const std::vector> &coords, const std::vector &R, double sigma) + : Cavity(coords, R, std::vector(R.size(), 1.0), std::vector(R.size(), 0.0), std::vector(R.size(), sigma)) {} double evalf(const mrcpp::Coord<3> &r) const override; auto getGradVector() const { return this->gradvector; } - std::vector> getCoordinates() const { return this->centers; } //!< Returns #centers. - std::vector getOriginalRadii() const { return this->radii_0; } //!< Returns #radii_0. - std::vector getRadii() const { return this->radii; } //!< Returns #radii. - std::vector getRadiiScalings() const { return this->alphas; } //!< Returns #alphas. - std::vector getWidths() const { return this->sigmas; } //!< Returns #sigmas. - std::vector getWidthScalings() const { return this->betas; } //!< Returns #betas. + + /** @brief Print parameters */ + void printParameters() const; protected: std::vector radii_0; //!< Contains the *unscaled* radius of each sphere in #Center. diff --git a/src/environment/Permittivity.cpp b/src/environment/Permittivity.cpp index b0cf01716..f62d0a38b 100644 --- a/src/environment/Permittivity.cpp +++ b/src/environment/Permittivity.cpp @@ -24,11 +24,12 @@ */ #include "Permittivity.h" -#include "Cavity.h" + #include -namespace mrchem { +#include "Cavity.h" +namespace mrchem { Permittivity::Permittivity(std::shared_ptr cavity, double epsilon_in, double epsilon_out, std::string formulation) : StepFunction(cavity, epsilon_in, epsilon_out) , formulation(formulation) {} diff --git a/src/environment/Permittivity.h b/src/environment/Permittivity.h index 2d04927fd..09e73a337 100644 --- a/src/environment/Permittivity.h +++ b/src/environment/Permittivity.h @@ -25,12 +25,12 @@ #pragma once -#include "Cavity.h" -#include "StepFunction.h" -#include "utils/print_utils.h" #include #include +#include "StepFunction.h" +#include "utils/print_utils.h" + namespace mrchem { /** @class Permittivity * @@ -65,7 +65,7 @@ class Permittivity final : public StepFunction { */ double evalf(const mrcpp::Coord<3> &r) const override; - void printHeader() const override { detail::print_header("Solvation Cavity", this->formulation, getValueIn(), getValueOut()); } + void printParameters() const override { detail::print_header("Permittivity", this->formulation, getValueIn(), getValueOut()); } private: std::string formulation{"exponential"}; diff --git a/src/environment/SCRF.cpp b/src/environment/SCRF.cpp index 8aaa1eff2..308cc90ee 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/SCRF.cpp @@ -30,6 +30,9 @@ #include #include +#include "Permittivity.h" +#include "qmfunctions/Density.h" +#include "qmfunctions/Orbital.h" #include "qmfunctions/density_utils.h" #include "qmoperators/two_electron/ReactionPotential.h" #include "scf_solver/KAIN.h" @@ -42,11 +45,10 @@ using mrcpp::Timer; using PoissonOperator_p = std::shared_ptr; using DerivativeOperator_p = std::shared_ptr>; -using OrbitalVector_p = std::shared_ptr; namespace mrchem { -SCRF::SCRF(const Permittivity &e, const Density &rho_nuc, PoissonOperator_p P, DerivativeOperator_p D, int kain_hist, int max_iter, bool dyn_thrs, const std::string &density_type) +SCRF::SCRF(const Permittivity &e, const Density &rho_nuc, PoissonOperator_p P, DerivativeOperator_p D, int kain_hist, int max_iter, bool dyn_thrs, SCRFDensityType density_type) : dynamic_thrs(dyn_thrs) , density_type(density_type) , max_iter(max_iter) @@ -73,19 +75,23 @@ double SCRF::setConvergenceThreshold(double prec) { return this->conv_thrs; } -void SCRF::computeDensities(OrbitalVector &Phi, Density &rho_out) { +void SCRF::computeDensities(const Density &rho_el, Density &rho_out) { Timer timer; - Density rho_el(false); - density::compute(this->apply_prec, rho_el, Phi, DensityType::Total); - rho_el.rescale(-1.0); - - if (this->density_type == "electronic") { - mrcpp::cplxfunc::deep_copy(rho_out, rho_el); - } else if (this->density_type == "nuclear") { - mrcpp::cplxfunc::deep_copy(rho_out, this->rho_nuc); - } else { - mrcpp::cplxfunc::add(rho_out, 1.0, rho_el, 1.0, this->rho_nuc, -1.0); + switch (this->density_type) { + case SCRFDensityType::TOTAL: + mrcpp::cplxfunc::add(rho_out, 1.0, rho_el, 1.0, this->rho_nuc, -1.0); + break; + case SCRFDensityType::ELECTRONIC: + // using const_cast is absolutely EVIL here and in principle shouldn't be needed! + // however MRCPP is not everywhere const-correct, so here we go! + mrcpp::cplxfunc::deep_copy(rho_out, const_cast(rho_el)); + break; + case SCRFDensityType::NUCLEAR: + mrcpp::cplxfunc::deep_copy(rho_out, this->rho_nuc); + break; + default: + MSG_ABORT("Invalid density type"); } print_utils::qmfunction(3, "Vacuum density", rho_out, timer); } @@ -112,7 +118,7 @@ void SCRF::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunctio mrcpp::clear(d_V, true); } -mrcpp::ComplexFunction SCRF::solvePoissonEquation(const mrcpp::ComplexFunction &in_gamma, mrchem::OrbitalVector &Phi) { +mrcpp::ComplexFunction SCRF::solvePoissonEquation(const mrcpp::ComplexFunction &in_gamma, const Density &rho_el) { mrcpp::ComplexFunction Poisson_func; mrcpp::ComplexFunction rho_eff; mrcpp::ComplexFunction first_term; @@ -121,7 +127,7 @@ mrcpp::ComplexFunction SCRF::solvePoissonEquation(const mrcpp::ComplexFunction & auto eps_inv_func = mrcpp::AnalyticFunction<3>([this](const mrcpp::Coord<3> &r) { return 1.0 / this->epsilon.evalf(r); }); Density rho_tot(false); - computeDensities(Phi, rho_tot); + computeDensities(rho_el, rho_tot); mrcpp::cplxfunc::multiply(first_term, rho_tot, eps_inv_func, this->apply_prec); @@ -154,7 +160,7 @@ void SCRF::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFu dPhi_n.clear(); } -void SCRF::nestedSCRF(const mrcpp::ComplexFunction &V_vac, std::shared_ptr Phi_p) { +void SCRF::nestedSCRF(const mrcpp::ComplexFunction &V_vac, const Density &rho_el) { KAIN kain(this->history); kain.setLocalPrintLevel(10); @@ -171,7 +177,7 @@ void SCRF::nestedSCRF(const mrcpp::ComplexFunction &V_vac, std::shared_ptrVr_n, 1.0, V_vac, -1.0); computeGamma(V_tot, gamma_n); - auto Vr_np1 = solvePoissonEquation(gamma_n, *Phi_p); + auto Vr_np1 = solvePoissonEquation(gamma_n, rho_el); norm = Vr_np1.norm(); // use a convergence accelerator @@ -223,10 +229,10 @@ void SCRF::printConvergenceRow(int i, double norm, double update, double time) c println(3, o_txt.str()); } -mrcpp::ComplexFunction &SCRF::setup(double prec, const OrbitalVector_p &Phi) { +mrcpp::ComplexFunction &SCRF::setup(double prec, const Density &rho_el) { this->apply_prec = prec; Density rho_tot(false); - computeDensities(*Phi, rho_tot); + computeDensities(rho_el, rho_tot); Timer t_vac; mrcpp::ComplexFunction V_vac; V_vac.alloc(NUMBER::Real); @@ -241,13 +247,13 @@ mrcpp::ComplexFunction &SCRF::setup(double prec, const OrbitalVector_p &Phi) { mrcpp::ComplexFunction gamma_0; mrcpp::ComplexFunction V_tot; computeGamma(V_vac, gamma_0); - this->Vr_n = solvePoissonEquation(gamma_0, *Phi); + this->Vr_n = solvePoissonEquation(gamma_0, rho_el); } // update the potential/gamma before doing anything with them Timer t_scrf; - nestedSCRF(V_vac, Phi); + nestedSCRF(V_vac, rho_el); print_utils::qmfunction(3, "Reaction potential", this->Vr_n, t_scrf); return this->Vr_n; } @@ -256,7 +262,7 @@ auto SCRF::computeEnergies(const Density &rho_el) -> std::tuple auto Er_nuc = 0.5 * mrcpp::cplxfunc::dot(this->rho_nuc, this->Vr_n).real(); auto Er_el = 0.5 * mrcpp::cplxfunc::dot(rho_el, this->Vr_n).real(); - return std::make_tuple(Er_el, Er_nuc); + return {Er_el, Er_nuc}; } void SCRF::resetComplexFunction(mrcpp::ComplexFunction &function) { diff --git a/src/environment/SCRF.h b/src/environment/SCRF.h index df9312841..573edbae7 100644 --- a/src/environment/SCRF.h +++ b/src/environment/SCRF.h @@ -25,12 +25,23 @@ #pragma once +#include +#include + +#include +#include + #include "Permittivity.h" #include "qmfunctions/Density.h" -#include "qmfunctions/Orbital.h" namespace mrchem { class KAIN; +class ReactionPotential; +class ReactionPotentialD1; +class ReactionPotentialD2; + +enum class SCRFDensityType : int { TOTAL = 0, ELECTRONIC = 1, NUCLEAR = 2 }; + /** @class SCRF * @brief class that performs the computation of the ReactionPotential, named Self Consistent Reaction Field. */ @@ -43,26 +54,29 @@ class SCRF final { int kain_hist, int max_iter, bool dyn_thrs, - const std::string &density_type); + SCRFDensityType density_type); ~SCRF(); double setConvergenceThreshold(double prec); - mrcpp::ComplexFunction &getCurrentReactionPotential() { return this->Vr_n; } - Permittivity &getPermittivity() { return this->epsilon; } void updateMOResidual(double const err_t) { this->mo_residual = err_t; } - friend class ReactionPotential; auto computeEnergies(const Density &rho_el) -> std::tuple; + auto getDensityType() const -> SCRFDensityType { return this->density_type; } + + friend class ReactionPotential; + friend class ReactionPotentialD1; + friend class ReactionPotentialD2; + protected: void clear(); private: bool dynamic_thrs; - std::string density_type; + SCRFDensityType density_type; int max_iter; int history; @@ -81,15 +95,15 @@ class SCRF final { std::shared_ptr> derivative; std::shared_ptr poisson; - void computeDensities(OrbitalVector &Phi, Density &rho_out); + void computeDensities(const Density &rho_el, Density &rho_out); void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma); - mrcpp::ComplexFunction solvePoissonEquation(const mrcpp::ComplexFunction &ingamma, mrchem::OrbitalVector &Phi); + mrcpp::ComplexFunction solvePoissonEquation(const mrcpp::ComplexFunction &ingamma, const Density &rho_el); void accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain); - void nestedSCRF(const mrcpp::ComplexFunction &V_vac, std::shared_ptr Phi_p); - mrcpp::ComplexFunction &setup(double prec, const std::shared_ptr &Phi_p); + void nestedSCRF(const mrcpp::ComplexFunction &V_vac, const Density &rho_el); + mrcpp::ComplexFunction &setup(double prec, const Density &rho_el); void resetComplexFunction(mrcpp::ComplexFunction &function); diff --git a/src/environment/StepFunction.cpp b/src/environment/StepFunction.cpp index 43b628122..7d5a35e7f 100644 --- a/src/environment/StepFunction.cpp +++ b/src/environment/StepFunction.cpp @@ -22,6 +22,7 @@ * For information on the complete list of contributors to MRChem, see: * */ + #include "StepFunction.h" #include @@ -35,6 +36,7 @@ void print_header(const std::string &header, const std::string &formulation, dou print_utils::text(0, "Formulation", formulation, true); print_utils::scalar(0, "Value inside Cavity", in_value, "(in)", 6); print_utils::scalar(0, "Value outside Cavity", out_value, "(out)", 6); + mrcpp::print::separator(0, '=', 2); } } // namespace detail @@ -42,67 +44,4 @@ StepFunction::StepFunction(std::shared_ptr cavity, double val_in : in(val_in) , out(val_out) , cavity{std::move(cavity)} {} - -void StepFunction::printParameters() const { - // Collect relevant quantities - auto c_pin = this->cavity; - - auto coords = c_pin->getCoordinates(); - auto radii = c_pin->getRadii(); - auto radii_0 = c_pin->getOriginalRadii(); - auto alphas = c_pin->getRadiiScalings(); - auto sigmas = c_pin->getWidths(); - auto betas = c_pin->getWidthScalings(); - - // Set widths - auto w0 = mrcpp::Printer::getWidth() - 1; - auto w1 = 5; - auto w2 = 9; - auto w3 = 6; - auto w4 = 10; - auto w5 = w0 - w1 - w2 - 3 * w3 - 3 * w4; - - // Build table column headers - std::stringstream o_head; - o_head << std::setw(w1) << "N"; - o_head << std::setw(w2) << "R_0"; - o_head << std::setw(w3 + 1) << "Alpha"; - o_head << std::setw(w3 - 1) << "Beta"; - o_head << std::setw(w3) << "Sigma"; - o_head << std::setw(w5) << "Radius"; - o_head << std::setw(w4) << "x"; - o_head << std::setw(w4) << "y"; - o_head << std::setw(w4) << "z"; - - // Print - printHeader(); - mrcpp::print::separator(0, '-'); - println(0, o_head.str()); - mrcpp::print::separator(0, '-'); - for (auto i = 0; i < coords.size(); i++) { - auto coord = coords[i]; - auto x = coord[0]; - auto y = coord[1]; - auto z = coord[2]; - auto r = radii[i]; - auto r_0 = radii_0[i]; - auto alpha = alphas[i]; - auto beta = betas[i]; - auto sigma = sigmas[i]; - - std::stringstream o_coord; - o_coord << std::setw(w1) << i; - o_coord << std::setw(w2) << std::setprecision(4) << std::fixed << r_0; - o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << alpha; - o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << beta; - o_coord << std::setw(w3) << std::setprecision(2) << std::fixed << sigma << " ->"; - o_coord << std::setw(w5 - 4) << std::setprecision(4) << std::fixed << r; - o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << x; - o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << y; - o_coord << std::setw(w4) << std::setprecision(6) << std::fixed << z; - println(0, o_coord.str()); - } - mrcpp::print::separator(0, '=', 2); -} - } // namespace mrchem diff --git a/src/environment/StepFunction.h b/src/environment/StepFunction.h index 00bfa58e9..fb4ced35a 100644 --- a/src/environment/StepFunction.h +++ b/src/environment/StepFunction.h @@ -24,6 +24,7 @@ */ #pragma once + #include #include @@ -67,9 +68,7 @@ class StepFunction : public mrcpp::RepresentableFunction<3> { std::shared_ptr getCavity() const { return this->cavity; } - void printParameters() const; - - virtual void printHeader() const { detail::print_header("Step function of Cavity values", "Standard", getValueIn(), getValueOut()); } + virtual void printParameters() const { detail::print_header("Step function", "Standard", getValueIn(), getValueOut()); } protected: double in; //!< Value of the function inside the #cavity. diff --git a/src/qmoperators/two_electron/CMakeLists.txt b/src/qmoperators/two_electron/CMakeLists.txt index e0451ff06..62a94ab26 100644 --- a/src/qmoperators/two_electron/CMakeLists.txt +++ b/src/qmoperators/two_electron/CMakeLists.txt @@ -1,4 +1,5 @@ -target_sources(mrchem PRIVATE +target_sources(mrchem + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/CoulombPotential.cpp ${CMAKE_CURRENT_SOURCE_DIR}/CoulombPotentialD1.cpp ${CMAKE_CURRENT_SOURCE_DIR}/CoulombPotentialD2.cpp @@ -10,4 +11,6 @@ target_sources(mrchem PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/XCPotentialD1.cpp ${CMAKE_CURRENT_SOURCE_DIR}/XCPotentialD2.cpp ${CMAKE_CURRENT_SOURCE_DIR}/ReactionPotential.cpp - ) + ${CMAKE_CURRENT_SOURCE_DIR}/ReactionPotentialD1.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ReactionPotentialD2.cpp + ) diff --git a/src/qmoperators/two_electron/FockBuilder.cpp b/src/qmoperators/two_electron/FockBuilder.cpp index 8a5a512e1..60909b62f 100644 --- a/src/qmoperators/two_electron/FockBuilder.cpp +++ b/src/qmoperators/two_electron/FockBuilder.cpp @@ -173,10 +173,7 @@ SCFEnergy FockBuilder::trace(OrbitalVector &Phi, const Nuclei &nucs) { Density rho_el(false); density::compute(this->prec, rho_el, Phi, DensityType::Total); rho_el.rescale(-1.0); - std::tuple Er = this->Ro->getHelper()->computeEnergies(rho_el); - - Er_el = std::get<0>(Er); - Er_nuc = std::get<1>(Er); + std::tie(Er_el, Er_nuc) = this->Ro->getHelper()->computeEnergies(rho_el); Er_tot = Er_nuc + Er_el; } diff --git a/src/qmoperators/two_electron/ReactionOperator.h b/src/qmoperators/two_electron/ReactionOperator.h index d5ab309cc..a3d07a0e5 100644 --- a/src/qmoperators/two_electron/ReactionOperator.h +++ b/src/qmoperators/two_electron/ReactionOperator.h @@ -27,23 +27,32 @@ #include "tensor/RankZeroOperator.h" -#include "ReactionPotential.h" +#include "ReactionPotentialD1.h" +#include "ReactionPotentialD2.h" +#include "environment/SCRF.h" /** @class ReactionOperator * * @brief Operator containing a single ReactionPotential * * This class is a simple TensorOperator realization of @class ReactionPotential. - * */ namespace mrchem { -class SCRF; - class ReactionOperator final : public RankZeroOperator { public: - ReactionOperator(std::unique_ptr scrf_p, std::shared_ptr Phi_p) { - potential = std::make_shared(std::move(scrf_p), Phi_p); + ReactionOperator(std::unique_ptr scrf_p, std::shared_ptr Phi_p, bool mpi_share = false) { + potential = std::make_shared(std::move(scrf_p), Phi_p, mpi_share); + // Invoke operator= to assign *this operator + RankZeroOperator &V = (*this); + V = potential; + V.name() = "V_r"; + } + + ReactionOperator(std::unique_ptr scrf_p, std::shared_ptr Phi, std::shared_ptr X, std::shared_ptr Y, bool mpi_share = false) { + // check that the SCRF object uses the electronic density only + if (scrf_p->getDensityType() != SCRFDensityType::ELECTRONIC) MSG_ERROR("Invalid SCRF object passed: only electronic density in response"); + potential = std::make_shared(std::move(scrf_p), Phi, X, Y, mpi_share); // Invoke operator= to assign *this operator RankZeroOperator &V = (*this); V = potential; @@ -56,10 +65,7 @@ class ReactionOperator final : public RankZeroOperator { std::shared_ptr getPotential() { return this->potential; } void updateMOResidual(double const err_t) { this->potential->updateMOResidual(err_t); } - mrcpp::ComplexFunction &getCurrentReactionPotential() { return this->potential->getCurrentReactionPotential(); } - private: std::shared_ptr potential{nullptr}; }; - } // namespace mrchem diff --git a/src/qmoperators/two_electron/ReactionPotential.cpp b/src/qmoperators/two_electron/ReactionPotential.cpp index bffc3e654..7a28bb44d 100644 --- a/src/qmoperators/two_electron/ReactionPotential.cpp +++ b/src/qmoperators/two_electron/ReactionPotential.cpp @@ -27,35 +27,43 @@ #include #include "ReactionPotential.h" +#include "qmfunctions/density_utils.h" + +using mrcpp::Printer; +using mrcpp::Timer; using SCRF_p = std::unique_ptr; using OrbitalVector_p = std::shared_ptr; namespace mrchem { -ReactionPotential::ReactionPotential(SCRF_p scrf_p, OrbitalVector_p Phi_p) - : QMPotential(1, false) - , helper(std::move(scrf_p)) - , Phi(Phi_p) {} +ReactionPotential::ReactionPotential(SCRF_p scrf, OrbitalVector_p Phi, bool mpi_share) + : QMPotential(1, mpi_share) + , helper(std::move(scrf)) + , orbitals(Phi) {} void ReactionPotential::setup(double prec) { + if (isSetup(prec)) return; setApplyPrec(prec); auto thrs = this->helper->setConvergenceThreshold(prec); - mrcpp::Timer timer; - auto plevel = mrcpp::Printer::getPrintLevel(); + Timer timer; + auto plevel = Printer::getPrintLevel(); mrcpp::print::separator(3, '='); print_utils::centered_text(3, "Building Reaction operator"); this->helper->printParameters(); mrcpp::print::value(3, "Precision", prec, "(rel)", 5); mrcpp::print::value(3, "Threshold", thrs, "(abs)", 5); mrcpp::print::separator(3, '-'); - auto potential = this->helper->setup(prec, this->Phi); + + auto potential = this->computePotential(prec); + mrcpp::cplxfunc::deep_copy(*this, potential); if (plevel == 2) print_utils::qmfunction(2, "Reaction operator", *this, timer); mrcpp::print::footer(3, timer, 2); } void ReactionPotential::clear() { + mrcpp::ComplexFunction::free(NUMBER::Total); // delete FunctionTree pointers clearApplyPrec(); this->helper->clear(); } diff --git a/src/qmoperators/two_electron/ReactionPotential.h b/src/qmoperators/two_electron/ReactionPotential.h index 698fc123c..ffbbb678a 100644 --- a/src/qmoperators/two_electron/ReactionPotential.h +++ b/src/qmoperators/two_electron/ReactionPotential.h @@ -39,32 +39,34 @@ namespace mrchem { * where \f$\rho\f$ is the total molecular density of a solute molecule, \f$\epsilon\f$ is * the Permittivity function of the continuum and \f$\gamma_s\f$ is the surface charge distribution. */ -class ReactionPotential final : public QMPotential { +class ReactionPotential : public QMPotential { public: /** @brief Initializes the ReactionPotential class. - * @param scrf_p A SCRF instance which contains the parameters needed to compute the ReactionPotential. - * @param Phi_p A pointer to a vector which contains the orbitals optimized in the SCF procedure. */ - ReactionPotential(std::unique_ptr scrf_p, std::shared_ptr Phi_p); - ~ReactionPotential() override { free(NUMBER::Total); } + * @param scrf A SCRF instance which contains the parameters needed to compute the ReactionPotential. + * @param Phi A pointer to a vector which contains the orbitals optimized in the SCF procedure. + */ + explicit ReactionPotential(std::unique_ptr scrf, std::shared_ptr Phi = nullptr, bool mpi_share = false); + ~ReactionPotential() override = default; SCRF *getHelper() { return this->helper.get(); } - /** @brief Updates the helper.mo_residual member variable. This variable is used to set the convergence criterion in - * the dynamic convergence method. */ + /** @brief Updates the helper.mo_residual member variable. + * + * This variable is used to set the convergence criterion in the dynamic convergence method. + */ void updateMOResidual(double const err_t) { this->helper->mo_residual = err_t; } - mrcpp::ComplexFunction &getCurrentReactionPotential() { return this->helper->getCurrentReactionPotential(); } - friend class ReactionOperator; protected: - void clear(); + std::unique_ptr helper; ///< A SCRF instance used to compute the ReactionPotential. + std::shared_ptr orbitals; ///< Unperturbed orbitals defining the ground-state electron density for the SCRF procedure. -private: - std::unique_ptr helper; //!< A SCRF instance used to compute the ReactionPotential. - std::shared_ptr Phi; //!< holds the Orbitals needed to compute the electronic density for the SCRF procedure. + void setup(double prec) override; + void clear() override; - void setup(double prec); +private: + virtual mrcpp::ComplexFunction &computePotential(double prec) const = 0; }; } // namespace mrchem diff --git a/src/qmoperators/two_electron/ReactionPotentialD1.cpp b/src/qmoperators/two_electron/ReactionPotentialD1.cpp new file mode 100644 index 000000000..349209fb4 --- /dev/null +++ b/src/qmoperators/two_electron/ReactionPotentialD1.cpp @@ -0,0 +1,51 @@ +/* + * MRChem, a numerical real-space code for molecular electronic structure + * calculations within the self-consistent field (SCF) approximations of quantum + * chemistry (Hartree-Fock and Density Functional Theory). + * Copyright (C) 2023 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors. + * + * This file is part of MRChem. + * + * MRChem is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRChem is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRChem. If not, see . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + +#include "MRCPP/Printer" +#include "MRCPP/Timer" + +#include "ReactionPotentialD1.h" +#include "environment/SCRF.h" +#include "qmfunctions/density_utils.h" +#include "utils/print_utils.h" + +using mrcpp::Printer; +using mrcpp::Timer; + +using PoissonOperator = mrcpp::PoissonOperator; +using PoissonOperator_p = std::shared_ptr; + +namespace mrchem { +mrcpp::ComplexFunction &ReactionPotentialD1::computePotential(double prec) const { + // construct electronic density from the orbitals + OrbitalVector &Phi = *this->orbitals; + Density rho_el(false); + density::compute(this->apply_prec, rho_el, Phi, DensityType::Total); + // change sign, because it's the electronic density + rho_el.rescale(-1.0); + + return this->helper->setup(prec, rho_el); +} +} // namespace mrchem \ No newline at end of file diff --git a/src/qmoperators/two_electron/ReactionPotentialD1.h b/src/qmoperators/two_electron/ReactionPotentialD1.h new file mode 100644 index 000000000..f51da4e7e --- /dev/null +++ b/src/qmoperators/two_electron/ReactionPotentialD1.h @@ -0,0 +1,42 @@ +/* + * MRChem, a numerical real-space code for molecular electronic structure + * calculations within the self-consistent field (SCF) approximations of quantum + * chemistry (Hartree-Fock and Density Functional Theory). + * Copyright (C) 2023 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors. + * + * This file is part of MRChem. + * + * MRChem is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRChem is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRChem. If not, see . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + +#pragma once + +#include "ReactionPotential.h" +#include "environment/SCRF.h" + +namespace mrchem { + +class ReactionPotentialD1 final : public ReactionPotential { +public: + ReactionPotentialD1(std::unique_ptr scrf, std::shared_ptr Phi, bool mpi_share = false) + : ReactionPotential(std::move(scrf), Phi, mpi_share) {} + +private: + mrcpp::ComplexFunction &computePotential(double prec) const override; +}; + +} // namespace mrchem diff --git a/src/qmoperators/two_electron/ReactionPotentialD2.cpp b/src/qmoperators/two_electron/ReactionPotentialD2.cpp new file mode 100644 index 000000000..c8dea2f51 --- /dev/null +++ b/src/qmoperators/two_electron/ReactionPotentialD2.cpp @@ -0,0 +1,56 @@ +/* + * MRChem, a numerical real-space code for molecular electronic structure + * calculations within the self-consistent field (SCF) approximations of quantum + * chemistry (Hartree-Fock and Density Functional Theory). + * Copyright (C) 2023 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors. + * + * This file is part of MRChem. + * + * MRChem is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRChem is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRChem. If not, see . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + +#include "MRCPP/Printer" +#include "MRCPP/Timer" + +#include "ReactionPotentialD2.h" +#include "qmfunctions/density_utils.h" +#include "utils/print_utils.h" + +using mrcpp::Printer; +using mrcpp::Timer; + +namespace mrchem { +mrcpp::ComplexFunction &ReactionPotentialD2::computePotential(double prec) const { + // construct perturbed density from the orbitals + if (this->orbitals == nullptr) MSG_ERROR("Orbitals not initialized"); + if (this->orbitals_x == nullptr) MSG_ERROR("Perturbed X orbitals not initialized"); + if (this->orbitals_y == nullptr) MSG_ERROR("Perturbed Y orbitals not initialized"); + + Density rho(false); + OrbitalVector &Phi = *this->orbitals; + OrbitalVector &X = *this->orbitals_x; + OrbitalVector &Y = *this->orbitals_y; + + Timer timer; + density::compute(prec, rho, Phi, X, Y, DensityType::Total); + print_utils::qmfunction(3, "Compute global density", rho, timer); + // change sign, because it's the electronic density + rho.rescale(-1.0); + + return this->helper->setup(prec, rho); +} +} // namespace mrchem diff --git a/src/qmoperators/two_electron/ReactionPotentialD2.h b/src/qmoperators/two_electron/ReactionPotentialD2.h new file mode 100644 index 000000000..8cc846d2f --- /dev/null +++ b/src/qmoperators/two_electron/ReactionPotentialD2.h @@ -0,0 +1,45 @@ +/* + * MRChem, a numerical real-space code for molecular electronic structure + * calculations within the self-consistent field (SCF) approximations of quantum + * chemistry (Hartree-Fock and Density Functional Theory). + * Copyright (C) 2023 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors. + * + * This file is part of MRChem. + * + * MRChem is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * MRChem is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with MRChem. If not, see . + * + * For information on the complete list of contributors to MRChem, see: + * + */ + +#pragma once + +#include "ReactionPotential.h" +#include "environment/SCRF.h" + +namespace mrchem { +class ReactionPotentialD2 final : public ReactionPotential { +public: + ReactionPotentialD2(std::unique_ptr scrf, std::shared_ptr Phi, std::shared_ptr X, std::shared_ptr Y, bool mpi_share = false) + : ReactionPotential(std::move(scrf), Phi, mpi_share) + , orbitals_x(X) + , orbitals_y(Y) {} + +private: + std::shared_ptr orbitals_x; ///< Perturbed orbitals + std::shared_ptr orbitals_y; ///< Perturbed orbitals + + mrcpp::ComplexFunction &computePotential(double prec) const override; +}; +} // namespace mrchem diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index e1f40990d..cc8929505 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -31,3 +31,4 @@ add_subdirectory(h2_scf_cube) add_subdirectory(li_solv) add_subdirectory(he_zora_scf_lda) add_subdirectory(cavity_input_parser) +add_subdirectory(h2_pol_solv) diff --git a/tests/cavity_input_parser/world_unit=angstrom-default.inp b/tests/cavity_input_parser/world_unit=angstrom-default.inp index 7807bdbdc..e451f93a4 100644 --- a/tests/cavity_input_parser/world_unit=angstrom-default.inp +++ b/tests/cavity_input_parser/world_unit=angstrom-default.inp @@ -19,7 +19,9 @@ WaveFunction { PCM { Permittivity { - epsilon_out = 80.0 + outside = { + epsilon_static = 80.0 + } } } diff --git a/tests/cavity_input_parser/world_unit=angstrom-mode=atoms-append.inp b/tests/cavity_input_parser/world_unit=angstrom-mode=atoms-append.inp index 14b71213c..ab2a50307 100644 --- a/tests/cavity_input_parser/world_unit=angstrom-mode=atoms-append.inp +++ b/tests/cavity_input_parser/world_unit=angstrom-mode=atoms-append.inp @@ -19,7 +19,9 @@ WaveFunction { PCM { Permittivity { - epsilon_out = 80.0 + outside = { + epsilon_static = 80.0 + } } Cavity { mode = atoms diff --git a/tests/cavity_input_parser/world_unit=angstrom-mode=atoms-replace.inp b/tests/cavity_input_parser/world_unit=angstrom-mode=atoms-replace.inp index a17fb8615..9a8b14791 100644 --- a/tests/cavity_input_parser/world_unit=angstrom-mode=atoms-replace.inp +++ b/tests/cavity_input_parser/world_unit=angstrom-mode=atoms-replace.inp @@ -19,7 +19,9 @@ WaveFunction { PCM { Permittivity { - epsilon_out = 80.0 + outside = { + epsilon_static = 80.0 + } } Cavity { mode = atoms diff --git a/tests/cavity_input_parser/world_unit=angstrom-mode=explicit.inp b/tests/cavity_input_parser/world_unit=angstrom-mode=explicit.inp index ba1e986fc..e8dae2340 100644 --- a/tests/cavity_input_parser/world_unit=angstrom-mode=explicit.inp +++ b/tests/cavity_input_parser/world_unit=angstrom-mode=explicit.inp @@ -19,7 +19,9 @@ WaveFunction { PCM { Permittivity { - epsilon_out = 80.0 + outside = { + epsilon_static = 80.0 + } } Cavity { mode = explicit diff --git a/tests/cavity_input_parser/world_unit=bohr-default.inp b/tests/cavity_input_parser/world_unit=bohr-default.inp index 1f476088a..65fe0526d 100644 --- a/tests/cavity_input_parser/world_unit=bohr-default.inp +++ b/tests/cavity_input_parser/world_unit=bohr-default.inp @@ -19,7 +19,9 @@ WaveFunction { PCM { Permittivity { - epsilon_out = 80.0 + outside = { + epsilon_static = 80.0 + } } } diff --git a/tests/cavity_input_parser/world_unit=bohr-mode=atoms-append.inp b/tests/cavity_input_parser/world_unit=bohr-mode=atoms-append.inp index 2c9b1ff7b..a84cc0139 100644 --- a/tests/cavity_input_parser/world_unit=bohr-mode=atoms-append.inp +++ b/tests/cavity_input_parser/world_unit=bohr-mode=atoms-append.inp @@ -19,7 +19,9 @@ WaveFunction { PCM { Permittivity { - epsilon_out = 80.0 + outside = { + epsilon_static = 80.0 + } } Cavity { mode = atoms diff --git a/tests/cavity_input_parser/world_unit=bohr-mode=atoms-replace.inp b/tests/cavity_input_parser/world_unit=bohr-mode=atoms-replace.inp index 1e769be17..776ff230d 100644 --- a/tests/cavity_input_parser/world_unit=bohr-mode=atoms-replace.inp +++ b/tests/cavity_input_parser/world_unit=bohr-mode=atoms-replace.inp @@ -19,7 +19,9 @@ WaveFunction { PCM { Permittivity { - epsilon_out = 80.0 + outside = { + epsilon_static = 80.0 + } } Cavity { mode = atoms diff --git a/tests/cavity_input_parser/world_unit=bohr-mode=explicit.inp b/tests/cavity_input_parser/world_unit=bohr-mode=explicit.inp index 07ffd2dcd..f1618596d 100644 --- a/tests/cavity_input_parser/world_unit=bohr-mode=explicit.inp +++ b/tests/cavity_input_parser/world_unit=bohr-mode=explicit.inp @@ -19,7 +19,9 @@ WaveFunction { PCM { Permittivity { - epsilon_out = 80.0 + outside = { + epsilon_static = 80.0 + } } Cavity { mode = explicit diff --git a/tests/h2_pol_solv/CMakeLists.txt b/tests/h2_pol_solv/CMakeLists.txt new file mode 100644 index 000000000..62fb9b009 --- /dev/null +++ b/tests/h2_pol_solv/CMakeLists.txt @@ -0,0 +1,10 @@ +if(ENABLE_MPI) + set(_h2_pol_solv_launcher "${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 1") +endif() + +add_integration_test( + NAME "H2_polarizability_solvent_effect" + LABELS "mrchem;h2_pol_solv;solvent;scf;energy;polarizability" + COST 150 + LAUNCH_AGENT ${_h2_pol_solv_launcher} + ) diff --git a/tests/h2_pol_solv/h2-eq.inp b/tests/h2_pol_solv/h2-eq.inp new file mode 100644 index 000000000..88c2556d8 --- /dev/null +++ b/tests/h2_pol_solv/h2-eq.inp @@ -0,0 +1,45 @@ +world_prec = 1.0e-3 +world_unit = angstrom + +Molecule { +$coords +H 0.0000 0.0000 -0.3705 +H 0.0000 0.0000 0.3705 +$end +} + +WaveFunction { + method = HF + environment = PCM +} + +Properties { + polarizability = true +} + +PCM { + Permittivity { + epsilon_out { + static = 4.0 + dynamic = 2.0 + nonequilibrium = false + } + } +} + +SCF { + run = false + guess_type = SAD_DZ +} + +Polarizability { + frequency = [0.0, 0.05] +} + +Response { + kain = 3 + max_iter = 10 + orbital_thrs = 0.01 + run = [false, false, true] +} +~ diff --git a/tests/h2_pol_solv/h2-neq.inp b/tests/h2_pol_solv/h2-neq.inp new file mode 100644 index 000000000..7af73eced --- /dev/null +++ b/tests/h2_pol_solv/h2-neq.inp @@ -0,0 +1,45 @@ +world_prec = 1.0e-3 +world_unit = angstrom + +Molecule { +$coords +H 0.0000 0.0000 -0.3705 +H 0.0000 0.0000 0.3705 +$end +} + +WaveFunction { + method = HF + environment = PCM +} + +Properties { + polarizability = true +} + +PCM { + Permittivity { + epsilon_out { + static = 4.0 + dynamic = 2.0 + nonequilibrium = true + } + } +} + +SCF { + run = false + guess_type = SAD_DZ +} + +Polarizability { + frequency = [0.0, 0.05] +} + +Response { + kain = 3 + max_iter = 10 + orbital_thrs = 0.01 + run = [false, false, true] +} +~ diff --git a/tests/h2_pol_solv/reference/README.md b/tests/h2_pol_solv/reference/README.md new file mode 100644 index 000000000..ae3475275 --- /dev/null +++ b/tests/h2_pol_solv/reference/README.md @@ -0,0 +1,12 @@ +# Checking the polarizability with finite differences + +This folder contains input files for checking the ZZ component of the static +polarizability by finite differences: + +- `h2_m.inp`: electric field with -0.01 field strength. +- `h2_p.inp`: electric field with +0.01 field strength. +- `h2.inp`: polarizability. + +After running the 3 MRChem calculations, run the `check-pol.py` script to +extract and compare the finite difference and analytical polarizabilities. + diff --git a/tests/h2_pol_solv/reference/check-pol.py b/tests/h2_pol_solv/reference/check-pol.py new file mode 100644 index 000000000..22f0fce91 --- /dev/null +++ b/tests/h2_pol_solv/reference/check-pol.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python3 + +from json import load +from pathlib import Path + +dipoles = [] +ef_strength = 0.01 + +curdir = Path(__file__).parent.resolve() + +for out in ["h2_m", "h2_p"]: + with (curdir / f"{out}.json").open("r") as fh: + data = load(fh) + dipoles.append( + data["output"]["properties"]["dipole_moment"]["dip-1"]["vector_el"][-1] + ) + +fd_pol = (dipoles[1] - dipoles[0]) / (2 * ef_strength) + +with (curdir / "h2.json").open("r") as fh: + data = load(fh) + pol = data["output"]["properties"]["polarizability"]["pol-0.000000"]["tensor"][-1] + +print(f"Finite difference: {fd_pol}") +print(f"Analytical: {pol}") diff --git a/tests/h2_pol_solv/reference/h2-eq.json b/tests/h2_pol_solv/reference/h2-eq.json new file mode 100644 index 000000000..0f553f789 --- /dev/null +++ b/tests/h2_pol_solv/reference/h2-eq.json @@ -0,0 +1,736 @@ +{ + "input": { + "constants": { + "angstrom2bohrs": 1.8897261246257702, + "dipmom_au2debye": 2.5417464739297717, + "electron_g_factor": -2.00231930436256, + "fine_structure_constant": 0.0072973525693, + "hartree2ev": 27.211386245988, + "hartree2kcalmol": 627.5094740630558, + "hartree2kjmol": 2625.4996394798254, + "hartree2simagnetizability": 78.9451185, + "hartree2wavenumbers": 219474.6313632, + "light_speed": 137.035999084 + }, + "molecule": { + "cavity": { + "spheres": [ + { + "alpha": 1.1, + "beta": 0.5, + "center": [ + 0.0, + 0.0, + -0.7001435291738478 + ], + "radius": 2.267671349550924, + "sigma": 0.2 + }, + { + "alpha": 1.1, + "beta": 0.5, + "center": [ + 0.0, + 0.0, + 0.7001435291738478 + ], + "radius": 2.267671349550924, + "sigma": 0.2 + } + ] + }, + "charge": 0, + "coords": [ + { + "atom": "h", + "r_rms": 2.6569547399e-05, + "xyz": [ + 0.0, + 0.0, + -0.7001435291738478 + ] + }, + { + "atom": "h", + "r_rms": 2.6569547399e-05, + "xyz": [ + 0.0, + 0.0, + 0.7001435291738478 + ] + } + ], + "multiplicity": 1 + }, + "mpi": { + "bank_size": -1, + "numerically_exact": false, + "shared_memory_size": 10000 + }, + "mra": { + "basis_order": 5, + "basis_type": "interpolating", + "boxes": [ + 2, + 2, + 2 + ], + "corner": [ + -1, + -1, + -1 + ], + "max_scale": 20, + "min_scale": -4 + }, + "printer": { + "file_name": "h2-eq", + "print_constants": false, + "print_level": 0, + "print_mpi": false, + "print_prec": 6, + "print_width": 75 + }, + "rsp_calculations": { + "ext_el-0.000000": { + "components": [ + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_0_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_0_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_0_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_0_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_0_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_0_vector.json", + "file_chk_x": "checkpoint/X_rsp_0", + "file_chk_y": "checkpoint/Y_rsp_0", + "file_x_a": "initial_guess/X_a_rsp_0", + "file_x_b": "initial_guess/X_b_rsp_0", + "file_x_p": "initial_guess/X_p_rsp_0", + "file_y_a": "initial_guess/Y_a_rsp_0", + "file_y_b": "initial_guess/Y_b_rsp_0", + "file_y_p": "initial_guess/Y_p_rsp_0", + "prec": 0.001, + "type": "none" + } + }, + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_1_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_1_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_1_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_1_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_1_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_1_vector.json", + "file_chk_x": "checkpoint/X_rsp_1", + "file_chk_y": "checkpoint/Y_rsp_1", + "file_x_a": "initial_guess/X_a_rsp_1", + "file_x_b": "initial_guess/X_b_rsp_1", + "file_x_p": "initial_guess/X_p_rsp_1", + "file_y_a": "initial_guess/Y_a_rsp_1", + "file_y_b": "initial_guess/Y_b_rsp_1", + "file_y_p": "initial_guess/Y_p_rsp_1", + "prec": 0.001, + "type": "none" + } + }, + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_2_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_2_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_2_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_2_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_2_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_2_vector.json", + "file_chk_x": "checkpoint/X_rsp_2", + "file_chk_y": "checkpoint/Y_rsp_2", + "file_x_a": "initial_guess/X_a_rsp_2", + "file_x_b": "initial_guess/X_b_rsp_2", + "file_x_p": "initial_guess/X_p_rsp_2", + "file_y_a": "initial_guess/Y_a_rsp_2", + "file_y_b": "initial_guess/Y_b_rsp_2", + "file_y_p": "initial_guess/Y_p_rsp_2", + "prec": 0.001, + "type": "none" + }, + "rsp_solver": { + "checkpoint": false, + "file_chk_x": "checkpoint/X_rsp_2", + "file_chk_y": "checkpoint/Y_rsp_2", + "final_prec": 0.001, + "helmholtz_prec": -1.0, + "kain": 3, + "max_iter": 10, + "method": "Hartree-Fock", + "orbital_thrs": 0.01, + "orth_prec": 1e-14, + "property_thrs": -1.0, + "start_prec": 0.001 + } + } + ], + "dynamic": false, + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "reaction_operator": { + "density_type": 1, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": false, + "poisson_prec": 0.001 + } + }, + "frequency": 0.0, + "perturbation": { + "operator": "h_e_dip", + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + }, + "properties": { + "polarizability": { + "pol-0.000000": { + "frequency": 0.0, + "operator": "h_e_dip", + "precision": 0.001, + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + } + } + }, + "unperturbed": { + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "kinetic_operator": { + "derivative": "abgv_55" + }, + "nuclear_operator": { + "nuclear_model": "point_like", + "proj_prec": 0.001, + "shared_memory": false, + "smooth_prec": 0.001 + }, + "reaction_operator": { + "density_type": 0, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": false, + "poisson_prec": 0.001 + } + }, + "localize": false, + "precision": 0.001 + } + }, + "ext_el-0.050000": { + "components": [ + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_0_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_0_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_0_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_0_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_0_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_0_vector.json", + "file_chk_x": "checkpoint/X_rsp_0", + "file_chk_y": "checkpoint/Y_rsp_0", + "file_x_a": "initial_guess/X_a_rsp_0", + "file_x_b": "initial_guess/X_b_rsp_0", + "file_x_p": "initial_guess/X_p_rsp_0", + "file_y_a": "initial_guess/Y_a_rsp_0", + "file_y_b": "initial_guess/Y_b_rsp_0", + "file_y_p": "initial_guess/Y_p_rsp_0", + "prec": 0.001, + "type": "none" + } + }, + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_1_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_1_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_1_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_1_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_1_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_1_vector.json", + "file_chk_x": "checkpoint/X_rsp_1", + "file_chk_y": "checkpoint/Y_rsp_1", + "file_x_a": "initial_guess/X_a_rsp_1", + "file_x_b": "initial_guess/X_b_rsp_1", + "file_x_p": "initial_guess/X_p_rsp_1", + "file_y_a": "initial_guess/Y_a_rsp_1", + "file_y_b": "initial_guess/Y_b_rsp_1", + "file_y_p": "initial_guess/Y_p_rsp_1", + "prec": 0.001, + "type": "none" + } + }, + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_2_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_2_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_2_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_2_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_2_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_2_vector.json", + "file_chk_x": "checkpoint/X_rsp_2", + "file_chk_y": "checkpoint/Y_rsp_2", + "file_x_a": "initial_guess/X_a_rsp_2", + "file_x_b": "initial_guess/X_b_rsp_2", + "file_x_p": "initial_guess/X_p_rsp_2", + "file_y_a": "initial_guess/Y_a_rsp_2", + "file_y_b": "initial_guess/Y_b_rsp_2", + "file_y_p": "initial_guess/Y_p_rsp_2", + "prec": 0.001, + "type": "none" + }, + "rsp_solver": { + "checkpoint": false, + "file_chk_x": "checkpoint/X_rsp_2", + "file_chk_y": "checkpoint/Y_rsp_2", + "final_prec": 0.001, + "helmholtz_prec": -1.0, + "kain": 3, + "max_iter": 10, + "method": "Hartree-Fock", + "orbital_thrs": 0.01, + "orth_prec": 1e-14, + "property_thrs": -1.0, + "start_prec": 0.001 + } + } + ], + "dynamic": true, + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "reaction_operator": { + "density_type": 1, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": false, + "poisson_prec": 0.001 + } + }, + "frequency": 0.05, + "perturbation": { + "operator": "h_e_dip", + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + }, + "properties": { + "polarizability": { + "pol-0.050000": { + "frequency": 0.05, + "operator": "h_e_dip", + "precision": 0.001, + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + } + } + }, + "unperturbed": { + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "kinetic_operator": { + "derivative": "abgv_55" + }, + "nuclear_operator": { + "nuclear_model": "point_like", + "proj_prec": 0.001, + "shared_memory": false, + "smooth_prec": 0.001 + }, + "reaction_operator": { + "density_type": 0, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": false, + "poisson_prec": 0.001 + } + }, + "localize": false, + "precision": 0.001 + } + } + }, + "scf_calculation": { + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "kinetic_operator": { + "derivative": "abgv_55" + }, + "nuclear_operator": { + "nuclear_model": "point_like", + "proj_prec": 0.001, + "shared_memory": false, + "smooth_prec": 0.001 + }, + "reaction_operator": { + "density_type": 0, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": false, + "poisson_prec": 0.001 + } + }, + "initial_guess": { + "environment": "PCM", + "external_field": "None", + "file_CUBE_a": "cube_vectors/CUBE_a_vector.json", + "file_CUBE_b": "cube_vectors/CUBE_b_vector.json", + "file_CUBE_p": "cube_vectors/CUBE_p_vector.json", + "file_basis": "initial_guess/mrchem.bas", + "file_chk": "checkpoint/phi_scf", + "file_gto_a": "initial_guess/mrchem.moa", + "file_gto_b": "initial_guess/mrchem.mob", + "file_gto_p": "initial_guess/mrchem.mop", + "file_phi_a": "initial_guess/phi_a_scf", + "file_phi_b": "initial_guess/phi_b_scf", + "file_phi_p": "initial_guess/phi_p_scf", + "localize": false, + "method": "Hartree-Fock", + "prec": 0.001, + "relativity": "None", + "restricted": true, + "screen": 12.0, + "type": "sad", + "zeta": 2 + }, + "properties": { + "dipole_moment": { + "dip-1": { + "operator": "h_e_dip", + "precision": 0.001, + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + } + } + } + }, + "schema_name": "mrchem_input", + "schema_version": 1 + }, + "output": { + "properties": { + "center_of_mass": [ + 0.0, + 0.0, + 0.0 + ], + "charge": 0, + "dipole_moment": { + "dip-1": { + "magnitude": 0.001436455029867261, + "r_O": [ + 0.0, + 0.0, + 0.0 + ], + "vector": [ + 0.0, + 0.0, + 0.001436455029867261 + ], + "vector_el": [ + 0.0, + 0.0, + 0.0014364550298491643 + ], + "vector_nuc": [ + 0.0, + 0.0, + 0.0 + ] + } + }, + "geometry": [ + { + "symbol": "H", + "xyz": [ + 0.0, + 0.0, + -0.7001435291738478 + ] + }, + { + "symbol": "H", + "xyz": [ + 0.0, + 0.0, + 0.7001435291738478 + ] + } + ], + "multiplicity": 1, + "orbital_energies": { + "energy": [ + -0.5985110862667959 + ], + "occupation": [ + 2.0 + ], + "spin": [ + "p" + ], + "sum_occupied": -1.1970221725335919 + }, + "polarizability": { + "pol-0.000000": { + "frequency": 0.0, + "isotropic_average": null, + "r_O": [ + 0.0, + 0.0, + 0.0 + ], + "tensor": [ + null, + null, + null, + null, + null, + null, + 0.0, + 0.0, + 9.607625115716449 + ] + }, + "pol-0.050000": { + "frequency": 0.05, + "isotropic_average": null, + "r_O": [ + 0.0, + 0.0, + 0.0 + ], + "tensor": [ + null, + null, + null, + null, + null, + null, + 0.0, + 0.0, + 9.72197933526629 + ] + } + }, + "scf_energy": { + "E_ee": 1.1798069065653403, + "E_eext": 0.0, + "E_el": -1.7997950754459169, + "E_en": -3.2910350186754416, + "E_kin": 0.8884670403157098, + "E_next": 0.0, + "E_nn": 0.7141392859689609, + "E_nuc": 0.7030989718207105, + "E_tot": -1.0966961036252063, + "E_x": -0.5877514493006317, + "E_xc": 0.0, + "Er_el": 0.01071744564910627, + "Er_nuc": -0.01104031414825042, + "Er_tot": -0.0003228684991441501 + } + }, + "provenance": { + "creator": "MRChem", + "mpi_processes": 1, + "nthreads": 1, + "routine": "mrchem.x", + "total_cores": 1, + "version": "1.2.0-alpha" + }, + "rsp_calculations": { + "ext_el-0.000000": { + "components": [ + null, + null, + { + "rsp_solver": { + "converged": true, + "cycles": [ + { + "mo_residual": 1.1840373878577934, + "property_update": -5.230056480276378, + "symmetric_property": -5.230056480276378, + "wall_time": 1.589366237 + }, + { + "mo_residual": 0.5039164949499658, + "property_update": -3.896545272260962, + "symmetric_property": -9.12660175253734, + "wall_time": 10.365826822 + }, + { + "mo_residual": 0.06312576948642024, + "property_update": -0.47827471818918177, + "symmetric_property": -9.604876470726522, + "wall_time": 9.840587724 + }, + { + "mo_residual": 0.01373468159902534, + "property_update": 0.020830591804632093, + "symmetric_property": -9.58404587892189, + "wall_time": 8.427857813 + }, + { + "mo_residual": 0.0033793658719311593, + "property_update": -0.023579236794558867, + "symmetric_property": -9.607625115716449, + "wall_time": 8.481558802 + } + ], + "wall_time": 38.70527487 + } + } + ], + "frequency": 0.0, + "perturbation": "h_e_dip", + "success": true + }, + "ext_el-0.050000": { + "components": [ + null, + null, + { + "rsp_solver": { + "converged": true, + "cycles": [ + { + "mo_residual": 1.7320741496375414, + "property_update": -5.425174312715582, + "symmetric_property": -5.425174312715582, + "wall_time": 3.228155599 + }, + { + "mo_residual": 0.7407765228467756, + "property_update": -4.061079791293976, + "symmetric_property": -9.486254104009557, + "wall_time": 13.582917248 + }, + { + "mo_residual": 0.08287169934740206, + "property_update": -0.22335149625176598, + "symmetric_property": -9.709605600261323, + "wall_time": 13.522357354 + }, + { + "mo_residual": 0.02379890307059201, + "property_update": 0.02381557327733752, + "symmetric_property": -9.685790026983986, + "wall_time": 12.259150599 + }, + { + "mo_residual": 0.007573242416437982, + "property_update": -0.03618930828230482, + "symmetric_property": -9.72197933526629, + "wall_time": 12.61330521 + } + ], + "wall_time": 55.205957909 + } + } + ], + "frequency": 0.05, + "perturbation": "h_e_dip", + "success": true + } + }, + "scf_calculation": { + "initial_energy": { + "E_ee": 1.1798069065653403, + "E_eext": 0.0, + "E_el": -1.7997950754459169, + "E_en": -3.2910350186754416, + "E_kin": 0.8884670403157098, + "E_next": 0.0, + "E_nn": 0.7141392859689609, + "E_nuc": 0.7030989718207105, + "E_tot": -1.0966961036252063, + "E_x": -0.5877514493006317, + "E_xc": 0.0, + "Er_el": 0.01071744564910627, + "Er_nuc": -0.01104031414825042, + "Er_tot": -0.0003228684991441501 + }, + "success": true + }, + "schema_name": "mrchem_output", + "schema_version": 1, + "success": true + } +} diff --git a/tests/h2_pol_solv/reference/h2-neq.json b/tests/h2_pol_solv/reference/h2-neq.json new file mode 100644 index 000000000..74be4ceee --- /dev/null +++ b/tests/h2_pol_solv/reference/h2-neq.json @@ -0,0 +1,736 @@ +{ + "input": { + "constants": { + "angstrom2bohrs": 1.8897261246257702, + "dipmom_au2debye": 2.5417464739297717, + "electron_g_factor": -2.00231930436256, + "fine_structure_constant": 0.0072973525693, + "hartree2ev": 27.211386245988, + "hartree2kcalmol": 627.5094740630558, + "hartree2kjmol": 2625.4996394798254, + "hartree2simagnetizability": 78.9451185, + "hartree2wavenumbers": 219474.6313632, + "light_speed": 137.035999084 + }, + "molecule": { + "cavity": { + "spheres": [ + { + "alpha": 1.1, + "beta": 0.5, + "center": [ + 0.0, + 0.0, + -0.7001435291738478 + ], + "radius": 2.267671349550924, + "sigma": 0.2 + }, + { + "alpha": 1.1, + "beta": 0.5, + "center": [ + 0.0, + 0.0, + 0.7001435291738478 + ], + "radius": 2.267671349550924, + "sigma": 0.2 + } + ] + }, + "charge": 0, + "coords": [ + { + "atom": "h", + "r_rms": 2.6569547399e-05, + "xyz": [ + 0.0, + 0.0, + -0.7001435291738478 + ] + }, + { + "atom": "h", + "r_rms": 2.6569547399e-05, + "xyz": [ + 0.0, + 0.0, + 0.7001435291738478 + ] + } + ], + "multiplicity": 1 + }, + "mpi": { + "bank_size": -1, + "numerically_exact": false, + "shared_memory_size": 10000 + }, + "mra": { + "basis_order": 5, + "basis_type": "interpolating", + "boxes": [ + 2, + 2, + 2 + ], + "corner": [ + -1, + -1, + -1 + ], + "max_scale": 20, + "min_scale": -4 + }, + "printer": { + "file_name": "h2-neq", + "print_constants": false, + "print_level": 0, + "print_mpi": false, + "print_prec": 6, + "print_width": 75 + }, + "rsp_calculations": { + "ext_el-0.000000": { + "components": [ + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_0_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_0_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_0_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_0_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_0_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_0_vector.json", + "file_chk_x": "checkpoint/X_rsp_0", + "file_chk_y": "checkpoint/Y_rsp_0", + "file_x_a": "initial_guess/X_a_rsp_0", + "file_x_b": "initial_guess/X_b_rsp_0", + "file_x_p": "initial_guess/X_p_rsp_0", + "file_y_a": "initial_guess/Y_a_rsp_0", + "file_y_b": "initial_guess/Y_b_rsp_0", + "file_y_p": "initial_guess/Y_p_rsp_0", + "prec": 0.001, + "type": "none" + } + }, + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_1_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_1_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_1_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_1_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_1_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_1_vector.json", + "file_chk_x": "checkpoint/X_rsp_1", + "file_chk_y": "checkpoint/Y_rsp_1", + "file_x_a": "initial_guess/X_a_rsp_1", + "file_x_b": "initial_guess/X_b_rsp_1", + "file_x_p": "initial_guess/X_p_rsp_1", + "file_y_a": "initial_guess/Y_a_rsp_1", + "file_y_b": "initial_guess/Y_b_rsp_1", + "file_y_p": "initial_guess/Y_p_rsp_1", + "prec": 0.001, + "type": "none" + } + }, + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_2_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_2_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_2_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_2_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_2_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_2_vector.json", + "file_chk_x": "checkpoint/X_rsp_2", + "file_chk_y": "checkpoint/Y_rsp_2", + "file_x_a": "initial_guess/X_a_rsp_2", + "file_x_b": "initial_guess/X_b_rsp_2", + "file_x_p": "initial_guess/X_p_rsp_2", + "file_y_a": "initial_guess/Y_a_rsp_2", + "file_y_b": "initial_guess/Y_b_rsp_2", + "file_y_p": "initial_guess/Y_p_rsp_2", + "prec": 0.001, + "type": "none" + }, + "rsp_solver": { + "checkpoint": false, + "file_chk_x": "checkpoint/X_rsp_2", + "file_chk_y": "checkpoint/Y_rsp_2", + "final_prec": 0.001, + "helmholtz_prec": -1.0, + "kain": 3, + "max_iter": 10, + "method": "Hartree-Fock", + "orbital_thrs": 0.01, + "orth_prec": 1e-14, + "property_thrs": -1.0, + "start_prec": 0.001 + } + } + ], + "dynamic": false, + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "reaction_operator": { + "density_type": 1, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": true, + "poisson_prec": 0.001 + } + }, + "frequency": 0.0, + "perturbation": { + "operator": "h_e_dip", + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + }, + "properties": { + "polarizability": { + "pol-0.000000": { + "frequency": 0.0, + "operator": "h_e_dip", + "precision": 0.001, + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + } + } + }, + "unperturbed": { + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "kinetic_operator": { + "derivative": "abgv_55" + }, + "nuclear_operator": { + "nuclear_model": "point_like", + "proj_prec": 0.001, + "shared_memory": false, + "smooth_prec": 0.001 + }, + "reaction_operator": { + "density_type": 0, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": true, + "poisson_prec": 0.001 + } + }, + "localize": false, + "precision": 0.001 + } + }, + "ext_el-0.050000": { + "components": [ + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_0_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_0_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_0_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_0_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_0_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_0_vector.json", + "file_chk_x": "checkpoint/X_rsp_0", + "file_chk_y": "checkpoint/Y_rsp_0", + "file_x_a": "initial_guess/X_a_rsp_0", + "file_x_b": "initial_guess/X_b_rsp_0", + "file_x_p": "initial_guess/X_p_rsp_0", + "file_y_a": "initial_guess/Y_a_rsp_0", + "file_y_b": "initial_guess/Y_b_rsp_0", + "file_y_p": "initial_guess/Y_p_rsp_0", + "prec": 0.001, + "type": "none" + } + }, + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_1_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_1_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_1_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_1_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_1_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_1_vector.json", + "file_chk_x": "checkpoint/X_rsp_1", + "file_chk_y": "checkpoint/Y_rsp_1", + "file_x_a": "initial_guess/X_a_rsp_1", + "file_x_b": "initial_guess/X_b_rsp_1", + "file_x_p": "initial_guess/X_p_rsp_1", + "file_y_a": "initial_guess/Y_a_rsp_1", + "file_y_b": "initial_guess/Y_b_rsp_1", + "file_y_p": "initial_guess/Y_p_rsp_1", + "prec": 0.001, + "type": "none" + } + }, + { + "initial_guess": { + "file_CUBE_x_a": "cube_vectors/CUBE_x_a_2_vector.json", + "file_CUBE_x_b": "cube_vectors/CUBE_x_b_2_vector.json", + "file_CUBE_x_p": "cube_vectors/CUBE_x_p_2_vector.json", + "file_CUBE_y_a": "cube_vectors/CUBE_y_a_2_vector.json", + "file_CUBE_y_b": "cube_vectors/CUBE_y_b_2_vector.json", + "file_CUBE_y_p": "cube_vectors/CUBE_y_p_2_vector.json", + "file_chk_x": "checkpoint/X_rsp_2", + "file_chk_y": "checkpoint/Y_rsp_2", + "file_x_a": "initial_guess/X_a_rsp_2", + "file_x_b": "initial_guess/X_b_rsp_2", + "file_x_p": "initial_guess/X_p_rsp_2", + "file_y_a": "initial_guess/Y_a_rsp_2", + "file_y_b": "initial_guess/Y_b_rsp_2", + "file_y_p": "initial_guess/Y_p_rsp_2", + "prec": 0.001, + "type": "none" + }, + "rsp_solver": { + "checkpoint": false, + "file_chk_x": "checkpoint/X_rsp_2", + "file_chk_y": "checkpoint/Y_rsp_2", + "final_prec": 0.001, + "helmholtz_prec": -1.0, + "kain": 3, + "max_iter": 10, + "method": "Hartree-Fock", + "orbital_thrs": 0.01, + "orth_prec": 1e-14, + "property_thrs": -1.0, + "start_prec": 0.001 + } + } + ], + "dynamic": true, + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "reaction_operator": { + "density_type": 1, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": true, + "poisson_prec": 0.001 + } + }, + "frequency": 0.05, + "perturbation": { + "operator": "h_e_dip", + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + }, + "properties": { + "polarizability": { + "pol-0.050000": { + "frequency": 0.05, + "operator": "h_e_dip", + "precision": 0.001, + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + } + } + }, + "unperturbed": { + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "kinetic_operator": { + "derivative": "abgv_55" + }, + "nuclear_operator": { + "nuclear_model": "point_like", + "proj_prec": 0.001, + "shared_memory": false, + "smooth_prec": 0.001 + }, + "reaction_operator": { + "density_type": 0, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": true, + "poisson_prec": 0.001 + } + }, + "localize": false, + "precision": 0.001 + } + } + }, + "scf_calculation": { + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.001 + }, + "kinetic_operator": { + "derivative": "abgv_55" + }, + "nuclear_operator": { + "nuclear_model": "point_like", + "proj_prec": 0.001, + "shared_memory": false, + "smooth_prec": 0.001 + }, + "reaction_operator": { + "density_type": 0, + "dynamic_thrs": true, + "epsilon_dynamic": 2.0, + "epsilon_in": 1.0, + "epsilon_static": 4.0, + "formulation": "exponential", + "kain": 5, + "max_iter": 100, + "nonequilibrium": true, + "poisson_prec": 0.001 + } + }, + "initial_guess": { + "environment": "PCM", + "external_field": "None", + "file_CUBE_a": "cube_vectors/CUBE_a_vector.json", + "file_CUBE_b": "cube_vectors/CUBE_b_vector.json", + "file_CUBE_p": "cube_vectors/CUBE_p_vector.json", + "file_basis": "initial_guess/mrchem.bas", + "file_chk": "checkpoint/phi_scf", + "file_gto_a": "initial_guess/mrchem.moa", + "file_gto_b": "initial_guess/mrchem.mob", + "file_gto_p": "initial_guess/mrchem.mop", + "file_phi_a": "initial_guess/phi_a_scf", + "file_phi_b": "initial_guess/phi_b_scf", + "file_phi_p": "initial_guess/phi_p_scf", + "localize": false, + "method": "Hartree-Fock", + "prec": 0.001, + "relativity": "None", + "restricted": true, + "screen": 12.0, + "type": "sad", + "zeta": 2 + }, + "properties": { + "dipole_moment": { + "dip-1": { + "operator": "h_e_dip", + "precision": 0.001, + "r_O": [ + 0.0, + 0.0, + 0.0 + ] + } + } + } + }, + "schema_name": "mrchem_input", + "schema_version": 1 + }, + "output": { + "properties": { + "center_of_mass": [ + 0.0, + 0.0, + 0.0 + ], + "charge": 0, + "dipole_moment": { + "dip-1": { + "magnitude": 0.001436455029867261, + "r_O": [ + 0.0, + 0.0, + 0.0 + ], + "vector": [ + 0.0, + 0.0, + 0.001436455029867261 + ], + "vector_el": [ + 0.0, + 0.0, + 0.0014364550298491643 + ], + "vector_nuc": [ + 0.0, + 0.0, + 0.0 + ] + } + }, + "geometry": [ + { + "symbol": "H", + "xyz": [ + 0.0, + 0.0, + -0.7001435291738478 + ] + }, + { + "symbol": "H", + "xyz": [ + 0.0, + 0.0, + 0.7001435291738478 + ] + } + ], + "multiplicity": 1, + "orbital_energies": { + "energy": [ + -0.5985110862667959 + ], + "occupation": [ + 2.0 + ], + "spin": [ + "p" + ], + "sum_occupied": -1.1970221725335919 + }, + "polarizability": { + "pol-0.000000": { + "frequency": 0.0, + "isotropic_average": null, + "r_O": [ + 0.0, + 0.0, + 0.0 + ], + "tensor": [ + null, + null, + null, + null, + null, + null, + 0.0, + 0.0, + 9.607625115716449 + ] + }, + "pol-0.050000": { + "frequency": 0.05, + "isotropic_average": null, + "r_O": [ + 0.0, + 0.0, + 0.0 + ], + "tensor": [ + null, + null, + null, + null, + null, + null, + 0.0, + 0.0, + 9.053197959077423 + ] + } + }, + "scf_energy": { + "E_ee": 1.1798069065653403, + "E_eext": 0.0, + "E_el": -1.7997950754459169, + "E_en": -3.2910350186754416, + "E_kin": 0.8884670403157098, + "E_next": 0.0, + "E_nn": 0.7141392859689609, + "E_nuc": 0.7030989718207105, + "E_tot": -1.0966961036252063, + "E_x": -0.5877514493006317, + "E_xc": 0.0, + "Er_el": 0.01071744564910627, + "Er_nuc": -0.01104031414825042, + "Er_tot": -0.0003228684991441501 + } + }, + "provenance": { + "creator": "MRChem", + "mpi_processes": 1, + "nthreads": 1, + "routine": "mrchem.x", + "total_cores": 1, + "version": "1.2.0-alpha" + }, + "rsp_calculations": { + "ext_el-0.000000": { + "components": [ + null, + null, + { + "rsp_solver": { + "converged": true, + "cycles": [ + { + "mo_residual": 1.1840373878577934, + "property_update": -5.230056480276378, + "symmetric_property": -5.230056480276378, + "wall_time": 1.595014405 + }, + { + "mo_residual": 0.5039164949499658, + "property_update": -3.896545272260962, + "symmetric_property": -9.12660175253734, + "wall_time": 10.362402619 + }, + { + "mo_residual": 0.06312576948642024, + "property_update": -0.47827471818918177, + "symmetric_property": -9.604876470726522, + "wall_time": 9.831763404 + }, + { + "mo_residual": 0.01373468159902534, + "property_update": 0.020830591804632093, + "symmetric_property": -9.58404587892189, + "wall_time": 8.421202646 + }, + { + "mo_residual": 0.0033793658719311593, + "property_update": -0.023579236794558867, + "symmetric_property": -9.607625115716449, + "wall_time": 8.485069895 + } + ], + "wall_time": 38.695528888 + } + } + ], + "frequency": 0.0, + "perturbation": "h_e_dip", + "success": true + }, + "ext_el-0.050000": { + "components": [ + null, + null, + { + "rsp_solver": { + "converged": true, + "cycles": [ + { + "mo_residual": 1.7320741496375414, + "property_update": -5.425174312715582, + "symmetric_property": -5.425174312715582, + "wall_time": 3.235084757 + }, + { + "mo_residual": 0.6815580986102592, + "property_update": -3.5252732385174532, + "symmetric_property": -8.950447551233035, + "wall_time": 11.921463333 + }, + { + "mo_residual": 0.0727349479097502, + "property_update": -0.08121816117184899, + "symmetric_property": -9.031665712404884, + "wall_time": 13.243448843 + }, + { + "mo_residual": 0.0204759712344416, + "property_update": -0.03845013785340079, + "symmetric_property": -9.070115850258285, + "wall_time": 12.193242469 + }, + { + "mo_residual": 0.006279219472324256, + "property_update": 0.016917891180861844, + "symmetric_property": -9.053197959077423, + "wall_time": 12.424313607 + } + ], + "wall_time": 53.017622849 + } + } + ], + "frequency": 0.05, + "perturbation": "h_e_dip", + "success": true + } + }, + "scf_calculation": { + "initial_energy": { + "E_ee": 1.1798069065653403, + "E_eext": 0.0, + "E_el": -1.7997950754459169, + "E_en": -3.2910350186754416, + "E_kin": 0.8884670403157098, + "E_next": 0.0, + "E_nn": 0.7141392859689609, + "E_nuc": 0.7030989718207105, + "E_tot": -1.0966961036252063, + "E_x": -0.5877514493006317, + "E_xc": 0.0, + "Er_el": 0.01071744564910627, + "Er_nuc": -0.01104031414825042, + "Er_tot": -0.0003228684991441501 + }, + "success": true + }, + "schema_name": "mrchem_output", + "schema_version": 1, + "success": true + } +} diff --git a/tests/h2_pol_solv/reference/h2.inp b/tests/h2_pol_solv/reference/h2.inp new file mode 100644 index 000000000..e1254f65d --- /dev/null +++ b/tests/h2_pol_solv/reference/h2.inp @@ -0,0 +1,42 @@ +world_prec = 1.0e-4 +world_unit = angstrom + +Molecule { +$coords +H 0.0000 0.0000 -0.3705 +H 0.0000 0.0000 0.3705 +$end +} + +WaveFunction { + method = HF + environment = PCM +} + +PCM { + Permittivity { + epsilon_out { + static = 4.0 + dynamic = 2.0 + nonequilibrium = false + } + } +} + +SCF { + run = true + guess_type = SAD_DZ +} + +Properties { + polarizability = true +} + +Polarizability { + frequency = [0.0] +} + +Response { + run = [false, false, true] +} +~ diff --git a/tests/h2_pol_solv/reference/h2_m.inp b/tests/h2_pol_solv/reference/h2_m.inp new file mode 100644 index 000000000..6d39c74a2 --- /dev/null +++ b/tests/h2_pol_solv/reference/h2_m.inp @@ -0,0 +1,38 @@ +world_prec = 1.0e-4 +world_unit = angstrom + +Molecule { +$coords +H 0.0000 0.0000 -0.3705 +H 0.0000 0.0000 0.3705 +$end +} + +WaveFunction { + method = HF + environment = PCM +} + +PCM { + Permittivity { + epsilon_out { + static = 4.0 + dynamic = 2.0 + nonequilibrium = false + } + } +} + +SCF { + run = true + guess_type = SAD_DZ +} + +ExternalFields { + electric_field = [0.0, 0.0, -0.01] +} + +Properties { + dipole_moment = true +} + diff --git a/tests/h2_pol_solv/reference/h2_p.inp b/tests/h2_pol_solv/reference/h2_p.inp new file mode 100644 index 000000000..595508c4b --- /dev/null +++ b/tests/h2_pol_solv/reference/h2_p.inp @@ -0,0 +1,38 @@ +world_prec = 1.0e-4 +world_unit = angstrom + +Molecule { +$coords +H 0.0000 0.0000 -0.3705 +H 0.0000 0.0000 0.3705 +$end +} + +WaveFunction { + method = HF + environment = PCM +} + +PCM { + Permittivity { + epsilon_out { + static = 4.0 + dynamic = 2.0 + nonequilibrium = false + } + } +} + +SCF { + run = true + guess_type = SAD_DZ +} + +ExternalFields { + electric_field = [0.0, 0.0, 0.01] +} + +Properties { + dipole_moment = true +} + diff --git a/tests/h2_pol_solv/test b/tests/h2_pol_solv/test new file mode 100755 index 000000000..b2d43ad9b --- /dev/null +++ b/tests/h2_pol_solv/test @@ -0,0 +1,31 @@ +#!/usr/bin/env python3 + +import sys +from pathlib import Path +import json + +sys.path.append(str(Path(__file__).resolve().parents[1])) + +from tester import * # isort:skip + +options = script_cli() + +filters = { + SUM_OCCUPIED: rel_tolerance(1.0e-6), + E_KIN: rel_tolerance(1.0e-6), + E_EN: rel_tolerance(1.0e-6), + E_EE: rel_tolerance(1.0e-6), + E_X: rel_tolerance(1.0e-6), + E_XC: rel_tolerance(1.0e-6), + E_EEXT: rel_tolerance(1.0e-6), + E_NEXT: rel_tolerance(1.0e-6), + E_EL: rel_tolerance(1.0e-6), + POLARIZABILITY(0.0): rel_tolerance(1.0e-6), + POLARIZABILITY(0.05): rel_tolerance(1.0e-6), +} + +ierr = 0 +for inp in ["eq", "neq"]: + ierr += run(options, input_file=f"h2-{inp}", filters=filters) + +sys.exit(ierr) diff --git a/tests/li_solv/li.inp b/tests/li_solv/li.inp index fcc9f4217..0b3ad0f49 100644 --- a/tests/li_solv/li.inp +++ b/tests/li_solv/li.inp @@ -23,8 +23,10 @@ }, "Permittivity": { "epsilon_in": 1.0, - "epsilon_out": 2.0, - "formulation": "exponential" + "formulation": "exponential", + "epsilon_out": { + "static": 2.0 + } } }, "SCF": { diff --git a/tests/li_solv/reference/li.json b/tests/li_solv/reference/li.json index fde26d561..6c6b2ced4 100644 --- a/tests/li_solv/reference/li.json +++ b/tests/li_solv/reference/li.json @@ -94,8 +94,9 @@ "reaction_operator": { "density_type": "total", "dynamic_thrs": false, + "epsilon_dynamic": 2.0, "epsilon_in": 1.0, - "epsilon_out": 2.0, + "epsilon_static": 2.0, "formulation": "exponential", "kain": 5, "max_iter": 100, @@ -165,7 +166,7 @@ "charge": 1, "dipole_moment": { "dip-1": { - "magnitude": 7.400866885254787e-14, + "magnitude": 2.7058089635853275e-13, "r_O": [ 0.0, 0.0, @@ -201,7 +202,7 @@ "multiplicity": 1, "orbital_energies": { "energy": [ - -2.204016950767741 + -2.204016950767725 ], "occupation": [ 2.0 @@ -209,23 +210,23 @@ "spin": [ "p" ], - "sum_occupied": -4.408033901535482 + "sum_occupied": -4.40803390153545 }, "scf_energy": { - "E_ee": 3.3327993838943106, + "E_ee": 3.332799383894344, "E_eext": 0.0, - "E_el": -7.094572904570907, - "E_en": -16.530755287695293, - "E_kin": 7.675564301900204, + "E_el": -7.094572904570924, + "E_en": -16.53075528769541, + "E_kin": 7.67556430190028, "E_next": 0.0, "E_nn": 0.0, - "E_nuc": -0.19217053480332283, - "E_tot": -7.28674343937423, - "E_x": -0.4164844953175972, - "E_xc": -1.2837729977609416, - "Er_el": 0.12807619040840895, - "Er_nuc": -0.19217053480332283, - "Er_tot": -0.06409434439491388 + "E_nuc": -0.19217053480331467, + "E_tot": -7.286743439374239, + "E_x": -0.4164844953176016, + "E_xc": -1.283772997760944, + "Er_el": 0.12807619040840695, + "Er_nuc": -0.19217053480331467, + "Er_tot": -0.06409434439490772 } }, "provenance": { @@ -239,20 +240,20 @@ "rsp_calculations": null, "scf_calculation": { "initial_energy": { - "E_ee": 3.3327993838943106, + "E_ee": 3.332799383894344, "E_eext": 0.0, - "E_el": -7.094572904570907, - "E_en": -16.530755287695293, - "E_kin": 7.675564301900204, + "E_el": -7.094572904570924, + "E_en": -16.53075528769541, + "E_kin": 7.67556430190028, "E_next": 0.0, "E_nn": 0.0, - "E_nuc": -0.19217053480332283, - "E_tot": -7.28674343937423, - "E_x": -0.4164844953175972, - "E_xc": -1.2837729977609416, - "Er_el": 0.12807619040840895, - "Er_nuc": -0.19217053480332283, - "Er_tot": -0.06409434439491388 + "E_nuc": -0.19217053480331467, + "E_tot": -7.286743439374239, + "E_x": -0.4164844953176016, + "E_xc": -1.283772997760944, + "Er_el": 0.12807619040840695, + "Er_nuc": -0.19217053480331467, + "Er_tot": -0.06409434439490772 }, "success": true }, diff --git a/tests/li_solv/reference/li.out b/tests/li_solv/reference/li.out index 4f21e9988..0085b9fd4 100644 --- a/tests/li_solv/reference/li.out +++ b/tests/li_solv/reference/li.out @@ -11,10 +11,10 @@ *** *** *** VERSION 1.2.0-alpha *** *** *** -*** Git branch vecmult *** -*** Git commit hash 07cbcd62711b3df44973-dirty *** -*** Git commit author Gabriel Gerez *** -*** Git commit date Mon Oct 30 13:15:37 2023 +0100 *** +*** Git branch pcm-terms-response *** +*** Git commit hash a1cfb2ea38e3362137ec-dirty *** +*** Git commit author Roberto Di Remigio Eikås *** +*** Git commit date Fri Nov 3 11:53:50 2023 +0100 *** *** *** *** Contact: luca.frediani@uit.no *** *** *** @@ -53,7 +53,7 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s Git commit date : Sat Oct 21 13:32:25 2023 +0200 Linear algebra : EIGEN v3.4.0 - Parallelization : OpenMP (1 threads) + Parallelization : NONE --------------------------------------------------------------------------- @@ -99,10 +99,6 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s =========================================================================== Solvation Cavity ---------------------------------------------------------------------------- - Formulation : exponential - Dielectric constant : (in) 1.000000 - : (out) 2.000000 --------------------------------------------------------------------------- N R_0 Alpha Beta Sigma Radius x y z --------------------------------------------------------------------------- @@ -110,6 +106,15 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s =========================================================================== +=========================================================================== + Solvation Permittivity +--------------------------------------------------------------------------- + Formulation : exponential + Dielectric constant : (in) 1.000000 + : (out) 2.000000 +=========================================================================== + + *************************************************************************** *** *** @@ -271,7 +276,7 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s *** *** *** Exiting MRChem *** *** *** -*** Wall time : 0h 0m 53s *** +*** Wall time : 0h 1m 13s *** *** *** *************************************************************************** diff --git a/tests/solventeffect/reaction_operator.cpp b/tests/solventeffect/reaction_operator.cpp index d381e46c5..40dfdda49 100644 --- a/tests/solventeffect/reaction_operator.cpp +++ b/tests/solventeffect/reaction_operator.cpp @@ -90,7 +90,7 @@ TEST_CASE("ReactionOperator", "[reaction_operator]") { auto rho_nuc = chemistry::compute_nuclear_density(prec, molecule, 100); int kain = 4; - auto scrf_p = std::make_unique(dielectric_func, rho_nuc, P_p, D_p, kain, 100, false, "total"); + auto scrf_p = std::make_unique(dielectric_func, rho_nuc, P_p, D_p, kain, 100, false, SCRFDensityType::TOTAL); auto Reo = std::make_shared(std::move(scrf_p), Phi_p); Reo->setup(prec);