Lattice symmetries

ext/DFTKJLD2Ext.jl

@@ -42,6 +42,7 @@ function load_basis(jld)
                       jld["variational"],
                       jld["kgrid"],
                       jld["symmetries_respect_rgrid"],
+                      jld["symmetries_respect_basis"],
                       jld["use_symmetries_for_kpoint_reduction"])
     else
         basis_args = nothing

src/PlaneWaveBasis.jl

@@ -106,10 +106,13 @@ struct PlaneWaveBasis{T,
     ## Symmetry operations that leave the discretized model (k and r grids) invariant.
     # Subset of model.symmetries.
     symmetries::Vector{SymOp{VT}}
-    # Whether the symmetry operations leave the rgrid invariant
+    # Whether the symmetry operations leave the rgrid invariant.
     # If this is true, the symmetries are a property of the complete discretized model.
     # Therefore, all quantities should be symmetric to machine precision
     symmetries_respect_rgrid::Bool
+    # Whether the symmetry operations leave the basis invariant. Otherwise the symmetry
+    # operations will leave the lattice invariant.
+    symmetries_respect_basis::Bool
     # Whether symmetry is used to reduce the number of explicit k-points to the
     # irreducible BZMesh. This is a debug option, useful when a part in the code does
     # not yet implement symmetry. See `unfold_bz` as a convenient way to use this.
@@ -175,6 +178,7 @@ end
 function PlaneWaveBasis(model::Model{T}, Ecut::Real, fft_size::Tuple{Int, Int, Int},
                         variational::Bool, kgrid::AbstractKgrid,
                         symmetries_respect_rgrid::Bool,
+                        symmetries_respect_basis::Bool,
                         use_symmetries_for_kpoint_reduction::Bool,
                         comm_kpts, architecture::Arch
                        ) where {T <: Real, Arch <: AbstractArchitecture}
@@ -209,7 +213,11 @@ function PlaneWaveBasis(model::Model{T}, Ecut::Real, fft_size::Tuple{Int, Int, I
     if symmetries_respect_rgrid
         symmetries = symmetries_preserving_rgrid(symmetries, fft_size)
     end
-    symmetries = symmetries_preserving_kgrid(symmetries, kgrid)
+    # If false, symmetries will respect the model (leave the lattice invariant), as is the
+    # default behaviour for Abinit (nsym 0).
+    if symmetries_respect_basis
+        symmetries = symmetries_preserving_kgrid(symmetries, kgrid)
+    end
 
     # Build the irreducible k-point coordinates
     if use_symmetries_for_kpoint_reduction
@@ -317,7 +325,7 @@ function PlaneWaveBasis(model::Model{T}, Ecut::Real, fft_size::Tuple{Int, Int, I
         kpoints, kweights, kgrid,
         kcoords_global, kweights_global,
         comm_kpts, krange_thisproc, krange_allprocs, krange_thisproc_allspin,
-        architecture, symmetries, symmetries_respect_rgrid,
+        architecture, symmetries, symmetries_respect_rgrid, symmetries_respect_basis,
         use_symmetries_for_kpoint_reduction, terms)
 
     # Instantiate the terms with the basis
@@ -342,6 +350,7 @@ with a maximal spacing of `2π * 0.022` per Bohr.
                                 kshift=[0, 0, 0],
                                 variational=true, fft_size=nothing,
                                 symmetries_respect_rgrid=isnothing(fft_size),
+                                symmetries_respect_basis=true,
                                 use_symmetries_for_kpoint_reduction=true,
                                 comm_kpts=MPI.COMM_WORLD, architecture=CPU()) where {T <: Real}
     if isnothing(fft_size)
@@ -372,8 +381,8 @@ with a maximal spacing of `2π * 0.022` per Bohr.
     end
 
     PlaneWaveBasis(model, austrip(Ecut), fft_size, variational, kgrid_inner,
-                   symmetries_respect_rgrid, use_symmetries_for_kpoint_reduction,
-                   comm_kpts, architecture)
+                   symmetries_respect_rgrid, symmetries_respect_basis,
+                   use_symmetries_for_kpoint_reduction, comm_kpts, architecture)
 end
 
 """
@@ -383,7 +392,7 @@ e.g. an [`MonkhorstPack`](@ref) or a [`ExplicitKpoints`](@ref) grid.
 @timing function PlaneWaveBasis(basis::PlaneWaveBasis, kgrid::AbstractKgrid)
     PlaneWaveBasis(basis.model, basis.Ecut,
                    basis.fft_size, basis.variational,
-                   kgrid, basis.symmetries_respect_rgrid,
+                   kgrid, basis.symmetries_respect_rgrid, basis.symmetries_respect_basis,
                    basis.use_symmetries_for_kpoint_reduction,
                    basis.comm_kpts, basis.architecture)
 end
@@ -544,6 +553,7 @@ function gather_kpts(basis::PlaneWaveBasis)
                        basis.variational,
                        basis.kgrid,
                        basis.symmetries_respect_rgrid,
+                       basis.symmetries_respect_basis,
                        basis.use_symmetries_for_kpoint_reduction,
                        MPI.COMM_SELF,
                        basis.architecture)

src/input_output.jl

@@ -157,6 +157,7 @@ function todict!(dict, basis::PlaneWaveBasis)
     dict["Ecut"]         = basis.Ecut
     dict["variational"]  = basis.variational
     dict["symmetries_respect_rgrid"] = basis.symmetries_respect_rgrid
+    dict["symmetries_respect_basis"] = basis.symmetries_respect_basis
     dict["use_symmetries_for_kpoint_reduction"] = basis.use_symmetries_for_kpoint_reduction
 
     # Update the symmetry as discretisation might have broken some symmetries

src/postprocess/stresses.jl

@@ -23,6 +23,7 @@ for details. In Voigt notation one would use the vector
         new_basis = PlaneWaveBasis(new_model,
                                    basis.Ecut, basis.fft_size, basis.variational,
                                    basis.kgrid, basis.symmetries_respect_rgrid,
+                                   basis.symmetries_respect_basis,
                                    basis.use_symmetries_for_kpoint_reduction,
                                    basis.comm_kpts, basis.architecture)
         ρ = compute_density(new_basis, scfres.ψ, scfres.occupation)

src/supercell.jl

@@ -44,9 +44,11 @@ function cell_to_supercell(basis::PlaneWaveBasis{T}) where {T}
     supercell_model = Model(model; supercell.lattice, supercell.atoms, supercell.positions,
                             n_electrons=n_electrons_supercell, symmetries=false)
     symmetries_respect_rgrid = true
+    symmetries_respect_basis = true
     use_symmetries_for_kpoint_reduction = true
     PlaneWaveBasis(supercell_model, basis.Ecut, supercell_fft_size,
-                   basis.variational, supercell_kgrid, symmetries_respect_rgrid,
+                   basis.variational, supercell_kgrid,
+                   symmetries_respect_rgrid, symmetries_respect_basis,
                    use_symmetries_for_kpoint_reduction,
                    basis.comm_kpts, basis.architecture)
 end

src/symmetry.jl

@@ -387,6 +387,7 @@ function unfold_bz(basis::PlaneWaveBasis)
         return PlaneWaveBasis(basis.model, basis.Ecut, basis.fft_size,
                               basis.variational, basis.kgrid,
                               basis.symmetries_respect_rgrid,
+                              basis.symmetries_respect_basis,
                               use_symmetry_for_kpoint_reduction,
                               basis.comm_kpts, basis.architecture)
     end

src/workarounds/forwarddiff_rules.jl

@@ -180,6 +180,7 @@ function construct_value(basis::PlaneWaveBasis{T}) where {T <: ForwardDiff.Dual}
                    basis.variational,
                    basis.kgrid,
                    basis.symmetries_respect_rgrid,
+                   basis.symmetries_respect_basis,
                    basis.use_symmetries_for_kpoint_reduction,
                    basis.comm_kpts,
                    basis.architecture)

test/todict.jl

@@ -53,6 +53,7 @@ function test_agreement_bands(band_data, dict; explicit_reshape=false, test_ψ=t
         @test dict["symmetries_translations"] ≈ translations atol=1e-12
         @test dict["use_symmetries_for_kpoint_reduction"] == basis.use_symmetries_for_kpoint_reduction
         @test dict["symmetries_respect_rgrid"] == basis.symmetries_respect_rgrid
+        @test dict["symmetries_respect_basis"] == basis.symmetries_respect_basis
 
         lattice_resh = condreshape(dict["lattice"], 3, 3)
         rotations_resh = [condreshape(rot, 3, 3) for rot in dict["symmetries_rotations"]]