The Pressure in Systems with Frozen Atoms

W.Smith ^$\dagger$ and P.M. Rodger ^$\ddagger$,
^$\dagger$Computation Science and Engineering Department,
C.C.L.R.C. Daresbury Laboratory, Daresbury, Warrington WA4 4AD, U.K.
^$\ddagger$Department of Chemistry, University of Warwick, Coventry, U.K.

1 Introduction

In molecular dynamics it is often an advantage to fix the position of some atoms, so that they are unaffected by the normal dynamics of the system. Such `frozen' atoms then form a fixed background to the events occuring in the simulation, but nevertheless exert an influence. This is of use, for example, in studies of catalysis, where atoms remote from the active site exert background forces on the process, but otherwise are not dynamically involved.

However, questions arise concerning the pressure in such systems. How are we to treat the forces derived from these frozen atoms? Do they contribute to the system virial and therefore to the system pressure? In this short note we take a look at these questions.

2 The Virial Theorem

It is instructive first to remind ourselves how pressure is calculated in systems where all the atoms are free to move. The basis of the pressure calculation comes from the Virial Theorem, which is due to Clausius.

We begin by defining a quantity W for a system of N atoms.

$$\sum_{i}^{N} \underline{r}_{i}^{} \underline{\dot{r}}_{i}^{}$$ (1)

The derivative with respect to time is

$${\frac{d}{dt}} \sum_{i}^{N} \underline{\dot{r}}_{i}^{2} \sum_{i}^{N} \underline{r}_{i}^{} \underline{\ddot{r}}_{i}^{}$$ (2)

which, according to Newton, may be written in terms of the force $\underline{f}_{i}^{}$ acting on each atom i as

$${\frac{d}{dt}} \sum_{i}^{N} \underline{\dot{r}}_{i}^{2} \sum_{i}^{N} \underline{r}_{i}^{} \underline{f}_{i}^{}$$ (3)

The force $\underline{f}_{i}^{}$ may be split into internal and external contributions as follows,

$$\underline{f}_{i}^{} \sum_{j}^{N} \prime \underline{f}_{ij}^{} \underline{f}_{i}^{o}$$ (4)

where $\underline{f}_{ij}^{}$ is the force on atom i due to atom j and $\underline{f}_{i}^{o}$ is the external force (i.e. one that does not arise from mutual atom-atom interactions). The dash $\prime$ on the summation indicates that the condition i$\ne$j applies in the sum.

Now we have

$${\frac{d}{dt}} \sum_{i}^{N} \underline{\dot{r}}_{i}^{2} \sum_{i}^{N} \sum_{j}^{N} \prime \underline{r}_{i}^{} \underline{f}_{ij}^{} \sum_{i}^{N} \underline{r}_{i}^{} \underline{f}_{i}^{o}$$ (5)

The ensemble average of this equation may be written as

$$\left\langle\vphantom{ \frac{d}{dt}W }\right. {\frac{d}{dt}} \left.\vphantom{ \frac{d}{dt}W }\right\rangle \Psi \left\langle\vphantom{ \sum_i^N \underline{r}_i \cdot \underline{f}_i^o }\right. \sum_{i}^{N} \underline{r}_{i}^{} \underline{f}_{i}^{o} \left.\vphantom{ \sum_i^N \underline{r}_i \cdot \underline{f}_i^o }\right\rangle$$ (6)

in which we have defined the the Kinetic Energy K as

$$\left\langle\vphantom{ \frac{1}{2}\sum_i^N m_i \underline{\dot{r}}_i^2 }\right. {\textstyle\frac{1}{2}} \sum_{i}^{N} \underline{\dot{r}}_{i}^{2} \left.\vphantom{ \frac{1}{2}\sum_i^N m_i \underline{\dot{r}}_i^2 }\right\rangle$$ (7)

and the Virial of Clausius $\Psi$ as

$$\Psi \left\langle\vphantom{ \sum_i \sum_{j>i}^N \underline{r}_{ij} \cdot \underline{f}_{ij} }\right. \sum_{i}^{} \sum_{j>i}^{N} \underline{r}_{ij}^{} \underline{f}_{ij}^{} \left.\vphantom{ \sum_i \sum_{j>i}^N \underline{r}_{ij} \cdot \underline{f}_{ij} }\right\rangle$$ (8)

in which we have used the well-known equivalence

$$\sum_{i}^{N} \sum_{j}^{N} \prime \underline{r}_{i}^{} \underline{f}_{ij}^{} \sum_{i}^{} \sum_{j>i}^{N} \underline{r}_{ij}^{} \underline{f}_{ij}^{}$$ (9)

Now, for a system confined within a cubic box, the external force of interest is the containing force supplied by the walls, which must balance the pressure of the system. Furthermore it only acts when the colliding atoms make contact with the walls. From these arguments we deduce that

$$\left\langle\vphantom{ \sum_i^N \underline{r}_i \cdot \underline{f}_i^o }\right. \sum_{i}^{N} \underline{r}_{i}^{} \underline{f}_{i}^{o} \left.\vphantom{ \sum_i^N \underline{r}_i \cdot \underline{f}_i^o }\right\rangle$$ (10)

In which P is the pressure, V the cube volume, L the width of the cube and A the area of one face of the cube. The minus sign on the right of this relation arises from the fact that the wall forces act inwards on the system. (It is noteworthy that in this argument it is assumed that the confining cube is sufficiently large that the size of the atoms is negligible by comparison.)

By the Ergodic Hypothesis we may replace an ensemble average by an average over time, thus we may write

$$\left\langle\vphantom{ \frac{d}{dt}W }\right. {\frac{d}{dt}} \left.\vphantom{ \frac{d}{dt}W }\right\rangle {\frac{lim}{\tau \rightarrow \infty}} {\frac{1}{\tau}} \int_{0}^{\tau} \left(\vphantom{ \frac{dW}{dt} }\right. {\frac{dW}{dt}} \left.\vphantom{ \frac{dW}{dt} }\right)$$ =

$${\frac{lim}{\tau \rightarrow \infty}} {\frac{1}{\tau}} \int_{0}^{\tau} \left(\vphantom{ \sum_i^N m_i \underline{r}_i \cdot \underline{\dot{r}}_i }\right. \sum_{i}^{N} \underline{r}_{i}^{} \underline{\dot{r}}_{i}^{} \left.\vphantom{ \sum_i^N m_i \underline{r}_i \cdot \underline{\dot{r}}_i }\right)$$ =

$${\frac{lim}{\tau \rightarrow \infty}} {\frac{1}{\tau}} \left[\vphantom{ \sum_i^N m_i \underline{r}_i \cdot \underline{\dot{r}}_i }\right. \sum_{i}^{N} \underline{r}_{i}^{} \underline{\dot{r}}_{i}^{} \left.\vphantom{ \sum_i^N m_i \underline{r}_i \cdot \underline{\dot{r}}_i }\right]^{\tau}_{0}$$ =

= 0. (11)

The last equality comes from the fact that the sum in square brackets must always be finite, but the denominator $\tau$ becomes infinite. Thus finally we have:

$$\Psi$$ (12)

or in the more usual form

$$\Psi$$ (13)

This is the equation we normally use to calculate the pressure in our simulations.

3 Adaptation to Systems with `Frozen' Atoms

