Zhonghua Ma and Mark Tuckerman, J. Chem. Phys 133 184110 (2010)
Three dimensional FBR basis, defined in terms of general cell matrix $\mathbf{h}$.
Consider an invertible transformation between:
The columns of $\mathbf{h}$ are three box vectors $\mathbf{a}$, $\mathbf{b}$ and $\mathbf{c}$.
The above coordinate transformation maps a simulation cell of shape $\mathbf{h}$ into a unit cubic box.
Let's demonstrate this equation for several lattice constant.
We will generate points in scaled coordinate $s_{\beta}$ and get the actual grid points by multiplying them with cell matrix $h_{\alpha\beta}$
In [1]:
function gen_lattice_hexagonal(a; coa=8.0/3.0)
LL = zeros(3,3)
LL[1,:] = [1.0, 0.0, 0.0]
LL[2,:] = [-cos(pi/3.0), sin(pi/3.0), 0.0]
LL[3,:] = [0.0, 0.0, coa]
return a*LL
end
function gen_lattice_fcc(a)
LL = zeros(3,3)
LL[1,:] = [ 0.0, 1.0, 1.0]
LL[2,:] = [ 1.0, 0.0, 1.0]
LL[3,:] = [ 1.0, 1.0, 0.0]
return 0.5*a*LL
end
Out[1]:
In [3]:
LL = gen_lattice_fcc(1.0)
a = LL[1,:]
b = LL[2,:]
c = LL[3,:]
LL
Out[3]:
Define general cell matrix $\mathbf{h}$
In [11]:
hh = zeros(3,3)
hh[:,1] = a
hh[:,2] = b
hh[:,3] = c
hh
Out[11]:
An old function function that I used to generate directly real space grid points
In [2]:
function init_r_grids( Ns, LatVecs )
#
Npoints = prod(Ns)
#
r = Array{Float64}(3,Npoints)
ip = 0
for k in 0:Ns[3]-1
for j in 0:Ns[2]-1
for i in 0:Ns[1]-1
ip = ip + 1
r[1,ip] = LatVecs[1,1]*i/Ns[1] + LatVecs[2,1]*j/Ns[2] + LatVecs[3,1]*k/Ns[3]
r[2,ip] = LatVecs[1,2]*i/Ns[1] + LatVecs[2,2]*j/Ns[2] + LatVecs[3,2]*k/Ns[3]
r[3,ip] = LatVecs[1,3]*i/Ns[1] + LatVecs[2,3]*j/Ns[2] + LatVecs[3,3]*k/Ns[3]
end
end
end
return r
end
Out[2]:
In [3]:
function write_xsf( filnam, LL, atpos; atsymb=nothing, molecule=false )
#
f = open(filnam, "w")
Natoms = size(atpos)[2]
#
if molecule
@printf(f, "MOLECULE\n")
else
@printf(f, "CRYSTAL\n")
end
@printf(f, "PRIMVEC\n")
@printf(f, "%18.10f %18.10f %18.10f\n", LL[1,1], LL[1,2], LL[1,3])
@printf(f, "%18.10f %18.10f %18.10f\n", LL[2,1], LL[2,2], LL[2,3])
@printf(f, "%18.10f %18.10f %18.10f\n", LL[3,1], LL[3,2], LL[3,3])
@printf(f, "PRIMCOORD\n")
@printf(f, "%8d %8d\n", Natoms, 1)
#
if atsymb == nothing
for ia = 1:Natoms
@printf(f, "X %18.10f %18.10f %18.10f\n", atpos[1,ia], atpos[2,ia], atpos[3,ia])
end
else
for ia = 1:Natoms
@printf(f, "%s %18.10f %18.10f %18.10f\n", atsymb[ia], atpos[1,ia], atpos[2,ia], atpos[3,ia])
end
end
close(f)
end
Out[3]:
Let's do this first using my old function init_r_grids
In [25]:
Ns = [11,11,11]
r = init_r_grids(Ns,10*LL)
write_xsf("r_fcc.xsf", 10*LL, r)
Using the new way, first we need to generate grid of $s$:
In [14]:
function init_s_grid(Ns)
if any( (Ns .% 2) .== 0 )
error("Ns must be odd numbers")
end
NN = round.(Int, (Ns-1)/2)
Npoints = prod(Ns)
s = zeros(3,Npoints)
ip = 0
for n = -NN[3]:NN[3]
for m = -NN[2]:NN[2]
for l = -NN[1]:NN[1]
ip = ip + 1
s[1,ip] = l/(2*NN[1] + 1)
s[2,ip] = m/(2*NN[2] + 1)
s[3,ip] = n/(2*NN[3] + 1)
end
end
end
return s
end
Out[14]:
Finally, a function to do the transformation from $s$ to $r$
In [15]:
function s_to_r(s::Array{Float64,2}, hh)
Npoints = size(s)[2]
r = zeros(3,Npoints)
for ip = 1:Npoints
for α = 1:3
r[α,ip] = 0.0
for β = 1:3
r[α,ip] = r[α,ip] + hh[α,β]*s[β,ip]
end
end
end
return r
end
Out[15]:
Test for fcc and hexagonal grid:
In [ ]:
println(10*hh)
r = s_to_r(s, 10*hh)
write_xsf("r_fcc_v2.xsf", 10*LL, r)
In [58]:
LL = gen_lattice_hexagonal(10)
hh = LL'
# Remember that Ns must be odd numbers
Ns = [11,11,15]
s = init_s_grid(Ns)
r = s_to_r(s,hh)
# don't forget to transpose hh because our convention is different from the one
# used in Ma and Tuckerman
write_xsf("r_hcp_v1.xsf", hh', r)
Periodic boundary conditions are used, plane-wave like functions $\mathbf{\phi}_{\hat{\mathbf{k}}}(\mathbf{s})$ can be defined as functions of $s_{\alpha} \in [ -\dfrac{1}{2},\dfrac{1}{2} ]$ as
where $\hat{\mathbf{k}} = (\hat{k}_1, \hat{k}_2, \hat{k}_3)$ is a vector of integers.
The basis $\{ \phi_{\hat{\mathbf{k}}}(\mathbf{s}) \}$ are orthonormal:
We can define a FBR from the curvilinear plane-wave functions $\{\phi_{\hat{\mathbf{k}}}\}$, restricting:
where the size of DVR quadrature grid is $2\mathcal{N}+1$
The corresponding DVR functions $\mu_{lmn}(\mathbf{s})$ can be defined as
For periodic system, the scalar product $\left\langle \left.\phi_{\hat{\mathbf{k}}}\right|\mu_{lmn}\right\rangle$ can be evaluated by $\mathcal{N}$-point Chebyshev-Gaussian quadrature of the first kind. An equally spaced grid in the curvilinear coordinates ensures the accuracy of the quadrature.
The DVR points $\mathbf{s}_{lmn}$ are
The DVR grid points in Cartesian coordinates are given by the matrix-vector product:
The corresponding DVR functions $\mu_{lmn}(\mathbf{s})$ then become
In [24]:
# Function to evaluate LF
function eval_μ_lmn(Ns, hh, s_grid, ip, s)
#
NN = round.(Int, (Ns-1)/2)
metric = 1/sqrt(det(hh))
#
s_lmn = s_grid[:,ip]
#
f = 0.0 + im*0.0
# loop structure should be the same as the one used in init_s_grid
for kz = -NN[3]:NN[3]
for ky = -NN[2]:NN[2]
for kx = -NN[1]:NN[1]
ks = dot( [kx,ky,kz], (s - s_lmn) )
f = f + cos(2π*ks)
end
end
end
return metric/sqrt(prod(Ns))*f
end
Out[24]:
In [30]:
# Sampling points
Ns = [21,21,11]
#LL = gen_lattice_hexagonal(10)
LL = 10.0*eye(3,3) # cubic
hh = LL'
s_grid = init_s_grid(Ns)
ip = 2
#r = [5.0, 5.0, 5.0]
#s = inv(hh)*r
s = s_grid[:,1]
println( eval_μ_lmn(Ns, hh, s_grid, ip, s) )
s = s_grid[:,2]
println( eval_μ_lmn(Ns, hh, s_grid, ip, s) )
With the above definition of the DVR functions, the overlap integrals of the DVR functions are equal to the Kronecker delta, and the values of $\mu_{lmn}(\mathbf{s})$ are zero at all DVR points except at $\mathbf{s}_{lmn}$
Laplacian operator: $$ \begin{equation} \nabla^2 = \sum_{\alpha=1}^{3}\sum_{\beta=1}^{3}\sum_{\gamma=1}^{3} \frac{\partial}{\partial s_{\beta}} \mathbf{h}^{-1}_{\beta\alpha} \frac{\partial}{\partial s_{\gamma}} \mathbf{h}^{-1}_{\gamma\alpha} \end{equation} $$
The kinetic energy can be computed as $$ \begin{align} E_{\mathrm{kin}} & = -\frac{1}{2}\sum_{i=1}^{N_{s}} f_{i} \left. \left\langle \psi_{i} \right| \nabla^2 \left| \psi_{i} \right\rangle \right. \\ & = -\frac{1}{2}\sum_{i=1}^{N_{s}}f_{i} \sum_{l,l'=-\mathcal{N}}^{\mathcal{N}} \sum_{m,m'=-\mathcal{N}}^{\mathcal{N}} \sum_{n,n'=-\mathcal{N}}^{\mathcal{N}} C_{lmn}^{i*} C_{l'm'n'}^{i} \mathbb{L}_{lmn}^{l'm'n'} \end{align} $$
The matrix $\mathbb{L}_{lmn}^{l'm'n'}$ is given by $$ \begin{align} \mathbb{L}_{lmn}^{l'm'n'} & = \sum_{\alpha=1}^{3} \mathbf{h}_{1\alpha}^{-1}\mathbf{h}_{1\alpha}^{-1} D_{ll'}\delta_{mm'}\delta_{nn'} \\ & + \sum_{\alpha=1}^{3} \mathbf{h}_{2\alpha}^{-1}\mathbf{h}_{2\alpha}^{-1} D_{mm'}\delta_{ll'}\delta_{nn'} \\ & + \sum_{\alpha=1}^{3} \mathbf{h}_{3\alpha}^{-1}\mathbf{h}_{3\alpha}^{-1} D_{nn'}\delta_{ll'}\delta_{mm'} \end{align} $$
where $$ \begin{equation} D_{jj'} = -4\pi^2 \left[ \frac{\mathcal{N}(\mathcal{N}-1)}{3} \delta_{jj'} + \frac{(-1)^{j-j'} \cos\left[ \dfrac{\pi(j-j')}{2\mathcal{N}+1} \right] } {2\sin^2\left[ \dfrac{\pi(j-j')}{2\mathcal{N}+1} \right]} (1-\delta_{jj'}) \right] \end{equation} $$
In [ ]: