View source for Lammps and plumed

This wiki serves as an overview of how to use the PLUMED<ref>The Plumed Consortium. (2019). Promoting transparency and reproducibility in enhanced molecular simulations. Nature Methods, 16(8), 670–673. https://doi.org/10.1038/s41592-019-0506-8</ref> <ref>Tribello, G. A., Bonomi, M., Branduardi, D., Camilloni, C., & Bussi, G. (2014). PLUMED 2: New feathers for an old bird. Computer Physics Communications, 185(2), 604–613. https://doi.org/10.1016/j.cpc.2013.09.018</ref> software package to perform biased molecular dynamics simulations in LAMMPS.

=Using LAMMPS=
The very first point of call for everything involving LAMMPS is the official [https://lammps.sandia.gov/doc/Manual.html LAMMPS manual]. One thing is crucial to know: LAMMPS does '''not''' create/contain force fields! You need to define the topology and the force field, LAMMPS only does the simulation itself. For starters, you can find a few example input files here: [[:File:Simple_example.zip]]. Here I use the classical (and classic) CL&P force field for ionic liquids.<ref>Canongia Lopes, J. N., & Padua, A. A. H. (2012). CL&P: A generic and systematic force field for ionic liquids modeling. Theoretical Chemistry Accounts, 131(3), 1129. https://doi.org/10.1007/s00214-012-1129-7</ref><ref>Canongia Lopes, J. N., Deschamps, J., & Padua, A. A. H. (2004). Modeling Ionic Liquids Using a Systematic All-Atom Force Field. The Journal of Physical Chemistry B, 108(6), 2038–2047. https://doi.org/10.1021/jp0362133</ref><ref>Canongia Lopes, J. N., & Padua, A. A. H. (2004). Molecular Force Field for Ionic Liquids Composed of Triflate or Bistriflylimide Anions. The Journal of Physical Chemistry B, 108(43), 16893–16898. https://doi.org/10.1021/jp0476545</ref> The content of the files and what you should watch out for is described in the following.
; run_2D.pbs : This is the runscript. It runs the [[https://github.com/z-gong/lammps.git ttdamp branch of LAMMPS]], compiled with PLUMED 2.6.0. If you want to reproduce my simulations with Drude particles further down: The required thermostat should now be implemented in the USER-DRUDE package.<ref>Goloviznina, K., Gong, Z., Costa Gomes, M. F., & Pádua, A. A. H. (2021). Extension of the CL&Pol Polarizable Force Field to Electrolytes, Protic Ionic Liquids, and Deep Eutectic Solvents. Journal of Chemical Theory and Computation, 17(3), 1606–1617. https://doi.org/10.1021/acs.jctc.0c01002</ref> Have a look beforehand what packages you need etc, this is very different depending on the application.
; METAD.inp : The LAMMPS input file. Line 37 invokes the fix plumed, with a separate input file named "plumedfile.dat". it runs for 10 ns which is ok, since there is only one ion. It needs the topology given as a data file.
; start.dat : This is the LAMMPS data file, containing information about the simulation box, connectivity, bonded force field parameters, cartesian coordinates and velocities.
; plumedfile.dat : Contains the instructions for PLUMED. This one is an example for 2D well-tempered Metadynamics, the two "collective" variables are the two C-S-N-S dihedral angles in the [NTf2]- anion.
After you run the simulation, you can process the results with the following tools:
; sumhills.pbs : You can use this to sum the gaussians from the Metadynamics simulation to produce a two dimensional free energy surface (fes.dat)
; plothills.py : Shows the hill height (of deposited gaussians) as a function of time, as well as the evolution of one of the collective variables.
; plotcolvar_2D.py : Plots the collective variable values as dots on a 2d map.

If you only want to test LAMMPS without PLUMED, you can change the runscript to:
<pre>
module load gcc
module load intel-suite mpi
module load fftw
module load lammps/11Aug17
cd $PBS_O_WORKDIR
mpiexec lammps -in METAD.inp
</pre>
However, if you do so, you also need to comment out the "fix plumed" in METAD.inp, otherwise LAMMPS terminates with the following error:
<pre>
ERROR: Unknown fix style plumed (../modify.cpp:854)
Last command: fix 1 all plumed plumedfile plumedfile.dat outfile plumed.out
</pre>

=Using PLUMED (2.6.0)=
PLUMED is a powerful software package, and I strongly recommend having a look at the [[http://www.plumed.org/doc-v2.5/user-doc/html/_colvar.html PLUMED documentation]]. There are a lot of tutorials on there, and it is worth having a look at the built-in colvars (and multicolvars). Similarly, there are a few good articles describing how nucleation and phase transitions can be modelled with PLUMED.<ref>Giberti, F., Salvalaglio, M., Mazzotti, M., & Parrinello, M. (2015). Insight into the nucleation of urea crystals from the melt. Chemical Engineering Science, 121, 51–59. https://doi.org/10.1016/j.ces.2014.08.032</ref><ref>Bussi, G., & Tribello, G. A. (2019). Analyzing and Biasing Simulations with PLUMED. In Methods in Molecular Biology (Vol. 2022, pp. 529–578). https://doi.org/10.1007/978-1-4939-9608-7_21</ref><ref>Tribello, G. A., Giberti, F., Sosso, G. C., Salvalaglio, M., & Parrinello, M. (2017). Analyzing and Driving Cluster Formation in Atomistic Simulations. Journal of Chemical Theory and Computation, 13(3), 1317–1327. https://doi.org/10.1021/acs.jctc.6b01073</ref> Here I will only show an example which is typical for me, and you can find a stripped-down version in [[:File:BMIMOTf.zip]].

The system in the example is [C4C1Im][OTf], an ionic liquid with the triflate anion. The force field is the CL&Pol force field developed by Goloviznina ''et al.''.<ref>Goloviznina, K., Canongia Lopes, J. N., Costa Gomes, M., & Pádua, A. A. H. (2019). Transferable, Polarizable Force Field for Ionic Liquids. Journal of Chemical Theory and Computation, 15(11), 5858–5871. https://doi.org/10.1021/acs.jctc.9b00689</ref> This anion, F3C-SO3, has an inherent directionality along its C3 axis. Thus we could hypothesise that there might be a liquid-liquid phase transition to an ordered state where all the OTf anions are aligned.<ref>Anaredy, R. S., & Shaw, S. K. (2018). Developing Distinct Chemical Environments in Ionic Liquid Films. The Journal of Physical Chemistry C, 122(34), 19731–19737. https://doi.org/10.1021/acs.jpcc.8b06608</ref> Now have a look at the '''plumedfile.dat''' in the [[:File:BMIMOTf.zip]] example pack. It has the following sections:
*  A header, defining units and file i/o flush, also [[http://www.plumed.org/doc-v2.5/user-doc/html/_r_e_s_t_a_r_t.html "RESTART"]] if the simulation is to be restarted - important not to overwrite existing files but rather to continue!
*  A section where groups are defined. Here I define all the CF3 groups as their centre of mass as c1,c2,...,c512 (because I have 512 ion pairs in the simulation). Similarly, the centre of mass of the SO3 groups is defined as s1-s512. Note that this simulation contains drude particles, if you wanted to be perfect you would include those in the definition. However, the Drude particles react pretty much instantaneously, so if you need to include them in your definition then your biasing is probably unphysical anyway.
<pre>
c1: COM ATOMS=12801-12804
(...all the other ones...)
c512: COM ATOMS=16889-16892
s1: COM ATOMS=12805-12808
(...all the other ones...)
s512: COM ATOMS=16893-16896
</pre>
*  Then we define [[http://www.plumed.org/doc-v2.5/user-doc/html/_m_o_l_e_c_u_l_e_s.html "molecules"]], here we use the centre of mass as the position of the anion. Note that this is a multicolvar that contains all the vectors from the CF3 groups to the SO3 groups!
<pre>
MOLECULES ...
MOL1=c1,s1
(...all the other ones...)
MOL512=c512,s512
LABEL=CSvectors
... MOLECULES
</pre>
* We condense the multicolvar to the actual collective variable which we later want to bias. This collective variable is the [[http://www.plumed.org/doc-v2.5/user-doc/html/_s_m_a_c.html SMAC variable]]. The switching function SWITCH selects molecules whose centres of mass are within a certain distance. The kernel KERNEL1 selects pairs of anions whose CS vectors are parallel. If we wanted to include antiparallel orientations as well, then we would need to add another kernel, etc. So far, SMAC is another multicolvar; each anion gets assigned a real number which is zero if there are no ordered anions nearby, and 1 if all nearby anions are aligned. Now we want to count those anions which are "ordered enough", which is done by another switching function in MORE_THAN1. It is possible to of course define more than 1 MORE_THAN, and to print them all - this is useful to see which values are reasonable for the switching function.
<pre>
SMAC ...
   SPECIES=CSvectors # this is the "DATA" from the manual
   SWITCH={RATIONAL R_0=9.0 NN=12 D_MAX=13.0} #This should be alright for BMIMOTf
   SWITCH_COORD={RATIONAL R_0=0.001} #That should be the "psi". we don't want to exclude neighbour-less molecules at this stage.
   KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} # the kernelfunction for 0° angle / parallel
   MORE_THAN1={RATIONAL R_0=0.5 NN=6 D_MAX=1.0}
   LABEL=smac
... SMAC
</pre>
*  Now we have a collective variable, smac.morethan-1, which takes values between 0 (no alignment whatsoever) and 512 (all anions are perfectly aligned). We can now bias this collective variable with the [[http://www.plumed.org/doc-v2.5/user-doc/html/_m_e_t_a_d.html METAD keyword]]:
<pre>
METAD ...
 LABEL=metad
 ARG=smac.morethan-1
 HEIGHT=5.0
 BIASFACTOR=5000
 TEMP=298.15
 PACE=200
 GRID_MIN=0.0
 GRID_MAX=512.0
 GRID_BIN=2048
 SIGMA=1.0
 FILE=HILLS
 WALKERS_N=25
 WALKERS_ID=XXXX
 WALKERS_DIR=../HILLS
... METAD
</pre>
HEIGHT is the maximum height of the deposited gaussians (decreases over time because we're using well-tempered metadynamics, see also the BIASFACTOR).<ref>Barducci, A., Bussi, G., & Parrinello, M. (2008). Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. Physical Review Letters, 100(2), 020603. https://doi.org/10.1103/PhysRevLett.100.020603</ref> TEMP is the temperature in K of the simulation, and PACE defines how often gaussians are deposited. Here we deposit on a grid which makes the calculation more efficient, since PLUMED doesn't need to keep track of thousands of gaussians, but rather only the cumulative grid file. SIGMA defines the width of the gaussians. What is interesting here are the WALKERS keywords. I am using 25 walkers here, which means that 25 simulations are running in parallel, each of them depositing gaussians on a common grid. The Walkers ID, XXXX, goes from 0 to 24 accordingly. Each of those should be run in their own directory to avoid a mess, and the common hills file directory should also be a clean separate directory.

If you look at the LAMMPS input file, producerestart.inp, you will see that it is slightly different than before. We are here reading a restart file restart1.lmp rather than the data file, and during the simulation restart files are written as well. Hence we can just resubmit the whole array job with the 25 walkers when it is finished but we need some more sampling. The pair coefficients take up quite a bit of space here, which is why they are in their own LAMMPS input file, pair-p-sc.lmp. The 25 restart files should not be the same even to begin with, otherwise that region gets oversamples. I personally like driving the simulation from the minimum colvar to the maximum colvar with a [[http://www.plumed.org/doc-v2.5/user-doc/html/_m_o_v_i_n_g_r_e_s_t_r_a_i_n_t.html movingrestraint]], meanwhile dumping restart files with LAMMPS so I start roughly evenly spaced in collective variable space.
=Using PLUMED (hack-the-tree)=
Some people are also interested in more complex anions, such as the [NTf2]- anion, which we have looked at in isolation in the first simple example.<ref>Anaredy, R. S., & Shaw, S. K. (2016). Long-Range Ordering of Ionic Liquid Fluid Films. Langmuir, 32(20), 5147–5154. https://doi.org/10.1021/acs.langmuir.6b00304</ref> This anion can exist in two states, cis and trans.<ref>Philippi, F., Pugh, D., Rauber, D., Welton, T., & Hunt, P. A. (2020). Conformational design concepts for anions in ionic liquids. Chemical Science, 11(25), 6405–6422. https://doi.org/10.1039/D0SC01379J</ref> The cis conformer has a larger charge arm (the "dipole moment" when the centre of mass is chosen as the origin of the coordinate system). Thus, the collective variable needs to be large only when the anion is in cis, when other ions in the neighbourhood are in cis as well, if the sulphur-sulphur vectors are aligned, and if the carbon-sulphur vectors are aligned. This is not possible with the current official version of PLUMED, thus the following refers to the experimental [[https://github.com/plumed/plumed2/tree/hack-the-tree hack-the-tree]] version. See also the files in [[:File:BMIMTFSI.zip]].

The definition of groups is similar to [C4C1Im][OTf], only that more vectors are required since the molecular anion is more complex. We use the SS vectors (from one SO2 group to the other, each defined as centre of mass) and the CC-SS vectors (from the centre of mass of both CF3 groups to the centre of mass of both SO2 groups).
What is interesting is the end of the plumedfile:

<pre>
osl: MORE_THAN ARG1=ccssdistance SWITCH={RATIONAL R_0=0.8 NN=6 D_MAX=1.7}
smac_cmap: CONTACT_MATRIX GROUP=ssvectors SWITCH={RATIONAL R_0=9.0 D_0=2.0 NN=12 D_MAX=13.0}
smac_dmat: DOTPRODUCT_MATRIX GROUP1=osl
smac_tor1: TORSIONS_MATRIX GROUP1=ssvectors.x GROUP2=ssvectors.y GROUP3=ssvectors.z POSITIONS=ssvectors
smac_ktor1: KERNEL ARG1=smac_tor1 CENTER=0 SIGMA=0.480 TYPE=gaussian
smac_ktor2: KERNEL ARG1=smac_tor1 CENTER=pi SIGMA=0.480 TYPE=gaussian
smac_ksum1: COMBINE ARG1=smac_ktor1 ARG2=smac_ktor2 PERIODIC=NO
smac_tor2: TORSIONS_MATRIX GROUP1=ccssvectors.x GROUP2=ccssvectors.y GROUP3=ccssvectors.z POSITIONS=ssvectors
smac_ktor3: KERNEL ARG1=smac_tor2 CENTER=0 SIGMA=0.480 TYPE=gaussian
smac_prod: MATHEVAL ARG1=smac_cmap.w ARG2=smac_dmat ARG3=smac_ksum1 ARG4=smac_ktor3 VAR=a,b,c,d FUNC=a*b*c*d PERIODIC=NO
smac: COORDINATIONNUMBER WEIGHT=smac_prod
smac_denom: COORDINATIONNUMBER WEIGHT=smac_cmap.w
smac_avg: MATHEVAL ARG1=smac ARG2=smac_denom FUNC=x/y PERIODIC=NO
smac_mt: MORE_THAN ARG=smac_avg SWITCH={RATIONAL R_0=0.5 NN=6 D_MAX=1.0} 
final: COMBINE ARG=smac_mt PERIODIC=NO
</pre>

I will explain the commands in reverse order, since I believe that this is easier to understand. the "final" colvar is the sum of all "smac_mt", which in turn is obtained from "smac_avg" by applying a switching function. This switching function essentially amplifies changes in "smac_avg". Both "smac_avg" and "smac_mt" are already multicolvars, i.e. one value per anion. "smac_avg" smoothly varies between 0 (disordered surroundings) and 1 (perfectly ordered surroundings) for every anion, whereas "smac_mt" should be 0 if the anion surroundings are too disordered (even if "smac_avg" is not 0) and 1 if they are ordered enough (even if "smac_avg" is not 1). "smac_avg" Itself is obtained by dividing "smac" by "smac_denom", which has the purpose of normalising smac to be between 0 and 1. Correspondingly, "smac_denom" is the number of neighbours. Thus, "smac" is the variable (more precisely, multicolvar) that counts how many other anions in the neighbourhood of each anion fulfill the above criteria. To this end, for each pair of anions, the product of "smac_cmap.w", "smac_dmat", "smac_ksum1" and "smac_ktor3" is calculated, and summed up per anion (done by "COORDINATIONNUMBER WEIGHT=smac_prod"). Here, "smac_cmap" decides whether a pair of anions is close enough by means of a switching function on the distance, "smac_dmat" selects only anions which are cis by means of a switching function on the distance between the centre of mass of the CF3 groups and the centre of mass of the SO2 groups, "smac_ksum1" selects anions whose SS vectors are either parallel ("smac_ktor1") or antiparallel ("smac_ktor2"), and finally "smac_ktor3" selects anions whose CC-SS vectors are aligned.

=References=
<References>