crystal.jl

using LinearAlgebra

###########################################
# Types and basis 
abstract type AbstractCrystalGroup end

struct CrystalGroup{K} <: AbstractCrystalGroup
    R::Array{K, 2}
    T::Array{K, 1}
end

CrystalGroup(R::Array{K, 2}) where K = CrystalGroup{K}(R, zeros(K, 3))
CrystalGroup(T::Array{K, 1}) where K = CrystalGroup{K}(Matrix{K}(I, 3, 3), T)

Base.:*(g1::CrystalGroup, g2::CrystalGroup) = CrystalGroup(g1.R*g2.R, g1.R*g2.T + g1.T)
Base.:*(λ::K1, g::CrystalGroup) where {K1 <: Number} = CrystalGroup(λ * g.R, λ * g.T)
Base.:*(g::CrystalGroup, v::Vector) = g.R * v + g.T

r1 = CrystalGroup(Matrix{Float64}(I, 3, 3))
r2 = CrystalGroup([0.0 -1 0; 1 -1 0; 0 0 1])
r3 = CrystalGroup([-1.0 1 0; -1 0 0; 0 0 1])
r4 = CrystalGroup([0.0 1 0; 1 0 0; 0 0 -1])
r5 = CrystalGroup([1.0 -1 0; 0 -1 0; 0 0 -1])
r6 = CrystalGroup([-1.0 0 0; -1 1 0; 0 0 -1])

t1 = CrystalGroup([0.0, 0, 0])
t2 = CrystalGroup([0, 0, 0.5])
t3 = CrystalGroup(([0, 0, 0] + [2/3, 1/3, 1/3]) - floor.([0, 0, 0] + [2/3, 1/3, 1/3]))
t4 = CrystalGroup(([0, 0, 0.5] + [2/3, 1/3, 1/3]) - floor.([0, 0, 0.5] + [2/3, 1/3, 1/3]))
t5 = CrystalGroup(([0, 0, 0] + [1/3, 2/3, 2/3]) - floor.([0, 0, 0] + [1/3, 2/3, 2/3]))
t6 = CrystalGroup(([0, 0, 0.5] + [1/3, 2/3, 2/3]) - floor.([0, 0, 0.5] + [1/3, 2/3, 2/3]))

rbasis = [[r1, r2, r3], [r4, r5, r6]]
tbasis = [[t1, t3, t5], [t2, t4, t6]]

basis = [[[t * r, t * (-1 * r)] for (r, t) in Iterators.product(rs, ts)] for (rs, ts) in zip(rbasis, tbasis)]
basis = [x for x in Iterators.flatten(basis)]
basis = [x for x in Iterators.flatten(basis)]

############################################################
# Constants and auxiliary functions

is_constrained(v, a, b, c) = (-a < v[1] < a) && (-b < v[2] < b) && (-c < v[3] < c)
a = 1
b = 1
c = 1
alpha = 90
beta = alpha
gamma = 120
a1 = 4.7606
b1 = a1
c1 = 12.994
abc = [a1 b1 c1]

M = [a b*cosd(gamma) c*cosd(beta); 0 b*sind(gamma) (c/sind(gamma))*(cosd(alpha)-cosd(beta)*cosd(gamma)); 0 0 (c*((1-cosd(alpha)^2 - cosd(beta)^2 - cosd(gamma)^2 + (2)*cosd(alpha)* cosd(beta)* cosd(gamma))^0.5))/(sind(gamma))]
M1 = M .* abc

Al = [0.0 0.0 0.35217]
O = [0.69365 0.0 0.25]
x = Al[1,1]
y = Al[1,2]
z = Al[1,3]
v = [x, y, z]

#################################################################
# Calculations

function mult_em_all(v, basis, a1, b1, c1, M, n)
    vs = [M * (g * v) for g in basis]
    atoms = Set{Vector{Float64}}()
    for i in -n:n, j in -n:n, k in -n:n
        for v1 in vs
            p = v1 + M * [i, j, k]
            is_constrained(p, a1, b1, c1) && push!(atoms, p)
        end
    end

    atoms
end

mult_em_all(v, basis, a1, b1, c1, M1, 50)