 |
 |
 |
| |
|
|
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.
Once this tutorial is completed students will have used PLUMED coupled with LAMMPS to perform an investigation of their own design.
The TARBALL for this project contains the following files:
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.
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.
You can find the instructions about how to install LAMMPS patched with PLUMED here:
https://github.com/plumed/conda#install-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:
touch plumed.dat /Users/gareth/Plumed-tutorial/LAMMPS/lammps-5Jun19/build/lmp < lammps-equil.in
(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:
# Specify that we are in natural units UNITS NATURAL( default=off ) use natural units # 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 ARGthe input for this action is the scalar output from one or more other actions. =cell.ax,cell.ay,cell.az POWERScompulsory keyword ( default=1.0 ) the powers to which you are raising each of the arguments in your function =2,2,2 PERIODICcompulsory keyword if the output of your function is periodic then you should specify the periodicity of the function. =NO bbb2: COMBINE ARGthe input for this action is the scalar output from one or more other actions. =cell.bx,cell.by,cell.bz POWERScompulsory keyword ( default=1.0 ) the powers to which you are raising each of the arguments in your function =2,2,2 PERIODICcompulsory keyword if the output of your function is periodic then you should specify the periodicity of the function. =NO ccc2: COMBINE ARGthe input for this action is the scalar output from one or more other actions. =cell.cx,cell.cy,cell.cz POWERScompulsory keyword ( default=1.0 ) the powers to which you are raising each of the arguments in your function =2,2,2 PERIODICcompulsory keyword if the output of your function is periodic then you should specify the periodicity of the function. =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 ARGthe input for this action is the scalar output from one or more other actions. =ee,aaa,bbb,ccc FILEthe name of the file on which to output these quantities =colvar STRIDEcompulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output =10
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.
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:
origin: FIXEDATOM ATcompulsory keyword coordinates of the virtual atom =0,0,0 dx: XDISTANCES ... ATOMS1the atoms involved in each of the distances you wish to calculate. =1,origin ATOMS2the atoms involved in each of the distances you wish to calculate. =2,origin BETWEENcalculate the number of values that are within a certain range. ={GAUSSIAN LOWER=0 UPPER=5.4373815 SMEAR=0.1} ... da: XDISTANCES ... ATOMS1the atoms involved in each of the distances you wish to calculate. =1,origin ATOMS2the atoms involved in each of the distances you wish to calculate. =2,origin ATOMS3the atoms involved in each of the distances you wish to calculate. =3,origin ATOMS4the atoms involved in each of the distances you wish to calculate. =4,origin BETWEEN1calculate the number of values that are within a certain range. ={GAUSSIAN LOWER=0 UPPER=5.4373815 SMEAR=0.1} BETWEEN2calculate the number of values that are within a certain range. ={GAUSSIAN LOWER=-5.4373815 UPPER=0 SMEAR=0.1} ...
(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:
UNITS NATURAL( default=off ) use natural units # Include the file that defines our collective variable INCLUDE FILEcompulsory keyword file to be included =__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 ARGthe input for this action is the scalar output from one or more other actions. =__FILL__ FUNCcompulsory keyword the function you wish to evaluate =x-y PERIODICcompulsory keyword if the output of your function is periodic then you should specify the periodicity of the function. =NO RESTRAINT ARGthe input for this action is the scalar output from one or more other actions. =ccc ATcompulsory keyword the position of the restraint =0.0 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 =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 ARGthe input for this action is the scalar output from one or more other actions. =__FILL__ AT0compulsory keyword AT\f$x\f$ is equal to the position of the restraint at time STEP\f$x\f$. =230 AT1compulsory keyword AT\f$x\f$ is equal to the position of the restraint at time STEP\f$x\f$. =450 AT2compulsory keyword AT\f$x\f$ is equal to the position of the restraint at time STEP\f$x\f$. =230 KAPPA0compulsory keyword KAPPA\f$x\f$ is equal to the value of the force constants at time STEP\f$x\f$. =100 STEP0compulsory keyword This keyword appears multiple times as STEP\f$x\f$ with x=0,1,2,...,n. =0 STEP1compulsory keyword This keyword appears multiple times as STEP\f$x\f$ with x=0,1,2,...,n. =50000 STEP2compulsory keyword This keyword appears multiple times as STEP\f$x\f$ with x=0,1,2,...,n. =100000 PRINT ARGthe input for this action is the scalar output from one or more other actions. =da.*,dx.between,res.* FILEthe name of the file on which to output these quantities =colvar STRIDEcompulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output =10
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.
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:
1.8.17