PLUMED Masterclass 22.8: Modelling Concentration-driven processes with PLUMED
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
tested on master
# LIQUID-like atoms
lq: COORDINATIONNUMBER 
SPECIES
this keyword is used for colvars such as coordination number.
=1-10000
SWITCH
This 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_THAN
calculate the number of variables more than a certain target value.
={RATIONAL R_0=5.0 D_MAX=10.0} PRINT
ARG
the input for this action is the scalar output from one or more other actions.
=lq.morethan
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=1
FILE
the 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

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
tested on master
# Identify liquid-like atoms
lq: COORDINATIONNUMBER 
SPECIES
this keyword is used for colvars such as coordination number.
=1-10000
SWITCH
This 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_THAN
calculate 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
ATOMS
The list of atoms for which you would like to calculate the contact matrix.
=lq
SWITCH
This 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
MATRIX
compulsory 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
CLUSTERS
compulsory keyword the label of the action that does the clustering
=dfs
CLUSTER
compulsory 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
CLUSTERS
compulsory keyword the label of the action that does the clustering
=dfs
CLUSTER
compulsory 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
CLUSTERS
compulsory keyword the label of the action that does the clustering
=dfs
CLUSTER
compulsory 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
CLUSTERS
compulsory keyword the label of the action that does the clustering
=dfs
CLUSTER
compulsory 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
CLUSTERS
compulsory keyword the label of the action that does the clustering
=dfs
MORE_THAN
calculate 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
ARG
the 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
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=1
FILE
the name of the file on which to output these quantities
=clusters.dat FLUSH
STRIDE
compulsory keyword the frequency with which all the open files should be flushed
=1

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
tested on master
# Define groups of atoms
LJ: GROUP 
ATOMS
the 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
ARG
the input for this action is the scalar output from one or more other actions.
=left
AT
compulsory keyword the position of the restraint
=1.4
KAPPA
compulsory 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
ARG
the input for this action is the scalar output from one or more other actions.
=right
AT
compulsory keyword the position of the restraint
=0.05
KAPPA
compulsory 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 ...
ARG
the input for this action is the scalar output from one or more other actions.
=left,right,rleft.bias,rright.bias
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=10
FILE
the 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

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
tested on master
# Define groups of atoms
LJ: GROUP 
ATOMS
the 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
ARG
the input for this action is the scalar output from one or more other actions.
=left
AT
compulsory keyword the position of the restraint
=3.
KAPPA
compulsory 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
ARG
the input for this action is the scalar output from one or more other actions.
=right
AT
compulsory keyword the position of the restraint
=3.
KAPPA
compulsory 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 ...
ARG
the input for this action is the scalar output from one or more other actions.
=left,right,rleft.bias,rright.bias
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=10
FILE
the name of the file on which to output these quantities
=CMUMD_log ...

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.