Equivalence of Atomic and Molecular Pressure

The volume scaling defined in Eq. (3.1) is not unique. Note that only the equation of motion for the center of mass momentum, Eq. (3.25), has a velocity dependent term that depends on the coordinates of the barostat through the matrix ${\bf G}$ defined in Eq. (3.11). The atomic momenta, Eq. (3.24), on the contrary, are not coupled to the barostat. This fact is also reflected in the equations of motion for the barostat momenta, Eq. (3.26), which is driven by the internal pressure due only to the molecular or group center of masses. In defining the extended Lagrangian one could as well have defined an atomic scaling of the form

$$r_{ik\alpha} = \sum_{\beta} h_{\alpha \beta}s_{i\alpha k}.$$ (3.33)

Atomic scaling might be trivially implemented by eliminating the kinetic energy, which depends on the $\dot l_{ik\alpha}$ velocities, from the starting Lagrangian (3.6) and replacing the term ${1 \over 2} \sum_{i}^{N} M_{i} s^{2} {\bf\dot S}_{i}^{t}{\bf h^{t} h} {\bf\dot S}_{i}$ with ${1 \over 2} \sum_{ik} m_{ik} s^{2} {\bf\dot s}_{ik}^{t}{\bf h^{t} h} {\bf\dot s}_{ik}$. The corresponding equations of motions for atomic scaling are then

$$\dot {\bf r}_{ik} = \frac{{\bf p}_{ik}}{m_{ik}}, ~~~~~ \qquad \Dot{\bf h}= \frac{{\bf P}_h}{W}, ~~~~~ \qquad \Dot{\eta} = \frac{p_{\eta}}{Q}$$ (3.34)

$$\dot {\bf p}_{ik} = {\bf h}^{-1} \dot {\bf p}_{ik}- \overleftrightarrow{\bf G}^{-1}\d... ...verleftrightarrow{\bf G}\dot {\bf p}_{ik}- \frac{p_{\eta}}{Q}\dot {\bf p}_{ik},$$ (3.35)

$$\dot {\overleftrightarrow{\bf P}}_h = \left ( {\boldmath\overleftrightarrow{\mathcal{V}}}+ \boldmath\ov... ...{K}}- {\bf h}^{-1} P_{ext} \det {\bf h}\right ) - \frac{p_{\eta}}{Q} {\bf P}_h,$$ (3.36)

$$\Dot{p}_{\eta} = {\mathcal{F}}_{\eta}$$ (3.37)

where the quantities ${\cal V}, {\cal K}, {\cal F}_{\eta}$ depend now on the atomic coordinates

$${\boldmath\overleftrightarrow{\mathcal{V}}} = \sum_{i=1k}^{N} {\bf f}_{ik}{\bf s}_{ik}^{t}$$

$$\boldmath\overleftrightarrow{\mathcal{K}} = \sum_{i=1}^{N} M_i \left ( {\bf h}{\bf s}_{ik} \right ) {\bf s}_{ik}^{t}$$ (3.38)

$${\cal F}_{\eta} = \frac{1}{2} \sum_{i=1}^{N} \sum_{k=1}^{n_i} \frac{ {\bf p}_{ik}^t {\bf p}_{ik}}{m_{ik}} - gk_B T .$$ (3.39)

In case of atomic, Eq. (3.34), or molecular scaling, Eq. (3.1), the internal pressure entering in Eqs. (3.26,3.37) is then

$$P_{int} = \langle P_{atom} \rangle = \left \langle { {1 \over 3V} } \sum_i ... ...ver { m_{ik} } } + {\bf r}_{ ik } \bullet {\bf f}_{ ik } \right) \right \rangle$$ (3.40)

$$P_{int} = \langle P_{mol} \rangle = \left \langle { {1 \over 3V} } \sum_i \... ...f P}^2_{i} } \over M }_i + {\bf R}_{i} \bullet {\bf F}_i \right) \right \rangle$$ (3.41)

respectively. Where the molecular quantities can be written in term of the atomic counterpart according to:

$${\bf R}_{i} = { {1 \over M}_{i} } \sum_k m_{ik} {\bf r}_{ ik }$$ (3.42)

$${\bf P}_{i} = \sum_k {\bf p}_{ ik }$$ (3.43)

$${\bf F}_{i} = \sum_k {\bf f}_{ ik }$$ (3.44)

The equation of motion for the barostat in the two cases, Eqs.(3.37, 3.26), has the same form whether atomic or molecular scaling is adopted. The internal pressure in the former case is given by Eq. (3.41) and in the latter is given by Eq. (3.42). The two pressures, Eqs. (3.41,3.42), differ instantaneously. Should the difference persist after averaging, then it would be obvious that the equilibrium thermodynamic state in the $N{\bf P}T$ ensemble depends on the scaling method. The two formulas (3.41,3.42) are fortunately equivalent. To prove this statement, we closely follow the route proposed by H. Berendsen and reported by Ciccotti and Ryckaert [98] and use Eqs. (3.43-3.45) to rearrange Eq. (3.42). We obtain

