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)$ and investigate a variety of ensembles.
You are expected to hand in the short report via OLAT, 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 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, learn how to organize and edit 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 L-J potential energy expression in use 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, $\sigma$ is related to the minimum of the Lennard-Jones potential, and $r_{ij}$ is the distance between atoms i and j.
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 for $r_{ij}$ of the L-J potential. For realistic results, the cut-off should be less than half of the simulation box size and larger than $\sigma$. The radial distribution function, (or pair correlation function) $g(r)$, in a system of particles (atoms, molecules, colloids, etc.), describes how the density varies as a function of distance from a reference particle.
Part I: Set up MD simulation
In this section, we provide you with an example CP2K input for an MD calculation. Extensive comments have been added to the file, which start with a has symbol '#'.
1. Step
Load the CP2K module as explained in Exercise 0, 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):
- argon.inp
## It's highly recommended to go to ## https://manual.cp2k.org/ ## and learn how to set up a CP2K ## calculation correctly, using the 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
& 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
$ cp2k.popt -i argon.inp | tee argon.out
which will write the output simultaneously to a file argon.out and show it in the terminal.
- 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 steps (100, 1000, 5000), and inspect the geometry in each case. Plot the total energy, temperature, and potential energy. Try to comment and explain what you observe.
Part II: Force Field Parameter
In this section we investigate the dependence of the L-J potential on the its different parameters. To do se, we investigate a small system of only two Ar atoms. The following code snippet shows the changes that you need to make to the input file in order to run such a calculation.
&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
In order to investigate the effect of the force-field parameters on the L-J potential you need to vary multiple parameters of your calculation:
- you need to actually vary the L-J force-field parameters, $\epsilon$ and $\sigma$
- you need to run each calculation at different distances between the two Ar atoms
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|" energy.inp > energy_${d}A.inp
cp2k.popt -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 4generates the numbers2.0,2.1,2.2, … ,4.0(try it out!) - With
for d in $(seq 2 0.1 4); dowe 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.inplooks for4 0 0in the fileenergy.inp(the original file from above) and replaces4 0 0by$d 0 0(that is:2.0,2.1,2.2, …)- … and using
> energy_${d}A.outwe redirect the output of thesedcommand to new filesenergy_2.0A.out,energy_2.1A.out, etc. - Then we run
cp2k.poptas shown in Exercise 0 on those new input files and write the output to new output files as well - Using
awkwe extract the energy from the output file and print it
- Plot the Lennard-Jones potential against the Ar-Ar distance $r$ (2-4 Å) using different values for $\epsilon$ and $\sigma$.
- Repeat the L-J MD calculation with different $\epsilon$ and $\sigma$, compare the potential energy and temperature evolution, and explain the relation between the calculated properties and force field parameters.
Part III: Radial distribution function
In this section we analyze the dependence of the radial distribution function (rdf), $g(r)$, on the temperature of the system. To do so, you should plot $g(r)$ against various temperatures and examine the effects. You can use VMD (as explained below) or write your own program (Fortran, C, C++, Python etc.) to calculate the rdf.
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 tell VMD the size of 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 the “File” menu to save the plot data.
- 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/how does the temperature contribute to the differences?
- Compare your results to the experimental data taken at 84 K (given in
exp_gr.dat, below). What does this say about the structure of the liquid and did you expect this?
- 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: Other Ensembles
In the previous sections, you have already run NVE ensemble molecular dynamics calculations for liquid Ar. In this section, we will focus on the NVT and NPT ensembles.
To set up an NVT calculation, change the settings in the &MD section as shown below:
&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
To set up an NPT calculation, change the settings in the &MD section as shown below:
&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
- Run a calculation using the NVT ensemble at 300K. Check the temperature and energy of the whole system, and compare the result to an NVE ensemble (300K). Rationalize and discuss the difference.
- Run a calculation using the NVT ensemble (300K) until the system is equilibrated, then run an NVE ensemble calculation. Check the temperature and energy of the whole system, and compare to the previous NVE simulation. Discuss your observations.
- Remove some atoms and run an NPT ensemble simulation. Then, check the size of the simulation box. Discuss your observations.
For the purposes of this example, you should see a file called ar108-1.restart (at least for Part 1 of this exercise). Note: there will also be multiple backup files (.bak-# suffixes) which you do not need to care about.
These files are nothing but another input file. However, their parameters are set such that they continue a CP2K calculation from the last step of the simulation which generated the RESTART file.
Here are two options for how you can use these RESTART-files:
1. Directly using the RESTART as an input.
- You can copy the RESTART file to a new input file:
cp ar108-1.restart argon_follow_up.inp
- Now you can change the input to your liking (e.g. change the ensemble, etc.)
- And finally simply run CP2K with the new input file:
cp2k -i argon_followup.inp -o argon_followup.out
2. You can also tell CP2K to load a specific RESTART-file.
- Write a new input file as usual:
argon_followup.inp
- Add an EXT_RESTART section:
&EXT_RESTART RESTART_FILE_NAME ar108-1.restart &END EXT_RESTART
- And now, again, simply run CP2K:
cp2k -i argon_followup.inp -o argon_followup.out