Action: TD_MULTITHERMAL_MULTIBARIC

Module	ves
Description	Usage
Multithermal-multibaric target distribution (dynamic).

Details and examples

Multithermal-multibaric target distribution (dynamic).

Use the target distribution to sample the multithermal-multibaric ensemble (see first paper cited below). In this way, in a single molecular dynamics simulation one can obtain information about the simulated system in a range of temperatures and pressures. This range is determined through the keywords MIN_TEMP, MAX_TEMP, MIN_PRESSURE, and MAX_PRESSURE. One should also specified the target pressure of the barostat with the keyword PRESSURE.

The collective variables (CVs) used to construct the bias potential must be: 1. the potential energy and the volume or, 2. the potential energy, the volume, and an order parameter.

Other choices of CVs or a different order of the above mentioned CVs are nonsensical. The third CV, the order parameter, must be used when the region of the phase diagram under study is crossed by a first order phase transition (for details see the first paper cited below).

The algorithm will explore the free energy at each temperature and pressure up to a predefined free energy threshold $\epsilon$ specified through the keyword THRESHOLD (in kT units). If only the energy and the volume are being biased, i.e. no phase transition is considered, then THRESHOLD can be set to around 5. If also an order parameter is used then the THRESHOLD should be greater than the barrier for the transformation in kT. For small systems undergoing a freezing transition THRESHOLD is typically between 20 and 50.

It is also important to specify the number of intermediate temperatures and pressures to consider. This is done through the keywords STEPS_TEMP and STEPS_PRESSURE. If the number of intermediate temperature and pressures is too small, then holes might appear in the target distribution. If it is too large, the performance will deteriorate with no additional advantage.

We now describe the algorithm more rigorously. The target distribution is given by

$p(E,\mathcal{V},s)= \begin{cases} 1/\Omega_{E,\mathcal{V},s} & \text{if there is at least one } \beta',P' \text{ such} \\ & \text{that } \beta' F_{\beta',P'}(E,\mathcal{V},s)<\epsilon \text{ with} \\ & \beta_1>\beta'>\beta_2 \text{ and } P_1<P'<P_2 \\ 0 & \text{otherwise} \end{cases}$

with $F_{\beta',P'}(E,\mathcal{V},s)$ the free energy as a function of energy $E$ and volume $\mathcal{V}$ (and optionally the order parameter $s$ ) at temperature $\beta'$ and pressure $P'$ , $\Omega_{E,\mathcal{V},s}$ is a normalization constant, and $\epsilon$ is the THRESHOLD. In practice the condition $\beta' F_{\beta',P'}(E,\mathcal{V},s)<\epsilon$ is checked in equally spaced points in each dimension $\beta'$ and $P'$ . The number of points is determined with the keywords STEPS_TEMP and STEPS_PRESSURE. In practice the target distribution is never set to zero but rather to a small value controlled by the keyword EPSILON. The small value is e^-EPSILON.

Much like in the Wang-Landau algorithm that is introduced in the second to last paper cited below or in the multicanonical ensemble that is discussed in the last paper cited below, a flat histogram is targeted. The idea behind this choice of target distribution is that all regions of potential energy and volume (and optionally order parameter) that are relevant at all temperatures $\beta_1<\beta'<\beta_2$ and pressure $P_1<P'<P_2$ are included in the distribution.

The free energy at temperature $\beta'$ and pressure $P'$ is calculated from the free energy at $\beta$ and $P$ using:

$\beta' F_{\beta',P'}(E,\mathcal{V},s) = \beta F_{\beta,P}(E,\mathcal{V},s) + (\beta' - \beta) E + (\beta' P' - \beta P ) \mathcal{V} + C$

with $C$ such that $F_{\beta',P'}(E_m,\mathcal{V}_m,s_m)=0$ with $E_{m},\mathcal{V}_m,s_m$ the position of the free energy minimum. $\beta F_{\beta,P}(E,\mathcal{V},s)$ is not know from the start and is instead found during the simulation. Therefore $p(E,\mathcal{V},s)$ is determined iteratively as done in the well tempered target distribution in the second paper cited below.

The output of these simulations can be reweighted in order to obtain information at all temperatures and pressures in the targeted region of Temperature-Pressure plane. The reweighting can be performed using the action REWEIGHT_TEMP_PRESS.

The multicanonical ensemble (fixed volume) can be targeted using TD_MULTICANONICAL.

Examples

The following input can be used to run a simulation in the multithermal-multibaric ensemble. The region of the temperature-pressure plane that will be explored is 260-350 K and 1 bar- 300 MPa. The energy and the volume are used as collective variables. Legendre polynomials are used to construct the two dimensional bias potential. The averaged stochastic gradient descent algorithm is chosen to optimize the VES functional. The target distribution is updated every 100 optimization steps (200 ps here) using the last estimation of the free energy.

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

