In this exercise we will compare linear scaling SCF 10.1021/ct200897x with normal SCF methods. We will run simulations on a 2D polymer system containing up to 2400 atoms. The 2D polymer system is described in more detail in this recent publication 10.1002/adma.201304705. To make this exercise computationally feasible, will us the efficient Density Functional based Tight Binding (DFTB) method 10.1103/PhysRevB.58.7260. It requires only a minima basis, but delivers nevertheless reasonable results due to an empirical correction term called repulsion potential.
In CP2K, Linear scaling (LS) methods are relatively new. Whereas LS is generally efficient with methods that require a minimal basis, the cross-over between LS and OT might require 10000s of atoms for DFT, dense 3D systems with DZVP-MOLOPT basis sets.
Calculate the energy of the 2D-polymer for various different system sizes, by using NREPX=1, 2, 3, …. Perform each calculation with the linear scaling and the normal SCF method. Plot the time required for the last SCF cycle of each calculation vs. the system size. Extrapolate the time required for an SCF cycle of a system containing 10×10 unit cells.
The number of seconds required of each SCF cycle is printed as part of the regular output.
For the normal SCF it's the 4th column in the following table:
Step Update method Time Convergence Total energy Change
------------------------------------------------------------------------------
1 P_Mix/Diag. 0.40E+00 3.7 0.15459443 -957.7477897329 -9.58E+02
2 P_Mix/Diag. 0.40E+00 5.6 0.06222675 -957.8604155443 -1.13E-01
3 DIIS/Diag. 0.11E-01 5.6 0.02839814 -957.9170257078 -5.66E-02
4 DIIS/Diag. 0.11E-01 5.6 0.02966521 -957.9973636669 -8.03E-02
5 DIIS/Diag. 0.22E-01 5.6 0.01746285 -957.9863069368 1.11E-02
For the linear scaling SCF it's the last column in lines like these:
SCF 1 -957.747786827 -957.747786827 12.559000 ... SCF 2 -957.970635063 -0.222848235 11.204000 ... SCF 3 -957.994490590 -0.023855528 11.233000
Calculate the energy of the single unit-cell with the linear scaling method at different filter thresholds, by changing EPS_FILTER. Plot the time required for the calculation and the obtained total energy vs the threshold. How does this compare with the total energy obtained from the normal SCF method? Note that this requires well converged SCF calculations (EPS_SCF). Also have a look at the sparsity report, which matrix is more dense (S or P) ?
Finally, compare the efficiency of three different algorithms (using NREP 2×2). The following three options can be tried:
sign matrix iterations 10.1021/ct200897x
! linear scaling SCF
&LS_SCF
! SIGN_MATRIX, only effective if the chemical potential is known in advance
PURIFICATION_METHOD SIGN_MATRIX
! threshold used to determine sparsity and thus speed and accuracy
EPS_FILTER 1E-7
! convergence for the SCF
EPS_SCF 1E-5
! chemical potential, an energy from within in the HOMO-LUMO gap
MU -0.15
S_PRECONDITIONER ATOMIC
&END
trs4 iterations 10.1063/1.1559913
! linear scaling SCF
&LS_SCF
! TRS4, does not need an estimate for the chemical potential
PURIFICATION_METHOD TRS4
! threshold used to determine sparsity and thus speed and accuracy
EPS_FILTER 1E-7
! convergence for the SCF
EPS_SCF 1E-5
! chemical potential, an energy from within in the HOMO-LUMO gap
MU -0.15
S_PRECONDITIONER ATOMIC
&END
curvy steps 10.1063/1.1558476
! linear scaling SCF
&LS_SCF
! TRS4, does not need an estimate for the chemical potential
PURIFICATION_METHOD TRS4
! threshold used to determine sparsity and thus speed and accuracy
EPS_FILTER 1E-7
! convergence for the SCF
EPS_SCF 1E-5
! chemical potential, an energy from within in the HOMO-LUMO gap
MU -0.15
S_PRECONDITIONER NONE
&CURVY_STEPS
&END CURVY_STEPS
&END
In making your judgement, consider the fact that the rate of convergence might not be the same for all algorithms.
From a practical point of view, the curvy_steps algorithm seems very robust in converging the SCF of electronically difficult systems, while the sign matrix algorithm is only worthwhile considering if the chemical potential is known (or fixed).
! Number of unit cell replications, increase this to create a large polymer.
@SET NREPX 1
@SET NREPY 1
&GLOBAL
PROJECT 2d_polymer
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
&DFT
! linear scaling SCF
&LS_SCF
! TRS4, does not need an estimate for the chemical potential
PURIFICATION_METHOD TRS4
! threshold used to determine sparsity and thus speed and accuracy
EPS_FILTER 1E-5
! convergence for the SCF
EPS_SCF 1E-5
! chemical potential, an energy from within in the HOMO-LUMO gap
MU -0.15
S_PRECONDITIONER ATOMIC
&END
! traditional SCF, diagonalization by default
&SCF
&END
&QS
LS_SCF ! Turns on linear scaling DFT. Comment this line to run with standard SCF.
! =============================================================================================
! === Description of Hamiltonian (Density Functional Tight Binding + Smooth Particle Meshed Ewald) ===
METHOD DFTB
&DFTB
HB_SR_GAMMA
SELF_CONSISTENT T
DO_EWALD T
DISPERSION T
&PARAMETER
! these parameters are in cp2k/data
PARAM_FILE_PATH ./DFTB/scc
PARAM_FILE_NAME scc_parameter
UFF_FORCE_FIELD uff_table
&END PARAMETER
&END DFTB
&END QS
&POISSON
&EWALD
EWALD_TYPE SPME
! scale the number of points with the system size, roughly 1pt/A.
GMAX 40*${NREPX} 40*${NREPY} 40
O_SPLINE 5
&END EWALD
&END POISSON
! ==============================End of Hamiltonian Description ================================
! =============================================================================================
&END DFT
&SUBSYS
&CELL
! a 2D system, z-direction is non-periodic
ABC 23.577999 40.8383 30
MULTIPLE_UNIT_CELL ${NREPX} ${NREPY} 1
&END CELL
&TOPOLOGY
MULTIPLE_UNIT_CELL ${NREPX} ${NREPY} 1
&END
&COORD
@INCLUDE '2d_polymer.coord'
&END COORD
&END SUBSYS
&END FORCE_EVAL
Download the following files into your project directory:
You can unpack them with the following commands:
$ gunzip 2d_polymer.coord.gz