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.
====== How to run Geometry Optimisation ======

===== Introduction =====
This tutorial is designed to illustrate how to relax the structure
of a system (without changing the cell dimensions) using ''CP2K''. We
use the relaxation of a water (H\(_2\)O) molecule as an example.

The example files are contained in {{:geometry_optimisation.tgz|geometry_optimisation.tgz}} that comes with this tutorial. The calculation was carried out with ''CP2K'' version
2.4.

It should be noted that before running the geometry optimisation,
the reader should have already know how to perform a simple
Kohn-Sham Density Functional Theory energy and force calculation
(this is covered in tutorial [[static_calculation|Calculating Energy and Forces using
QUICKSTEP]]), and they should also know how how to find a sufficient
grid cutoff for the static energy calculations (this is covered in
tutorial [[converging_cutoff|Converging the CUTOFF and REL_CUTOFF]]).

DIIS (direct inversion in the iterative subspace or direct inversion of the iterative subspace), also known as Pulay mixing, is an extrapolation technique. DIIS was developed by Peter Pulay in the field of computational quantum chemistry with the intent to accelerate and stabilize the convergence of the Hartree–Fock self-consistent field method.
===== Input Files =====
The input file for a geometry calculation is shown below:

<code cp2k>
&GLOBAL
  PROJECT H2O
  RUN_TYPE GEO_OPT
  PRINT_LEVEL LOW
&END GLOBAL
&FORCE_EVAL
  METHOD QS
  &SUBSYS
    &CELL
      ABC 12.4138 12.4138 12.4138
    &END CELL
    &COORD
      O      12.235322       1.376642      10.869880
      H      12.415139       2.233125      11.257611
      H      11.922476       1.573799       9.986994
    &END COORD
    &KIND H
      BASIS_SET DZVP-GTH-PADE
      POTENTIAL GTH-PADE-q1
    &END KIND
    &KIND O
      BASIS_SET DZVP-GTH-PADE
      POTENTIAL GTH-PADE-q6
    &END KIND
  &END SUBSYS
  &DFT
    BASIS_SET_FILE_NAME ./BASIS_SET
    POTENTIAL_FILE_NAME ./POTENTIAL
    &QS
      EPS_DEFAULT 1.0E-7
    &END QS
    &MGRID
      CUTOFF 200
      NGRIDS 4
      REL_CUTOFF 30
    &END MGRID
    &SCF
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-05
      MAX_SCF 200
      &DIAGONALIZATION T
        ALGORITHM STANDARD
      &END DIAGONALIZATION
      &MIXING T
        ALPHA 0.5
        METHOD PULAY_MIXING
        NPULAY 5
      &END MIXING
      &PRINT
        &RESTART OFF
        &END RESTART
      &END PRINT
    &END SCF
    &XC
      &XC_FUNCTIONAL PADE
      &END XC_FUNCTIONAL
    &END XC
  &END DFT
&END FORCE_EVAL
&MOTION
  &GEO_OPT
    TYPE MINIMIZATION
    MAX_DR    1.0E-03
    MAX_FORCE 1.0E-03
    RMS_DR    1.0E-03
    RMS_FORCE 1.0E-03
    MAX_ITER 200
    OPTIMIZER CG
    &CG
      MAX_STEEP_STEPS  0
      RESTART_LIMIT 9.0E-01
    &END CG
  &END GEO_OPT
  &CONSTRAINT
    &FIXED_ATOMS
      COMPONENTS_TO_FIX XYZ
      LIST 1
    &END FIXED_ATOMS
  &END CONSTRAINT
&END MOTION
</code>