# Use energy and volume as CVs
energyThe ENERGY action with label energy calculates the following quantities: Quantity    Type    Description  
energy scalar the internal energy: ENERGYCalculate the total potential energy of the simulation box. More details
volThe VOLUME action with label vol calculates the following quantities: Quantity    Type    Description  
vol scalar the volume of simulation box: VOLUMECalculate the volume the simulation box. More details

# Basis functions
bf1: BF_LEGENDRELegendre polynomials basis functions. More details ORDERThe order of the basis function expansion=10 MINIMUMThe minimum of the interval on which the basis functions are defined=-14750 MAXIMUMThe maximum of the interval on which the basis functions are defined=-12250
The BF_LEGENDRE action with label bf1 calculates somethingbf2: BF_LEGENDRELegendre polynomials basis functions. More details ORDERThe order of the basis function expansion=10 MINIMUMThe minimum of the interval on which the basis functions are defined=6.5 MAXIMUMThe maximum of the interval on which the basis functions are defined=8.25

# Target distribution - 1 bar = 0.06022140857 and 300 MPa = 180.66422571
The BF_LEGENDRE action with label bf2 calculates somethingTD_MULTITHERMAL_MULTIBARICMultithermal-multibaric target distribution (dynamic). More details ...
 MIN_TEMPMinimum energy=260
 MAX_TEMPMaximum energy=350
 MAX_PRESSUREMaximum pressure=180.66422571
 MIN_PRESSUREMinimum pressure=0.06022140857
 PRESSURETarget pressure of the barostat used in the MD engine=0.06022140857
 LABELa label for the action so that its output can be referenced in the input to other actions=td_multi
... TD_MULTITHERMAL_MULTIBARIC

# Bias expansion
The TD_MULTITHERMAL_MULTIBARIC action with label td_multi calculates somethingVES_LINEAR_EXPANSIONLinear basis set expansion bias. This action has hidden defaults. More details ...
 ARGthe labels of the scalars on which the bias will act=energy,vol
 BASIS_FUNCTIONSthe label of the one dimensional basis functions that should be used=bf1,bf2
 TEMPthe system temperature - this is needed if the MD code does not pass the temperature to PLUMED=300.0
 GRID_BINSthe number of bins used for the grid=200,200
 TARGET_DISTRIBUTIONthe label of the target distribution to be used=td_multi
 LABELa label for the action so that its output can be referenced in the input to other actions=b1The VES_LINEAR_EXPANSION action with label b1 calculates the following quantities: Quantity    Type    Description  
b1.bias scalar the instantaneous value of the bias potential
b1.force2 scalar the instantaneous value of the squared force due to this bias potential.
... VES_LINEAR_EXPANSION
VES_LINEAR_EXPANSIONLinear basis set expansion bias. This action uses the defaults shown here. More details ...
 ARGthe labels of the scalars on which the bias will act=energy,vol
 BASIS_FUNCTIONSthe label of the one dimensional basis functions that should be used=bf1,bf2
 TEMPthe system temperature - this is needed if the MD code does not pass the temperature to PLUMED=300.0
 GRID_BINSthe number of bins used for the grid=200,200
 TARGET_DISTRIBUTIONthe label of the target distribution to be used=td_multi
 LABELa label for the action so that its output can be referenced in the input to other actions=b1
 GRID_FMT the format to use when outputting the numbers in the grids=%14.9f
... VES_LINEAR_EXPANSION

# Optimization algorithm
OPT_AVERAGED_SGDAveraged stochastic gradient decent with fixed step size. This action has hidden defaults. More details ...
  BIASthe label of the VES bias to be optimized=b1
  STRIDEthe frequency of updating the coefficients given in the number of MD steps=500
  LABELa label for the action so that its output can be referenced in the input to other actions=o1
  STEPSIZEthe step size used for the optimization=1.0
  FES_OUTPUThow often the FES(s) should be written out to file=500
  BIAS_OUTPUThow often the bias(es) should be written out to file=500
  TARGETDIST_OUTPUThow often the dynamic target distribution(s) should be written out to file=500
  COEFFS_OUTPUT how often the coefficients should be written to file=100
  TARGETDIST_STRIDEstride for updating a target distribution that is iteratively updated during the optimization=100
