CP2K 2024.1: HSE/ADMM UKS SCF cycling for neutral Be_Ga C2v-like state in periodic GaN

[Asif Billah]

Dear CP2K community,

I am trying to get results similar to those in a paper titled "Koopmans-tuned Heyd-Scuseria-Ernzerhof hybrid functional calculations of acceptors in GaN." They are doing VASP HSE hybrid functional calculations to obtain both deep localized polaronic states and shallow delocalized states of Be in GaN. Be acts as an AX center in GaN. Now, I want to do the same calculations using CP2K. I also want to use the HSE hybrid functional calculation. My immediate goal is to demonstrate the existence of shallow and deep Be states using CP2K HSE hybrid functional simulations.

Now, I am trying to converge a neutral Be substitutional acceptor in wurtzite GaN, targeting a metastable C2v-like configuration with an anisotropically delocalized basal plane hole.

My immediate goal is state discovery and structural relaxation. I am not yet using this small cell for quantitative formation energies or final transition levels.

System and method:
- CP2K version: 2024.1
- Periodic 3x3x2 wurtzite GaN supercell, 72 atoms
- One Ga substituted by Be
- Gamma point calculation
- Neutral open shell defect:
    CHARGE 0
    UKS TRUE
    MULTIPLICITY 2
  Total electrons: N_alpha = 320, N_beta = 319

- HSE-type screened hybrid setup:
    HF FRACTION 0.31
    POTENTIAL_TYPE SHORTRANGE
    OMEGA 0.10584 bohr^-1  (= 0.20 Angstrom^-1)

- GTH-PBE0 pseudopotentials with q13 Ga, q5 N, q4 Be
- DZVP-MOLOPT-PBE0-GTH basis sets with ADMM auxiliary basis sets
- Four-center HFX with ADMM; no RI-HFX approximation
- MGRID CUTOFF 600 Ry, REL_CUTOFF 60 Ry

The C3v-like axial polaron seed converges with the same general setup. It has one long axial Be-N bond and the unpaired spin localizes primarily on the axial N atom, as expected for the deep state. The simulation has converged when I changed ALPHA from 0.20 to 0.05. At first, it was cycling for neutral Be_Ga C3v-like state with ALPHA 0.20.

The problem appears specifically for my intended shallow C2v-like seed.

C2v-like initial geometry:
I perturb Be laterally within the basal plane so that two adjacent basal Be-N bonds are changed symmetrically. For the seed:
- selected basal Be-N bonds: 1.9033 and 1.9033 Angstrom
- third basal Be-N bond:     2.0433 Angstrom
- axial Be-N bond:           1.9573 Angstrom

The two selected basal N atoms are assigned a separate N kind only to provide an initial shared spin bias:

    &KIND N_PAIR
      ELEMENT N
      MAGNETIZATION 0.5
      BASIS_SET DZVP-MOLOPT-PBE0-GTH-q5
      BASIS_SET AUX_FIT admm-dz-q5
      POTENTIAL GTH-PBE0-q5
    &END KIND

There are no geometry, charge, spin, or orbital constraints after the initial atomic guess.

SCF issue:
I use diagonalization plus Broyden mixing:

    &SCF
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-6
      MAX_SCF 160
      &DIAGONALIZATION
        ALGORITHM STANDARD
      &END DIAGONALIZATION
      &MIXING
        METHOD BROYDEN_MIXING
        ALPHA 0.05
        BETA 1.50
        NBUFFER 8
      &END MIXING
    &END SCF

The first run reached residuals in the low 10^-4 range repeatedly, but did not converge after ~160 iterations. Its final residual was approximately 1.8E-4.

I then restarted from the saved WFN with:

    WFN_RESTART_FILE_NAME seed_from_previous-RESTART.wfn
    SCF_GUESS RESTART
    ALPHA 0.01

Here, I changed ALPHA from 0.05 to 0.01. The restart confirms:

    SCF PARAMETERS Density guess: RESTART

but it continues to oscillate. It occasionally reaches residuals around 1.6E-4 to 1.7E-4, then returns to 10^-3 or occasionally 10^-2. The continuation has not reached EPS_SCF = 1.0E-6.

The calculation is otherwise technically stable:
- no OOM, MPI, or node errors
- peak memory is about 5.7 GiB per MPI task
- 8 MPI ranks, 128 GB requested
- four-center HFX setup completes normally

Questions:
1. For this type of near degenerate open shell defect state, would you recommend a different SCF strategy than diagonalization plus Broyden mixing?

2. Would orbital transformation, a temporary small electronic smearing, level shifting, a different mixing method, or occupation/orbital following control be the most appropriate next step in CP2K?

3. Is there a recommended way in CP2K to stabilize a targeted metastable spin-polarized defect state without artificially constraining the final state?

4. Does the paired basal bond geometry and shared initial spin seed seem like a reasonable way to search for a C2v-like basin, or is there a more appropriate initialization strategy?

5. CP2K also prints the following type of warning:

   "*** Periodic Hartree Fock calculation
   *** requested with use of a truncated or shortrange potential. The cutoff
   *** radius is larger than half the minimal cell dimension..."

   The output reports:

       Interaction Potential: SHORTRANGE
       Omega: 0.10584
       Cutoff Radius [angstrom]: 14.8112

I understand this may be less problematic for an HSE-type screened potential than for an explicitly truncated Coulomb calculation, but I would appreciate guidance on whether this warning should affect the reliability or SCF stability of this small cell state discovery test.

Please address my approach and the issues I am facing. Any suggestions for improvement regarding my approach or other details would be appreciated. Thank you for your time and help!