Converging magnetic systems in CP2K
Magnetic systems are common in nature, and present unique challenges for electronic structure calculations. In this tutorial I will explain how to calculate magnetic structures in CP2K, and share some techniques to accelerate SCF convergence. The techniques discussed are applicable to all magnetic systems, and may be helpful for other challenging convergence problems such as defects in crystals.
I will start with a ferromagnetic bulk Ni 2x2x2 supercell, and then consider the same structure instead as a slab in vacuum. These calculation will run in a few minutes on 20 cores on the HPC YOUNG, and therefore should be easy to follow. I would like to note that this is a particularly difficult system chosen deliberately as a learning example. For example, all magnetic systems that I studied during my PhD [1] including hematite, lepidocrocite, goethite and white rust [2] were significantly easier to converge. All files are available for reference on [3] A presentation given on 25/11/2022 which was the basis of this tutorial is available [4]. Note that some minor differences are present between the calculations shown in the presentation and the present tutorial.
Please note that I have now modified CP2K to allow for a non-integer number of electrons in the alpha and beta spin channels, significantly improving convergence.
This tutorial will be updated in the future, for now please refer to: [5] [6]
Bulk Ni
Non-magnetic structure
Let us start by considering the following input file:
&GLOBAL
PROJECT Ni_bulk
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
STRESS_TENSOR ANALYTICAL
&DFT
BASIS_SET_FILE_NAME BASIS_MOLOPT
POTENTIAL_FILE_NAME GTH_POTENTIALS
CHARGE 0
MULTIPLICITY 1
UKS .TRUE.
&QS
EPS_DEFAULT 1.0E-12
EXTRAPOLATION ASPC
&END QS
&MGRID
NGRIDS 5
CUTOFF 600.0
REL_CUTOFF 60.0
&END MGRID
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-8
MAX_SCF 500
ADDED_MOS 1000
MAX_DIIS 10
&DIAGONALIZATION
ALGORITHM STANDARD
&END DIAGONALIZATION
&MIXING
METHOD BROYDEN_MIXING
ALPHA 0.1
NBUFFER 10
&END MIXING
&SMEAR
ELECTRONIC_TEMPERATURE [K] 500.0
METHOD FERMI_DIRAC
&END SMEAR
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&PRINT
&E_DENSITY_CUBE ON
&END E_DENSITY_CUBE
&HIRSHFELD
SHAPE_FUNCTION DENSITY
&END HIRSHFELD
&END PRINT
&END DFT
&SUBSYS
&CELL
ABC 3.50579800 3.50579800 3.50579800
ALPHA_BETA_GAMMA 90.000 90.000 90.000
MULTIPLE_UNIT_CELL 2 2 2
&END CELL
&COORD
Ni 0.00000 0.00000 0.00000
Ni 0.00000 1.75290 1.75290
Ni 1.75290 0.00000 1.75290
Ni 1.75290 1.75290 0.00000
&END COORD
&TOPOLOGY
MULTIPLE_UNIT_CELL 2 2 2
&END TOPOLOGY
&KIND Ni
ELEMENT Ni
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q18
&END KIND
&END SUBSYS
&END FORCE_EVAL
We will be modelling bulk Ni with PBE and a DZVP basis set, with the use of Diagonalisation and Fermi-Dirac smearing as bulk Ni is a metal with no band gap. Importantly, we are performing unrestricted Kohn-Sham (UKS TRUE), with a multiplicity of 1 [7].
Part of the CP2K output file related to the initial guess is copied below. Recall that in CP2K the initial guess is obtained by superposition of atomic densities, and that a good starting guess is very important to ensure reliable convergence.
Atomic guess: The first density matrix is obtained in terms of atomic orbitals
and electronic configurations assigned to each atomic kind
Guess for atomic kind: Ni
Electronic structure
Total number of core electrons 10.00
Total number of valence electrons 18.00
Total number of electrons 28.00
Multiplicity not specified
S [ 2.00 2.00] 2.00 2.00
P [ 6.00] 6.00
D 8.00
Spin 1
Re-scaling the density matrix to get the right number of electrons for spin 1
# Electrons Trace(P) Scaling factor
288 288.033 1.000
Spin 2
Re-scaling the density matrix to get the right number of electrons for spin 2
# Electrons Trace(P) Scaling factor
288 288.033 1.000
We can see that CP2K is using an atomic guess composed of closed shell Ni atoms with 8 electrons in the d orbital, 4 electrons in each spin channel.
*** SCF run converged in 25 steps *** ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.): -5417.715937573439078
#Atom Element Kind Ref Charge Population Spin moment Net charge
1 Ni 1 18.000 9.000 9.000 -0.000 -0.000
2 Ni 1 18.000 9.000 9.000 -0.000 -0.000
Our SCF cycle converges in 25 steps, resulting a spin moment of 0.00 on each Ni atom. As such, we have started with a closed shell atomic guess and CP2K has converged to a closed shell density. This is the expected behaviour, as generally CP2K is not able to converge to a magnetic structure without specifying this manually as the starting guess.
Magnetic structure
To calculate magnetic Ni we must provide a suitable value for the keyword MULTIPLICITY. If we consider that the electronic configuration for Ni is [Ar] 3d⁸ 4s², we calculate a total multiplicity for 32 atoms of: 2S+1=2*2*1/2*32+1=65. Replacing MULTIPLICITY 1 with MULTIPLICITY 65 in our input file, we find the following in the CP2K output file:
Spin 1
Re-scaling the density matrix to get the right number of electrons for spin 1
# Electrons Trace(P) Scaling factor
320 288.033 1.111
Spin 2
Re-scaling the density matrix to get the right number of electrons for spin 2
# Electrons Trace(P) Scaling factor
256 288.033 0.889
As we have not changed the initial atomic guess the trace of the density matrix has equal number of electrons in spin channel 1 and 2, while the number of electrons calculated from the multiplicity is different. CP2K therefore imposes a scaling factor on the initial density matrix, rescaling it to match the specified multiplicity.
*** SCF run converged in 36 steps *** ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.): -5417.832272529466536
#Atom Element Kind Ref Charge Population Spin moment Net charge
1 Ni 1 18.000 9.344 8.656 0.687 -0.000
2 Ni 1 18.000 9.344 8.656 0.687 -0.000
We have therefore confirmed that the ferromagnetic solution for bulk Ni is -5417.715937573439078--5417.832272529466536 = 0.1 Ha = 3 eV lower in energy than the non-magnetic solution.
However, we also notice that the number of SCF steps has increased from 25 to 36. In order to accelerate the SCF convergence, we can consider instead a starting guess with multiplicity of 2S+1=2*1/2*32+1=33
Spin 1
Re-scaling the density matrix to get the right number of electrons for spin 1
# Electrons Trace(P) Scaling factor
304 288.033 1.055
Spin 2
Re-scaling the density matrix to get the right number of electrons for spin 2
# Electrons Trace(P) Scaling factor
272 288.033 0.944
*** SCF run converged in 19 steps *** ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.): -5417.832268755088080
#Atom Element Kind Ref Charge Population Spin moment Net charge
1 Ni 1 18.000 9.344 8.656 0.687 -0.000
2 Ni 1 18.000 9.344 8.656 0.687 -0.000
CP2K converges to the same final energy, however with 19 SCF steps instead of 36. This is clearly a better initial guess than the multiplicity of 65 from chemical intuition, and highlights how the results of a calculation can be used to accelerate future calculations. As all we are doing is changing the initial density matrix, we are not imposing any constraint on the system.
Improving the atomic guess
While changing the multiplicity leads to successful convergence of the magnetic structure, we have noted that CP2K is rescaling the density matrix. For this system with only Ni atoms this is not a problem, however if there were multiple atomic species then we would require a greater degree of control. This is performed by adding either a &BS section or the MAGNETIZATION keyword to the &KIND section of the atomic species that we are interested in.
&KIND Ni
ELEMENT Ni
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q18
&BS
&ALPHA
N 4 3
L 0 2
NEL 0 2
&END ALPHA
&BETA
N 4 3
L 0 2
NEL 0 -2
&END BETA
&END BS
&END KIND
The &BS contains &ALPHA and &BETA sections to control spin channel 1 (alpha) and spin channel 2 (beta). Keywords N, L and NEL refer to the principal quantum number N, the angular momentum L and 2x the number of electrons NEL. Therefore N=4, L=2 and NEL=2 will add 1 electron to the 3d orbital of the alpha spin channel. Recall that CP2K requires the number of electrons in the alpha spin channel to be greater than the number of electrons in the beta spin channel.
&KIND Ni
ELEMENT Ni
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q18
MAGNETIZATION 2.0
&END KIND
An alternative to the &BS section is to simply use the keyword MAGNETIZATION. The difference of 2 electrons between the alpha and beta spin channels corresponds to MAGNETIZATION 2. While the &BS section allows for a greater degree of control over which orbitals CP2K adds or removes electrons, the values for N, L and NEL are integers while the MAGNETIZATION keyword allows for floating point numbers. As such, fractional electrons can be specified in the initial atomic guess by using the MAGNETIZATION keyword.
Ni slab in vacuum
With an understanding of the multiplicity, &BS section and MAGNETIZATION keyword we are now ready to look at the more complicated system of a Ni slab in vacuum. The same input file is used, the only change being an increase in the z component of the cell from 7 Å to 20 Å which has the effect of adding 13 Å of vacuum. Converging this system is considerably more challenging the bulk Ni, and as such I have chosen to decrease the value of EPS_SCF from 1.0E-8 to 1.0E-5. The error introduced is <10 meV, reasonable for most calculations.
Non-magnetic structure
*** SCF run converged in 89 steps *** ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.): -5417.244330897219697
#Atom Element Kind Ref Charge Population Spin moment Net charge
1 Ni 1 18.000 9.002 9.002 0.000 -0.003
2 Ni 1 18.000 8.998 8.998 0.000 0.003
Starting with a non-magnetic structure, we can converge this in 89 SCF steps with a spin moment of 0.00 as found for bulk Ni.
Magnetic structure
#Atom Element Kind Ref Charge Population Spin moment Net charge
1 Ni 1 18.000 9.375 8.628 0.746 -0.003
2 Ni 1 18.000 9.339 8.658 0.680 0.003
For the magnetic structure, converging with a multiplicity of 33 or 65 is possible however requires many hundreds of SCF steps in additional to multiple restarts of the calculation. An example can be found in Gitlab for a multiplicity of 33 requiring a total of 584 SCF steps, restarting the calculation with EPS_SCF 1e-3 -> 1e-4 -> 5e-5 -> 1e-5. The results of such a brute force convergence reveal there are two different atomic environments, atoms on the two surface layers with a spin moment of 0.746 and atoms in the two bulk-like layers with a spin moment of 0.680 which is very similar to the spin moment of 0.687 observed for bulk Ni. As such, we can define two different atomic kinds in the input file with different multiplicities:
&KIND Ni_b
ELEMENT Ni
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q18
MAGNETIZATION 0.7
&END KIND
&KIND Ni_i
ELEMENT Ni
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q18
MAGNETIZATION 0.8
&END KIND
Recall that the MAGNETIZATION keyword is a float, while MULTIPLICITY keyword must be an integer. If the MAGNETIZATION does not sum to give the correct MULTIPLICITY, then the density matrix is rescaled. As such, we must carefully choose our MAGNETIZATION. I have chosen 0.7 and 0.8 as these are similar to the calculated spin moments of 0.680 and 0.746, giving integer multiplicity (16*0.7+16*0.8)+1=25.
ENERGY| Total FORCE_EVAL ( QS ) energy (a.u.): -5417.320042237174675
#Atom Element Kind Ref Charge Population Spin moment Net charge
1 Ni 1 18.000 9.375 8.629 0.746 -0.003
2 Ni 2 18.000 9.339 8.658 0.680 0.003
With the improved atomic guess we reach the same ground state electronic structure, but with 77 SCF steps instead of 616. The convergence is still challenging, requiring restarting the calculation with EPS_SCF 1e-3 -> 1e-4 -> 1e-5. I would expect that with further refinement of the MAGNETIZATION and MULTIPLICITY even faster convergence may be found.
Conclusion
I have demonstrated how a particular magnetic structure may be found by choosing an appropriate value of MULTIPLICITY, and that the &BS section and MAGNETIZATION keywords can be used to accelerate SCF convergence. In general I would suggest a workflow outlined below:
- Converge bulk crystal with multiplicity according to Hund’s rule
- Converge interface with multiplicity according to Hund’s rule
- Define MAGNETIZATION for each unique atomic environment
- Optimise MAGNETIZATION (float) consistently with MULTIPLICITY (integer)
Finally, I would like to stress again that this is a particularly difficult system chosen deliberately as a learning example. For example, all magnetic systems that I studied during my PhD including hematite, lepidocrocite, goethite and white rust were significantly easier to converge.
References
<References>
- ↑ PhD thesis of Christian Ahart [1]
- ↑ Polaronic structure of excess electrons and holes for a series of bulk iron oxides [2]
- ↑ Gitlab repository for this tutorial [3]
- ↑ Video of presentation which is the basis of this tutorial [4]
- ↑ Github repository containing my CP2K branch [5]
- ↑ Slides showing proof of concept for my CP2K code changes for Ni bulk and slab models [6]
- ↑ Wikipedia page on multiplicity [7]