Constant pressure ab initio molecular dynamics with discrete variable representation basis sets

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:

a set of curvilinear coordinates $s_{\alpha}$, $\alpha = 1,2,3$, which has three axes along three box vectors of a general cell matrix, and
a set of Cartesian coordinates $r_{\alpha}$, $\alpha = 1,2,3$

$$ \begin{equation} r_{\alpha} = \sum_{j=1}^{3} h_{\alpha\beta}s_{\beta} \end{equation} $$

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]:





gen_lattice_fcc (generic function with 1 method)



In [3]:

    
LL = gen_lattice_fcc(1.0)
a = LL[1,:]
b = LL[2,:]
c = LL[3,:]
LL









    Out[3]:





3×3 Array{Float64,2}:
 0.0  0.5  0.5
 0.5  0.0  0.5
 0.5  0.5  0.0

Define general cell matrix $\mathbf{h}$



In [11]:

    
hh = zeros(3,3)
hh[:,1] = a
hh[:,2] = b
hh[:,3] = c
hh









    Out[11]:





3×3 Array{Float64,2}:
 0.0  0.5  0.5
 0.5  0.0  0.5
 0.5  0.5  0.0

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]:





init_r_grids (generic function with 1 method)



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]:





write_xsf (generic function with 1 method)

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]:





init_s_grid (generic function with 1 method)

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]:





s_to_r (generic function with 1 method)

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

$$ \begin{equation} \phi_{\hat{\mathbf{k}}}(\mathbf{s}) = \left| \frac{\partial\mathbf{s}}{\partial\mathbf{r}} \right|^{1/2} \exp\left[ 2\pi\imath\hat{\mathbf{k}}\cdot\mathbf{s} \right] = \dfrac{1}{\sqrt{\mathrm{det}(\mathbf{h})}} \exp\left[ 2\pi\imath\hat{\mathbf{k}}\cdot\mathbf{s} \right] \end{equation} $$

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:

$$ \begin{equation} \int \, \mathrm{d}\mathbf{s} \, \left|\frac{\partial\mathbf{r}}{\partial\mathbf{s}}\right| \end{equation} \phi^{*}_{\hat{\mathbf{k}}}(\mathbf{s}) \phi_{\hat{\mathbf{k}}'}(\mathbf{s}) = \int \, \mathrm{d}\mathbf{s} \, \exp[ -2\pi\imath\hat{\mathbf{k}}\cdot\mathbf{s} ] \, \exp[ 2\pi\imath\hat{\mathbf{k}}'\cdot\mathbf{s} ] = \delta_{\hat{\mathbf{k}},\hat{\mathbf{k}}'} $$

We can define a FBR from the curvilinear plane-wave functions $\{\phi_{\hat{\mathbf{k}}}\}$, restricting:

$$ \begin{equation} \hat{\mathbf{k}}_{\alpha} \in \left[-\mathcal{N},\mathcal{N}\right] \end{equation} $$

where the size of DVR quadrature grid is $2\mathcal{N}+1$

The corresponding DVR functions $\mu_{lmn}(\mathbf{s})$ can be defined as

$$ \begin{equation} \mu_{lmn}(\mathbf{s}) = \sum_{\hat{\mathbf{k}}} \phi_{\hat{\mathbf{k}}}(\mathbf{s}) \left\langle \left.\phi_{\hat{\mathbf{k}}}\right| \mu_{lmn} \right\rangle \end{equation} $$

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

$$ \begin{equation} \mathbf{s}_{lmn} = \left( \frac{l}{2\mathcal{N}+1}, \frac{m}{2\mathcal{N}+1}, \frac{n}{2\mathcal{N}+1} \right) \end{equation} $$

The DVR grid points in Cartesian coordinates are given by the matrix-vector product:

$$ \begin{equation} \mathbf{r}_{lmn} = \mathbf{h}\mathbf{s}_{lmn} \end{equation} $$

The corresponding DVR functions $\mu_{lmn}(\mathbf{s})$ then become

$$ \begin{equation} \mu_{lmn}(\mathbf{s}) = \sum_{\hat{\mathbf{k}}} \frac{1}{\sqrt{\left(2\mathcal{N}+1\right)^3}} \left| \frac{\partial\mathbf{s}}{\partial\mathbf{r}} \right|^{1/2} \cos\left[ 2\pi\hat{\mathbf{k}}\cdot \left( \mathbf{s} - \mathbf{s}_{lmn} \right) \right] \end{equation} $$

$$ \begin{equation} \left| \frac{\partial\mathbf{s}}{\partial\mathbf{r}} \right|^{1/2} = \dfrac{1}{\sqrt{\mathrm{det}(\mathbf{h})}} \end{equation} $$



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]:





eval_μ_lmn (generic function with 2 methods)



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) )









    



8.003188339268475e-16 + 0.0im
2.202498581157318 + 0.0im

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}$

$$ \begin{equation} \int\,\mathrm{d}\mathbf{s}\, \left| \frac{\partial\mathbf{r}}{\partial\mathbf{s}} \right| \mu_{lmn}(\mathbf{s}) \mu_{l'm'n'}(\mathbf{s}) = \delta_{ll'}\delta_{mm'}\delta_{nn'} \end{equation} $$

$$ \begin{equation} \mu_{lmn}\left(\mathbf{s}_{l'm'n'}\right) = \delta_{ll'}\delta_{mm'}\delta_{nn'} \sqrt{\frac{\left(2\mathcal{N}+1\right)^3} {\left|\partial\mathbf{r}/\partial\mathbf{s}\right|}} \end{equation} $$

Kinetic energy term

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 [ ]: