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
157
158
159
160
161
162
163
164
165
166
|
! Copyright (C) 2002-2007 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU Lesser General Public
! License. See the file COPYING for license details.
!BOP
! !ROUTINE: genppts
! !INTERFACE:
subroutine genppts(tfbz,nsym,sym,ngridp,npptnr,epslat,bvec,boxl,nppt,ipvip, &
ipvipnr,ivp,vpl,vpc,wppt,wpptnr)
! !INPUT/OUTPUT PARAMETERS:
! tfbz : .true. if vpl and vpc should be mapped to the first Brillouin zone
! (in,logical)
! nsym : number of point group symmetries used for reduction, set to 1 for
! no reduction (in,integer)
! sym : symmetry matrices in lattice coordinates (in,integer(3,3,*))
! ngridp : p-point grid sizes (in,integer(3))
! npptnr : number of non-reduced p-points: ngridp(1)*ngridp(2)*ngridp(3)
! (in,integer)
! epslat : tolerance for determining identical vectors (in,real)
! bvec : reciprocal lattice vectors (in,real(3,3))
! boxl : corners of box containing p-points in lattice coordinates, the
! zeroth vector is the origin (in,real(3,0:3))
! nppt : total number of p-points (out,integer)
! ipvip : map from (i1,i2,i3) to reduced p-point index
! (out,integer(0:ngridp(1)-1,0:ngridp(2)-1,0:ngridp(3)-1))
! ipvipnr : map from (i1,i2,i3) to non-reduced p-point index
! (out,integer(0:ngridp(1)-1,0:ngridp(2)-1,0:ngridp(3)-1))
! ivp : integer coordinates of the p-points
! (out,integer(3,ngridp(1)*ngridp(2)*ngridp(3)))
! vpl : lattice coordinates of each p-point
! (out,real(3,ngridp(1)*ngridp(2)*ngridp(3)))
! vpc : Cartesian coordinates of each p-point
! (out,real(3,ngridp(1)*ngridp(2)*ngridp(3)))
! wppt : weights of each reduced p-point
! (out,real(ngridp(1)*ngridp(2)*ngridp(3)))
! wpptnr : weight of each non-reduced p-point (out,real)
! !DESCRIPTION:
! This routine is used for generating $k$-point or $q$-point sets. Since these
! are stored in global arrays, the points passed to this and other routines
! are referred to as $p$-points. In lattice coordinates, the ${\bf p}$ vectors
! are given by
! $$ {\bf p}=\left(\begin{matrix} & & \\
! {\bf B}_2-{\bf B}_1 & {\bf B}_3-{\bf B}_1 & {\bf B}_4-{\bf B}_1 \\
! & & \end{matrix}\right)
! \left(\begin{matrix}i_1/n_1 \\ i_2/n_2 \\ i_3/n_3 \end{matrix}\right)
! +{\bf B}_1 $$
! where $i_j$ runs from 0 to $n_j-1$, and the ${\bf B}$ vectors define the
! corners of a box with ${\bf B}_1$ as the origin. If {\tt tfbz} is
! {\tt .true.} then each {\tt vpl} vector is mapped to the first Brillouin
! zone. If {\tt tfbz} is {\tt .false.} and then the coordinates of each
! {\tt vpl} are mapped to the $[0,1)$ interval. The $p$-point weights are
! stored in {\tt wppt} and the array {\tt ipvip} contains the map from the
! integer coordinates to the reduced index.
!
! !REVISION HISTORY:
! Created August 2002 (JKD)
! Updated April 2007 (JKD)
! Added mapping to the first Brillouin zone, September 2008 (JKD)
! Made independent of modmain, February 2010 (JKD)
!EOP
!BOC
implicit none
! arguments
logical, intent(in) :: tfbz
integer, intent(in) :: nsym,sym(3,3,*),ngridp(3),npptnr
real(8), intent(in) :: epslat,bvec(3,3),boxl(3,0:3)
integer, intent(out) :: nppt
integer, intent(out) :: ipvip(0:ngridp(1)-1,0:ngridp(2)-1,0:ngridp(3)-1)
integer, intent(out) :: ipvipnr(0:ngridp(1)-1,0:ngridp(2)-1,0:ngridp(3)-1)
integer, intent(out) :: ivp(3,npptnr)
real(8), intent(out) :: vpl(3,npptnr),vpc(3,npptnr)
real(8), intent(out) :: wppt(npptnr),wpptnr
! local variables
integer i1,i2,i3,i
integer isym,ip,jp
real(8) v1(3),v2(3),v3(3)
real(8) b(3,3),t1
if ((ngridp(1) < 1).or.(ngridp(2) < 1).or.(ngridp(3) < 1)) then
write(*,*)
write(*,'("Error(genppts): invalid ngridp :",3(X,I0))') ngridp
write(*,*)
stop
end if
if (npptnr /= ngridp(1)*ngridp(2)*ngridp(3)) then
write(*,*)
write(*,'("Error(genppts): mismatched npptnr and ngridp :",4(X,I0))') npptnr,&
ngridp
write(*,*)
stop
end if
! box vector matrix
b(1:3,1)=boxl(1:3,1)-boxl(1:3,0)
b(1:3,2)=boxl(1:3,2)-boxl(1:3,0)
b(1:3,3)=boxl(1:3,3)-boxl(1:3,0)
! weight of each non-reduced p-point
wpptnr=1.d0/dble(ngridp(1)*ngridp(2)*ngridp(3))
ip=0
jp=npptnr
do i3=0,ngridp(3)-1
v1(3)=dble(i3)/dble(ngridp(3))
do i2=0,ngridp(2)-1
v1(2)=dble(i2)/dble(ngridp(2))
do i1=0,ngridp(1)-1
v1(1)=dble(i1)/dble(ngridp(1))
call r3mv(b,v1,v2)
v2(1:3)=v2(1:3)+boxl(1:3,0)
! map vector components to [0,1)
call r3frac(epslat,v2)
if (nsym > 1) then
! determine if this point is equivalent to one already in the set
do isym=1,nsym
call i3mtrv(sym(:,:,isym),v2,v3)
call r3frac(epslat,v3)
do i=1,ip
t1=abs(vpl(1,i)-v3(1))+abs(vpl(2,i)-v3(2))+abs(vpl(3,i)-v3(3))
if (t1 < epslat) then
! equivalent p-point found so add to existing weight
ipvip(i1,i2,i3)=i
wppt(i)=wppt(i)+wpptnr
! add new point to back of set
ipvipnr(i1,i2,i3)=jp
ivp(1,jp)=i1; ivp(2,jp)=i2; ivp(3,jp)=i3
vpl(1:3,jp)=v2(1:3)
wppt(jp)=0.d0
jp=jp-1
goto 10
end if
end do
end do
end if
! add new point to front of set
ip=ip+1
ipvip(i1,i2,i3)=ip
ipvipnr(i1,i2,i3)=ip
ivp(1,ip)=i1; ivp(2,ip)=i2; ivp(3,ip)=i3
vpl(1:3,ip)=v2(1:3)
wppt(ip)=wpptnr
10 continue
end do
end do
end do
nppt=ip
do ip=1,npptnr
! map vpl to the first Brillouin zone if required
if (tfbz) call vecfbz(epslat,bvec,vpl(:,ip))
! determine the Cartesian coordinates of the p-points
call r3mv(bvec,vpl(:,ip),vpc(:,ip))
end do
contains
pure subroutine i3mtrv(a,x,y)
implicit none
! arguments
integer, intent(in) :: a(3,3)
real(8), intent(in) :: x(3)
real(8), intent(out) :: y(3)
y(1)=a(1,1)*x(1)+a(2,1)*x(2)+a(3,1)*x(3)
y(2)=a(1,2)*x(1)+a(2,2)*x(2)+a(3,2)*x(3)
y(3)=a(1,3)*x(1)+a(2,3)*x(2)+a(3,3)*x(3)
end subroutine
end subroutine
!EOC
|
