Symmetry Issue with HCP Lattice

[umbriquse]

The following code is attempting to obtain the self consistent field for an HCP lattice for Platinum. The error that is produced can be suppressed by setting the symmetry option for the Model to :off.

I don't understand why an HCP lattice throws this type of error, or why a lattice has to be a Unitary rotation matrix in order to use symmetry operations.

Code

using DFTK
using Unitful, UnitfulAtomic

latticeConstant = 10
α = [1, 1, 1.633]
HCPLattice = 0.5 * [
    [2 0 0];
    [1 √3 0];
    [0 0 2]
] .* α

element = ElementPsp(:Pt, psp = load_psp(list_psp(:Pt; functional = "lda")[2].identifier))
temp = austrip(300u"K")
atoms = [element => [[0., 0, 0], ones(3)/3]]

cutoffEnergy = 10.0
kGrid = ones(Int,3) * 4

lattice = latticeConstant * HCPLattice

model = model_LDA(lattice, atoms; temperature = temp)
basis = PlaneWaveBasis(model, cutoffEnergy; kgrid = kGrid)
scfGroundState = self_consistent_field(basis)
@info scfGroundState.energies

Error

ERROR: LoadError: spglib returned non-unitary rotation matrix
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] spglib_get_symmetry(::SArray{Tuple{3,3},Float64,2,9}, ::Array{Pair{ElementPsp,Array{Array{Float64,1},1}},1}; tol_symmetry::Float64) at /Users/jasonlehto/.julia/packages/DFTK/dtzTz/src/external/spglib.jl:82
 [3] symmetry_operations(::SArray{Tuple{3,3},Float64,2,9}, ::Array{Pair{ElementPsp,Array{Array{Float64,1},1}},1}; tol_symmetry::Float64, kcoords::Nothing) at /Users/jasonlehto/.julia/packages/DFTK/dtzTz/src/symmetry.jl:51
 [4] symmetry_operations at /Users/jasonlehto/.julia/packages/DFTK/dtzTz/src/symmetry.jl:49 [inlined]
 [5] Model(::Array{Float64,2}; n_electrons::Nothing, atoms::Array{Pair{ElementPsp,Array{Array{Float64,1},1}},1}, terms::Array{Any,1}, temperature::Float64, smearing::Nothing, spin_polarization::Symbol, symmetry::Symbol) at /Users/jasonlehto/.julia/packages/DFTK/dtzTz/src/Model.jl:129
 [6] model_atomic(::Array{Float64,2}, ::Array{Pair{ElementPsp,Array{Array{Float64,1},1}},1}; extra_terms::Array{Any,1}, kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:temperature,),Tuple{Float64}}}) at /Users/jasonlehto/.julia/packages/DFTK/dtzTz/src/Model.jl:154
 [7] model_DFT(::Array{Float64,2}, ::Array{Pair{ElementPsp,Array{Array{Float64,1},1}},1}, ::Symbol; extra_terms::Array{Any,1}, kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:temperature,),Tuple{Float64}}}) at /Users/jasonlehto/.julia/packages/DFTK/dtzTz/src/Model.jl:162
 [8] #model_LDA#12 at /Users/jasonlehto/.julia/packages/DFTK/dtzTz/src/Model.jl:170 [inlined]
 [9] top-level scope at /Users/jasonlehto/git/Dino.jl/test/issue.jl:21
 [10] include(::String) at ./client.jl:457
 [11] top-level scope at REPL[25]:1
in expression starting at /Users/jasonlehto/git/Dino.jl/test/issue.jl:21

[antoine-levitt]

Aaaaaaaaaaargh (sorry, it's a reflex when it comes to symmetries and spglib issues). Thanks for the report! I'm indeed not sure what's going on, looks like a tolerance set too low.

[antoine-levitt]

Hm no that's not it; either we messed up in the equations somewhere or it's an spglib bug.

[mfherbst]

Yikes ... symmetry again. Well thanks for the report!

So I took a brief look into this. So basically it happens directly after spglib has returned and while I could be wrong, I'm quite doubtful this is us screwing up, since it literally is the step where we check the returned rotation matrices to be unitary. It's not a simple transpose issue in the matrix (that I checked).

Also using spglib_standardize_cell on this lattice returns only garbage (two atoms on top of each other), which points to me that spglib is not able to digest this lattice properly. So regarding your problem, did you try to define the hcp lattice in a different way ... spglib seems to sometimes misbehave if you define lattices in an unconventional way. For example if you use ASE to construct an hcp Pt lattice (ase.build.bulk("Pt", crystalstructure="hcp", a=5)) works fine with DFTK.

Regarding the bug in DFTK, I think what we can do is to explicitly remove the rotation matrices, which we detect to not be unitary. What do you think @antoine-levitt.

[antoine-levitt]

Also using spglib_standardize_cell on this lattice returns only garbage (two atoms on top of each other), which points to me that spglib is not able to digest this lattice properly

Indeed, that's a good sign we're OK and spglib is not. Should we just open a bug in spglib? I'd rather avoid monkey-patching DFTK further to avoid upstream issues.

[mfherbst]

We can. I expect it goes into a similar direction than spglib issue 101, however.

[antoine-levitt]

Could it also be that the 1.633 is badly rounded, resulting in an almost-symmetric lattice?

[mfherbst]

I went down to a symmetry tolerance 1e-2 and that did not help. It's more substantial than that form looking at the non-orthogonality numbers.

[antoine-levitt]

OK, I'll open an issue at spglib and see what they say. At the very least it's a simpler reproducer than 101.

In the meantime @umbriquse you can just disable symmetry or, if that is too slow, discard the non-unitary rotations.

[antoine-levitt]

Seems to be a transpose issue again?

import numpy as np
import spglib

lattice = np.array([[2.0, 0.0, 0.0],
                    [1.0, 1.7320508075688772, 0.0],
                    [0.0, 0.0, 2.0]])
positions = [[0., 0., 0.]             ]
numbers = [1, ] * 1
cell = (lattice, positions, numbers)

syms = spglib.get_symmetry(cell)
# lattice = np.transpose(lattice)
for r in syms['rotations']:
    rcart = np.dot(np.dot(lattice, r), np.linalg.inv(lattice))
    print(np.dot(np.transpose(rcart), rcart))

works if I put the transpose in. Did we mess up?

[mfherbst]

Are you sure this is the same setup? You only have one atom, but the failing example has two.

[antoine-levitt]

Doesn't change anything

[mfherbst]

Did you take care of all the transposes in python versus C versus julia? You remembered that spglib changes the ordering between the C and python interfaces?

[antoine-levitt]

No, I'm always very confused by this, I was hoping you'd remember :-)

