View source for 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 <ref>PhD thesis of Christian Ahart [https://discovery.ucl.ac.uk/id/eprint/10153463/2/thesis.pdf]</ref> including hematite, lepidocrocite, goethite and white rust <ref>Polaronic structure of excess electrons and holes for a series of bulk iron oxides [https://doi.org/10.1039/C9CP06482F]
</ref> were significantly easier to converge.  All files are available for reference on <ref>Gitlab repository for this tutorial [https://gitlab.doc.ic.ac.uk/rgc/converging-magnetic-systems-in-cp2k]</ref> A presentation given on 25/11/2022 which was the basis of this tutorial is available <ref>Video of presentation which is the basis of this tutorial [https://imperiallondon-my.sharepoint.com/:v:/g/personal/cahart_ic_ac_uk/EY2b8Mkqgx5Kgb0Ij0hOUm8B3RerDQ52ptRNc9_dlEDWeQ]</ref>. 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:'''
<ref>Github repository containing my CP2K branch [https://github.com/chrisahart/cp2k/tree/nelectron-specify]</ref>
<ref>Slides showing proof of concept for my CP2K code changes for Ni bulk and slab models [https://gitlab.doc.ic.ac.uk/rgc/converging-magnetic-systems-in-cp2k/-/blob/master/slides/magnetic_cp2k_changes.pdf]</ref>

[[File:NiSlab.png |thumb|200px|32 atom Ni slab in vacuum.]]

=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 <ref>Wikipedia page on multiplicity [https://en.wikipedia.org/wiki/Multiplicity_(chemistry)]</ref>. 

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 [https://gitlab.doc.ic.ac.uk/rgc/converging-magnetic-systems-in-cp2k 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 [https://discovery.ucl.ac.uk/id/eprint/10153463/2/thesis.pdf PhD] including [https://doi.org/10.1039/C9CP06482F hematite, lepidocrocite, goethite and white rust] were significantly easier to converge.


=References=
<References>


__FORCETOC__