TD_MULTITHERMAL_MULTIBARIC

	This is part of the ves module
	It is only available if you configure PLUMED with ./configure –enable-modules=ves . Furthermore, this feature is still being developed so take care when using it and report any problems on the mailing list.

Multithermal-multibaric target distribution (dynamic).

Use the target distribution to sample the multithermal-multibaric ensemble [86] [81]. 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 [87] .

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 [111] or in the multicanonical ensemble [13] , 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 [109].

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
energy: ENERGY The ENERGY action with label energy calculates a single scalar value
vol: VOLUME The VOLUME action with label vol calculates a single scalar value
# Basis functions
bf1: BF_LEGENDRE 
ORDER
compulsory keyword 
The order of the basis function expansion. 
=10 
MINIMUM
compulsory keyword 
The minimum of the interval on which the basis functions are defined. 
=-14750 
MAXIMUM
compulsory keyword 
The maximum of the interval on which the basis functions are defined. 
=-12250 The BF_LEGENDRE action with label bf1
bf2: BF_LEGENDRE 
ORDER
compulsory keyword 
The order of the basis function expansion. 
=10 
MINIMUM
compulsory keyword 
The minimum of the interval on which the basis functions are defined. 
=6.5 
MAXIMUM
compulsory keyword 
The maximum of the interval on which the basis functions are defined. 
=8.25 The BF_LEGENDRE action with label bf2
# Target distribution - 1 bar = 0.06022140857 and 300 MPa = 180.66422571
td_multi: TD_MULTITHERMAL_MULTIBARIC ...
   
MIN_TEMP
compulsory keyword 
Minimum energy. 
=260 
   
MAX_TEMP
compulsory keyword 
Maximum energy. 
=350 
   
MAX_PRESSURE
compulsory keyword 
Maximum pressure. 
=180.66422571 
   
MIN_PRESSURE
compulsory keyword 
Minimum pressure. 
=0.06022140857 
   
PRESSURE
compulsory keyword 
Target pressure of the barostat used in the MD engine. 
=0.06022140857 
   
...The TD_MULTITHERMAL_MULTIBARIC action with label td_multi
# Bias expansion
b1: VES_LINEAR_EXPANSION ...
   
ARG
the input for this action is the scalar output from one or more other actions. 
=energy,vol 
   
BASIS_FUNCTIONS
compulsory keyword 
the label of the one dimensional basis functions that should be used. 
=bf1,bf2 
   
TEMP
the system temperature - this is needed if the MD code does not pass the temperature
to PLUMED. 
=300.0 
   
GRID_BINS
the number of bins used for the grid. 
=200,200 
   
TARGET_DISTRIBUTION
the label of the target distribution to be used. 
=td_multi 
   
...The VES_LINEAR_EXPANSION action with label b1 calculates the following quantities:

 Quantity   
 Description  
b1.bias
the instantaneous value of the bias potential
b1.force2
the instantaneous value of the squared force due to this bias potential.


# Optimization algorithm
o1: OPT_AVERAGED_SGD ...
   
BIAS
compulsory keyword 
the label of the VES bias to be optimized 
=b1 
   
STRIDE
compulsory keyword 
the frequency of updating the coefficients given in the number of MD steps. 
=500 
   
   
STEPSIZE
compulsory keyword 
the step size used for the optimization 
=1.0 
   
FES_OUTPUT
how often the FES(s) should be written out to file. 
=500 
   
BIAS_OUTPUT
how often the bias(es) should be written out to file. 
=500 
   
TARGETDIST_OUTPUT
how often the dynamic target distribution(s) should be written out to file. 
=500 
   
COEFFS_OUTPUT
compulsory keyword ( default=100 )
how often the coefficients should be written to file. 
=100 
   
TARGETDIST_STRIDE
stride for updating a target distribution that is iteratively updated during the
optimization. 
=100 
   The OPT_AVERAGED_SGD action with label o1

... 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

contract shortcuts

# Use energy, volume and FCCUBIC as CVs
energy: ENERGY The ENERGY action with label energy calculates a single scalar value
vol: VOLUME The VOLUME action with label vol calculates a single scalar value
fcc: FCCUBIC 
SPECIES
this keyword is used for colvars such as coordination number. 
=1-256 
SWITCH
the switching function that it used in the construction of the contact matrix 
={CUBIC D_0=0.4 D_MAX=0.5}  
MORE_THAN
calculate the number of variables that are more than a certain target value. 
={RATIONAL R_0=0.45 NN=12 MM=24}  The FCCUBIC action with label fcc calculates a single scalar value
# Basis functions
bf1: BF_LEGENDRE 
ORDER
compulsory keyword 
The order of the basis function expansion. 
=8 
MINIMUM
compulsory keyword 
The minimum of the interval on which the basis functions are defined. 
=-26500 
MAXIMUM
compulsory keyword 
The maximum of the interval on which the basis functions are defined. 
=-23500 The BF_LEGENDRE action with label bf1
bf2: BF_LEGENDRE 
ORDER
compulsory keyword 
The order of the basis function expansion. 
=8 
MINIMUM
compulsory keyword 
The minimum of the interval on which the basis functions are defined. 
=8.0 
MAXIMUM
compulsory keyword 
The maximum of the interval on which the basis functions are defined. 
=11.5 The BF_LEGENDRE action with label bf2
bf3: BF_LEGENDRE 
ORDER
compulsory keyword 
The order of the basis function expansion. 
=8 
MINIMUM
compulsory keyword 
The minimum of the interval on which the basis functions are defined. 
=0.0 
MAXIMUM
compulsory keyword 
The maximum of the interval on which the basis functions are defined. 
=256.0 The BF_LEGENDRE action with label bf3
# Target distribution
td_multitp: TD_MULTITHERMAL_MULTIBARIC ...
   
   
MIN_TEMP
compulsory keyword 
Minimum energy. 
=350.0 
   
