CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
howto:resp

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
howto:resp [2015/11/11 11:32] – [Restraints] dgolzehowto:resp [2024/01/15 09:24] (current) – oschuett
Line 1: Line 1:
-====== How to fit RESP charges with CP2K ====== +This page has been moved to: https://manual.cp2k.org/trunk/methods/properties/resp_charges.html
- +
-===== 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: +
-* [[  http://dx.doi.org/10.1021/j100142a004 | 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  +
-* [[  http://dx.doi.org/10.1039/C4CP04638B | 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 [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES.html|''PROPERTIES'']] section. +
- +
-<code cp2k> +
-&PROPERTIES +
-  &RESP +
-     &SPHERE_SAMPLING +
-     &END +
-  &END RESP +
-&END PROPERTIES +
-</code> +
- +
-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 [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT.html#CHARGE|''CHARGE'']] is   retained, i.e. [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP.html#INTEGER_TOTAL_CHARGE|''INTEGER_TOTAL_CHARGE'']] is set to ''.TRUE.'' by default. +
-  * All atoms except the hydrogens are weakly restrained to zero, i.e. [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP.html#RESTRAIN_HEAVIES_TO_ZERO|''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: +
-<code cp2k> +
-&RESP +
-  .... +
-  &PRINT +
-    &COORD_FIT_POINTS +
-    &END +
-  &END  +
-&END RESP +
-</code> +
-For better visualization it is recommended to center the coordinates of the systems using [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/TOPOLOGY/CENTER_COORDINATES.html|''CENTER_COORDINATES'']]. +
-==== Sphere sampling ==== +
-<imgcaption fig:sphere_samp |Methanol molecule and fitting points (gray) sampled.>{{ ch3oh_fitpoints.png?200x200}}</imgcaption> +
-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 <imgref fig:sphere_samp>. 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 [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP/SPHERE_SAMPLING.html#AUTO_VDW_RADII_TABLE|''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}}$. +
-<code cp2k> +
-&SPHERE_SAMPLING +
-  AUTO_VDW_RADII_TABLE CAMBRIDGE +
-  AUTO_RMIN_SCALE 1.0 +
-  AUTO_RMAX_SCALE 10.0 +
-  RMIN_KIND 2.1 C +
-&END +
-</code> +
- +
-==== 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 <imgref fig:monolayer_samp>, is for example: +
- +
-<imgcaption fig:monolayer_samp |Flat monolayer and fitting points (gray).>{{ fitpoints_monolayer.png?280}}</imgcaption> +
-<code cp2k> +
-&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 +
-</code> +
-<imgcaption fig:1atom_samp | Fit points sampled over surface atom.>{{fitpoints_1atom.png?90 }}</imgcaption> +
-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 <imgref fig:1atom_samp>. 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  <imgref fig:corrugated_samp>, is given below. ''ATOM_LIST'' lists in this case the indexes of the carbon atoms. +
- +
- +
-<imgcaption fig:corrugated_samp |Corrugated graphene layer on ruthenium and fitting points (gray)>{{ fittingpoints_corrugated.png?450}}</imgcaption> +
- +
-<code cp2k> +
-&SLAB_SAMPLING +
-  RANGE 2.0 4.0 +
-  LENGTH 2.0 +
-  ATOM_LIST 1..1250 +
-  SURF_DIRECTION Z +
-&END +
-</code> +
- +
-===== Constraints ===== +
-A constraint on the total charge of the system is introduced by the keyword [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP.html#INTEGER_TOTAL_CHARGE|''INTEGER_TOTAL_CHARGE'']], which is set by default to  ''.TRUE.''. Further explicit constraints can be given via the subsection [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP/CONSTRAINT.html|''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:  +
-<code cp2k> +
-&CONSTRAINT +
-  EQUAL_CHARGES +
-  ATOM_LIST 1 2 3 +
-&END +
-</code> +
-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 [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP/CONSTRAINT.html#TARGET|''TARGET'']]. The coefficients $\{c_i\}$ are defined by the keyword [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP/CONSTRAINT.html#ATOM_COEF|''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$. +
-<code cp2k> +
-&CONSTRAINT +
-  ATOM_LIST 3 5 +
-  ATOM_COEF 1.0 -2.0 +
-  TARGET 0.0 +
-&END +
-</code> +
- +
-===== 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 [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP.html#RESTRAIN_HEAVIES_TO_ZERO|''RESTRAIN_HEAVIES_TO_ZERO'']] is set to ''.TRUE.'' by default. The strength of this restraint is controlled by [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP.html#list_RESTRAIN_HEAVIES_STRENGTH|''RESTRAIN_HEAVIES_STRENGTH'']]. Instead, restraints can be also defined explicitly via the subsection [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP/RESTRAINT.html| ''RESTRAINT'']]: +
-<code cp2k> +
-&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 +
-</code>  +
-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 [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/PROPERTIES/RESP/RESTRAINT.html#STRENGTH|''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) ===== +
-===== Check the quality of the fit ===== +
-===== Example input files =====+
howto/resp.txt · Last modified: 2024/01/15 09:24 by oschuett