Table of Contents

Lennard-Jones liquids

In this exercise, you will simulate a fluid of monoatomic particles that interact with a Lennard-Jones potential. The method to be used is molecular dynamics (MD) with periodic boundary conditions using CP2K. The aim is to explore the method, calculate the radial distribution function $g(r)$

You are expected to hand in the respective plots by email, ONLY PDF format.

Run your first simulation using CP2K

When you check CP2K's features and the outline of the lecture you will notice that there are many levels of theory, methods, and possibilities to combine them. This results in a large number of possible options and coefficients to setup, control and tune a specific simulation. Together with the parameters for the numerical solvers, this means that an average CP2K configuration file will contain quite a number of options (even though for many others the default value will be applied) and not all of them will be discussed in the lecture or the exercises.

The CP2K Manual is the complete reference for all configuration options. Where appropriate you will find a reference to the respective paper when looking up a specific keyword/option.

To get you started, we will do a simple exercise using Molecular Mechanics (that is: a classical approach). The point is to get familiar with the options, organizing and editing the input file and analyze the output.

Background

You are expected to carry out an MD simulation of Lennard-Jones fluid containing mono-atomic particles. The potential energy expression used is

$U(x)=4\epsilon \left [\left ( \frac{\sigma }{r_{ij}} \right )^{12}- \left ( \frac{\sigma }{r_{ij}} \right )^{6} \right ]$

where $\epsilon$ is the well depth, σ is related to the minimum of Lennard-Jones potential, rij is the distance between atoms i and j.

Periodic boundary conditions should be used in this simulation. The atom near the ”edge” of the simulation box interacts with atoms contained in the neighboring image of the box. In computer simulations, one of these is the original simulation box, and others are copies called images. During the simulation, only the properties of the original simulation box need to be recorded and propagated. The minimum-image convention is a common form of PBC particle bookkeeping in which each individual particle in the simulation interacts with the closest image of the remaining particles in the system. To prevent artifacts, it requires a cut-off value of rij for the L-J potential. For realistic results, the cut-off should be less than the half of the simulation box size and over σ Radial distribution function, (or pair correlation function) $g(r)$ in a system of particles (atoms, molecules, colloids, etc.), describes how density varies as a function of distance from a reference particle.

Part I: Lennard-Jones potential energy

In the first part of this set of exercises you are asked to perform some single point calculation, where only the energy is computed for a set of atomic configurations. Specifically, you have to compute the energy of systems of noble gas dimers, modeled as Lennard-Jones particles, and determine the potential energy as a function of the distance between the two atoms. You are asked to do this for Ar, Ne, Kr.

1. Step

We suggest to create a directory ex1 and change to it:

mkdir ex1
cd ex1

And inside this directory make three others, one for each element (Ar, Ne, Kr).

2. Step

Now you will run single point calculations on the specified systems. We provide you with the input file for the Ar dimers. Copy the input file below into the Ar subdirectory and run a CP2K calculation with the following command:

$ cp2k.sopt -i argon.inp -o argon.out 
&GLOBAL
  PROJECT ar108           #Project Name
  RUN_TYPE energy         #Calculation Type : MD (molecular dynamics), GEO_OPT (Geometry Optimization), Energy (Energy Calculation)
&END GLOBAL