The reader should already be familiar with the [[http://manual.cp2k.org/trunk/CP2K_INPUT/GLOBAL.html|''GLOBAL'']] and
[[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL.html|''FORCE_EVAL'']] sections. For geometry optimisation calculations, we
must set [[http://manual.cp2k.org/trunk/CP2K_INPUT/GLOBAL.html#desc_RUN_TYPE|''RUN_TYPE'']] in ''GLOBAL'' section to ''GEO_OPT'':

<code cp2k>
RUN_TYPE GEO_OPT
</code>

In this example, we note that we have chosen diagonalisation of the
Kohn-Sham Hamiltonian for the evaluation of wavefunctions, and used
Pulay mixing for the self-consistency loops. 5 histories are used
for Pulay mixing.

The important section for geometry optimisation settings are
contained in subsection [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/GEO_OPT.html|''GEO_OPT'']] of [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION.html|''MOTION'']] section. Note that
''GEO_OPT'' subsection only applies to the calculation where the cell
dimensions do not change. Calculations which allows the relaxation
of the cell are covered in a separate tutorial.

<code cp2k>
&GEO_OPT
  TYPE MINIMIZATION
  MAX_DR    1.0E-03
  MAX_FORCE 1.0E-03
  RMS_DR    1.0E-03
  RMS_FORCE 1.0E-03
  MAX_ITER 200
  OPTIMIZER CG
  &CG
    MAX_STEEP_STEPS  0
    RESTART_LIMIT 9.0E-01
  &END CG
&END GEO_OPT
</code>

The [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/GEO_OPT.html#desc_TYPE|''TYPE'']] keyword sets whether the geometry optimisation is for
finding the local minima (''MINIMIZATION'') or for finding the saddle
point transition state (''TRANSITION_STATE''). The keywords [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/GEO_OPT.html#desc_MAX_DR|''MAX_DR'']],
[[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/GEO_OPT.html#desc_MAX_FORCE|''MAX_FORCE'']], [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/GEO_OPT.html#desc_RMS_DR|''RMS_DR'']] and [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/GEO_OPT.html#desc_RMS_FORCE|''RMS_FORCE'']] set the criteria of whether an
optimised geometry is reached. ''MAX_DR'' and ''RMS_DR'' (in Bohr) are
the tolerance on the maximum and root-mean-square of atomic
displacements from the previous geometry optimisation iteration;
''MAX_FORCE'' and ''RMS_FORCE'' (in Bohr/Hartree) are the tolerance on
the maximum and root-mean-square of atomic forces. The geometry is
considered to be optimised //only when all four criteria are
satisfied//. The keyword [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/GEO_OPT.html#desc_MAX_ITER|''MAX_ITER'']] sets the maximum number of
geometry optimisation iterations. [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/GEO_OPT.html#desc_OPTIMIZER|''OPTIMIZER'']] sets the algorithm for
finding the stationary points; in this example we have chosen the
conjugate gradients (''CG'') method.

The [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/GEO_OPT/CG.html|''CG'']] subsection sets options for the conjugate gradients
algorithm. In this case, we have configured it so that no steepest
descent steps are to be performed before the start of the conjugate
gradients algorithm; and the CG algorithm should be reset (and one
steepest descent step is performed) if the cosine of the angles
between two consecutive searching directions is less than 0.9.

<code cp2k>
&CONSTRAINT
  &FIXED_ATOMS
    COMPONENTS_TO_FIX XYZ
    LIST 1
  &END FIXED_ATOMS
&END CONSTRAINT
</code>

We can add constraints to atomic movements by using the [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/CONSTRAINT.html|''CONSTRAINT'']]
subsection in ''MOTION'' section. In this example, we choose to fix
particular atoms using the [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/CONSTRAINT/FIXED_ATOMS.html|''FIXED_ATOMS'']] subsection. The keyword
[[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/CONSTRAINT/FIXED_ATOMS.html#desc_COMPONENTS_TO_FIX|''COMPONENTS_TO_FIX'']] sets which of the ''X'' ''Y'' ''Z'' directions are to
be fixed, and in this case, the atoms will be completely pinned in
all directions (''XYZ''). The list of atoms to be constrained are
given by the [[http://manual.cp2k.org/trunk/CP2K_INPUT/MOTION/CONSTRAINT/FIXED_ATOMS.html#desc_LIST|''LIST'']] keyword:

<code cp2k>
LIST 1 2 3 ... N
</code>

The numbers to the right of ''LIST'' are the list of atomic indices,
and correspond to the order (from top to bottom) of the atoms given
in the [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS/COORD.html|''COORD'']] subsection of [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/SUBSYS.html|''SUBSYS'']] (of [[http://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL.html|''FORCE_EVAL'']]). In our
example, we have fixed the oxygen atom during geometry optimisation,
so that the water molecule will not move around while its structure
is being relaxed.


===== Results =====
The example is run using the serial version of the ''CP2K'' binaries:

<code>
cp2k.sopt -o H2O.out H2O.inp &
</code>

After the job has finished, you should obtain the following files:

  * ''H2O.out''
  * ''H2O-pos-1.xyz''
  * ''H2O-1.restart''
  * ''H2O-1.restart.bak-1''
  * ''H2O-1.restart.bak-2''
  * ''H2O-1.restart.bak-3''

Again, the file ''H2O.out'' contains the main output of the
job. ''H2O-pos-1.xyz'' contains the trace of atomic coordinates at
each geometry optimisation step in the ''xyz'' file format. The last
set of atomic coordinates corresponds to the relaxed
structure. ''H2O-1.restart'' is a ''CP2K'' input file, similar to
''H2O.inp'', which contains the latest atomic coordinates of the water
molecule. Should the job die for some reason, you can continue the
job using the latest atomic coordinates by using command:

<code>
cp2k.sopt -o H2O.out H2O-1.restart &
</code>

You can of course also use ''H2O-1.restart'' as a template for writing
an input for further calculations using the relaxed atomic
structures.

The files ''H2O-1.restart.bak-*'' are backup restart files with
atomic coordinates obtained from the previous 1, 2 and 3 geometry
optimisation iterations. ''H2O-1.restart.bak-1'' should be the same as
''H2O-1.restart''.

In the main output file ''H2O.out'', at the end of each geometry
optimisation step, we will have the following information:

<code>
--------  Informations at step =     1 ------------
 Optimization Method        =                   SD
 Total Energy               =       -17.1643447508
 Real energy change         =        -0.0006776683
 Decrease in energy         =                  YES
 Used time                  =               90.837

 Convergence check :
 Max. step size             =         0.0336570168
 Conv. limit for step size  =         0.0010000000
 Convergence in step size   =                   NO
 RMS step size              =         0.0168136889
 Conv. limit for RMS step   =         0.0010000000
 Convergence in RMS step    =                   NO
 Max. gradient              =         0.0182785685
 Conv. limit for gradients  =         0.0010000000
 Conv. for gradients        =                   NO
 RMS gradient               =         0.0091312361
 Conv. limit for RMS grad.  =         0.0010000000
 Conv. for gradients        =                   NO
---------------------------------------------------
</code>

The above output segment states that at the end of geometry
optimisation step 1, the total energy of the system is
-17.1643447508 (Ha) and none of the criteria for optimised geometry
has been reached. The iteration therefore will carry on, until all
criteria becomes "''YES''".

At the end of geometry optimisation, one should obtain something
like:

<code>
--------  Informations at step =    11 ------------
 Optimization Method        =                   SD
 Total Energy               =       -17.1646204766
 Real energy change         =        -0.0000000529
 Decrease in energy         =                  YES
 Used time                  =               49.893

 Convergence check :
 Max. step size             =         0.0003393150
 Conv. limit for step size  =         0.0010000000
 Convergence in step size   =                  YES
 RMS step size              =         0.0001493298
 Conv. limit for RMS step   =         0.0010000000
 Convergence in RMS step    =                  YES
 Max. gradient              =         0.0001787448
 Conv. limit for gradients  =         0.0010000000
 Conv. in gradients         =                  YES
 RMS gradient               =         0.0000786642
 Conv. limit for RMS grad.  =         0.0010000000
 Conv. in RMS gradients     =                  YES
---------------------------------------------------
</code>

which clearly shows all criteria have been satisfied.

The final Kohn-Sham energies can be obtained at the end of the
output:

<code>
*******************************************************************************
***                    GEOMETRY OPTIMIZATION COMPLETED                      ***
*******************************************************************************

                   Reevaluating energy at the minimum

Number of electrons:                                                          8
Number of occupied orbitals:                                                  4
Number of molecular orbitals:                                                 4

Number of orbital functions:                                                 23
Number of independent orbital functions:                                     23

 Parameters for the always stable predictor-corrector (ASPC) method:

 ASPC order: 3

 B(1) =   3.000000
 B(2) =  -3.428571
 B(3) =   1.928571
 B(4) =  -0.571429
 B(5) =   0.071429

Extrapolation method: ASPC


SCF WAVEFUNCTION OPTIMIZATION

 Step     Update method      Time    Convergence         Total energy    Change
 ------------------------------------------------------------------------------
    1 Pulay/Diag. 0.50E+00    0.5     0.00005615       -17.1646204762 -1.72E+01
    2 Pulay/Diag. 0.50E+00    1.0     0.00000563       -17.1646347711 -1.43E-05

 *** SCF run converged in     2 steps ***


 Electronic density on regular grids:         -8.0000016293       -0.0000016293
 Core density on regular grids:                7.9999992554       -0.0000007446
 Total charge density on r-space grids:       -0.0000023739
 Total charge density g-space grids:          -0.0000023739

 Overlap energy of the core charge distribution:               0.00000004555422
 Self energy of the core charge distribution:                -43.83289054591484
 Core Hamiltonian energy:                                     12.82175605770555
 Hartree energy:                                              17.97395116120845
 Exchange-correlation energy:                                 -4.12745148966141

 Total energy:                                               -17.16463477110803

ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.):              -17.164634771108034
</code>
Open SourceMolecular Dynamics

User Tools

Site Tools

Page Tools

Open Source
Molecular Dynamics
