diff --git a/.circleci/config.yml b/.circleci/config.yml index f21425fe6..d44d3263d 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -1,9 +1,9 @@ version: 2 variables: - ubuntu-2004: &ubuntu-2004 + ubuntu-2204: &ubuntu-2204 docker: - - image: ghcr.io/mrchemsoft/metamr/circleci_ubuntu-20.04:sha-343e011 + - image: ghcr.io/mrchemsoft/metamr/circleci_ubuntu-22.04:sha-9f6ecd4 name: tsubame user: merzbow working_directory: ~/mrchem @@ -52,14 +52,14 @@ variables: jobs: serial-py3: - <<: *ubuntu-2004 + <<: *ubuntu-2204 steps: - checkout - *configure-serial - *build - *tests omp-py3: - <<: *ubuntu-2004 + <<: *ubuntu-2204 environment: - OMP_NUM_THREADS: '2' steps: @@ -68,7 +68,7 @@ jobs: - *build - *tests mpi-py3: - <<: *ubuntu-2004 + <<: *ubuntu-2204 steps: - checkout - *configure-mpi diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index cbafdc6d8..1198d5c85 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -42,19 +42,19 @@ jobs: activate-environment: mrchem-gha environment-file: .github/mrchem-gha.yml channel-priority: true - python-version: 3.6 + python-version: 3.9 use-only-tar-bz2: true # IMPORTANT: This needs to be set for caching to work properly! - name: Configure shell: bash -l {0} run: | python ./setup --type=$BUILD_TYPE --omp --arch-flags=false --generator=Ninja --prefix=$GITHUB_WORKSPACE/Software/MRChem build - + - name: Build shell: bash -l {0} run: | cmake --build build --config $BUILD_TYPE --target install -- -v -d stats - + - name: Test shell: bash -l {0} run: | diff --git a/.github/workflows/code-coverage.yml b/.github/workflows/code-coverage.yml index 4efe22654..7dbdec202 100644 --- a/.github/workflows/code-coverage.yml +++ b/.github/workflows/code-coverage.yml @@ -26,18 +26,18 @@ jobs: activate-environment: mrchem-codecov environment-file: .github/mrchem-codecov.yml channel-priority: true - python-version: 3.6 + python-version: 3.9 - name: Configure shell: bash -l {0} run: | python ./setup --type=$BUILD_TYPE --arch-flags=false --coverage --generator=Ninja --prefix=$GITHUB_WORKSPACE/Software/MRChem build - + - name: Build shell: bash -l {0} run: | cmake --build build --config $BUILD_TYPE --target install -- -v -d stats - + - name: Test MRChem and generate coverage report shell: bash -l {0} run: | diff --git a/doc/index.rst b/doc/index.rst index 004f30cd6..400ccc45c 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -42,6 +42,8 @@ Features in MRChem-1.1: - Electric field + Solvent effects - Cavity-free PCM + - Poisson-Boltzmann PCM + - Linearized Poisson-Boltzmann PCM * Properties: + Ground state energy + Dipole moment diff --git a/doc/installation.rst b/doc/installation.rst index 1c267a197..83f07899b 100644 --- a/doc/installation.rst +++ b/doc/installation.rst @@ -6,7 +6,7 @@ Installation Build prerequisites ------------------- -- Python-3.7 (or later) +- Python-3.9 (or later) - CMake-3.14 (or later) - GNU-5.4 or Intel-17 (or later) compilers (C++14 standard) diff --git a/doc/programmers/code_reference/environment.rst b/doc/programmers/code_reference/environment.rst index 74f69f74d..7eb0c3381 100644 --- a/doc/programmers/code_reference/environment.rst +++ b/doc/programmers/code_reference/environment.rst @@ -21,11 +21,38 @@ Permittivity :protected-members: :private-members: -SCRF +DHScreening ------------ -.. doxygenclass:: mrchem::SCRF +.. doxygenclass:: mrchem::DHScreening :project: MRChem :members: :protected-members: :private-members: + +GPESolver +------------ + +.. doxygenclass:: mrchem::GPESolver + :project: MRChem + :members: + :protected-members: + :private-members: + +PBESolver +------------ + +.. doxygenclass:: mrchem::PBESolver + :project: MRChem + :members: + :protected-members: + :private-members: + +LPBESolver +------------ + +.. doxygenclass:: mrchem::LPBESolver + :project: MRChem + :members: + :protected-members: + :private-members: \ No newline at end of file diff --git a/doc/programming.bib b/doc/programming.bib index 010b47457..e284b31b5 100644 --- a/doc/programming.bib +++ b/doc/programming.bib @@ -333,7 +333,6 @@ @article{Fosso-Tande2013 abstract = {We describe and present results of the implementation of the surface and volume polarization for electrostatics (SVPE) solvation model. Unlike most other implementations of the solvation model where the solute and the solvent are described with multiple numerical representations, our implementation uses a multiresolution, adaptive multiwavelet basis to describe both the solute and the solvent. This requires reformulation to use integral equations throughout as well as a conscious management of numerical properties of the basis. {\textcopyright} 2013 Elsevier B.V. All rights reserved.}, author = {Fosso-Tande, Jacob and Harrison, Robert J.}, doi = {10.1016/j.cplett.2013.01.065}, -file = {:home/ggerez/.local/share/data/Mendeley Ltd./Mendeley Desktop/Downloaded/Fosso-Tande, Harrison - 2013 - Implicit solvation models in a multiresolution multiwavelet basis.pdf:pdf}, issn = {00092614}, journal = {Chem. Phys. Lett.}, pages = {179--184}, @@ -343,3 +342,20 @@ @article{Fosso-Tande2013 year = {2013} } +@article{gerez2023, +author = {Gerez S, Gabriel A. and Di Remigio Eikås, Roberto and Jensen, Stig Rune and Bjørgve, Magnar and Frediani, Luca}, +title = {Cavity-Free Continuum Solvation: Implementation and Parametrization in a Multiwavelet Framework}, +journal = {Journal of Chemical Theory and Computation}, +volume = {19}, +number = {7}, +pages = {1986-1997}, +year = {2023}, +doi = {10.1021/acs.jctc.2c01098}, + note ={PMID: 36933225}, +URL = { + https://doi.org/10.1021/acs.jctc.2c01098 +}, +eprint = { + https://doi.org/10.1021/acs.jctc.2c01098 +} +} diff --git a/doc/users/user_inp.rst b/doc/users/user_inp.rst index c3a02d239..711a092fd 100644 --- a/doc/users/user_inp.rst +++ b/doc/users/user_inp.rst @@ -194,7 +194,7 @@ as the ``world_origin`` is the true origin). WaveFunction ------------ -Here we give the wavefunction method and whether we run spin restricted (alpha +Here we give the wavefunction method, environment used (for solvent models) and whether we run spin restricted (alpha and beta spins are forced to occupy the same spatial orbitals) or not (method must be specified, otherwise defaults are shown): @@ -203,6 +203,7 @@ must be specified, otherwise defaults are shown): WaveFunction { method = # Core, Hartree, HF or DFT restricted = true # Spin restricted/unrestricted + environment = pcm # Environment (pcm, pcm-pb, pcm-lpb) defaults to none } There are currently four methods available: Core Hamiltonian, Hartree, @@ -212,6 +213,11 @@ B3LYP``), *or* you can set ``method = DFT`` and specify a "non-standard" functional in the separate DFT section (see below). See :ref:`User input reference` for a list of available default functionals. +The solvent model implemented is a cavity free PCM, described in :cite:`gerez2023`. +In this model we have implemented the Generalized Poisson equation solver, keyword ``pcm``, a +Poisson-Boltzmann solver, keyword ``pcm-pb`` and a Linearized Poisson-Boltzmann solver, keyword ``pcm-lpb``. +Further details for the calculation have to be included in the ``PCM`` section, see :ref: `User input reference` for details. + .. note:: Restricted open-shell wavefunctions are not supported. diff --git a/doc/users/user_ref.rst b/doc/users/user_ref.rst index 41fc52b62..75ac097cd 100644 --- a/doc/users/user_ref.rst +++ b/doc/users/user_ref.rst @@ -329,14 +329,14 @@ User input reference **Predicates** - ``value.lower() in ['none', 'zora', 'nzora']`` - :environment: Set method for treatment of environment. ``none`` for vacuum calculation. ``PCM`` for Polarizable Continuum Model, which will activate the ``PCM`` input section for further parametrization options. + :environment: Set method for treatment of environment. ``none`` for vacuum calculation. ``PCM`` for Polarizable Continuum Model, which will activate the ``PCM`` input section for further parametrization options. The ``PB`` and ``LPB`` variants add the Poisson-Boltzmann and Linearized Poisson-Boltzmann solvers, respectively. **Type** ``str`` **Default** ``none`` **Predicates** - - ``value.lower() in ['none', 'pcm']`` + - ``value.lower() in ['none', 'pcm', 'pcm_pb', 'pcm_lpb']`` :nuclear_model: Type of nucleus model. Point-like (numerical smoothing): HFYGB (default), parabola or minimal. Finite models (physical smoothing): Gaussian or Homogeneous sphere Finite models are derived from nuclear RMS radius, Visscher (1997) @@ -860,6 +860,79 @@ User input reference **Default** ``user['SCF']['kain']`` + :Solvent: Parameters for the Self-Consistent Reaction Field optimization. + + :red:`Sections` + :Permittivity: Parameters for the permittivity function. + + :red:`Keywords` + :epsilon_in: Permittivity inside the cavity. 1.0 is the permittivity of free space, anything other than this is undefined behaviour. + + **Type** ``float`` + + **Default** ``1.0`` + + :formulation: Formulation of the Permittivity function. Currently only the exponential is available. + + **Type** ``str`` + + **Default** ``exponential`` + + **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']['Solvent']['Permittivity']['epsilon_out']['static']`` + + :DebyeHuckelScreening: Parameters for the Debye-Huckel screening factor + + :red:`Keywords` + :ion_strength: Ionic strength of the electrolyte in mol/L. This represents the concentration of the ions in the bulk solvent. + + **Type** ``float`` + + **Default** ``1.0`` + + :ion_radius: Amount with which the vdw-radius of the atoms will be increased. The screening factor will have an area of effect that is often going to be larger than the vdw-cavity, but centered in the same atoms. + + **Type** ``float`` + + **Default** ``0.0`` + + :ion_width: Width of the transition between the solute and the ion accessible part. + + **Type** ``float`` + + **Default** ``0.2`` + + :formulation: formulation of the debye-huckel screening factor. Currently only the variable factor is implemented. ``variable``: implement the screening functions as k = (1-C_ion)k_out + + **Type** ``str`` + + **Default** ``variable`` + + **Predicates** + - ``value.lower() in ['variable']`` + :Cavity: Define the interlocking spheres cavity. :red:`Keywords` @@ -896,47 +969,7 @@ User input reference **Default** ``0.2`` - :Permittivity: Parameters for the permittivity function. - - :red:`Keywords` - :epsilon_in: Permittivity inside the cavity. 1.0 is the permittivity of free space, anything other than this is undefined behaviour. - - **Type** ``float`` - - **Default** ``1.0`` - - :formulation: Formulation of the Permittivity function. Currently only the exponential is available. - - **Type** ``str`` - - **Default** ``exponential`` - - **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']`` - - :GeometryOptimizer: Includes parameters related to the internal geometry optimization using the SQNM (Stabilized Quasi-Newton Method) for noisy PES. + :GeometryOptimizer: Includes parameters related to the internal geometry optimization using the SQNM (Stabilized Quasi-Newton Method) for noisy PES. Geometry optimizations require accurate forces. Consider setting world_prec to 1e-5 to 1e-7. Convergence issues can usually be solved by increasing the precision of the SCF calculation. If that does not work, try setting the initial step size manually. :red:`Keywords` :run: Run optimizer. Otherwise single point energy/properties are computed. @@ -951,37 +984,37 @@ User input reference **Default** ``False`` - :init_step_size: Initial step size. + :init_step_size: Initial step size. For systems with hard bonds (e.g. C-C) use a value between and 1.0 and 2.5. If a system only contains weaker bonds a value up to 5.0 may speed up the convergence. Use a small negative number (should be between -0.1 and -0.5) for an automatic guess. The optimal step size is the inverse of the largest eigenvalue of the hessian matrix. **Type** ``float`` **Default** ``-0.5`` - :minimal_step_size: Minimal step size. + :minimal_step_size: Minimal step size. It rarely makes sense to change it. **Type** ``float`` **Default** ``0.01`` - :max_history_length: Maximum length of history. + :max_history_length: Maximum length of history list. Energies and forces from the previous n geometry optimization iterations are used to estimate the hessian matrix. Use a value between 2 and 20. A lower value makes the SQNM algorithm behave more like steepest descent and slows down convergence. But it can handle more noise in the energies and forces. It rarely makes sense to change it. **Type** ``int`` **Default** ``10`` - :subspace_tolerance: Subspace tolerance. + :subspace_tolerance: Lower limit on linear dependencies of basis vectors in history listSubspace tolerance. Use a number between 1e-9 and 1e-1. A high subspace tolerance slows down convergence but improves numerical stability when the energies and forces contain a lot of noise. It rarely makes sense to change it. **Type** ``float`` **Default** ``0.001`` - :max_iter: Maximum number of iterations. + :max_iter: Maximum number of geometry optimization iterations. **Type** ``int`` **Default** ``100`` - :max_force_component: Maximum force component. + :max_force_component: The geometry optimization stopps when the absolute value of all force components is smaller than this keyword. A value between 1e-3 and 1e-4 is tight enough for most applications. **Type** ``float`` @@ -1050,6 +1083,36 @@ User input reference **Default** ``2.5417464739297717`` + :boltzmann_constant: | Boltzmann constant in (unit: J K^-1). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``1.380649e-23`` + + :elementary_charge: | Elementary charge in (unit: C). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``1.602176634e-19`` + + :e0: | Permittivity of free space (unit: F m^-1). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``8.8541878128e-12`` + + :N_a: | Avogadro constant (unit: mol^-1). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``6.02214076e+23`` + + :meter2bohr: | conversion factor from meter to Bohr radius (unit: m^-1). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``18897261246.2577`` + :Elements: list of elements with data :red:`Sections` diff --git a/python/mrchem/api.py b/python/mrchem/api.py index 0ac8c2c23..f0cad11e1 100644 --- a/python/mrchem/api.py +++ b/python/mrchem/api.py @@ -25,16 +25,11 @@ import math -from .helpers import ( - write_scf_fock, - write_scf_guess, - write_scf_solver, - write_scf_properties, - write_scf_plot, - write_rsp_calc, - parse_wf_method, -) -from .periodictable import PeriodicTable as PT, PeriodicTableByZ as PT_Z +from .helpers import (parse_wf_method, write_rsp_calc, write_scf_fock, + write_scf_guess, write_scf_plot, write_scf_properties, + write_scf_solver) +from .periodictable import PeriodicTable as PT +from .periodictable import PeriodicTableByZ as PT_Z from .validators import MoleculeValidator @@ -128,7 +123,8 @@ def write_molecule(user_dict, origin): "charge": mol.charge, "coords": mol.get_coords_in_program_syntax(), } - if user_dict["WaveFunction"]["environment"].lower() == "pcm": + + if "pcm" in user_dict["WaveFunction"]["environment"].lower(): mol_dict["cavity"] = { "spheres": mol.get_cavity_in_program_syntax(), } diff --git a/python/mrchem/helpers.py b/python/mrchem/helpers.py index 616c7dd3c..cff57844e 100644 --- a/python/mrchem/helpers.py +++ b/python/mrchem/helpers.py @@ -23,6 +23,7 @@ # # +from math import sqrt from pathlib import Path from .CUBEparser import parse_files @@ -72,7 +73,7 @@ def write_scf_fock(user_dict, wf_dict, origin): } # Reaction - if user_dict["WaveFunction"]["environment"].lower() == "pcm": + if user_dict["WaveFunction"]["environment"].lower() != "none": fock_dict["reaction_operator"] = _reaction_operator_handler(user_dict) # Coulomb @@ -128,22 +129,43 @@ def _reaction_operator_handler(user_dict, rsp=False): else: density_type = 2 - return { + # reaction field operator settings common to all continuum models + reo_dict = { + "solver_type": "Generalized_Poisson", "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"][ + "epsilon_in": user_dict["PCM"]["Solvent"]["Permittivity"]["epsilon_in"], + "epsilon_static": user_dict["PCM"]["Solvent"]["Permittivity"]["epsilon_out"]["static"], + "epsilon_dynamic": user_dict["PCM"]["Solvent"]["Permittivity"]["epsilon_out"]["dynamic"], + "nonequilibrium": user_dict["PCM"]["Solvent"]["Permittivity"]["epsilon_out"][ "nonequilibrium" ], - "formulation": user_dict["PCM"]["Permittivity"]["formulation"], + "formulation": user_dict["PCM"]["Solvent"]["Permittivity"]["formulation"], + "kappa_out": 0.0, + "ion_radius": user_dict["PCM"]["Solvent"]["DebyeHuckelScreening"]["ion_radius"], + "ion_width": user_dict["PCM"]["Solvent"]["DebyeHuckelScreening"]["ion_width"], + "DHS-formulation": user_dict["PCM"]["Solvent"]["DebyeHuckelScreening"]["formulation"], } + # ionic solvent continuum model + ionic_model = user_dict["WaveFunction"]["environment"].lower().split("_")[-1] + if ionic_model in ("pb", "lpb"): + permittivity = user_dict["PCM"]["Solvent"]["Permittivity"]["epsilon_out"]["static"] + ionic_strength = user_dict["PCM"]["Solvent"]["DebyeHuckelScreening"]["ion_strength"] + kappa_out = compute_kappa(user_dict["Constants"], permittivity, ionic_strength) + reo_dict |= { + "kappa_out": kappa_out, + "solver_type": "Poisson-Boltzmann" + if ionic_model == "pb" + else "Linearized_Poisson-Boltzmann", + } + + return reo_dict + def write_scf_guess(user_dict, wf_dict): guess_str = user_dict["SCF"]["guess_type"].lower() @@ -422,7 +444,7 @@ def write_rsp_fock(user_dict, wf_dict): } # Reaction - if user_dict["WaveFunction"]["environment"].lower() == "pcm": + if user_dict["WaveFunction"]["environment"].lower() != "none": fock_dict["reaction_operator"] = _reaction_operator_handler(user_dict, rsp=True) return fock_dict @@ -547,3 +569,19 @@ def parse_wf_method(user_dict): "dft_funcs": dft_funcs, } return wf_dict + + +def compute_kappa(constants, eps, I): + kb = constants["boltzmann_constant"] + e = constants["elementary_charge"] + e_0 = constants["e0"] + N_a = constants["N_a"] + m2au = constants["meter2bohr"] + T = 298.15 + + numerator = e_0 * eps * kb * T + denominator = 2.0 * (e**2) * N_a * 1000.0 * I + + debye_length = sqrt(numerator / denominator) * m2au + + return 1.0 / debye_length diff --git a/python/mrchem/input_parser/api.py b/python/mrchem/input_parser/api.py index a9abc9a33..3be32d7d7 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 2024-01-29 +# This file was automatically generated by parselglossy on 2024-02-08 # Editing is *STRONGLY DISCOURAGED* from copy import deepcopy @@ -300,7 +300,9 @@ def stencil() -> JSONDict: 'name': 'environment', 'predicates': [ 'value.lower() ' "in ['none', " - "'pcm']"], + "'pcm', " + "'pcm_pb', " + "'pcm_lpb']"], 'type': 'str'}, { 'default': 'point_like', 'name': 'nuclear_model', @@ -584,6 +586,43 @@ def stencil() -> JSONDict: 'name': 'kain', 'type': 'int'}], 'name': 'SCRF'}, + { 'name': 'Solvent', + 'sections': [ { 'keywords': [ { 'default': 1.0, + 'name': 'epsilon_in', + 'type': 'float'}, + { 'default': 'exponential', + 'name': 'formulation', + 'predicates': [ 'value.lower() ' + 'in ' + "['exponential']"], + 'type': 'str'}], + 'name': 'Permittivity', + 'sections': [ { 'keywords': [ { 'default': False, + 'name': 'nonequilibrium', + 'type': 'bool'}, + { 'default': 1.0, + 'name': 'static', + 'type': 'float'}, + { 'default': "user['PCM']['Solvent']['Permittivity']['epsilon_out']['static']", + 'name': 'dynamic', + 'type': 'float'}], + 'name': 'epsilon_out'}]}, + { 'keywords': [ { 'default': 1.0, + 'name': 'ion_strength', + 'type': 'float'}, + { 'default': 0.0, + 'name': 'ion_radius', + 'type': 'float'}, + { 'default': 0.2, + 'name': 'ion_width', + 'type': 'float'}, + { 'default': 'variable', + 'name': 'formulation', + 'predicates': [ 'value.lower() ' + 'in ' + "['variable']"], + 'type': 'str'}], + 'name': 'DebyeHuckelScreening'}]}, { 'keywords': [ { 'default': 'atoms', 'name': 'mode', 'predicates': [ 'value.lower() ' @@ -603,27 +642,7 @@ def stencil() -> JSONDict: { 'default': 0.2, 'name': 'sigma', 'type': 'float'}], - 'name': 'Cavity'}, - { 'keywords': [ { 'default': 1.0, - 'name': 'epsilon_in', - 'type': 'float'}, - { 'default': 'exponential', - 'name': 'formulation', - 'predicates': [ 'value.lower() ' - 'in ' - "['exponential']"], - 'type': 'str'}], - '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'}]}]}, + 'name': 'Cavity'}]}, { 'keywords': [ { 'default': False, 'name': 'run', 'type': 'bool'}, @@ -678,6 +697,21 @@ def stencil() -> JSONDict: 'type': 'float'}, { 'default': 2.5417464739297717, 'name': 'dipmom_au2debye', + 'type': 'float'}, + { 'default': 1.380649e-23, + 'name': 'boltzmann_constant', + 'type': 'float'}, + { 'default': 1.602176634e-19, + 'name': 'elementary_charge', + 'type': 'float'}, + { 'default': 8.8541878128e-12, + 'name': 'e0', + 'type': 'float'}, + { 'default': 6.02214076e+23, + 'name': 'N_a', + 'type': 'float'}, + { 'default': 18897261246.2577, + 'name': 'meter2bohr', 'type': 'float'}], 'name': 'Constants'}, { 'name': 'Elements', diff --git a/python/mrchem/input_parser/docs/user_ref.rst b/python/mrchem/input_parser/docs/user_ref.rst index 41fc52b62..75ac097cd 100644 --- a/python/mrchem/input_parser/docs/user_ref.rst +++ b/python/mrchem/input_parser/docs/user_ref.rst @@ -329,14 +329,14 @@ User input reference **Predicates** - ``value.lower() in ['none', 'zora', 'nzora']`` - :environment: Set method for treatment of environment. ``none`` for vacuum calculation. ``PCM`` for Polarizable Continuum Model, which will activate the ``PCM`` input section for further parametrization options. + :environment: Set method for treatment of environment. ``none`` for vacuum calculation. ``PCM`` for Polarizable Continuum Model, which will activate the ``PCM`` input section for further parametrization options. The ``PB`` and ``LPB`` variants add the Poisson-Boltzmann and Linearized Poisson-Boltzmann solvers, respectively. **Type** ``str`` **Default** ``none`` **Predicates** - - ``value.lower() in ['none', 'pcm']`` + - ``value.lower() in ['none', 'pcm', 'pcm_pb', 'pcm_lpb']`` :nuclear_model: Type of nucleus model. Point-like (numerical smoothing): HFYGB (default), parabola or minimal. Finite models (physical smoothing): Gaussian or Homogeneous sphere Finite models are derived from nuclear RMS radius, Visscher (1997) @@ -860,6 +860,79 @@ User input reference **Default** ``user['SCF']['kain']`` + :Solvent: Parameters for the Self-Consistent Reaction Field optimization. + + :red:`Sections` + :Permittivity: Parameters for the permittivity function. + + :red:`Keywords` + :epsilon_in: Permittivity inside the cavity. 1.0 is the permittivity of free space, anything other than this is undefined behaviour. + + **Type** ``float`` + + **Default** ``1.0`` + + :formulation: Formulation of the Permittivity function. Currently only the exponential is available. + + **Type** ``str`` + + **Default** ``exponential`` + + **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']['Solvent']['Permittivity']['epsilon_out']['static']`` + + :DebyeHuckelScreening: Parameters for the Debye-Huckel screening factor + + :red:`Keywords` + :ion_strength: Ionic strength of the electrolyte in mol/L. This represents the concentration of the ions in the bulk solvent. + + **Type** ``float`` + + **Default** ``1.0`` + + :ion_radius: Amount with which the vdw-radius of the atoms will be increased. The screening factor will have an area of effect that is often going to be larger than the vdw-cavity, but centered in the same atoms. + + **Type** ``float`` + + **Default** ``0.0`` + + :ion_width: Width of the transition between the solute and the ion accessible part. + + **Type** ``float`` + + **Default** ``0.2`` + + :formulation: formulation of the debye-huckel screening factor. Currently only the variable factor is implemented. ``variable``: implement the screening functions as k = (1-C_ion)k_out + + **Type** ``str`` + + **Default** ``variable`` + + **Predicates** + - ``value.lower() in ['variable']`` + :Cavity: Define the interlocking spheres cavity. :red:`Keywords` @@ -896,47 +969,7 @@ User input reference **Default** ``0.2`` - :Permittivity: Parameters for the permittivity function. - - :red:`Keywords` - :epsilon_in: Permittivity inside the cavity. 1.0 is the permittivity of free space, anything other than this is undefined behaviour. - - **Type** ``float`` - - **Default** ``1.0`` - - :formulation: Formulation of the Permittivity function. Currently only the exponential is available. - - **Type** ``str`` - - **Default** ``exponential`` - - **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']`` - - :GeometryOptimizer: Includes parameters related to the internal geometry optimization using the SQNM (Stabilized Quasi-Newton Method) for noisy PES. + :GeometryOptimizer: Includes parameters related to the internal geometry optimization using the SQNM (Stabilized Quasi-Newton Method) for noisy PES. Geometry optimizations require accurate forces. Consider setting world_prec to 1e-5 to 1e-7. Convergence issues can usually be solved by increasing the precision of the SCF calculation. If that does not work, try setting the initial step size manually. :red:`Keywords` :run: Run optimizer. Otherwise single point energy/properties are computed. @@ -951,37 +984,37 @@ User input reference **Default** ``False`` - :init_step_size: Initial step size. + :init_step_size: Initial step size. For systems with hard bonds (e.g. C-C) use a value between and 1.0 and 2.5. If a system only contains weaker bonds a value up to 5.0 may speed up the convergence. Use a small negative number (should be between -0.1 and -0.5) for an automatic guess. The optimal step size is the inverse of the largest eigenvalue of the hessian matrix. **Type** ``float`` **Default** ``-0.5`` - :minimal_step_size: Minimal step size. + :minimal_step_size: Minimal step size. It rarely makes sense to change it. **Type** ``float`` **Default** ``0.01`` - :max_history_length: Maximum length of history. + :max_history_length: Maximum length of history list. Energies and forces from the previous n geometry optimization iterations are used to estimate the hessian matrix. Use a value between 2 and 20. A lower value makes the SQNM algorithm behave more like steepest descent and slows down convergence. But it can handle more noise in the energies and forces. It rarely makes sense to change it. **Type** ``int`` **Default** ``10`` - :subspace_tolerance: Subspace tolerance. + :subspace_tolerance: Lower limit on linear dependencies of basis vectors in history listSubspace tolerance. Use a number between 1e-9 and 1e-1. A high subspace tolerance slows down convergence but improves numerical stability when the energies and forces contain a lot of noise. It rarely makes sense to change it. **Type** ``float`` **Default** ``0.001`` - :max_iter: Maximum number of iterations. + :max_iter: Maximum number of geometry optimization iterations. **Type** ``int`` **Default** ``100`` - :max_force_component: Maximum force component. + :max_force_component: The geometry optimization stopps when the absolute value of all force components is smaller than this keyword. A value between 1e-3 and 1e-4 is tight enough for most applications. **Type** ``float`` @@ -1050,6 +1083,36 @@ User input reference **Default** ``2.5417464739297717`` + :boltzmann_constant: | Boltzmann constant in (unit: J K^-1). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``1.380649e-23`` + + :elementary_charge: | Elementary charge in (unit: C). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``1.602176634e-19`` + + :e0: | Permittivity of free space (unit: F m^-1). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``8.8541878128e-12`` + + :N_a: | Avogadro constant (unit: mol^-1). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``6.02214076e+23`` + + :meter2bohr: | conversion factor from meter to Bohr radius (unit: m^-1). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation. + + **Type** ``float`` + + **Default** ``18897261246.2577`` + :Elements: list of elements with data :red:`Sections` diff --git a/python/mrchem/physical_constants.py b/python/mrchem/physical_constants.py index 7a150db61..e96c6b4ff 100644 --- a/python/mrchem/physical_constants.py +++ b/python/mrchem/physical_constants.py @@ -131,6 +131,41 @@ def __init__(self): docstring = f"| Conversion factor for dipoles from atomic units to Debye (unit: {unit}). Affected code: Printed value of dipole moments." self.add_constant(name, unit, value, docstring) + # Boltzmann constant in J K^{-1} + name = "boltzmann_constant" + unit = "J K^-1" + value = self.qce.kb + docstring = f"| Boltzmann constant in (unit: {unit}). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation." + self.add_constant(name, unit, value, docstring) + + # elementary charge in C + name = "elementary_charge" + unit = "C" + value = self.qce.get("elementary charge") + docstring = f"| Elementary charge in (unit: {unit}). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation." + self.add_constant(name, unit, value, docstring) + + # Permittivity of free space in F m^-1 + name = "e0" + unit = "F m^-1" + value = self.qce.e0 + docstring = f"| Permittivity of free space (unit: {unit}). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation." + self.add_constant(name, unit, value, docstring) + + # Avogadro constant in mol^-1 + name = "N_a" + unit = "mol^-1" + value = self.qce.na + docstring = f"| Avogadro constant (unit: {unit}). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation." + self.add_constant(name, unit, value, docstring) + + # meter to Bohr radius + name = "meter2bohr" + unit = "m^-1" + value = 1 / self.qce.bohr2m + docstring = f"| conversion factor from meter to Bohr radius (unit: {unit}). Affected code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation." + self.add_constant(name, unit, value, docstring) + # Set our constants to instance attributes for easy access for name, _, value, _ in self.data: self.__setattr__(name.lower(), float(value)) diff --git a/python/mrchem/validators.py b/python/mrchem/validators.py index a6b4f60c4..98c9c2a02 100644 --- a/python/mrchem/validators.py +++ b/python/mrchem/validators.py @@ -140,7 +140,7 @@ def __init__(self, user_dict, origin): self.atomic_coords = self.ang2bohr_array(self.atomic_coords) # Cavity related data - if user_dict["WaveFunction"]["environment"].lower() == "pcm": + if "pcm" in user_dict["WaveFunction"]["environment"].lower(): self.cavity_dict = user_dict["PCM"]["Cavity"] self.cavity_mode = self.cavity_dict["mode"] self.spheres_raw = self.cavity_dict["spheres"] diff --git a/python/template.yml b/python/template.yml index 4dca4880c..75e0fd55c 100644 --- a/python/template.yml +++ b/python/template.yml @@ -295,10 +295,11 @@ sections: type: str default: none predicates: - - value.lower() in ['none', 'pcm'] + - value.lower() in ['none', 'pcm', 'pcm_pb', 'pcm_lpb'] docstring: | Set method for treatment of environment. ``none`` for vacuum calculation. ``PCM`` for Polarizable Continuum Model, - which will activate the ``PCM`` input section for further parametrization options. + which will activate the ``PCM`` input section for further parametrization options. The ``PB`` and ``LPB`` variants add the + Poisson-Boltzmann and Linearized Poisson-Boltzmann solvers, respectively. - name: nuclear_model type: str default: point_like @@ -844,6 +845,83 @@ sections: ``total`` uses the total charge density. ``nuclear`` uses only the nuclear part of the total charge density. ``electronic`` uses only the electronic part of the total charge density. + - name: Solvent + docstring: | + Parameters for the Self-Consistent Reaction Field optimization. + sections: + - name: Permittivity + docstring: | + Parameters for the permittivity function. + keywords: + - name: epsilon_in + type: float + default: 1.0 + docstring: | + Permittivity inside the cavity. 1.0 is the permittivity of free + space, anything other than this is undefined behaviour. + - name: formulation + type: str + default: exponential + predicates: + - value.lower() in ['exponential'] + docstring: | + Formulation of the Permittivity function. Currently only the + 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']['Solvent']['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: DebyeHuckelScreening + docstring: | + Parameters for the Debye-Huckel screening factor + keywords: + - name: ion_strength + type: float + default: 1.0 + docstring: | + Ionic strength of the electrolyte in mol/L. This represents the concentration of the ions in the bulk solvent. + - name: ion_radius + type: float + default: 0.0 + docstring: | + Amount with which the vdw-radius of the atoms will be increased. + The screening factor will have an area of effect that is often going + to be larger than the vdw-cavity, but centered in the same atoms. + - name: ion_width + type: float + default: 0.2 + docstring: | + Width of the transition between the solute and the ion accessible part. + - name: formulation + type: str + default: variable + predicates: + - value.lower() in ['variable'] + docstring: | + formulation of the debye-huckel screening factor. Currently only the variable factor is implemented. + ``variable``: implement the screening functions as k = (1-C_ion)k_out - name: Cavity docstring: | Define the interlocking spheres cavity. @@ -926,56 +1004,12 @@ sections: Width of cavity boundary, smaller value means sharper transition. **This quantity has dimensions of length. The default value is in atomic units**. - - name: Permittivity - docstring: | - Parameters for the permittivity function. - keywords: - - name: epsilon_in - type: float - default: 1.0 - docstring: | - Permittivity inside the cavity. 1.0 is the permittivity of free - space, anything other than this is undefined behaviour. - - name: formulation - type: str - default: exponential - predicates: - - value.lower() in ['exponential'] - docstring: | - Formulation of the Permittivity function. Currently only the - 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: GeometryOptimizer docstring: | Includes parameters related to the internal geometry optimization using the SQNM (Stabilized Quasi-Newton Method) for noisy PES. Geometry optimizations require accurate forces. Consider setting world_prec to 1e-5 to 1e-7. - Convergence issues can usually be solved by increasing the precision of the SCF calculation. If that does + Convergence issues can usually be solved by increasing the precision of the SCF calculation. If that does not work, try setting the initial step size manually. keywords: - name: run @@ -1085,6 +1119,32 @@ sections: type: float docstring: '| Conversion factor for dipoles from atomic units to Debye (unit: ?). Affected code: Printed value of dipole moments.' + - name: boltzmann_constant + default: 1.380649e-23 + type: float + docstring: '| Boltzmann constant in (unit: J K^-1). Affected code: Value of + the Debye-Huckel screening parameter in the Poisson-Boltzmann equation.' + - name: elementary_charge + default: 1.602176634e-19 + type: float + docstring: '| Elementary charge in (unit: C). Affected code: Value of the + Debye-Huckel screening parameter in the Poisson-Boltzmann equation.' + - name: e0 + default: 8.8541878128e-12 + type: float + docstring: '| Permittivity of free space (unit: F m^-1). Affected code: Value + of the Debye-Huckel screening parameter in the Poisson-Boltzmann equation.' + - name: N_a + default: 6.02214076e+23 + type: float + docstring: '| Avogadro constant (unit: mol^-1). Affected code: Value of the + Debye-Huckel screening parameter in the Poisson-Boltzmann equation.' + - name: meter2bohr + default: 18897261246.2577 + type: float + docstring: '| conversion factor from meter to Bohr radius (unit: m^-1). Affected + code: Value of the Debye-Huckel screening parameter in the Poisson-Boltzmann + equation.' - name: Elements sections: - name: h diff --git a/recipes/Singularity.nompi b/recipes/Singularity.nompi index ab3d05c9b..1f5e7258e 100644 --- a/recipes/Singularity.nompi +++ b/recipes/Singularity.nompi @@ -45,7 +45,7 @@ Stage: build %post apt-get update -y DEBIAN_FRONTEND=noninteractive apt-get install -y --no-install-recommends \ - python3 + python3.9 rm -rf /var/lib/apt/lists/* @@ -78,7 +78,7 @@ From: ubuntu:20.04 %post apt-get update -y DEBIAN_FRONTEND=noninteractive apt-get install -y --no-install-recommends \ - python3 + python3.9 rm -rf /var/lib/apt/lists/* diff --git a/recipes/recipe_nompi.py b/recipes/recipe_nompi.py index 8b5c56ed1..cadfc6c69 100644 --- a/recipes/recipe_nompi.py +++ b/recipes/recipe_nompi.py @@ -5,7 +5,7 @@ Ubuntu 20.04 GNU compilers (upstream) CMake 3.20.6 - Python3 + Python3.9 MRChem (current source version) Generating recipe (stdout): diff --git a/recipes/recipe_openmpi4.0.py b/recipes/recipe_openmpi4.0.py index c3b6882df..b5fbb1863 100644 --- a/recipes/recipe_openmpi4.0.py +++ b/recipes/recipe_openmpi4.0.py @@ -9,7 +9,7 @@ UCX OpenMPI 4.0.5 CMake 3.20.6 - Python3 + Python3.9 MRChem (current source version) Generating recipe (stdout): @@ -111,7 +111,7 @@ (--dryrun) the container on the main input file, here called molecule.inp: $ ./.sif --dryrun molecule - + This will produce a new file molecule.json in the current directory which can be passed to the mrchem.x program inside the container using the singularity exec command: diff --git a/src/driver.cpp b/src/driver.cpp index 59c63c79a..3acffc882 100644 --- a/src/driver.cpp +++ b/src/driver.cpp @@ -77,8 +77,10 @@ #include "scf_solver/LinearResponseSolver.h" #include "environment/Cavity.h" +#include "environment/GPESolver.h" +#include "environment/LPBESolver.h" +#include "environment/PBESolver.h" #include "environment/Permittivity.h" -#include "environment/SCRF.h" #include "mrdft/Factory.h" @@ -1053,6 +1055,7 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild auto cavity_p = mol.getCavity_p(); cavity_p->printParameters(); + auto solver_type = json_fock["reaction_operator"]["solver_type"]; 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"]; @@ -1069,6 +1072,12 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild auto noneq = json_fock["reaction_operator"]["nonequilibrium"]; auto formulation = json_fock["reaction_operator"]["formulation"]; + // info for the poisson-boltzmann solver, this will anyways be ignored if we are not using it + double ion_radius = json_fock["reaction_operator"]["ion_radius"]; + auto kappa_o = json_fock["reaction_operator"]["kappa_out"]; + auto width_ion = json_fock["reaction_operator"]["ion_width"]; + auto kformulation = json_fock["reaction_operator"]["formulation"]; + // compute nuclear charge density Density rho_nuc(false); rho_nuc = chemistry::compute_nuclear_density(poisson_prec, nuclei, 100); @@ -1091,7 +1100,34 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild 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); + std::unique_ptr scrf_p; + + std::shared_ptr cavity_ion; + if (ion_radius != 0.0) { + auto radii_0 = cavity_p->getOriginalRadii(); + auto radii_ion = std::vector(radii_0.size()); + + for (int i = 0; i < radii_0.size(); i++) { radii_ion[i] = radii_0[i] + ion_radius; } + auto cavity_centers = cavity_p->getCoordinates(); + cavity_ion = std::make_shared(cavity_centers, radii_ion, width_ion); + } else { + cavity_ion = cavity_p; + } + if (solver_type == "Poisson-Boltzmann") { + DHScreening dhscreening(cavity_ion, kappa_o, kformulation); // this is now deciding the pb formulation, but it really shouldn't, the formulation here is for the DHScreening where we + // have 4 different parametrizations, not all implemented yet. + dhscreening.printParameters(); + scrf_p = std::make_unique(dielectric_func, dhscreening, rho_nuc, P_p, D_p, kain, max_iter, dynamic_thrs, density_type); + } else if (solver_type == "Linearized_Poisson-Boltzmann") { + DHScreening dhscreening(cavity_ion, kappa_o, kformulation); // this is now deciding the pb formulation, but it really shouldn't, the formulation here is for the DHScreening where we + // have 4 different parametrizations, not all implemented yet. + dhscreening.printParameters(); + scrf_p = std::make_unique(dielectric_func, dhscreening, rho_nuc, P_p, D_p, kain, max_iter, dynamic_thrs, density_type); + } else if (solver_type == "Generalized_Poisson") { + scrf_p = std::make_unique(dielectric_func, rho_nuc, P_p, D_p, kain, max_iter, dynamic_thrs, density_type); + } else { + MSG_ERROR("Solver type not implemented"); + } // initialize reaction potential object auto V_R = [&] { diff --git a/src/environment/CMakeLists.txt b/src/environment/CMakeLists.txt index a61dc42ec..5e7f108c5 100644 --- a/src/environment/CMakeLists.txt +++ b/src/environment/CMakeLists.txt @@ -1,6 +1,9 @@ target_sources(mrchem PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/Cavity.cpp ${CMAKE_CURRENT_SOURCE_DIR}/Permittivity.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/SCRF.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/DHScreening.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/GPESolver.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/PBESolver.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/LPBESolver.cpp ${CMAKE_CURRENT_SOURCE_DIR}/StepFunction.cpp ) diff --git a/src/environment/Cavity.cpp b/src/environment/Cavity.cpp index b7fda58a5..d2b9d128d 100644 --- a/src/environment/Cavity.cpp +++ b/src/environment/Cavity.cpp @@ -37,8 +37,7 @@ namespace mrchem { namespace detail { /** @relates mrchem::Cavity * @brief Constructs a single element of the gradient of the Cavity. - * - * This constructs the analytical partial derivative of the Cavity \f$C\f$ with respect to \f$x\f$, \f$y\f$ or \f$z\f$ + * @details This constructs the analytical partial derivative of the Cavity \f$C\f$ with respect to \f$x\f$, \f$y\f$ or \f$z\f$ * coordinates and evaluates it at a point \f$\mathbf{r}\f$. This is given for \f$x\f$ by * \f[ * \frac{\partial C\left(\mathbf{r}\right)}{\partial x} = \left(1 - C{\left(\mathbf{r} \right)}\right) diff --git a/src/environment/Cavity.h b/src/environment/Cavity.h index 25e299a35..495d6dd4a 100644 --- a/src/environment/Cavity.h +++ b/src/environment/Cavity.h @@ -34,7 +34,7 @@ namespace mrchem { /** @class Cavity * @brief Interlocking spheres cavity centered on the nuclei of the molecule. - * The Cavity class represents the following function \cite Fosso-Tande2013 + * @details The Cavity class represents the following function \cite Fosso-Tande2013 * * \f[ * C(\mathbf{r}) = 1 - \prod^N_{i=1} (1-C_i(\mathbf{r})) \\ @@ -59,12 +59,13 @@ namespace mrchem { * - \f$R_{0,i}\f$ is the atomic radius. By default, the van der Waals radius. * - \f$\alpha_{i}\f$ is a scaling factor. By default, 1.1 * - \f$\beta_{i}\f$ is a width scaling factor. By default, 0.5 - * - \f$\sigma_{i}\f$ is the width. By default, 0.2 + * - \f$\sigma_{i}\f$ is the width. By default, 0.2 bohr */ 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); + /** @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: @@ -77,9 +78,18 @@ class Cavity final : public mrcpp::RepresentableFunction<3> { */ 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; diff --git a/src/environment/DHScreening.cpp b/src/environment/DHScreening.cpp new file mode 100644 index 000000000..5d6b53dd9 --- /dev/null +++ b/src/environment/DHScreening.cpp @@ -0,0 +1,40 @@ +/* + * 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 "DHScreening.h" + +namespace mrchem { + +DHScreening::DHScreening(std::shared_ptr cavity_ion, double kappa_out, const std::string &formulation) + : StepFunction(cavity_ion, 0.0, kappa_out) + , formulation(formulation) {} + +double DHScreening::evalf(const mrcpp::Coord<3> &r) const { + auto c_pin = this->cavity_p; + auto kappa_squared = std::pow(this->out, 2) * (1 - c_pin->evalf(r)); + return kappa_squared; +} + +} // namespace mrchem diff --git a/src/environment/DHScreening.h b/src/environment/DHScreening.h new file mode 100644 index 000000000..5f5a5a575 --- /dev/null +++ b/src/environment/DHScreening.h @@ -0,0 +1,81 @@ +/* + * 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 + +#include +#include + +#include "StepFunction.h" +#include "utils/print_utils.h" + +namespace mrchem { +class Cavity; +/** @class DHScreening + * + * @brief Square of the Debye-Huckel Screening parameter. + * @details + * This is used for the Poisson-Boltzmann solver. The DHScreening function is defined as + * \f[ + * \kappa^2(\mathbf{r}) = \begin{cases} + * \kappa^2_{out} & \text{if } \mathbf{r} \notin \Omega_{ion} \\ + * 0.0 & \text{if } \mathbf{r} \in \Omega_{ion} + * \end{cases} + * \f] + * This can be parametrized a number of ways. The one used here is + * \f[ + * \kappa^2(\mathbf{r}) = (1 - C_{ion}(\mathbf{r})) \kappa^2_{out} + * \f] + * Where \f$C_{ion}(\mathbf{r})\f$ is the ion accessible Cavity function. + */ +class DHScreening final : public StepFunction { +public: + /** @brief Standard constructor. Initializes the #cavity_ion and #kappa_out with the input parameters. + * @param cavity_ion interlocking spheres of Cavity class. + * @param kappa_out value of the screening function outside the #cavity_ion. + * @param formulation Decides which formulation of the #DHScreening function to implement, only continuous screening function available. + * @details #kappa_out is given by + * \f[ + * \kappa = \sqrt{\frac{2000 I_{0} e^2 N_a I_0}{\epsilon_{out} \epsilon_{in} k_B T}} + * \f] + * where \f$N_a\f$ is the Avogadro constant, e is the elementary charge, \f$I_0\f$ is the concentration of the ions, + * \f$k_B\f$ is the Boltzmann constant, \f$T\f$ is the temperature, \f$\epsilon_{out}\f$ is the permittivity of the solvent and \f$\epsilon_{in}\f$ is the permittivity of free space. + */ + DHScreening(std::shared_ptr cavity_ion, double kappa_out, const std::string &formulation); + + /** @brief Evaluates DHScreening at a point in 3D space. + * @param r coordinates of a 3D point in space. + * @return Value at point r. + */ + double evalf(const mrcpp::Coord<3> &r) const override; + void printParameters() const override { detail::print_header("Ion-accessible Cavity", this->formulation, getValueIn(), getValueOut()); } + +private: + std::string formulation{"Continuous Screening Function"}; //!< Formulation of the DHScreening function. Only linear variable is used now. +}; + +} // namespace mrchem diff --git a/src/environment/SCRF.cpp b/src/environment/GPESolver.cpp similarity index 82% rename from src/environment/SCRF.cpp rename to src/environment/GPESolver.cpp index 308cc90ee..df7141c04 100644 --- a/src/environment/SCRF.cpp +++ b/src/environment/GPESolver.cpp @@ -23,7 +23,7 @@ * */ -#include "SCRF.h" +#include "GPESolver.h" #include #include @@ -48,7 +48,7 @@ using DerivativeOperator_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, SCRFDensityType density_type) +GPESolver::GPESolver(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) @@ -59,23 +59,23 @@ SCRF::SCRF(const Permittivity &e, const Density &rho_nuc, PoissonOperator_p P, D , derivative(D) , poisson(P) {} -SCRF::~SCRF() { +GPESolver::~GPESolver() { this->rho_nuc.free(NUMBER::Real); clear(); } -void SCRF::clear() { +void GPESolver::clear() { this->apply_prec = -1.0; } -double SCRF::setConvergenceThreshold(double prec) { +double GPESolver::setConvergenceThreshold(double prec) { // converge_thrs should be in the interval [prec, 1.0] this->conv_thrs = prec; if (this->dynamic_thrs and this->mo_residual > 10 * prec) this->conv_thrs = std::min(1.0, this->mo_residual); return this->conv_thrs; } -void SCRF::computeDensities(const Density &rho_el, Density &rho_out) { +void GPESolver::computeDensities(const Density &rho_el, Density &rho_out) { Timer timer; switch (this->density_type) { @@ -96,13 +96,13 @@ void SCRF::computeDensities(const Density &rho_el, Density &rho_out) { print_utils::qmfunction(3, "Vacuum density", rho_out, timer); } -void SCRF::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma) { +void GPESolver::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma) { auto d_V = mrcpp::gradient(*derivative, potential.real()); // FunctionTreeVector resetComplexFunction(out_gamma); for (int d = 0; d < 3; d++) { - auto C_pin = this->epsilon.getCavity(); + auto C_pin = this->epsilon.getCavity_p(); mrcpp::AnalyticFunction<3> d_cav(C_pin->getGradVector()[d]); mrcpp::ComplexFunction cplxfunc_prod; mrcpp::cplxfunc::multiply(cplxfunc_prod, get_func(d_V, d), d_cav, this->apply_prec, 1); @@ -118,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, const Density &rho_el) { +mrcpp::ComplexFunction GPESolver::solvePoissonEquation(const mrcpp::ComplexFunction &in_gamma, const Density &rho_el) { mrcpp::ComplexFunction Poisson_func; mrcpp::ComplexFunction rho_eff; mrcpp::ComplexFunction first_term; @@ -139,7 +139,7 @@ mrcpp::ComplexFunction SCRF::solvePoissonEquation(const mrcpp::ComplexFunction & return Vr_np1; } -void SCRF::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain) { +void GPESolver::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain) { OrbitalVector phi_n(0); OrbitalVector dPhi_n(0); phi_n.push_back(Orbital(SPIN::Paired)); @@ -160,14 +160,14 @@ void SCRF::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFu dPhi_n.clear(); } -void SCRF::nestedSCRF(const mrcpp::ComplexFunction &V_vac, const Density &rho_el) { +void GPESolver::runMicroIterations(const mrcpp::ComplexFunction &V_vac, const Density &rho_el) { KAIN kain(this->history); kain.setLocalPrintLevel(10); mrcpp::print::separator(3, '-'); + auto update = 10.0, norm = 1.0; - double update = 10.0, norm = 1.0; - int iter = 1; + auto iter = 1; while (update >= this->conv_thrs && iter <= max_iter) { Timer t_iter; // solve the poisson equation @@ -197,15 +197,19 @@ void SCRF::nestedSCRF(const mrcpp::ComplexFunction &V_vac, const Density &rho_el Vr_np1.free(NUMBER::Real); printConvergenceRow(iter, norm, update, t_iter.elapsed()); + + // compute new PB term + if (update < this->conv_thrs) break; iter++; } + if (iter > max_iter) println(0, "Reaction potential failed to converge after " << iter - 1 << " iterations, residual " << update); mrcpp::print::separator(3, '-'); kain.clear(); } -void SCRF::printConvergenceRow(int i, double norm, double update, double time) const { +void GPESolver::printConvergenceRow(int i, double norm, double update, double time) const { auto pprec = Printer::getPrecision(); auto w0 = Printer::getWidth() - 1; auto w1 = 9; @@ -229,7 +233,7 @@ void SCRF::printConvergenceRow(int i, double norm, double update, double time) c println(3, o_txt.str()); } -mrcpp::ComplexFunction &SCRF::setup(double prec, const Density &rho_el) { +mrcpp::ComplexFunction &GPESolver::solveEquation(double prec, const Density &rho_el) { this->apply_prec = prec; Density rho_tot(false); computeDensities(rho_el, rho_tot); @@ -253,51 +257,34 @@ mrcpp::ComplexFunction &SCRF::setup(double prec, const Density &rho_el) { // update the potential/gamma before doing anything with them Timer t_scrf; - nestedSCRF(V_vac, rho_el); + runMicroIterations(V_vac, rho_el); print_utils::qmfunction(3, "Reaction potential", this->Vr_n, t_scrf); return this->Vr_n; } -auto SCRF::computeEnergies(const Density &rho_el) -> std::tuple { +auto GPESolver::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 {Er_el, Er_nuc}; } -void SCRF::resetComplexFunction(mrcpp::ComplexFunction &function) { +void GPESolver::resetComplexFunction(mrcpp::ComplexFunction &function) { if (function.hasReal()) function.free(NUMBER::Real); if (function.hasImag()) function.free(NUMBER::Imag); function.alloc(NUMBER::Real); } -void SCRF::printParameters() const { - std::stringstream o_iter; - if (this->max_iter > 0) { - o_iter << this->max_iter; - } else { - o_iter << "Off"; - } - - std::stringstream o_kain; - if (this->history > 0) { - o_kain << this->history; - } else { - o_kain << "Off"; - } - +void GPESolver::printParameters() const { nlohmann::json data = { - {"Method ", "SCRF"}, - {"Optimizer ", "Potential"}, - {"Max iterations ", max_iter}, - {"KAIN solver ", o_kain.str()}, - {"Density type ", density_type}, - {"Dynamic threshold ", (dynamic_thrs) ? "On" : "Off"}, + {"Method ", this->solver_name}, + {"Density ", this->density_type}, + {"Max iterations ", this->max_iter}, + {"KAIN solver ", (this->history > 0) ? std::to_string(this->history) : "Off"}, + {"Dynamic threshold ", (this->dynamic_thrs) ? "On" : "Off"}, }; mrcpp::print::separator(3, '~'); print_utils::json(3, data, false); mrcpp::print::separator(3, '~'); } - } // namespace mrchem diff --git a/src/environment/GPESolver.h b/src/environment/GPESolver.h new file mode 100644 index 000000000..acd8adf8e --- /dev/null +++ b/src/environment/GPESolver.h @@ -0,0 +1,223 @@ +/* + * 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 +#include + +#include +#include + +#include "DHScreening.h" +#include "Permittivity.h" +#include "qmfunctions/Density.h" + +namespace mrchem { +class KAIN; +class ReactionPotential; +class ReactionPotentialD1; +class ReactionPotentialD2; + +enum class SCRFDensityType : int { TOTAL = 0, ELECTRONIC = 1, NUCLEAR = 2 }; + +/** @class GPESolver + * @brief Solves the Generalized Poisson equation iteratively + * @details The Generalized Poisson equation is given by + * \f[ + * \nabla \cdot \left( \epsilon(\mathbf{r}) \nabla V(\mathbf{r}) \right) = -4\pi \rho(\mathbf{r}) + * \f] + * where \f$\epsilon(\mathbf{r})\f$ is the permittivity, \f$V(\mathbf{r})\f$ is the total electrostatic potential and \f$\rho(\mathbf{r})\f$ is the molecular charge density defined as: + * \f[ + * \rho(\mathbf{r}) = \rho_{el}(\mathbf{r}) + \rho_{nuc}(\mathbf{r}) + * \f] + * where \f$\rho_{el}\f$ is the electronic charge density and \f$\rho_{nuc}\f$ is the nuclear charge density. + * The Generalized Poisson equation is solved iteratively through a set of micro-iteration on each SCF-iteration by appliation of the Poisson operator \f$\mathcal{P}\f$ :cite:`Fosso-Tande2013` + * \f[ + * V_R(\mathbf{r}) = \mathcal{P} \star \left[ \rho_{eff}(\mathbf{r}) - \rho(\mathbf{r}) + \gamma_s(\mathbf{r}) \right] + * \f] + * where \f$\gamma_s(\mathbf{r})\f$ is the surface charge distribution describing the polarization at the surface, \f$\rho_{eff}(\mathbf{r})\f$ is the effective charge density given by + * \f$\frac{\rho(\mathbf{r})}{\epsilon(\mathbf{r})}\f$ and \f$V_R(\mathbf{r})\f$ is the reaction potential. + * + * We utilize a so-called dynamic threshold to more easily converge the reaction potential. This is done by setting the convergence threshold of the micro-iterations to + * the MO update of the previous SCF iteration, unless the MO update is small enough (once the quality of the MOs is good enough, we use the default convergence threshold). + * Another optimization used is that we utilize the previous SCF converged Reaction potential as an initial guess for the next micro-iterations. These procedures are + * investigated and explained in :cite:`gerez2023` + */ +class GPESolver { +public: + GPESolver(const Permittivity &e, + const Density &rho_nuc, + std::shared_ptr P, + std::shared_ptr> D, + int kain_hist, + int max_iter, + bool dyn_thrs, + SCRFDensityType density_type); + ~GPESolver(); + + /** @brief Sets the convergence threshold for the micro-iterations, used with dynamic thresholding. + * @param prec value to set the convergence threshold to + * @return the current convergence threshold. + * @details will check if the MO update is small enough (ten times as big) wrt. to the scf convergence threshold, if so, it will use the default convergence threshold. + * If not, it will use the MO update as the convergence threshold. + */ + double setConvergenceThreshold(double prec); + + Permittivity &getPermittivity() { return this->epsilon; } + + void updateMOResidual(double const err_t) { this->mo_residual = err_t; } + + /** @brief Computes the energy contributions from the reaction potential + * @param rho_el the electronic charge density + * @return a tuple containing the electronic and nuclear energy contributions + * @details We compute the reaction energy through the following integral: + * \f[ + * E_{R} = \frac{1}{2}\int \rho_{el}(\mathbf{r}) V_R(\mathbf{r}) d\mathbf{r} + \frac{1}{2} \int \rho_{nuc}(\mathbf{r}) V_R(\mathbf{r}) d\mathbf{r} + * \f] + * Each term represents the electronic and nuclear contributions to the reaction energy, respectively. + * We compute each term separately, and return a tuple containing both. + */ + 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: + bool dynamic_thrs; + SCRFDensityType density_type; //!< Decides which density we will use for computing the reaction potential, options are ``total``, ``electronic`` and ``nuclear``. + + int max_iter; + int history; + double apply_prec{-1.0}; + double conv_thrs{1.0}; + double mo_residual{1.0}; + std::string solver_name{"Generalized Poisson"}; + + Permittivity epsilon; + + Density rho_nuc; // As of right now, this is the biggest memory hog. + // Alternative could be to precompute its contributions, as a potential is not as heavy as a density (maybe) + // another one could be to define a representable function which only has the exact analytical form of the nuclear contribution. + // Since we already have \nabla^2 V_nuc = -4*pi*rho_nuc (nuclear potential) we could use this as a way to bypass computing rho_nuc at all + // Same with the coulomb potential, which is basically what is left from V_vac after subtracting V_nuc. in one way we could just precompute both and + // just iterate through V_R only. Only issue here is (1 -1\epsilon)/\epsilon * \rho_nuc as I am not sure how to represent this as an analytitcal function, + // maybe after applying the Poisson operator? + + mrcpp::ComplexFunction Vr_n; + + std::shared_ptr> derivative; + std::shared_ptr poisson; + + void clear(); + + /** @brief computes density wrt. the density_type variable + * @param Phi the molecular orbitals + * @param rho_out Density function in which the density will be computed. + * @details The total charge density is given by the sum of the electronic and nuclear charge densities: + * \f[ + * \rho(\mathbf{r}) = \rho_{el}(\mathbf{r}) + \rho_{nuc}(\mathbf{r}) + * \f] + * where \f$\rho_{el}\f$ is the electronic charge density and \f$\rho_{nuc}\f$ is the nuclear charge density. + * The nuclear charge density is stored in the class variable #rho_nuc, while we compute the electronic charge density + * from the molecular orbitals. + * The class variable #density_type decides the density which will be computed in rho_out, options are ``total``, ``electronic`` and ``nuclear``. + */ + void computeDensities(const Density &rho_el, Density &rho_out); + + /** @brief Computes the surface charge distibution due to polarization at the solute-solvent boundary + * @param potential Potential used to compute \f$\nabla V(\mathbf{r})\f$ + * @param out_gamma ComplexFunction in which the surface charge distribution will be computed. + * @details The surface charge distribution is given by + * \f[ + * \gamma_s(\mathbf{r}) = \frac{\log \frac{\epsilon_{in}}{\epsilon_{out}}}{4 \pi } \left( \nabla C(\mathbf{r}) \cdot \nabla V(\mathbf{r})\right) + * \f] + * where \f$\epsilon_{in}\f$ is the permittivity inside the cavity and \f$\epsilon_{out}\f$ is the permittivity outside the cavity. + */ + virtual void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma); + + /** @brief Iterates once through the Generalized Poisson equation to compute the reaction potential + * @param ingamma the surface charge distribution + * @param Phi the molecular orbitals + * @return the reaction potential + * @details Constructs the effective charge density \f$\rho_{eff}(\mathbf{r})\f$ and the Poisson operator \f$\mathcal{P}\f$ as: + * \f[ + * V_R(\mathbf{r}) = \mathcal{P} \star \left[ \rho_{eff}(\mathbf{r}) - \rho(\mathbf{r}) + \gamma_s(\mathbf{r}) \right] + * \f] + * where \f$\gamma_s(\mathbf{r})\f$ is the surface charge distribution describing the polarization at the surface, \f$\rho_{eff}(\mathbf{r})\f$ is the effective charge density given by + * \f$\frac{\rho(\mathbf{r})}{\epsilon(\mathbf{r})}\f$ and \f$V_R(\mathbf{r})\f$ is the reaction potential. + */ + mrcpp::ComplexFunction solvePoissonEquation(const mrcpp::ComplexFunction &ingamma, const Density &rho_el); + + /** @brief Uses KAIN to accelerate convergece of the reaction potential + * @param dfunc the current update of the reaction potential + * @param func the current reaction potential + * @param kain the KAIN object + */ + void accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain); + + /** @brief Iterates through the application of the Poisson operator to Solve the Generalized Poisson equation + * @param V_vac the vacuum potential + * @param Phi_p the molecular orbitals + * @details Iterating through the application of the Poisson operator is done through a set of micro-iterations, + * where the convergence threshold is set to the MO update of the previous SCF iteration. + * The micro-iterations are done through the following steps: + * -# Compute the total potential as \f$V(\mathbf{r}) = V_{vac}(\mathbf{r}) + V_R(\mathbf{r})\f$ + * -# Compute the surface charge distribution \f$\gamma_s(\mathbf{r})\f$ with #computeGamma + * -# Compute a new reaction potential \f$V_R(\mathbf{r})\f$ by application of the Poisson operator with #solvePoissonEquation + * -# Calculate the update of the reaction potential as \f$\Delta V_R(\mathbf{r}) = V_R(\mathbf{r}) - V_R^{old}(\mathbf{r})\f$ + * -# Accelerate convergence of the reaction potential through KAIN + * -# Update the reaction potential as \f$V_R(\mathbf{r}) = V_R^{old}(\mathbf{r}) + \Delta V_R(\mathbf{r})\f$ + * -# Check if the reaction potential has converged, if not, repeat from step 1. + */ + void runMicroIterations(const mrcpp::ComplexFunction &V_vac, const Density &rho_el); + + /** @brief Setups and computes the reaction potential through the microiterations + * @param V_vac the vacuum potential + * @param Phi_p the molecular orbitals + * @return The converged reaction potential for the current SCF iteration + * @details An initial guess of the reaction potential is computed with the following steps: + * -# Set the total potential as \f$V(\mathbf{r}) = V_{vac}(\mathbf{r})\f$ + * -# Compute the surface charge distribution \f$\gamma_s(\mathbf{r})\f$ from this potential + * -# Iterate once through the application of the Poisson operator to return the initial guess of the reaction potential \f$V_R(\mathbf{r})\f$ + * + * the method then runs the micro-iterations through #runMicroIterations and returns the converged reaction potential. + * If this is not the first SCF iteration, the previous converged reaction potential is used as an initial guess for the micro-iterations. + */ + mrcpp::ComplexFunction &solveEquation(double prec, const Density &rho_el); + + /** @brief Frees the memory used by the FunctionTrees of the input Complexfunction and reallocates them. + * @param function the ComplexFunction to reset + * @details This is done to avoid memory leaks when the ComplexFunction is used in the micro-iterations. + */ + void resetComplexFunction(mrcpp::ComplexFunction &function); + + virtual void printParameters() const; + void printConvergenceRow(int i, double norm, double update, double time) const; +}; +} // namespace mrchem diff --git a/src/environment/LPBESolver.cpp b/src/environment/LPBESolver.cpp new file mode 100644 index 000000000..6d6bb6868 --- /dev/null +++ b/src/environment/LPBESolver.cpp @@ -0,0 +1,65 @@ +/* + * 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 "LPBESolver.h" + +#include +#include +#include +#include + +#include "GPESolver.h" +#include "chemistry/PhysicalConstants.h" +#include "chemistry/chemistry_utils.h" +#include "qmfunctions/density_utils.h" +#include "qmoperators/two_electron/ReactionPotential.h" +#include "utils/print_utils.h" + +using mrcpp::Printer; +using mrcpp::Timer; + +using PoissonOperator_p = std::shared_ptr; +using DerivativeOperator_p = std::shared_ptr>; +using OrbitalVector_p = std::shared_ptr; + +namespace mrchem { +LPBESolver::LPBESolver(const Permittivity &e, + const DHScreening &k, + const Density &rho_nuc, + std::shared_ptr P, + std::shared_ptr> D, + int kain_hist, + int max_iter, + bool dyn_thrs, + SCRFDensityType density_type) + : PBESolver(e, k, rho_nuc, P, D, kain_hist, max_iter, dyn_thrs, density_type) {} +// TODO separate this for the linear and non-linear solver +void LPBESolver::computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) { + resetComplexFunction(pb_term); + mrcpp::cplxfunc::multiply(pb_term, V_tot, this->kappa, this->apply_prec); + pb_term.rescale(salt_factor / (4.0 * mrcpp::pi)); +} + +} // namespace mrchem diff --git a/src/environment/LPBESolver.h b/src/environment/LPBESolver.h new file mode 100644 index 000000000..25f617d3c --- /dev/null +++ b/src/environment/LPBESolver.h @@ -0,0 +1,70 @@ +/* + * 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 "DHScreening.h" +#include "PBESolver.h" +#include "Permittivity.h" +#include "qmfunctions/Density.h" +#include "qmfunctions/Orbital.h" + +namespace mrchem { +/** @class LPBESolver + * @brief Solves the Linearized Poisson-Boltzmann equation iteratively + * @details The Linearized Poisson-Boltzmann equation is solved iteratively using the SCRF procedure outlined in GPESolver and PBESolver. + * The linearized Poisson-Boltzmann equation is a further simplification of the Poisson-Boltzmann equation, outlined in PBESolver, where the PB term is expanded and + * only the linear term is included. This is a good approximation for low ionic strength solutions. + * The linearized Poisson-Boltzmann equation is given by + * \f[ + * \nabla^2 V_{R} = -4\pi\frac{1-\epsilon}{\epsilon}\left(\rho_{el} + \rho_{nuc}\right) + \gamma_s - \kappa^2 V_{tot} + * \f] + * where \f$\gamma_s\f$ is the surface charge density, \f$\kappa\f$ is obtained from the DHScreening class and \f$V_{R}\f$ is the reaction potential. + */ +class LPBESolver final : public PBESolver { +public: + LPBESolver(const Permittivity &e, + const DHScreening &k, + const Density &rho_nuc, + std::shared_ptr P, + std::shared_ptr> D, + int kain_hist, + int max_iter, + bool dyn_thrs, + SCRFDensityType density_type); + + friend class ReactionPotential; + +protected: + /** @brief Computes the PB term + * @param[in] V_tot the total potential + * @param[in] salt_factor the salt factor deciding how much of the total concentration to include in the PB term + * @param[out] pb_term the ComplexFunction in which to store the result + * @details The PB term is computed as \f$ \kappa^2 V_{tot} \f$ and returned. + */ + void computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) override; + std::string solver_name{"Linearized Poisson-Boltzmann"}; +}; +} // namespace mrchem diff --git a/src/environment/PBESolver.cpp b/src/environment/PBESolver.cpp new file mode 100644 index 000000000..6c5e22fd6 --- /dev/null +++ b/src/environment/PBESolver.cpp @@ -0,0 +1,99 @@ +/* + * 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 "PBESolver.h" + +#include +#include +#include +#include + +#include "GPESolver.h" +#include "chemistry/PhysicalConstants.h" +#include "chemistry/chemistry_utils.h" +#include "qmfunctions/density_utils.h" +#include "utils/print_utils.h" + +using mrcpp::Printer; +using mrcpp::Timer; + +using PoissonOperator_p = std::shared_ptr; +using DerivativeOperator_p = std::shared_ptr>; +using OrbitalVector_p = std::shared_ptr; + +namespace mrchem { + +PBESolver::PBESolver(const Permittivity &e, + const DHScreening &k, + const Density &rho_nuc, + std::shared_ptr P, + std::shared_ptr> D, + int kain_hist, + int max_iter, + bool dyn_thrs, + SCRFDensityType density_type) + : GPESolver(e, rho_nuc, P, D, kain_hist, max_iter, dyn_thrs, density_type) + , kappa(k) {} + +void PBESolver::computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) { + // create a lambda function for the sinh(V) term and multiply it with kappa and salt factor to get the PB term + auto sinh_f = [salt_factor](const double &V) { return (salt_factor / (4.0 * mrcpp::pi)) * std::sinh(V); }; + resetComplexFunction(pb_term); + mrcpp::ComplexFunction sinhV; + sinhV.alloc(NUMBER::Real); + mrcpp::map(this->apply_prec / 100, sinhV.real(), V_tot.real(), sinh_f); + + mrcpp::cplxfunc::multiply(pb_term, sinhV, this->kappa, this->apply_prec); +} + +void PBESolver::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma) { + + auto d_V = mrcpp::gradient(*derivative, potential.real()); // FunctionTreeVector + resetComplexFunction(out_gamma); + + for (int d = 0; d < 3; d++) { + auto C_pin = this->epsilon.getCavity_p(); + mrcpp::AnalyticFunction<3> d_cav(C_pin->getGradVector()[d]); + mrcpp::ComplexFunction cplxfunc_prod; + mrcpp::cplxfunc::multiply(cplxfunc_prod, get_func(d_V, d), d_cav, this->apply_prec, 1); + // add result into out_gamma + if (d == 0) { + mrcpp::cplxfunc::deep_copy(out_gamma, cplxfunc_prod); + } else { + out_gamma.add(1.0, cplxfunc_prod); + } + } + + out_gamma.rescale(std::log((epsilon.getValueIn() / epsilon.getValueOut())) * (1.0 / (4.0 * mrcpp::pi))); + mrcpp::clear(d_V, true); + + // add PB term + mrcpp::ComplexFunction pb_term; + auto salt_factor = 1.0; // placeholder for now, want to change it wrt to convergence in future + computePBTerm(potential, salt_factor, pb_term); + out_gamma.add(-1.0, pb_term); +} + +} // namespace mrchem diff --git a/src/environment/PBESolver.h b/src/environment/PBESolver.h new file mode 100644 index 000000000..c5215d55b --- /dev/null +++ b/src/environment/PBESolver.h @@ -0,0 +1,94 @@ +/* + * 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 "DHScreening.h" +#include "GPESolver.h" +#include "Permittivity.h" +#include "qmfunctions/Density.h" +#include "qmfunctions/Orbital.h" + +namespace mrchem { +/** @class PBESolver + * @brief Solves the Poisson-Boltzmann equation iteratively + * @details The Poisson-Boltzmann equation is solved iteratively using the SCRF procedure outlined in GPESolver. + * The Poisson-Boltzmann equation models the electrostatic potential in a solvent with electrolytes. + * The general equation for electrolyte solutions is given by + * \f[ + * \nabla \cdot \epsilon \nabla V_{tot} = -4\pi\left(\rho_{el} + \rho_{nuc} + \rho_{ext}\right) + * \f] + * where \f$ V_{tot} \f$ is the total electrostatic potential, \f$ \epsilon\f$ is the permittivity function of the solvent, + * \f$\rho_{el}\f$ is the electronic charge density, \f$\rho_{nuc}\f$ is the nuclear charge density and \f$\rho_{ext}\f$ is the external charge density. + * In the general form for the Poisson-Boltzmann equation, the external charge density is approximated by assuming a boltzmann distribution of the ions. + * \f[ + * \rho_{ext} = \sum_{i}^{N_{ion}} q_i e I_{0, i}\exp\left(-\frac{q_i e V_{tot}}{k_B T}\right) + * \f] + * where \f$I_{0, i}\f$ is the concentration of the i-th ion species, \f$q_i\f$ is the charge of the i-th ion species, \f$k_B\f$ is the Boltzmann constant and \f$T\f$ is the temperature. + * In this implementation we assume a 1:1 (\f$I_{0, 0} = I_{0, 1}\f$) electrolyte soluttion of ions of same opposite charges (\f$z_i = +1, -1\f$). This simplifies the external density to + * \f[ + * \rho_{ext} = -2 e I_{0} \sinh\left(\frac{e V_{tot}}{2 k_B T}\right) + * \f] + * where \f$I_{0}\f$ is the concentration of the ions. + * We can plug this into the first equation (and massage terms a bit) to arrive at the Poisson-Boltzmann equation for 1:1 electrolyte solution + * \f[ + * \nabla^2 V_{R} = -4\pi\frac{1-\epsilon}{\epsilon}\left(\rho_{el} + \rho_{nuc}\right) + \gamma_s - \kappa^2 \sinh\left(V_{tot}\right) + * \f] + * where \f$\gamma_s\f$ is the surface charge density, \f$\kappa\f$ is obtained from the DHScreening class and \f$V_{R}\f$ is the reaction potential. + */ +class PBESolver : public GPESolver { +public: + PBESolver(const Permittivity &e, + const DHScreening &k, + const Density &rho_nuc, + std::shared_ptr P, + std::shared_ptr> D, + int kain_hist, + int max_iter, + bool dyn_thrs, + SCRFDensityType density_type); + + friend class ReactionPotential; + +protected: + DHScreening kappa; ///< the DHScreening object used to compute the PB term \f$\kappa\f$ + std::string solver_name{"Poisson-Boltzmann"}; + + /** @brief constructs the surface chage distribution and adds it to the PB term + * @param[in] potential the potential to compute \f$\nabla V\f$ from + * @param[out] out_gamma the ComplexFunction in which to store the result + * @details Method follows the implementation in GPESolver::computeGamma, but adds the PB term to the surface charge distribution. + */ + void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma) override; + + /** @brief Computes the PB term + * @param[in] V_tot the total potential + * @param[in] salt_factor the salt factor deciding how much of the total concentration to include in the PB term + * @param[out] pb_term the ComplexFunction in which to store the result + * @details The PB term is computed as \f$ \kappa^2 \sinh(V_{tot}) \f$ and returned. + */ + virtual void computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term); +}; +} // namespace mrchem diff --git a/src/environment/Permittivity.cpp b/src/environment/Permittivity.cpp index f62d0a38b..852c85aff 100644 --- a/src/environment/Permittivity.cpp +++ b/src/environment/Permittivity.cpp @@ -35,7 +35,7 @@ Permittivity::Permittivity(std::shared_ptr cavity, double epsilo , formulation(formulation) {} double Permittivity::evalf(const mrcpp::Coord<3> &r) const { - auto c_pin = this->cavity; + auto c_pin = this->cavity_p; return this->in * std::exp(std::log(this->out / this->in) * (1 - c_pin->evalf(r))); } diff --git a/src/environment/Permittivity.h b/src/environment/Permittivity.h index 09e73a337..bbac87199 100644 --- a/src/environment/Permittivity.h +++ b/src/environment/Permittivity.h @@ -25,6 +25,8 @@ #pragma once +#include + #include #include diff --git a/src/environment/SCRF.h b/src/environment/SCRF.h deleted file mode 100644 index 573edbae7..000000000 --- a/src/environment/SCRF.h +++ /dev/null @@ -1,113 +0,0 @@ -/* - * 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 -#include - -#include -#include - -#include "Permittivity.h" -#include "qmfunctions/Density.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. - */ -class SCRF final { -public: - SCRF(const Permittivity &e, - const Density &rho_nuc, - std::shared_ptr P, - std::shared_ptr> D, - int kain_hist, - int max_iter, - bool dyn_thrs, - SCRFDensityType density_type); - ~SCRF(); - - double setConvergenceThreshold(double prec); - - Permittivity &getPermittivity() { return this->epsilon; } - - void updateMOResidual(double const err_t) { this->mo_residual = err_t; } - - 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; - SCRFDensityType density_type; - - int max_iter; - int history; - double apply_prec{-1.0}; - double conv_thrs{1.0}; - double mo_residual{1.0}; - - Permittivity epsilon; - - Density rho_nuc; // As of right now, this is the biggest memory hog. - // Alternative could be to precompute its contributions, as a potential is not as heavy as a density (maybe) - // another one could be to define a representable function which only has the exact analytical form of the nuclear contribution. - - mrcpp::ComplexFunction Vr_n; - - std::shared_ptr> derivative; - std::shared_ptr poisson; - - 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, const Density &rho_el); - - void accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain); - - 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); - - void printParameters() const; - void printConvergenceRow(int i, double norm, double update, double time) const; -}; -} // namespace mrchem diff --git a/src/environment/StepFunction.cpp b/src/environment/StepFunction.cpp index 7d5a35e7f..f967f7907 100644 --- a/src/environment/StepFunction.cpp +++ b/src/environment/StepFunction.cpp @@ -43,5 +43,5 @@ void print_header(const std::string &header, const std::string &formulation, dou StepFunction::StepFunction(std::shared_ptr cavity, double val_in, double val_out) : in(val_in) , out(val_out) - , cavity{std::move(cavity)} {} + , cavity_p{std::move(cavity)} {} } // namespace mrchem diff --git a/src/environment/StepFunction.h b/src/environment/StepFunction.h index fb4ced35a..406cd4792 100644 --- a/src/environment/StepFunction.h +++ b/src/environment/StepFunction.h @@ -66,14 +66,14 @@ class StepFunction : public mrcpp::RepresentableFunction<3> { auto getValueIn() const { return this->in; } auto getValueOut() const { return this->out; } - std::shared_ptr getCavity() const { return this->cavity; } + std::shared_ptr getCavity_p() const { return this->cavity_p; } virtual void printParameters() const { detail::print_header("Step function", "Standard", getValueIn(), getValueOut()); } protected: - double in; //!< Value of the function inside the #cavity. - double out; //!< value of the function outside the #cavity. - std::shared_ptr cavity; //!< A Cavity class instance. + double in; //!< Value of the function inside the #cavity. + double out; //!< value of the function outside the #cavity. + std::shared_ptr cavity_p; //!< A shared pointer to a Cavity class instance. }; } // namespace mrchem diff --git a/src/qmoperators/two_electron/FockBuilder.cpp b/src/qmoperators/two_electron/FockBuilder.cpp index 60909b62f..ad9aaa5d9 100644 --- a/src/qmoperators/two_electron/FockBuilder.cpp +++ b/src/qmoperators/two_electron/FockBuilder.cpp @@ -173,7 +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::tie(Er_el, Er_nuc) = this->Ro->getHelper()->computeEnergies(rho_el); + std::tie(Er_el, Er_nuc) = this->Ro->getSolver()->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 a3d07a0e5..3f859c671 100644 --- a/src/qmoperators/two_electron/ReactionOperator.h +++ b/src/qmoperators/two_electron/ReactionOperator.h @@ -29,30 +29,31 @@ #include "ReactionPotentialD1.h" #include "ReactionPotentialD2.h" -#include "environment/SCRF.h" +#include "environment/GPESolver.h" /** @class ReactionOperator * * @brief Operator containing a single ReactionPotential * * This class is a simple TensorOperator realization of @class ReactionPotential. + * */ namespace mrchem { class ReactionOperator final : public RankZeroOperator { public: - 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); + ReactionOperator(std::unique_ptr gpesolver_p, std::shared_ptr Phi_p, bool mpi_share = false) { + potential = std::make_shared(std::move(gpesolver_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); + ReactionOperator(std::unique_ptr gpesolver_p, std::shared_ptr Phi, std::shared_ptr X, std::shared_ptr Y, bool mpi_share = false) { + // check that the GPESolver object uses the electronic density only + if (gpesolver_p->getDensityType() != SCRFDensityType::ELECTRONIC) MSG_ERROR("Invalid SCRF object passed: only electronic density in response"); + potential = std::make_shared(std::move(gpesolver_p), Phi, X, Y, mpi_share); // Invoke operator= to assign *this operator RankZeroOperator &V = (*this); V = potential; @@ -61,7 +62,7 @@ class ReactionOperator final : public RankZeroOperator { ComplexDouble trace(OrbitalVector &Phi) { return RankZeroOperator::trace(Phi); } - SCRF *getHelper() { return this->potential->getHelper(); } + GPESolver *getSolver() { return this->potential->getSolver(); } std::shared_ptr getPotential() { return this->potential; } void updateMOResidual(double const err_t) { this->potential->updateMOResidual(err_t); } diff --git a/src/qmoperators/two_electron/ReactionPotential.cpp b/src/qmoperators/two_electron/ReactionPotential.cpp index 7a28bb44d..1a6b671c0 100644 --- a/src/qmoperators/two_electron/ReactionPotential.cpp +++ b/src/qmoperators/two_electron/ReactionPotential.cpp @@ -29,34 +29,29 @@ #include "ReactionPotential.h" #include "qmfunctions/density_utils.h" -using mrcpp::Printer; -using mrcpp::Timer; - -using SCRF_p = std::unique_ptr; +using GPESolver_p = std::unique_ptr; using OrbitalVector_p = std::shared_ptr; namespace mrchem { -ReactionPotential::ReactionPotential(SCRF_p scrf, OrbitalVector_p Phi, bool mpi_share) +ReactionPotential::ReactionPotential(GPESolver_p gpesolver_p, OrbitalVector_p Phi, bool mpi_share) : QMPotential(1, mpi_share) - , helper(std::move(scrf)) + , solver(std::move(gpesolver_p)) , orbitals(Phi) {} void ReactionPotential::setup(double prec) { if (isSetup(prec)) return; setApplyPrec(prec); - auto thrs = this->helper->setConvergenceThreshold(prec); - Timer timer; - auto plevel = Printer::getPrintLevel(); + auto thrs = this->solver->setConvergenceThreshold(prec); + mrcpp::Timer timer; + auto plevel = mrcpp::Printer::getPrintLevel(); mrcpp::print::separator(3, '='); print_utils::centered_text(3, "Building Reaction operator"); - this->helper->printParameters(); + this->solver->printParameters(); mrcpp::print::value(3, "Precision", prec, "(rel)", 5); mrcpp::print::value(3, "Threshold", thrs, "(abs)", 5); mrcpp::print::separator(3, '-'); - 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); @@ -65,7 +60,7 @@ void ReactionPotential::setup(double prec) { void ReactionPotential::clear() { mrcpp::ComplexFunction::free(NUMBER::Total); // delete FunctionTree pointers clearApplyPrec(); - this->helper->clear(); + this->solver->clear(); } } // namespace mrchem diff --git a/src/qmoperators/two_electron/ReactionPotential.h b/src/qmoperators/two_electron/ReactionPotential.h index ffbbb678a..6e37fc6ea 100644 --- a/src/qmoperators/two_electron/ReactionPotential.h +++ b/src/qmoperators/two_electron/ReactionPotential.h @@ -25,10 +25,9 @@ #pragma once +#include "environment/GPESolver.h" #include "qmoperators/QMPotential.h" -#include "environment/SCRF.h" - namespace mrchem { /** @class ReactionPotential * @brief class containing the solvent-substrate interaction reaction potential @@ -42,24 +41,22 @@ namespace mrchem { class ReactionPotential : public QMPotential { public: /** @brief Initializes the ReactionPotential class. - * @param scrf A SCRF instance which contains the parameters needed to compute the ReactionPotential. + * @param scrf A GPESolver 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); + explicit ReactionPotential(std::unique_ptr scrf, std::shared_ptr Phi = nullptr, bool mpi_share = false); ~ReactionPotential() override = default; - SCRF *getHelper() { return this->helper.get(); } + GPESolver *getSolver() { return this->solver.get(); } - /** @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; } + /** @brief Updates the solver.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->solver->mo_residual = err_t; } friend class ReactionOperator; protected: - std::unique_ptr helper; ///< A SCRF instance used to compute the ReactionPotential. + std::unique_ptr solver; //!< A GPESolver instance used to compute the ReactionPotential. std::shared_ptr orbitals; ///< Unperturbed orbitals defining the ground-state electron density for the SCRF procedure. void setup(double prec) override; diff --git a/src/qmoperators/two_electron/ReactionPotentialD1.cpp b/src/qmoperators/two_electron/ReactionPotentialD1.cpp index 349209fb4..2b12e30fd 100644 --- a/src/qmoperators/two_electron/ReactionPotentialD1.cpp +++ b/src/qmoperators/two_electron/ReactionPotentialD1.cpp @@ -23,11 +23,11 @@ * */ -#include "MRCPP/Printer" -#include "MRCPP/Timer" - #include "ReactionPotentialD1.h" -#include "environment/SCRF.h" + +#include +#include + #include "qmfunctions/density_utils.h" #include "utils/print_utils.h" @@ -46,6 +46,6 @@ mrcpp::ComplexFunction &ReactionPotentialD1::computePotential(double prec) const // change sign, because it's the electronic density rho_el.rescale(-1.0); - return this->helper->setup(prec, rho_el); + return this->solver->solveEquation(prec, rho_el); } -} // namespace mrchem \ No newline at end of file +} // namespace mrchem diff --git a/src/qmoperators/two_electron/ReactionPotentialD1.h b/src/qmoperators/two_electron/ReactionPotentialD1.h index f51da4e7e..e470fb170 100644 --- a/src/qmoperators/two_electron/ReactionPotentialD1.h +++ b/src/qmoperators/two_electron/ReactionPotentialD1.h @@ -26,14 +26,14 @@ #pragma once #include "ReactionPotential.h" -#include "environment/SCRF.h" +#include "environment/GPESolver.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) {} + ReactionPotentialD1(std::unique_ptr gpesolver, std::shared_ptr Phi, bool mpi_share = false) + : ReactionPotential(std::move(gpesolver), Phi, mpi_share) {} private: mrcpp::ComplexFunction &computePotential(double prec) const override; diff --git a/src/qmoperators/two_electron/ReactionPotentialD2.cpp b/src/qmoperators/two_electron/ReactionPotentialD2.cpp index c8dea2f51..03328b0ed 100644 --- a/src/qmoperators/two_electron/ReactionPotentialD2.cpp +++ b/src/qmoperators/two_electron/ReactionPotentialD2.cpp @@ -23,10 +23,11 @@ * */ -#include "MRCPP/Printer" -#include "MRCPP/Timer" - #include "ReactionPotentialD2.h" + +#include +#include + #include "qmfunctions/density_utils.h" #include "utils/print_utils.h" @@ -51,6 +52,6 @@ mrcpp::ComplexFunction &ReactionPotentialD2::computePotential(double prec) const // change sign, because it's the electronic density rho.rescale(-1.0); - return this->helper->setup(prec, rho); + return this->solver->solveEquation(prec, rho); } } // namespace mrchem diff --git a/src/qmoperators/two_electron/ReactionPotentialD2.h b/src/qmoperators/two_electron/ReactionPotentialD2.h index 8cc846d2f..c77aafcbc 100644 --- a/src/qmoperators/two_electron/ReactionPotentialD2.h +++ b/src/qmoperators/two_electron/ReactionPotentialD2.h @@ -26,13 +26,13 @@ #pragma once #include "ReactionPotential.h" -#include "environment/SCRF.h" +#include "environment/GPESolver.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) + ReactionPotentialD2(std::unique_ptr gpesolver, std::shared_ptr Phi, std::shared_ptr X, std::shared_ptr Y, bool mpi_share = false) + : ReactionPotential(std::move(gpesolver), Phi, mpi_share) , orbitals_x(X) , orbitals_y(Y) {} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 10d08d523..e7a1417fd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -30,6 +30,8 @@ add_subdirectory(hf_grad_lda) add_subdirectory(cube_parser) add_subdirectory(h2_scf_cube) add_subdirectory(li_solv) +add_subdirectory(h_pb) +add_subdirectory(h_lpb) 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 e451f93a4..8489eedb0 100644 --- a/tests/cavity_input_parser/world_unit=angstrom-default.inp +++ b/tests/cavity_input_parser/world_unit=angstrom-default.inp @@ -18,9 +18,11 @@ WaveFunction { } PCM { - Permittivity { - outside = { - epsilon_static = 80.0 + Solvent { + Permittivity { + 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 ab2a50307..e38cb7f31 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 @@ -18,9 +18,11 @@ WaveFunction { } PCM { - Permittivity { - outside = { - epsilon_static = 80.0 + Solvent { + Permittivity { + outside = { + epsilon_static = 80.0 + } } } Cavity { 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 9a8b14791..95c0f0a8f 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 @@ -18,10 +18,12 @@ WaveFunction { } PCM { - Permittivity { - outside = { - epsilon_static = 80.0 - } + Solvent { + Permittivity { + 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 e8dae2340..9e6686b33 100644 --- a/tests/cavity_input_parser/world_unit=angstrom-mode=explicit.inp +++ b/tests/cavity_input_parser/world_unit=angstrom-mode=explicit.inp @@ -18,9 +18,11 @@ WaveFunction { } PCM { - Permittivity { - outside = { - epsilon_static = 80.0 + Solvent { + Permittivity { + outside = { + epsilon_static = 80.0 + } } } Cavity { diff --git a/tests/cavity_input_parser/world_unit=bohr-default.inp b/tests/cavity_input_parser/world_unit=bohr-default.inp index 65fe0526d..7d35425cc 100644 --- a/tests/cavity_input_parser/world_unit=bohr-default.inp +++ b/tests/cavity_input_parser/world_unit=bohr-default.inp @@ -18,9 +18,11 @@ WaveFunction { } PCM { - Permittivity { - outside = { - epsilon_static = 80.0 + Solvent { + Permittivity { + 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 a84cc0139..b917f907c 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 @@ -18,9 +18,11 @@ WaveFunction { } PCM { - Permittivity { - outside = { - epsilon_static = 80.0 + Solvent { + Permittivity { + outside = { + epsilon_static = 80.0 + } } } Cavity { 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 776ff230d..446e7d9e8 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 @@ -18,9 +18,11 @@ WaveFunction { } PCM { - Permittivity { - outside = { - epsilon_static = 80.0 + Solvent { + Permittivity { + outside = { + epsilon_static = 80.0 + } } } Cavity { 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 f1618596d..219b574c7 100644 --- a/tests/cavity_input_parser/world_unit=bohr-mode=explicit.inp +++ b/tests/cavity_input_parser/world_unit=bohr-mode=explicit.inp @@ -18,9 +18,11 @@ WaveFunction { } PCM { - Permittivity { - outside = { - epsilon_static = 80.0 + Solvent { + Permittivity { + outside = { + epsilon_static = 80.0 + } } } Cavity { diff --git a/tests/h2_pol_solv/h2-eq.inp b/tests/h2_pol_solv/h2-eq.inp index 88c2556d8..802ae4571 100644 --- a/tests/h2_pol_solv/h2-eq.inp +++ b/tests/h2_pol_solv/h2-eq.inp @@ -18,11 +18,13 @@ Properties { } PCM { - Permittivity { - epsilon_out { - static = 4.0 - dynamic = 2.0 - nonequilibrium = false + Solvent { + Permittivity { + epsilon_out { + static = 4.0 + dynamic = 2.0 + nonequilibrium = false + } } } } diff --git a/tests/h2_pol_solv/h2-neq.inp b/tests/h2_pol_solv/h2-neq.inp index 7af73eced..d0f8adc46 100644 --- a/tests/h2_pol_solv/h2-neq.inp +++ b/tests/h2_pol_solv/h2-neq.inp @@ -18,11 +18,13 @@ Properties { } PCM { - Permittivity { - epsilon_out { - static = 4.0 - dynamic = 2.0 - nonequilibrium = true + Solvent { + Permittivity { + epsilon_out { + static = 4.0 + dynamic = 2.0 + nonequilibrium = true + } } } } diff --git a/tests/h_lpb/CMakeLists.txt b/tests/h_lpb/CMakeLists.txt new file mode 100644 index 000000000..ff600cf17 --- /dev/null +++ b/tests/h_lpb/CMakeLists.txt @@ -0,0 +1,10 @@ +if(ENABLE_MPI) + set(_h_lpb_launcher "${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 1") +endif() + +add_integration_test( + NAME "H_linear_poisson_boltzmann" + LABELS "mrchem;h_lpb;linear_poisson_boltzmann;scf;energy" + COST 100 + LAUNCH_AGENT ${_h_lpb_launcher} + ) diff --git a/tests/h_lpb/h.inp b/tests/h_lpb/h.inp new file mode 100644 index 000000000..35985d29f --- /dev/null +++ b/tests/h_lpb/h.inp @@ -0,0 +1,41 @@ +{ +"world_prec": 1.0e-4, +"world_size": 5, +"MPI": { + "numerically_exact": true +}, +"Molecule": { + "charge": -1, + "coords": "H 0.0 0.0 0.0" +}, +"WaveFunction": { + "method": "pbe0", + "environment": "pcm_lpb" +}, +"PCM": { + "SCRF": { + "kain": 6, + "max_iter": 100, + "dynamic_thrs": false + }, + "Cavity": { + "spheres": "0 2.645616384 1.0 0.0 0.2" + }, + "Solvent": { + "Permittivity": { + "epsilon_in": 1.0, + "epsilon_out": { "static": 78.4}, + "formulation": "exponential" + }, + "DebyeHuckelScreening": { + "ion_strength": 0.25, + "ion_radius": 0.0, + "ion_width": 0.2 + } + } +}, +"SCF": { + "run": false, + "guess_type": "sad_gto" +} +} diff --git a/tests/h_lpb/reference/h.json b/tests/h_lpb/reference/h.json new file mode 100644 index 000000000..cd9e7a311 --- /dev/null +++ b/tests/h_lpb/reference/h.json @@ -0,0 +1,273 @@ +{ + "input": { + "constants": { + "N_a": 6.02214076e+23, + "angstrom2bohrs": 1.8897261246257702, + "boltzmann_constant": 1.380649e-23, + "dipmom_au2debye": 2.5417464739297717, + "e0": 8.8541878128e-12, + "electron_g_factor": -2.00231930436256, + "elementary_charge": 1.602176634e-19, + "fine_structure_constant": 0.0072973525693, + "hartree2ev": 27.211386245988, + "hartree2kcalmol": 627.5094740630558, + "hartree2kjmol": 2625.4996394798254, + "hartree2simagnetizability": 78.9451185, + "hartree2wavenumbers": 219474.6313632, + "light_speed": 137.035999084, + "meter2bohr": 18897261246.2577 + }, + "molecule": { + "cavity": { + "spheres": [ + { + "alpha": 1.0, + "beta": 0.0, + "center": [ + 0.0, + 0.0, + 0.0 + ], + "radius": 2.645616384, + "sigma": 0.2 + } + ] + }, + "charge": -1, + "coords": [ + { + "atom": "h", + "r_rms": 2.6569547399e-05, + "xyz": [ + 0.0, + 0.0, + 0.0 + ] + } + ], + "multiplicity": 1 + }, + "mpi": { + "bank_size": -1, + "numerically_exact": true, + "shared_memory_size": 10000 + }, + "mra": { + "basis_order": 6, + "basis_type": "interpolating", + "boxes": [ + 2, + 2, + 2 + ], + "corner": [ + -1, + -1, + -1 + ], + "max_scale": 20, + "min_scale": -4 + }, + "printer": { + "file_name": "h", + "print_constants": false, + "print_level": 0, + "print_mpi": false, + "print_prec": 6, + "print_width": 75 + }, + "rsp_calculations": {}, + "scf_calculation": { + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.0001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.0001 + }, + "kinetic_operator": { + "derivative": "abgv_55" + }, + "nuclear_operator": { + "nuclear_model": "point_like", + "proj_prec": 0.0001, + "shared_memory": false, + "smooth_prec": 0.0001 + }, + "reaction_operator": { + "DHS-formulation": "variable", + "density_type": "total", + "dynamic_thrs": false, + "epsilon_in": 1.0, + "epsilon_out": 78.4, + "formulation": "exponential", + "ion_radius": 0.0, + "ion_width": 0.2, + "kain": 6, + "kappa_out": 0.08703231499578493, + "max_iter": 100, + "poisson_prec": 0.0001, + "solver_type": "Linearized_Poisson-Boltzmann" + }, + "xc_operator": { + "shared_memory": false, + "xc_functional": { + "cutoff": 0.0, + "functionals": [ + { + "coef": 1.0, + "name": "pbe0" + } + ], + "spin": false + } + } + }, + "initial_guess": { + "environment": "None", + "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": "DFT (PBE0)", + "prec": 0.001, + "relativity": "None", + "restricted": true, + "screen": 12.0, + "type": "sad_gto", + "zeta": 0 + }, + "properties": { + "dipole_moment": { + "dip-1": { + "operator": "h_e_dip", + "precision": 0.0001, + "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": -1, + "dipole_moment": { + "dip-1": { + "magnitude": 7.715572311852292e-13, + "r_O": [ + 0.0, + 0.0, + 0.0 + ], + "vector": [ + 0.0, + 0.0, + 0.0 + ], + "vector_el": [ + 0.0, + 0.0, + 0.0 + ], + "vector_nuc": [ + 0.0, + 0.0, + 0.0 + ] + } + }, + "geometry": [ + { + "symbol": "H", + "xyz": [ + 0.0, + 0.0, + 0.0 + ] + } + ], + "multiplicity": 1, + "orbital_energies": { + "energy": [ + -0.11654513076836014 + ], + "occupation": [ + 2.0 + ], + "spin": [ + "p" + ], + "sum_occupied": -0.23309026153672027 + }, + "scf_energy": { + "E_ee": 1.220280348193976, + "E_eext": 0.0, + "E_el": -0.7847284165934681, + "E_en": -1.9145538948639653, + "E_kin": 0.9258603759485317, + "E_next": 0.0, + "E_nn": 0.0, + "E_nuc": 0.19089592768192742, + "E_tot": -0.5938324889115407, + "E_x": -0.15252084470528926, + "E_xc": -0.4885331476131019, + "Er_el": -0.3752612535536192, + "Er_nuc": 0.19089592768192742, + "Er_tot": -0.1843653258716918 + } + }, + "provenance": { + "creator": "MRChem", + "mpi_processes": 1, + "nthreads": 1, + "routine": "mrchem.x", + "total_cores": 1, + "version": "1.2.0-alpha" + }, + "rsp_calculations": null, + "scf_calculation": { + "initial_energy": { + "E_ee": 1.220280348193976, + "E_eext": 0.0, + "E_el": -0.7847284165934681, + "E_en": -1.9145538948639653, + "E_kin": 0.9258603759485317, + "E_next": 0.0, + "E_nn": 0.0, + "E_nuc": 0.19089592768192742, + "E_tot": -0.5938324889115407, + "E_x": -0.15252084470528926, + "E_xc": -0.4885331476131019, + "Er_el": -0.3752612535536192, + "Er_nuc": 0.19089592768192742, + "Er_tot": -0.1843653258716918 + }, + "success": true + }, + "schema_name": "mrchem_output", + "schema_version": 1, + "success": true + } +} diff --git a/tests/h_lpb/reference/h.out b/tests/h_lpb/reference/h.out new file mode 100644 index 000000000..517379217 --- /dev/null +++ b/tests/h_lpb/reference/h.out @@ -0,0 +1,290 @@ + + +*************************************************************************** +*** *** +*** *** +*** __ __ ____ ____ _ *** +*** | \/ | _ \ / ___| |__ ___ _ __ ___ *** +*** | |\/| | |_) | | | '_ \ / _ \ '_ ` _ \ *** +*** | | | | _ <| |___| | | | __/ | | | | | *** +*** |_| |_|_| \_\\____|_| |_|\___|_| |_| |_| *** +*** *** +*** VERSION 1.2.0-alpha *** +*** *** +*** Git branch unknown *** +*** Git commit hash unknown *** +*** Git commit author unknown *** +*** Git commit date unknown *** +*** *** +*** Contact: luca.frediani@uit.no *** +*** *** +*** Radovan Bast Magnar Bjorgve *** +*** Roberto Di Remigio Antoine Durdek *** +*** Luca Frediani Gabriel Gerez *** +*** Stig Rune Jensen Jonas Juselius *** +*** Rune Monstad Peter Wind *** +*** *** +*************************************************************************** + +--------------------------------------------------------------------------- + + MPI processes : (no bank) 1 + OpenMP threads : 1 + Total cores : 1 + +--------------------------------------------------------------------------- + +XCFun DFT library Copyright 2009-2020 Ulf Ekstrom and contributors. +See http://dftlibs.org/xcfun/ for more information. + +This is free software; see the source code for copying conditions. +There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. For details see the documentation. +Scientific users of this library should cite +U. Ekstrom, L. Visscher, R. Bast, A. J. Thorvaldsen and K. Ruud; +J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s + +--------------------------------------------------------------------------- + + MRCPP version : 1.6.0-alpha + Git branch : unknown + Git commit hash : unknown + Git commit author : unknown + Git commit date : unknown + + Linear algebra : EIGEN v3.4.0 + Parallelization : OpenMP (1 threads) + +--------------------------------------------------------------------------- + + + +=========================================================================== + MultiResolution Analysis +--------------------------------------------------------------------------- + polynomial order : 6 + polynomial type : Interpolating +--------------------------------------------------------------------------- + total boxes : 8 + boxes : [ 2 2 2 ] + unit lengths : [ 16.00000 16.00000 16.00000 ] + scaling factor : [ 1.00000 1.00000 1.00000 ] + lower bounds : [ -16.00000 -16.00000 -16.00000 ] + upper bounds : [ 16.00000 16.00000 16.00000 ] + total length : [ 32.00000 32.00000 32.00000 ] +=========================================================================== + + + +*************************************************************************** +*** *** +*** Initializing Molecule *** +*** *** +*************************************************************************** + + +=========================================================================== + Molecule +--------------------------------------------------------------------------- + Charge : -1 + Multiplicity : 1 +--------------------------------------------------------------------------- + N Atom : x y z +--------------------------------------------------------------------------- + 0 H : 0.000000 0.000000 0.000000 +--------------------------------------------------------------------------- + Center of mass : 0.000000 0.000000 0.000000 +=========================================================================== + + +=========================================================================== + Solvation Cavity +--------------------------------------------------------------------------- + Formulation : exponential + Dielectric constant : (in) 1.000000 + : (out) 78.400000 +--------------------------------------------------------------------------- + N R_0 Alpha Beta Sigma Radius x y z +--------------------------------------------------------------------------- + 0 2.6456 1.00 0.00 0.20 -> 2.6456 0.000000 0.000000 0.000000 +=========================================================================== + + +=========================================================================== + Square of the Debye-Huckel screening parameter +--------------------------------------------------------------------------- + Formulation : exponential + Screening function value: (in) 0.000000 + : (out) 0.087032 +--------------------------------------------------------------------------- + N R_0 Sigma Radius x y z +--------------------------------------------------------------------------- + 0 2.6456 0.20 -> 2.6456 0.000000 0.000000 0.000000 +=========================================================================== + + + +*************************************************************************** +*** *** +*** Computing Initial Guess Wavefunction *** +*** *** +*************************************************************************** + + +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Calculation : Compute initial orbitals + Method : Diagonalize SAD Hamiltonian + Precision : 1.00000e-03 + Screening : 1.20000e+01 StdDev + Restricted : True + Functional : LDA (SVWN5) + AO basis : 3-21G +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +=========================================================================== + Molecular Orbitals +--------------------------------------------------------------------------- + Alpha electrons : 1 + Beta electrons : 1 + Total electrons : 2 +--------------------------------------------------------------------------- + n Occ Spin : Norm +--------------------------------------------------------------------------- + 0 2 p : 9.999999944384e-01 +=========================================================================== + + +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Calculation : Compute initial energy + Method : DFT (PBE0) + Relativity : None + Environment : None + External fields : None + Precision : 1.00000e-03 + Localization : Off +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +=========================================================================== + Molecular Energy (initial) +--------------------------------------------------------------------------- + Kinetic energy : (au) 0.925860375949 + E-N energy : (au) -1.914553894864 + Coulomb energy : (au) 1.220280348194 + Exchange energy : (au) -0.152520844705 + X-C energy : (au) -0.488533147613 + N-N energy : (au) 0.000000000000 +--------------------------------------------------------------------------- + Reaction energy (el) : (au) -0.375261253554 + Reaction energy (nuc) : (au) 0.190895927682 + Reaction energy (tot) : (au) -0.184365325872 +--------------------------------------------------------------------------- + Electronic energy : (au) -0.784728416593 + Nuclear energy : (au) 0.190895927682 +--------------------------------------------------------------------------- + Total energy : (au) -5.938324889115e-01 + : (kcal/mol) -3.726355127984e+02 + : (kJ/mol) -1.559106985549e+03 + : (eV) -1.615900522119e+01 +=========================================================================== + + +=========================================================================== + Orbital Energies (initial) +--------------------------------------------------------------------------- + n Occ Spin : Epsilon +--------------------------------------------------------------------------- + 0 2 p : (au) -0.116545130768 +--------------------------------------------------------------------------- + Sum occupied : (au) -0.233090261537 +=========================================================================== + + + +*************************************************************************** +*** *** +*** Printing Molecular Properties *** +*** *** +*************************************************************************** + + +=========================================================================== + Molecule +--------------------------------------------------------------------------- + Charge : -1 + Multiplicity : 1 +--------------------------------------------------------------------------- + N Atom : x y z +--------------------------------------------------------------------------- + 0 H : 0.000000 0.000000 0.000000 +--------------------------------------------------------------------------- + Center of mass : 0.000000 0.000000 0.000000 +=========================================================================== + + +=========================================================================== + Molecular Energy (final) +--------------------------------------------------------------------------- + Kinetic energy : (au) 0.925860375949 + E-N energy : (au) -1.914553894864 + Coulomb energy : (au) 1.220280348194 + Exchange energy : (au) -0.152520844705 + X-C energy : (au) -0.488533147613 + N-N energy : (au) 0.000000000000 +--------------------------------------------------------------------------- + Reaction energy (el) : (au) -0.375261253554 + Reaction energy (nuc) : (au) 0.190895927682 + Reaction energy (tot) : (au) -0.184365325872 +--------------------------------------------------------------------------- + Electronic energy : (au) -0.784728416593 + Nuclear energy : (au) 0.190895927682 +--------------------------------------------------------------------------- + Total energy : (au) -5.938324889115e-01 + : (kcal/mol) -3.726355127984e+02 + : (kJ/mol) -1.559106985549e+03 + : (eV) -1.615900522119e+01 +=========================================================================== + + +=========================================================================== + Orbital Energies (final) +--------------------------------------------------------------------------- + n Occ Spin : Epsilon +--------------------------------------------------------------------------- + 0 2 p : (au) -0.116545130768 +--------------------------------------------------------------------------- + Sum occupied : (au) -0.233090261537 +=========================================================================== + + +=========================================================================== + Dipole Moment (dip-1) +--------------------------------------------------------------------------- + r_O : 0.000000 0.000000 0.000000 +--------------------------------------------------------------------------- + Electronic vector : 0.000000 0.000000 0.000000 + Magnitude : (au) 0.000000 + : (Debye) 0.000000 +--------------------------------------------------------------------------- + Nuclear vector : -0.000000 -0.000000 -0.000000 + Magnitude : (au) 0.000000 + : (Debye) 0.000000 +--------------------------------------------------------------------------- + Total vector : -0.000000 -0.000000 -0.000000 + Magnitude : (au) 0.000000 + : (Debye) 0.000000 +=========================================================================== + + + + +*************************************************************************** +*** *** +*** Exiting MRChem *** +*** *** +*** Wall time : 0h 0m 33s *** +*** *** +*************************************************************************** + + diff --git a/tests/h_lpb/test b/tests/h_lpb/test new file mode 100755 index 000000000..8fe8c91be --- /dev/null +++ b/tests/h_lpb/test @@ -0,0 +1,29 @@ +#!/usr/bin/env python3 + +import sys +from pathlib import Path + +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), + ER_TOT: rel_tolerance(1.0e-6), + ER_EL: rel_tolerance(1.0e-6), + ER_NUC: rel_tolerance(1.0e-6), +} + +ierr = run(options, input_file="h", filters=filters, extra_args=['--json']) + +sys.exit(ierr) diff --git a/tests/h_pb/CMakeLists.txt b/tests/h_pb/CMakeLists.txt new file mode 100644 index 000000000..a4f360299 --- /dev/null +++ b/tests/h_pb/CMakeLists.txt @@ -0,0 +1,10 @@ +if(ENABLE_MPI) + set(_h_pb_launcher "${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 1") +endif() + +add_integration_test( + NAME "H_poisson_boltzmann" + LABELS "mrchem;h_pb;poisson_boltzmann;scf;energy" + COST 100 + LAUNCH_AGENT ${_h_pb_launcher} + ) diff --git a/tests/h_pb/h.inp b/tests/h_pb/h.inp new file mode 100644 index 000000000..1fd09d87e --- /dev/null +++ b/tests/h_pb/h.inp @@ -0,0 +1,41 @@ +{ +"world_prec": 1.0e-4, +"world_size": 5, +"MPI": { + "numerically_exact": true +}, +"Molecule": { + "charge": -1, + "coords": "H 0.0 0.0 0.0" +}, +"WaveFunction": { + "method": "pbe0", + "environment": "pcm_pb" +}, +"PCM": { + "SCRF": { + "kain": 6, + "max_iter": 100, + "dynamic_thrs": false + }, + "Cavity": { + "spheres": "0 2.645616384 1.0 0.0 0.2" + }, + "Solvent":{ + "Permittivity": { + "epsilon_in": 1.0, + "epsilon_out": { "static": 78.4}, + "formulation": "exponential" + }, + "DebyeHuckelScreening": { + "ion_strength": 0.25, + "ion_radius": 0.0, + "ion_width": 0.2 + } + } +}, +"SCF": { + "run": false, + "guess_type": "sad_gto" +} +} diff --git a/tests/h_pb/reference/h.json b/tests/h_pb/reference/h.json new file mode 100644 index 000000000..0ed1227a7 --- /dev/null +++ b/tests/h_pb/reference/h.json @@ -0,0 +1,273 @@ +{ + "input": { + "constants": { + "N_a": 6.02214076e+23, + "angstrom2bohrs": 1.8897261246257702, + "boltzmann_constant": 1.380649e-23, + "dipmom_au2debye": 2.5417464739297717, + "e0": 8.8541878128e-12, + "electron_g_factor": -2.00231930436256, + "elementary_charge": 1.602176634e-19, + "fine_structure_constant": 0.0072973525693, + "hartree2ev": 27.211386245988, + "hartree2kcalmol": 627.5094740630558, + "hartree2kjmol": 2625.4996394798254, + "hartree2simagnetizability": 78.9451185, + "hartree2wavenumbers": 219474.6313632, + "light_speed": 137.035999084, + "meter2bohr": 18897261246.2577 + }, + "molecule": { + "cavity": { + "spheres": [ + { + "alpha": 1.0, + "beta": 0.0, + "center": [ + 0.0, + 0.0, + 0.0 + ], + "radius": 2.645616384, + "sigma": 0.2 + } + ] + }, + "charge": -1, + "coords": [ + { + "atom": "h", + "r_rms": 2.6569547399e-05, + "xyz": [ + 0.0, + 0.0, + 0.0 + ] + } + ], + "multiplicity": 1 + }, + "mpi": { + "bank_size": -1, + "numerically_exact": true, + "shared_memory_size": 10000 + }, + "mra": { + "basis_order": 6, + "basis_type": "interpolating", + "boxes": [ + 2, + 2, + 2 + ], + "corner": [ + -1, + -1, + -1 + ], + "max_scale": 20, + "min_scale": -4 + }, + "printer": { + "file_name": "h", + "print_constants": false, + "print_level": 0, + "print_mpi": false, + "print_prec": 6, + "print_width": 75 + }, + "rsp_calculations": {}, + "scf_calculation": { + "fock_operator": { + "coulomb_operator": { + "poisson_prec": 0.0001, + "shared_memory": false + }, + "exchange_operator": { + "exchange_prec": -1.0, + "poisson_prec": 0.0001 + }, + "kinetic_operator": { + "derivative": "abgv_55" + }, + "nuclear_operator": { + "nuclear_model": "point_like", + "proj_prec": 0.0001, + "shared_memory": false, + "smooth_prec": 0.0001 + }, + "reaction_operator": { + "DHS-formulation": "variable", + "density_type": "total", + "dynamic_thrs": false, + "epsilon_in": 1.0, + "epsilon_out": 78.4, + "formulation": "exponential", + "ion_radius": 0.0, + "ion_width": 0.2, + "kain": 6, + "kappa_out": 0.08703231499578493, + "max_iter": 100, + "poisson_prec": 0.0001, + "solver_type": "Poisson-Boltzmann" + }, + "xc_operator": { + "shared_memory": false, + "xc_functional": { + "cutoff": 0.0, + "functionals": [ + { + "coef": 1.0, + "name": "pbe0" + } + ], + "spin": false + } + } + }, + "initial_guess": { + "environment": "None", + "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": "DFT (PBE0)", + "prec": 0.001, + "relativity": "None", + "restricted": true, + "screen": 12.0, + "type": "sad_gto", + "zeta": 0 + }, + "properties": { + "dipole_moment": { + "dip-1": { + "operator": "h_e_dip", + "precision": 0.0001, + "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": -1, + "dipole_moment": { + "dip-1": { + "magnitude": 7.715572311852292e-13, + "r_O": [ + 0.0, + 0.0, + 0.0 + ], + "vector": [ + 0.0, + 0.0, + 0.0 + ], + "vector_el": [ + 0.0, + 0.0, + 0.0 + ], + "vector_nuc": [ + 0.0, + 0.0, + 0.0 + ] + } + }, + "geometry": [ + { + "symbol": "H", + "xyz": [ + 0.0, + 0.0, + 0.0 + ] + } + ], + "multiplicity": 1, + "orbital_energies": { + "energy": [ + -0.1165723464459209 + ], + "occupation": [ + 2.0 + ], + "spin": [ + "p" + ], + "sum_occupied": -0.2331446928918418 + }, + "scf_energy": { + "E_ee": 1.220280348193976, + "E_eext": 0.0, + "E_el": -0.7847556322707149, + "E_en": -1.9145538948639653, + "E_kin": 0.9258603759485317, + "E_next": 0.0, + "E_nn": 0.0, + "E_nuc": 0.19090981420562303, + "E_tot": -0.5938458180650918, + "E_x": -0.15252084470528926, + "E_xc": -0.4885331476131019, + "Er_el": -0.37528846923086595, + "Er_nuc": 0.19090981420562303, + "Er_tot": -0.18437865502524292 + } + }, + "provenance": { + "creator": "MRChem", + "mpi_processes": 1, + "nthreads": 1, + "routine": "mrchem.x", + "total_cores": 1, + "version": "1.2.0-alpha" + }, + "rsp_calculations": null, + "scf_calculation": { + "initial_energy": { + "E_ee": 1.220280348193976, + "E_eext": 0.0, + "E_el": -0.7847556322707149, + "E_en": -1.9145538948639653, + "E_kin": 0.9258603759485317, + "E_next": 0.0, + "E_nn": 0.0, + "E_nuc": 0.19090981420562303, + "E_tot": -0.5938458180650918, + "E_x": -0.15252084470528926, + "E_xc": -0.4885331476131019, + "Er_el": -0.37528846923086595, + "Er_nuc": 0.19090981420562303, + "Er_tot": -0.18437865502524292 + }, + "success": true + }, + "schema_name": "mrchem_output", + "schema_version": 1, + "success": true + } +} diff --git a/tests/h_pb/reference/h.out b/tests/h_pb/reference/h.out new file mode 100644 index 000000000..0cc156321 --- /dev/null +++ b/tests/h_pb/reference/h.out @@ -0,0 +1,290 @@ + + +*************************************************************************** +*** *** +*** *** +*** __ __ ____ ____ _ *** +*** | \/ | _ \ / ___| |__ ___ _ __ ___ *** +*** | |\/| | |_) | | | '_ \ / _ \ '_ ` _ \ *** +*** | | | | _ <| |___| | | | __/ | | | | | *** +*** |_| |_|_| \_\\____|_| |_|\___|_| |_| |_| *** +*** *** +*** VERSION 1.2.0-alpha *** +*** *** +*** Git branch unknown *** +*** Git commit hash unknown *** +*** Git commit author unknown *** +*** Git commit date unknown *** +*** *** +*** Contact: luca.frediani@uit.no *** +*** *** +*** Radovan Bast Magnar Bjorgve *** +*** Roberto Di Remigio Antoine Durdek *** +*** Luca Frediani Gabriel Gerez *** +*** Stig Rune Jensen Jonas Juselius *** +*** Rune Monstad Peter Wind *** +*** *** +*************************************************************************** + +--------------------------------------------------------------------------- + + MPI processes : (no bank) 1 + OpenMP threads : 1 + Total cores : 1 + +--------------------------------------------------------------------------- + +XCFun DFT library Copyright 2009-2020 Ulf Ekstrom and contributors. +See http://dftlibs.org/xcfun/ for more information. + +This is free software; see the source code for copying conditions. +There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. For details see the documentation. +Scientific users of this library should cite +U. Ekstrom, L. Visscher, R. Bast, A. J. Thorvaldsen and K. Ruud; +J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s + +--------------------------------------------------------------------------- + + MRCPP version : 1.6.0-alpha + Git branch : unknown + Git commit hash : unknown + Git commit author : unknown + Git commit date : unknown + + Linear algebra : EIGEN v3.4.0 + Parallelization : OpenMP (1 threads) + +--------------------------------------------------------------------------- + + + +=========================================================================== + MultiResolution Analysis +--------------------------------------------------------------------------- + polynomial order : 6 + polynomial type : Interpolating +--------------------------------------------------------------------------- + total boxes : 8 + boxes : [ 2 2 2 ] + unit lengths : [ 16.00000 16.00000 16.00000 ] + scaling factor : [ 1.00000 1.00000 1.00000 ] + lower bounds : [ -16.00000 -16.00000 -16.00000 ] + upper bounds : [ 16.00000 16.00000 16.00000 ] + total length : [ 32.00000 32.00000 32.00000 ] +=========================================================================== + + + +*************************************************************************** +*** *** +*** Initializing Molecule *** +*** *** +*************************************************************************** + + +=========================================================================== + Molecule +--------------------------------------------------------------------------- + Charge : -1 + Multiplicity : 1 +--------------------------------------------------------------------------- + N Atom : x y z +--------------------------------------------------------------------------- + 0 H : 0.000000 0.000000 0.000000 +--------------------------------------------------------------------------- + Center of mass : 0.000000 0.000000 0.000000 +=========================================================================== + + +=========================================================================== + Solvation Cavity +--------------------------------------------------------------------------- + Formulation : exponential + Dielectric constant : (in) 1.000000 + : (out) 78.400000 +--------------------------------------------------------------------------- + N R_0 Alpha Beta Sigma Radius x y z +--------------------------------------------------------------------------- + 0 2.6456 1.00 0.00 0.20 -> 2.6456 0.000000 0.000000 0.000000 +=========================================================================== + + +=========================================================================== + Square of the Debye-Huckel screening parameter +--------------------------------------------------------------------------- + Formulation : exponential + Screening function value: (in) 0.000000 + : (out) 0.087032 +--------------------------------------------------------------------------- + N R_0 Sigma Radius x y z +--------------------------------------------------------------------------- + 0 2.6456 0.20 -> 2.6456 0.000000 0.000000 0.000000 +=========================================================================== + + + +*************************************************************************** +*** *** +*** Computing Initial Guess Wavefunction *** +*** *** +*************************************************************************** + + +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Calculation : Compute initial orbitals + Method : Diagonalize SAD Hamiltonian + Precision : 1.00000e-03 + Screening : 1.20000e+01 StdDev + Restricted : True + Functional : LDA (SVWN5) + AO basis : 3-21G +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +=========================================================================== + Molecular Orbitals +--------------------------------------------------------------------------- + Alpha electrons : 1 + Beta electrons : 1 + Total electrons : 2 +--------------------------------------------------------------------------- + n Occ Spin : Norm +--------------------------------------------------------------------------- + 0 2 p : 9.999999944384e-01 +=========================================================================== + + +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + Calculation : Compute initial energy + Method : DFT (PBE0) + Relativity : None + Environment : None + External fields : None + Precision : 1.00000e-03 + Localization : Off +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +=========================================================================== + Molecular Energy (initial) +--------------------------------------------------------------------------- + Kinetic energy : (au) 0.925860375949 + E-N energy : (au) -1.914553894864 + Coulomb energy : (au) 1.220280348194 + Exchange energy : (au) -0.152520844705 + X-C energy : (au) -0.488533147613 + N-N energy : (au) 0.000000000000 +--------------------------------------------------------------------------- + Reaction energy (el) : (au) -0.375288469231 + Reaction energy (nuc) : (au) 0.190909814206 + Reaction energy (tot) : (au) -0.184378655025 +--------------------------------------------------------------------------- + Electronic energy : (au) -0.784755632271 + Nuclear energy : (au) 0.190909814206 +--------------------------------------------------------------------------- + Total energy : (au) -5.938458180651e-01 + : (kcal/mol) -3.726438769686e+02 + : (kJ/mol) -1.559141981237e+03 + : (eV) -1.615936792593e+01 +=========================================================================== + + +=========================================================================== + Orbital Energies (initial) +--------------------------------------------------------------------------- + n Occ Spin : Epsilon +--------------------------------------------------------------------------- + 0 2 p : (au) -0.116572346446 +--------------------------------------------------------------------------- + Sum occupied : (au) -0.233144692892 +=========================================================================== + + + +*************************************************************************** +*** *** +*** Printing Molecular Properties *** +*** *** +*************************************************************************** + + +=========================================================================== + Molecule +--------------------------------------------------------------------------- + Charge : -1 + Multiplicity : 1 +--------------------------------------------------------------------------- + N Atom : x y z +--------------------------------------------------------------------------- + 0 H : 0.000000 0.000000 0.000000 +--------------------------------------------------------------------------- + Center of mass : 0.000000 0.000000 0.000000 +=========================================================================== + + +=========================================================================== + Molecular Energy (final) +--------------------------------------------------------------------------- + Kinetic energy : (au) 0.925860375949 + E-N energy : (au) -1.914553894864 + Coulomb energy : (au) 1.220280348194 + Exchange energy : (au) -0.152520844705 + X-C energy : (au) -0.488533147613 + N-N energy : (au) 0.000000000000 +--------------------------------------------------------------------------- + Reaction energy (el) : (au) -0.375288469231 + Reaction energy (nuc) : (au) 0.190909814206 + Reaction energy (tot) : (au) -0.184378655025 +--------------------------------------------------------------------------- + Electronic energy : (au) -0.784755632271 + Nuclear energy : (au) 0.190909814206 +--------------------------------------------------------------------------- + Total energy : (au) -5.938458180651e-01 + : (kcal/mol) -3.726438769686e+02 + : (kJ/mol) -1.559141981237e+03 + : (eV) -1.615936792593e+01 +=========================================================================== + + +=========================================================================== + Orbital Energies (final) +--------------------------------------------------------------------------- + n Occ Spin : Epsilon +--------------------------------------------------------------------------- + 0 2 p : (au) -0.116572346446 +--------------------------------------------------------------------------- + Sum occupied : (au) -0.233144692892 +=========================================================================== + + +=========================================================================== + Dipole Moment (dip-1) +--------------------------------------------------------------------------- + r_O : 0.000000 0.000000 0.000000 +--------------------------------------------------------------------------- + Electronic vector : 0.000000 0.000000 0.000000 + Magnitude : (au) 0.000000 + : (Debye) 0.000000 +--------------------------------------------------------------------------- + Nuclear vector : -0.000000 -0.000000 -0.000000 + Magnitude : (au) 0.000000 + : (Debye) 0.000000 +--------------------------------------------------------------------------- + Total vector : -0.000000 -0.000000 -0.000000 + Magnitude : (au) 0.000000 + : (Debye) 0.000000 +=========================================================================== + + + + +*************************************************************************** +*** *** +*** Exiting MRChem *** +*** *** +*** Wall time : 0h 0m 51s *** +*** *** +*************************************************************************** + + diff --git a/tests/h_pb/test b/tests/h_pb/test new file mode 100755 index 000000000..8fe8c91be --- /dev/null +++ b/tests/h_pb/test @@ -0,0 +1,29 @@ +#!/usr/bin/env python3 + +import sys +from pathlib import Path + +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), + ER_TOT: rel_tolerance(1.0e-6), + ER_EL: rel_tolerance(1.0e-6), + ER_NUC: rel_tolerance(1.0e-6), +} + +ierr = run(options, input_file="h", filters=filters, extra_args=['--json']) + +sys.exit(ierr) diff --git a/tests/li_solv/li.inp b/tests/li_solv/li.inp index 0b3ad0f49..863f532da 100644 --- a/tests/li_solv/li.inp +++ b/tests/li_solv/li.inp @@ -21,11 +21,13 @@ "Cavity": { "spheres": "0 4.0 1.0 0.0 0.5" }, - "Permittivity": { - "epsilon_in": 1.0, - "formulation": "exponential", - "epsilon_out": { - "static": 2.0 + "Solvent":{ + "Permittivity": { + "epsilon_in": 1.0, + "formulation": "exponential", + "epsilon_out": { + "static": 2.0 + } } } }, diff --git a/tests/solventeffect/CMakeLists.txt b/tests/solventeffect/CMakeLists.txt index 93eabf748..c1a857595 100644 --- a/tests/solventeffect/CMakeLists.txt +++ b/tests/solventeffect/CMakeLists.txt @@ -1,7 +1,8 @@ target_sources(mrchem-tests PRIVATE - ${CMAKE_CURRENT_SOURCE_DIR}/reaction_operator.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/cavity_function.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/reaction_operator.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/cavity_function.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/PB_solver.cpp ) add_Catch_test( @@ -13,3 +14,8 @@ add_Catch_test( NAME reaction_operator LABELS "reaction_operator" ) + + add_Catch_test( + NAME PB_solver + LABELS "PB_solver" + ) diff --git a/tests/solventeffect/PB_solver.cpp b/tests/solventeffect/PB_solver.cpp new file mode 100644 index 000000000..a53ab9fe5 --- /dev/null +++ b/tests/solventeffect/PB_solver.cpp @@ -0,0 +1,169 @@ +/* + * 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 "catch.hpp" + +#include + +#include "mrchem.h" + +#include "analyticfunctions/HydrogenFunction.h" +#include "chemistry/Nucleus.h" +#include "chemistry/PeriodicTable.h" +#include "chemistry/chemistry_utils.h" +#include "environment/Cavity.h" +#include "environment/DHScreening.h" +#include "environment/GPESolver.h" +#include "environment/LPBESolver.h" +#include "environment/PBESolver.h" +#include "environment/Permittivity.h" +#include "qmfunctions/density_utils.h" +#include "qmfunctions/orbital_utils.h" +#include "qmfunctions/qmfunction_fwd.h" +#include "qmoperators/two_electron/ReactionOperator.h" + +using namespace mrchem; + +namespace PB_solver { +/* this test is the zeroth case of cases shown in 10.1137/18M119553X + * case 0 is just a positive charge, the rho_el should be zero and the energy should be (in the exact solution) -0.1373074208 Hartree + * this is meant to be lightweight test of the solver for theoretical correctness. + */ + +TEST_CASE("Poisson Boltzmann equation solver standard", "[PB_solver][pb_standard]") { + const double prec = 1.0e-5; + const double thrs = 1.0e-6; + + auto dyn_thrs = false; + auto kain = 7; + auto max_iter = 100; + auto eps_in = 1.0; + auto eps_out = 78.54; + auto kappa_out = 0.054995; + auto slope = 0.2; + + auto R = std::vector({3.7794522509156563}); + auto sph_coords = std::vector>({{0.0, 0.0, 0.0}}); + // initialize spherical cavity + auto sphere = std::make_shared(sph_coords, R, slope); + // initialize dielectric function + auto dielectric_func = Permittivity(sphere, eps_in, eps_out, "exponential"); + // initialize DHScreening object + auto kappa_sq = DHScreening(sphere, kappa_out, "continuous"); + + auto P_p = std::make_shared(*MRA, prec); + auto D_p = std::make_shared>(*MRA, 0.0, 0.0); + + PeriodicTable PT; + SECTION("case 0: one positive charge in center of sphere (born model)", "[PB_solver][pb_standard][case_0]") { + // initialize Nuclei in center of sphere + auto q_coords = std::vector>({{0.0, 0.0, 0.0}}); + Nucleus Q(PT.getElement(0), q_coords[0]); + Nuclei molecule; + molecule.push_back(Q); + + // initialize orbital vector + auto Phi_p = std::make_shared(); + auto &Phi = *Phi_p; + Phi.push_back(Orbital(SPIN::Paired)); + Phi.distribute(); + + HydrogenFunction f(1, 0, 0); + if (mrcpp::mpi::my_orb(Phi[0])) mrcpp::cplxfunc::project(Phi[0], f, NUMBER::Real, prec); + + auto rho_nuc = chemistry::compute_nuclear_density(prec, molecule, 100); + + auto scrf_p = std::make_unique(dielectric_func, kappa_sq, rho_nuc, P_p, D_p, kain, max_iter, dyn_thrs, SCRFDensityType::NUCLEAR); + auto Reo = std::make_shared(std::move(scrf_p), Phi_p); + Reo->setup(prec * 10); + + Density rho_el(false); + + auto [Er_el, Er_nuc] = Reo->getSolver()->computeEnergies(rho_el); + + Reo->clear(); + REQUIRE((Er_nuc) == Approx(-1.358726143734e-01).epsilon(thrs)); // exact is -0.1373074208 Hartree, though ours is close, i think we are a bit too far away, some parameterization issue + } +} + +TEST_CASE("Poisson Boltzmann equation solver linearized", "[PB_solver][pb_linearized]") { + const double prec = 1.0e-5; + const double thrs = 1.0e-6; + + auto dyn_thrs = false; + auto kain = 5; + auto max_iter = 200; + auto eps_in = 1.0; + auto eps_out = 78.54; + auto kappa_out = 0.054995; + auto slope = 0.2; + + auto R = std::vector({3.7794522509156563}); + auto sph_coords = std::vector>({{0.0, 0.0, 0.0}}); + // initialize spherical cavity + auto sphere = std::make_shared(sph_coords, R, slope); + // initialize dielectric function + auto dielectric_func = Permittivity(sphere, eps_in, eps_out, "exponential"); + // initialize DHScreening object + auto kappa_sq = DHScreening(sphere, kappa_out, "continuous"); + + auto P_p = std::make_shared(*MRA, prec); + auto D_p = std::make_shared>(*MRA, 0.0, 0.0); + + PeriodicTable PT; + SECTION("case 0: one positive charge in center of sphere (born model)", "[PB_solver][pb_linearized][case_0]") { + + // initialize Nuclei in center of sphere + auto q_coords = std::vector>({{0.0, 0.0, 0.0}}); + Nucleus Q(PT.getElement(0), q_coords[0]); + Nuclei molecule; + molecule.push_back(Q); + + // initialize orbital vector + auto Phi_p = std::make_shared(); + auto &Phi = *Phi_p; + Phi.push_back(Orbital(SPIN::Paired)); + Phi.distribute(); + + HydrogenFunction f(1, 0, 0); + if (mrcpp::mpi::my_orb(Phi[0])) mrcpp::cplxfunc::project(Phi[0], f, NUMBER::Real, prec); + + auto rho_nuc = chemistry::compute_nuclear_density(prec, molecule, 100); + + auto scrf_p = std::make_unique(dielectric_func, kappa_sq, rho_nuc, P_p, D_p, kain, max_iter, dyn_thrs, SCRFDensityType::NUCLEAR); + + auto Reo = std::make_shared(std::move(scrf_p), Phi_p); + Reo->setup(prec * 10); + + Density rho_el(false); + + auto [Er_el, Er_nuc] = Reo->getSolver()->computeEnergies(rho_el); + + Reo->clear(); + REQUIRE(Er_nuc == Approx(-1.358725427728e-01).epsilon(thrs)); // what we get in standard GPESolver is -1.455145361712e-01, while with PB we get -1.329978908155e-01 + } +} + +} // namespace PB_solver diff --git a/tests/solventeffect/reaction_operator.cpp b/tests/solventeffect/reaction_operator.cpp index 40dfdda49..323b579a0 100644 --- a/tests/solventeffect/reaction_operator.cpp +++ b/tests/solventeffect/reaction_operator.cpp @@ -40,8 +40,8 @@ #include "chemistry/PeriodicTable.h" #include "chemistry/chemistry_utils.h" #include "environment/Cavity.h" +#include "environment/GPESolver.h" #include "environment/Permittivity.h" -#include "environment/SCRF.h" #include "qmfunctions/Orbital.h" #include "qmfunctions/density_utils.h" #include "qmfunctions/orbital_utils.h" @@ -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, SCRFDensityType::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); @@ -99,7 +99,7 @@ TEST_CASE("ReactionOperator", "[reaction_operator]") { density::compute(prec, rho_el, Phi, DensityType::Total); rho_el.rescale(-1.0); - auto [Er_nuc, Er_el] = Reo->getHelper()->computeEnergies(rho_el); + auto [Er_nuc, Er_el] = Reo->getSolver()->computeEnergies(rho_el); auto total_energy = Er_nuc + Er_el; Reo->clear(); REQUIRE(total_energy == Approx(-1.022729683846e-01).epsilon(thrs));