3.18 Continuum Solvation Methods

The simulation of a solvent in a quantum mechanical calculation can, in principle, be done in two ways. One way is to include explicit solvent molecules in the calculation. This straightforward approach usually requires molecular dynamics (MD) simulations in order to yield thermodynamically meaningful observables as e.g. solvation free energies.

The second way is to average the effect of the solvent and treat it as a continuum which responds to the electrostatic potential created by the solute, i.e. the entity that is to be solvated. There are several flavors to this comparatively inexpensive approximation, e.g. the polarizable continuum model (PCM) [218], the conductor like screening model (COSMO) [166], the self-consistent continuum solvation (SCCS) model [11], the “SMx” models [67, 210, 211], or CMIRSv1.1 [322] to name some of the more popular ones. Statistical sampling then only needs to be performed for the degrees of freedom of the solute which obviously makes it computationally much cheaper.

In general, the necessary integration of the solvent’s degrees of freedom beforehand leads to a problem where one now needs to solve a generalized Poisson’s equation,

\displaystyle\nabla\left(\varepsilon_{0}\varepsilon(\bm{r})\nabla\Phi(\bm{r})% \right)=-4\pi\varrho(\bm{r}),

(3.49)

to obtain the electrostatic potential $\Phi(\bm{r})$ created by the total charge density $\varrho(\bm{r})$ which now accounts for the electrostatic polarization potential of the solvent (often called “reaction field”). Eq. 3.49 contains a spatially dependent dielectric permittivity function $\varepsilon(\bm{r})$ in contrast to the regular Poisson’s equation,

\displaystyle\nabla\left(\varepsilon_{0}\nabla\Phi_{\mathrm{H}}(\bm{r})\right)% =-4\pi\varrho(\bm{r}),

(3.50)

which is solved in a regular DFT calculation in every step of the SCF cycle to get the Hartree potential $\Phi_{\mathrm{H}}$ . $\varepsilon(\bm{r})$ is defined in such a way that it has the experimentally accessible bulk value of the solvent at large distance to the solute. Inside the region in space occupied by the solute, the so-called ‘cavity’, it has a value of 1, corresponding to vacuum permittivity. Any polarization inside that region is assumed to be accounted for at the DFT level. The definition of the cavity and the transition of $\varepsilon(\bm{r})$ at its boundary differ between different implicit solvation models.

Currently, FHI-aims intrinsically supports three models which have different strengths and capabilities. The Multipole Expansion (MPE) implicit solvation method uses a sharp, step-like transition of $\varepsilon(\bm{r})$ at the cavity boundary. It separates the FHI-aims grid into two (original MPE method[281]) or multiple (MPE- $n$ c[91], MPE- $n$ cps) domains and couples them via electrostatic boundary conditions. It therefore in fact solves multiple coupled Poisson equations with different dielectric permittivities. The Finite ion-size and Stern layer modified Poisson-Boltzmann (SMPB) method solves a single Poisson equation on the full FHI-aims integration grid by defining a smooth dielectric permittivity function. On top of this ion-free implicit solvation model, the SMPB approach also supports the modeling of finite ionic strengths in the solution. The available features for MPE and SMPE are summarized in Table 3.1. The conductor-like screening model (COSMO) initially approximates the solvent as a conductor and later applies a correction to account for the finite permittivity of the actual solvent.