[mfherbst]

That's why I would go for an example in C ... one less transpose to worry about ...

[antoine-levitt]

wouldn't be too sure about that, C has row major ordering?

[mfherbst]

yes, ah ... I get your point. Ok so the lattice definition in the spglib python interface takes lattice vectors as rows that is important to keep in mind. But of course compared to Julia python is row-major. So when you actually do matmuls you need to transpose the lattice to get the same thing we do.

[antoine-levitt]

But then it works?

[mfherbst]

Let's be clear about this:

DFTK

• Column-major
• Lattice vectors as columns

Spglib

• Row-major
• Lattice vectors as rows

So you do not need to take a transpose when using our lattice vectors in spglib and to get the same algebra we do for comparing things you also need to take a transpose when using the lattice matrix in the comparison to get the same algebra we do.

[mfherbst]

(Reason I don't believe it's just a transpose missing is that we have used spglib on non-symmetric lattice matrices)

[antoine-levitt]

I think I'm more confused than before. The correct snippet is

import numpy as np
import spglib

lattice = np.array([[2.0, 0.0, 0.0],
                    [1.0, 1.7320508075688772, 0.0],
                    [0.0, 0.0, 2.0]])
print(lattice)
positions = [[.0,.0,.0]             ]
numbers = [1, ] * 1
cell = (np.transpose(lattice), positions, numbers)

syms = spglib.get_symmetry(cell)
for r in syms['rotations']:
    rcart = np.dot(np.dot(lattice, r), np.linalg.inv(lattice))
    print(np.dot(np.transpose(rcart), rcart))

where the lattice gets transposed in the call to spglib_get_symmetry to accomodate spglib's python requirements that the vectors a b c are given as rows. Then there's no problem. So the problem is in DFTK, no?

[mfherbst]

Actually, I think you are right. I looked at the code again and it seems to me the transpose might be missing here:

https://github.com/JuliaMolSim/DFTK.jl/blob/master/src/external/spglib.jl#L56

in the lattice matrix. From the spglib docs it seems that spglib does the lattice vectors in columns, same as we do, but of course you need to transpose once for column major -> row major.

If that is indeed the case, it needs fixing here, too: https://github.com/JuliaMolSim/DFTK.jl/blob/master/src/external/spglib.jl#L119

I wonder why we did not spot this before ... we definitely used lattices which are not symmetric.

[antoine-levitt]

https://spglib.github.io/spglib/variable.html#id2 what's the memory layout of this in C?

[mfherbst]

[a_x, b_x, c_x, a_y, b_y, c_y, a_z, b_z, c_z ] no?

[antoine-levitt]

Then yes it needs to be transposed

[mfherbst]

@antoine-levitt Are you fixing it or should I?

[mfherbst]

• Implement bugfix
• Once we found the bug, we should add this example as a testcase. At least as an abinit testcase.

[antoine-levitt]

Not right now, feel free if you have time

[mfherbst]

Me neither, but I put it on the pile.

[antoine-levitt]

Fun fact: Julia's transpose is a lazy wrapper, so if you pass it to a C function it'll get ignored

[antoine-levitt]

There were actually two bugs here: the lattice was passed untransposed, and the rotations were passed untransposed in get_stabilized_reciprocal_mesh. Maybe they compensated each other...

[antoine-levitt]

Should be fixed by #314. Thanks a lot @umbriquse ! It's wonderful having users to actually find bugs for you :-)

[mfherbst]

I can only agree 😄

[antoine-levitt]

This was a pretty serious bug, tag a new release when the current round of PRs are in?

[mfherbst]

of course!