This Masterclass aims to introduce users to the tools available via PLUMED to analyze and perform simulations of concentration-driven processes such as nucleation, growth, and diffusion.
Once this Masterclass is completed, users will be able to:
We assume that you are familiar with PLUMED. However, the 2021 PLUMED Masterclass is a great place to start if you are not.
Follow the instructions provided here to install gromacs+plumed2.8 compiled with CmuMD.
PLUMED facilitates the calculation of functions of atom coordinates and the introduction of bias potentials in atomistic simulations.
In this tutorial, we review how to use PLUMED functionalities for two complementary tasks:
The data needed to complete the exercises of this Masterclass can be found on GitHub. You can clone this repository locally on your machine using the following command:
git clone https://github.com/mme-ucl/PLUMED_MClass_22-08
The solution of this masterclass is available on GitHub as a jupyter notebook.
We can start by using plumed to analyze a trajectory and identify the total number of atoms belonging to a liquid phase condensing from a supersaturated vapour in a long unbiased MD simulation.
To this aim, a typical choice is the implementation of a criterion inspired by the work of Ten Wolde and Frenkel: all the atoms with a threshold number of nearest neighbours larger than a threshold is considered as part of the liquid phase.
In PLUMED this criterion can be readily implemented using the multicolvar COORDINATIONNUMBER:
# LIQUID-like atoms lq: COORDINATIONNUMBER SPECIESthis keyword is used for colvars such as coordination number. =1-10000 SWITCHThis keyword is used if you want to employ an alternative to the continuous switching function defined above. ={CUBIC D_0=0.45 D_MAX=0.55} MORE_THANcalculate the number of variables more than a certain target value. ={RATIONAL R_0=5.0 D_MAX=10.0} PRINT ARGthe input for this action is the scalar output from one or more other actions. =lq.morethan STRIDEcompulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output =1 FILEthe name of the file on which to output these quantities =nliquid.dat
The time evolution of the number of atoms in the liquid phase can be computed by analyzing the trajectory using the driver functionality (where liquid.dat is the plumed file for this analysis).
plumed driver --mf_xtc LJvap_to_liq_red.xtc --plumed liquid.dat
Using other MultiColvar keywords can you count the number of atoms on the surface of the clusters forming in this trajectory? Can you estimate how the surface to volume ratio of the droplets forming in the system changes directly within PLUMED?
The MD trajectory that we will analyze can be found in the folder called data
:
LJvap.gro
: reference conformation of 10000 LJ particles in a metastable vapour phase.LJvap_to_liq.xtc
: trajectory.In certain contexts, it can be helpful to characterize the phase transition by following analyzing the population of clusters that form during the nucleation process. In PLUMED, this is possible by taking advantage of the DFSCLUSTERING algorithm, which builds on the construction of a CONTACTMATRIX to identify clusters of atoms with specific properties.
Here we analyse the clusters distribution emerging in the trajectory LJvap_to_liq.xtc
, focussing on the number of clusters, and following the evolution of the four largest clusters in the system:
# Identify liquid-like atoms lq: COORDINATIONNUMBER SPECIESthis keyword is used for colvars such as coordination number. =1-10000 SWITCHThis keyword is used if you want to employ an alternative to the continuous switching function defined above. ={CUBIC D_0=0.45 D_MAX=0.55} MORE_THANcalculate the number of variables more than a certain target value. ={RATIONAL R_0=5.0 D_MAX=10.0} # Define a contact matrix & perfom DFS clustering cm: CONTACT_MATRIX ATOMSThe list of atoms for which you would like to calculate the contact matrix. =lq SWITCHThis keyword is used if you want to employ an alternative to the continuous switching function defined above. ={CUBIC D_0=0.45 D_MAX=0.55} dfs: DFSCLUSTERING MATRIXcompulsory keyword the action that calculates the adjacency matrix vessel we would like to analyze =cm # Compute the size of the four largest clusters cluster_1: CLUSTER_NATOMS CLUSTERScompulsory keyword the label of the action that does the clustering =dfs CLUSTERcompulsory keyword ( default=1 ) which cluster would you like to look at 1 is the largest cluster, 2 is the second largest, 3 is the the third largest and so on. =1 cluster_2: CLUSTER_NATOMS CLUSTERScompulsory keyword the label of the action that does the clustering =dfs CLUSTERcompulsory keyword ( default=1 ) which cluster would you like to look at 1 is the largest cluster, 2 is the second largest, 3 is the the third largest and so on. =2 cluster_3: CLUSTER_NATOMS CLUSTERScompulsory keyword the label of the action that does the clustering =dfs CLUSTERcompulsory keyword ( default=1 ) which cluster would you like to look at 1 is the largest cluster, 2 is the second largest, 3 is the the third largest and so on. =3 cluster_4: CLUSTER_NATOMS CLUSTERScompulsory keyword the label of the action that does the clustering =dfs CLUSTERcompulsory keyword ( default=1 ) which cluster would you like to look at 1 is the largest cluster, 2 is the second largest, 3 is the the third largest and so on. =4 # Compute the number of clusters nclust: CLUSTER_DISTRIBUTION CLUSTERScompulsory keyword the label of the action that does the clustering =dfs MORE_THANcalculate the number of variables more than a certain target value. ={GAUSSIAN D_0=D_0=1.95 R_0=0.01 D_MAX=1.99} # PRINT to file PRINT ARGthe input for this action is the scalar output from one or more other actions. =lq.morethan,nclust.*,cluster_1,cluster_2,cluster_3,cluster_4 STRIDEcompulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output =1 FILEthe name of the file on which to output these quantities =clusters.dat FLUSH STRIDEcompulsory keyword the frequency with which all the open files should be flushed =1
Modifying the setup above, can you compute the number of clusters with a size larger than 20 LJ particles?
The MD trajectory that we will analyze can be found in the folder called data
:
LJvap.gro
: reference conformation of 10000 LJ particles in a metastable vapour phase.LJvap_to_liq.xtc
: trajectory.In the two examples discussed in Exercises 1 and 2 it is apparent that the process of condensation reaches steady state due to finite-size effects. However, looking at the process more carefully, one can notice that the driving force leading to phase separation is far from constant. C$\mu$MD allows using PLUMED to apply ad-hoc forces and keep constant the composition of spatial regions of the simulation box called control regions. This feature allows to model processes driven by concentration in steady, out-of-equilibrium conditions.
The first example we will focus on is a purely diffusive process in a LJ vapour.
A PLUMED file that allows imposing a steady concentration difference with CmuMD looks like:
# Define groups of atoms LJ: GROUP ATOMSthe numerical indexes for the set of atoms in the group. =1-1000:1 # Provide parameters for the CV left: CMUMD GROUP could not find this keyword =lj NSV could not find this keyword =1 FIXED could not find this keyword =0.5 DCR could not find this keyword =0.25 CRSIZE could not find this keyword =0.1 WF could not find this keyword =0.0001 ASYMM could not find this keyword =-1 NINT could not find this keyword =0.1 NZ could not find this keyword =291 right: CMUMD GROUP could not find this keyword =lj NSV could not find this keyword =1 FIXED could not find this keyword =0.5 DCR could not find this keyword =0.25 CRSIZE could not find this keyword =0.1 WF could not find this keyword =0.0001 ASYMM could not find this keyword =1 NINT could not find this keyword =0.1 NZ could not find this keyword =291 # CmuMD is implemented as a restraint on the densities of species in CR rleft: RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =left ATcompulsory keyword the position of the restraint =1.4 KAPPAcompulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are =1000.0 rright: RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =right ATcompulsory keyword the position of the restraint =0.05 KAPPAcompulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are =2000.0 # Report the densities and bias PRINT ... ARGthe input for this action is the scalar output from one or more other actions. =left,right,rleft.bias,rright.bias STRIDEcompulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output =10 FILEthe name of the file on which to output these quantities =CMUMD_log ...
We can run a CMUMD simulation under these conditions as:
gmx_mpi mdrun --plumed cmumd_diff.dat
Perform simulations at varying concentration differences. Compute the concentration gradient across the simulation box as a function of the concentration difference between left and right controlled volumes.
In the folder data
:
LJ_diffusion.gro
: reference conformation of 1000 LJ particles.LJ_diffusion.xtc
: trajectory evolving under the effect of the stationary concentration gradient.The fourth exercise combines aspects of the previous three. We will use CmuMD to control the driving force associated with the growth of a dense phase represented by a slab at the center of a simulation box.
Similarly to the diffusion case the plumed file that can be used to perform this type of simulation reads:
# Define groups of atoms LJ: GROUP ATOMSthe numerical indexes for the set of atoms in the group. =1001-2000:1 # Provide parameters for the CV left: CMUMD GROUP could not find this keyword =lj NSV could not find this keyword =1 FIXED could not find this keyword =0.4 DCR could not find this keyword =0.25 CRSIZE could not find this keyword =0.1 WF could not find this keyword =0.0001 ASYMM could not find this keyword =-1 NINT could not find this keyword =0.1 NZ could not find this keyword =291 right: CMUMD GROUP could not find this keyword =lj NSV could not find this keyword =1 FIXED could not find this keyword =0.6 DCR could not find this keyword =0.25 CRSIZE could not find this keyword =0.1 WF could not find this keyword =0.0001 ASYMM could not find this keyword =1 NINT could not find this keyword =0.1 NZ could not find this keyword =291 # CmuMD is implemented as a restraint on the densities of species in CR left: RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =left ATcompulsory keyword the position of the restraint =3. KAPPAcompulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are =2000.0 right: RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =right ATcompulsory keyword the position of the restraint =3. KAPPAcompulsory keyword ( default=0.0 ) specifies that the restraint is harmonic and what the values of the force constants on each of the variables are =2000.0 # Report the densities and bias PRINT ... ARGthe input for this action is the scalar output from one or more other actions. =left,right,rleft.bias,rright.bias STRIDEcompulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output =10 FILEthe name of the file on which to output these quantities =CMUMD_log ...
In the folder data
, a CmuMD trajectory ready for analysis can be found together with the gromacs input files necessary to setup and run it.
LJ_slab.gro
: reference conformation of 1000 LJ particles.LJ_slab.xtc
: trajectory evolving under the effect of the stationary concentration gradient.md_input_LJ_slab
: MD input files