... OPT_AVERAGED_SGD
The OPT_AVERAGED_SGD action with label o1 calculates the following quantities: Quantity    Description  
o1.value a scalarOPT_AVERAGED_SGDAveraged stochastic gradient decent with fixed step size. This action uses the defaults shown here. More details ...
  BIASthe label of the VES bias to be optimized=b1
  STRIDEthe frequency of updating the coefficients given in the number of MD steps=500
  LABELa label for the action so that its output can be referenced in the input to other actions=o1
  STEPSIZEthe step size used for the optimization=1.0
  FES_OUTPUThow often the FES(s) should be written out to file=500
  BIAS_OUTPUThow often the bias(es) should be written out to file=500
  TARGETDIST_OUTPUThow often the dynamic target distribution(s) should be written out to file=500
  COEFFS_OUTPUT how often the coefficients should be written to file=100
  TARGETDIST_STRIDEstride for updating a target distribution that is iteratively updated during the optimization=100
 COEFFS_FILE the name of output file for the coefficients=coeffs.data
... OPT_AVERAGED_SGD

The multithermal-multibaric target distribution can also be used to explore regions of the phase diagram crossed by first order phase transitions. Consider a system of 250 atoms that crystallizes in the FCC crystal structure. The region of the temperature-pressure plane that will be explored is 350-450 K and 1bar-1GPa. We assume that inside this region we can find the liquid-FCC coexistence line that we would like to obtain. In this case in addition to the energy and volume, an order parameter must also be biased. The energy, volume, and an order parameter are used as collective variables to construct the bias potential. We choose as order parameter the FCCUBIC. Legendre polynomials are used to construct the three dimensional bias potential. The averaged stochastic gradient descent algorithm is chosen to optimize the VES functional. The target distribution is updated every 100 optimization steps (200 ps here) using the last estimation of the free energy.

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

