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