Separation of variables in molecular-dynamics simulations: A criterion to estimate the quality of the approximation

We propose a simple method to evaluate the approximation of separation of variables in molecular dynamics (MD) and related fields. It is based on a point-by-point evaluation of the difference between the true potential and the corresponding potential where the separation of variables is applied. The major advantage of such an approach is the fact that it requires only the analytical form of the potential as provided in most of the MD codes. We provide two examples of application, namely, a diatomic molecule adsorbing on a flat surface and an alkane (aliphatic) chain.

Figure 1

Schematic representation of the system considered. In (I) the mapping from $(z_a,z_b)$ to $(r,\theta )$ is shown pictorially. If $\theta$ and $r$ are separable one can represent the interaction of the diatomic molecule with the surface as an effective center-of-mass–surface interaction and disregard the evolution with respect to $\theta$, as shown in II.

Figure 2

(Color online) Quality factor $Q$ studied as a function of the angle $\theta$ with $r$ varying parametrically from $3.0\text{to}1.0\mathrm{nm}$. The fixed point used in this example is ${\theta }_0=\pi ∕6$ and $r_0=2.0\mathrm{nm}$. Here $\sigma =0.5\mathrm{nm}$, $l=0.5\mathrm{nm}$, and $ϵ=10kT$. The negative value of $Q$ indicates that the potential term where the separation is applied overestimates the true potential. As the distance of the center of mass from the surface decreases, the dependence on $\theta$ becomes stronger, and approximately for $r⩽1.2\mathrm{nm}$ the ASV starts to be questionable since the error induced can be larger than $2kT$.

Figure 3

(Color online) Variation of $Q$ with respect to $r_0$ $\left({\delta }_{r_0}\right)$ for some example values of ${\theta }_0$ and $\theta$. The symmetric curves represent the two extreme cases of values for ${\theta }_0$ and $\theta$ while in between are reported two more examples. The message of this plot is that there is a general trend according to which the dependence of $Q$ on $x_0$ becomes crucial for, approximately, $r_0⩽1.4\mathrm{nm}$. The vertical dashed line indicates that the region $r_0⩽1.4\mathrm{nm}$ should be considered as the critical region for the variable $r$.

Figure 4

(Color online) Variation of $Q$ with respect to ${\theta }_0$ $\left({\delta }_{{\theta }_0}\right)$ for some example values of $r_0$ and $r$. For values of $r_0$ and $r$ outside the critical region defined by Fig. 3, the dependence of $Q$ on ${\theta }_0$ is negligible; however, for values within the critical region of $r$, $r_0⩽1.4$, the dependence becomes stronger and identifies the critical region of the variable $\theta$. The vertical dashed line indicates that the region below ${\theta }_0⩾\pi ∕4$ is the critical one.

Figure 5

Chemical structure of the alkane chain in which each vertex of the backbone represent a carbon atom.

Figure 6

(Color online) CG structure of the backbone for an AC and the underlying atomistic (carbon group) structure. Each carbon group (e.g., $\mathrm{C}{\mathrm{H}}_2$ or $\mathrm{C}{\mathrm{H}}_3$) is represented by a large dark (blue) circle close to a letter. (a) shows the MS where there are two carbon groups per bead (1:2), and (b) three per bead (1:3). The “atomistic” variables ${\theta }_i$, $\forall i=1,2,3,4,5$, and the CG variables $R^{\left(1\right)}$, $R^{\left(2\right)}$, and $\Omega$ are also shown.

Figure 7

(Color online) (a) Quality factor $Q_1=Q\left({\overline{R}}^{\left(1\right)}\right)={\Delta }_{R_0^{\left(2\right)},{\Omega }_0}\left({\overline{R}}^{\left(1\right)}\right)∕kT$ for the mapping 1:2 with ${\overline{R}}^{\left(1\right)}=2.4\mathrm{\AA }$ as a parametric value and the fixed points equal to $R_0^{\left(2\right)}=2.4\mathrm{\AA }$ and ${\Omega }_0=140°$. In order to determine (a) we first chose an arbitrary fixed point $(R_0^{\left(2\right)},{\Omega }_0)$ and then analyzed $Q$ and its dependence on the fixed point as shown in (b). Next, according to the procedure shown for the example of a diatomic molecule, we can identify two critical regions which are given $\forall$ $R_0^{\left(2\right)}∊[2.425,2.525]$ and $\forall$ ${\Omega }_0^{\left(2\right)}∊[145°,180°]$. Accordingly, we have then chosen the fixed point of (a) outside the critical region. In (a), the light gray region shows the vicinity of the fixed point $(R_0^{\left(2\right)},{\Omega }_0)$ where the CG potential is decoupled (by definition) and thus of no interest in this context. The real interest is in the critical regions determined by the plot shown in (b) [${\delta }_{{\Omega }_0}$, not reported, shows a behavior similar to that of (b)]. Here the potential is no longer decoupled and the variables are highly correlated. Note also that the white regions are those where the CG variables are not defined and thus $Q$ is not calculated. The same applies for all the following figures.

