File: genkpakq.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (156 lines) | stat: -rw-r--r-- 4,079 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156

! Copyright (C) 2017 T. Mueller, J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine genkpakq
use modmain
use modulr
implicit none
! local variables
integer i1,i2,i3,j1,j2,j3
integer ikpa,ik0,ik,iq,ir
real(8) v(3)
! allocatable arrays
real(8), allocatable :: vkl0(:,:),vkc0(:,:)
! make sure the kappa-point grid is odd
do i1=1,3
  if (mod(ngridkpa(i1),2).eq.0) ngridkpa(i1)=ngridkpa(i1)+1
end do
! number of kappa-points
nkpa=ngridkpa(1)*ngridkpa(2)*ngridkpa(3)
! integer grid intervals for the kappa-points
intkpa(1,:)=ngridkpa(:)/2-ngridkpa(:)+1
intkpa(2,:)=ngridkpa(:)/2
! Q-point grid should be twice the kappa-point grid
ngridq(:)=2*ngridkpa(:)-1
! find next largest FFT-compatible grid size
call nfftifc(ngridq(1))
call nfftifc(ngridq(2))
call nfftifc(ngridq(3))
! total number of Q-points
nqpt=ngridq(1)*ngridq(2)*ngridq(3)
! integer grid intervals for the Q-points
intq(1,:)=ngridq(:)/2-ngridq(:)+1
intq(2,:)=ngridq(:)/2
! allocate global Q-point arrays
if (allocated(ivq)) deallocate(ivq)
allocate(ivq(3,nqpt))
if (allocated(ivqiq)) deallocate(ivqiq)
allocate(ivqiq(intq(1,1):intq(2,1),intq(1,2):intq(2,2),intq(1,3):intq(2,3)))
if (allocated(iqfft)) deallocate(iqfft)
allocate(iqfft(nqpt))
if (allocated(vql)) deallocate(vql)
allocate(vql(3,nqpt))
if (allocated(vqc)) deallocate(vqc)
allocate(vqc(3,nqpt))
if (allocated(qc)) deallocate(qc)
allocate(qc(nqpt))
! store the kappa-points as the first nkpa entries in the Q-point arrays
iq=0
do i1=intkpa(1,1),intkpa(2,1)
  do i2=intkpa(1,2),intkpa(2,2)
    do i3=intkpa(1,3),intkpa(2,3)
      iq=iq+1
      ivq(1,iq)=i1
      ivq(2,iq)=i2
      ivq(3,iq)=i3
    end do
  end do
end do
! store the remaining Q-points
do i1=intq(1,1),intq(2,1)
  do i2=intq(1,2),intq(2,2)
    do i3=intq(1,3),intq(2,3)
      if ((i1.lt.intkpa(1,1)).or.(i1.gt.intkpa(2,1)).or. &
          (i2.lt.intkpa(1,2)).or.(i2.gt.intkpa(2,2)).or. &
          (i3.lt.intkpa(1,3)).or.(i3.gt.intkpa(2,3))) then
        iq=iq+1
        ivq(1,iq)=i1
        ivq(2,iq)=i2
        ivq(3,iq)=i3
      end if
    end do
  end do
end do
! ensure the first point is the zero vector
do iq=1,nkpa
  if ((ivq(1,iq).eq.0).and.(ivq(2,iq).eq.0).and.(ivq(3,iq).eq.0)) then
    ivq(:,iq)=ivq(:,1)
    ivq(:,1)=0
    exit
  end if
end do
do iq=1,nqpt
  i1=ivq(1,iq)
  i2=ivq(2,iq)
  i3=ivq(3,iq)
! map from (i1,i2,i3) to Q-vector index
  ivqiq(i1,i2,i3)=iq
! Fourier transform index
  if (i1.ge.0) then
    j1=i1
  else
    j1=ngridq(1)+i1
  end if
  if (i2.ge.0) then
    j2=i2
  else
    j2=ngridq(2)+i2
  end if
  if (i3.ge.0) then
    j3=i3
  else
    j3=ngridq(3)+i3
  end if
  iqfft(iq)=j3*ngridq(2)*ngridq(1)+j2*ngridq(1)+j1+1
! Q-vector in Cartesian coordinates
  vqc(:,iq)=dble(i1)*bvecu(:,1) &
           +dble(i2)*bvecu(:,2) &
           +dble(i3)*bvecu(:,3)
! length of Q-vector
  qc(iq)=sqrt(vqc(1,iq)**2+vqc(2,iq)**2+vqc(3,iq)**2)
! Q-vector in (unit cell) lattice coordinates
  call r3mv(binv,vqc(:,iq),vql(:,iq))
end do
! store the R-vectors in Cartesian coordinates spanning the ultracell
if (allocated(vrcu)) deallocate(vrcu)
allocate(vrcu(3,nqpt))
ir=0
do i3=0,ngridq(3)-1
  v(3)=dble(i3)/dble(ngridq(3))
  do i2=0,ngridq(2)-1
    v(2)=dble(i2)/dble(ngridq(2))
    do i1=0,ngridq(1)-1
      v(1)=dble(i1)/dble(ngridq(1))
      ir=ir+1
      call r3mv(avecu,v,vrcu(:,ir))
    end do
  end do
end do
! store the existing k-point and weight arrays
allocate(vkl0(3,nkpt),vkc0(3,nkpt))
if (allocated(wkpt0)) deallocate(wkpt0)
allocate(wkpt0(nkpt))
vkl0(:,1:nkpt)=vkl(:,1:nkpt)
vkc0(:,1:nkpt)=vkc(:,1:nkpt)
wkpt0(1:nkpt)=wkpt(1:nkpt)
! number of k+kappa-points
nkpt0=nkpt
nkpt=nkpt0*nkpa
! deallocate and reallocate k-point and weight arrays
deallocate(vkl,vkc,wkpt)
allocate(vkl(3,nkpt),vkc(3,nkpt),wkpt(nkpt))
ik=0
do ik0=1,nkpt0
  do ikpa=1,nkpa
    ik=ik+1
    vkl(:,ik)=vkl0(:,ik0)+vql(:,ikpa)
    vkc(:,ik)=vkc0(:,ik0)+vqc(:,ikpa)
    wkpt(ik)=wkpt0(ik0)/dble(nkpa)
  end do
end do
deallocate(vkl0,vkc0)
return
end subroutine