In addition to these models implemented directly in FHI-aims, we also provide a communication interface[90] for the continuum solvation package Environ. Environ implements a wide range of continuum embedding methods for molecules and surfaces, most notably the self-consistent continuum solvation (SCCS) [11] model (on which SMPB is based) and the soft-sphere continuum solvation (SSCS) model.[92] The full documentation of Environ can be found at (https://environ.readthedocs.io).

The MPE solvation model the faster compared to SMPB with only a small overhead with respect to vacuum calculations. The overhead of implicit solvation methods is reduced, when performing expensive hybrid calculations, since the actual time for the implicit solvation calculations does not vary with the functional. The overhead of Environ is typically higher than that of the intrinsically supported methods for isolated molecules, as its Fast Fourier Transform based formalism requires the formal introduction of a periodic cell (although corrections which remove interactions with periodic images are available) which is filled with a regular grid. Essentially, this means that in Environ, a significant amount of computational time is spent on regions in space that are free of charge, which is not necessary in FHI-aims. For periodic surfaces where a larger portion of the cell is filled with charge, the relative overhead of Environ becomes smaller.

	MPE	SMPB
Solvent Parametrizations	H₂O (N,C)	H₂O (N,C), CH₃OH (N,C)
		(more in work)
Dissolved ions/salt	no	yes (SMPB/LPB)
Salt Parametrizations	–	Aq. Monoval. Salt Solutions
CPU speed	fast	moderate
Forces	no (in work)	yes
PBCs	yes	no (in work)
Developers	Markus Sinstein	Stefan Ringe
	Jakob Filser, Pavel Stishenko	Christoph Muschielok, Marvin H. Lechner

Table 3.1: Comparison of the two implicit solvation methods in FHI-aims. Parameter sets for the single cavity MPE method are available in ref. [281], for MPE-

n

c in ref. [91], for the SMPB method in ref. [11] and [81] (parameters for methanol as more solvents in current work). N and C indicate parameter sets fitted for neutral and charged solutes, respectively.

In the following, both models are summarized and the key input parameters presented.

3.18.1 MPE Implicit Solvent Model

The periodic systems are supported only by the MPE- $n$ cps flavour.

The current implementation does not have analytical forces yet.

Generally, when combining MPE with other functionality, you should know what you are doing. No specific interactions with other methods beyond single point DFT are implemented, so only methods which do not interfere with MPE are safe to use.

Default behaviour: If solvent mpe-nc is set, an MPE- $n$ c calculation in water will be performed with tight expansion orders and isodensity and non-electrostatic parameter optimized for the respective xc functional (or optimized for HSE06, if no parameterization for the chosen functional exists), see ref. [91], table 2. If solvent mpe2017 is set and the solute charge is 0 or +1, then the solvation parameters for neutrals and cations (SPANC) are used.[281]

As outlined in more detail in Ref. [281], the multipole expansion (MPE) implicit solvent model offers an efficient way of solving Eq. 3.49 based on the knowledge of the Hartree potential $\Phi_{\mathrm{H}}$ readily available from a splined representation in FHI-aims (cf. Sec. 3.7) via least-squares fitting instead of integration. The central Ansatz is

\Phi(\bm{r})=\varepsilon^{-1}(\bm{r})\Phi_{\text{H}}(\bm{r})+\Phi_{\text{MPE}}% (\bm{r})

(3.51)

With the potential that the solute would generate in vacuum $\Phi_{\text{H}}(\bm{r})$ and the inverse permittivity $\varepsilon^{-1}(\bm{r})$ . The remaining part $\Phi_{\text{MPE}}(\bm{r})$ of the potential is defined separately for the cavity and the dielectric region. The dielectric function $\varepsilon(\bm{r})$ here needs to be a step-function in 3D-space where the following boundary conditions apply at the step ( $\bm{n}$ denotes the normal direction to the interface):


	$\displaystyle\Phi_{+}=\Phi_{-}$		(3.52a)

	$\displaystyle\bm{n}\cdot\varepsilon_{+}\nabla\Phi_{+}=\bm{n}\cdot\varepsilon_{% -}\nabla\Phi_{-}$		(3.52b)

Then, the above equations are discretized in two ways:

•

The potentials $\Phi_{\text{MPE}}$ in- and outside the cavity are expressed in a truncated multipole series with expansion orders $l_{\mathrm{max,R}}$ and $l_{\mathrm{max,O}}$ , and
•

equations 3.52a and 3.52b are evaluated at $N$ points on the interface manifold.

Thereby, $N$ is chosen such that the resulting system of linear equations (SLE) is overdetermined (typically by a factor of two to three).

MPE- $n$ c additionally separates the cavity into atom-centered subregions, each with its own multipole series. This ensures fast convergence of the total (DFT) energy with expansion order. In contrast to the earlier single-cavity version, it does not rely on a specific kind of solute size dependent error cancellation, improving transferability (although different kinds of error cancellation are likely still present).[91] Although slightly more costly than the single-cavity version, MPE- $n$ c is typically not the computational bottleneck of any calculation. It is therefore more advisable in most cases.

One noteable exception may be small (fewer than $\approx 10$ non-hydrogen atoms) neutral or cationic solutes, for which a parameterization of single-cavity MPE with isc_cavity_type rho_multipole_static exists with an overall smaller reported error than any of the MPE- $n$ c parameter sets.[281] Note that this cavity type was simply never parameterized for MPE- $n$ c, because with this cavity definition, parameters are not transferable to anionic solutes. If desired, a parameterization of MPE- $n$ c with isc_cavity_type rho_multipole_static for neutral and cationic solutes is certainly possible and would likely be at least equally as accurate as the single-cavity MPE version.

MPE- $n$ cps uses piecewise representation for both solute cavity and for the solvent region (Piecewise Solvent). The solvent region separation is based on the vicinity to heavy (non-hydorgen) atoms, in the same way as it is done for MPE- $n$ c method. The Voronoi tesselation algorithm implemented in the Voro++ library [268] by Chris H. Rycroft is used for solvent region separation. Piecewise solvent representation allows simulation of systems with periodic boundary conditions and may accelerate calculations for large systems. Only orthogonal lattices are currently supported. Large Ewald_radius may be necessary - do convergence testing.

If needed for this special case, or for consistency with previous work, the original implementation can still be accessed by solvent mpe2017.

Tags for general subsection of control.in:

The keywords controlling the MPE module are divided into six categories:

elementary

These are the most important keywords which define the physical model.
convergence

Here, the most important convergence parameters are collected which should be checked before doing (large scale) production runs.
performance

These settings should only influence computational time and memory requirements.
expert

These settings should only be modified by an experienced user as they allow quite profound modifications.
debug

Debug settings are intended to give valuable insight for developers into intermediate results.
deprecated

These settings are no longer actively maintained and should not be used except by developers who know what is actually implemented.

elementary

Tag: solvent(control.in)

Usage: solvent method
Purpose: Specifies the desired implicit solvent model.
Options: method is a string which specifies the implicit solvent method

•

mpe-ncps: MPE- $n$ cps is an extension of MPE- $n$ c that uses Piecewise Solvent representation to support periodic systems. Only orthogonal lattices currently supported. Large Ewald_radius may be necessary - do convergence testing..
•

mpe-nc: MPE- $n$ c[91]
•

mpe2017: Single cavity MPE.[281] Note that this is not entirely equivalent to the deprecated option mpe, as the defaults for mpe_solvent_permittivity, isc_cavity_type and mpe_nonelectrostatic_model have changed.
•

mpb (cf. Sec. 3.18.2)
•

cosmo (cf. Sec. 3.18.3)
•

Environ (cf. Sec. 3.18.4)

All further options described in this section apply to all mpe-ncps, mpe-nc and mpe2017. Different choices of isc_cavity_type and mpe_nonelectrostatic_model are required for the two models.

Tag: mpe_solvent_permittivity(control.in)

Usage: mpe_solvent_permittivity epsilon
Purpose: Specifies the dielectric constant of the bulk solvent.
epsilon is a positive real number equal to the macroscopic dielectric constant of the solvent.
Default: 78.36 (water)

Tag: isc_cavity_type(control.in)

Usage: isc_cavity_type type
Purpose: This keyword controls the model used to sample the implicit solvent cavity for the MPE method. Depending on type, further flags (or even lines) might be necessary. Those are explained below.
Options: Currently supported options are rho_free, rho_multipole_static, overlapping_spheres, and rho_multipole_dynamic. Parameterizations exist only for the first two options. The latter two are thus described in the expert section.
Default: If solvent mpe-nc, then rho_free with isodensity optimized for water and the chosen xc functional, cf. ref. [91], table 2.
If solvent mpe2017, and the solute charge is 0 or +1, then rho_multipole_static with the solvation parameters for neutrals and cations (SPANC).[281]

isc_cavity_type sub-tag: rho_free(control.in)

Usage: isc_cavity_type rho_free rho_iso
Purpose: Constructs the cavity as an iso-density surface of the superposed electron density of the neutral, free atoms in the solute.
rho_iso is a positive real number specifying the desired iso-density value in units of e/Å³ .

isc_cavity_type sub-tag: rho_multipole_static(control.in)

Usage: isc_cavity_type rho_multipole_static rho_iso
Purpose: Constructs the cavity as an iso-density surface of the (multipole-expanded) converged electron density of the solute in vacuum, keeping the cavity static throughout the actual MPE calculation.

First, a vacuum calculation (SCF with any MPE related keywords turned off) will be performed until self-consistency is achieved, then SCF will be restarted, with specified MPE related keywords turned on and the isodensity cavity constructed from the converged electron density in vacuum.

rho_iso is a positive real number specifying the desired iso-density value in units of e/Å³ .

Tag: mpe_nonelectrostatic_model(control.in)

Usage: mpe_nonelectrostatic_model model
Purpose: This keyword controls any additional, “non-electrostatic” term not included in the purely electrostatic treatment of the solvent.
Options:

•

linear_OV: see below
•

none: compute no non-electrostatic terms

Default: If solvent mpe-nc, then linear_OV with $\beta$ =0 and $\alpha$ optimized for water and the chosen xc functional, cf. ref. [91], table 2.
If solvent mpe2017, and the solute charge is 0 or +1, then the solvation parameters for neutrals and cations (SPANC).[281]

mpe_nonelectrostatic_model sub-tag: linear_OV(control.in)

Usage: mpe_nonelectrostatic_model linear_OV $\alpha$ $\beta$
Purpose: Corrects the total energy term by $\alpha O+\beta V$ where $O$ is the surface area of the cavity and $V$ its volume.
$\alpha$ is a real number in units of eV/Å².
$\beta$ is a real number in units of eV/Å³.

This non-electrostatic model is in principle identical to the one proposed by Andreussi et al. [11]. Note, however, that the surface tension of the solvent is here included in the parameter $\alpha$ .

convergence

Tag: mpe_lmax_rf(control.in)

Usage: mpe_lmax_rf lmax
Purpose: Specifies the expansion order of the polarization potential aka reaction field inside the cavity.
lmax is a non-negative integer number.
Default: 8

This is a critical convergence parameter of the MPE model. For MPE- $n$ c, a value of 8 should be sufficent in most cases. 12 is a safe choice. If MPE ever becomes the bottleneck of the calculation (which usually isn’t the case), then 6 can be used for a “quick and dirty” approximate calculation.[91]

solvent mpe2017 generally needs higher expansion orders for larger solutes, but also relies to some degree on error cancellation of an underconverged electrostatic energy with an overestimated non-electrostatic attraction, as discussed in ref. [91].

Tag: mpe_lmax_ep(control.in)

Usage: mpe_lmax_ep lmax
Purpose: Specifies the expansion order of the polarization potential aka reaction field outside of the cavity.
lmax is a non-negative integer number.
Default: With solvent mpe-nc, default is 6. With solvent mpe2017, default is maximum value of l_hartree for all species.

This parameter is similar to but usually less critical than mpe_lmax_rf. For MPE- $n$ c, a value of 6 should be sufficent in most cases. 8 is a safe choice. If MPE ever becomes the bottleneck of the calculation (which usually isn’t the case), then 4 can be used for a “quick and dirty” approximate calculation.[91]

Note that mpe2017 places expansion centers for the external potentials on solute hydrogen atoms, whereas MPE- $n$ c does not. Both place centers on non-hydrogen atoms.

Tag: mpe_degree_of_determination(control.in)

Usage: mpe_degree_of_determination dod
Purpose: Defines the desired ratio of number of rows to columns in left-hand side matrix of the MPE equation.
dod is a real number $\geq 1.0$ .
Default: 5.0

For the very limited (!) number of applications of the MPE method so far, the default value of 5.0 has been a save choice. Note, that the requested degree of determination can only very approximately be reached. This can lead to an under- determination of the MPE equations and a subsequent termination of the program when values for dod very close to 1 are chosen.

performance

Tag: mpe_solver(control.in)

Usage: mpe_solver type
Purpose: Defines the numerical method used to solver the MPE equations.
Options: type can be chosen from: direct, lsqr.
Default: direct

The default behavior is to factorize the left hand side of the MPE equations using any of the methods described in mpe_factorization_type. This allows to robustly solve the MPE equation via the pseudo-inverse of the left-hand side.
lsqr uses the iterative LSQR algorithm with a sparse representation of the left hand side matrix. This only makes sense together with solvent mpe-nc, as there will be no sparsity otherwise. While the direct solvers are generally faster, LSQR requires less memory and can be used as a fallback option for very large systems. Note also the caveat at mpe_factorization_type.

Tag: mpe_factorization_type(control.in)

Usage: mpe_factorization_type type
Purpose: Defines the numerical method used to factorize the left-hand side of the MPE equations as the first step to the numerical solution. Has no effect if mpe_solver lsqr is chosen.
Options: type can be chosen from: qr, qr+svd, and svd.
Default: If solvent mpe-nc is set, and mpe_lmax_rf is 8 or smaller and mpe_lmax_ep is 6 or smaller, then the default is to perform a QR factorization and solve the linear system directly from there.

If solvent mpe2017 is used, or if one or both of the two expansion orders is higher than the above values, the default behavior is to perform a QR factorization with a singular value decomposition (SVD) on top. This allows to robustly solve the MPE equation via the pseudo-inverse of the left-hand side. This is the safest option, and sometimes necessary to deal with ill-conditioning at high expansion orders.

qr is generally significantly faster than the other options. However:

Caveat! Using the (non rank-revealing) QR factorization alone can fail when the left-hand side is rank deficient which can easily happen—especially for large expansion orders mpe_lmax_rf and/or mpe_lmax_ep! On the other hand, svd does not necessarily mean that no QR factorization is performed as this (at least for the parallel implementation) depends on the (Sca)LAPACK driver routine used.

Tag: mpe_sc_block_size(control.in)

Usage: mpe_sc_block_size size
Purpose: Defines the ScaLAPACK block size used for solving the linear system of MPE. This may have an influence on computational time in some cases.
size is a positive integer number. If chosen too large to distribute the matrix over the given number of processes, the size will be reduced automatically.
Default: With mpe_solver direct: 32
With mpe_solver lsqr: 16

Tag: mpe_tol_load_imbalance(control.in)

Usage: mpe_tol_load_imbalance tol
Purpose: Works only together with mpe_solver lsqr. Defines how uneven the matrix elements are allowed to be distributed over processes. Example: If the MPE matrix has 1000 non-zero elements, and you are running on 10 cores, then it would ideally be possible to give exactly 100 elements to each process. With mpe_tol_load_imbalance 1.05, each process can be given up to 105 elements.
tol is a real number > 1.
Default: 1.05

Tag: mpe_tol_overdistribution(control.in)

Usage: mpe_tol_overdistribution tol
Purpose: Works only together with mpe_solver lsqr. Defines the factor by which a matrix block can be distributed over more cores than would ideally be necessary. Example: If the MPE matrix has 1000 non-zero elements, and you are running on 10 cores, then a matrix block with 200 elements could ideally be distributed over 2 cores. With mpe_tol_overdistribution 4, it can be distributed over up to 8 cores.
tol is a real number > 1.
Default: 4.0

The load balancing of mpe_solver lsqr tries to distribute matrix blocks over processes in such a way that the maximum number of matrix elements on any process is minimized while not distributing each individual block over more processes than is necessary. If no compromise is found, both mpe_tol_load_imbalance and mpe_tol_overdistribution will automatically be gradually increased until one is found.

expert

isc_cavity_type sub-tag: rho_multipole_dynamic(control.in)

Usage: isc_cavity_type rho_multipole_dynamic rho_iso
Purpose: Constructs the cavity as an iso-density surface of the self-consistent (multipole-expanded) electron density of the solute, updating the cavity in each SCF step.
rho_iso is a positive real number specifying the desired iso-density value in units of e/Å³ .

With this method, the cavity is updated to the current electron density in every SCF step. This also means that the MPE equations have to be solved in every SCF step making it computationally more expensive.

Caveat: Although it should in principle be implemented, the combination of this option together with solvent mpe-nc is not being actively maintained.

isc_cavity_type sub-tag: overlapping_spheres(control.in)

Usage: isc_cavity_type overlapping_spheres type value
Purpose: Constructs the cavity as a superposition of overlapping spheres around all atoms.
type specifies how the radii of the atomic spheres are determined.

•

radius: all spheres have the same radius given by value in units of Å;
•

rho: the atomic spheres are iso-density surfaces based on the electron density of the isolated, neutral atom with an iso-value of value in units of e/Å³ .

value is a real number whose meaning and units depend on the choice of radius (see above).

WARNING: The usage of this cavity type is strongly discouraged! It has been helpful in the development to analyze the cavity sampling process itself. The resulting cavities, however, are almost certainly not smooth and were never intended to be used in production calculations. When used with the MPE model, the whole calculation is prone to numerical problems and the results are very often unphysical. Instead, use the rho_free type that builds the cavity based on the superposition of atomic densities (which is again smooth) or use other types based on the (self-consistent) electron density of the solute.

Tag: mpe_tol_adjR2(control.in)

Usage: mpe_tol_adjR2 tol
Purpose: Defines the tolerance for the adjusted coefficient of determination $\bar{R}^{2}$ of the solved MPE equations. Will abort if $\bar{R}^{2}<1-tol$ .
tol is a real number between $0.0$ and $1.0$ . Default: 0.05

The MPE equations are sometimes not solvable in the regular solid harmonic basis used for the reaction field. This is the case especially for large molecules. A low $\bar{R}^{2}$ indicates such a bad solution. In some cases increasing mpe_lmax_rf helps, but there are pathologic cases where increasing mpe_lmax_rf leads to a perpetual decrease in $\Delta_{\text{solv}}^{\text{el}}G$ without ever converging.

For $\bar{R}^{2}<0.95$ it is likely that less than 90% of $\Delta_{\text{solv}}^{\text{el}}G$ are captured. This is however based on experience from a limited number of cases. Feedback to the developers (filser@fhi-berlin.mpg.de) will be appreciated!

Tag: mpe_tol_adjR2_wait_scf(control.in)

Usage: mpe_tol_adjR2_wait_scf bool
Purpose: If .true., will wait until the SCF cycle is converged before it is checked whether $\bar{R}^{2}<1-tol$ .
Default: .false.

Although MPE does not actively try to converge, $\bar{R}^{2}$ tends to improve during the SCF procedure. Setting mpe_tol_adjR2_wait_scf can thus help borderline cases converge, at the cost of spending the full computation time of the SCF procedure on a calculation that might ultimately fail.

Tag: isc_surface_curvature_correction(control.in)

Usage: isc_surface_curvature_correction bool
Purpose: When this flag is turned on, the calculated surface area (and volume) of the cavity is approximately corrected for the cavity curvature.
bool is of Boolean type. Default: .true.

The effect of this keyword is usually rather negligible. For more details regarding the correction, please consult Ref. [281].

Tag: isc_rho_rel_deviation_threshold(control.in)

Usage: isc_rho_rel_deviation_threshold threshold
Purpose: Defines the convergence criterion of the cavity generation process: The walker dynamics simulation is run until the density values for all walkers deviate from the chosen iso value by at most threshold.
threshold is a small, positive real number. Default: $10^{-3}$

This keyword is only applicable for an isc_cavity_type defined by an iso-density value.

Tag: isc_max_dyn_steps(control.in)

Usage: isc_max_dyn_steps num
Purpose: Determines the maximum number of allowed steps to reach convergence of the walker dynamics simulation in the cavity creation process.
num is a positive integer number. Default: 300

Tag: isc_try_restore_convergence(control.in)

Usage: isc_try_restore_convergence bool
Purpose: When convergence of the cavity creation dynamics run could not be achieved within the number of allowed steps specified by isc_max_dyn_steps, this flag allows to enforce convergence by simply deleting all walkers not satisfying the convergence criterion given by isc_rho_rel_deviation_threshold.
bool is of Boolean type. Default: .false.

Although a simple check is done to stop the calculation when too many walkers do not satisfy the convergence criterion, one should always manually checking the resulting cavity for larger holes that might result from the deletion of walkers which can lead to a bad estimate of the cavity’s surface area and volume and maybe also have an impact on the quality of the polarization potential.

Tag: isc_kill_ratio(control.in)

Usage: isc_kill_ratio fraction
Purpose: This keyword can be helpful when the walker dynamics run does not converge due to trapped walkers by killing the worst fraction of walkers at each neighbor list update step (also see isc_update_nlist_interval).
fraction is a non-negative real number much smaller than 1. Default: 0.0

As the number of possibly trapped walkers depends a lot on the shape of the electron density, it is rather difficult to give a recommendation about a sensible value for fraction. In case walkers get stuck, we propose to use a rather conservative kill ratio of $\sim 10^{-3}$ and only increase it if necessary.

Tag: isc_update_nlist_interval(control.in)

Usage: isc_update_nlist_interval num
Purpose: This keyword triggers a re-evaluation of the neighbor lists in the density walkers dynamics simulation after every num steps.
num is a positive integer number. Default: 50

Tag: isc_dynamics_friction(control.in)

Usage: isc_dynamics_friction fric
Purpose: The value of fric determines how much “kinetic energy” is removed from the walkers in every step of the simulations via a simple velocity scaling.
fric is a real number between 0 and 1. Default: 0.1

A value of 0 for fric means that no energy is removed from the system which may lead to a bad convergence behavior. On the other hand, a value of 1 means that all kinetic energy is removed at each step which tends to slow down the rate of convergence.

Tag: isc_dt(control.in)

Usage: isc_dt delta
Purpose: Determines the “time” step of the walker dynamics simulation.
delta is a positive real number. Default: 0.1

Note: Since this is no actual physical quantity, arbitrary time units are used.

Tag: isc_rho_k(control.in)

Usage: isc_rho_k k
Purpose: Determines the force constant k for the “density” force, i.e. the harmonic force that pulls the walkers along the density gradient to the specified iso-density value.
k is a positive real number. Default: 1.0

Note: Since this is no actual physical quantity, arbitrary time units are used.

Tag: isc_rep_k(control.in)

Usage: isc_rep_k k
Purpose: Determines the force constant k for the repulsive interaction between walkers perpendicular to the density gradient.
k is a positive real number. Default: 0.01

Note: Since this is no actual physical quantity, arbitrary time units are used.

Tag: isc_g_k(control.in)

Usage: isc_g_k k
Purpose: Determines the force constant k for the “gravitational” force that drags walkers to the center of gravity of the solute in case the local density gradient is too small.
k is a positive real number. Default: 2.0

Usually this should not happen, but when walkers move too far away from the solute, the density gradient becomes very small and its direction is unreliable due to numerical noise (see isc_gradient_threshold). In this case, the walker is dragged to the center of the solute until the density gradient is again large enough. Note: Since this is no actual physical quantity, arbitrary time units are used.

Tag: isc_gradient_threshold(control.in)

Usage: isc_gradient_threshold thsq
Purpose: When the squared norm of the electron density gradient at the position of a walker is less than thsq, this gradient is considered unreliable. Instead, a simple “gravitational” force towards the center of the solute is applied.
thsq is a positive real number. Default: $10^{-8}$

The force constant of the “gravitational” force is determined by isc_g_k.

debug

Tag: mpe_timing_output_level(control.in)

Usage: mpe_timing_output_level int
Purpose: Defines how deep timing analysis of the MPE method will regress into subroutines. Default: 2

Tag: isc_cavity_restart(control.in)

Usage: isc_cavity_restart filename
Purpose: Read the solvation cavity from restart file (if available) and write new cavity to same file.
filename is the name of the restart file.

Specifying this keyword is almost equivalent to specifying both isc_cavity_restart_read and isc_cavity_restart_write with the same filename option except that with isc_cavity_restart the program does not abort when there is no restart file to read from.
Note: This keyword is intended for debugging purposes. Do not rely on the current structure of the cavity restart file as it might change in the future.

Tag: isc_cavity_restart_read(control.in)

Usage: isc_cavity_restart_read filename
Purpose: Read the solvation cavity from restart file instead of constructing a new one.
filename is the name of the file (in .xyz format) containing the cavity points and normal vectors. Additionally, a file <filename>.bin is written which contains the entire cavity information in not human-readable form. While the primary purpose of the former is visualization, the latter is the actual restart file.

Note: This keyword is intended for debugging purposes. Do not rely on the current structure of the cavity restart file as it might change in the future.

Tag: isc_cavity_restart_write(control.in)

Usage: isc_cavity_restart_write filename
Purpose: Write the cavity to the specified restart file once created.
filename is the name of the restart file (in .xyz format). If additionally <filename>.bin is present, the cavity is read from the latter instead. Note that the .xyz file has to be present in both cases. While this file itself is sufficient to create a cavity, only the .bin file allows for a fully deterministic restart.

Note: This keyword is intended for debugging purposes. Do not rely on the current structure of the cavity restart file as it might change in the future.

Tag: isc_record_cavity_creation(control.in)

Usage: isc_record_cavity_creation filename num
Purpose: Controls the output of snapshots during the cavity generation process. When num is positive, every num steps an XYZ snapshot of the cavity is written to file filename. For other choices of num, no output will be generated.
filename is of type string.
num is of type integer. Default: 0

This keyword is intended for debugging purposes. Note that the size of the output file can become very large!

deprecated

Tag: mpe_n_centers_ep(control.in)

Usage: mpe_n_centers_ep n
Purpose: Defines the number of centers used for the expansion of the polarization potential outside of the cavity.
n is a positive integer number. Default: number of centers for the Hartree potential expansion (cf. 3.7)

The first n centers defined in geometry.in are used as expansion centers. The default is to use all of them. Only change this value if you fully understand what you are doing and why you want to do this!

Tag: mpe_n_boundary_conditions(control.in)

Usage: mpe_n_boundary_conditions nbc
Purpose: Determines the number of boundary conditions imposed at every point on the cavity interface.
Valid choices for nbc are 2 and 4. Default: 2

As outlined in Ref. [281], there are at least two more boundary conditions other than Eqns. 3.52a and 3.52b that can be imposed on the electrostatic potential / field / flux at a dielectric interface. The default is to enforce continuity of the potential and continuity of the dielectric flux perpendicular to the interface, i.e. nbc equals 2. Furthermore, continuity of the electric field parallel to the interface can be imposed, i.e. nbc equals 4. However, this should automatically be satisfied by the former two boundary conditions and—in the best case—only leads to a higher order correction of the fit. Warning: The non-default has not been tested thoroughly. Verify your results carefully when using it!

Tag: isc_calculate_surface_and_volume(control.in)

Usage: isc_calculate_surface_and_volume bool
Purpose: Determines whether the surface area and volume of the cavity are calculated.
bool is of Boolean type. Default: .true.

As the only currently implemented mpe_nonelectrostatic_model linear_OV requires the calculated measures, this flag is automatically turned on when it has been turned off but is needed.

Tag: mpe_xml_logging(control.in)

Usage: mpe_xml_logging filename level
Purpose: Controls the MPE module’s internal XML logging output.
filename specifies the name of the log file to be written. Default: mpe_interface.xml
level defines the detail of the output. Supported log levels are: off, basic, medium, and detailed. Default: off

This keyword is intended for debugging purposes. Note: Depending on the log level, the size of the output can become quite large.

3.18.2 SMPB Implicit Electrolyte Model

In FHI-aims, implicit solvation effects or electrolyte effects (z:z electrolytes) can be included by solving the Stern- and finite ion-size Modified Poisson-Boltzmann equation ((S)MPBE) in each SCF step:

\nabla\cdot\left[\varepsilon[n_{\rm el}(\bm{r})]\nabla v(\bm{r})\right]\;=\;-4% \pi n_{\rm sol}(\bm{r})-4\pi n_{\mathrm{ion}}^{\mathrm{MPB}}(\bm{r})\quad,

(3.53)

with

n_{\mathrm{ion}}^{\mathrm{MPB}}(\bm{r})\;=\;z\left[c^{\rm s}_{+}(\bm{r})-c^{% \rm s}_{-}(\bm{r})\right]\quad,

(3.54)

where $\varepsilon[n_{\mathrm{el}}]$ is a parameterized function of the electron density, $v$ is the electrostatic potential, $n_{\mathrm{sol}}$ is the solute charge density consisting of electrons and nuclei and $n_{\mathrm{ion}}^{\mathrm{MPB}}$ is the ionic charge density modeled as a function of the exclusion function $\alpha_{\mathrm{ion}}[n_{\mathrm{el}}]$ being parameterized via the electron density and the electrostatic potential $v$ . The implementation so far supports different kind of models for the ionic charge density, that is the modified, the linearized or the standard PBE. All models include a model for the Stern layer by a repulsion of the ions from the solute modeled via $\alpha_{\mathrm{ion}}[n_{\mathrm{el}}]$ and the size-modified version also a finite ion size $a$ . Parameterizations are needed for the dielectric function ( $n_{\mathrm{min}}$ and $n_{\mathrm{max}}$ ) and nonmean-field interaction of solvent with solute ( $(\alpha+\gamma)$ and $\beta$ ) which are readily available for water solvents but have to be obtained first for other solvents. Ionic parameters (ion size $a$ and Stern layer defining parameters $d_{\alpha_{\mathrm{ion}}}$ and $\xi_{\alpha_{\mathrm{ion}}}$ ) are not known so far and we are currently working on deriving them.

The energies are outputted in the end of FHI-aims under the header MPBE Solvation Additional Free Energies:

•

Total energy = Electrostatic part of the energy. This does NOT consider yet any non-electrostatic corrections (see next term)
•

Free Energy in Electrolyte = $\Omega_{\circ}$ in ref. [262]. Free energy of solute in electrolytic environment, which is Total energy + Nonelectrostatic Free Energy + Additional Nonelstatic MPBE Solvation Energy, where Total energy is the normally outputted energy in Aims (electrostatic part)
•

Surface Area of Cavity = quantum surface of solvation cavity
•

Volume of Cavity = quantum volume of solvation cavity
•

Nonelectrostatic Free Energy = non-electrostatic part of solvation energy due to solute-solvent interactions, $\Omega^{\mathrm{non-mf}}$ in the publication
•

Additional Nonelstatic MPBE Solvation Energy = non-electrostatic part of free energy due to ions. For ion-free calculations this is zero.

For more details see [262, 263]. If you want to do any calculations considering solvent or ion effects, please contact the authors, we are happy to help and cooperate.

The keywords listed here are the main part of all keywords. Some of the keywords were left out because they are highly experimental, if one is interested in more options, please contact the authors.

Tags for general subsection of control.in:

Tag: solvent mpb(control.in)

Usage: solvent mpb
Purpose: Switches MPB solvent effects on.
Restriction: Only for cluster systems (no periodic systems).

solvent mpb sub-tag: dielec_func(control.in)

Usage: dielec_func type parameters
Purpose: Define the dielectric function.
type integer describes the type of dielectric function used, type=0 Fattebert & Gygi[89] or type=1 Andreussi & Marzari[11]
parameters settings for dielectric function, separated by space:
type=0: bulk dielectric constant $\epsilon^{\mathrm{s,bulk}}$ , $\beta$ , $n_{0}$
type=1: bulk dielectric constant $\epsilon^{\mathrm{s,bulk}}$ , $n_{\mathrm{min}}$ , $n_{\mathrm{max}}$
Default: 1 78.36 0.0001 0.005[11]

solvent mpb sub-tag: ions_{parameter}(control.in)

Usage: ions_{parameter} parameter
Purpose: Set the parameters defining the ions in the electrolyte. In our recent publication[263] we explain how to choose these for different monovalent salt solutions.
parameter {parameter} =

•

temp (temperature (K))
•

conc (bulk concentration $c^{\mathrm{s,bulk}}$ (mol/L))
•

charge ( $z$ )
•

size (lenght of lattice cell a (Å))
•

kind (0 for sharp step function, 1 for smooth function)
•

mod_alpha ( $d_{\alpha_{\mathrm{ion}}}$ , $\xi_{\alpha_{\mathrm{ion}}}$ )

Defaults: T = 300K, $c^{\mathrm{s,bulk}}=$ 1M,z=1,a=5, kind = 1, $d_{\alpha_{\mathrm{ion}}}=$ 0.5, $\xi_{\alpha_{\mathrm{ion}}}=$ 1.0
Remarks: The inclusion of a second $\alpha_{\mathrm{ion}}$ function for the anions is experimental and should not be used. The use of a sharp cutoff function for $\alpha_{\mathrm{ion}}$ is not recommended, not properly implemented and just there for testing purposes.

solvent mpb sub-tag: SPE_{setting}(control.in)

Usage: SPE_{setting} parameter
Purpose: Change numerical parameters of the SPE solver.
parameter {setting} =

•

lmax (maximum angular momentum $l_{\mathrm{max}}$ of multipole expansion and of all species)
•

conv ( $\tau_{\mathrm{MERM}}$ , $\eta$ , separated by space)
•

cut_and_lmax_ff (distance from atom centers at which far field is turned on – multipole_radius_SPE, $l_{\mathrm{max}}^{\mathrm{ff}}$ – maximum angular momentum in the far field, separated by space)

Defaults: $l_{\mathrm{max}}$ = max(l_hartree), $\tau_{\mathrm{MERM}}=$ 1e-10, $\eta=$ 0.5, ${l_{\mathrm{max}}^{\mathrm{ff}}=l_{\mathrm{max}}}$ , multipole_radius_SPE is per default not used and the species dependent multipole_radius_free + 2.0 is used as far field cutoff radius.
Remarks: Due to our present tests, we do not recommend to use ${l_{\mathrm{max}}^{\mathrm{ff}}<l_{\mathrm{max}}}$ , the errors in the energies at the normal cutoff radius are too big. $\tau_{\mathrm{MERM}}=$ 1e-8 can be enough in most cases and speed up the calculation. The species dependend l_hartree can be by implementation not larger than $l_{\mathrm{max}}$ , so it is reduced to $l_{\mathrm{max}}$ if higher for the SPE solver.

solvent mpb sub-tag: dynamic_{quantity}_off(control.in)

Usage: dynamic_{quantity}_off
Purpose: If these keywords are used, {quantity} is parameterized before the SCF cycle from the superposition of free energy densities.
{quantity} =

•

cavity dielectric function $\varepsilon$
•

ions exclusion function $\alpha_{\mathrm{ion}}$

Default: both keywords not used by default, so both quantities are calculated self-consistently by parameterizing it with the full electron density.

solvent mpb sub-tag: delta_rho_in_merm(control.in)

Usage: delta_rho_in_merm
Purpose: Setting this keyword, evaluates the change of the source term ${q-\frac{1}{4\pi}\hat{L}_{1}\delta v_{n+1}}$ during the MERM iteration and solves the SPE for this change rather than the full source density.
Default: Not used. This keyword is under development and experimental, do not use it, yet.

solvent mpb sub-tag: nonsc_Gnonmf(control.in)

Usage: nonsc_Gnonmf
Purpose: Setting this keyword, calculates the free energy term $\Omega^{\mathrm{non-mf}}$ as a post-correction after the convergence of the SCF cycle, so no Kohn-Sham correction is added which would normally arise from this term. This has been proven to give very similar results for solvation energies like the fully self-consistent calculation of this term. Since people observed numerical instabilities due to this term, sometimes it might be better to set this flag.
Default: Not used. Fully self-consistent evaluation of $\Omega^{\mathrm{non-mf}}$

solvent mpb sub-tag: Gnonmf_FD_delta(control.in)

Usage: Gnonmf_FD_delta parameter
parameter $\Delta$ parameter defining the thickness of the cavity
Purpose: Used to calculated the quantum surface $S$ and volume $V$ to evaluate the free energy contribution $\Omega^{\mathrm{non-mf}}$
Default: 1e-8

solvent mpb sub-tag: not_converge_rho_mpb(control.in)

Usage: not_converge_rho_mpb
Purpose: Setting this keyword, runs a vacuum calculation first and then subsequently solves the MPBE once with the vacuum electron density and then outputs all energetics.
Default: Not used. This could be of interest for either very big systems to get first approximations without running the Newton method in each SCF step but only once, but of course then does not involve any self-consistent solution of the coupled Kohn-Sham and MPB equations. Originally, this feature was introduced to evaluate electrostatic potentials and compare them to other codes, like e.g. FEM codes.

solvent mpb sub-tag: solve_lpbe_only(control.in)

Usage: solve_lpbe_only logical
Purpose: Instead of solving the MPBE, solve the linearized version of this, also called the LPBE. For neutral molecules electrostatic fields are often small, so the LPBE electrostatic potential is often a good approximation to the true MPBE potential. The solution of the LPBE can be done directly using the MERM without the Newton method and is therefore faster for most cases.
logical if .True., use the LPB electrostatic potential, but the MPB free energy expression which contains additional entropic terms compared to the LPB expression.
Default: By default the MPBE is solved, so this is not used.

solvent mpb sub-tag: MERM_in_SPE_solver(control.in)

Usage: MERM_in_SPE_solver logical
Purpose: Do the MERM iterations inside the SPE_solver.f90 routine without updating $\delta v_{n+1}$ on the full integration grid at each step, but only on the points where we actually need it to form the source term. By this, we can gain speed, especially for $c^{\mathrm{s,bulk}}=0$ .
logical
Default: .true. Remark: In general the both options should give exactly the same result at convergence. If any difficulties arise, one is however recommended to try the .False. options, since it should be the more stable version of the solver.

solvent mpb sub-tag: MERM_atom_wise(control.in)

Usage: MERM_atom_wise logical
Purpose: Do the MERM iterations for each atom separately, i.e. we write eq. (33)[262] (Generalized Poisson or LPB-kind of equation) as:

	$\displaystyle\left(\nabla\left[\varepsilon\nabla\right]-h^{2}[v_{n}]\right)% \delta v_{n+1,\mathrm{at}}=-4\pi\varepsilon p_{\mathrm{at}}q[v_{n}]$		(3.55)
	$\displaystyle\delta v_{n+1}=\sum_{\mathrm{at}}\delta v_{n+1,\mathrm{at}}$		(3.56)
	$\displaystyle q[v_{n}]=\sum_{\mathrm{at}}p_{\mathrm{at}}q[v_{n}]$		(3.57)

In order to perform the MERM iterations for each atom, the full grid of the respective atom has to be used, i.e. also the electron density needs to be updated on points where commonly the partition_tab is vanishing. However, by this we avoid the cross-update of atomic potentials on the atomic grid of other atoms as needed in the original method and this is usually most costly in particular for larger systems. In terms of convergence with the maximum angular momentum $l_{\mathrm{max}}$ , this method performs a bit worse than the original method, which is why we recommend to use $l_{\mathrm{max}}=8$ for production runs. Still this method should be faster also with this higher accuracy in the multipole expansion.
logical
Default: .false.
Remark: Using this flag will automatically set MERM_in_SPE_solver = .True..

solvent mpb sub-tag: set_nonelstat_params(control.in)

Usage: set_nonelstat_params value value
Purpose: Set the parameters for the nonelectrostatic solvent-solute interactions.
value value two real numbers, $\alpha+\gamma$ (dyn/cm) and $\beta$ (GPa), separated by space.
Default: $\alpha+\gamma=50$ dyn/cm, $\beta=-0.35$ GPa

3.18.3 COSMO Implicit Solvent Model

In FHI-aims, we implemented a widely used conductor-like screening model (COSMO) model[167, 321] for the implicit solvation effects. The solvent are modeled as a conductor with $\epsilon=\infty$ , thus within COSMO, we solve a Poisson equation with a modified boundary condition. Given, real solvents have finite dielectric constants, a correction coefficient of this form is usually multiplied onto the calculated solute-solvent interaction

f(\epsilon)=\frac{\epsilon-1}{\epsilon+k},

(3.58)

where k is an empirical parameter usually takes 0.5. The boundary is defined by a user provided surface grid. Our implementation closely follows the smooth version of COSMO developed by York and Karplus[321]. For further details, users are encouraged to refer to their paper.

Tags for general subsection of control.in:

Tag: solvent cosmo(control.in)

Usage: solvent cosmo
Purpose: Switches COSMO solvent effects on.
Restriction: Only for cluster systems (no periodic systems).

Tag: cosmo_epsilon(control.in)

Usage: cosmo_epsilon epsilon
Purpose: Define the dielectric constant of the solvent.

Tag: cosmo_grid(control.in)

Usage: cosmo_grid file:cavity.bin
Purpose: Read the surface grid from the file cavity.bin. Please see the next page for detailed descrption of the grid file.

Description of the surface grid file.

All used references to equations rely on the paper of York and Karplus (DOI:10.1021/jp992097l). The very first line of the binary file contains a flexible version stamp "cosmo_version_yymmdd". The following lines will provide information in order to compute the A-matrix and B-vector (eq. 66, 67,72 ) unambiguously. For this, the following quantities are mandatory:
$\rightarrow$ Atomic information:
atom coordinates (already provided in input file)
atom radii $R_{i}$
atomic switching widths $R_{sw,i}$ (lhs of eq. 87)
atomic exponential terms $\zeta_{i}$ (lhs of eq. 69)
number of molecular cavity grid points $N_{i}$ associated with a particular atom
$\rightarrow$ Molecular cavity point information:
point coordinates $r_{ik}=\left(x_{ik},y_{ik},z_{ik}\right)$
contributing neighbors $n_{\text{neighbors,ik }}$ to compute associated switch value
We keep this separated scheme for the atomic and cavity information. So the first $n_{\text{atom }}$ rows contain the atomic information with each row providing the mentioned information above:

R_{i}

R_{sw,i}

\zeta_{i}

N_{i}

Note that the provided exponential terms $\zeta_{i}$ are not grid point dependent (different from eq. 61). After $n_{\text{atom }}$ rows the information of the molecular cavity points is provided as:

$r_{ik}=\left(x_{ik},y_{ik},z_{ik}\right)$	$n_{\text{neighbors,ik }}$ (e.g. 4)	$7,1,3,2$
$r_{jl}=\left(x_{jl},y_{jl},z_{jl}\right)$	$n_{\text{neighors, }jl}$ (e.g. 0)	empty

The order of those rows strictly follows the given atomic input, that is the first $N_{i}$ rows correspond to the molecular cavity points associated with the first atom specified, the following $N_{j}$ rows contain information of those points associated with the second atom and so on. Note that we provide the integer $n_{\text{neighbors }}$ after the cavity point coordinates. This value is a counter for the following integers that correspond to neighbor atoms contributing to the switch value associated to the given cavity point. The numbering of neighbor atoms follow the order of the atomic input and starts at 1 for the first atom. Since all provided cavity points represent the molecular cavity (i.e. the solvent accessible surface) the switch value $S_{ik}$ (lhs eq. 70) of a cavity point with $n_{\text{neighbor }}=0$ is equal to 1 . On the other hand, if $n_{\text{neighbor }}\neq 0$ it holds for the switch value $0<S_{ik}<1$ .

The actual switch value associated with a cavity point is a product (eq. 70) where each neighbor atom contributes with a specific factor evaluated by the switching function $S_{wf}(r)$ in eq. 64 . We follow eq. 71 for the computation of the argument of $S_{wf}$ where we apply eq. 62 for $R_{in,i}$ and make use of eq. 87 with our provided $R_{\text{sw, }i}$ .

3.18.4 Interface for Continuum Embedding Package Environ

FHI-aims features an interface for the continuum solvation package Environ.[11] The main quantities communicated from FHI-aims to Environ are the geometry and nuclear charges, the unit cell, and the electron density. In return, Environ communicates to FHI-aims the solvent contributions to the Kohn-Sham potential, the free energy, and the atomic forces. Additionally, various metadata and grid information is communicated. Communication happens at each SCF step, resulting in a self-consistent solution. The final energy output of FHI-aims contains all solvent contributions and double counting corrections for the solvent-electron interaction, and it should be interpreted as the free energy of the system in solution. The final force output are the nuclear derivatives of this free energy.

The electron density that is passed to Environ is NOT the actual one used in FHI-aims. To allow for a stable electrostatic solution on Environ’s evenly spaced grid, the core regions are smoothened in a way that conserves the electrostatic potential in the solvent region. This conservation results in the correct solvent potential, energy, and force contributions. Furthermore, Environ treats the nuclei as Gaussian charge distributions rather than point charges. All of this is typically fully automated and hidden from the user. However, when using Environ to write quantities like the charge density or the full electrostatic potential to cube files, one should be aware that these are not the actual quantities used in FHI-aims. The full procedure is described in ref. [90].

Currently, only orthorhombic cells are implemented. For surfaces with non-orthogonal primitive cells, orthogonal supercells should be used if possible (e.g., the conventional cell instead of a rhombic primitive cell). If no such supercell can be constructed, the surface is, unfortunately, outside the scope of this interface at the current time.

Input

To use both programs together, FHI-aims requires its usual control.in and geometry.in files. Environ requires an input file called environ.in. Its contents are described in the Environ documentation at https://environ.readthedocs.io. To activate the communication on the side of FHI-aims, the solvent keyword needs to be set to Environ.

Additionally, when running an isolated system, a unit cell with three lattice vectors (see lattice_vector) and a k_grid (typically 1 1 1) needs to be defined. Interactions with periodic images are removed with Environ’s pbc_correction = ’parabolic’ and pbc_dim = 0 options. The unit cell should be cubic and large enough to allow for some space between the solute and the cell boundary. However, a larger cell also increases Environ’s computational cost. If in doubt, it is advisable to test the convergence of energy and forces with cell size.

In surface slab calculations with Environ, the net dipole should NOT be compensated using use_dipole_correction in FHI-aims. Instead, pbc_correction = ’parabolic’ with pbc_dim = 2 and an appropriate pbc_axis (= 3 for surface normal in z direction) should be used in Environ. This compensates the field at the cell boundary while taking the solvent contribution to the field into account.

In contrast to plane-wave DFT, FHI-aims has no regular grid which Environ could inherit. Therefore, the cutoff parameter env_ecut which defines the grid density always needs to be given in the &ENVIRON namelist. The exact value which is required to achieve energy and force convergence differs between solutes and between solvent models. Some examples are given in ref. [90], but we encourage you to perform convergence tests on some representative examples of the kind of system that you are interested in.

When using FHI-aims together with the SCCS[11] model as implemented in Environ, some systems may require excessively dense integration grids in Environ. This is alleviated by using the deriv_lowpass_p1 and deriv_lowpass_p2 parameters. When using SCCS, it is generally recommended to set these parameters to 7 and 4, respectively. Evaluating the electron density from FHI-aims on a regular grid and taking the Fourier transform generally exhibits high frequency components which produce artifacts in the derivatives. While this is generally unproblematic in SSCS,[92] it does have an impact on the derivatives of the cavity boundary in SCCS due to its dependence on the electron density. Setting the lowpass parameters will filter out high frequency components when taking any derivatives (gradient, Laplacian, Hessian, or forces) in Fourier space.

Additionally, FHI-aims may also require the really_tight integration grids (see sec. 2.2) to produce converged results. When using SCCS, it should generally be tested if this is necessary, or if tight is sufficient.

Lastly, solver = ’fixed-point’ generally yields higher force accuracy than solver = ’cg’ (conjugate gradient). With SCCS, the latter may produce inconsistent forces for some systems.

Unless linked to the respective entry in this manual, the parameters mentioned here are defined in the input file of Environ, not FHI-aims, and are accordingly documented in Environ’s documentation. We mention them here only to make the user aware of these issues which are specific to the combination of FHI-aims and Environ and to point out possible remedies.

Running FHI-aims with Environ

Communication between FHI-aims and Environ is entirely file based. FHI-aims and Environ run as separate executables in the same working directory. FHI-aims knows to communicate with Environ when solvent Environ is set. Environ knows to communicate with FHI-aims when the -n with_aims command line flag is used. It is required that both programs use the same number of MPI processes. Additionally, it is recommended to redirect the STDERR of FHI-aims to AIMS_CRASH and that of Environ to CRASH (both in the working directory). Each program will look for error messages from the other in the respective file. If it finds any, it will stop waiting for input and terminate. This helps prevent runaway processes which wait indefinitely for input from another program which is no longer running. It is also worth mentioning that if the outputs of both programs should be gathered into a single output file, this needs to be done in append mode. It should be checked beforehand that the file does not exist, as it will otherwise be appended instead of overwritten. However, all relevant results are written by FHI-aims, not Environ, so simply redirecting FHI-aims’ output as usual is often enough. An example command to run FHI-aims with Environ could look like this:

mpiexec -np 4 $environ_binary -n with_aims 2>CRASH | tee -a aims.out & mpiexec -np 4 $aims_binary 2>AIMS_CRASH | tee -a aims.out

with $aims_binary pointing to the usual FHI-aims executable and $environ_binary pointing the Environ’s driver executable.

FHI-aims with Environ on HPC clusters

When running FHI-aims with Environ on an HPC cluster with a job scheduler, you need to make sure that both programs are allowed to use the full number of physical cores. By default, workload managers commonly block resources allocated to one MPI job from being used by any other MPI job. This prevents FHI-aims from running on the same physical cores as Environ. However, ideally, FHI-aims and Environ should actually use the same physical cores because they never perform actual work simultaneously. While one is running its calculations, the other is idle waiting for input from the former. The workload manager should be instructed to allow MPI jobs to share resources. For example, in the SLURM scheduler, this is achieved with the --overlap flag for the srun command, see https://slurm.schedmd.com/srun.html#OPT_overlap.

Dealing with crashed calculations

If a calculation ever crashes, temporary files may not be properly cleaned up and need to be removed manually before rerunning the calculation. Apart from the error message files CRASH, AIMS_CRASH, and AIMS_ENV_CRASH, the names of these temporary files typically start with aims_to_env, env_to_aims or env_pot. Additionally, it should be checked before rerunning that no runaway processes exist anymore. When redirecting STDERR as described above, processes should typically terminate, but this may take some time. If, for example, FHI-aims is manually terminated with CTRL+C while Environ is running its solver, Environ will first complete its current calculation before checking for an error message from FHI-aims. If the error file has been manually deleted in the meantime, Environ will continue to wait for input indefinitely and needs to be cancelled manually.

Optional keyword

Tag: environ_rcut(control.in)

Usage: environ_rcut $r^{\text{cut}}_{1}$ $r^{\text{cut}}_{2}$ ...
Purpose: Overrides the default cutoff radii for the core regions in which the electron density which is passed to Environ is smoothened.
$r^{\text{cut}}_{1}$ $r^{\text{cut}}_{2}$ ... is a list of real numbers containing the cutoff radius of each species in Å.
Default: Determine automatically, see ref. [90].