Computational Physics

256 Molecular dynamics simulations

att+handt−hafter such a translation can cause severe errors: this should be done

beforemoving the molecule back into the cell.

(a) Implement this algorithm for liquid nitrogen.

The program can be checked by verifying whether the kinetic energies

associated with translational and vibrational degrees of freedom satisfy

equipartition. The total kinetic energyKtotcan be determined as in the argon case

by taking allatomicvelocities into account. From this, the temperature can be

determined asNkBT= 4 / 5 KtotwhereNis the number of molecules. The

translational kinetic energyKtranscan be calculated by taking into account the

molecular velocities (sums of velocities of the two atoms) and the temperature

can be found from this asNkBT= 3 / 2 Ktrans. The average temperatures should be

the same for both procedures.

Check whether this requirement is satisfied.

(b) The virial theorem applies as usual: molecular forces should be used and the

separation occurring in this theorem is the separation between the centres of mass

of the molecules. The correction term is evaluated usingg≡1 for the correlation

function beyond the cut-off distance, where it is assumed thatgis independent of

distance but also of the angular configuration of the molecular pairs.

(c) Calculate the pressure also using the atomic forces (including the constraint

forces), and compare the result with (b).

(d) Calculate the pressure for various temperatures and densities. Cheung and

Powles give extensive data on thermodynamic quantities [46]. The table below

gives some of the data (in reduced units) obtained by Cheung and Powles.

ρ TPU

0.6964 2.86 8.35 −17.16

0.6964 1.72 1.29 −18.68

0.6220 2.70 2.50 −15.82

0.6220 2.17 0.27 −16.30

8.8 [C] In this problem, we consider the implementation of the Andersen method for

simulating a system in the canonical ensemble. Remember that the preferred energy

estimator for the Verlet/leap-frog algorithm is

E=

∑

i

[pi(t+h/ 2 )+pi(t−h/ 2 )]^2

8

+V[R(t)],

whereRis the combined position coordinate of the system which consists of particles

of massm=1. In view of the form of this estimator, it seems sensible to update the

momenta at the same time instances for which we calculate the positions, and it is

convenient to define theith component of the momentum coordinate at timet:

pi(t)=[pi(t+h/ 2 )+pi(t−h/ 2 )]/2.