====== Calculations with Periodic Boundary Conditions ======

This short exercise illustrates the property of periodic boundary conditions and at the same time gives you some ideas/skeletons on how to create useful bash scripts. The files you need for this are:

A geometry file for a <chem>H2O</chem> molecule:
<code - h2o.xyz>
3
Atom    X   Y   Z
O   0.0000000   0.0000000   0.1194180
H   0.0000000   0.7654990   -0.4776700
H   0.0000000   -0.7654990  -0.4776700
</code>

A CP2K input file to run a calculation with it:
<code - h2o_pbc.inp>
&GLOBAL
  PROJECT h2o_pbc
  RUN_TYPE ENERGY 
  PRINT_LEVEL MEDIUM
&END GLOBAL

&FORCE_EVAL
  METHOD Quickstep              ! Electronic structure method (DFT,...)
  &DFT
    BASIS_SET_FILE_NAME  BASIS_MOLOPT
    POTENTIAL_FILE_NAME  POTENTIAL

    &POISSON                    ! Solver requested for non periodic calculations
      PERIODIC XYZ
    &END POISSON
    &SCF                        ! Parameters controlling the convergence of the SCF. This section should not be changed. 
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-6
      MAX_SCF 300
    &END SCF
    &XC                         ! Parameters needed to compute the electronic exchange potential 
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT

  &SUBSYS
    &CELL
      ABC 10. 10. 10.
      PERIODIC XYZ
    &END CELL
    &TOPOLOGY                    ! Section used to center the atomic coordinates in the given box. Useful for big molecules
      COORD_FILE_FORMAT xyz
      COORD_FILE_NAME  ./h2o.xyz
      &CENTER_COORDINATES
        CENTER_POINT 5. 5. 5.
      &END
    &END
    &KIND H
      ELEMENT H
      BASIS_SET TZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE-q1
    &END KIND
    &KIND O
      ELEMENT O
      BASIS_SET TZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE-q6
    &END KIND
  &END SUBSYS
&END FORCE_EVAL
</code>

And a script to run it:

<code - run.sh>
#!/bin/bash

set -o errexit
set -o nounset
set -o pipefail

xseq=$(seq 5 0.1 6)

for x in $xseq; do
    rm -f "h2o_pbc_x-${x}.out"

    sed \
        -e "s|h2o_pbc|h2o_pbc_x-$x|" \
        -e "s|CENTER_POINT .*|CENTER_POINT $x 5. 5.|" \
        h2o_pbc.inp > "h2o_pbc_x-${x}.inp"

    cp2k.sopt -i "h2o_pbc_x-${x}.inp" -o "h2o_pbc_x-${x}.out" &
done

for job in $(jobs -p) ; do
    wait $job
done

# create or truncate the file named "energies"
:>| energies

for x in $xseq; do
    energy=$(awk '/Total FORCE_EVAL/ { print $9 }' "h2o_pbc_x-${x}.out")
    echo "$x $energy" >> energies
done
</code>

You have to mark the script as an executable after creating it, using:

<code>
$ chmod +x run.sh
</code>

===== Create a plot =====

Run the script after loading the CP2K module:

<code>
$ ./run.sh
</code>

Besides various input and output files, you should get a file named ''energies''.

  * Plot it, including meaningful axis labels. To that end, try to understand the provided script. Looking at the difference between the initially provided input file ''h2o_pbc.inp'' and the newly created input files may give a hint. The command ''$ diff h2o_pbc.inp SOMEOTHERFILE'' may also help.
  * What would you expect from such a plot and why? What order is the deviation from what you would expect? What could be the source for this deviation?

===== Make the script your own =====

  * Read the script and try to understand the commands. Look them up using your favorite search engine, the man pages on the system (''$ man sed'' for example) or try them out by hand on the command line.
  * Try to document the script in such a way that you will still understand it when looking at it again in a year or two (so, not every detail, but make sure to cover the obscure stuff and the big blocks)
  * Change the script to calculate energies for ''x=z=5.'' and ''y=2,2.2,2.4,..,7.6,7.8,8.0'' and plot them again.