CP2K Open Source Molecular Dynamics

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Exercise 1: Electronic energy of the L-alanine crystal ======

Most parts of this exercise are adapted from the cubic Si example on the CP2K "How to" website:
[[https://www.cp2k.org/howto:static_calculation|How to Calculate Energy and Forces]]

As a first step, we perform a static self-consistent Kohn–Sham density functional theory (DFT) calculation to obtain the electronic energy. Our model system is the L-alanine crystal, consisting of 13 atoms per molecule and 4 molecules in the unit cell.

===== Input files =====

We first look at the files required for this calculation, that can be downloaded [[https://github.com/JWilhelm/CP2K_Computational_Methods_in_Crystallography/tree/main/1_Electronic_energy_of_the_L-alanine_crystal|here]]:

  * ''Lalanine.inp'': the main input file for the calculation, which defines the system and the job parameters
  * ''Lalanine.cif'': the coordinate file when defined externally

Parameter files:

  * ''BASIS_SET'': file containing parameters for the basis sets
  * ''GTH_POTENTIALS'': file containing parameters for the pseudopotentials

A list of basis set and pseudopotential files may be found in ''cp2k/data/'', which comes with a CP2K source release.

===== Input structure =====

Let us look at the main input: ''Lalanine.inp'' in detail. The filename and extension can be chosen freely. The input file is structured into ordered blocks and keywords, the order of which are unimportant. Each input block is referred to as a 'section' in this tutorial, and some sections are 'subsections' of other sections. Full definitions of the input file format and keywords are available via the CP2K input reference manual.

The main sections in the input file are:

  * ''GLOBAL'': contains general options for the CP2K run, such as the name of the job, the type of run etc.
  * ''FORCE_EVAL'': contains all parameters associated with the evaluation of forces on atoms, including the initial atomic coordinates.

We now look at each section in detail, starting with the ''GLOBAL'' section from ''Lalanine.inp'':
<code>
&GLOBAL
  PROJECT LALANINE
  RUN_TYPE ENERGY
  PRINT_LEVEL MEDIUM
&END GLOBAL
</code>
Here we perform a static energy calculation; in this case, we must set ''RUN_TYPE'' to ''ENERGY''. The keyword ''PROJECT'' sets the root name of the calculation, in this case ''Lalanine''. Any output files automatically generated by CP2K will have the name prefixed by ''Lalanine''. ''PRINT_LEVEL'' controls the default verbosity of the main output of CP2K, in this example it is set to 'medium'.


Next, we examine the ''FORCE_EVAL'' section line-by-line.

''METHOD Quickstep''

The keyword ''METHOD'' specifies ''QUICKSTEP'' as the approach for evaluating the forces on the atoms, i.e.
density functional theory using the Gaussian and Plane Waves (GPW) method.

The subsection ''SUBSYS'' defines the simulation unit cell and the initial coordinates of the atoms in the
calculation.
<code>
&SUBSYS
  &CELL
    ABC 5.94 12.274 5.806
    ALPHA_BETA_GAMMA 90.0 90.0 90.0
  &END CELL
  &COORD
    H 2.524500 3.041497 1.962428
    ...
  &END COORD
  &KIND H
    ELEMENT H
    BASIS_SET DZVP-GTH-PBE
    POTENTIAL GTH-PBE-q1
  &END KIND
  &KIND C
    ELEMENT C
    BASIS_SET DZVP-GTH-PBE
    POTENTIAL GTH-PBE-q4
  &END KIND
  &KIND N
    ELEMENT N
    BASIS_SET DZVP-GTH-PBE
    POTENTIAL GTH-PBE-q5
  &END KIND
  &KIND O
    ELEMENT O
    BASIS_SET DZVP-GTH-PBE
    POTENTIAL GTH-PBE-q6
  &END KIND
&END SUBSYS
</code>

The subsection ''CELL'' defines the simulation cell used in a calculation. In this example, we have a tetragonal cell with $a$=5.94 Å, $b$=12.274 Å, $c$=5.806 Å.

The initial atomic coordinates are specified in the ''COORD'' subsection. The default input format for atomic coordinates in CP2K is:
''<ATOM_KIND> X Y Z''
where $X$, $Y$, $Z$ are Cartesian coordinates in Å. This can be changed by configuring keyword ''SCALED'' to ''.TRUE.'' in the ''COORD''subsection, which makes the coordinate input $X$ $Y$ $Z$ fractional with respect to the lattice vectors. One can also change the unit for the Cartesian coordinates by setting the keyword ''UNIT'' within the subsection. ''<ATOM_KIND>'' should be a label that corresponds to the definition of the elements in the ''KIND'' subsections.

<note classic>
**Task:** Obtain the CIF file from the WebCSD (Refcode LALNIN23) and load it into CP2K in place of manual coordinates. You can compare your input to the example [[https://github.com/JWilhelm/CP2K_Computational_Methods_in_Crystallography/tree/main/1_Electronic_energy_of_the_L-alanine_crystal|here]].
</note>

The subsection ''KIND'' gives definitions of elements in the calculation. There must be one ''KIND'' subsection per element. In this example, we have defined the basis set to be ''DZVP-GTH-PBE'' (double-ζ with polarisation basis optimised for Goedecker–Teter–Hutter PBE pseudopotentials); and the pseudopotentials such as ''GTH-PBE-q4'' (Goedecker–Teter–Hutter PBE pseudopotential with 4 valence electrons).

The basis set and pseudopotential names must correspond to an existing entry in the corresponding basis set and pseudopotential files defined by the ''BASIS_SET_FILE_NAME'' and ''POTENTIAL_FILE_NAME'' keywords in the ''DFT'' subsection of the ''FORCE_EVAL'' section.
<note classic>
**Task:** Identify where these files are located.
</note>

<note important>
**Hint:** The path where CP2K is located is printed in the beginning of the output file.
</note>

After the ''SUBSYS'' section in ''Lalanine.inp'' comes the ''DFT'' subsection, which controls all aspects of the self-consistent Kohn–Sham density functional theory calculation. This subsection is relevant only if the ''METHOD'' keyword in ''FORCE_EVAL'' is set to ''QUICKSTEP''.

<code>
BASIS_SET_FILE_NAME  BASIS_SET
POTENTIAL_FILE_NAME  GTH_POTENTIALS
</code>

As already mentioned above, the keywords ''BASIS_SET_FILE_NAME'' and ''POTENTIAL_FILE_NAME'' set the files that contain basis set and pseudopotential parameters.
<code>
&QS
  EPS_DEFAULT 1.0E-10
&END QS
</code>

General control parameters for QUICKSTEP are specified in the ''QS'' subsection. The parameter ''EPS_DEFAULT'' defines the global default tolerance applied throughout the calculation. More specific tolerance parameters (''EPS_*''), if present, override this global value.

<code>
&MGRID
  NGRIDS 4
  CUTOFF 300
  REL_CUTOFF 60
&END MGRID
</code>

The ''MGRID'' subsection defines how the integration grid used in QUICKSTEP calculations should be set up. QUICKSTEP employs a multigrid scheme to represent Gaussian functions numerically: Narrow and sharp Gaussians are mapped onto finer grids than wider and smoother Gaussians. In this case, we are telling the code to set up 4 levels of multigrids, with the planewave cutoff of the finest grid set to 300 Ry, and with the grid spacing underneath any Gaussian functions to be finer than the equivalent planewave cutoff of 60 Ry. The users should read the tutorial [[https://manual.cp2k.org/trunk/methods/dft/cutoff.html|How to Converge the CUTOFF and REL_CUTOFF]] for details on how these parameters affect the grid constructed, and how to define a sufficient grid for their calculation.

<note classic>
**Task:** Increase the grid parameters and compare how this affects both the energy and the computational time.
</note>

<code>
&XC
  &XC_FUNCTIONAL PBE
  &END XC_FUNCTIONAL
&END XC
</code>

This subsection defines which exchange–correlation density functional we want to use. In this case we choose the PBE functional, which is consistent with the basis set and pseudopotential we have chosen.
<code>
&SCF
  SCF_GUESS ATOMIC
  EPS_SCF 1.0E-7
  MAX_SCF 300
  &DIAGONALIZATION
    ALGORITHM STANDARD
  &END DIAGONALIZATION
  &MIXING
    METHOD BROYDEN_MIXING
    ALPHA 0.4
    NBROYDEN 8
  &END MIXING
&END SCF
</code>

The ''SCF'' subsection defines all the settings related to methods used to find a self-consistent solution of the Kohn–Sham DFT formalism.

''SCF_GUESS'' specifies how the initial trial electron density function $\vec{r}$ is generated. In this example (''ATOMIC''), the initial density is obtained from overlapping atomic charge densities. A good starting point for the electron density in the self-consistency loop is important for achieving convergence efficiently. ''EPS_SCF'' sets the tolerance of the charge density residual. This overrides the ''EPS_DEFAULT'' value set in the ''QS'' subsection. ''MAX_SCF'' sets the maximum number of self-consistency loops QUICKSTEP is allowed to perform for each ground-state energy calculation.

<code>
&DIAGONALIZATION ON
  ALGORITHM STANDARD
&END DIAGONALIZATION
</code>

The ''DIAGONALIZATION'' subsection instructs CP2K to use the traditional diagonalisation method for obtaining the ground-state Kohn–Sham energy and electron density. The subsection heading also takes an argument; in this case it is set to ''ON'', which is equivalent to ''.TRUE.'' or ''T'', and enables diagonalisation. If the argument is omitted, the default is also ''.TRUE.''. The alternative to diagonalisation is the Orbital Transform (OT) method. In that case, the user should either delete the ''DIAGONALIZATION'' block or set the heading to ''OFF'' (or ''.FALSE.''), and instead add the ''OT'' subsection. The keyword ''ALGORITHM'' specifies the diagonalisation algorithm; here ''STANDARD'' means the standard LAPACK/SCALAPACK routines are used. 

<code>
&MIXING T
  METHOD BROYDEN_MIXING
  ALPHA 0.4
  NBROYDEN 8
&END MIXING
</code>

The ''MIXING'' subsection defines the parameters for charge mixing in a self-consistency calculation. It can take a value of either ''.TRUE.'' (''T'') or ''.FALSE.'' (''F''), which switches charge mixing on or off. The default is ''.TRUE.''. Note that this subsection only applies when using the diagonalisation method; the OT method uses a different approach to charge mixing, which is explained in other tutorials. The keyword ''ALPHA'' sets the mixing parameter: in this example, 40% of the output density is mixed with 60% of the input density to form the new input density in the next SCF iteration. The keyword ''METHOD'' specifies the mixing method; here Broyden mixing is used. The keyword ''NBROYDEN'' (an alias for ''NBUFFER'') sets the number of previous iterations to be stored in the Broyden mixing algorithm.

===== Running the calculation and analysing the output files =====

To run the calculation, the user needs to have a working CP2K package compiled in the system ''PATH''. The calculation can be started with the command:

<code>
cp2k.psmp -o Lalanine.out Lalanine.inp &
</code>

This launches CP2K in the background. The option ''-o'' redirects the CP2K output to the file ''Lalanine.out''. Note that the ''-o'' option appends output from successive runs to the same file. If you want to start from scratch, you should delete ''Lalanine.out'' before launching a new calculation.

After the job has finished, the following files should be present:
  - ''Lalanine.out''
  - ''Lalanine-RESTART.wfn''

The file ''Lalanine.out'' contains the main output of the job. The file ''Lalanine-RESTART.wfn'' contains the final Kohn–Sham wavefunctions from the calculation. Files of the form ''Lalanine-RESTART.wfn.bak-<n>'' store the wavefunctions obtained from the $n$-th previous SCF step; in this case, ''Lalanine-RESTART.wfn.bak-1'' contains the wavefunctions from the last SCF step. 

If you want to start a new calculation from previously calculated wavefunctions, change the keyword ''SCF_GUESS'' in the ''SCF'' subsection of the ''FORCE_EVAL'' section to:
<code>
SCF_GUESS RESTART
</code>

This requires that the new calculation uses the same ''PROJECT_NAME'' as the one that generated the restart files. Otherwise, the restart wavefunction files need to be renamed accordingly.

<note classic>
**Task:** Try this and check how it affects the computational time.
</note>

In the output file, CP2K reports quantities such as the number of electrons, the number of occupied orbitals, the SCF convergence behaviour, the total energy, and the atomic forces. During the class we will examine these sections of the output in detail.