# Open SourceMolecular Dynamics

### Sidebar

#### For Developers

events:2018_summer_school:converging_cutoff

This is an old revision of the document!

# Converging the cutoff for a more difficult problem

## Input files

The complete set of files for this exercise can be found here.

This exercise is similar to the previous one, but uses a setup and system more typical of CP2K usage. We will use a system of 32 H2O water molecules within a periodic box. Here is the input template:

&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT cuttoff-test
RUN_TYPE ENERGY_FORCE
&END GLOBAL

&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME GTH_POTENTIALS
WFN_RESTART_FILE_NAME ../cuttoff-test-RESTART.wfn
CHARGE 0
MULTIPLICITY 1
&MGRID
NGRIDS 4
CUTOFF LT_cutoff
REL_CUTOFF LT_rel_cutoff
&END
&QS
EPS_DEFAULT 1.0E-12
METHOD GPW
&END

&SCF
SCF_GUESS RESTART
EPS_SCF 5.e-7
MAX_SCF 15
&OT
PRECONDITIONER FULL_ALL
MINIMIZER DIIS
&END OT
&OUTER_SCF
EPS_SCF 5.0E-7
MAX_SCF 1
&END OUTER_SCF
&END SCF

&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&XC_GRID
! defaults
XC_SMOOTH_RHO NONE
XC_DERIV PW
&END XC_GRID
&END XC

&END DFT
&SUBSYS
&CELL
ABC 9.8528 9.8528 9.8528
PERIODIC XYZ
&END CELL

&KIND H
BASIS_SET DZVP-MOLOPT-SR-GTH-q1
POTENTIAL GTH-PBE-q1
&END

&KIND O
BASIS_SET DZVP-MOLOPT-SR-GTH-q6
POTENTIAL GTH-PBE-q6
&END KIND

&TOPOLOGY
COORDINATE XYZ
COORD_FILE_NAME ../structure.xyz
CONNECTIVITY OFF
&END TOPOLOGY
&END SUBSYS

&PRINT
&FORCES
&END
&END
&END FORCE_EVAL

Compared to the Si example, this is a larger system, we are using the OT optimizer in a good setup for a small to medium insulating system:

    &SCF
SCF_GUESS RESTART
EPS_SCF 5.e-7
MAX_SCF 15
&OT
PRECONDITIONER FULL_ALL
MINIMIZER DIIS
&END OT
&OUTER_SCF
EPS_SCF 5.0E-7
MAX_SCF 1
&END OUTER_SCF
&END SCF

and we are also saving the forces on the atoms

  &PRINT
&FORCES
&END
&END

We save the forces as for many purposes (MD) converging the forces reasonably is more important than the total energy of the system.

## Running the system

The runcutoff file is a shell script as before to generate the different input files:

#!/bin/bash

cutoffs="100 200 300 400 500 600 700 800 900 1000 1100 1200"

template_file=input_template.inp
input_file=input.inp

rel_cutoff=60

for ii in $cutoffs ; do work_dir=cutoff_${ii}Ry
if [ ! -d $work_dir ] ; then mkdir$work_dir
else
rm -r $work_dir/* fi sed -e "s/LT_rel_cutoff/${rel_cutoff}/g" \
-e "s/LT_cutoff/${ii}/g" \$template_file > $work_dir/$input_file
done
remember to make it executable

When you run the shell script you should get a series of directories, cutoff_\${cutoff}Ry. Run the input files in each directory (you may want to setup a script to do this).

At the end you should have a set of output files that contain the total energy of the system and the forces on each atom.

1. Extract and plot the total energy of the system as a function of cutoff
2. Extract and plot the total force on the system as a function of cutoff (search for 'SUM OF ATOMIC FORCES')
3. Extract and plot the force on some chosen atoms from the system as a function of cutoff

## Analysis

What is converged?

Compare the convergence of forces to the default convergence criteria for geometry optimization.

What sets the required cutoff? It is the basis set (which is dictated by the pseudopotentials). You will need to be able to represent the Gaussian with largest exponent well on the realspace grids. Oxygen, being very electronegative (on the right of the period table with many protons) has very contracted 2s states. You can see in the output

       Normalised Cartesian orbitals:

Set   Shell   Orbital            Exponent    Coefficient

1       1    2s               10.389228       0.396646
3.849621       0.208811
1.388401      -0.301641
0.496955      -0.274061
0.162492      -0.033677

That there is a Gaussian with an exponent of 10.4 Bohr-2. If we compare to the basis set for Silicon

 Si DZVP-MOLOPT-GTH DZVP-MOLOPT-GTH-q4
1
2 0 2 6 2 2 1
2.693604434572  0.015333179500 -0.073325401200 -0.005800105400  0.023996406700  0.043919650100
1.359613855428 -0.283798205000  0.484815594600 -0.059172026000  0.055459199900  0.134639409600
0.513245176029 -0.228939692700 -0.276015880000  0.121487149900 -0.269559268100  0.517732111300
0.326563011394  0.728834000900 -0.228394679700  0.423382421100 -0.259506329000  0.282311245100
0.139986977410  0.446205299300 -0.018311553000  0.474592116300  0.310318217600  0.281350794600
0.068212286977  0.122025292800  0.365245476200  0.250129397700  0.647414251100  0.139066843800

we see that the largest exponent is only 2.7 Bohr-2, so can be represented on a much coarser grid.

Task If you like, have a look at the BASIS_MOLOPT file (in the data directory, or online here) and see how the exponents change across the periodic table

The convergence is largely dominated by the calculation of the gradient terms in a GGA functional (compare a simulation with LDA to the PBE used here). The evaluation of these terms on the grids are demanding, and very dependent on the functional.

    &XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&XC_GRID
! defaults
XC_SMOOTH_RHO NONE
XC_DERIV PW
&END XC_GRID
&END XC

For BLYP functional some smoothing needs to be applied. The smoothing may also converge forces more rapidly than the default settings, but at the expense of modifying the functional slightly.

compare to the previous calculation, but using a smoothing section in the XC section.

    &XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&XC_GRID
XC_SMOOTH_RHO NN50
XC_DERIV NN50_SMOOTH
&END
&END XC

compare the convergence of LDA and BLYP to PBE.

&XC_FUNCTIONAL PADE # or BLYP
&END XC_FUNCTIONAL
Also change the psuedo potential to the appropriate functional.
    &KIND O
BASIS_SET DZVP-MOLOPT-SR-GTH-q6
&END KIND