Python Simulation Chains
The program mc_chain_nvt_cbmc_lj.py simulates a single Lennard-Jones chain,
where the atoms are linked by harmonic springs.
There are, therefore,
both bonded and non-bonded interactions,
the former being used to select atom positions,
and the latter appearing in Rosenbluth weights,
which govern the acceptance/rejection of moves.
For comparison with the paper of Calvo, Doye and Wales, J Chem Phys, 116, 2642 (2002),
test runs were carried out using N=13 atoms, a bond length of 1.122462σ
(prepared using initialize.py with "molecules":"chain" to give random non-overlapping atom positions)
and a rather low spring potential kspring=20.
We only use CBMC moves in this code: for a practical application it would be advisable
to include other kinds of move, for example crankshaft, pivot, and bridging moves.
Replica exchange (as used by Calvo et al) would also improve the sampling at low temperature.
Below we report the excess, potential, energy per atom PE,
and the excess heat capacity per atom Cv(ex),
as well as the radius of gyration Rg.
The program default run length is 10 blocks of 5000 steps,
but the results below were obtained with blocks of 100000 steps,
like the Fortran examples.
For lower temperatures (below 0.40), longer runs (10 blocks of 1000000
steps) were used for the Fortran examples, but these were not attempted here.
For temperatures higher than 1.0,
the bond length fluctuations become unphysically large for this value of kspring.
| T | PE | Rg | Cv(ex) |
|---|---|---|---|
| 0.40 | -1.515(4) | 1.174(2) | 2.39(8) |
| 0.45 | -1.405(3) | 1.197(1) | 2.19(4) |
| 0.50 | -1.297(3) | 1.225(1) | 2.00(2) |
| 1.00 | -0.437(2) | 1.538(1) | 1.38(1) |
Here we give analogous results for the program default spring constant of kspring=400.
| T | PE | Rg | Cv(ex) |
|---|---|---|---|
| 0.40 | -1.538(7) | 1.176(2) | 3.6(2) |
| 0.45 | -1.384(3) | 1.212(1) | 2.44(6) |
| 0.5 | -1.268(1) | 1.239(1) | 2.03(2) |
| 1.00 | -0.463(2) | 1.511(1) | 1.228(6) |
| 2.00 | 0.384(3) | 1.853(2) | 0.580(4) |
| 5.00 | 1.999(3) | 2.035(1) | 0.472(2) |
Similar models were employed in md_chain_nve_lj and md_chain_mts_lj:
N=13 atoms and equilibrium bond length of 1.122462σ.
Here we report results for constrained bond lengths, using the first program,
and for kspring=400 and 10000 (the program default value), using the second program.
Averages below were computed over 10 blocks of 100000 steps of length δt=0.002
(same as for the Fortran examples),
whereas the program default is a more modest 10×10000 steps;
otherwise program default input values were used.
The primary indicator of a correctly-functioning program is energy conservation,
and this was checked in all cases.
Energies were chosen to give average temperatures close to the values used in
the MC simulations above.
Results for constrained system (columns 2:4 RATTLE, columns 5:7 MILC-SHAKE):
| E | T | Rg | Cv | T | Rg | Cv |
|---|---|---|---|---|---|---|
| -1.3495 | 0.401(2) | 1.184(1) | 2.47(3) | 0.402(4) | 1.183(3) | 2.42(2) |
| -1.2195 | 0.449(2) | 1.209(2) | 2.35(1) | 0.449(2) | 1.208(1) | 2.35(2) |
| -1.0968 | 0.504(1) | 1.230(1) | 2.30(2) | 0.504(2) | 1.227(1) | 2.29(2) |
| -0.1244 | 1.016(5) | 1.459(5) | 2.03(2) | 1.010(7) | 1.466(5) | 2.02(2) |
| 1.0456 | 1.999(6) | 1.752(7) | 1.653(2) | 2.006(5) | 1.760(7) | 1.650(2) |
| 3.6459 | 4.996(3) | 1.897(7) | 1.534(1) | 4.993(2) | 1.890(4) | 1.534(1) |
Results for kspring=10000 system using MTS:
| E | T | Rg | Cv |
|---|---|---|---|
| -0.9494 | 0.401(2) | 1.199(2) | 3.27(4) |
| -0.7694 | 0.440(1) | 1.240(1) | 3.24(4) |
| -0.5943 | 0.510(2) | 1.244(2) | 2.90(4) |
| 0.7857 | 1.003(4) | 1.415(7) | 2.63(3) |
| 2.8858 | 2.00(2) | 1.75(2) | 2.02(3) |
| 8.3859 | 4.99(2) | 1.903(5) | 1.98(2) |
Results for kspring=400 system using MTS:
| E | T | Rg | Cv |
|---|---|---|---|
| -0.9942 | 0.399(1) | 1.181(1) | 3.6(1) |
| -0.8042 | 0.449(1) | 1.210(1) | 3.12(5) |
| -0.6558 | 0.496(1) | 1.229(2) | 3.07(5) |
| 0.7565 | 0.996(3) | 1.450(5) | 2.54(2) |
| 2.9036 | 2.014(5) | 1.75(1) | 2.08(2) |
| 8.3488 | 4.999(5) | 1.932(7) | 1.973(5) |
When comparing results with the MC program, several points should be remembered.
- Constraining the bond lengths affects average potential energy, kinetic energy, and heat capacity.
- While we use kspring=10000 to highlight the multiple timestep method, it is quite likely that energy flow between bond vibrations and other degrees of freedom will be inefficient, due to the timescale separation.
- The constant-NVE and constant-NVT ensembles are expected to yield different behaviour around the collapse transition. We do not focus on this region, in the Python tests.
- Molecular dynamics is not expected to thoroughly explore the energy landscape at low temperatures, giving instead (typically) quasi-harmonic vibrations in a single basin. We do not focus on this region, in the Python tests.
For the hard-sphere square-well chain, the aim was to show the operation of the Wang-Landau method.
In mc_chain_wl_sw.py we use pivot and crankshaft moves as well as CBMC regrowth.
In a practical application it would be advisable to include some bridging moves as well.
Reasonably long chains, N=128, have been studied by this method,
and exact results are available for very short chains;
see, for example,
- MP Taylor, J Chem Phys, 118, 883 (2003),
- JE Magee, L Lue, RA Curtis, Phys Rev E, 78, 031803 (2008),
- MP Taylor, W Paul, K Binder, J Chem Phys, 131, 114907 (2009),
who provide references to other simulation work.
For testing purposes our aims are quite modest:
we choose N=6, bond length equal to σ, and a nonbonded interaction range of 1.5σ.
The starting chain configuration can be prepared using initialize.py in the usual way
(note the non-default value of the bond length).
Default parameters are used in mc_chain_wl_sw.py,
including a flatness criterion of 0.9.
The entropy modification constant ds is halved at each stage,
and there are 20 stages.
For this system, the energy range (in units of the well depth) is
E = 0 … -10.
The principal result is the histogram of entropies S(E) produced at the final stage.
For convenience we (arbitrarily) define S(0)=0.
We conduct a set of nine independent WL runs,
and report the results from the two runs with the highest and lowest values of S(-10),
which bracket all the other results in the set,
as a rough indication of the errors.
We compare with the exact values calculated from the density of states
of Taylor (2003), normalized in the same way to make S(0)=0.
| E | S(E) (exact) | S(E) (WL) | S(E) (WL) |
|---|---|---|---|
| 0.0 | 0.0000 | 0.0000 | 0.0000 |
| -1.0 | 0.7521 | 0.7561 | 0.7506 |
| -2.0 | 0.6661 | 0.6655 | 0.6510 |
| -3.0 | 0.2108 | 0.2067 | 0.1984 |
| -4.0 | -0.4433 | -0.4352 | -0.4560 |
| -5.0 | -1.3484 | -1.3837 | -1.3385 |
| -6.0 | -2.4438 | -2.5051 | -2.4393 |
| -7.0 | -3.6832 | -3.7044 | -3.6725 |
| -8.0 | -5.8548 | -5.8702 | -5.8534 |
| -9.0 | -8.4766 | -8.4457 | -8.5045 |
| -10.0 | -14.9981 | -14.5896 | -15.2831 |
As a further check, we ran a set of canonical ensemble calculations for the same system
with mc_chain_nvt_sw at selected temperatures.
The program default is to run for 10 blocks, each of 10000 steps,
but the results here were obtained with 10×500000 steps
for temperatures below 0.25,
and 10×100000 steps for higher temperatures.
The results may be compared with values reconstructed using the
wl_hist program from the simulation histograms.
Below we show the heat capacity per atom from the WL runs (shaded region, bounded by the two runs mentioned above),
from the exact density of states of Taylor (2003) (line),
and from the canonical ensemble calculations (points).
It is also straightforward to compare average energies and radii of gyration, but we do not do that here.
Because of the slow speed of the Python codes, we make no attempt to reproduce the tests conducted with the Fortran versions for the N=13 chains.