Skip to content

Python Simulation Polyatomic Molecules

Jump to bottom
LonelyProf edited this page Jun 18, 2019 · 4 revisions

Polyatomic Lennard-Jones program

The program mc_nvt_poly_lj.py conducts Monte Carlo simulations of a system of rigid molecules composed of Lennard-Jones interaction sites. For simplicity the sites are taken to be identical, although the program is easily generalized. Molecular orientations are represented by quaternions, which are used to calculate the rotation matrix and hence the interaction site positions.

We test this with the three-site model of orthoterphenyl, a fragile glassformer, described in the following publications amongst others

We compare with the results of Mossa et al (2002). The sites are arranged at the vertices of an isosceles triangle with bond angle 75 degrees, LJ parameters ε = 5.276 kJ mol-1, σ=0.483nm, and two equal bonds of length σ. The program employs the usual reduced units based on ε and σ and in these units the potential cutoff of Mossa et al (2002) is Rc=2.612; the pair potential is Lennard-Jones with a shifted-force correction term, linear in r, to make the potential and its derivative vanish at r=Rc. Apart from the temperatures and much longer runs (10 blocks of 20000 steps each, the same as the Fortran tests), default program parameters were used throughout.

Tests were performed at ρ=0.32655 which is equivalent to ρ4=1.108g cm-3 in Mossa et al (2002). Comparisons of potential energy (PE=E-3 T converted to kJ/mol with a factor 5.276) were made with the fit given by eqn (23) of that paper. Note that ε/kB≈635 K.

T E P T (K) PE (kJ/mol) PE (kJ/mol) eqn (23)
0.5 -12.995(2) 1.763(9) 317 -76.48(2) -76.945
1.0 -9.813(15) 5.86(4) 635 -67.60(8) -67.634
1.5 -6.87(1) 9.33(2) 952 -59.99(5) -60.011
2.0 -4.10(1) 12.30(3) 1270 -53.29(5) -53.265

A second set of tests was performed at T=0.6≈380K at the specified densities ρ1, … ρ5 of Mossa et al (2002). A set of starting configurations is provided in the Data repository. Here the excess pressure (P(ex)=P-ρT converted to MPa with a factor 77.75 based on the values of ε and σ) is compared with the fit given by eqn (28) and the coefficients in Table III of Mossa et al (2002). NB the volumes to insert into the equation are those of their Table I, which are specific to their system size. In addition their eqn (29) with coefficients in Table V is a fit to their potential energy, which we calculate from the simulation as described above.

Id ρ E P P(ex) (MPa) P(ex) (MPa) eqn (28) PE (kJ/mol) PE (kJ/mol) eqn (29)
1 0.30533 -11.617(7) 0.47(2) 22(2) 19.077 -70.79(4) -70.818
2 0.31240 -11.927(4) 1.00(2) 63(2) 60.143 -72.42(2) -72.289
3 0.31918 -12.16(1) 1.66(3) 114(2) 112.798 -73.65(5) -73.601
4 0.32655 -12.393(6) 2.52(2) 181(2) 177.222 -74.88(3) -74.886
5 0.33451 -12.498(3) 3.82(1) 281(1) 253.510 -75.44(2) -75.825

In making these comparisons, our limited run length (10 blocks of 20000 sweeps each) should be borne in mind, since this system can show sluggish behaviour. The MD simulations of Mossa et al (2002) are reported to extend to several hundred nanoseconds (of order 107 MD timesteps) at the lowest temperatures. It should be noted that this Python code is quite slow compared to the simpler atomic LJ examples.

For comparison we provide a molecular dynamics code md_nvt_poly_lj.py for the same model. The program takes the molecular mass M to be unity. Mossa et al (2002) ascribe a notional mass of 78u to each of the three LJ sites, so M≈3.9×10-25kg. Combined with the above values of ε and σ, this gives a time scale (M/ε)1/2σ ≈ 3.22 ps. The timestep of δt=0.01 ps used by Mossa et al (2002) corresponds to the default value in the program dt=0.003 in these units. By default, the program simulates the constant-NVE ensemble, but there is an option to simulate at constant NVT by velocity randomization (Andersen thermostat). If the latter option is selected, the program will read configurations in the same format as mc_nvt_poly_lj.py (positions and quaternions only), selecting random initial velocities and angular momenta, which can be convenient.

By default the program calculates the inertia tensor from the LJ site bond vectors, assuming equal masses. For simplicity it is assumed that the bond vectors are defined such that the principal axes of the inertia tensor coincide with the xyz axes of the molecular coordinate system, with the centre of mass at the origin; it is always possible to arrange this. In general, the three principal moments of inertia will all be different, so the molecule is an asymmetric top. The MD algorithm for rotation is a symplectic one in which a kick propagator advances the space-fixed angular momenta, using the torque on each molecule, and a succession of drift steps implement free rotation about each of the principal axes. This is described in the text, section 3.3; see

The results below are for test runs in both constant-NVE and constant-NVT ensembles, at (approximately) the same state points as those given above. All runs were 10×5000 steps in length and used program defaults, except for t_interval=1 and the specified temperature in the NVT case. Because of the slow execution of the Python code, these runs are significantly shorter than the comparable Fortran examples, and very much shorter than the runs of Mossa et al (2002). For constant-NVE runs we report RMS energy fluctuations, and T is the average translational temperature.

ρ T E P E(RMS)
0.32655 0.5 -13.035(7) 1.66(2)
0.32655 0.5025(9) -13.0355 1.606(9) 1.14×10-8
0.32655 1.0 -9.74(1) 6.03(3)
0.32655 1.013(1) -9.7372 5.975(7) 1.01×10-7
0.32655 1.5 -6.81(1) 9.41(4)
0.32655 1.508(1) -6.8136 9.407(5) 3.64×10-7
0.32655 2.0 -4.16(3) 12.26(7)
0.32655 1.986(2) -4.1570 12.291(8) 9.04×10-7
ρ T E P E (RMS)
0.30533 0.6 -11.624(6) 0.48(2)
0.30533 0.6018(7) -11.6239 0.453(6) 1.33×10-8
0.31240 0.6 -11.902(9) 1.13(2)
0.31240 0.604(1) -11.9018 1.049(7) 1.49×10-8
0.31918 0.6 -12.206(7) 1.59(2)
0.31918 0.596(1) -12.2065 1.63(1) 1.64×10-8
0.32655 0.6 -12.379(8) 2.50(2)
0.32655 0.599(1) -12.3793 2.58(1) 1.93×10-8
0.33451 0.6 -12.541(4) 3.66(1)
0.33451 0.601(1) -12.5411 3.70(1) 2.40×10-8
Clone this wiki locally