Run your first simulation using CP2K
Computation of the Lennard Jones curve
- Part I: Single Point (Energy) calculation
- Part II: Computation of the LJ energy curve
Tips & Tricks
- Parsing the output
- Generating input files

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 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:

mkdir ex1
cd ex1

Save the following input to a file named energy.inp (for example using $ vim energy.inp):

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:

$ cp2k.sopt -i energy.inp -o energy.out

Alternatively you can run it using

$ 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:

[...]
  **** **** ******  **  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:

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:

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:

$ 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:

$ 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:

$ 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:

         &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:

      &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:

         &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:

$ grep "Total FORCE_EVAL" energy.out

which gives you:

 ENERGY| Total FORCE_EVAL ( FIST ) energy (a.u.):             -0.000250281091139

Using the awk tool:
```
$ awk '/Total FORCE_EVAL/ { print $9; }' energy.out
```
which returns:
```
-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):

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

Table of Contents

Run your first simulation using CP2K

Computation of the Lennard Jones curve

Part I: Single Point (Energy) calculation

1. Step

2. Step

3. Step

Part II: Computation of the LJ energy curve

1. Step

2. Step

3. Step

4. Step

Tips & Tricks

Parsing the output

Generating input files