Author: Matteo Salvalaglio

Date: May 20, 2022

Aims

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.

Objectives

Once this Masterclass is completed, users will be able to:

Write a PLUMED input file to analyze a trajectory with multicolvars.
Use the Depth-First-Search tools available in PLUMED to characterize clusters of atoms in an existing vapour condensation trajectory.
Write a PLUMED input file to analyze the condensation of a simple liquid phase from vapour.
Write a PLUMED input file for CmuMD, and perform a CmuMD simulation of liquid condensation at constant driving force.
Write a PLUMED input file for CmuMD, and perform a CmuMD simulation to model steady-state diffusive flux across a simulation box.

Prerequisites

We assume that you are familiar with PLUMED. However, the 2021 PLUMED Masterclass is a great place to start if you are not.

Setting up the software

Follow the instructions provided here to install gromacs+plumed2.8 compiled with CmuMD.

Background Overview

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:

Calculating collective variables that capture processes of assembly/disassembly such as the nucleation, condensation, dissolution etc. (Ex. 1, 2)
Introducing biasing forces through CmuMD, a method that mimics open-boundary conditions and allows for mitigation of finite-size effects. (Ex. 3, 4)

Resources

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

Solution

The solution of this masterclass is available on GitHub as a jupyter notebook.

Exercise 1: analysis of a nucleation trajectory I

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:

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

# 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}   You cannot view the components that are calculated by each action for this input file. Sorry 
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  You cannot view the components that are calculated by each action for this input file. Sorry

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

Try on your own

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?

Available data

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.

Exercise 2: analysis of a nucleation trajectory II

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:

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

# 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} You cannot view the components that are calculated by each action for this input file. Sorry
# 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} You cannot view the components that are calculated by each action for this input file. Sorry
dfs: DFSCLUSTERING MATRIXcompulsory keyword
the action that calculates the adjacency matrix vessel we would like to analyze
=cm You cannot view the components that are calculated by each action for this input file. Sorry
# 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 You cannot view the components that are calculated by each action for this input file. Sorry
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 You cannot view the components that are calculated by each action for this input file. Sorry
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 You cannot view the components that are calculated by each action for this input file. Sorry
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 You cannot view the components that are calculated by each action for this input file. Sorry
# 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} You cannot view the components that are calculated by each action for this input file. Sorry
# 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 You cannot view the components that are calculated by each action for this input file. Sorry
FLUSH STRIDEcompulsory keyword
the frequency with which all the open files should be flushed
=1 You cannot view the components that are calculated by each action for this input file. Sorry

Try on your own

Modifying the setup above, can you compute the number of clusters with a size larger than 20 LJ particles?

Available data

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.

Exercise 3: Steady-state diffusive flux with CmuMD

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:

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

# Define groups of atoms
LJ: GROUP ATOMSthe numerical indexes for the set of atoms in the group. 
=1-1000:1  You cannot view the components that are calculated by each action for this input file. Sorry 
# 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  You cannot view the components that are calculated by each action for this input file. Sorry 
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  You cannot view the components that are calculated by each action for this input file. Sorry 
# 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  You cannot view the components that are calculated by each action for this input file. Sorry 
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  You cannot view the components that are calculated by each action for this input file. Sorry 
# 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 
... You cannot view the components that are calculated by each action for this input file. Sorry

We can run a CMUMD simulation under these conditions as:

gmx_mpi mdrun --plumed cmumd_diff.dat

Try on your own

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.

Available data

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.

Exercise 4: Steady-state condensation process with C$\mu$MD

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:

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

# Define groups of atoms
LJ: GROUP ATOMSthe numerical indexes for the set of atoms in the group.
=1001-2000:1 You cannot view the components that are calculated by each action for this input file. Sorry
# 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 You cannot view the components that are calculated by each action for this input file. Sorry
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 You cannot view the components that are calculated by each action for this input file. Sorry
# 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 You cannot view the components that are calculated by each action for this input file. Sorry
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 You cannot view the components that are calculated by each action for this input file. Sorry
# 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
... You cannot view the components that are calculated by each action for this input file. Sorry

Available data

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

Try on your own

Setup and run CMUMD simulations for the LJ slab system at varying CR concentrations.
Use the multicolvar-based approach discussed in Ex1 to follow the condensation process.
Use the graph-based approach discussed in Ex2 to follow the condensation process.
Monitor the density profile across the simulation box.