User Tools

Site Tools


howto:resp

This is an old revision of the document!


How to fit RESP charges with CP2K

Introduction

In CP2K, Restrained Electrostatic Potential (RESP) charges can be fitted for periodic and nonperiodic systems. It is automatically decided by the program whether a periodic or nonperiodic RESP fit is carried out. If the electrostatic (Hartree) potential is periodic (i.e. a periodic Poisson solver is used), a periodic RESP fit is performed. If the Hartree potential is computed using a nonperiodic Poisson solver, the nonperiodic fitting is employed. In any case, a least squares fitting procedure at defined grid points $\mathbf{r}_k$ is carried out,

\begin{equation} R_{\mathrm{esp}}=\frac{1}{N}\sum_k^N{(V_{\mathrm{QM}}(\mathbf{r}_k)-V_{\mathrm{RESP}}(\mathbf{r}_k))^2}, \end{equation}

where $V_{\mathrm{QM}}$ is the quantum mechanical (QM) potential and $V_{\mathrm{RESP}}$ the potential generated by the RESP charges. $N$ is the number of selected fit points.

Nonperiodic fitting

The fitted potential is obtained from a set of point charges $\{q_a\}$ center at atom $a$ according to \begin{equation} V_{\mathrm{RESP}}(\mathbf{r}_k) = \sum_a\frac{q_a}{|\mathbf{R}_a-\mathbf{r}_k|}. \end{equation} For more details see: * J. Phys. Chem., 97 , 10269-10280 (1993).

Periodic fitting

The fitted potential is generated from the charge distribution $\rho$, \begin{equation} \rho(\mathbf{r})=\sum_a{q_a g_a(\mathbf{r},\mathbf{R}_a)}, \end{equation} where $g_a$ is a Gaussian function centered at atom $a$. The periodic fitting is embedded in a Gaussian and plane waves (GPW) framework and described in detail in * Phys. Chem. Chem. Phys., 17 , 14307-14316 (2015).
In the periodic case, CP2K offers also the possibility to fit the variance of the potential instead of the absolute values, see below.

Basic input

The RESP fitting is a post-SCF step and included as a subsection of the ''PROPERTIES'' section.

&PROPERTIES
  &RESP
     &SPHERE_SAMPLING
     &END
  &END RESP
&END PROPERTIES

With this basis setup, the following defaults are employed:

  • All points outside the van der Waals radii (taken from the Cambridge database) of the atoms are included
  • The total charge of the system as defined in ''CHARGE'' is retained, i.e. ''INTEGER_TOTAL_CHARGE'' is set to .TRUE. by default.
  • All atoms except the hydrogens are weakly restrained to zero, i.e. ''RESTRAIN_HEAVIES_TO_ZERO'' is set to .TRUE. by default.

Sampling of fit points

There are different options to sample the fit points. The fit points can be sampled in shells/spheres around the atoms or, for slab-like systems, in a certain range above the surface. In any case, the systems and the sampled fitting points can be printed as .xyz file and visualized with, e.g. VMD, by enabling:

&RESP
  ....
  &PRINT
    &COORD_FIT_POINTS
    &END
  &END 
&END RESP

For better visualization it is recommended to center the coordinates of the systems using ''CENTER_COORDINATES''.

Sphere sampling

Fig. ##: Methanol molecule and fitting points (gray) sampled. This type of sampling is employed for isolated molecules and porous periodic structures suchs as metal-organic frameworks (MOFs). All grid points within a given spherical shell around the atom are included in the fitting, see figure ##. The spherical shells are defined by a minimal radius r$_{\mathrm{min}}$ and a maximal radius r$_{\mathrm{max}}$. The paramters r$_{\mathrm{min}}$ and r$_{\mathrm{max}}$ can be defined specifically for each element and are, by default, based on the van der Waals (vdW) radii. For the vdW radii, the values from the Cambridge Structural Database CAMBRIDGE or the Universal Force Field UFF can be specified via ''AUTO_VDW_RADII_TABLE''. Using the keywords AUTO_RMIN_SCALE and AUTO_RMAX_SCALE, r$_{\mathrm{min}}$ and r$_{\mathrm{max}}$ are then calculated as follows:

  • r$_{\mathrm{min}}$ = AUTO_RMIN_SCALE $\cdot$ vdW_radius
  • r$_{\mathrm{max}}$ = AUTO_RMAX_SCALE $\cdot$ vdW_radius

These settings can be overwritten for all atoms using the keywords RMIN and RMAX or only for specific atoms using RMIN_KIND and RMAX_KIND. In the following example, r$_{\mathrm{min}}$ is overwritten for all carbon atoms to 2.1$\,\mathring{\mathrm{A}}$.

&SPHERE_SAMPLING
  AUTO_VDW_RADII_TABLE CAMBRIDGE
  AUTO_RMIN_SCALE 1.0
  AUTO_RMAX_SCALE 10.0
  RMIN_KIND 2.1 C
&END

Slab sampling

RESP charges can be also fitted for slab-like systems. In this case, the potential should be well reproduced above the surface where, e.g., adsorption processes take place. The input for a flat monolayer, see figure ##, is for example:

Fig. ##: Flat monolayer and fitting points (gray).

&SLAB_SAMPLING
  RANGE 1.0 3.0
  ATOM_LIST 1..32
  SURF_DIRECTION Z
&END
&SLAB_SAMPLING
  RANGE 1.0 3.0
  ATOM_LIST 1..32
  SURF_DIRECTION -Z
&END

