diff --git a/docs/methods/properties/index.md b/docs/methods/properties/index.md index 5261482deb..634af480c3 100644 --- a/docs/methods/properties/index.md +++ b/docs/methods/properties/index.md @@ -13,5 +13,5 @@ bandstructure_gw nmr STM Images Electron Transport -RESP Charges +resp_charges ``` diff --git a/docs/methods/properties/resp_1atom_fitpoints.png b/docs/methods/properties/resp_1atom_fitpoints.png new file mode 100644 index 0000000000..213fd81b5b Binary files /dev/null and b/docs/methods/properties/resp_1atom_fitpoints.png differ diff --git a/docs/methods/properties/resp_ch3oh_fitpoints.png b/docs/methods/properties/resp_ch3oh_fitpoints.png new file mode 100644 index 0000000000..dbd1fd4ba6 Binary files /dev/null and b/docs/methods/properties/resp_ch3oh_fitpoints.png differ diff --git a/docs/methods/properties/resp_charges.md b/docs/methods/properties/resp_charges.md new file mode 100644 index 0000000000..30191500be --- /dev/null +++ b/docs/methods/properties/resp_charges.md @@ -0,0 +1,387 @@ +# RESP Charges + +## 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, + +$$ +R_\mathrm{esp}=\frac{1}{N}\sum_k^N{(V_\mathrm{QM}(\mathbf{r}_k)-V_\mathrm{RESP}(\mathbf{r}_k))^2}, +$$ + +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\}$ centered at atom $a$ +according to + +$$ +V_\mathrm{RESP}(\mathbf{r}_k) = \sum_a\frac{q_a}{|\mathbf{R}_a-\mathbf{r}_k|}. +$$ + +For more details see: +[J. Phys. Chem., 97 , 10269-10280 (1993).](https://dx.doi.org/10.1021/j100142a004) + +### Periodic fitting + +The fitted potential is generated from the charge distribution $\rho$, + +$$ +\rho(\mathbf{r})=\sum_a{q_a g_a(\mathbf{r},\mathbf{R}_a)}, +$$ + +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).](https://dx.doi.org/10.1039/C4CP04638B) In the +periodic case, CP2K offers also the possibility to fit the variance of the potential instead of the +absolute values, see [below](#repeat_method). + +## Basic input + +The RESP fitting is a post-SCF step and included as a subsection of the +[PROPERTIES](#CP2K_INPUT.FORCE_EVAL.PROPERTIES) section. + +```none +&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](#CP2K_INPUT.FORCE_EVAL.DFT.CHARGE) is + retained, i.e. [INTEGER_TOTAL_CHARGE](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.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: + +```none +&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](#CP2K_INPUT.FORCE_EVAL.SUBSYS.TOPOLOGY.CENTER_COORDINATES). + +### Sphere sampling + +![](resp_ch3oh_fitpoints.png){align=center width=350px} **Figure 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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SPHERE_SAMPLING.AUTO_VDW_RADII_TABLE). +Using the keywords +[AUTO_RMIN_SCALE](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SPHERE_SAMPLING.AUTO_RMIN_SCALE) and +[AUTO_RMAX_SCALE](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SPHERE_SAMPLING.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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SPHERE_SAMPLING.RMIN) and +[RMAX](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SPHERE_SAMPLING.RMAX) or only for specific atoms using +[RMIN_KIND](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SPHERE_SAMPLING.RMIN_KIND) and +[RMAX_KIND](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SPHERE_SAMPLING.RMAX_KIND). In the following +example, $r_\mathrm{min}$ is overwritten for all carbon atoms to 2.1$\,\mathring{\mathrm{A}}$. + +```none +&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 2, is for example: + +![](resp_monolayer_fitpoints.png){align=center width=450px} **Figure 2.** Flat monolayer and fitting +points (gray). + +```none +&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 +``` + +![](resp_1atom_fitpoints.png){align=center width=150px} **Figure 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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SLAB_SAMPLING.SURF_DIRECTION). With the +keyword [ATOM_LIST](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SLAB_SAMPLING.ATOM_LIST) the atoms that +constitute the surface are defined. This list should contain indexes of atoms of the first surface +layer. [RANGE](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SLAB_SAMPLING.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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SLAB_SAMPLING.RANGE) . The length $a=b$ are given by +the keyword [LENGTH](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SLAB_SAMPLING.LENGTH). For flat surface +layers, [LENGTH](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SLAB_SAMPLING.LENGTH) is set to a +sufficiently large values (3 $\mathring{\mathrm{A}}$ or more). For corrugated surface layers, +[LENGTH](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SLAB_SAMPLING.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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SLAB_SAMPLING.ATOM_LIST) lists in this case the +indexes of the carbon atoms. + +![](resp_corrugated_fitpoints.png){align=center width=550px} **Figure 4.** Corrugated graphene layer +on ruthenium and fitting points (gray). + +```none +&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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.INTEGER_TOTAL_CHARGE), which is set by +default to `.TRUE.`. Further explicit constraints can be given via the subsection +[CONSTRAINT](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.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: + +```none +&CONSTRAINT + EQUAL_CHARGES + ATOM_LIST 1 2 3 +&END +``` + +where [ATOM_LIST](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.CONSTRAINT.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_\text{list}}c_i q_i=t$. The sum is running over the atoms given in +[ATOM_LIST](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.CONSTRAINT.ATOM_LIST) and $t$ is the target value +of the constraint given by [TARGET](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.CONSTRAINT.TARGET). The +coefficients $\{c_i\}$ are defined by the keyword +[ATOM_COEF](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.CONSTRAINT.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$. + +```none +&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, + +$$ +R_\mathrm{rest} = \beta \sum_j (q_j-t_j)^2, +$$ + +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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.RESTRAIN_HEAVIES_TO_ZERO) is set +to `.TRUE.` by default. The strength of this restraint is controlled by +[RESTRAIN_HEAVIES_STRENGTH](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.RESTRAIN_HEAVIES_STRENGTH). +Restraints can be also defined explicitly via the subsection +[RESTRAINT](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.RESTRAINT): + +```none +&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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.RESTRAINT.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](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.RESTRAIN_HEAVIES_TO_ZERO) must be +switched to `.FALSE.`. + +(repeat_method)= + +## Fitting the variance (REPEAT method) + +CP2K offers also the possibility to fit the variance of the potential as proposed in +[](#Campana2009). This is only valid for periodic systems, since the reference state of the ESP is +arbitrary in the periodic case. The modified residual reads: + +$$ +R_\mathrm{repeat}=\frac{1}{N}\sum_k^N{(V_\mathrm{QM}(\mathbf{r}_k)-V_\mathrm{RESP}(\mathbf{r}_k)-\delta)^2}, +$$ + +where + +$$ +\delta = \frac{1}{N}\sum_k^N(V_\mathrm{QM}(\mathbf{r}_k)-V_\mathrm{RESP}(\mathbf{r}_k)). +$$ + +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 \ne 0$. 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\sim0$. However, minimizing $R_\mathrm{repeat}$ will partly +also yield improvements for such systems. + +To enable this option, add the keyword +[USE_REPEAT_METHOD](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.USE_REPEAT_METHOD): + +```none +&RESP + ... + USE_REPEAT_METHOD +&END RESP +``` + +Note that +[RESTRAIN_HEAVIES_TO_ZERO](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.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](https://bitbucket.org/ccampana/repeat_atomic_charges)*, 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 + +```none +&RESP + USE_REPEAT_METHOD + &SPHERE_SAMPLING + AUTO_VDW_RADII_TABLE UFF + &END +&END RESP +``` + +Use the keyword +[AUTO_RMIN_SCALE](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SPHERE_SAMPLING.AUTO_RMIN_SCALE) and +[AUTO_RMAX_SCALE](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.RESP.SPHERE_SAMPLING.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)](https://dx.doi.org/10.1039/C4CP04638B) , whereas +the *REPEAT code* uses Ewald summation. + +## Check the quality of the fit + +A measure for the quality of the fit are the root-mean square (RMS) error + +$$ +\mathrm{RMS}=\sqrt{\frac{\sum_k^N~(V_{\mathrm{QM}}(\mathbf{r}_k)-V_\mathrm{RESP}(\mathbf{r}_k))^2}{N}} +$$ + +and the relative root-mean square (RRMS) error + +$$ +\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}}. +$$ + +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).](https://dx.doi.org/10.1021/j100142a004) +- [J. Phys. Chem. A, 114 , 10225-10233 (2010).](https://dx.doi.org/10.1021/jp103944q) +- [Phys. Chem. Chem. Phys., 17 , 14307-14316 (2015).](https://dx.doi.org/10.1039/C4CP04638B) + +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: + +```none +&RESP + .... + &PRINT + &V_RESP_CUBE + &END + &END +&END RESP +``` + +The QM potential can be as well printed as cube file using +[V_HARTREE_CUBE](#CP2K_INPUT.FORCE_EVAL.DFT.PRINT.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. + +## Example input files + +- Isolated methanol molecule - nonperiodic fit: + [methanol_nonperiodic.tar.gz](https://www.cp2k.org/_media/howto:methanol_nonperiodic.tar.gz) +- Metal-organic framework IRMOF-1 - periodic fit using REPEAT: + [irmof-1_REPEAT.tar.gz](https://www.cp2k.org/_media/howto:irmof-1_repeat.tar.gz) +- Metal-organic framework MIL-53-Al - periodic fit using REPEAT: + [mil-53-al-repeat.tar.gz](https://www.cp2k.org/_media/howto:mil-53-al-repeat.tar.gz) +- Graphene on Ru(0001) - periodic fit: + [graphene_Ru.tar.gz](https://www.cp2k.org/_media/howto:graphene_ru.tar.gz) diff --git a/docs/methods/properties/resp_corrugated_fitpoints.png b/docs/methods/properties/resp_corrugated_fitpoints.png new file mode 100644 index 0000000000..a8822b1702 Binary files /dev/null and b/docs/methods/properties/resp_corrugated_fitpoints.png differ diff --git a/docs/methods/properties/resp_monolayer_fitpoints.png b/docs/methods/properties/resp_monolayer_fitpoints.png new file mode 100644 index 0000000000..f02f522b57 Binary files /dev/null and b/docs/methods/properties/resp_monolayer_fitpoints.png differ