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.9
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 input for this action is the scalar output from one or more other actions
=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.9
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 input for this action is the scalar output from one or more other actions
=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 arguments on which the bias is acting
=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

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.

Original result Result with PLUMED Effect of peturbation % Difference
-26662.0234 20.5064 -26662.0234 20.5064 0.0000 0.0000 0.0000 0.0000
-26632.8451 20.5065 -32231.8700 20.5065 4.2609 0.0000 131403.3915 0.0000
-26571.8418 20.5067 -33340.5345 20.5067 17.1626 0.0000 39438.6459 0.0000
-26484.6485 20.5070 -33486.1608 20.5070 38.7570 0.0001 18065.1732 0.0000
-26378.6384 20.5075 -33410.7323 20.5075 69.0190 0.0002 10188.6288 0.0000
-26263.1010 20.5081 -33279.3548 20.5081 107.6309 0.0002 6518.8106 0.0000
-26147.5261 20.5088 -33136.1905 20.5088 153.6553 0.0003 4548.2754 0.0000
-26041.4056 20.5095 -33000.0056 20.5095 206.2019 0.0004 3374.6538 0.0000
-25953.4726 20.5103 -32883.4737 20.5103 263.0785 0.0005 2634.1956 0.0000
-25889.9420 20.5112 -32795.4714 20.5112 321.9729 0.0006 2144.7547 0.0000
-25854.7954 20.5121 -32741.8445 20.5121 379.6169 0.0006 1814.2104 0.0000
-25849.5326 20.5129 -32725.3200 20.5129 433.2041 0.0007 1587.1936 0.0000
-25872.5350 20.5138 -32744.8523 20.5138 479.5418 0.0007 1433.1008 0.0000
-25918.8984 20.5145 -32795.3171 20.5145 514.9060 0.0006 1335.4707 0.0000
-25982.6266 20.5152 -32869.6430 20.5152 537.2445 0.0006 1281.9148 0.0000
-26055.6035 20.5158 -32958.2286 20.5158 544.7532 0.0005 1267.1105 0.0000
-26129.2387 20.5162 -33050.4665 20.5162 536.8589 0.0003 1289.2079 0.0000
-26195.4679 20.5166 -33136.0659 20.5166 514.7496 0.0002 1348.3444 0.0000
-26247.4342 20.5167 -33206.0079 20.5167 480.3899 0.0001 1448.5261 0.0000
-26279.6508 20.5167 -33252.9125 20.5167 437.2449 0.0004 1594.8184 0.0000