&FORCE_EVAL
  METHOD FIST             #Method to calculate force: Fist (Molecular Mechanics), QS or QUICKSTEP (Electronic structure methods, like DFT)
  &MM
    &FORCEFIELD
    &SPLINE
     EMAX_SPLINE 100000
    &END
      &CHARGE             #charge of the MM atoms
        ATOM Ar           #Defines the atomic kind of the charge
        CHARGE 0.0        #Defines the charge of the MM atom in electron charge unit
      &END
      &NONBONDED
        &LENNARD-JONES    #LENNARD-JONES potential type.Functional form: V(r) = 4.0 * EPSILON * [(SIGMA/r)^12-(SIGMA/r)^6]
          atoms Ar Ar     #Defines the atomic kind involved in the nonbonded potential
          EPSILON 119.8   #Defines the EPSILON parameter of the LJ potential (K_e)
          SIGMA 3.405     #Defines the SIGMA parameter of the LJ potential (Angstrom)
          RCUT 8.4        #Defines the cutoff parameter of the LJ potential
        &END LENNARD-JONES
      &END NONBONDED
    &END FORCEFIELD
    &POISSON              # Poisson solver
      &EWALD
        EWALD_TYPE none
      &END EWALD
    &END POISSON
  &END MM

      &SUBSYS                 ! system description
       &CELL
         ABC [angstrom] 40 40 40
         PERIODIC NONE
       &END CELL
      &COORD
        UNIT angstrom
        Ar  0 0 0
        Ar  4 0 0
      &END COORD
  &END SUBSYS

&END FORCE_EVAL
Instructions starting with an ampersand & start a section and must be terminated with an &END SECTION-NAME. Other instructions are called keywords. The indentation is ignored but recommended for readability.

After having performed the calculation above for a given distance, you will have to run this type of calculation at different distances and to check how the potential changes.

A simple way to generate the different input files is using shell scripting in combination with sed (the stream editor):

for d in $(seq 3.0 0.1 9); do
  sed -e "s|4 0 0|${d} 0 0|" argon.inp > energy_${d}A.inp
  cp2k.sopt -i energy_${d}A.inp -o energy_${d}A.out
  awk '/Total FORCE_EVAL/ { print "'"${d}"'", $9; }' energy_${d}A.out >> Pot_En_vs_distance.dat
done

Plot distance vs. potential energy and find the minimum in energy, which corresponds to the equilibrium distance. After having done it, you can calculate the minimum analytically as well.

3. Step

Now we will change the atoms and we will repeat the same procedure we have already seen. To do this you need to change the L-J parameters. Change directory to the one related to Ne or Kr, copy the input file into the directory and remember to change the atom type in the coordinate section as well.

