Refactor symmetries & symmetrization (JuliaMolSim#244)

docs/src/advanced/symmetries.md

@@ -70,7 +70,7 @@ and a `[4, 4, 4]` Monkhorst-Pack grid.
 First we perform the DFT calculation disabling symmetry handling
 ```@example symmetries
 model = model_LDA(lattice, atoms)
-basis_nosym = PlaneWaveBasis(model, Ecut; kgrid=kgrid, enable_bzmesh_symmetry=false)
+basis_nosym = PlaneWaveBasis(model, Ecut; kgrid=kgrid, use_symmetry=false)
 scfres_nosym = @time self_consistent_field(basis_nosym, tol=1e-8)
 nothing  # hide
 ```

src/Model.jl

@@ -1,10 +1,8 @@
 # Contains the physical specification of the model
 
-"""
-A physical specification of a model. 
-Contains the geometry information, but no discretization parameters.
-The exact model used is defined by the list of terms.
-"""
+# A physical specification of a model. 
+# Contains the geometry information, but no discretization parameters.
+# The exact model used is defined by the list of terms.
 struct Model{T <: Real}
     # Lattice and reciprocal lattice vectors in columns
     lattice::Mat3{T}
@@ -14,7 +12,7 @@ struct Model{T <: Real}
     dim::Int # Dimension of the system; 3 unless `lattice` has zero columns
 
     # Electrons, occupation and smearing function
-    n_electrons::Int # not necessarily consistent with `atoms` field
+    n_electrons::Int # usually consistent with `atoms` field, but doesn't have to
 
     # spin_polarization values:
     #     :none       No spin polarization, αα and ββ density identical,
@@ -24,26 +22,59 @@ struct Model{T <: Real}
     #     :full       Generic magnetization, non-uniform direction.
     #                 αβ, βα, αα, ββ all nonzero, different
     #     :spinless   No spin at all ("spinless fermions", "mathematicians' electrons").
-    #                 Difference with :none is that the occupations are 1 instead of 2
-    spin_polarization::Symbol  # :none, :collinear, :full, :spinless
+    #                 The difference with :none is that the occupations are 1 instead of 2
+    spin_polarization::Symbol
 
     # If temperature=0, no fractional occupations are used.
     # If temperature is nonzero, the occupations are
     # `fn = max_occ*smearing((εn-εF) / temperature)`
     temperature::T
     smearing::Smearing.SmearingFunction # see Smearing.jl for choices
 
-    atoms::Vector{Pair} # Vector of pairs Element => vector of vec3 (positions, fractional coordinates)
-    # Possibly empty. Right now, the consistency of `atoms` with the different terms is *not* checked
+    # Vector of pairs Element => vector of vec3 (positions, fractional coordinates)
+    # Possibly empty.
+    # `atoms` just contain the information on the atoms, it's up to the `terms` to make use of it (or not)
+    atoms::Vector{Pair}
 
     # each element t must implement t(basis), which instantiates a
     # term in a given basis and gives back a term (<: Term)
     # see terms.jl for some default terms
     term_types::Vector
+
+    # list of symmetries of the model
+    symops::Vector{SymOp}
 end
 
-function Model(lattice::AbstractMatrix{T}; n_electrons=nothing, atoms=[], terms=[], temperature=0.0,
-               smearing=nothing, spin_polarization=:none) where {T <: Real}
+"""
+    Model(lattice; n_electrons, atoms, terms, temperature,
+                   smearing, spin_polarization, symmetry)
+
+Creates the physical specification of a model (without any
+discretization information).
+
+`n_electrons` is taken from `atoms` if not specified
+
+`spin_polarization` is :none by default (paired electrons)
+
+`smearing` is Fermi-Dirac if `temperature` is non-zero, none otherwise
+
+The `symmetry` kwarg can be:
+- auto: determine from terms if they respect the symmetry of the lattice/atoms.
+- off: no symmetries at all
+- force: force all the symmetries of the lattice/atoms.
+Careful that in this last case, wrong results can occur if the
+external potential breaks symmetries (this is not checked).
+"""
+function Model(lattice::AbstractMatrix{T};
+               n_electrons=nothing,
+               atoms=[],
+               terms=[],
+               temperature=T(0.0),
+               smearing=nothing,
+               spin_polarization=:none,
+               symmetry=:auto
+               ) where {T <: Real}
+
     lattice = Mat3{T}(lattice)
 
     if n_electrons === nothing
@@ -85,11 +116,24 @@ function Model(lattice::AbstractMatrix{T}; n_electrons=nothing, atoms=[], terms=
         error("Having several terms of the same name is not supported.")
     end
 
+    @assert symmetry in (:auto, :force, :off)
+    if symmetry == :auto
+        # ask the terms if they break symmetry
+        compute_symmetry = !(any(breaks_symmetries, terms))
+    else
+        compute_symmetry = (symmetry == :force)
+    end
+
+    if compute_symmetry
+        symops = symmetry_operations(lattice, atoms)
+    else
+        symops = [identity_symop()]
+    end
+
     Model{T}(lattice, recip_lattice, unit_cell_volume, recip_cell_volume, d, n_electrons,
-             spin_polarization, T(temperature), smearing, atoms, terms)
+             spin_polarization, T(temperature), smearing, atoms, terms, symops)
 end
 
-
 """
 Convenience constructor, which builds a standard atomic (kinetic + atomic potential) model.
 Use `extra_terms` to add additional terms.

src/PlaneWaveBasis.jl

@@ -6,28 +6,6 @@ include("fft.jl")
 # products of orbitals. This also defines the real-space grid
 # (as the dual of the cubic basis set).
 
-# The full (reducible) Brillouin zone is implicitly represented by
-# a set of (irreducible) kpoints (see explanation in docs/). Each
-# irreducible kpoint k comes with a list of symmetry operations
-# (S,τ) (containing at least the trivial operation (I,0)), where S
-# is a rotation matrix (/!\ not unitary in reduced coordinates)
-# and τ a translation vector. The kpoint is then used to represent
-# implicitly the information at all the kpoints Sk. The
-# relationship between the Hamiltonians is
-# H_{Sk} = U H_k U*, with
-# (Uu)(x) = u(S^-1(x-τ))
-# or in Fourier space
-# (Uu)(G) = e^{-i G τ} u(S^-1 G)
-# In particular, we can choose the eigenvectors at Sk as u_{Sk} = U u_k
-
-# We represent then the BZ as a set of irreducible points `kpoints`, a
-# list of symmetry operations `ksymops` allowing the reconstruction of
-# the full (reducible) BZ, and a set of weights `kweights` (summing to
-# 1). The value of an observable (eg energy) per unit cell is given as
-# the value of that observable at each irreducible kpoint, weighted by
-# kweight
-
-
 """
 Discretization information for kpoint-dependent quantities such as orbitals.
 More generally, a kpoint is a block of the Hamiltonian;
@@ -70,11 +48,12 @@ struct PlaneWaveBasis{T <: Real}
 
     # irreducible kpoints
     kpoints::Vector{Kpoint}
-    # Brillouin zone integration weights,
+    # BZ integration weights, summing up to 1
     # kweights[ik] = length(ksymops[ik]) / sum(length(ksymops[ik]) for ik=1:Nk)
     kweights::Vector{T}
     # ksymops[ikpt] is a list of symmetry operations (S,τ)
-    ksymops::Vector{Vector{Tuple{Mat3{Int}, Vec3{T}}}}
+    # mapping to points in the reducible BZ
+    ksymops::Vector{Vector{SymOp}}
 
     # fft_size defines both the G basis on which densities and
     # potentials are expanded, and the real-space grid
@@ -89,7 +68,13 @@ struct PlaneWaveBasis{T <: Real}
     # Instantiated terms (<: Term), that contain a backreference to basis.
     # See Hamiltonian for high-level usage
     terms::Vector{Any}
+
+    # symmetry operations that leave the reducible Brillouin zone invariant.
+    # Subset of model.symops, and superset of all the ksymops.
+    # Independent of the `use_symmetry` option
+    symops::Vector{SymOp}
 end
+
 # Default printing is just too verbose
 Base.show(io::IO, basis::PlaneWaveBasis) =
     print(io, "PlaneWaveBasis (Ecut=$(basis.Ecut), $(length(basis.kpoints)) kpoints)")
@@ -124,8 +109,9 @@ end
 build_kpoints(basis::PlaneWaveBasis, kcoords) =
     build_kpoints(basis.model, basis.fft_size, kcoords, basis.Ecut)
 
+# This is the "internal" constructor; the higher-level one below should be preferred
 function PlaneWaveBasis(model::Model{T}, Ecut::Number,
-                        kcoords::AbstractVector, ksymops;
+                        kcoords::AbstractVector, ksymops, symops=nothing;
                         fft_size=determine_grid_size(model, Ecut)) where {T <: Real}
     @assert Ecut > 0
     fft_size = Tuple{Int, Int, Int}(fft_size)
@@ -158,9 +144,16 @@ function PlaneWaveBasis(model::Model{T}, Ecut::Number,
 
     terms = Vector{Any}(undef, length(model.term_types))
 
+    if symops === nothing
+        # TODO instead compute the group generated by with ksymops, or
+        # just retire this constructor. Not critical because this
+        # should not be used in this context anyway...
+        symops = vcat(ksymops...)
+    end
+
     basis = PlaneWaveBasis{T}(
         model, Ecut, build_kpoints(model, fft_size, kcoords, Ecut),
-        kweights, ksymops, fft_size, opFFT, ipFFT, opIFFT, ipIFFT, terms)
+        kweights, ksymops, fft_size, opFFT, ipFFT, opIFFT, ipIFFT, terms, symops)
 
     # Instantiate terms
     for (it, t) in enumerate(model.term_types)
@@ -173,8 +166,8 @@ end
 Creates a new basis identical to `basis`, but with a different set of kpoints
 """
 function PlaneWaveBasis(basis::PlaneWaveBasis, kcoords::AbstractVector,
-                        ksymops::AbstractVector)
-    PlaneWaveBasis(basis.model, basis.Ecut, kcoords, ksymops;
+                        ksymops::AbstractVector, symops=nothing)
+    PlaneWaveBasis(basis.model, basis.Ecut, kcoords, ksymops, symops;
                    fft_size=basis.fft_size)
 end
 
@@ -183,20 +176,23 @@ end
 Creates a `PlaneWaveBasis` using the kinetic energy cutoff `Ecut` and a Monkhorst-Pack
 kpoint grid `kgrid` shifted by `kshift` (0 or 1/2 in each direction).
 
-If `enable_bzmesh_symmetry` is `true` (default) the symmetries of the
+If `use_symmetry` is `true` (default) the symmetries of the
 crystal are used to reduce the number of ``k``-Points which are
 treated explicitly. In this case all guess densities and potential
 functions must agree with the crystal symmetries or the result is
 undefined.
 """
 function PlaneWaveBasis(model::Model, Ecut::Number; kgrid=[1, 1, 1], kshift=[0, 0, 0],
-                        enable_bzmesh_symmetry=true, kwargs...)
-    if enable_bzmesh_symmetry
-        kcoords, ksymops = bzmesh_ir_wedge(kgrid, model.lattice, model.atoms, kshift=kshift)
+                        use_symmetry=true, kwargs...)
+    if use_symmetry
+        kcoords, ksymops, symops = bzmesh_ir_wedge(kgrid, model.symops, kshift=kshift)
     else
