howto:resp
Differences
This shows you the differences between two versions of the page.
Both sides previous revisionPrevious revisionNext revision | Previous revision | ||
howto:resp [2015/12/02 16:53] – [Example input files] dgolze | howto: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:// |
- | + | ||
- | ===== 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$ | + | |
- | + | ||
- | \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\}$ centered at atom $a$ according to | + | |
- | \begin{equation} | + | |
- | | + | |
- | \end{equation} | + | |
- | For more details see: | + | |
- | * [[ http:// | + | |
- | + | ||
- | ==== 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}, | + | |
- | \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:// | + | |
- | 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:// | + | |
- | + | ||
- | <code cp2k> | + | |
- | & | + | |
- | &RESP | + | |
- | & | + | |
- | & | + | |
- | &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 [[http:// | + | |
- | * All atoms except the hydrogens are weakly restrained to zero, i.e. [[http:// | + | |
- | + | ||
- | ===== Sampling of fit points ===== | + | |
- | There are different options to sample the fit points. The fit points can be sampled in shells/ | + | |
- | <code cp2k> | + | |
- | &RESP | + | |
- | .... | + | |
- | & | + | |
- | & | + | |
- | &END | + | |
- | &END | + | |
- | &END RESP | + | |
- | </ | + | |
- | For better visualization it is recommended to center the coordinates of the systems using [[http:// | + | |
- | ==== Sphere sampling ==== | + | |
- | < | + | |
- | 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, | + | |
- | For the vdW radii, the values from the Cambridge Structural Database '' | + | |
- | * 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 '' | + | |
- | <code cp2k> | + | |
- | & | + | |
- | 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 <imgref fig: | + | |
- | + | ||
- | < | + | |
- | <code cp2k> | + | |
- | & | + | |
- | RANGE 1.0 3.0 | + | |
- | ATOM_LIST 1..32 | + | |
- | SURF_DIRECTION Z | + | |
- | &END | + | |
- | & | + | |
- | RANGE 1.0 3.0 | + | |
- | ATOM_LIST 1..32 | + | |
- | SURF_DIRECTION -Z | + | |
- | &END | + | |
- | </ | + | |
- | < | + | |
- | 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 '' | + | |
- | 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: | + | |
- | An input example for a corrugated graphene layer on ruthenium, see <imgref fig: | + | |
- | + | ||
- | + | ||
- | < | + | |
- | + | ||
- | <code cp2k> | + | |
- | & | + | |
- | 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 [[http:// | + | |
- | <code cp2k> | + | |
- | & | + | |
- | EQUAL_CHARGES | + | |
- | ATOM_LIST 1 2 3 | + | |
- | &END | + | |
- | </ | + | |
- | where '' | + | |
- | 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 '' | + | |
- | <code cp2k> | + | |
- | & | + | |
- | 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} | + | |
- | | + | |
- | \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:// | + | |
- | <code cp2k> | + | |
- | &RESP | + | |
- | ... | + | |
- | & | + | |
- | ATOM_LIST 1..3 | + | |
- | TARGET -0.18 | + | |
- | STRENGTH 0.0001 | + | |
- | &END | + | |
- | & | + | |
- | 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 '' | + | |
- | ===== Fitting the variance (REPEAT method) ===== | + | |
- | CP2K offers also the possibility to fit the variance of the potential as proposed in * [[ http:// | + | |
- | + | ||
- | \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} | + | |
- | | + | |
- | \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. | + | |
- | To enable this option, add the keyword '' | + | |
- | + | ||
- | <code cp2k> | + | |
- | &RESP | + | |
- | ... | + | |
- | USE_REPEAT_METHOD | + | |
- | &END RESP | + | |
- | </ | + | |
- | Note that '' | + | |
- | To obtain REPEAT charges in a stricter sense, i.e. as computed by the [[https:// | + | |
- | <code cp2k> | + | |
- | &RESP | + | |
- | USE_REPEAT_METHOD | + | |
- | & | + | |
- | | + | |
- | & | + | |
- | &END RESP | + | |
- | </ | + | |
- | Use the keyword '' | + | |
- | + | ||
- | + | ||
- | ===== Check the quality of the fit ===== | + | |
- | + | ||
- | 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}}. | + | |
- | | + | |
- | + | ||
- | Both errors are printed to the output file. They should be as small as possible. Typical values can be found here:\\ | + | |
- | * [[ http:// | + | |
- | * [[ http:// | + | |
- | * [[ http:// | + | |
- | + | ||
- | 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: | + | |
- | <code cp2k> | + | |
- | &RESP | + | |
- | .... | + | |
- | & | + | |
- | & | + | |
- | &END | + | |
- | &END | + | |
- | &END RESP | + | |
- | </ | + | |
- | The QM potential can be as well printed as cube file using | + | |
- | [[https:// | + | |
- | ===== Example input files ===== | + | |
- | + | ||
- | * Isolated methanol molecule - nonperiodic fit: {{: | + | |
- | * Metal organic framework IRMOF-1 - periodic fit using REPEAT : {{: | + |
howto/resp.1449075233.txt.gz · Last modified: 2020/08/21 10:15 (external edit)