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.

The fitted potential is obtained from a set of point charges $\{q_a\}$ centered 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).

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.

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.

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.

Fig. 1: 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 1. 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

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 2, is for example:

Fig. 2: 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. 3: 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 3. 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 4, is given below. ATOM_LIST lists in this case the indexes of the carbon atoms.

Fig. 4: 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

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

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. 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. 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. If only the explicitly given restraints should be used, RESTRAIN_HEAVIES_TO_ZERO must be switched to .FALSE..

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 GPW 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. A stabilization of the fitting procedure is thus also expected in CP2K when $\delta\ne0$. The value of $\delta$ depends on the sampling of fitting points. If all points are included, it strictly holds that $\delta=0$, since we have $\sum_k^{N_{all}}V(r_k)=0$ (in CP2K). In this case, the original and modified residuals are identical, i.e. $R_{\mathrm{esp}}= R_{\mathrm{repeat}}$. If the fitting points are sampled in spheres around the atom, which is done for molecular periodic structures like MOFs, $\delta$ will be non-zero and fitting the variance is strongly recommended. When sampling the fitting points above a surface, we often find that $\delta\sim 0$. However, minimizing $R_{\mathrm{repeat}}$ will partly also yield improvements for such systems.
To enable this option, add the keyword USE_REPEAT_METHOD:

&RESP
  ...
  USE_REPEAT_METHOD
&END RESP

Note that RESTRAIN_HEAVIES_TO_ZERO is then automatically switched to .FALSE.. Furthermore, the definition of explicit restraints is usually not necessary.
To obtain REPEAT charges in a stricter sense, i.e. as computed by the ''REPEAT code'', sphere sampling has to be enabled, the van der Waals radii must be retrieved from the Universal Force Field and the total charge must be retained. The corresponding input is

&RESP
  USE_REPEAT_METHOD
  &SPHERE_SAMPLING
     AUTO_VDW_RADII_TABLE UFF
   &END
&END RESP

Use the keyword AUTO_RMIN_SCALE and AUTO_RMAX_SCALE to scale the van der Waals radii as described above. Note that small numerical deviations compared to the REPEAT code are possible since the fitting is embedded in a GPW framwork as described in Phys. Chem. Chem. Phys., 17 , 14307-14316 (2015) , whereas the REPEAT code uses Ewald summation.

A measure for the quality of the fit are the root-mean square (RMS) error

\begin{equation} \mathrm{RMS}=\sqrt{\frac{\sum_k^N~(V_{\mathrm{QM}}(\mathbf{r}_k)-V_{\mathrm{RESP}}(\mathbf{r}_k))^2}{N}} \end{equation}

and the relative root-mean square (RRMS) error \begin{equation} \mathrm{RRMS}=\sqrt{\frac{\sum_k^N~(V_{\mathrm{QM}}(\mathbf{r}_k)-V_{\mathrm{RESP}}(\mathbf{r}_k))^2}{\sum_k^N~V_{{\mathrm{QM}}}(\mathbf{r}_k)^2}}. \end{equation}

Both errors are printed to the output file. They should be as small as possible. Typical values can be found here:

• J. Phys. Chem., 97 , 10269-10280 (1993).
  
• J. Phys. Chem. A, 114 , 10225-10233 (2010).
  
• Phys. Chem. Chem. Phys., 17 , 14307-14316 (2015).
  

When the variance is fitted, $V_{\mathrm{RESP}}$ is shifted by $\delta$ with respect to $V_{\mathrm{QM}}$. Thus, $V_{\mathrm{RESP}}$ is replaced by $\tilde{V}_{\mathrm{RESP}}=V_{\mathrm{RESP}}+\delta$ in the formulas for the RMS and RRMS values.
The RESP potential can be printed in cube file format using the following option:

&RESP
  ....
  &PRINT
    &V_RESP_CUBE
    &END
  &END 
&END RESP

The QM potential can be as well printed as cube file using V_HARTREE_CUBE. The cube files can be visualized with, e.g. VMD, and the RESP and the QM potential can be directly compared. Note that $\tilde{V}_{\mathrm{RESP}}$ is printed instead of $V_{\mathrm{RESP}}$ when the variance is fitted.

• Isolated methanol molecule - nonperiodic fit: methanol_nonperiodic.tar.gz
  
• Metal-organic framework IRMOF-1 - periodic fit using REPEAT : irmof-1_REPEAT.tar.gz
  
• Metal-organic framework MIL-53-Al - periodic fit using REPEAT: mil-53-al-repeat.tar.gz
  
• Graphene on Ru(0001) - periodic fit: graphene_Ru.tar.gz