-        kcoords, ksymops = bzmesh_uniform(kgrid, kshift=kshift)
+        kcoords, ksymops, _ = bzmesh_uniform(kgrid, kshift=kshift)
+        # even when not using symmetry to reduce computations, still
+        # store in symops the set of kgrid-preserving symops
+        symops = symops_preserving_kgrid(model.symops, kcoords)
     end
-    PlaneWaveBasis(model, Ecut, kcoords, ksymops; kwargs...)
+    PlaneWaveBasis(model, Ecut, kcoords, ksymops, symops; kwargs...)
 end
 
 """
@@ -361,8 +357,8 @@ Convert a `basis` into one that uses or doesn't use BZ symmetrization
 Mainly useful for debug purposes (e.g. in cases we don't want to
 bother with symmetry)
 """
-function PlaneWaveBasis(basis::PlaneWaveBasis; enable_bzmesh_symmetry)
-    enable_bzmesh_symmetry && error("Not implemented")
+function PlaneWaveBasis(basis::PlaneWaveBasis; use_symmetry)
+    use_symmetry && error("Not implemented")
     if all(s -> length(s) == 1, basis.ksymops)
         return basis
     end
@@ -373,6 +369,6 @@ function PlaneWaveBasis(basis::PlaneWaveBasis; enable_bzmesh_symmetry)
         end
     end
     new_basis = PlaneWaveBasis(basis.model, basis.Ecut, kcoords,
-                               [[(Mat3{Int}(I), Vec3(zeros(3)))] for _ in 1:length(kcoords)];
+                               [[identity_symop()] for _ in 1:length(kcoords)];
                                fft_size=basis.fft_size)
 end

