Linear response with solvent (#471)

* Use C++17

* Use std::tie to unpack tuple return value

* Remove some cruft from Cavity and Permittivity

* SCRF deals in densities only, not orbitals

* Refactor Reaction{Potential,Operator} following Coulomb example

In preparation for the computation of the response terms, we make ReactionPotential abstract and add
subclass ReactionPotentialD1 to deal with one single set of orbitals.

* Modify input to accept static and dynamic permittivities

* Handle nonequilibrium solvation more gracefully

* Sign change was necessary

* Update input parser

* Small fixes to include order

* Rename section for permittivity outside

* Add tests

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)

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`

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
 
 

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*

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*

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'},

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

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`

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

python/template.yml

@@ -936,20 +936,40 @@ 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
             predicates:
               - 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:

src/driver.cpp

@@ -104,7 +104,7 @@ namespace driver {
 DerivativeOperator_p get_derivative(const std::string &name);
 template <int I> RankOneOperator<I> get_operator(const std::string &name, const json &json_oper);
 template <int I, int J> RankTwoOperator<I, J> 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<PoissonOperator>(*MRA, poisson_prec);
         auto D_p = std::make_shared<mrcpp::ABGVOperator<3>>(*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<SCRF>(dielectric_func, rho_nuc, P_p, D_p, kain, max_iter, dynamic_thrs, density_type);
-        auto V_R = std::make_shared<ReactionOperator>(std::move(scrf_p), Phi_p);
+
+        // initialize reaction potential object
+        auto V_R = [&] {
+            if (order == 0) {
+                return std::make_shared<ReactionOperator>(std::move(scrf_p), Phi_p);
+            } else {
+                return std::make_shared<ReactionOperator>(std::move(scrf_p), Phi_p, X_p, Y_p);
+            }
+        }();
+
+        // set reaction potential in FockBuilder object
         F.getReactionOperator() = V_R;
     }
     ///////////////////////////////////////////////////////////