To adapt Clausius' theorem to this case, we begin by considering the first N_f atoms from the N atom system as frozen. What effect do they exert on the remaining system? Clearly forces arise from them which act on the other atoms, but they are not themselves meaningfully acted upon - inasmuch as they are not moved by these forces. The principal effect is to exclude the mobile atoms from the volume that surrounds them, and thus effectively reduce the system volume. Can we use this observation?

Let us write the virial equation for the full N atom system:

$$\left\langle\vphantom{ \sum_i^{N} m_i \underline{\dot{r}}_i^2 }\right. \sum_{i}^{N} \underline{\dot{r}}_{i}^{2} \left.\vphantom{ \sum_i^{N} m_i \underline{\dot{r}}_i^2 }\right\rangle \left\langle\vphantom{ \sum_i \sum_{j>i}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right. \sum_{i}^{} \sum_{j>i}^{N} \underline{r}_{ij}^{} \underline{f}_{ij}^{} \left.\vphantom{ \sum_i \sum_{j>i}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right\rangle$$ (14)

in which we assume the first N_f atoms are frozen and the remainder are free to move as Newton dictates. The first term on the left may now be written as

$$\left\langle\vphantom{ \sum_i^{N} m_i \underline{\dot{r}}_i^2 }\right. \sum_{i}^{N} \underline{\dot{r}}_{i}^{2} \left.\vphantom{ \sum_i^{N} m_i \underline{\dot{r}}_i^2 }\right\rangle \left\langle\vphantom{ \sum_{i>N_f}^{N} m_i \underline{\dot{r}}_i^2 }\right. \sum_{i>N_f}^{N} \underline{\dot{r}}_{i}^{2} \left.\vphantom{ \sum_{i>N_f}^{N} m_i \underline{\dot{r}}_i^2 }\right\rangle$$ (15)

where K is the kinetic energy of the free atoms only.

The virial term may be expanded into the following contributions

$$\left\langle\vphantom{ \sum_i \sum_{j>i}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right. \sum_{i}^{} \sum_{j>i}^{N} \underline{r}_{ij}^{} \underline{f}_{ij}^{} \left.\vphantom{ \sum_i \sum_{j>i}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right\rangle \left\langle\vphantom{ \sum_i \sum_{j>i}^{N_f} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right. \sum_{i}^{} \sum_{j>i}^{N_f} \underline{r}_{ij}^{} \underline{f}_{ij}^{} \left.\vphantom{ \sum_i \sum_{j>i}^{N_f} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right\rangle$$ =

$$\left\langle\vphantom{ \sum_{i>N_f} \sum_{j>i}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right. \sum_{i>N_f}^{} \sum_{j>i}^{N} \underline{r}_{ij}^{} \underline{f}_{ij}^{} \left.\vphantom{ \sum_{i>N_f} \sum_{j>i}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right\rangle$$

$$\left\langle\vphantom{ \sum_i^{N_f} \sum_{j>N_f}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right. \sum_{i}^{N_f} \sum_{j>N_f}^{N} \underline{r}_{ij}^{} \underline{f}_{ij}^{} \left.\vphantom{ \sum_i^{N_f} \sum_{j>N_f}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right\rangle$$ (16)

The first term on the right, contains contributions from the frozen atoms alone. If we insist that these atoms never change their positions it is difficult to see how this term can contribute to the pressure. In thermodynamic terms this subset of interactions does not scale with the system volume and thus cannot contribute to the pressure via the standard relation:

$$\left\{\vphantom{\frac{\partial A}{\partial V} }\right. {\frac{\partial A}{\partial V}} \left.\vphantom{\frac{\partial A}{\partial V} }\right\}_{T}^{}$$ (17)

Thus we fell justified in neglecting this term altogther.

The second term we may identify with the virial of the freely moving atoms:

$$\Psi{^\prime} \left\langle\vphantom{ \sum_{i>N_f} \sum_{j>i}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right. \sum_{i>N_f}^{} \sum_{j>i}^{N} \underline{r}_{ij}^{} \underline{f}_{ij}^{} \left.\vphantom{ \sum_{i>N_f} \sum_{j>i}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right\rangle$$ (18)

The third term involves both frozen and mobile atoms. Can this be identified with a reduction in system volume as suggested above? It turns out that it can.

Consider a single pair of atoms, one frozen and one free. The virial contribution is $\underline{r}_{ij}^{}$ ^. $\underline{f}_{ij}^{}$, which because the force and interatomic vectors have the same direction, reduces simply to r_ijf_ij, if we assume for the moment that the system is comprised of uniform hard spheres, this contribution if finite only when r_ij = $\sigma$, which is the atomic diameter. f_ij then corresponds to an impulse force. From the viewpoint of the frozen atom, the sum of such impulse forces gives rise to a pressure acting on the surface of a sphere of radius $\sigma$ centred on the frozen atom. This pressure is, of course, the system pressure P. Thus we may write for the hard sphere case (by analogy with the wall effect of the containing cube):

$$\left\langle\vphantom{ \sum_i^{N_f} \sum_{j>N_f}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right. \sum_{i}^{N_f} \sum_{j>N_f}^{N} \underline{r}_{ij}^{} \underline{f}_{ij}^{} \left.\vphantom{ \sum_i^{N_f} \sum_{j>N_f}^{N} \underline{r}_{ij} \cdot \underline{f}_{ij} }\right\rangle \sum_{i}^{N_f} \sigma_{i}^{} \pi \sigma_{i}^{2}$$ (19)

Where V_h is the excluded volume around a frozen atom. The pressure force in this case acts in the same direction as the vector r_ij and leads to a positive result. Thus we see that the additional term is effectively a correction to the system volume due to the finite size of the frozen atoms. i.e. the pressure equation is now

$$\Psi{^\prime}$$ (20)

In systems with continuous potentials however, we cannot express the correction so neatly, but at least we see what it means. In practice it is more straightforward to combine the volume correction terms with the virial for the free atoms in 18 (and compute them at the same time).

Thus, in summary, the only operational differences between calculating pressure for a system of mobile atoms and one containing frozen atoms are:

1. When calculating the kinetic energy, the contributions from the frozen atoms are zero.
2. When calculating the virial, assume all forces between frozen atoms are zero, but include all other pair forces.

This is, perhaps, a rather obvious result, but at least it now has some foundation!