Add and use cis2pi functions (JuliaMolSim#641)

src/DFTK.jl

@@ -27,6 +27,7 @@ include("common/split_evenly.jl")
 include("common/mpi.jl")
 include("common/threading.jl")
 include("common/printing.jl")
+include("common/cis2pi.jl")
 
 export PspHgh
 include("pseudo/NormConservingPsp.jl")

src/common/cis2pi.jl

@@ -0,0 +1,8 @@
+"""Function to compute exp(2π i x)"""
+cis2pi(x) = cispi(2x)
+
+"""Function to compute sin(2π x)"""
+sin2pi(x) = sinpi(2x)
+
+"""Function to compute cos(2π x)"""
+cos2pi(x) = cospi(2x)

src/fft.jl

@@ -112,7 +112,7 @@ function G_to_r_matrix(basis::PlaneWaveBasis{T}) where {T}
     ret = zeros(complex(T), prod(basis.fft_size), prod(basis.fft_size))
     for (iG, G) in enumerate(G_vectors(basis))
         for (ir, r) in enumerate(r_vectors(basis))
-            ret[ir, iG] = cis(2π * dot(r, G)) / sqrt(basis.model.unit_cell_volume)
+            ret[ir, iG] = cis2pi(dot(r, G)) / sqrt(basis.model.unit_cell_volume)
         end
     end
     ret
@@ -122,7 +122,7 @@ function r_to_G_matrix(basis::PlaneWaveBasis{T}) where {T}
     for (iG, G) in enumerate(G_vectors(basis))
         for (ir, r) in enumerate(r_vectors(basis))
             Ω = basis.model.unit_cell_volume
-            ret[iG, ir] = cis(-2π * dot(r, G)) * sqrt(Ω) / prod(basis.fft_size)
+            ret[iG, ir] = cis2pi(-dot(r, G)) * sqrt(Ω) / prod(basis.fft_size)
         end
     end
     ret

src/guess_density.jl

@@ -104,7 +104,7 @@ function gaussian_superposition(basis::PlaneWaveBasis{T}, gaussians) where {T}
         Gsq = sum(abs2, basis.model.recip_lattice * G)
         for (coeff, decay_length, r) in gaussians
             form_factor::T = exp(-Gsq * T(decay_length)^2)
-            ρ[iG] += T(coeff) * form_factor * cis(-2T(π) * dot(G, r))
+            ρ[iG] += T(coeff) * form_factor * cis2pi(-dot(G, r))
         end
     end
 

src/symmetry.jl

@@ -121,7 +121,7 @@ function apply_symop(symop::SymOp, basis, kpoint, ψk::AbstractVecOrMat)
         for (ig, G_full) in enumerate(Gs_full)
             igired = index_G_vectors(basis, kpoint, invS * G_full)
             @assert igired !== nothing
-            ψSk[ig, iband] = cis(-2π * dot(G_full, τ)) * ψk[igired, iband]
+            ψSk[ig, iband] = cis2pi(-dot(G_full, τ)) * ψk[igired, iband]
         end
     end
 
@@ -160,7 +160,7 @@ function accumulate_over_symmetries!(ρaccu, ρin, basis, symmetries)
         for (ig, G) in enumerate(G_vectors_generator(basis.fft_size))
             igired = index_G_vectors(basis, invS * G)
             if igired !== nothing
-                @inbounds ρaccu[ig] += cis(-2T(π) * T(dot(G, symop.τ))) * ρin[igired]
+                @inbounds ρaccu[ig] += cis2pi(-T(dot(G, symop.τ))) * ρin[igired]
             end
         end
     end  # symop

src/terms/ewald.jl

@@ -113,8 +113,8 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing, fo
                 continue
             end
 
-            cos_strucfac = sum(Z * cos(2T(π) * dot(r, G)) for (r, Z) in zip(positions, charges))
-            sin_strucfac = sum(Z * sin(2T(π) * dot(r, G)) for (r, Z) in zip(positions, charges))
+            cos_strucfac = sum(Z * cos2pi(dot(r, G)) for (r, Z) in zip(positions, charges))
+            sin_strucfac = sum(Z * sin2pi(dot(r, G)) for (r, Z) in zip(positions, charges))
             sum_strucfac = cos_strucfac^2 + sin_strucfac^2
 
             any_term_contributes = true
@@ -123,8 +123,8 @@ function energy_ewald(lattice, recip_lattice, charges, positions; η=nothing, fo
             if forces !== nothing
                 for (ir, r) in enumerate(positions)
                     Z = charges[ir]
-                    dc = -Z*2T(π)*G*sin(2T(π) * dot(r, G))
-                    ds = +Z*2T(π)*G*cos(2T(π) * dot(r, G))
+                    dc = -Z*2T(π)*G*sin2pi(dot(r, G))
+                    ds = +Z*2T(π)*G*cos2pi(dot(r, G))
                     dsum = 2cos_strucfac*dc + 2sin_strucfac*ds
                     forces_recip[ir] -= dsum * exp(-exponent)/Gsq
                 end

src/terms/local.jl

@@ -78,7 +78,7 @@ function (::AtomicLocal)(basis::PlaneWaveBasis{T}) where {T}
         pot = sum(model.atom_groups) do group
             element = model.atoms[first(group)]
             form_factor::T = local_potential_fourier(element, norm(model.recip_lattice * G))
-            form_factor * sum(cis(-2T(π) * dot(G, r)) for r in @view model.positions[group])
+            form_factor * sum(cis2pi(-dot(G, r)) for r in @view model.positions[group])
         end
         pot / sqrt(model.unit_cell_volume)
     end
@@ -116,7 +116,7 @@ function _force_local_internal(basis, ρ_fourier, form_factors, r)
     for (iG, G) in enumerate(G_vectors(basis))
         f -= real(conj(ρ_fourier[iG])
                   .* form_factors[iG]
-                  .* cis(-2T(π) * dot(G, r))
+                  .* cis2pi(-dot(G, r))
                   .* (-2T(π)) .* G .* im
                   ./ sqrt(basis.model.unit_cell_volume))
     end

src/terms/nonlocal.jl

@@ -66,7 +66,7 @@ end
             form_factors = build_form_factors(element.psp, qs_cart)
             for idx in group
                 r = model.positions[idx]
-                structure_factors = [cis(-2T(π) * dot(q, r)) for q in qs]
+                structure_factors = [cis2pi(-dot(q, r)) for q in qs]
                 P = structure_factors .* form_factors ./ sqrt(unit_cell_volume)
 
                 forces[idx] += map(1:3) do α
@@ -160,7 +160,7 @@ function build_projection_vectors_(basis::PlaneWaveBasis{T}, kpt::Kpoint,
         # Combine with structure factors
         for r in positions
             # k+G in this formula can also be G, this only changes an unimportant phase factor
-            structure_factors = map(q -> cis(-2T(π) * dot(q, r)), Gplusk_vectors(basis, kpt))
+            structure_factors = map(q -> cis2pi(-dot(q, r)), Gplusk_vectors(basis, kpt))
             @views for iproj = 1:count_n_proj(psp)
                 proj_vectors[:, offset+iproj] .= (
                     structure_factors .* form_factors[:, iproj] ./ sqrt(unit_cell_volume)

src/workarounds/intervals.jl

@@ -7,6 +7,11 @@ const Interval = IntervalArithmetic.Interval
 # should be done e.g. by changing  the rounding mode ...
 erfc(i::Interval) = Interval(prevfloat(erfc(i.lo)), nextfloat(erfc(i.hi)))
 
+# This is done to avoid using sincospi(x), called by cispi(x),
+# which has not been implemented in IntervalArithmetic
+# see issue #513 on IntervalArithmetic repository
+cis2pi(x::Interval) = exp(2 * (pi * (im * x)))
+
 function compute_Glims_fast(lattice::AbstractMatrix{<:Interval}, args...; kwargs...)
     # This is done to avoid a call like ceil(Int, ::Interval)
     # in the above implementation of compute_fft_size,