Fortran Simulation Chains
The program mc_chain_nvt_cbmc_lj 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 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 100000 steps.
For lower temperatures (below 0.40), longer runs (10 blocks of 1000000
steps) were used.
For temperatures higher than 1.0,
the bond length fluctuations become unphysically large for this value of kspring.
| T | PE | Rg | Cv(ex) |
|---|---|---|---|
| 0.26 | -2.044(7) | 1.074(1) | 3.3(2) |
| 0.28 | -1.967(5) | 1.086(1) | 4.16(9) |
| 0.29 | -1.905(4) | 1.096(1) | 4.6(1) |
| 0.30 | -1.893(7) | 1.098(1) | 4.71(9) |
| 0.31 | -1.819(3) | 1.1105(8) | 4.34(6) |
| 0.32 | -1.789(2) | 1.1162(5) | 4.2(1) |
| 0.33 | -1.742(3) | 1.1254(5) | 4.05(1) |
| 0.34 | -1.705(4) | 1.133(1) | 3.7(1) |
| 0.35 | -1.672(3) | 1.140(1) | 3.49(8) |
| 0.40 | -1.50(1) | 1.173(3) | 2.51(17) |
| 0.45 | -1.40(1) | 1.201(3) | 2.26(8) |
| 0.50 | -1.297(8) | 1.224(1) | 1.99(2) |
| 1.00 | -0.438(2) | 1.538(2) | 1.371(3) |
At the lowest temperatures, the acceptance rate of CBMC moves (with the default parameters) was around 2%, while at T=0.35 it was around 11%, increasing further at higher temperatures. The results are broadly in agreement with Calvo et al (2002) showing a similar sized peak in Cv, although at a somewhat lower temperature (0.30 as opposed to 0.35).
Here we give analogous results for the program default spring constant of kspring=400.
| T | PE | Rg | Cv(ex) |
|---|---|---|---|
| 0.26 | -2.076(3) | 1.069(1) | 2.13(5) |
| 0.28 | -2.027(5) | 1.075(1) | 2.91(15) |
| 0.29 | -1.993(4) | 1.081(1) | 3.41(8) |
| 0.30 | -1.958(4) | 1.087(1) | 3.79(8) |
| 0.31 | -1.915(4) | 1.094(1) | 4.40(6) |
| 0.32 | -1.859(6) | 1.105(1) | 4.74(8) |
| 0.33 | -1.816(4) | 1.114(1) | 4.94(6) |
| 0.34 | -1.765(5) | 1.125(1) | 4.90(9) |
| 0.35 | -1.719(4) | 1.135(1) | 4.75(7) |
| 0.40 | -1.505(5) | 1.183(1) | 3.06(14) |
| 0.45 | -1.379(3) | 1.212(1) | 2.32(6) |
| 0.50 | -1.266(3) | 1.240(1) | 2.03(3) |
| 1.00 | -0.459(1) | 1.512(1) | 1.233(6) |
| 2.00 | 0.387(2) | 1.850(1) | 0.591(3) |
| 5.00 | 1.986(3) | 2.035(2) | 0.465(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.
The default run lengths are fairly modest here: 10 blocks,
each consisting of 100000 steps of length δt=0.002.
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 |
|---|---|---|---|---|---|---|
| -2.0246 | 0.2485(2) | 1.06374(4) | 2.176(6) | 0.2475(1) | 1.06450(3) | 2.172(7) |
| -1.9145 | 0.296(2) | 1.073(1) | 2.38(8) | 0.2989(3) | 1.0724(1) | 2.27(2) |
| -1.6145 | 0.345(4) | 1.125(2) | 3.14(8) | 0.347(2) | 1.125(1) | 3.16(6) |
| -1.3495 | 0.404(1) | 1.182(2) | 2.39(1) | 0.404(2) | 1.183(2) | 2.50(2) |
| -1.2195 | 0.451(1) | 1.207(1) | 2.36(2) | 0.449(2) | 1.210(1) | 2.34(2) |
| -1.0968 | 0.499(2) | 1.234(1) | 2.28(1) | 0.503(2) | 1.231(2) | 2.28(2) |
| -0.1244 | 1.009(5) | 1.471(5) | 2.04(2) | 1.024(6) | 1.455(7) | 2.01(2) |
| 1.0456 | 2.008(5) | 1.754(9) | 1.653(3) | 2.009(7) | 1.753(8) | 1.652(3) |
| 3.6459 | 4.996(4) | 1.889(7) | 1.534(1) | 4.988(3) | 1.901(6) | 1.534(1) |
Results for kspring=10000 system using MTS:
| E | T | Rg | Cv |
|---|---|---|---|
| -1.7734 | 0.2496(2) | 1.0695(1) | 2.96(4) |
| -1.3444 | 0.301(2) | 1.144(2) | 5.5(3) |
| -1.1394 | 0.350(2) | 1.173(2) | 3.70(5) |
| -0.9494 | 0.399(1) | 1.199(1) | 3.19(5) |
| -0.7694 | 0.448(1) | 1.230(2) | 3.09(3) |
| -0.5943 | 0.497(2) | 1.262(3) | 3.04(5) |
| 0.7857 | 1.000(4) | 1.467(12) | 2.50(3) |
| 2.8858 | 1.98(2) | 1.752(14) | 2.08(6) |
| 8.3859 | 5.04(4) | 1.904(2) | 1.94(2) |
Results for kspring=400 system using MTS:
| E | T | Rg | Cv |
|---|---|---|---|
| -1.7934 | 0.2487(2) | 1.0646(1) | 2.89(5) |
| -1.6250 | 0.299(1) | 1.0753(4) | 3.24(7) |
| -1.2900 | 0.346(5) | 1.127(3) | 5.7(3) |
| -0.9942 | 0.401(2) | 1.179(2) | 3.58(8) |
| -0.8042 | 0.451(2) | 1.208(2) | 3.21(4) |
| -0.6558 | 0.497(2) | 1.224(2) | 3.03(3) |
| 0.7565 | 0.995(3) | 1.447(6) | 2.55(3) |
| 2.9036 | 2.006(6) | 1.757(9) | 2.08(1) |
| 8.3488 | 5.00(1) | 1.92(1) | 1.97(1) |
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.
- Molecular dynamics is not expected to thoroughly explore the energy landscape at low temperatures, giving instead (typically) quasi-harmonic vibrations in a single basin.
For the hard-sphere square-well chain, the aim was to show the operation of the Wang-Landau method.
In mc_chain_wl_sw 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 in the usual way
(note the non-default value of the bond length).
Default parameters are used in mc_chain_wl_sw,
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.7629 | 0.7518 |
| -2.0 | 0.6661 | 0.7014 | 0.6683 |
| -3.0 | 0.2108 | 0.2308 | 0.2108 |
| -4.0 | -0.4433 | -0.4152 | -0.4449 |
| -5.0 | -1.3484 | -1.3316 | -1.3444 |
| -6.0 | -2.4438 | -2.4256 | -2.4322 |
| -7.0 | -3.6832 | -3.6634 | -3.6733 |
| -8.0 | -5.8548 | -5.8440 | -5.8620 |
| -9.0 | -8.4766 | -8.4050 | -8.4733 |
| -10.0 | -14.9981 | -14.6824 | -15.0295 |
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 100000 steps;
this was increased to 10 blocks of 500000 steps for temperatures
below 0.25.
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.
As the chain length increases, the energy landscape becomes more challenging.
For N=13, with the same bond length of σ,
and nonbonded interaction range of 1.5σ,
sensible results may still be achieved
with the simple example program mc_chain_wl_sw.
Once more, as a reference for comparison,
we ran a set of canonical ensemble calculations with mc_chain_nvt_sw,
using the same parameters described above.
The results are shown on the left of the following table.
| T | PE | Rg | Cv(ex) | PE | Rg | Cv(ex) |
|---|---|---|---|---|---|---|
| method | NVT | NVT | NVT | WL | WL | WL |
| 0.15 | -2.81(1) | 1.070(2) | 1.1(2) | -2.814 | 1.068 | 2.053 |
| 0.18 | -2.759(8) | 1.072(2) | 2.2(2) | -2.744 | 1.073 | 2.498 |
| 0.20 | -2.699(8) | 1.077(2) | 2.4(1) | -2.694 | 1.078 | 2.366 |
| 0.22 | -2.645(4) | 1.082(1) | 2.16(7) | -2.649 | 1.082 | 2.226 |
| 0.25 | -2.586(8) | 1.090(2) | 2.15(9) | -2.584 | 1.089 | 2.121 |
| 0.30 | -2.482(6) | 1.104(2) | 1.97(6) | -2.481 | 1.103 | 2.010 |
| 0.50 | -2.127(2) | 1.161(1) | 1.50(2) | -2.128 | 1.161 | 1.491 |
| 1.00 | -1.547(2) | 1.318(1) | 1.024(6) | -1.543 | 1.319 | 1.020 |
| 2.00 | -0.885(1) | 1.658(1) | 0.330(2) | -0.883 | 1.660 | 0.332 |
| 5.00 | -0.566(1) | 1.889(1) | 0.0318(1) | -0.565 | 1.890 | 0.032 |
Results obtained from a run of the Wang-Landau program mc_chain_wl_sw,
using the same model, are given on the right of the table above.
The program was run with default parameters,
except that the flatness criterion was set at 80%.
The results are from the histograms produced in the 20th stage.
This analysis can also be performed (for any desired temperature) by the program wl_hist, after the run.
The results are generally in good agreement with the canonical ensemble test runs.
The most significant discrepancies are in the heat capacities at the lowest two temperatures,
which reflects the poor sampling of the canonical ensemble program (with these basic MC moves),
and the sampling problems about to be discussed.
This particular test run illustrated one drawback of the simplest Wang-Landau implementation:
two low-lying energies (corresponding to q=38 and 39 square-well interactions) were discovered
during the very last stage, in which ds is very small.
Accordingly, the system remained stuck in these low-energy states for a very long time,
until their tabulated entropy S(q) reached a high enough value to allow q to change;
even then, the final weight of the lowest state in the final "flat" histogram
could not be considered completely reliable.
The overall run length was of order 90000 blocks of 10000 steps each,
most of it spent in stage 20.
Repeating the run with the same parameters typically produces different results,
depending on whether, and when, these low-lying states are discovered.
This affects the canonical ensemble results reconstructed from the histograms
through wl_hist, particularly at the lower temperatures,
while the higher temperatures are largely unaffected.
From a single run, or a few runs, there might be no indication of anything wrong.
It is always a danger with any Monte Carlo method, including Wang-Landau,
that inaccessible configurations will not be sampled.
Various improvements of the method may be found in the literature.
(The reader is invited to work out the structure of the 13-atom chain configuration
with hard sphere diameter and bond length both equal to σ
having 39 nonbonded interactions within range 1.5σ.
Hint: it is not one of the highest-symmetry clusters such as the icosahedron;
however, it is based on a fairly symmetric cluster,
with small distortions due to the fixed bond length and the need to maximise interactions.)
To obtain more reliable results it would be advisable to add a
production run at the end of the mc_chain_wl_sw program,
during which the weights are no longer adjusted,
allowing averages to be generated using a proper Markov chain.