Fig. ##: Fit points sampled over surface atom. Here, the system is 2D-periodic in the xy-dimensions and the fitting points are sampled above (+z-direction) and below (-z-direction) which is specified by SURF_DIRECTION. With the keyword ATOM_LIST the atoms that constitute the surface are defined. This list should contain indexes of atoms of the first surface layer. RANGE defines that the points are sampled between 1-3 $\mathring{\mathrm{A}}$ above the surface layer.
The sampling technique is flexible enough to follow a corrugation of the surface. The sampling technique works as follows: An orthogonal box with box length $abc$ is constructed over each surface atom, see figure ##. All fitting points within this box are included in the fitting. The height $c$ of the box is defined by the keyword RANGE. The length $a=b$ are given by the keyword LENGTH. For flat surface layers, LENGTH is set to a sufficiently large values (3 $\mathring{\mathrm{A}}$ or more). For corrugated surface layers, LENGTH should be in the range of the distance between the surface atoms.
An input example for a corrugated graphene layer on ruthenium, see figure ##, is given below. ATOM_LIST lists in this case the indexes of the carbon atoms.

Fig. ##: Corrugated graphene layer on ruthenium and fitting points (gray)

&SLAB_SAMPLING
  RANGE 2.0 4.0
  LENGTH 2.0
  ATOM_LIST 1..1250
  SURF_DIRECTION Z
&END

Constraints

A constraint on the total charge of the system is introduced by the keyword ''INTEGER_TOTAL_CHARGE'', which is set by default to .TRUE.. Further explicit constraints can be given via the subsection ''CONSTRAINT''. It is possible to enforce the same charges for chemically equivalent atoms, e.g. for the hydrogen atoms of a methyl group. The corresponding input is:

&CONSTRAINT
  EQUAL_CHARGES
  ATOM_LIST 1 2 3
&END

where ATOM_LIST lists the indexes of the atoms that should have the same charge.
The definition of more elaborate constraints is also possible. The constraints are always linear following the formula $\sum_i^{n\_list}c_iq_i=t$. The sum is running over the atoms given in ATOM_LIST and $t$ is the target value of the constraint given by ''TARGET''. The coefficients $\{c_i\}$ are defined by the keyword ''ATOM_COEF''. With the following input it is achieved that the (absolute value) of the charge on atom 3 is twice as large as the charge on atom 5, i.e. $q_3=2q_5$.

&CONSTRAINT
  ATOM_LIST 3 5
  ATOM_COEF 1.0 -2.0
  TARGET 0.0
&END

Restraints

To avoid unphysical values for the fitted charges, restraints can be set. The restraints in CP2K are addressed by harmonic penalty functions, \begin{equation} R_{\mathrm{rest}} = \beta \sum_j (q_j-t_j)^2, \end{equation} where $t_j$ is the target value for charge $q_j$ and $\beta$ is the strength of the restraint. By default, all elements except hydrogen are weakly restrained to zero, i.e. the keyword ''RESTRAIN_HEAVIES_TO_ZERO'' is set to .TRUE. by default. The strength of this restraint is controlled by ''RESTRAIN_HEAVIES_STRENGTH''. Instead, restraints can be also defined explicitly via the subsection ''RESTRAINT'':

&RESP
  ...
  &RESTRAINT
    ATOM_LIST 1..3
    TARGET -0.18
    STRENGTH 0.0001
  &END
  &RESTRAINT
    ATOM_LIST 4
    TARGET 0.21
    STRENGTH 0.0001
  &END
  RESTRAIN_HEAVIES_TO_ZERO .FALSE.
&END RESP

In this example, charges on atoms with indexes 1..3 are restrained to -0.18 and the charge on atom 4 to 0.21. The target values $t_j$ of the restraints can be, e.g., inspired from DDAPC, Mulliken charges etc. Note the RESTRAIN_HEAVIES_TO_ZERO should be switched to .FALSE. when defining restraints explicitly. The strength $\beta$ of the restraint is defined by ''STRENGTH''. Large values for $\beta$ will limit increasingly the flexibility of the charge fitting and decrease the quality of the fit.

Fitting the variance (REPEAT method)

CP2K offers also the possibility to fit the variance of the potential as proposed in * J. Chem. Theory Comput., 5 , 2866–2878 (2009). This is only valid for periodic systems, since the reference state of the ESP is arbitrary in the periodic case. The modified residual reads:

\begin{equation} R_{\mathrm{repeat}}=\frac{1}{N}\sum_k^N{(V_{\mathrm{QM}}(\mathbf{r}_k)-V_{\mathrm{RESP}}(\mathbf{r}_k)-\delta)^2}, \end{equation} where \begin{equation} \delta = \frac{1}{N}\sum_k^N(V_{\mathrm{QM}}(\mathbf{r}_k)-V_{\mathrm{RESP}}(\mathbf{r}_k)) \end{equation}

When $V_{\mathrm{QM}}$ is obtained from, e.g., a plane wave code and the periodicity of $V_{\mathrm{RESP}}$ is later treated by, e.g., Ewald summation, both potentials will have different offsets. The modified residual $R_{\mathrm{repeat}}$ was introduced to overcome this problem. In CP2k, $V_{\mathrm{QM}}$ and $V_{\mathrm{RESP}}$ are both evaluated with the same method, the Gaussian and plane waves approach, and have thus the same offset. However, fitting the variance of the potential is an easier task than fitting the absolute values and avoids a strong fluctuation of the charges. Fitting the variance of the potential will thus also stabilize the fitting procedure in CP2K when $\delta\ne0$. The value of $\delta$ depends on the sampling of fitting points. If all points are included, it holds that $\delta=0$ since $\sum_k^{N_{all}}V(r_k)=0$

Check the quality of the fit

Example input files

howto/resp.1447843299.txt.gz · Last modified: 2020/08/21 10:15 (external edit)