CP2K Polaron DFT-MD with PBE+U

Author: Chris Ahart

---

The tutorial extends the previous tutorial 'CP2K Polaron Geometry Optimisation with PBE+U' to consider DFT-MD of polarons in anatase TiO2. Please work through part 1 before following this tutorial. Hopefully by the end of this tutorial you will be comfortable generating images such as the ones above.

CP2K is generally regarded as the best software for DFT-MD, as it was originally designed for MD simulations of large chemical and biological systems. The use of short-range basis sets (DZVP-MOLOPT-SR) in combination with efficient SCF minimization methods (Orbital Transform) and MD wavefunction propagators (Always Stable Predictor Corrector) result in highly scalable and stable DFT-MD calculations. Do note however that CP2K has generally very poor support for k-points, so you should use VASP if you would like to perform DFT-MD with k-points. Similarly, I recommend using VASP if you intend to study metallic or small band gap materials as Orbital Transform requires a finite band-gap (>0.5 eV).

All files corresponding to the tutorial are available on our group Github page. I have organised the files such that the examples you can run are in the folder /examples, and my origonal files are also included for reference in a seperate folder.

To begin, we perform PBE+U MD of a 2x2x1 supercell of anatase TiO2. Refer to your HPC provider of choice for a job submission script, and run CP2K using the provided input file. CP2K will generate a number of output files, including a log file and multiple MD specific output files: 'tio2-1.ener', 'tio2-frc-1.xyz', 'tio2-pos-1.xyz' and 'tio2-vel-1.xyz' which contain the energy, forces, positions and velocities of the atoms. You can load 'tio2-pos-1.xyz' in VMD and watch your atoms move.

CP2K does not have a good way to automatically increase the temperature during MD. While there is a simulated annealing keyword, the temperature profile is not linear and therefore I recommend performing equilibration using an external driving code such as MD or by doing so manually. To do this, we generally want to start at some low temperature and slowly increase the temperature. For example we can perform dynamics for some amount of time at 100 K, then restart at 200 K and run for another amount of time, and repeat this process until the desired temperature is reached.

For initial rapid equilibration you can use the following keywords

     ENSEMBLE NVE
     COMVEL_TOL 1e-10
     TEMP_TOL 10

This performs velocity rescaling whenever the temperature is 10 K above or below the target temperature, and will also re-scale velocities to ensure that the centre of mass remains constant.

Once you have performed a short initial equilibration at your target temperature, you should switch to the canonical sampling through velocity rescaling thermostat (CSVR)

     ENSEMBLE NVT
     &THERMOSTAT
       TYPE CSVR
       &CSVR
         TIMECON  1
       &END CSVR
     &END THERMOSTAT

These lines are added to the &MD subsection of the &MOTION section, following the structure that is visible in the current directory of the CP2K manual.

The keyword TIMECON is the time constant of the thermostat, the smaller the timecon the stronger the thermostatting. Ideally, you should use TIMECON 1 for equilibration and TIMECON 1000 for production runs however in practice I typically use TIMECON 1 for both. I generally perform DFT-MD with hybrid functionals, where equilibrating the system fully is not possible and therefore a small TIMECON is necessary. If you are performing PBE-MD, then you should be able to perform full equilibration.

Once you have performed equilibration with the CSVR thermostat for some appropriate amount of time, you may then choose to switch to ENSEMBLE NPT_F. Unfortunately CP2K will crash if you change the ensemble without making any further modifications to your input file, as you must explicitely tell CP2K to include a barostat to control the pressure

 &BAROSTAT 
 &END BAROSTAT

and to calculate the STRESS_TENSOR. In addition, you should add

 &CELL_REF 
  ABC A B C
 &END CELL_REF

This section 'CELL_REF' should define a cell with dimensions A B C slightly larger than the largest supercell size that will be sampled during DFT-MD. This reference cell is used to keep the FFT grid fixed during DFT-MD, and is also required when performing cell optimisations.

An alternative thermostat to CSVR is Nose-Hoover chains. Opinions differ on which of CSVR or Nose-Hoover chains is best, in my limited experience CSVR is faster for equilibration and therefore I use choose to use CSVR. I encourage you to try both and see which works best for your system.

Extensions:

You can apply what you learnt in the previous tutorial to set the charge and mulitplicity of the calculation to run PBE+U MD with an electron hole polaron. You can then plot images such as the ones at the beginning of this tutorial using the provided Python script as a guide.

Other remarks:

1. MD stability in CP2K is generally good. However, occasionally you will encounter bad geometries that may cause problems in the SCF convergence. You should check the column in 'tio2-1.ener' for 'UsedTime[s]' to make sure the value is roughly constant in time, and periodically check the CP2K output file to make sure the SCF is converging smoothly. There is a keyword IGNORE_CONVERGENCE_FAILURE with default value FALSE, that will cause CP2K to crash if the SCF convergence fails. You may set this to TRUE, so long as you are careful.