Skip to content
Snippets Groups Projects
Commit 237d8109 authored by Gareth Tribello's avatar Gareth Tribello
Browse files

Added tutorial on PLUMED + LAMMPS

parent 974d6931
No related branches found
No related tags found
No related merge requests found
......@@ -936,3 +936,12 @@ diagonalize
diagonalized
diagonalization
vmdrc
lammps
equilibration
xdistances
homebrew
ffmpeg
gsl
lmp
equilibrating
equilibrates
/**
\page lugano-6a Lugano tutorial: Using PLUMED with LAMMPS
\section lugano-6a-aim Aims
In this tutorial I will show you how you can use PLUMED to perform a simulation in which a mixture of
Lennard Jones particles are forced to separate and mix. You then will be encouraged to investigate
how the work done during this process depends on the simulation parameters or the parameters that are
used for the enhanced sampled algorithm.
\section lugano-6a-lo Objectives
Once this tutorial is completed students will
- Know how to compile LAMMPS together with PLUMED
- Have used PLUMED coupled with LAMMPS to perform an investigation of their own design
\section lugano-6a-resources Resources
The \tarball{lugano-4} for this project contains the following files:
- lammps-equilibration.in : The input for LAMMPS that we will use for our setup and equilibration steps of our investigation
- lammps-production.in : The input for LAMMPS that will be used for the production steps of our investigation
- xdistances.dat : A PLUMED file that defines the collective variables that we will use in the investigation. We will include this file in the PLUMED input files we write for this tutorial.
This tutorial has been tested on v2.5 but it should also work with other versions of PLUMED.
Also notice that the `.solutions` direction of the tarball contains correct input files for the exercises.
Please only look at these files once you have tried to solve the problems yourself. Similarly the tutorial
below contains questions for you to answer that are shown in bold. You can reveal the answers to these questions
by clicking on links within the tutorial but you should obviously try to answer things yourself before reading these
answers.
\section lugano-6a-intro Introduction
This tutorial is somewhat different to the others that you have done. The tutorial will show you how to run a simulation on a binary Lennard Jones system using PLUMED and
LAMMPS. In this simulation the two types of atoms will be forced to separate and mix. You will then be encouraged to design your investigation to further investigate this phenomenon
by means of simulation. There are a variety of different things you can investigate. You could investigate how the work done to separate the two components depends on the parameters of
the Lennard Jones model or temperature, you could investigate whether you can reduce the amount of hysteresis that is observed with the current simulation parameters. You may even decide
that the method proposed is not a particularly good way to investigate the phenomenon of phase separation and propose a different method for simulating the process.
\section lugano-6a-exercises Exercises
\subsection lugano-6a-installation Installing LAMMPS with PLUMED
You can get a copy of LAMMPS from this website:
https://lammps.sandia.gov/download.html
In preparing the material for this exercises I downloaded the latest stable release C++ source. You can try to download the pre-built executables. I did
not do this, however, as I have a Mac and have broken many things on it in the past by using homebrew. I therefore am, as a rule, distrustful of homebrew.
I thus instead made a directory called:
\verbatim
/Users/gareth/Plumed-tutorial/LAMMPS
\endverbatim
On my directory, I copied the downloaded LAMMPS tar ball into this directory and extracted it from the tar ball using the command:
\verbatim
tar -xvf lammps-5Jun19.tar
\endverbatim
Once this operation was completed I found that when I ls the LAMMPS directory I created I have the following files:
\verbatim
lammps-5Jun19 lammps-5Jun19.tar
\endverbatim
Before building LAMMPS I rebuilt a version of PLUMED to link with LAMMPS. I did this using the following commands:
\verbatim
./configure --enable-modules=all --prefix=/Users/gareth/Plumed-tutorial/LAMMPS
make
make install
\endverbatim
I set the --prefix flag here to the directory I created to hold my compiled version of LAMMPS here rather than using the default in order to be able to control exactly what version of PLUMED I was linking
LAMMPS with. Setting the prefix in this way ensures that once the above three commands have been run the contents of the LAMMPS directory I created in the first step are as follows:
\verbatim
bin include lammps-5Jun19 lammps-5Jun19.tar lib
\endverbatim
To build LAMMPS we now issue the following commands:
\verbatim
cd lammps-5Jun19
mkdir build
cd build
export PKG_CONFIG_PATH=$PKG_CONFIG_PATH:/Users/gareth/Plumed-tutorial/LAMMPS/lib/pkgconfig
cmake -D PKG_MANYBODY=yes -D PKG_KSPACE=yes -D PKG_MOLECULE=yes -D PKG_RIGID=yes -D PKG_USER-PLUMED=yes ../cmake
\endverbatim
LAMMPS requires cmake, ffmpeg, gsl and xdrfile, which for some reason were not installed on my computer. I therefore installed these packages using the commands:
\verbatim
sudo port install cmake
sudo port install ffmpeg
sudo port install gsl
sudo port install xdrfile
\endverbatim
Finally now that cmake was installed I could build LAMMPS using the command:
\verbatim
make
\endverbatim
This generates an executable called lmp in the build directory, which is used for the remainder of the exercises in this tutorial.
\subsection lugano-6a-setup-equilibration Setting up the simulation and equilibrating
Having compiled LAMMPS we can now run a simulation with LAMMPS. If we take the input lammps-equilibration.in that we downloaded from the tarball we can run this simulation using the following commands:
\verbatim
touch plumed.dat
/Users/gareth/Plumed-tutorial/LAMMPS/lammps-5Jun19/build/lmp < lammps-equil.in
\endverbatim
(In the second command here you obviously need to specify the location where you compiled LAMMPS instead of the location where I compiled LAMMPS).
As you can see if you visualize the trajectory output.xyz that is generated by this simulation this input creates an initial configuration containing 1000 atoms.
There are 500 Ar atoms and 500 Ne atoms and in the initial configuration all the Ne atoms are placed in one half of the simulation box while the Ar atoms are placed
in the other half of the box. During the simulation you have just run you should observe that these two types of atom mix together as the system equilibrates. We are now
going to run this simulation again but we will now use PLUMED to monitor the extent how the system comes to equilibrium. In particular, we are going to use a filled in
version of the following input file to monitor the energy and the cell parameters over the course of the simulation:
\plumedfile
# Specify that we are in natural units
UNITS NATURAL
# Retrieve the value of the potential energy
ee: ENERGY
# Retrieve the 3 cell vectors
cell: CELL
# Compute the square moduli of three cell vectors
aaa2: COMBINE ARG=cell.ax,cell.ay,cell.az POWERS=2,2,2 PERIODIC=NO
bbb2: COMBINE ARG=cell.bx,cell.by,cell.bz POWERS=2,2,2 PERIODIC=NO
ccc2: COMBINE ARG=cell.cx,cell.cy,cell.cz POWERS=2,2,2 PERIODIC=NO
# Compute the moduli of the three cell vectors
aaa: __FILL__
bbb: __FILL__
ccc: __FILL__
# Print the energy and the moduli of the three cell vectors
PRINT ARG=ee,aaa,bbb,ccc FILE=colvar STRIDE=10
\endplumedfile
Rerun the equilibration calculation that you just performed but with a plumed.dat file that now contains a filled-in version of the above input. You should observe that the
energy and cell lengths that are output to the colvar file initially change substantially, which is normal as we are running in the NPT ensemble. After a while, however, these
quantities equilibrate and then simply fluctuate around some fixed value.
\subsection lugano-6a-production Simulating phase separation
Having equilibrated our system lets now simulate phase separation. What we are going to do is force the Ne atoms to be in a different part of the simulation cell from the Ar atoms.
We will thus force the system to be a un-mixed, two phase system from being a mixed, single-phase system. In other words, we are going to force the system to go back to a configuration
much like the one that we started our simulation from. Before we get on to running the simulation, however, we are going to need to copy a few pieces to the directory. In particular, we
need the following files:
- We need to start our simulation from the equilibrated configuration that we arrived in at the end of the simulation run that we have just performed. We can start in this why by copying
the lammps restart file that was output from that simulation: lj.equilibration.restart.100000
- We need the lammps input file called lammps-production.in that we obtained from the tarball that we downloaded at the beginning of the exercise
- We need the PLUMED file called xdistances.dat that we downloaded from the tar ball. This file defines one virtual atom and three collective variables and looks something like the following:
\plumedfile
origin: FIXEDATOM AT=0,0,0
dx: XDISTANCES ...
ATOMS1=1,origin
ATOMS2=2,origin
BETWEEN={GAUSSIAN LOWER=0 UPPER=5.4373815 SMEAR=0.1}
...
da: XDISTANCES ...
ATOMS1=1,origin
ATOMS2=2,origin
ATOMS3=3,origin
ATOMS4=4,origin
BETWEEN1={GAUSSIAN LOWER=0 UPPER=5.4373815 SMEAR=0.1 LABEL=left}
BETWEEN2={GAUSSIAN LOWER=-5.4373815 UPPER=0 SMEAR=0.1 LABEL=right}
...
\endplumedfile
(In the actual file there are 500 ATOMS keywords in the first XDISTANCES command and 1000 ATOMS keywords in the second XDISTANCES command, however, which is why I have not shown the whole file.)
These XDISTANCES commands compute the x-component of the vector connecting the pair of specified atom. The BETWEEN keywords then specify that we want to calculate how many of these x-components
are between a particular lower and upper bound. In this particular case, the FIXEDATOM is placed at the origin (center) of the simulation cell and the CVs monitor how many of the x-components of
the vector connecting the origin and each atom are positive and how many are negative. We thus divide the box into equally sized regions and count how many atoms are in each of these regions.
To this correctly, however, you have to ensure that the parameters for the BETWEEN keywords are set correctly. Wherever you see 5.4373815 in the input file you need to replace this value by
half the simulation box length that you will use for these production calculations. You can get the simulation box for the initial configuration (the final configuration from the
equilibration simulation) from the colvar file that was output from your equilibration simulation. You simply need to look at the value that this quantity took in the last line of the colvar file.
This quantity will stay fixed throughout your simulation as in this production run we are going to use the NVT ensemble.
With these modifications in place lets now look at the PLUMED input file we are going to use for this calculation:
\plumedfile
UNITS NATURAL
# Include the file that defines our collective variable
INCLUDE FILE=__FILL__
# Add a restraint that ensures that the number of atoms in the left part of the box stays the same as the number of atoms
# in the right part of the box
ccc: CUSTOM ARG=__FILL__ FUNC=x-y PERIODIC=NO
RESTRAINT ARG=ccc AT=0.0 KAPPA=10
# Add a moving restraint that forces the Ar atoms (1-500) to separate from the Ne atoms and to then mix again with the Ne atoms
res: MOVINGRESTRAINT ARG=__FILL__ AT0=230 AT1=450 AT2=230 KAPPA0=100 STEP0=0 STEP1=50000 STEP2=100000
PRINT ARG=da.*,dx.between,res.* FILE=colvar STRIDE=10
\endplumedfile
Fill in the blanks in the input above and then try to run the simulation. You should see as the simulation progresses the Ne and Ar atoms un-mix into two separate phases. They then remix
and by the end of the simulation there is only one phase once more.
\section lugano-6a-conclusions Conclusions
This tutorial has shown you how you can run a simulation in which an interesting physical phenomenon is forced to occur. The phenomenon in question is phase separation - the simulation was
started from a mixed phase of Ar and Ne. These two chemical components were then forced to separate during the simulation and a two phase system was formed with all the Ne atoms in one
part of the simulation box and all the Ar atoms in the other part of the box. Your task is now to investigate this phenomenon further. Some questions you might want to ask yourself as you
design your investigation are the following:
- How do you get a quantitative estimate of the energetic cost for performing the phase separation reaction?
- How does this energetic cost depend on the simulation parameters e.g. ratio of Ar/Ne, temperature, pressure?
- How does this energetic cost depend on the ratio of force field parameters?
- Is the method that you were taught in this tutorial the best method for studying this process? What other methods might we use to investigate phase separation?
*/
link: @subpage lugano-6a
description: An exercise on running PLUMED with LAMMPS
additional-files: lugano-6a
# Specify that we are in natural units
UNITS NATURAL
# Retrieve the value of the potential energy
ee: ENERGY
# Retrieve the 3 cell vectors
cell: CELL
# Compute the square moduli of three cell vectors
aaa2: COMBINE ARG=cell.ax,cell.ay,cell.az POWERS=2,2,2 PERIODIC=NO
bbb2: COMBINE ARG=cell.bx,cell.by,cell.bz POWERS=2,2,2 PERIODIC=NO
ccc2: COMBINE ARG=cell.cx,cell.cy,cell.cz POWERS=2,2,2 PERIODIC=NO
# Compute the moduli of the three cell vectors
aaa: CUSTOM ARG=aaa2 FUNC=sqrt(x) PERIODIC=NO
bbb: CUSTOM ARG=bbb2 FUNC=sqrt(x) PERIODIC=NO
ccc: CUSTOM ARG=ccc2 FUNC=sqrt(x) PERIODIC=NO
# Print the energy and the moduli of the three cell vectors
PRINT ARG=ee,aaa,bbb,ccc FILE=colvar STRIDE=10
UNITS NATURAL
# Include the file that defines our collective variable
INCLUDE FILE=xdistances.dat
# Add a restraint that ensures that the number of atoms in the left part of the box stays the same as the number of atoms
# in the right part of the box
ccc: CUSTOM ARG=da.right,da.left FUNC=x-y PERIODIC=NO
RESTRAINT ARG=ccc AT=0.0 KAPPA=10
# Add a moving restraint that forces the Ar atoms (1-500) to separate from the Ne atoms and to then mix again with the Ne atoms
res: MOVINGRESTRAINT ARG=dx.between AT0=230 AT1=450 AT2=230 KAPPA0=100 STEP0=0 STEP1=50000 STEP2=100000
PRINT ARG=da.*,dx.between,res.* FILE=colvar STRIDE=10
#
# define units
#
units lj
#
# specify periodic boundary conditions
#
boundary p p p
#
# define atom_style
# full covers everything
#
atom_style full
#
# Setup the simulatioon box
#
variable blen equal 12.5
region boxid block 0.0 ${blen} 0.0 ${blen} 0.0 ${blen}
create_box 2 boxid
#
# Split the simulation into two regions
#
variable half_x equal ${blen}/2.0
region boxa block 0.0 ${half_x} 0.0 ${blen} 0.0 ${blen}
region boxb block ${half_x} ${blen} 0.0 ${blen} 0.0 ${blen}
#
# specify initial positions of atoms
# sc = simple cubic
# 0.5 = density in lj units
#
lattice sc 0.50
# place atoms of type 1 in boxa
# place atoms of type 2 in boxb
create_atoms 1 region boxa
create_atoms 2 region boxb
#
# define mass
#
# mass of atom type 1
mass 1 1.0
mass 2 1.0
#
# specify initial velocity of atoms
# group = all
# reduced temperature is T = 1.0 = lj-eps/kb
# seed for random number generator
# distribution is gaussian (e.g. Maxwell-Boltzmann)
#
velocity all create 1.0 87287 dist gaussian
#
# specify interaction potential
# pairwise interaction via the Lennard-Jones potential with a cut-off at 2.5 lj-sigma
#
pair_style lj/cut 2.5
# specify parameters between atoms of type 1 with an atom of type 1
# epsilon = 1.0, sigma = 1.0, cutoff = 2.5
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 2 2 1.2 0.9 2.5
#
# specify parameters for neighbor list
# rnbr = rcut + 0.3
#
neighbor 0.3 bin
#
# Specify timestep
#
timestep 0.001
#
# Write thermodynamic parameters to log
#
thermo_style custom step pe ke etotal temp press density
thermo 1000
thermo_modify norm no
#
# Setup plumed and that we are going to run npt equilibration
#
fix 1 all plumed plumedfile plumed.dat outfile p.log
fix 2 all npt temp 1.0 1.0 0.1 iso 0.0 0.0 0.4
#
# save trajectory and restart file
# dumpid = 1
# filename = output.xyz
#
dump 1 all xyz 1000 output.xyz
dump_modify 1 element Ar Ne
restart 100000 lj.equilibration.restart
#
# Run 5000000 of equilibration
#
run 100000
#
# define units
#
units lj
#
# specify periodic boundary conditions
#
boundary p p p
#
# define atom_style
# full covers everything
#
atom_style full
#
# Read the restart file and restart the counter on timesteps
#
read_restart lj.equilibration.restart.100000
reset_timestep 0
#
# define mass
#
# mass of atom type 1
mass 1 1.0
mass 2 1.0
#
# specify interaction potential
# pairwise interaction via the Lennard-Jones potential with a cut-off at 2.5 lj-sigma
#
pair_style lj/cut 2.5
# specify parameters between atoms of type 1 with an atom of type 1
# epsilon = 1.0, sigma = 1.0, cutoff = 2.5
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 2 2 1.2 0.9 2.5
#
# specify parameters for neighbor list
# rnbr = rcut + 0.3
#
neighbor 0.3 bin
#
# Specify timestep
#
timestep 0.001
#
# Write thermodynamic parameters to log
#
thermo_style custom step pe ke etotal temp press density
thermo 1000
thermo_modify norm no
#
# Setup plumed and that we are going to run npt equilibration
#
fix 1 all plumed plumedfile plumed.dat outfile p.log
fix 2 all nvt temp 1.0 1.0 0.1
#
# save trajectory and restart file
# dumpid = 1
# filename = output.xyz
#
dump 1 all xyz 1000 output.xyz
dump_modify 1 element Ar Ne
restart 100000 lj.prod.restart
#
# Run 5000000 of equilibration
#
run 100000
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment