The actual energy – and therefore also the stress tensor – depends on many parameters, like the selected functional. This means that geometrical parameters like the lattice constant may also vary and therefore need to be optimized first when building a new geometry. While this could be done using CP2K's CELL_OPT run type, optimizing both the lattice/cell constants and the geometry simultaneously, we are going to do it manually here, especially since we can assume that only the lattice constant will actually change.

What we are using to determine the center volume (the volume for which the energy is minimal) is the Birch–Murnaghan equation of state (to be precise: the BM equation integrated over pressure), which links the energy and the volume using the minimal energy $E_0$, the center volume $V_0$, the bulk modulus $B_0$ and its derivative $B_1$:

\begin{align*} E(V) = E_0 + \frac{9 V_0 B_0}{16} \Bigg\{ \left[ \left(\frac{V_0}{V}\right)^{2/3} - 1 \right]^3 B_1 \; + \left[ \left(\frac{V_0}{V}\right)^{2/3} - 1 \right]^2 \left[ 6 - 4 \left(\frac{V_0}{V}\right)^{2/3} \right] \Bigg\} \end{align*}

Use the following input file as a starting point, and an adapted version of the script you documented in a previous exercise to generate a number of input files for different lattice constants and run the respective calculation. A good interval for the fraction of the lattice constant is $0.90-1.10$ with a step size of $0.025$.

Extract the energies and fit $E_0$, $V_0$, $B_0$, $B_1$ using the Birch–Murnaghan EOS and using the new $V0$ to determine the lattice constant. Plot your fit.

graphene.inp

&GLOBAL
  PROJECT graphene
  RUN_TYPE ENERGY
  PRINT_LEVEL MEDIUM
&END GLOBAL
 
&FORCE_EVAL
  METHOD Quickstep
  &DFT
    BASIS_SET_FILE_NAME  BASIS_MOLOPT
    POTENTIAL_FILE_NAME  POTENTIAL
 
    &POISSON
      PERIODIC XYZ
    &END POISSON
    &SCF
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-6
      MAX_SCF 300
 
      # The following settings help with convergence:
      ADDED_MOS 100
      CHOLESKY INVERSE
      &SMEAR ON
        METHOD FERMI_DIRAC
        ELECTRONIC_TEMPERATURE [K] 300
      &END SMEAR
      &DIAGONALIZATION
        ALGORITHM STANDARD
        EPS_ADAPT 0.01
      &END DIAGONALIZATION
      &MIXING
        METHOD BROYDEN_MIXING
        ALPHA 0.2
        BETA 1.5
        NBROYDEN 8
      &END MIXING
    &END SCF
    &XC
      &XC_FUNCTIONAL PBE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
 
  &SUBSYS
    &CELL
      # create a hexagonal unit cell:
      ABC 2.4612 2.4612 15.0
      ALPHA_BETA_GAMMA 90. 90. 60.
      SYMMETRY HEXAGONAL
      PERIODIC XYZ
    &END CELL
    &COORD
      SCALED
      C  1./3.  1./3.  0.
      C  2./3.  2./3.  0.
    &END
    &KIND C
      ELEMENT C
      BASIS_SET DZVP-MOLOPT-GTH
      POTENTIAL GTH-PBE
    &END KIND
  &END SUBSYS
 
&END FORCE_EVAL

bc -l <<< "5.6 * 12.3"

# you can also use variables and capture the output again in a variable:
x=1.025
a=$(bc -l <<< "$x * 2.4612")

a=3.54
sed -e "s|2.4612|$a|g" graphene.inp > "graphene_V-${x}.inp"