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