In this exercise we are going to show how to obtain Phonon spectra using CP2K.

While some solid state code packages (like Quantum Espresso) contain functionality to calculate phonon spectra, CP2K does not. Using the Phonopy package it is nevertheless possible to calculate them using the Finite Displacement Method (FDM) with CP2K as the backend to calculate the forces.

You can install Phonopy in your home-directory on the server by running the following, which will download the package and its requirements:

pip install --user phonopy
# plus some tools required to handle CP2K input files:
pip install --user cp2k_tools

afterwards you will have the command phonopy available in ~/.local/bin. To make it known to the shell as any other command, run the following and logout and login back again:

echo 'export PATH="${HOME}/.local/bin:${PATH}"' >> ~/.bashrc

Now you should be able to run phonopy --help as follows and get some help text:

studentXY@tcopt2 ~ $ phonopy --help
Usage: phonopy [options]

Options:
  -h, --help            show this help message and exit
[...]

A typical Phonopy workflow looks as follows:

1. Phonopy takes a CP2K input file containing at least a cell specification plus coordinates for a crystal and generates a number of input file skeletons containing new cell definitions and new coordinates
2. You have to take these input file skeletons, complete them with the rest of the simulation parameters and run a energy+forces simulation. The input file skeleton will also contain an option for CP2K to write the forces into a separate file
3. Phonopy is then run a second time to parse these forces and generate a so-called force set out of it
4. The third time Phonopy is called, it is used to actually calculate and generate a plot of the phonon band structure for a given path

As a first example we are going to reproduce the text book band structure for bulk silicon.

Take the following CP2K input file:

Si.inp

&GLOBAL
   PROJECT Si
   RUN_TYPE ENERGY_FORCE
   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
         MAX_SCF 300
         ADDED_MOS 100
         &DIAGONALIZATION
            ALGORITHM STANDARD
            EPS_ADAPT 0.01
         &END DIAGONALIZATION
         &SMEAR  ON
            METHOD FERMI_DIRAC
            ELECTRONIC_TEMPERATURE [K] 300
         &END SMEAR
 
         &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
         ABC [angstrom] 5.4661639157319968 5.4661639157319968 5.4661639157319968
         PERIODIC XYZ
      &END CELL
      &COORD
         SCALED
         Si 0.875  0.875  0.875
         Si 0.875  0.375  0.375
         Si 0.375  0.875  0.375
         Si 0.375  0.375  0.875
         Si 0.125  0.125  0.125
         Si 0.125  0.625  0.625
         Si 0.625  0.125  0.625
         Si 0.625  0.625  0.125
      &END COORD
      &KIND Si
         ELEMENT Si
         BASIS_SET DZVP-MOLOPT-GTH
         POTENTIAL GTH-PBE
      &END KIND
 
   &END SUBSYS
&END FORCE_EVAL

Now let Phonopy parse the CP2K input file and have it generate a 2x2x2 supercell with one displacement, by running the following command:

phonopy --cp2k -c Si.inp -d --dim="2 2 2"

which should give you the following additional files:

$ ls
disp.yaml  phonopy_disp.yaml  Si.inp  supercell-001.inp  supercell.inp

supercell.inp is a skeleton with only the newly generated supercell, but without the displacement, while supercell-001.inp also contains the displacement of one Silicon atom and is therefore the calculation we need to run.

You therefore have to extend supercell-001.inp with the rest of the missing calculation parameters from the original input file Si.inp. Also make sure that you change project name to Si-supercell-001.

After running CP2K, you should get the following output file in addition to the usual CP2K output: Si-supercell-001-forces-1_0.xyz.

Now run Phonopy on this output file to generate the force set (afterwards found in a file called FORCE_SETS):

phonopy --cp2k -f Si-supercell-001-forces-1_0.xyz

… and run Phonopy again to get the band structure:

phonopy --cp2k -c Si.inp -p --dim="2 2 2" --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --band="1/2 1/2 1/2 0 0 0 1/2 0 1/2"

This command is supposed to open a window with the plot on your local machine. But this works only if you have setup the X11-Forwarding as described at the beginning of the lecture. Should this not work, then your best chance is to install the phonopy and cp2k-tools package on your local machine and run the last command on your local machine after having transferred all files in the current project folder to your local machine. Phonopy does not call CP2K directly so CP2K does not have to be installed.

Compare the generated plot to plots usually found in literature. Which special points have been specified in the above command? Change them to plot the path $\Gamma-X-K-\Gamma-L$ instead.