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