Check virial contribution

If you are running simulations at constant pressure then the virial forces cause the cell parameters to change with time. Any CVs calculated by PLUMED contribute to these virial forces and PLUMED must, therefore, have a mechanism to pass virial forces back to the MD code.

To debug this mechanism we run a constant pressure simulation at 1 bar using the MD code. During this simulation we use the following PLUMED input to monitor the cell volume:

Click on the labels of the actions for more information on what each action computes

vThe VOLUME action with label v calculates the following quantities:
 Quantity   
 Type   
 Description  
v
scalar
the volume of simulation box: 
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
=v 
FILE
the name of the file on which to output these quantities
=volume

We then run a second constant pressure MD simulation at a pressure of 1001 bar and the input above.

If the virial has been implemented correctly within PLUMED the following PLUMED restraint will apply a negative pressure of 1000bar, which should compensate the fact that the second calculation was run at higher pressure. We thus run a third MD calculation with the following input file:

Click on the labels of the actions for more information on what each action computes

vThe VOLUME action with label v calculates the following quantities:
 Quantity   
 Type   
 Description  
v
scalar
the volume of simulation box: 
VOLUME
Calculate the volume of the simulation box. More details
 
# slope should be just 10 times the Avogadro constant:
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
=v 
SLOPE
 specifies that the restraint is linear and what the values of the force constants on each of the variables are
=-60.2214129
The RESTRAINT action with label  calculates the following quantities:
 Quantity   
 Description  
.bias
the instantaneous value of the bias potential
.force2
the instantaneous value of the squared force due to this bias potential
PRINT
Print quantities to a file. More details
 
ARG
the labels of the values that you would like to print to the file
=v 
FILE
the name of the file on which to output these quantities
=volume2

The time series for the volumes that are output by the files volume and volume2 above should thus be close to identical.

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. The time series for the volume that was obtained from the simulation in that was performed at 1 bar, $x_{md}$.
2. The time series for the volume that was obtained from the simulation that was performed at 1001 bar and in which PLUMED applied a restraint on the volume, $x_{pl}$.
3. The absolute value of the difference between the time series of volumes that were obtained from the simulations running at 1001 bar and 1 bar, $\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.