# Open SourceMolecular Dynamics

### Site Tools

exercises:2017_uzh_cmest:first_simulation_run

# Differences

This shows you the differences between two versions of the page.

 — exercises:2017_uzh_cmest:first_simulation_run [2017/09/20 07:31] (current)tmueller created 2017/09/20 07:31 tmueller created 2017/09/20 07:31 tmueller created Line 1: Line 1: + ====== 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 quiet 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. + + ====== Computation of the Lennard Jones curve ====== + + In this exercise you will compute the Lennard-Jones energy curve for a system of two Krypton (Kr) atoms using a Molecular Mechanics simulation rather than the analytical form of the potential. + + In Part I you find the instructions for computing the energy of two Kr atoms at a distance $r=4.00 Å$. + + In Part II you find the instructions for getting the energy profile as a function of $r$. + + Additonal parameters for Neon (Ne) and combination rules to obtain new parameters are provided in  Part III and IV. + + You are expected to hand in the respective plots by email, either as one PDF or as one file per plot (PNG, JPEG, or PDF format). + ===== Part I:  Single Point (Energy) calculation ===== + + In this section a commented CP2K input example for a single point calculation is provided. + Comments are added, they start with an exclamation mark '​!'​. + + === 1. Step === + + Load the CP2K module as shown before, create a directory ''​ex1''​ and change to it: + + <​code>​ + mkdir ex1 + cd ex1 + ​ + + Save the following input to a file named ''​energy.inp''​ (for example using ''​$vim energy.inp''​):​ + + ​ + &​GLOBAL ​ ! section to select the kind of calculation + ​RUN_TYPE ENERGY ​ ! select type of calculation. In this case: ENERGY (=Single point calculation) + &END GLOBAL + &​FORCE_EVAL ​ ! section with parameters and system description + METHOD FIST ! Molecular Mechanics method + &​MM ​ ! specification of MM parameters ​ + &​FORCEFIELD ​ ! parameters needed to describe the potential ​ + &SPLINE + EMAX_SPLINE 10000 ! numeric parameter to ensure calculation stability. Should not be changed + &END SPLINE + &​NONBONDED ​ ! parameters for the non bonded interactions + &​LENNARD-JONES ! Lennard-Jones parameters + ATOMS Kr Kr + EPSILON ​ [K_e] 164.56 + SIGMA [angstrom] ​ 3.601 + RCUT [angstrom] ​ 25.0 + &END LENNARD-JONES + &END NONBONDED + &CHARGE + ATOM Kr + CHARGE 0.0 + &END CHARGE + &END FORCEFIELD + &​POISSON ​ ! solver for non periodic calculations + PERIODIC NONE + &EWALD + EWALD_TYPE none + &END EWALD + &END POISSON + &END MM + &​SUBSYS ​ ! system description ​ + &CELL + ABC [angstrom] 10 10 10 ​ + PERIODIC NONE + &END CELL + &​COORD ​ ​ + UNIT angstrom + Kr 0 0 0 + Kr 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.​ + + === 2. Step === + + Run a CP2K calculation with the following command: + + <​code>​ +$ cp2k.sopt -i energy.inp -o energy.out + ​ + + ​Alternatively you can run it using <​code>​$cp2k.sopt -i energy.inp | tee energy.out​ which will write the output simultaneously to a file ''​energy.out''​ and show it in the terminal.​ + + === 3. Step === + + The output (''​$ less energy.out''​) should look like this: + + <​code>​ + [...] + **** **** ****** ​ **  PROGRAM STARTED AT               ​2016-09-22 15:​15:​15.977 + ***** ** ***  *** **   ​PROGRAM STARTED ON                                tcopt3 + ​** ​   ****   ​****** ​   PROGRAM STARTED BY                             ​studentXX + ***** **    ** ** **   ​PROGRAM PROCESS ID                                112277 + **** **  ******* ​ **  PROGRAM STARTED IN          /​data/​students/​studentXX/​ex1 + [...] + ENERGY| Total FORCE_EVAL ( FIST ) energy (a.u.): ​         -0.000518941408898 + [...] + + The number of warnings for this run is : 0 + + ​------------------------------------------------------------------------------- + **** **** ****** ​ **  PROGRAM ENDED AT                 ​2016-09-22 15:​15:​16.027 + ***** ** ***  *** **   ​PROGRAM RAN ON                                    tcopt3 + ​** ​   ****   ​****** ​   PROGRAM RAN BY                                 ​studentXX + ***** **    ** ** **   ​PROGRAM PROCESS ID                                112277 + **** **  ******* ​ **  PROGRAM STOPPED IN          /​data/​students/​studentXX/​ex1 + ​ + + If you get the closing banner you know that CP2K finished. + + The following line tells you the result: + <​code>​ + ENERGY| Total FORCE_EVAL ( FIST ) energy (a.u.): ​            ​-0.000518941408898 + ​ + + This is the energy (in Hartree) for a system of 2 Kr atoms at distance $r=4.00 Å$ + + Note, that in the input-file ''​EPSILON''​ is given in units of //Kelvin//, whereas in the output the energy is printed in //​Hartree//,​ which is the unit of energy in the system of atomic units (a.u.). + + To convert from //Kelvin// to //Hartree// you have to multiply with the Boltzmann constant $k_\text{b} = 3.1668154 \cdot 10^{-6} \frac{E_\text{H}}{\text{K}}$ . + + ​Always check the number of warnings by looking at the line: <​code>​The number of warnings for this run is : ... ​ If that number is not zero you **must** check the rest of the output for warnings and act accordingly,​ otherwise you may work with wrong results.​ + ===== Part II: Computation of the LJ energy curve ===== + + In this section a few scripts to get the LJ energy profiles are presented. + + === 1. Step === + + In order to get a good profile, a set of energy values as a function of the interatomic distance is needed. You can use the ''​energy.inp''​ input file and change the Kr coordinates in order to get different starting distances. + + ​ + The output file will be overwritten every time you run a calculation,​ unless you change its name. + ​ + + To do so: + <​code>​ + $mv energy.out energy_dist4A.out + ​ + + + If you run multiple calculations,​ it is always good to keep track of what you have done by producing an input and an output for every distance you are planning to run. + ​ + + For doing so: + <​code>​ +$ cp energy.inp energy_dist2A.inp ​ + ​ + + then edit the input file to update to the new coordinates (e.g. 2 Å) and rerun CP2K to produce a new output file: + + <​code>​ + $cp2k.sopt -i energy_dist2A.inp -o energy_dist2A.out + ​ + + === 2. Step === + + When you have tested a few distances, you can produce a table looking like: + + ^ Input file ^ Distance (Å) ^ Energy ​ ​(Eh) ​ ^ + | energy_dist1A.inp ​ | 1 | ... | + | energy_dist2A.inp ​ | 2 | ... | + | energy_dist3A.inp ​ | 3 | ... | + | ... | ... | ... | + + This is the Lennard Jones energy curve for two Kr atoms. + + By using any plotting program you can now get a representation of the energy profile. + Choose a an appropriate minimum distance and step size. + + === 3. Step === + + Now we do the same for Ne atoms: use the previous input file as a template and update the parameters according to the following: + + <​code>​ + &​NONBONDED ​ + &​LENNARD-JONES ! Lennard-Jones Ne parameters + ATOMS Ne Ne + ​EPSILON ​ [K_e] 36.831 ​ + SIGMA [angstrom] ​ 2.775 + ​RCUT ​ [angstrom] 25.0 + &​END LENNARD-JONES + &​END NONBONDED + &​CHARGE + ATOM Ne + ​CHARGE 0.0 + &​END CHARGE + ​ + + Plot the energy curve again. + + === 4. Step === + + Finally we look at the curve for Kr-Ne. + + The epsilon and sigma for the Lennard-Jones potential between Kr and Ne can be calculated using the parameters from the Kr-Kr and Ne-Ne interaction:​ + + $$\sigma_{ij}= \sqrt{\sigma_i\sigma_j}$$ \\ + $$\epsilon_{ij}= \sqrt{\epsilon_i\epsilon_j}$$ + + Please note: + + * The ''​LENNARD-JONES''​ section must be present for all the three possible couples now: Kr-Kr, Ne-Ne and Ne-Kr: ​ + + <​code>​ + &​LENNARD-JONES ! Lennard-Jones parameters for Ar-Ar interaction + ATOMS Kr Kr + EPSILON ​ [K_e] 164.56 + SIGMA [angstrom] ​ 3.601 + RCUT [angstrom] ​ 25.0 + &END LENNARD-JONES + &​LENNARD-JONES ! Lennard-Jones Ne-Ne parameters + ATOMS Ne Ne + EPSILON ​ [K_e] 36.831 ​ + SIGMA [angstrom] ​ 2.775 + RCUT [angstrom] 25.0 + &END LENNARD-JONES + &​LENNARD-JONES ! Lennard-Jones parameters for Kr-Ne interaction + ATOMS Kr Ne + EPSILON ​ [K_e] YOUR EPSILON + SIGMA [angstrom] ​ YOUR SIGMA + RCUT [angstrom] ​ 25.0 + &END LENNARD-JONES ​ + ​ + + * The ''​CHARGE''​ section must be also duplicated: ​ + + <​code>​ + &​CHARGE + ATOM Ne + ​CHARGE 0.0 + &​END CHARGE + &​CHARGE + ATOM Kr + ​CHARGE 0.0 + &​END CHARGE + + ​ + + * one of the atom kinds in the ''&​COORD''​ section must be set to Kr, the other to Ne. + + Plot again the energy curve. + + ====== Tips & Tricks ====== + + ===== Parsing the output ===== + + Many times you will have to get some value out of a simulation output, in this case, the energy. + This can achieved in a number of ways: + + * Using the ''​grep''​ command:<​code>​ +$ grep "Total FORCE_EVAL"​ energy.out + ​which gives you:<​code>​ + ​ENERGY| Total FORCE_EVAL ( FIST ) energy (a.u.): ​            ​-0.000250281091139 + ​ + * Using the ''​awk''​ tool:<​code>​ + $awk '/​Total FORCE_EVAL/ { print$9; }' energy.out + ​which returns:<​code>​ + -0.000250281091139 + ''​awk''​ reads a given file line-by-line and splits a line into multiple fields (using whitespace as delimiter). The command above tells ''​awk''​ to look for the string ''​Total FORCE_EVAL''​ in a line and if found, print field number 9 of it, effectively returning the energy. + + ===== Generating input files ===== + + Many 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): + + <​code>​ + 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.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