MAX_TEMP
compulsory keyword 
Maximum energy. 
=450.0 
   
MIN_PRESSURE
compulsory keyword 
Minimum pressure. 
=0.06022140857 
   
MAX_PRESSURE
compulsory keyword 
Maximum pressure. 
=602.2140857 
   
PRESSURE
compulsory keyword 
Target pressure of the barostat used in the MD engine. 
=301.10704285 
   
SIGMA
The standard deviation parameters of the Gaussian kernels used for smoothing the
target distribution. 
=250.0,0.1,10.0 
   
THRESHOLD
compulsory keyword ( default=5 )
Maximum exploration free energy in kT. 
=15 
   
STEPS_TEMP
compulsory keyword ( default=20 )
Number of temperature steps. 
=20 
   
STEPS_PRESSURE
compulsory keyword ( default=20 )
Number of pressure steps. 
=20 
...The TD_MULTITHERMAL_MULTIBARIC action with label td_multitp
# Expansion
b1: VES_LINEAR_EXPANSION ...
   
ARG
the input for this action is the scalar output from one or more other actions. 
=energy,vol,fcc.morethan 
   
BASIS_FUNCTIONS
compulsory keyword 
the label of the one dimensional basis functions that should be used. 
=bf1,bf2,bf3 
   
TEMP
the system temperature - this is needed if the MD code does not pass the temperature
to PLUMED. 
=400.0 
   
GRID_BINS
the number of bins used for the grid. 
=40,40,40 
   
TARGET_DISTRIBUTION
the label of the target distribution to be used. 
=td_multitp 
   
...The VES_LINEAR_EXPANSION action with label b1 calculates the following quantities:

 Quantity   
 Description  
b1.bias
the instantaneous value of the bias potential
b1.force2
the instantaneous value of the squared force due to this bias potential.


# Optimization algorithm
o1: OPT_AVERAGED_SGD ...
   
BIAS
compulsory keyword 
the label of the VES bias to be optimized 
=b1 
   
STRIDE
compulsory keyword 
the frequency of updating the coefficients given in the number of MD steps. 
=500 
   
   
STEPSIZE
compulsory keyword 
the step size used for the optimization 
=1.0 
   
FES_OUTPUT
how often the FES(s) should be written out to file. 
=500 
   
BIAS_OUTPUT
how often the bias(es) should be written out to file. 
=500 
   
TARGETDIST_OUTPUT
how often the dynamic target distribution(s) should be written out to file. 
=500 
   
COEFFS_OUTPUT
compulsory keyword ( default=100 )
how often the coefficients should be written to file. 
=100 
   
TARGETDIST_STRIDE
stride for updating a target distribution that is iteratively updated during the
optimization. 
=500 
   The OPT_AVERAGED_SGD action with label o1

# Use energy, volume and FCCUBIC as CVs
energy: ENERGY The ENERGY action with label energy calculates a single scalar value
vol: VOLUME The VOLUME action with label vol calculates a single scalar value
fcc_grp: GROUP 
ATOMS
the numerical indexes for the set of atoms in the group. 
=1-256 The GROUP action with label fcc_grp defines a group of atoms so that they can be referred to later in the input
fcc_mat: CONTACT_MATRIX 
GROUP
specifies the list of atoms that should be assumed indistinguishable. 
=1-256 
SWITCH
specify the switching function to use between two sets of indistinguishable atoms.
=CUBIC,D_0=0.4,D_MAX=0.5 
COMPONENTS
( default=off ) also calculate the components of the vector connecting the atoms
in the contact matrix 
 The CONTACT_MATRIX action with label fcc_mat... OPT_AVERAGED_SGD
Glossary of keywords and components

Compulsory keywords



 THRESHOLD  
( default=5 ) Maximum exploration free energy in kT.   

 EPSILON  
( default=10 ) The zeros of the target distribution are changed to e^-EPSILON.   

 MIN_TEMP  
Minimum energy.   

 MAX_TEMP  
Maximum energy.   

 MIN_PRESSURE  
Minimum pressure.   

 MAX_PRESSURE  
Maximum pressure.   

 PRESSURE  
Target pressure of the barostat used in the MD engine.   

 STEPS_TEMP  
( default=20 ) Number of temperature steps.   

 STEPS_PRESSURE  
( default=20 ) Number of pressure steps.   

Options



 SIGMA  
The standard deviation parameters of the Gaussian kernels used for smoothing the target distribution. One value must be specified for each argument, i.e. one value per CV. A value of 0.0 means that no smoothing is performed, this is the default behavior.