My symmetry-equivalent sites don't have the same energy

[Kane Shenton]

Hello!

Help, my symmetry-equivalent sites don't have the same energy!

__The problem:__

My system is insulating and antiferromagnetic Cr2O3 with a positive interstitial H in a hexagonal supercell. There are 6 equivalent* sites in a 'doughnut' around a large interstitial region. When I compute the total energy of these sites using CP2K, I get some weird artifacts that I've not been able to get rid of, namely that pairs of sites opposite each other in this ring have exactly the same energy (to about 1e-6 eV) but these pairs have different energy to one another (by about 1e-3 eV).

Normally this sort of energy discrepancy is not so important, but in my case it's crucial and I would like to know how to get rid of it. I think it's an egg-box effect of some kind, but I've tried a bunch of things to no avail.

* equivalent in that I generated them using symmetry operations from one initial relaxed structure, confirmed using pymatgen's structure matcher tool that they are indeed structurally the same and with VASP that they have the same total energy to 1e-6 eV.

__What I've tried so far:__

• - Crank up EPS_DEFAULT to 1e-18, CUTOFF to 3000 Ry and REL_CUTOFF to 500 Ry (with 4 grids)
• - Go down to just NGRID=1 at 1400 Ry
• - Switch between GAPW and GPW
• - Uniform translations of the structure (made ~ no difference to total energy)
• - Rotated unit cell by 120 degrees, keeping fractional coordinates the same (this resulted in a significant change in energy ~1 meV).
• - Flip all magnetic moments (no change in energy, as expected)
  
• - Switch to diagonalisation (rather than OT) (slight change in total energy but if anything adjacent equivalent sites got further from each other in energy)
  
• - Larger supercell (2x2x1 -> 3x3x2)
  
• - Changed from LDA+U Mulliken to LDA+U Lowdin
• - Changed from LDA+U to hybrid
• - Increased basis set on O from DZVP to TZVP (had to switch to Lowdin PLUS_U_METHOD due to convergence issues).
• - CP2K versions 7.1, 8.0, 8.1 and 8.2

None of these things came close to getting rid of the ~1 meV difference :(

__Minimal working example:__

Attached are 1) a cif file containing all six sites and 2) a tarball of six sets of cp2k inputs and outputs these sites (with a adjacent sites having an difference in energy of about 1 meV) that should all be energetically equivalent. These are for the relatively cheap 2x2x1 supercell and parameters that hopefully serve as a minimal working example.

Can someone explain why I'm getting different energies for these equivalent sites and how I might reduce this discrepancy.

Thanks and best wishes,

Kane

[Krack Matthias (PSI)]

Hello Kane

 

Without any constraint the symmetry of the Cr 3d orbital occupations will break and Cr atoms equivalent by symmetry will not show exactly the same occupation pattern on convergence. The numerical noise introduced by the employed grids in CP2K will always trigger such a break in symmetry, even for insanely high cutoffs and already small fluctuations in the 3d occupation patterns of topological equivalent Cr atoms can easily cause energy differences up to a milliHartree on convergence as you observe. CP2K has currently no tools implemented to impose the required symmetry strictly. I tried to keep at least the 3d occupations of the Cr atoms spherical by using an AMF setting which results in the following energies for the system configuration you provided (GPW, 2000 Ry cutoff, EPS_DEFAULT 1.0E-16):

 

# System    Energy [a.u.]

  00       -5320.66059565

  01       -5320.66059828

  02       -5320.66059579

  03       -5320.66059679

  04       -5320.66059787

  05       -5320.66059655

 

The total 3d occupations of symmetry equivalent Cr atoms, however, can still vary slightly, which causes the microHartree fluctuations in the total energy. I am afraid, it will be difficult to get it much better.

 

Best regards

 

Matthias

 

--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/addd5835-db77-45f0-9c48-74a5dea43aafn%40googlegroups.com.

[Kane Shenton]

Hi Matthias,

Thanks a lot for looking into this and for your helpful explanation!

Your energies look just about converged enough for my purposes actually so I'd like to understand how you got them.

I've had a look through the CP2K manual, github repo and forums but can't find mention of an AMF setting. Could you maybe post your input file / point me to a resource on this setting? 

Thanks and best wishes,

Kane

[Krack Matthias (PSI)]

Hi Kane

 

I enforced a smearing of the Cr 3d orbital occupations by adding the following section in the &DFT_PLUS_U of each Cr kind:

 

    &ENFORCE_OCCUPATION on

     EPS_SCF 1.0E-7

     ORBITALS -2 -1 +0 +1 +2

     SMEAR

    &END ENFORCE_OCCUPATION

 

Note, with that, you do not converge to the ground state. This would require the implementation of a symmetry constraint which keeps the same 3d occupation pattern for all Cr atoms equivalent by symmetry.

 

Best regards

 

Matthias

 

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/7044e83b-444c-4637-a0fb-b0cca61ebf98n%40googlegroups.com.

[Kane Shenton]

Hi Matthias,

Many thanks again for your reply explaining the how you got the configurations to have the same energy.

For this specific case (Cr3+ --> d^3, octahedrally coordinated), I think that smearing the occupations over all five d states is probably too 'non-physical'. I tried smearing just across the three t2g states and I got something in between your result and the non-smeared result. I thought about trying to enforce a specific matrix occupancy for symmetry-equivalent Cr ions 'by hand' since I know the mapping, but I think it's not possible with CP2K (?).

In any case, when I now look at energy differences between a given structure before and after a perturbation such as relaxation under an electric field I seem* to be able to resolve the energy differences (i.e. energy with electric field - energy without electric field) to ~1e-6 eV at least without any occupation smearing.

So for these kinds of grid effects it seems that, even if you have 1e-3 eV differences the total energy of in symmetry-equivalent structures, you can still resolve ~1e-6 eV differences for a given structure before and after a small perturbation. Would you say that makes sense or am I being too optimistic?

Cheers,

Kane

* When  I do this independently for each of my six configurations I get consistent results on that level.

[Krack Matthias (PSI)]

Hi Kane

 

You can assign each group of Cr ions equivalent by symmetry to a specific atomic kind and then enforce the desired d orbital occupations for that atomic kind and hope that the initial symmetry won’t break during the wavefunction optimization. Currently, CP2K provides no measure to keep an equal number of d electrons (smeared within the selected subset of d orbitals) per symmetry equivalent Cr ion within such a group during the SCF procedure. You will need to hack the code to achieve that additional symmetry constraint. You can contact me directly when you need help for that.

 

It is always good to be optimistic. I agree with you that it might be well possible to achieve the desired energy resolution for relative energies with and without perturbation in static calculations (i.e. without an ongoing cell size or shape change or significant atomic movements).

 

Cheers

 

Matthias

 

To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/e0ec4b94-2d43-4cb1-9ad5-32bde50f39f7n%40googlegroups.com.