Neon: Ne, $\epsilon$ = 3.0840 meV and $\sigma$ = 0.2782 nm. Kripton: Kr, $\epsilon$ = 14.391 meV and $\sigma$ = 0.3633 nm.
In these calculations you will not use the CP2K default units for the L-J parameters ([K_e] and [angstrom]). CP2K allows you to add the unit you want inside squared brackets (here you have the list of the units that you can use https://manual.cp2k.org/cp2k-5_1-branch/units.html). Keep in mind that you have to transform meV in eV. Example:
        atoms Kr Kr
        EPSILON [eV] 0.014391
        SIGMA [nm] 0.3633

Repeat the procedure explained above, compare the potentials energy curves for the three dimers and determine the equilibrium distance and the stationary point.

Part II: Set up MD simulation

Now you have to run an actual MD simulation for liquid Ar modeled as Lennard-Jones fluid. Copy the input below to a new directory as done previously.

argon.inp
##  It's highly recommended to go 
##  https://manual.cp2k.org/ 
##  and learn how to set up CP2K 
##  calculation correctly using manual.

&GLOBAL
  PROJECT ar108           #Project Name
  RUN_TYPE md             #Calculation Type : MD (molecular dynamics), GEO_OPT (Geometry Optimization), Energy (Energy Calculation)
&END GLOBAL

&MOTION
  &MD
    ENSEMBLE NVE          #The ensemble for MD propagation, NVE (microcanonical), NVT (canonical), NPT_I (NPT with isotropic cell) 
    STEPS 30000           #The number of MD steps to perform
    TIMESTEP 5.           #The length of an integration step (fs)
    TEMPERATURE 85.       #The temperature in K used to initialize the velocities with init and pos restart, and in the NPT/NVT simulations
  &END MD
&END MOTION

&FORCE_EVAL
  METHOD FIST             #Method to calculate force: Fist (Molecular Mechanics), QS or QUICKSTEP (Electronic structure methods, like DFT)
  &MM
    &FORCEFIELD
      &CHARGE             #charge of the MM atoms 
        ATOM Ar           #Defines the atomic kind of the charge
        CHARGE 0.0        #Defines the charge of the MM atom in electron charge unit
      &END
      &NONBONDED
        &LENNARD-JONES    #LENNARD-JONES potential type.Functional form: V(r) = 4.0 * EPSILON * [(SIGMA/r)^12-(SIGMA/r)^6]
          atoms Ar Ar     #Defines the atomic kind involved in the nonbonded potential 
          EPSILON 119.8   #Defines the EPSILON parameter of the LJ potential (K_e) 
          SIGMA 3.405     #Defines the SIGMA parameter of the LJ potential (Angstrom)
          RCUT 8.4        #Defines the cutoff parameter of the LJ potential 
        &END LENNARD-JONES
      &END NONBONDED
    &END FORCEFIELD
    &POISSON              # Poisson solver
      &EWALD
        EWALD_TYPE none
      &END EWALD
    &END POISSON
  &END MM
  &SUBSYS
    &CELL
      ABC 17.1580 17.158 17.158    #Simulation box size
    &END CELL

    &COORD
   Ar -8.53869012951987116 -15.5816257770688615 2.85663672298278293
   Ar 1.53007304829383051 9.28528179040142554 11.1777824543317941
   Ar 11.9910225119590699 -7.48825329565798015 -9.96545306345559823
   Ar -12.6782400030290496 -3.34105872014234606 4.07471097818485806
   Ar -1.77046254278594462 -0.232459464264201887 13.2012946017273016
   Ar 8.01761371186688443 -2.57249587730733298 -4.12720554747711432
   Ar 8.57849517232300052 4.01396664624232002 5.57368821983998419
   Ar -3.89200679277030925 -10.2930917801117356 -6.98640232289045482
   Ar -3.35457160564444568 -16.1119619276890056 16.1358515626317427
   Ar 9.78957155103081966 -16.2628264194939263 -5.69790857071688350
   Ar 0.505143495414835719 -4.22978415759568183 12.4854171634357307
   Ar 15.5632243939617503 -7.98048905093276240 2.20994708545912832
   Ar -5.40741643995084953 -2.64764457113743079 -0.681485212640798199
   Ar -0.983719068448489081E-01 -1.73674004862212694 -7.11915545117132265
   Ar 7.52655781331927187 -5.52969969672439632 -12.8886150439489313
   Ar -5.45655410995716128 0.564445754429787061 2.03902510096247536
   Ar -11.8590998267164665 3.40407446386207724 3.72687933934436399
   Ar 16.7175362589401821 -7.47132377347522780 -1.02274476672697889
   Ar -20.4572129717055340 -5.73700807719791683 4.81845086375497811
   Ar 14.8485522289272627 -1.41608633045414667 -16.0839111490847451
   Ar 8.04379470511429595 -8.14033814842439263 -4.75543123809189261
   Ar 12.2738439612049568 -1.70589834674486429 12.9622486199573572
   Ar -0.421851806372696092 -11.1177490353157999 20.4545363332536283
   Ar 2.28194341698637571 5.92083917539752136 -11.1732449877738436
   Ar -13.9648466918215064 8.77923885764231926 8.07373370482465091
   Ar  -10.3147439499058429 6.38529561240966004 -15.3411964215061527
   Ar  -2.71899964647918457 -21.4890074469143855 10.8678096818980006
   Ar  -17.7923879123397271 -10.7840901151121251 -4.83954996524571968
   Ar  5.23494138507746420 -6.79222906792632841 -6.07187690814296133
   Ar  3.52448750638480446 -10.1225951872349782 2.96829048662758721
   Ar  -16.1586602901979361 -5.18274316385346445 8.57072694078649455
   Ar  -5.80982824422251287 4.32640193501643733 2.55599101868223322
   Ar  6.29160109084684382 27.5741337288405717 15.0246410590392632
   Ar  -3.18741711710350684 23.2996469099840624 -16.8034854143018748
   Ar  -4.20225755039435622 9.36037725943080190 16.5891306154890081
   Ar  -7.64392908749747946 -9.52432384411045341 -29.8228731471089645
   Ar  0.545352525792712428 13.9240554617015260 -0.383786780333776500
   Ar  -5.27432886808646906 -5.53813781787395865 -20.3014703747109415
   Ar  22.9921850152838871 6.78619371666398941 -1.98289905290632484
   Ar  19.7720034229251880 -10.2373337687313679 -3.33081818566269172
   Ar  0.156776902886395425 6.59630118110908725 8.90749062505743083
   Ar  5.57937381862174053 0.233106223140015806 1.02752287819280941
   Ar  -3.64343561800208793 3.96448881012491006 25.8752124557059595
   Ar  -0.248491698112870391 20.4489725648023182 -2.51220445353457666
   Ar  2.93626708600658270 0.859812213376437984 9.96743307236779508
   Ar  3.30384315693043895 -2.92421266591109408 -6.34927042371499883
   Ar  -6.15490235244551265 -6.84961480075890883 -6.46204144605644260
   Ar  -23.2388291761596619 -28.1213094673208666 7.13721047187827917
   Ar  4.11526291325474780 2.71564143367947342 -0.852030043744060328
   Ar  14.6194148692240713 2.80815182256426210 1.93601975975151541
   Ar  18.9667954753247869 16.5700888519293095 13.3423444868082761
   Ar  -28.6124161416877705 2.84353637083477562 -9.23601973326721648
   Ar  -5.97004594556101331 -16.2230172568109978 -9.22928061840017477
   Ar  10.0481077882725955 16.3854819569745231 5.12578711346205651
   Ar  -7.22508507825336643 6.34615422233080650 -0.680757463730119028
   Ar  -12.0138912984383506 -10.4653110276797570 -6.43434787584580103
   Ar  -8.53169926903037457 12.8976589212818862 -0.890361252446473683
   Ar  -25.3700692950848676 3.33119906434656077 -7.93917685683272722
   Ar  2.90163480643285920 -12.9668181360039672 -7.94907759259854707
   Ar  -13.5963940986222074 11.9896580951935974 -3.55068754869933789
   Ar  13.1416029517342476 4.97143783446568488 -3.50841252726170705
   Ar  -13.3295460955805254 16.0410015777677764 7.05282797577515375
   Ar  -4.15068335494176122 -19.5111913798076593 21.0255971827539376
   Ar  -8.42944270819351793 16.3065160593537009 -18.2887817284733885
   Ar  0.788636333898691255 9.59016836817029095 22.1772606194495872
   Ar  -2.92606778628861974 -7.97408054890791007 -21.3519900334304964
   Ar  -6.39959865978756426 -4.56280461803643256 2.75533571094951402
   Ar  -4.04423878093174860 -14.9275965394452506 -5.58561473738824965
   Ar  -27.8524912514281482 0.802052180719123098 -3.02663789713126441
   Ar  -2.83966529645897792 7.11627121253196915 6.18547332762273783
   Ar  -8.68327887612401739 -6.67088300493855879 -9.15815450219801264
   Ar  -11.9620847111198501 -2.20956249614563038 -1.83979975374852245
   Ar  22.6848553724304907 12.2047209420099971 1.01238797839832362
   Ar  6.29501012040417507 -0.769712471349173866 -6.91454332254278281
   Ar  3.49995546789933476 -8.00704920137973453 -0.426526631939732892
   Ar  0.385154812289867643 17.8769740351009112 -17.4065226240143041
   Ar  21.2288869131365736 10.2327102035561044 -13.0872200859088803
   Ar  1.22082587001210396 5.83597435065779457 16.8450099266840283
   Ar  -7.08754036219628425 6.03412971863339109 -22.3251445579668015
   Ar  -0.244265849036998037E-01 17.4693605251376454 7.37116730966604194
   Ar  15.0981822679441553 9.88940516251130397 -8.49382740142986670
   Ar  -6.57877688336587152 -15.0484532074656290 14.7230359830473887
   Ar  -2.22666666633409394 -4.18421900331013674 -2.47007887105670587
   Ar  5.20621069851729867 -22.6565181989138011 7.39475674805799521
   Ar  -8.85828800414884299 -2.47510661993999781 2.35441398531938617
   Ar  6.75202354538700167 0.430391383628436597 5.43492495261394382
   Ar  11.9263127546080856 8.13267254152258445 2.40081132956567966
   Ar  -14.5507562394484040 -0.471540677239574602 -13.7058431104765983
   Ar  14.1157692422228553 -2.98968593175088149 24.6842798176059546
   Ar  -3.35107336204723527 -0.681362546744063047 -7.37039916831594510
   Ar  7.79269876443546838 3.30687615091469800 -0.732378021069576002E-01
   Ar  -1.13289059102623746 -17.1672835835708497 -12.9126466371968966
   Ar  -9.21054349522787241 -10.4846510042527843 -8.38485797788161591
   Ar  -6.47848777956778044 -3.90736653076878993 -10.6499668409808841
   Ar  0.987874979233200667 13.7363585340729077 5.07209659800543733
   Ar  8.86097814789463278 9.96103887786039799 1.09373795795780060
   Ar  -6.58068766844202013 -20.4019345282015756 -7.28935608176262662
   Ar  -0.448977062720621045 19.5862520159664086 11.0351198968750293
   Ar  7.36056937465398153 -2.69594281683156067 7.26081874603436361
   Ar  13.8791344546872004 12.1903465249438128 1.24889885444881155
   Ar  -3.65782753722175302 19.7829061761924159 -11.8161510229542408
   Ar  1.49729450944005649 -5.39289977250827679 1.92445849672255198
   Ar  18.5861605633917577 3.00868366398259690 2.06440131010935168
   Ar  6.20730767975507014 -9.47418398815358032 5.54930507752316249
   Ar  3.65054837888884753 3.43181054126032858 -4.31160813615129435
   Ar  2.67862616463048520 2.29300605146530545 5.98502962150055051
   Ar 24.5113122275150914 4.00733170976478448 13.1412501215423774
   Ar -0.600233262008137092 3.62825631372324597 6.38411284716526772
    &END COORD
    &TOPOLOGY
    &END TOPOLOGY
  &END SUBSYS
&END FORCE_EVAL

TASK
  • Run the calculation and visualize the trajectories using VMD. Plot the total energies (conserved quantity), temperature, potential energies (you need to look into the *.ener file). Comment and explain what you observe.
  • Change the timestep to 40 fs and rerun, look into the ener file again and compare to the simulation with the smaller timestep. What happens to the energy conservation?
  • Change the temperature to 150 K (use the smaller timestep of 5 fs) and compare the temperatures that you set up in the input files with the average temperature. Are they the same? Comment.
See exercise 0 to remember the commands to define the unit cell in VMD and how to wrap the coordinates to the original cell.

Part III: Radial distribution functions in NVT ensemble

In this exercise you are asked to compute the radial distribution function of liquid Ar at different temperatures. First of all perform two simulations at 85 K and 150 K for liquid Ar in the NVT ensemble to ensure the simulations are equilibrated at the right temperatures. To perform simulations in NVT copy the relevant section in the input file as shown below.

  &MD
    ENSEMBLE NVT
    STEPS 10000
    TIMESTEP 5
    TEMPERATURE 85.0
    &THERMOSTAT
      &NOSE                    #Uses the Nose-Hoover thermostat
        TIMECON 100           #timeconstant of the thermostat chain, how often does thermostat adjust your system 
      &END NOSE
    &END
  &END MD

Use VMD or write your own program (Fortran, C, C++, Python etc.) to calculate radial distribution $g(r)$. Plot $g(r)$, and against various the temperatures to examine the effects. VMD comes with an extension for exactly this purpose: In the VMD Main window open “Extensions → Analysis” click on “Radial Pair Distribution function $g(r)$. In the appearing window use “Utilities → Set unit cell dimensions” to let VMD know the simulation box you used. After that use Selection 1 and 2 to define the atomic types that you want to calculate the rdf for, for example “element Ar”. In the plot window, use “File”, you can save the plot data.

TASK
  • Plot $g(r)$ at 85 and 150 K into the same graph.
  • What are the differences in the height of the first peak, and why does temperature contribute to the differences?
  • Compared to experimental data exp_gr.dat taken at 84 K, what does this say about the structure of the liquid and is this expected? .
exp_gr.dat
.0681      -.0789                                              
.1362      -.0545                                              
.2043      -.0263                                              
.2724      -.0062                                              
.3405      -.0008                                              
.4086      -.0088                                              
.4767      -.0222                                              
.5448      -.0313                                              
.6129      -.0305                                              
.6810      -.0206                                              
.7491      -.0082                                              
.8172      -.0010                                              
.8853      -.0031                                              
.9534      -.0123                                              
1.0215     -.0222                                              
1.0896     -.0259                                              
1.1577     -.0203                                              
1.2258     -.0085                                              
1.2939     .0026                                               
1.3620     .0065                                               
1.4301     .0011                                               
1.4982     -.0095                                              
1.5663     -.0180                                              
1.6344     -.0180                                              
1.7025     -.0086                                              
1.7706     .0053                                               
1.8387     .0155                                               
1.9068     .0157                                               
1.9749     .0058                                               
2.0430     -.0076                                              
2.1111     -.0153                                              
2.1792     .0111                                               
2.2473     .0038                                               
2.3154     .0210                                               
2.3835     .0295                                               
2.4516     .0232                                               
2.5197     .0052                                               
2.5878     -.0128                                              
2.6559     -.0176                                              
2.7240     .0030                                                
2.7921     .0241                                               
2.8602     .0455                                               
2.9283     .0427                                               
2.9964     .0112                                               
3.0645     -.0279                                              
3.1326     -.0276                                              
3.2007     .0721                                               
3.2688     .3212                                               
3.3369     .7356                                               
3.4050     1.2831                                              
3.4731     1.8857                                              
3.5412     2.4408                                              
3.6093     2.8510                                              
3.6774     3.0542                                              
3.7455     3.0403                                              
3.8136     2.8497                                              
3.8817     2.5549                                              
3.9498     2.2343                                              
4.0179     1.9474                                              
4.0860     1.7218                                              
4.1541     1.5545                                              
4.2222     1.4240                                              
4.2903     1.3059                                              
4.3584    1.1856                                              
4.4265     1.0624                                              
4.4946     .9458                                               
4.5627     .8471                                               
4.6308     .7730                                               
4.6989     .7219                                               
4.7670     .6862                                               
4.8351     .6571                                               
4.9032     .6293                                               
4.9713     .6023                                               
5.0394     .5796                                               
5.1075     .5651                                               
5.1756     .5604                                               
5.2437     .5640                                               
5.3118     .5725                                               
5.3799     .5829                                               
5.4480     .5945                                               
5.5161     .6091                                               
5.5842     .6294                                               
5.6523     .6569                                               
5.7204     .6911                                               
5.7885     .7291                                                
5.8566     .7675                                               
5.9247     .8043                                               
5.9928     .8395                                               
6.0609     .8751                                               
6.1290     .9137                                               
6.1971     .9569                                               
6.2652     1.0038                                              
6.3333     1.0520                                              
6.4014     1.0983                                              
6.4695     1.1398                                              
6.5376     1.1753                                              
6.6057     1.2048                                              
6.6738     1.2287                                              
6.7419     1.2475                                              
6.8100     1.2608                                              
6.8781     1.2685                                              
6.9462     1.2706                                              
7.0143     1.2679                                              
7.0824     1.2618                                              
7.1505     1.2536                                              
7.2186     1.2435                                              
7.2867     1.2304                                              
7.3548     1.2128                                              
7.4229     1.1891                                              
7.4910     1.1594                                              
7.5591     1.1254                                              
7.6272     1.0899                                              
7.6953     1.0561                                              
7.7634     1.0255                                              
7.8315     .9982                                               
7.8996     .9726                                               
7.9677     .9470                                               
8.0358     .9206                                               
8.1039     .8942                                               
8.1720     .8698                                               
8.2401     .8499                                               
8.3082     .8360                                               
8.3763     .8280                                               
8.4444     .8247                                               
8.5125     .8244                                               
8.5806     .8259                                               
8.6487     .8293                                               
8.7168     .8355                                               
8.7849     .8452                                               
8.8530     .8590                                               
8.9211     .8757                                               
8.9892     .8937                                               
9.0573     .9113                                               
9.1254     .9276                                               
9.1935     .9431                                               
9.2616     .9587                                               
9.3297     .9757                                               
9.3978     .9943                                               
9.4659     1.0138                                              
9.5340     1.0326                                              
9.6021     1.0492                                              
9.6702     1.0628                                              
9.7383     1.0738                                              
9.8064     1.0830                                              
9.8745     1.0913                                              
9.9426     1.0988                                              
10.0107    1.1048                                              
10.0788    1.1082                                              
10.1469    1.1079                                              
10.2150    1.1041                                              
10.2831    1.0977                                              
10.3512    1~0904                                              
10.4193    1.0835                                              
10.4874    1.0774                                              
10.5558    1.0716                                              
10.6236    1.0648                                              
10.6917    1.0559                                              
10.7598    1.0448                                              
10.8279    1.0322                                              
10.8960    1.0196                                              
10.9641    1.0083                                              

Part IV: Ensembles

In previous section, you have already run NVE and NVT ensemble molecular dynamics for Ar liquid. In this section, we will focus on the NPT ensembles and you will compare the results in different ensembles.

Set up NPT calculation, change the setting in &MD section.

&FORCE_EVAL
...
STRESS_TENSOR ANALYTICAL
...
&END FORCE_EVAL
...
...
&MD
  ENSEMBLE NPT_I                #constant temperature and pressure using an isotropic cell
  STEPS 10000
  TIMESTEP 5.
  TEMPERATURE 85.0
  &BAROSTAT
    PRESSURE 1.                 # PRESSURE, unit[bar]
    TIMECON 1000
  &END BAROSTAT
  &THERMOSTAT
    &NOSE
      TIMECON 1000
    &END NOSE
 &END MD
TASK
  • For the calculations in the NVT ensemble at 150 and 85 K, check the temperature, and energy of the whole system, and compare the result to the previous simulations run in NVE, and rationalize the difference.
  • It is a common practice to first perform a simulation in NVT and then run an NVE simulation. What is a possible reason for doing this?
  • It is often needed to perform MD simulations in the NPT ensemble. For the case of Argon, it is liquid at 85 K at atmospheric pressure. First perform an NPT simulation at 85 K and atmospheric pressure. Then, based on the phase diagram reported in this link, choose a possible value of the pressure for liquid Argon at 150 K, edit the input file accordingly and run the simulation in NPT.
  • For these two simulations plot the time evolution of the cell volume as well as the pressure. To print the cell to file you need to add the keyword &CELL in the print section. For the pressure you can use the following command from the terminal grep 'PRESSURE' NAME_OF_OUTPUT_FILE.out | awk '{i++; print i, $4}' Plot and comment the results.
 &END MD
...   
  &PRINT
   &CELL
     &EACH
       MD 1
     &END
   &END
  &END
  • Based on the differences between the results obtained from the simulations carried out in the different ensambles and on what you have learned during the lectures, explain some of the possible reasons to perform simulations in a given ensamble instead of another.