Figure 8

(Color online) (a) Quality $Q_2=Q\left({\overline{R}}^{\left(2\right)}\right)={\Delta }_{R_0^{\left(1\right)},{\Omega }_0}\left({\overline{R}}^{\left(2\right)}\right)∕kT$ for the mapping 1:2 with ${\overline{R}}^{\left(2\right)}=2.35\mathrm{\AA }$ as a parametric value and the fixed points equal to $R_0^{\left(1\right)}=2.5\mathrm{\AA }$ and ${\Omega }_0=162°$. The procedure used is the same as in the previous figure. Note that the value of the derivative shown in (b) is almost constant for each parametric curve (the same behavior holds for ${\delta }_{{\Omega }_0}$, not reported). This means that there exists no real definition of the noncritical region and in this case one can proceed by using several fixed points in the domain and averaging the quality factor resulting from each study.

Figure 9

(Color online) (a) Quality factor $Q_3=Q\left(\overline{\Omega }\right)={\Delta }_{R_0^{\left(1\right)},R_0^{\left(2\right)}}\left(\overline{\Omega }\right)∕kT$ for the mapping 1:2 with $\overline{\Omega }=162°$ as a parametric value and the fixed points equal to $R_0^{\left(1\right)}=2.6\mathrm{\AA }$ and $R_0^{\left(2\right)}=2.4\mathrm{\AA }$. The procedure for analyzing (b) [${\delta }_{R_0^{\left(1\right)}}$, not reported, displays a similar behavior to that reported in (b)] is the same as described in the previous figures.

Figure 10

(Color online) (a) Quality factor $Q_1=Q\left({\overline{R}}^{\left(1\right)}\right)={\Delta }_{R_0^{\left(2\right)},{\Omega }_0}\left({\overline{R}}^{\left(1\right)}\right)∕kT$ of the mapping 1:3 with ${\overline{R}}^{\left(1\right)}=2.7\mathrm{\AA }$ and ${\theta }_2=126°$ as parametric values. The corresponding fixed points are $R_0^{\left(2\right)}=3.25\mathrm{\AA }$ and ${\Omega }_0=176.5°$. As for Fig. 8, but in this case due to the rapidly varying dependence on the fixed point [see (b)] one has to use several fixed points and consider the $Q$ resulting from the average of each study, in order to have more valid quantitative information. ${\delta }_{R_0^{\left(1\right)}}$, not reported, shows a similar behavior.

Figure 11

(Color online) (a) Quality factor $Q_2=Q\left({\overline{R}}^{\left(2\right)}\right)={\Delta }_{R_0^{\left(1\right)},{\Omega }_0}\left({\overline{R}}^{\left(2\right)}\right)∕kT$ for the mapping 1:3 with ${\overline{R}}^{\left(2\right)}=3.25\mathrm{\AA }$ and ${\theta }_2=126°$ as parametric values. The corresponding fixed points are $R_0^{\left(1\right)}=2.7\mathrm{\AA }$ and ${\Omega }_0=176.5°$. As in the previous figure, there is no noncritical region [see (b)]. ${\delta }_{{\Omega }_0}$, not reported, shows the same behavior as in (b). This means that one must explore critical points over the whole domain in order to estimate $Q$ quantitatively.

Figure 12

(Color online) (a) Quality factor $Q_3=Q\left(\overline{\Omega }\right)={\Delta }_{R_0^{\left(1\right)},R_0^{\left(2\right)}}\left(\overline{\Omega }\right)∕kT$ for the mapping 1:3 with $\overline{\Omega }=174°$ and ${\theta }_2=126°$ as parametric values. The corresponding fixed points are $R_0^{\left(1\right)}=2.71\mathrm{\AA }$ and $R_0^{\left(2\right)}=3.4\mathrm{\AA }$. As in the other cases for $Q_1$, $Q_2$, and now for $Q_3$ noncritical regions within the domains of the CG variables cannot be found [see (b)] and thus the same considerations apply as in Figs. 10, 11. ${\delta }_{R_0^{\left(1\right)}}$, not reported, shows a similar behavior to that reported in (b).