====== 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 short report by email, ONLY in 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 [[https://manual.cp2k.org/|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 (L-J) fluid containing mono-atomic
particles. The [[https://en.wikipedia.org/wiki/Lennard-Jones_potential | L-J 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 is $\epsilon$ the well depth, σ is related to the minimum of Lennard-Jones potential, rij is the
distance between atoms i and j.
[[https://en.wikipedia.org/wiki/Periodic_boundary_conditions|Periodic boundary conditions (PBC)]] 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: Set up MD simulation =====
In this section, a commented CP2K input example for an MD calculation is provided.
Comments are added, they start with a hash mark '#'.
=== 1. Step ===
Load the CP2K module as shown before, create a directory ''ex1'' and change to it:
mkdir ex1
cd ex1
Save the following input to a file named ''argon.inp'' (for example using ''$ vim 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 100 #The number of MD steps to perform
TIMESTEP 5. #The length of an integration step (fs)
TEMPERATURE 85.0 #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
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.
=== 2. Step ===
Run a CP2K calculation with the following command:
$ cp2k.popt -i argon.inp -o argon.out
Alternatively you can run it using $ cp2k.sopt -i argon.inp | tee argon.out which will write the output simultaneously to a file ''argon.out'' and show it in the terminal.
**TASK**
*Run the calculation and visualize the trajectories using VMD
*Run calculations with different timesteps (0.5 2 5fs), different temperatures(84, 300, 400K), different densities($\rho$=0.25,0.5,1), and a different number of timesteps (100,1000,5000), and inspect geometry in each case, plot the total energies, temperature, potential energies. Try to comment and explain what you observe.
===== Part II: Force Field Parameter =====
You need to modify the L-J force field parameter. Run calculations with different $\sigma$ and $\epsilon$ and inspect the insight of L-J potential (Energies, velocities, temperatures etc ). Then use the force field parameter provided to plot the L-J potential curve.
&GLOBAL ! section to select the kind of calculation
RUN_TYPE ENERGY ! select type of calculation. In this case: ENERGY (=Single point calculation)
&END GLOBAL
...
&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
...
&SUBSYS ! system description
&CELL
ABC [angstrom] 10 10 10
PERIODIC NONE
&END CELL
&COORD
UNIT angstrom
Ar 0 0 0
Ar 4 0 0
&END COORD
&END SUBSYS
any times you will have to run the same simulation with different parameters (here the distance).
A simple way to generate the different input files is using shell scripting in combination with ''sed'' (the stream editor):
for d in $(seq 2 0.1 4); 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 $9; }' energy_${d}A.out
done
* The command ''seq 2 0.1 4'' generates the numbers ''2.0'', ''2.1'', ''2.2'', ... , ''4.0'' (try it out!)
* With ''for d in $(seq 2 0.1 4); do'' we use the shell to run all commands which follow once for every number (stored in ''$d'')
* ''sed -e "s|4 0 0|$d 0 0|" energy.inp'' looks for ''4 0 0'' in the file ''energy.inp'' (the original file from above) and replaces ''4 0 0'' by ''$d 0 0'' (that is: ''2.0'', ''2.1'', ''2.2'', ...)
* ... and using ''> energy_${d}A.out'' we redirect the output of the ''sed'' command to new files ''energy_2.0A.out'', ''energy_2.1A.out'', etc.
* Then we run ''cp2k.sopt'' as shown before on those new input files and write the output to new output files as well
* Using ''awk'' we extract the energy from the output file
**TASK**
*Plot the Lennard-Jones potential against Ar-Ar distance $r$ (2-4 Å) using different $\epsilon$ and $\sigma$.
*Repeat the L-J MD calculation with different $\epsilon$ and $\sigma$, and compare the potential energy, temperature evolution, and explain the relation between calculated properties and force field parameters.
===== Part III: Radial distribution functions =====
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 84, 300 and 400 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? .
.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 ensemble molecular dynamics for Ar liquid. In this section, we will focus on the NVT, NPT ensembles.
Step up NVT calculation, change the setting in &MD section.
&MD
ENSEMBLE NVT
STEPS 3000
TIMESTEP 5
TEMPERATURE 298
&THERMOSTAT
REGION MASSIVE
&NOSE #Uses the Nose-Hoover thermostat
TIMECON 1000 #timeconstant of the thermostat chain, how often does thermostat adjust your system
&END NOSE
&END
&END MD
Step 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 3000
TIMESTEP 5.
TEMPERATURE 85.0
&BAROSTAT
PRESSURE 0. # PRESSURE, unit[bar]
TIMECON 1000
&END BAROSTAT
&THERMOSTAT
&NOSE
TIMECON 1000
&END NOSE
&END THERMOSTAT
&END MD
**TASK**
*Run calculation using NVT at 300K, check the temperature, and energy of the whole system, and compare the result to NVE (300K), and rationalize the difference.
*Run calculation using NVT (300K) until the system is equilibrated then run NVE, check the temperature, and energy of the whole system, and compare to previous NVE simulation.
*Remove some atoms and run NPT simulation then check the size of the simulation box.