User Tools

Site Tools


exercises:2017_uzh_cmest:first_simulation_run

Differences

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

Link to this comparison view

exercises:2017_uzh_cmest:first_simulation_run [2017/09/20 07:31] (current)
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
 +</​code>​
 +
 +Save the following input to a file named ''​energy.inp''​ (for example using ''​$ vim energy.inp''​):​
 +     
 +<code - 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
 +</​code>​
 +
 +<note tip>​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.</​note>​
 +
 +=== 2. Step ===
 +
 +Run a CP2K calculation with the following command:
 +
 +<​code>​
 +$ cp2k.sopt -i energy.inp -o energy.out
 +</​code>​
 +
 +<note tip>​Alternatively you can run it using <​code>​$ cp2k.sopt -i energy.inp | tee energy.out</​code>​ which will write the output simultaneously to a file ''​energy.out''​ and show it in the terminal.</​note>​
 +
 +=== 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
 +</​code>​
 +
 +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
 +</​code>​
 +
 +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}} $ .
 +
 +<note warning>​Always check the number of warnings by looking at the line: <​code>​The number of warnings for this run is : ...</​code> ​ 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.</​note>​
 +===== 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.
 +
 +<note important>​
 +The output file will be overwritten every time you run a calculation,​ unless you change its name.
 +</​note>​
 +
 +To do so: 
 +<​code>​
 +$ mv energy.out energy_dist4A.out
 +</​code>​
 +
 +<note tip>
 +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. 
 +</​note>​
 +
 +For doing so:
 +<​code>​
 +$ cp energy.inp energy_dist2A.inp ​
 +</​code>​
 +
 +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
 +</​code>​
 +
 +=== 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
 +</​code>​
 +
 +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 ​
 +</​code>​
 +
 +  * The ''​CHARGE''​ section must be also duplicated: ​
 +
 +<​code>​
 +         &​CHARGE
 +           ATOM Ne
 +           ​CHARGE 0.0
 +         &​END CHARGE
 +         &​CHARGE
 +           ATOM Kr
 +           ​CHARGE 0.0
 +         &​END CHARGE
 +         
 +</​code>​
 +
 +  * 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
 +</​code>​which gives you:<​code>​
 + ​ENERGY| Total FORCE_EVAL ( FIST ) energy (a.u.): ​            ​-0.000250281091139
 +</​code>​
 +  * Using the ''​awk''​ tool:<​code>​
 +$ awk '/​Total FORCE_EVAL/ { print $9; }' energy.out
 +</​code>​which returns:<​code>​
 +-0.000250281091139
 +</​code>''​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
 +</​code>​
 +
 +  * 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
exercises/2017_uzh_cmest/first_simulation_run.txt · Last modified: 2017/09/20 07:31 by tmueller