src/bzmesh.jl

@@ -34,7 +34,7 @@ operations are the identity.
 """
 function bzmesh_uniform(kgrid_size; kshift=[0, 0, 0])
     kcoords = kgrid_monkhorst_pack(kgrid_size; kshift=kshift)
-    kcoords, [[(Mat3{Int}(I), Vec3(zeros(3)))] for _ in 1:length(kcoords)]
+    kcoords, [[identity_symop()] for _ in 1:length(kcoords)], [identity_symop()]
 end
 
 
@@ -48,8 +48,7 @@ the full mesh. `lattice` are the lattice vectors, column by column, `atoms` are
 representing a mapping from `Element` objects to a list of positions in fractional
 coordinates. `tol_symmetry` is the tolerance used for searching for symmetry operations.
 """
-function bzmesh_ir_wedge(kgrid_size, lattice, atoms;
-                         tol_symmetry=1e-5, kshift=[0, 0, 0])
+function bzmesh_ir_wedge(kgrid_size, symops; kshift=[0, 0, 0])
     all(isequal.(kgrid_size, 1)) && return bzmesh_uniform(kgrid_size, kshift=kshift)
 
     # Transform kshift to the convention used in spglib:
@@ -60,8 +59,8 @@ function bzmesh_ir_wedge(kgrid_size, lattice, atoms;
 
     kpoints_mp = kgrid_monkhorst_pack(kgrid_size, kshift=kshift)
 
-    # Get the list of symmetry operations (S,τ) that preserve the MP grid
-    symops = symmetry_operations(lattice, atoms; kcoords=kpoints_mp)
+    # Filter those symmetry operations (S,τ) that preserve the MP grid
+    symops = symops_preserving_kgrid(symops, kpoints_mp)
 
     # Give the remaining symmetries to spglib to compute an irreducible k-Point mesh
     # TODO implement time-reversal symmetry and turn the flag to true
@@ -132,7 +131,7 @@ function bzmesh_ir_wedge(kgrid_size, lattice, atoms;
     @assert all(findfirst(Sτ -> iszero(Sτ[1] - I) && iszero(Sτ[2]), ops) !== nothing
                 for ops in ksymops)
 
-    kirreds, ksymops
+    kirreds, ksymops, symops
 end
 
 

src/common/types.jl

@@ -4,3 +4,6 @@ using StaticArrays
 const Mat3{T} = SMatrix{3, 3, T, 9} where T
 const Vec3{T} = SVector{3, T} where T
 const AbstractArray3{T} = AbstractArray{T, 3}
+# Represents a symmetry (S,τ)
+const SymOp = Tuple{Mat3{Int}, Vec3{Float64}}
+identity_symop() = (Mat3{Int}(I), Vec3(zeros(3)))

src/external/spglib.jl

@@ -7,7 +7,7 @@ in the convention required to take the place of a `cell` datastructure used in s
 """
 function spglib_cell_atommapping_(lattice, atoms)
     lattice = Matrix{Float64}(lattice)  # spglib operates in double precision
-    n_attypes = sum(length(positions) for (type, positions) in atoms)
+    n_attypes = isempty(atoms) ? 0 : sum(length(positions) for (type, positions) in atoms)
     spg_numbers = Vector{Int}(undef, n_attypes)
     spg_positions = Matrix{Float64}(undef, n_attypes, 3)
 
@@ -37,14 +37,22 @@ function spglib_get_symmetry(lattice, atoms; tol_symmetry=1e-5)
     spglib = pyimport("spglib")
     lattice = Matrix{Float64}(lattice)  # spglib operates in double precision
 
+    if isempty(atoms)
+        # spglib doesn't like no atoms, so we default to
+        # no symmetries (even though there are lots)
+        return [Mat3{Int}(I)], [Vec3(zeros(3))]
+    end
+
     # Ask spglib for symmetry operations and for irreducible mesh
     spg_symops = spglib.get_symmetry(spglib_cell(lattice, atoms),
                                      symprec=tol_symmetry)
 
     # If spglib does not find symmetries give an error
-    spg_symops !== nothing || error(
-        "spglib failed to get the symmetries. Check your lattice, or use a uniform BZ mesh."
-    )
+    if spg_symops === nothing
+        err_message=spglib.get_error_message()
+        error("spglib failed to get the symmetries. Check your lattice, use a " *
+              "uniform BZ mesh or disable symmetries. Spglib reported : " * err_message)
+    end
 
     Stildes = [St for St in eachslice(spg_symops["rotations"]; dims=1)]
     τtildes = [rationalize.(τt, tol=tol_symmetry)

src/postprocess/band_structure.jl

@@ -23,7 +23,7 @@ end
 function compute_bands(basis, ρ, kcoords, n_bands;
                        eigensolver=lobpcg_hyper, tol=1e-3, show_progress=true, kwargs...)
     # Create basis with new kpoints, where we cheat by using any symmetry operations.
-    ksymops = [[(Mat3{Int}(I), Vec3(zeros(3)))] for _ in 1:length(kcoords)]
+    ksymops = [[identity_symop()] for _ in 1:length(kcoords)]
     # For some reason rationalize(2//3) isn't supported (julia 1.4)
     myrationalize(x::T) where {T <: AbstractFloat} = rationalize(x, tol=10eps(T))
     myrationalize(x) = x