 |
 |
 |
| |
|
|
| 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:
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 1. 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.
Much like in the Wang-Landau algorithm [113] or in the multicanonical ensemble [11] , 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 [111].
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.
| THRESHOLD | ( default=1 ) Maximum exploration free energy in kT. |
| 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. |
| 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. |
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.
# Use energy and volume as CVs energy: ENERGY vol: VOLUME # Basis functions bf1: BF_LEGENDRE ORDER=10 MINIMUM=-14750 MAXIMUM=-12250 bf2: BF_LEGENDRE ORDER=10 MINIMUM=6.5 MAXIMUM=8.25 # Target distribution TD_MULTITHERMAL_MULTIBARIC ... MIN_TEMP=260 MAX_TEMP=350 MAX_PRESSURE=180.66422571 # 300 MPa MIN_PRESSURE=0.06022140857 # 1 bar PRESSURE=0.06022140857 # 1 bar STEPS_PRESSURE=20 STEPS_TEMP=20 SIGMA=50.,0.05 THRESHOLD=1 LABEL=td_multi ... TD_MULTITHERMAL_MULTIBARIC # Bias expansion VES_LINEAR_EXPANSION ... ARG=energy,vol BASIS_FUNCTIONS=bf1,bf2 TEMP=300.0 GRID_BINS=200,200 TARGET_DISTRIBUTION=td_multi LABEL=b1 ... VES_LINEAR_EXPANSION # Optimization algorithm OPT_AVERAGED_SGD ... BIAS=b1 STRIDE=500 LABEL=o1 STEPSIZE=1.0 FES_OUTPUT=500 BIAS_OUTPUT=500 TARGETDIST_OUTPUT=500 COEFFS_OUTPUT=100 TARGETDIST_STRIDE=100 ... 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.
# Use energy, volume and FCCUBIC as CVs energy: ENERGY vol: VOLUME 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} # Basis functions bf1: BF_LEGENDRE ORDER=8 MINIMUM=-26500 MAXIMUM=-23500 bf2: BF_LEGENDRE ORDER=8 MINIMUM=8.0 MAXIMUM=11.5 bf3: BF_LEGENDRE ORDER=8 MINIMUM=0.0 MAXIMUM=256.0 # Target distribution TD_MULTITHERMAL_MULTIBARIC ... LABEL=td_multitp MIN_TEMP=350.0 MAX_TEMP=450.0 MIN_PRESSURE=0.06022140857 MAX_PRESSURE=602.2140857 PRESSURE=301.10704285 SIGMA=250.0,0.1,10.0 THRESHOLD=15 STEPS_TEMP=20 STEPS_PRESSURE=20 ... TD_MULTITHERMAL_MULTIBARIC # Expansion VES_LINEAR_EXPANSION ... ARG=energy,vol,fcc.morethan BASIS_FUNCTIONS=bf1,bf2,bf3 TEMP=400.0 GRID_BINS=40,40,40 TARGET_DISTRIBUTION=td_multitp LABEL=b1 ... VES_LINEAR_EXPANSION # Optimization algorithm OPT_AVERAGED_SGD ... BIAS=b1 STRIDE=500 LABEL=o1 STEPSIZE=1.0 FES_OUTPUT=500 BIAS_OUTPUT=500 TARGETDIST_OUTPUT=500 COEFFS_OUTPUT=100 TARGETDIST_STRIDE=500 ... OPT_AVERAGED_SGD
1.8.14