Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : module m_wann_kptsreduc2
8 : use m_juDFT
9 : contains
10 0 : subroutine wann_kptsreduc2(
11 : > mhp,
12 0 : > nop,mrot,bmat,tau,film,
13 : > l_nocosoc)
14 : c*****************************************************************
15 : c Apply the symmetries to reduce the number of k-points.
16 : c Frank Freimuth
17 : c*****************************************************************
18 :
19 : USE m_constants
20 :
21 : implicit none
22 :
23 : integer, intent(in) :: mhp(3)
24 : integer, intent(in) :: nop
25 : logical,intent(in) :: film
26 : real, intent(in) :: bmat(3,3)
27 : real,intent(in) :: tau(3,nop)
28 : integer, intent(in) :: mrot(3,3,nop)
29 : logical,intent(in) :: l_nocosoc
30 :
31 0 : real,allocatable :: weight(:)
32 : integer :: reduznumk,ikpt2
33 0 : integer,allocatable :: mapk(:),mapkoper(:)
34 : integer :: ikpt,nmat,k,j
35 : integer :: k1,k2,k3
36 : real :: bbkpt(3)
37 : real :: wk
38 : integer :: brot(3)
39 : integer :: bkpt(3),bkpt2(3)
40 : real :: kpoint(3)
41 : integer :: oper,iter
42 : logical :: l_file
43 : real :: scale
44 : integer :: nkpts
45 : integer :: sumweights
46 0 : integer,allocatable :: irreduc(:)
47 0 : integer,allocatable :: shiftkpt(:,:)
48 : real :: bbmat(3,3)
49 0 : real,allocatable :: kptslen(:)
50 : integer :: i,kk1,kk2,kk3
51 :
52 0 : call timestart("wann_kptsreduc2")
53 0 : bbmat=matmul(bmat,transpose(bmat))
54 0 : write(oUnit,*) "Apply the symmetries to w90kpts"
55 : c**************************************************************
56 : c The array 'mhp' specifies the Monkhorst-Pack mesh.
57 : c**************************************************************
58 0 : nkpts = mhp(1) * mhp(2) * mhp(3)
59 :
60 0 : allocate(mapk(nkpts),mapkoper(nkpts))
61 0 : allocate(weight(nkpts),irreduc(nkpts))
62 :
63 : c*********************************************************
64 : c determine lengths of kpoints
65 : c*********************************************************
66 0 : allocate( kptslen(nkpts) )
67 0 : iter=0
68 0 : do k3=0,mhp(3)-1
69 0 : do k2=0,mhp(2)-1
70 0 : do k1=0,mhp(1)-1
71 0 : iter=iter+1
72 0 : kpoint(1)=real(k1)/real(mhp(1))
73 0 : kpoint(2)=real(k2)/real(mhp(2))
74 0 : kpoint(3)=real(k3)/real(mhp(3))
75 : kptslen(iter)=dot_product( axmintx(kpoint(:)),
76 0 : & matmul(bbmat,axmintx(kpoint(:))) )
77 : enddo !k1
78 : enddo !k2
79 : enddo !k3
80 :
81 : c***********************************************
82 : c determine irreducible part
83 : c***********************************************
84 0 : weight(:)=1.0
85 0 : mapk(:)=0
86 0 : reduznumk=0
87 0 : mapkoper(:)=0
88 0 : irreduc(:)=0
89 0 : shiftkpt(:,:)=0
90 0 : ikpt=0
91 0 : do k3=0,mhp(3)-1
92 0 : do k2=0,mhp(2)-1
93 0 : do k1=0,mhp(1)-1
94 0 : ikpt=ikpt+1
95 0 : if(mapk(ikpt).ne.0) cycle
96 0 : reduznumk=reduznumk+1
97 0 : irreduc(ikpt)=reduznumk
98 : c bkpt(:)=kpoints(:,ikpt)
99 0 : bkpt(1)=k1*mhp(2)*mhp(3)
100 0 : bkpt(2)=k2*mhp(1)*mhp(3)
101 0 : bkpt(3)=k3*mhp(1)*mhp(2)
102 0 : do oper=1,nop
103 : do i=1,3
104 : brot(i)=0.0
105 0 : do k=1,3
106 : brot(i)=brot(i)+mrot(k,i,oper)*bkpt(k)
107 : enddo
108 : enddo
109 : ikpt2=0
110 0 : do kk3=0,mhp(3)-1
111 0 : do kk2=0,mhp(2)-1
112 0 : do kk1=0,mhp(1)-1
113 0 : ikpt2=ikpt2+1
114 0 : if(ikpt2.le.ikpt)cycle
115 0 : bkpt2(1)=kk1*mhp(2)*mhp(3)
116 0 : bkpt2(2)=kk2*mhp(1)*mhp(3)
117 0 : bkpt2(3)=kk3*mhp(1)*mhp(2)
118 0 : if(mapk(ikpt2).ne.0)cycle
119 0 : if(abs(kptslen(ikpt2)-kptslen(ikpt)).gt.1.e-8)cycle
120 0 : kpoint(1)=real(bkpt(1)-bkpt2(1))/real(mhp(1)*mhp(2)*mhp(3))
121 0 : kpoint(2)=real(bkpt(2)-bkpt2(2))/real(mhp(1)*mhp(2)*mhp(3))
122 0 : kpoint(3)=real(bkpt(3)-bkpt2(3))/real(mhp(1)*mhp(2)*mhp(3))
123 0 : if( all( axmintx(kpoint(:)) .lt.1e-6 ) ) then
124 0 : weight(ikpt)=weight(ikpt)+1.0
125 0 : mapk(ikpt2)=ikpt
126 0 : mapkoper(ikpt2)=oper
127 : endif
128 : enddo !kk1
129 : enddo !kk2
130 : enddo !kk3
131 : enddo !oper
132 0 : if(.not.l_nocosoc)then
133 0 : do oper=1,nop
134 : do i=1,3
135 : brot(i)=0.0
136 0 : do k=1,3
137 : brot(i)=brot(i)+mrot(k,i,oper)*bkpt(k)
138 : enddo
139 : enddo
140 : ikpt2=0
141 0 : do kk3=0,mhp(3)-1
142 0 : do kk2=0,mhp(2)-1
143 0 : do kk1=0,mhp(1)-1
144 0 : ikpt2=ikpt2+1
145 0 : if(ikpt2.le.ikpt)cycle
146 0 : bkpt2(1)=kk1*mhp(2)*mhp(3)
147 0 : bkpt2(2)=kk2*mhp(1)*mhp(3)
148 0 : bkpt2(3)=kk3*mhp(1)*mhp(2)
149 0 : if(mapk(ikpt2).ne.0)cycle
150 0 : if(abs(kptslen(ikpt2)-kptslen(ikpt)).gt.1.e-8)cycle
151 : kpoint(1)=
152 0 : & real(bkpt(1)-bkpt2(1))/real(mhp(1)*mhp(2)*mhp(3))
153 : kpoint(2)=
154 0 : & real(bkpt(2)-bkpt2(2))/real(mhp(1)*mhp(2)*mhp(3))
155 : kpoint(3)=
156 0 : & real(bkpt(3)-bkpt2(3))/real(mhp(1)*mhp(2)*mhp(3))
157 0 : if( all( axmintx(kpoint(:)).lt.1e-6 ) ) then
158 0 : weight(ikpt)=weight(ikpt)+1.0
159 0 : mapk(ikpt2)=ikpt
160 0 : mapkoper(ikpt2)=oper
161 : endif
162 : enddo !kk1
163 : enddo !kk2
164 : enddo !kk3
165 : enddo !oper
166 : endif !l_nocosoc
167 : enddo !k1
168 : enddo !k2
169 : enddo !k3
170 :
171 : c****************************************************
172 : c write results to files
173 : c w90kpts: whole Brillouin zone (for w90)
174 : c kpts: irreducible part (for fleur)
175 : c kptsmap: mapping from w90kpts to kpts
176 : c****************************************************
177 : c open(117,file='kptsmap',form='formatted',recl=1000)
178 : c do ikpt=1,nkpts
179 : c if(mapk(ikpt)==0)then
180 : c write(117,*)ikpt,irreduc(ikpt),1,0,0,0
181 : c else
182 : c write(117,*)ikpt,irreduc(mapk(ikpt)),mapkoper(ikpt),
183 : c & shiftkpt(:,ikpt)
184 : c endif
185 : c enddo
186 : c close(117)
187 :
188 0 : scale=1.000
189 :
190 0 : open(119,file='kpts',form='formatted')
191 0 : if (film)then
192 0 : write(119,'(i5,f20.10,3x,l1)')reduznumk,scale,.false.
193 : else
194 0 : write(119,'(i5,f20.10)')reduznumk,scale
195 : endif
196 0 : sumweights=0
197 :
198 0 : ikpt=0
199 0 : do k3=0,mhp(3)-1
200 0 : do k2=0,mhp(2)-1
201 0 : do k1=0,mhp(1)-1
202 0 : ikpt=ikpt+1
203 0 : kpoint(1)=real(k1)/real(mhp(1))
204 0 : kpoint(2)=real(k2)/real(mhp(2))
205 0 : kpoint(3)=real(k3)/real(mhp(3))
206 :
207 :
208 0 : if (mapk(ikpt)==0)then
209 0 : sumweights=sumweights+weight(ikpt)
210 0 : write(oUnit,*)"ikpt=",ikpt
211 0 : write(oUnit,*)"irreducible"
212 0 : write(oUnit,fmt='(a10,3f9.6)')"internal: ",kpoint(:)/scale
213 :
214 0 : if (film)then
215 0 : write(119,'(3f10.5)')kpoint(1:2),weight(ikpt)
216 : else
217 0 : write(119,'(4f10.5)')kpoint(:),weight(ikpt)
218 : endif
219 :
220 : elseif(mapkoper(ikpt).gt.0)then
221 : c write(oUnit,*)"ikpt=",ikpt
222 : c write(oUnit,*)"reducible"
223 : c write(oUnit,*)"map=",mapk(ikpt)
224 : c brot(:)=0.0
225 : c do k=1,3
226 : c brot(:)=brot(:)+
227 : c + mrot(k,:,mapkoper(ikpt))*kpoints(k,mapk(ikpt))
228 : c enddo
229 : c write(oUnit,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
230 : elseif(mapkoper(ikpt).lt.0)then
231 : c write(oUnit,*)"ikpt=",ikpt
232 : c write(oUnit,*)"reducible"
233 : c write(oUnit,*)"map=",mapk(ikpt)
234 : c brot(:)=0.0
235 : c do k=1,3
236 : c brot(:)=brot(:)-
237 : c + mrot(k,:,-mapkoper(ikpt))*kpoints(k,mapk(ikpt))
238 : c enddo
239 : c write(oUnit,'(a19,3f9.6)')"rotated internal: ",brot(:)/scale
240 : endif
241 : enddo !k1
242 : enddo !k2
243 : enddo !k3
244 0 : close(119)
245 :
246 0 : write(oUnit,*)"reduznumk=",reduznumk
247 0 : write(oUnit,*)"nkpts=",nkpts
248 0 : write(oUnit,*)"sumweights=",sumweights
249 :
250 0 : IF(sumweights/=nkpts) CALL juDFT_error
251 : + ("sum of weights differs from nkpts",calledby
252 0 : + ="wann_kptsreduc2")
253 0 : deallocate(mapk,mapkoper,weight,irreduc)
254 :
255 0 : call timestop("wann_kptsreduc2")
256 : contains
257 0 : real elemental function axmintx(x)
258 : implicit none
259 : real,intent(in) :: x
260 0 : axmintx = abs(x-nint(x))
261 : end function
262 :
263 : end subroutine wann_kptsreduc2
264 : end module m_wann_kptsreduc2
|
