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_mmnk_symm
8 : use m_juDFT
9 : USE m_types
10 : private
11 : public:: wann_mmnk_symm
12 : contains
13 : c******************************************************************
14 : c Find out minimal set of k-point-pairs that have to be
15 : c calculated; map symmetry-related k-point-pairs to this
16 : c minimal set.
17 : c Frank Freimuth
18 : c******************************************************************
19 2 : subroutine wann_mmnk_symm(input,kpts,
20 2 : > fullnkpts,nntot,bpt,gb,l_bzsym,
21 2 : > irreduc,mapkoper,l_p0,film,nop,
22 2 : > invtab,mrot,tau,
23 2 : < pair_to_do,maptopair,kdiff,l_q,param_file)
24 :
25 : USE m_constants
26 :
27 : implicit none
28 :
29 : TYPE(t_input),INTENT(IN) :: input
30 : TYPE(t_kpts),INTENT(IN) :: kpts
31 : integer,intent(in) :: nop
32 : integer,intent(in) :: mrot(3,3,nop)
33 : integer,intent(in) :: invtab(nop)
34 : integer,intent(in) :: fullnkpts
35 : integer,intent(in) :: nntot
36 : integer,intent(in) :: bpt(nntot,fullnkpts)
37 : integer,intent(in) :: gb(3,nntot,fullnkpts)
38 : logical,intent(in) :: l_bzsym
39 : integer,intent(in) :: irreduc(fullnkpts)
40 : integer,intent(in) :: mapkoper(fullnkpts)
41 : logical,intent(in) :: l_p0,film
42 : real,intent(in) :: tau(3,nop)
43 :
44 : integer,intent(out):: pair_to_do(fullnkpts,nntot)
45 : integer,intent(out):: maptopair(3,fullnkpts,nntot)
46 : real,intent(out) :: kdiff(3,nntot)
47 :
48 : logical, intent(in) :: l_q
49 : character(len=20), intent(in) :: param_file
50 :
51 : integer :: ikpt,ikpt_k,kptibz,ikpt_b
52 : integer :: num_pair,kptibz_b
53 2 : integer :: index(fullnkpts,nntot)
54 : integer :: oper,oper_b,num_rot,num_conj,repo,repo_b
55 : integer :: k,kx,ky,repkpt,repkpt_b,repkpt_bb
56 : integer :: sign,sign_b,ngis,ngis_b
57 : logical :: startloop,l_file
58 2 : real :: kpoints(3,fullnkpts)
59 : real :: kdiffvec(3)
60 2 : integer :: multtab(nop,nop)
61 : integer :: fullnkpts_tmp,kr
62 : real :: scale
63 : real :: brot(3)
64 : logical :: l_testnosymm,l_nosymm1,l_nosymm2
65 :
66 2 : call timestart("wann_mmnk_symm")
67 146 : index(:,:)=0
68 2 : num_pair=0
69 2 : num_rot=0
70 2 : num_conj=0
71 146 : pair_to_do(:,:)=0
72 :
73 2 : inquire(file='testnosymm',exist=l_testnosymm)
74 :
75 : c-----Test for nonsymmorphic space groups
76 2 : if(l_bzsym)then
77 0 : do ikpt=1,fullnkpts
78 0 : oper=abs(mapkoper(ikpt))
79 0 : if( any( abs(tau(:,oper)).gt.1.e-6 ) ) l_testnosymm=.true.
80 : enddo
81 : endif
82 :
83 2 : if(l_bzsym)then
84 0 : call close_pt(nop,mrot,multtab)
85 : endif
86 :
87 18 : do 10 ikpt = 1,fullnkpts ! loop by k-points starts
88 16 : l_nosymm1=.false.
89 16 : kptibz=ikpt
90 16 : if(l_bzsym) then
91 0 : kptibz=irreduc(ikpt)
92 0 : oper=mapkoper(ikpt)
93 0 : if(oper.lt.0)then
94 0 : oper=-oper
95 0 : sign=-1
96 : else
97 : sign=1
98 : endif
99 0 : if(
100 : & any( abs(tau(:,oper)).gt.1.e-6 )
101 0 : & )l_nosymm1=.true.
102 : endif
103 :
104 144 : do 15 ikpt_b = 1,nntot
105 128 : l_nosymm2=.false.
106 128 : if(index(ikpt,ikpt_b).eq.1)cycle
107 128 : kptibz_b=bpt(ikpt_b,ikpt)
108 128 : if(l_bzsym) then
109 0 : oper_b=mapkoper(kptibz_b)
110 0 : kptibz_b=irreduc(kptibz_b)
111 0 : if(oper_b.lt.0)then
112 0 : oper_b=-oper_b
113 0 : sign_b=-1
114 : else
115 : sign_b=1
116 : endif
117 0 : if(
118 : & any( abs(tau(:,oper_b)).gt.1.e-6 )
119 0 : & )l_nosymm2=.true.
120 : endif
121 :
122 128 : if(l_testnosymm)goto 33
123 128 : if(l_nosymm1.or.l_nosymm2)goto 33
124 :
125 : c***************************************************************
126 : c..the conjugation selection rule, which speeds up significantly
127 : c***************************************************************
128 896 : do ikpt_k = 1,nntot
129 : if((bpt(ikpt_k,bpt(ikpt_b,ikpt)).eq.ikpt).and.
130 : & gb(1,ikpt_b,ikpt).eq.(-gb(1,ikpt_k,bpt(ikpt_b,ikpt))).and.
131 : & gb(2,ikpt_b,ikpt).eq.(-gb(2,ikpt_k,bpt(ikpt_b,ikpt))).and.
132 832 : & gb(3,ikpt_b,ikpt).eq.(-gb(3,ikpt_k,bpt(ikpt_b,ikpt))).and.
133 64 : & (index(bpt(ikpt_b,ikpt),ikpt_k).eq.1))then
134 64 : index(ikpt,ikpt_b)=1
135 64 : maptopair(1,ikpt,ikpt_b)=bpt(ikpt_b,ikpt)
136 64 : maptopair(2,ikpt,ikpt_b)=ikpt_k
137 64 : maptopair(3,ikpt,ikpt_b)=1
138 : c print*,"conjugation"
139 64 : num_conj=num_conj+1
140 64 : goto 15
141 : endif
142 : enddo !ikpt_k
143 : c****************************************************************
144 : c check whether k-point pairs can be mapped onto each other
145 : c by rotation
146 : c****************************************************************
147 64 : if(l_bzsym)then
148 : c if(all(gb(:,ikpt_b,ikpt).eq.0))then
149 0 : do k=1,fullnkpts
150 0 : if(irreduc(k).eq.kptibz)then
151 0 : repkpt=k
152 0 : repo=mapkoper(k)
153 0 : if(repo.lt.0)then
154 0 : repo=-repo
155 0 : ngis=-1
156 : else
157 : ngis=1
158 : endif
159 :
160 0 : do kx=1,fullnkpts
161 0 : if(irreduc(kx).eq.kptibz_b)then
162 0 : repkpt_bb=kx
163 0 : repo_b=mapkoper(kx)
164 0 : if(repo_b.lt.0)then
165 0 : repo_b=-repo_b
166 0 : ngis_b=-1
167 : else
168 : ngis_b=1
169 : endif
170 0 : do ky=1,nntot
171 0 : if(bpt(ky,repkpt).eq.repkpt_bb)then
172 0 : repkpt_b=ky
173 0 : if (index(repkpt,repkpt_b).eq.1)then
174 0 : if(.not.all(gb(:,ikpt_b,ikpt).eq.0))then
175 0 : brot(:)=0.0
176 0 : do kr=1,3
177 : brot(:)=brot(:)
178 : & -sign_b*mrot(kr,:,oper)*gb(kr,ikpt_b,ikpt)
179 0 : & +ngis_b*mrot(kr,:,oper_b)*gb(kr,repkpt_b,repkpt)
180 : enddo
181 0 : if( any( abs(brot).gt.1e-6 ) )cycle
182 : endif
183 0 : if(sign*ngis*multtab(invtab(oper),repo).eq.
184 : & sign_b*ngis_b*multtab(invtab(oper_b),repo_b))then
185 0 : maptopair(1,ikpt,ikpt_b)=repkpt
186 0 : maptopair(2,ikpt,ikpt_b)=repkpt_b
187 0 : maptopair(3,ikpt,ikpt_b)=2+(1-ngis*sign)/2
188 0 : index(ikpt,ikpt_b)=1
189 0 : num_rot=num_rot+1
190 0 : goto 15
191 : endif
192 : endif
193 : endif
194 : enddo
195 : endif
196 : enddo
197 : endif
198 : enddo
199 : c endif !gb=0
200 : endif
201 :
202 : 33 continue
203 :
204 64 : index(ikpt,ikpt_b)=1
205 64 : num_pair=num_pair+1
206 64 : pair_to_do(ikpt,ikpt_b)=num_pair
207 :
208 16 : 15 continue !loop over nearest neighbor k-points
209 2 : 10 continue ! end of cycle by the k-points
210 :
211 2 : if(l_p0)then
212 1 : write(oUnit,*)"pairs to calculate: ",num_pair
213 1 : write(oUnit,*)"maps by conjugation: ",num_conj
214 1 : write(oUnit,*)"maps by rotation:", num_rot
215 1 : write(oUnit,*)"num_pair+num_rot+num_conj:",
216 2 : + num_pair+num_conj+num_rot
217 1 : if(.not.l_q) then
218 1 : write(oUnit,*)"fullnkpts*nntot:", fullnkpts*nntot
219 : else
220 0 : write(oUnit,*)"fullnqpts*nntot:", fullnkpts*nntot
221 : endif
222 : endif !l_p0
223 :
224 : c*****************************************************************
225 : c determine difference vectors that occur on the k-mesh
226 : c*****************************************************************
227 2 : if (l_bzsym) then
228 0 : l_file=.false.
229 0 : IF(.NOT.l_q)THEN
230 0 : inquire(file='w90kpts',exist=l_file)
231 0 : IF(.NOT.l_file) CALL juDFT_error
232 : + ("w90kpts not found, needed if bzsym",calledby
233 0 : + ="wann_mmnk_symm")
234 0 : open(412,file='w90kpts',form='formatted')
235 : ELSE
236 0 : inquire(file='w90qpts',exist=l_file)
237 0 : IF(.NOT.l_file) CALL juDFT_error
238 : + ("w90qpts not found, needed if bzsym",calledby
239 0 : + ="wann_mmnk_symm")
240 0 : open(412,file='w90qpts',form='formatted')
241 : ENDIF
242 0 : read(412,*)fullnkpts_tmp,scale
243 0 : do k=1,fullnkpts
244 0 : read(412,*)kpoints(:,k)
245 : enddo
246 0 : kpoints=kpoints/scale
247 0 : close(412)
248 : else
249 2 : IF(.not.l_q)THEN
250 2 : fullnkpts_tmp = kpts%nkpt
251 18 : do k=1,fullnkpts
252 66 : kpoints(:,k) = kpts%bk(:,k)
253 : enddo
254 : ELSE
255 0 : open(412,file=param_file,form='formatted')
256 0 : read(412,*)fullnkpts_tmp,scale
257 0 : do k=1,fullnkpts
258 0 : read(412,*)kpoints(:,k)
259 : enddo
260 0 : kpoints(:,:)=kpoints/scale
261 : ENDIF
262 :
263 2 : if (film) kpoints(3,:)=0.0
264 2 : close(412)
265 : endif
266 :
267 2 : if(l_p0)then
268 1 : IF(.not.l_q) THEN
269 1 : print*,"vectors combining nearest neighbor k-points:"
270 : ELSE
271 0 : print*,"vectors combining nearest neighbor q-points:"
272 : ENDIF
273 : endif
274 2 : ky=1
275 18 : do k=1,fullnkpts
276 146 : do kx=1,nntot
277 512 : kdiffvec=kpoints(:,bpt(kx,k))+gb(:,kx,k)-kpoints(:,k)
278 576 : do ikpt=1,ky-1
279 1056 : if(all(abs(kdiff(:,ikpt)-kdiffvec).le.0.0001))goto 200
280 : enddo
281 16 : IF(ky>nntot) CALL juDFT_error("problem in wann_mmnk_symm"
282 0 : + ,calledby ="wann_mmnk_symm")
283 64 : kdiff(:,ky)=kdiffvec(:)
284 16 : if(l_p0)then
285 8 : print*,ky,k,kx,kdiff(:,ky)
286 : endif
287 :
288 128 : ky=ky+1
289 16 : 200 continue
290 :
291 : enddo
292 : enddo
293 :
294 2 : call timestop("wann_mmnk_symm")
295 2 : end subroutine
296 :
297 0 : SUBROUTINE close_pt(
298 0 : > nops,mrot,
299 0 : < mtable)
300 :
301 : USE m_constants
302 :
303 : IMPLICIT NONE
304 :
305 : INTEGER, INTENT (IN) :: nops,mrot(3,3,nops)
306 : INTEGER, INTENT (OUT) :: mtable(nops,nops) ! table(i,j) = {R_i|0}{R_j|0}
307 :
308 0 : INTEGER :: i,j,k,mp(3,3),map(nops)
309 :
310 0 : call timestart("close_pt")
311 : !---> loop over all operations
312 0 : DO j=1,nops
313 :
314 0 : map(1:nops) = 0
315 :
316 : !---> multiply {R_j|0}{R_i|0}
317 0 : DO i=1,nops
318 0 : mp = matmul( mrot(:,:,j) , mrot(:,:,i) )
319 :
320 : !---> determine which operation this is
321 0 : DO k = 1, nops
322 0 : IF ( all( mp(:,:)==mrot(:,:,k) ) ) THEN
323 0 : IF ( map(i) .eq. 0 ) THEN
324 0 : map(i) = k
325 : ELSE
326 0 : WRITE (oUnit,'(" Symmetry error : multiple ops")')
327 : CALL juDFT_error("close_pt: Multiple ops (Bravais)"
328 0 : + ,calledby ="wann_mmnk_symm")
329 : ENDIF
330 : ENDIF
331 : ENDDO
332 :
333 0 : IF (map(i).eq.0) THEN
334 0 : WRITE (oUnit,'(" Group not closed (Bravais lattice)")')
335 : WRITE (oUnit,'(" operation j=",i2," map=",12i4,:/,
336 0 : & (21x,12i4))') j, map(1:nops)
337 : CALL juDFT_error("close_pt: Not closed",calledby
338 0 : + ="wann_mmnk_symm")
339 : ENDIF
340 : ENDDo
341 0 : mtable(j,1:nops) = map(1:nops)
342 : ENDDO
343 :
344 0 : call timestop("close_pt")
345 0 : END SUBROUTINE close_pt
346 :
347 : end module m_wann_mmnk_symm
|
