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