$$\sum_i \langle {\bf R}_{i} \bullet {\bf F}_{i} \rangle = \sum_i { {1 \over M}_i \sum_{ kl } \langle m_{ik} {\bf r}_{ ik } \bullet {\bf f}_{ il } \rangle} .$$ (3.45)

Adding and subtracting $m_{ik} {\bf r}_{ il } \bullet {\bf f}_{ il }$, we get

$$= \sum_i { {1 \over M}_i } \sum_{ kl } \langle m_{ik} ({\bf r}_{ ik... ...) \bullet {\bf f}_{ il } + m_{ik} {\bf r}_{ il } \bullet {\bf f}_{ il } \rangle$$ (3.46)

which can be rearranged as

$$= \sum_i {1 \over M_i} \left\{ \sum_{ kl } \left[ {1 \over 2} \lang... ... \right ] + \sum_l \langle {\bf r}_{ il } \bullet {\bf f}_{il} \rangle \right\}$$ (3.47)

using the newton law ${\bf f}_{ ik } = m_{ik} {\bf a}_{ ik }$, where ${\bf a}_{ik}$ is the acceleration, we obtain

$$= \sum_i {1 \over M_i} \left \{ \sum_{ kl } \left[ {1 \over 2} \lan... ...ght ] + \sum_l \langle {\bf r}_{ il } \bullet {\bf f}_{ il } \rangle \right \}.$$ (3.48)

The first term in the above equation can be decomposed according to:

$$({\bf r}_{ ik } - {\bf r}_{ il } ) \bullet ({\bf a}_{ il } - {\bf a}_{ ik } ) = { d \over { dt } } \left[ ({\bf r}_{ ik } - {\bf r}_{ il } ) \bul... ...\bf v}_{ il } - {\bf v}_{ ik } ) \right] + ({\bf v}_{ il } - {\bf v}_{ ik } )^2$$ (3.49)

The first derivative term on the right hand side is zero rigorously for rigid molecules or rigid group and is zero on average for flexible molecules or groups, assuming that the flexible molecules or groups do not dissociate. This can be readily seen in case of ergodic systems, by evaluating directly the average of this derivatives as

$$\langle { d \over { dt } } \left[ ({\bf r}_{ ik } - {\bf r}_{ il } ) \bullet ({\bf v}_{ il } - {\bf v}_{ ik } ) \right] \rangle = \lim_{\tau \rightarrow \infty} {1 \over \tau} \int_{0}^{\infty} {... ...{ ik } - {\bf r}_{ il } ) \bullet ({\bf v}_{ il } - {\bf v}_{ ik } ) \right] dt$$ (3.50)

$$= \lim_{\tau \rightarrow \infty} {1 \over \tau} \left [ ({\bf r}_{ ... ...l }(\tau )) \bullet ({\bf v}_{ il }(\tau) - {\bf v}_{ ik } (\tau)) + C \right ]$$ (3.51)

So if the quantity ${\bf r}_{ikl}(\tau) \bullet {\bf v}_{ilk}(\tau)$ remains bounded (which is true if the potential is not dissociative, since $k,l$ refers to the same molecule $i$), the average in Eq. (3.52) is zero^3.8. Thus, we can rewrite the average of Eq. (3.49) as

$$\langle \sum_i {\bf R}_{i} \bullet {\bf F}_{i} \rangle = \sum_i { 1 \over { 2 M_i } } \sum_{ kl} m_{ik} m_{il} \langle ({\... ...2 \rangle + \sum_{ ik } \langle {\bf r}_{ ik } \bullet {\bf f}_{ ik } \rangle .$$ (3.52)

The first term on the right hand side of the above equation can be further developed obtaining the trivial identity:

$$\sum_{kl} m_{ik} m_{il} \langle ({\bf v}_{ il } - {\bf v}_{ ik } )^2 \rangle = \sum_{ kl } m_{ik} m_{il} \langle {\bf v}_{ ik }^2 \rangle+ \sum_{ kl } m_{ik} m_{il} \langle {\bf v}_{ il }^2 \rangle -$$

$$- \sum_{ kl } 2 m_{ik} m_{il} \langle {\bf v}_{ il } \bullet {\bf v}_{ ik } \rangle$$ (3.53)

$$= 2 M_i \sum_k m_{ik} \langle {\bf v}_{ ik }^2 \rangle - 2 \langle {\bf P}_{i}^2 \rangle$$ (3.54)

Substituting Eq. (3.55) in Eq. (3.53) we get

$$\langle \sum_i {\bf R}_{i} \bullet {\bf F}_{i} \rangle = \sum_{ ik } m_{ik} \langle {\bf v}_{ ik }^2 \rangle - { 1 \over {... ...}^2 \rangle + \sum_{ ik } \langle {\bf r}_{ ik } \bullet {\bf f}_{ ik } \rangle$$ (3.55)

Substituting Eq. (3.56) into Eq. (3.42) leads speedily to (3.41) which completes the proof. As a consequence of the above discussion, it seems likely that both the equilibrium and non equilibrium properties of the MD system are not affected by coordinate scaling. We shall see later that this is actually the case.