# Use energy, volume and FCCUBIC as CVs
energyThe ENERGY action with label energy calculates the following quantities: Quantity    Type    Description  
energy scalar the internal energy: ENERGYCalculate the total potential energy of the simulation box. More details
volThe VOLUME action with label vol calculates the following quantities: Quantity    Type    Description  
vol scalar the volume of simulation box: VOLUMECalculate the volume the simulation box. More details
fccThe FCCUBIC action with label fcc calculates the following quantities: Quantity    Type    Description  
fcc vector the symmetry function for each of the specified atoms
fcc_morethan scalar the number of colvars that have a value more than a threshold: FCCUBICMeasure how similar the environment around atoms is to that found in a FCC structure. This action is a shortcut and it has hidden defaults. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-256 SWITCHthe switching function that it used in the construction of the contact matrix. Options for this keyword are explained in the documentation for LESS_THAN.={CUBIC D_0=0.4 D_MAX=0.5} MORE_THANcalculate the number of variables that are more than a certain target value. Options for this keyword are explained in the documentation for MORE_THAN.={RATIONAL R_0=0.45 NN=12 MM=24}
fcc: FCCUBICMeasure how similar the environment around atoms is to that found in a FCC structure. This action is a shortcut and uses the defaults shown here. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-256 SWITCHthe switching function that it used in the construction of the contact matrix. Options for this keyword are explained in the documentation for LESS_THAN.={CUBIC D_0=0.4 D_MAX=0.5} MORE_THANcalculate the number of variables that are more than a certain target value. Options for this keyword are explained in the documentation for MORE_THAN.={RATIONAL R_0=0.45 NN=12 MM=24}  PHI The Euler rotational angle phi=0.0 THETA The Euler rotational angle theta=0.0 PSI The Euler rotational angle psi=0.0 ALPHA The alpha parameter of the angular function that is used for FCCUBIC=3.0
# PLUMED interprets the command:
# fcc: FCCUBIC SPECIES=1-256 SWITCH={CUBIC D_0=0.4 D_MAX=0.5} MORE_THAN={RATIONAL R_0=0.45 NN=12 MM=24}
# as follows (Click the red comment above to revert to the short version of the input):
fcc_grpThe GROUP action with label fcc_grp calculates the following quantities: Quantity    Type    Description  
fcc_grp atoms indices of atoms specified in GROUP: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=1-256
fcc_matThe CONTACT_MATRIX action with label fcc_mat calculates the following quantities: Quantity    Type    Description  
fcc_mat.w matrix a matrix containing the weights for the bonds between each pair of atoms
fcc_mat.x matrix the projection of the bond on the x axis
fcc_mat.y matrix the projection of the bond on the y axis
fcc_mat.z matrix the projection of the bond on the z axis: CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More details GROUPspecifies the list of atoms that should be assumed indistinguishable=1-256 SWITCHthe input for the switching function that acts upon the distance between each pair of atoms. Options for this keyword are explained in the documentation for LESS_THAN.={CUBIC D_0=0.4 D_MAX=0.5} COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix
fcc_vfuncThe FCCUBIC_FUNC action with label fcc_vfunc calculates the following quantities: Quantity    Type    Description  
fcc_vfunc matrix the matrix obtained by doing an element-wise application of a function that measures the similarity with an fcc environment to the input matrix: FCCUBIC_FUNCMeasure how similar the environment around atoms is to that found in a FCC structure. More details ARGthe values input to this function=fcc_mat.x,fcc_mat.y,fcc_mat.z ALPHA The alpha parameter of the angular function=3.0
fcc_wvfuncThe CUSTOM action with label fcc_wvfunc calculates the following quantities: Quantity    Type    Description  
fcc_wvfunc matrix the matrix obtained by doing an element-wise application of an arbitrary function to the input matrix: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=fcc_vfunc,fcc_mat.w FUNCthe function you wish to evaluate=x*y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
fcc_onesThe CONSTANT action with label fcc_ones calculates the following quantities: Quantity    Type    Description  
fcc_ones vector the constant value that was read from the plumed input: ONESCreate a constant vector with all elements equal to one More details SIZEthe number of ones that you would like to create=256
fccThe MATRIX_VECTOR_PRODUCT action with label fcc calculates the following quantities: Quantity    Type    Description  
fcc vector the vector that is obtained by taking the product between the matrix and the vector that were input: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=fcc_wvfunc,fcc_ones
fcc_denomThe MATRIX_VECTOR_PRODUCT action with label fcc_denom calculates the following quantities: Quantity    Type    Description  
fcc_denom vector the vector that is obtained by taking the product between the matrix and the vector that were input: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=fcc_mat.w,fcc_ones
fcc_nThe CUSTOM action with label fcc_n calculates the following quantities: Quantity    Type    Description  
fcc_n vector the vector obtained by doing an element-wise application of an arbitrary function to the input vectors: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=fcc,fcc_denom FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
fcc_mtThe MORE_THAN action with label fcc_mt calculates the following quantities: Quantity    Type    Description  
fcc_mt vector the vector obtained by doing an element-wise application of a function that is one if the if the input is more than a threshold to the input vectors: MORE_THANUse a switching function to determine how many of the input variables are more than a certain cutoff. More details ARGthe values input to this function=fcc_n SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above={RATIONAL R_0=0.45 NN=12 MM=24}
fcc_morethanThe SUM action with label fcc_morethan calculates the following quantities: Quantity    Type    Description  
fcc_morethan scalar the SUM of the elements in the input value: SUMCalculate the sum of the arguments More details ARGthe vector/matrix/grid whose elements shuld be added together=fcc_mt PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# --- End of included input --- 
# Basis functions
bf1: BF_LEGENDRELegendre polynomials basis functions. More details ORDERThe order of the basis function expansion=8 MINIMUMThe minimum of the interval on which the basis functions are defined=-26500 MAXIMUMThe maximum of the interval on which the basis functions are defined=-23500
The BF_LEGENDRE action with label bf1 calculates somethingbf2: BF_LEGENDRELegendre polynomials basis functions. More details ORDERThe order of the basis function expansion=8 MINIMUMThe minimum of the interval on which the basis functions are defined=8.0 MAXIMUMThe maximum of the interval on which the basis functions are defined=11.5
The BF_LEGENDRE action with label bf2 calculates somethingbf3: BF_LEGENDRELegendre polynomials basis functions. More details ORDERThe order of the basis function expansion=8 MINIMUMThe minimum of the interval on which the basis functions are defined=0.0 MAXIMUMThe maximum of the interval on which the basis functions are defined=256.0

# Target distribution
The BF_LEGENDRE action with label bf3 calculates somethingTD_MULTITHERMAL_MULTIBARICMultithermal-multibaric target distribution (dynamic). More details ...
 LABELa label for the action so that its output can be referenced in the input to other actions=td_multitp
 MIN_TEMPMinimum energy=350.0
 MAX_TEMPMaximum energy=450.0
 MIN_PRESSUREMinimum pressure=0.06022140857
 MAX_PRESSUREMaximum pressure=602.2140857
 PRESSURETarget pressure of the barostat used in the MD engine=301.10704285
 SIGMAThe standard deviation parameters of the Gaussian kernels used for smoothing the target distribution=250.0,0.1,10.0
 THRESHOLD Maximum exploration free energy in kT=15
 STEPS_TEMP Number of temperature steps=20
 STEPS_PRESSURE Number of pressure steps=20
... TD_MULTITHERMAL_MULTIBARIC

