Check force on energy

It is common practise to use the potential energy as a collective energy. Some MD codes thus pass the potential energy to PLUMED and PLUMED can then apply forces on this collective variable. We test that any forces that PLUMED applies on the potential energy are correctly passed back to the MD code by doing the following test. We first run a short simulation at $T$ K with a timestep of $\tau$ ps. During the course of this simulation we monitor the potential energy using the following PLUMED input:

Click on the labels of the actions for more information on what each action computes
tested on2.10
tested onmaster
e: 
ENERGY
Calculate the total potential energy of the simulation box. More details
v:
VOLUME
Calculate the volume of the simulation box. More details
PRINT
Print quantities to a file. More details
ARG
the labels of the values that you would like to print to the file
=e,v
FILE
the name of the file on which to output these quantities
=energy1

We then run a second simulation (starting from identical conditions) at a temperature of $T\alpha$ and with a timestep of $\tau/\sqrt(\alpha)$. The thermostat and barostat relaxation times are similarly divided by $\sqrt(\alpha)$. In the tests that are run on this website we set $\sqrt(\alpha)=1.1$. The PLUMED file above is used when this test is run but a different time series of energy values is recorded as the MD parameters in this second simulation are different.

If PLUMED is working correctly we should be able to recapture the time series of energy values for the first simulation by running an MD simulation with the modified parameters that were used in the second simulation and the following PLUMED input file:

Click on the labels of the actions for more information on what each action computes
tested on2.10
tested onmaster
e: 
ENERGY
Calculate the total potential energy of the simulation box. More details
v:
VOLUME
Calculate the volume of the simulation box. More details
# slope is such that
PRINT
Print quantities to a file. More details
ARG
the labels of the values that you would like to print to the file
=e
FILE
the name of the file on which to output these quantities
=energy2 # slope should be (alpha-1)=0.21
RESTRAINT
Adds harmonic and/or linear restraints on one or more variables. More details
AT
the position of the restraint
=0.0
ARG
the values the harmonic restraint acts upon
=e
SLOPE
specifies that the restraint is linear and what the values of the force constants on each of the variables are
=0.21

In other words, when forces are passed correctly the time series for the energies and volumes from the first and third of these calculations should be identical.

To determine if PLUMED passes this test we calculate the difference between the time series that were observed in the first and third simulations described above. We then divide this by the difference between the first and second time series.

An NPT version of this calculation is performed as well as an NVT calculation if the virial is passed to PLUMED.

Trajectories

  1. Input and output files for the unpeturbed calculation are available in this zip archive

  2. Input and output files for the peturbed calculation are available in this zip archive

  3. Input and output files for the peturbed calculation in which a PLUMED restraint is used to undo the effect of the changed MD parameters are available in this zip archive

Results

Original With PLUMED Effect of peturbation % Difference
-18174.8223 11.6346 -18174.8223 11.6346 0.0000 0.0000 0.0000 0.0000
-18171.8203 11.6346 -18193.2891 11.6346 6.8730 0.0000 312.3615 0.0000
-18195.5430 11.6475 -18232.3770 11.6321 24.0508 0.0054 153.1509 284.9140
-18170.8691 11.6475 -18282.1035 11.6321 23.5449 0.0054 472.4347 284.9140
-18161.4902 11.6475 -18317.6230 11.6321 30.8164 0.0054 506.6548 284.9140
-18147.2266 11.6475 -18336.5059 11.6321 36.4531 0.0054 519.2403 284.9140
-18132.9375 11.6475 -18339.7812 11.6321 41.5020 0.0054 498.3952 284.9140
-18124.1738 11.6475 -18334.1953 11.6321 47.0840 0.0054 446.0572 284.9140
-18124.1035 11.6475 -18327.7793 11.6321 53.1387 0.0054 383.2911 284.9140
-18131.8105 11.6475 -18325.6777 11.6321 58.0391 0.0054 334.0288 284.9140
-18142.9453 11.6475 -18328.7461 11.6321 59.7559 0.0054 310.9332 284.9140
-18152.4590 11.6475 -18334.7852 11.6321 57.4980 0.0054 317.0998 284.9140
-18158.1934 11.6475 -18341.8340 11.6321 53.4453 0.0054 343.6047 284.9140
-18162.9258 11.6475 -18350.7266 11.6321 52.4082 0.0054 358.3423 284.9140
-18172.3887 11.6475 -18364.4219 11.6321 58.3555 0.0054 329.0749 284.9140
-18190.2090 11.6475 -18384.2441 11.6321 70.3828 0.0054 275.6854 284.9140
-18213.3125 11.6475 -18406.3926 11.6321 81.2031 0.0054 237.7742 284.9140
-18232.8691 11.6475 -18422.7383 11.6321 81.5977 0.0054 232.6895 284.9140
-18239.9102 11.6475 -18425.9766 11.6321 67.3926 0.0054 276.0933 284.9140
-18231.1348 11.6475 -18414.6250 11.6321 42.8105 0.0054 428.6099 284.9140

The table below includes some of the results from the calculation. The columns contain:

  1. Time series for the energy and volume that were obtained from the simulation at $T$ K, $x_{md}$.
  2. Time series for the energy and volume that were obtained from the simulation at $\alpha T$ K and in which PLUMED applied a restraint on the energy, $x_{pl}$.
  3. The absolute value of the difference between the time series of energies and volumes that were obtained from the simulations running at $T$ K and $\alpha T$ K, $\vert x_{md}’-x_{md} \vert$. No PLUMED restraints were applied in either of these simulations.
  4. The values of $100\frac{\vert x_{md} - x_{pl}\vert }{ \vert x_{md}’-x_{md} \vert}$.

If the PLUMED interface is working correctly the first two sets of numbers should be identical and the final column should be filled with zeros.

Graphical representation (beta)

A visualization of the table above:
engvir_master