# Expansion
The TD_MULTITHERMAL_MULTIBARIC action with label td_multitp calculates somethingVES_LINEAR_EXPANSIONLinear basis set expansion bias. This action has hidden defaults. More details ...
 ARGthe labels of the scalars on which the bias will act=energy,vol,fcc.morethan
 BASIS_FUNCTIONSthe label of the one dimensional basis functions that should be used=bf1,bf2,bf3
 TEMPthe system temperature - this is needed if the MD code does not pass the temperature to PLUMED=400.0
 GRID_BINSthe number of bins used for the grid=40,40,40
 TARGET_DISTRIBUTIONthe label of the target distribution to be used=td_multitp
 LABELa label for the action so that its output can be referenced in the input to other actions=b1The VES_LINEAR_EXPANSION action with label b1 calculates the following quantities: Quantity    Type    Description  
b1.bias scalar the instantaneous value of the bias potential
b1.force2 scalar the instantaneous value of the squared force due to this bias potential.
... VES_LINEAR_EXPANSION
VES_LINEAR_EXPANSIONLinear basis set expansion bias. This action uses the defaults shown here. More details ...
 ARGthe labels of the scalars on which the bias will act=energy,vol,fcc.morethan
 BASIS_FUNCTIONSthe label of the one dimensional basis functions that should be used=bf1,bf2,bf3
 TEMPthe system temperature - this is needed if the MD code does not pass the temperature to PLUMED=400.0
 GRID_BINSthe number of bins used for the grid=40,40,40
 TARGET_DISTRIBUTIONthe label of the target distribution to be used=td_multitp
 LABELa label for the action so that its output can be referenced in the input to other actions=b1
 GRID_FMT the format to use when outputting the numbers in the grids=%14.9f
... VES_LINEAR_EXPANSION

# Optimization algorithm
OPT_AVERAGED_SGDAveraged stochastic gradient decent with fixed step size. This action has hidden defaults. More details ...
  BIASthe label of the VES bias to be optimized=b1
  STRIDEthe frequency of updating the coefficients given in the number of MD steps=500
  LABELa label for the action so that its output can be referenced in the input to other actions=o1
  STEPSIZEthe step size used for the optimization=1.0
  FES_OUTPUThow often the FES(s) should be written out to file=500
  BIAS_OUTPUThow often the bias(es) should be written out to file=500
  TARGETDIST_OUTPUThow often the dynamic target distribution(s) should be written out to file=500
  COEFFS_OUTPUT how often the coefficients should be written to file=100
  TARGETDIST_STRIDEstride for updating a target distribution that is iteratively updated during the optimization=500
... OPT_AVERAGED_SGD
The OPT_AVERAGED_SGD action with label o1 calculates the following quantities: Quantity    Description  
o1.value a scalarOPT_AVERAGED_SGDAveraged stochastic gradient decent with fixed step size. This action uses the defaults shown here. More details ...
  BIASthe label of the VES bias to be optimized=b1
  STRIDEthe frequency of updating the coefficients given in the number of MD steps=500
  LABELa label for the action so that its output can be referenced in the input to other actions=o1
  STEPSIZEthe step size used for the optimization=1.0
  FES_OUTPUThow often the FES(s) should be written out to file=500
  BIAS_OUTPUThow often the bias(es) should be written out to file=500
  TARGETDIST_OUTPUThow often the dynamic target distribution(s) should be written out to file=500
  COEFFS_OUTPUT how often the coefficients should be written to file=100
  TARGETDIST_STRIDEstride for updating a target distribution that is iteratively updated during the optimization=500
 COEFFS_FILE the name of output file for the coefficients=coeffs.data
... OPT_AVERAGED_SGD

Full list of keywords

The following table describes the keywords and options that can be used with this action

Keyword	Type	Default	Description
THRESHOLD	compulsory	5	Maximum exploration free energy in kT
EPSILONThis keyword do not have examples	compulsory	10	The zeros of the target distribution are changed to e^-EPSILON
MIN_TEMP	compulsory	none	Minimum energy
MAX_TEMP	compulsory	none	Maximum energy
MIN_PRESSURE	compulsory	none	Minimum pressure
MAX_PRESSURE	compulsory	none	Maximum pressure
PRESSURE	compulsory	none	Target pressure of the barostat used in the MD engine
STEPS_TEMP	compulsory	20	Number of temperature steps
STEPS_PRESSURE	compulsory	20	Number of pressure steps
SIGMA	optional	not used	The standard deviation parameters of the Gaussian kernels used for smoothing the target distribution

References

More information about how this action can be used is available in the following articles:

Quantity	Type	Description
energy	scalar	the internal energy