TABLE OF CONTENTS
ABINIT/irrzg [ Functions ]
NAME
irrzg
FUNCTION
Find the irreducible zone in reciprocal space under the symmetry group with real space rotations in symrel(3,3,nsym). The (integer) rotation matrices symrel(3,3,nsym) express the new real space positions (e.g. rotated atom positions) in REDUCED coordinates, i.e. in coordinates expressed as fractions of real space primitive translations (atomic coordinates xred). tnons(3,nsym) express the associated nonsymmorphic translations, again in reduced coordinates. Special data structure created in irrzon. First half holds mapping from irr zone to full zone; part of second half holds repetition number info. work1 is a work array to keep track of grid points found so far. In case nspden=2 and nsppol=1, one has to take care of antiferromagnetic operations. The subgroup of non-magnetic operations is used to generate irrzon(:,:,2) and phnons(:,:,2), while the full group is used to generate irrzon(:,:,1) and phnons(:,:,1)
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
nspden=number of spin-density components nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetry elements in group n1,n2,n3=box dimensions of real space grid (or fft box) symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry matrices in real space (integers) tnons(3,nsym)=reduced nonsymmorphic translations (symrel and tnons are in terms of real space primitive translations)
OUTPUT
irrzon(n1*n2*n3,2+(nspden/4),(nspden/nsppol)-3*(nspden/4))=integer array which contains the locations of related grid points and the repetition number for each symmetry class. phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4))=phases associated with nonsymmorphic translations
PARENTS
m_ab7_kpoints,setsym,wfd_mkrho
CHILDREN
sort_int,wrtout
SOURCE
54 #if defined HAVE_CONFIG_H 55 #include "config.h" 56 #endif 57 58 #include "abi_common.h" 59 60 subroutine irrzg(irrzon,nspden,nsppol,nsym,n1,n2,n3,phnons,symafm,symrel,tnons) 61 62 use defs_basis 63 use m_profiling_abi 64 use m_errors 65 use m_sort 66 67 !This section has been created automatically by the script Abilint (TD). 68 !Do not modify the following lines by hand. 69 #undef ABI_FUNC 70 #define ABI_FUNC 'irrzg' 71 use interfaces_14_hidewrite 72 !End of the abilint section 73 74 implicit none 75 76 !Arguments ------------------------------------ 77 !scalars 78 integer,intent(in) :: n1,n2,n3,nspden,nsppol,nsym 79 !arrays 80 integer,intent(in) :: symafm(nsym),symrel(3,3,nsym) 81 integer,intent(out) :: irrzon(n1*n2*n3,2,(nspden/nsppol)-3*(nspden/4)) 82 real(dp),intent(in) :: tnons(3,nsym) 83 real(dp),intent(out) :: phnons(2,n1*n2*n3,(nspden/nsppol)-3*(nspden/4)) 84 85 !Local variables------------------------------- 86 !scalars 87 integer :: i1,i2,i3,id1,id2,id3,ifft,imagn,ind1,ind2,ipt,irep,isym,izone 88 integer :: izonemax,j1,j2,j3,jj,k1,k2,k3,l1,l2,l3,nfftot,npt,nsym_used 89 integer :: nzone,setzer,sppoldbl 90 real(dp) :: arg,ph1i,ph1r,ph2i,ph2r,tau1,tau2,tau3 91 logical,parameter :: afm_noncoll=.true. ! TRUE if antiferro symmetries are used in non-collinear magnetism 92 character(len=500) :: message 93 !arrays 94 integer,allocatable :: class(:),iperm(:),symafm_used(:),symrel_used(:,:,:) 95 integer,allocatable :: work1(:) 96 real(dp),allocatable :: tnons_used(:,:),work2(:,:) 97 98 ! ************************************************************************* 99 100 ABI_ALLOCATE(class,(nsym)) 101 ABI_ALLOCATE(iperm,(nsym)) 102 ABI_ALLOCATE(work1,(n1*n2*n3)) 103 ABI_ALLOCATE(work2,(2,n1*n2*n3)) 104 105 nfftot=n1*n2*n3 106 107 id1=n1/2+2 108 id2=n2/2+2 109 id3=n3/2+2 110 111 sppoldbl=nspden/nsppol;if (nspden==4) sppoldbl=1 112 113 do imagn=1,sppoldbl 114 115 ! Treat in a similar way the case of the full group and the non-magnetic subgroup 116 nsym_used=0 117 do isym=1,nsym 118 if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. & 119 & ((nspden==4).and.afm_noncoll) )then 120 nsym_used=nsym_used+1 121 end if 122 end do 123 124 if(imagn==2 .and. nsym_used/=nsym/2)then 125 write(message, '(a,a,a,a,a,i4,a,i0)' )& 126 & ' The number of ferromagnetic symmetry operations must be',ch10,& 127 & ' half the total number of operations, while it is observed that',ch10,& 128 & ' nsym=',nsym,' and nsym_magn=',nsym_used 129 MSG_BUG(message) 130 end if 131 132 ABI_ALLOCATE(symafm_used,(nsym_used)) 133 ABI_ALLOCATE(symrel_used,(3,3,nsym_used)) 134 ABI_ALLOCATE(tnons_used,(3,nsym_used)) 135 136 nsym_used=0 137 do isym=1,nsym 138 if( (imagn==1 .and. sppoldbl==2) .or. symafm(isym)==1 .or. & 139 & ((nspden==4).and.afm_noncoll) ) then 140 nsym_used=nsym_used+1 141 symrel_used(:,:,nsym_used)=symrel(:,:,isym) 142 tnons_used(:,nsym_used)=tnons(:,isym) 143 symafm_used(nsym_used)=symafm(isym) 144 end if 145 end do 146 if ((nspden/=4).or.(.not.afm_noncoll)) symafm_used=1 147 148 149 ! Zero out work array--later on, a zero entry will mean that 150 ! a given grid point has not yet been assigned to an ibz point 151 work1(1:nfftot)=0 152 irrzon(:,2,imagn)=0 153 154 ! Initialize at G=0 (always in irreducible zone) 155 nzone=1 156 irrzon(1,1,imagn)=1 157 irrzon(1,2,imagn)=nsym_used 158 ! Set phase exp(2*Pi*I*G dot tnons) for G=0 159 phnons(1,1,imagn)=one 160 phnons(2,1,imagn)=zero 161 npt=1 162 ! setting work1(1)=1 indicates that first grid point (G=0) is 163 ! in the iz (irreducible zone) 164 work1(1)=1 165 166 ind1=0 167 168 ! Loop over reciprocal space grid points: 169 do i3=1,n3 170 do i2=1,n2 171 do i1=1,n1 172 173 ind1=ind1+1 174 175 ! Check to see whether present grid point is equivalent to 176 ! any previously identified ibz point--if not, a new ibz point 177 ! has been found 178 179 if (work1(ind1)==0) then 180 181 ! A new class has been found. 182 183 ! Get location of G vector (grid point) centered at 0 0 0 184 l3=i3-(i3/id3)*n3-1 185 l2=i2-(i2/id2)*n2-1 186 l1=i1-(i1/id1)*n1-1 187 188 do isym=1,nsym_used 189 190 ! Get rotated G vector Gj for each symmetry element 191 ! -- here we use the TRANSPOSE of symrel_used; assuming symrel_used expresses 192 ! the rotation in real space, the transpose is then appropriate 193 ! for G space symmetrization (p. 1172d,e of notes, 2 June 1995). 194 j1=symrel_used(1,1,isym)*l1+& 195 & symrel_used(2,1,isym)*l2+symrel_used(3,1,isym)*l3 196 j2=symrel_used(1,2,isym)*l1+& 197 & symrel_used(2,2,isym)*l2+symrel_used(3,2,isym)*l3 198 j3=symrel_used(1,3,isym)*l1+& 199 & symrel_used(2,3,isym)*l2+symrel_used(3,3,isym)*l3 200 201 ! Map into [0,n-1] and then add 1 for array index in [1,n] 202 k1=1+mod(n1+mod(j1,n1),n1) 203 k2=1+mod(n2+mod(j2,n2),n2) 204 k3=1+mod(n3+mod(j3,n3),n3) 205 ! k1=1+map(j1,n1) 206 ! k2=1+map(j2,n2) 207 ! k3=1+map(j3,n3) 208 209 ! Get linear index of rotated point Gj 210 ind2=k1+n1*((k2-1)+n2*(k3-1)) 211 212 ! Store info for new class: 213 class(isym)=ind2 214 iperm(isym)=isym 215 216 ! Setting work array element to 1 indicates grid point has been 217 ! identified with iz point 218 work1(ind2)=1 219 220 ! End of loop on isym 221 end do 222 223 ! Sort integers into ascending order in each class 224 ! (this lumps together G vectors with the same linear index, i.e. 225 ! groups together symmetries which land on the same Gj) 226 call sort_int(nsym_used,class,iperm) 227 228 ! Check repetition factor (how many identical copies of Gj occur 229 ! from all symmetries applied to G) 230 irep=0 231 do isym=1,nsym_used 232 if (class(isym)==class(1)) then 233 irep=irep+1 234 end if 235 end do 236 ipt=nsym_used/irep 237 238 ! Repetition factor must be divisor of nsym_used: 239 if (nsym_used/=(ipt*irep)) then 240 write(message, '(a,i5,a,i6,a,a,a,a,a,a)' )& 241 & ' irep=',irep,' not a divisor of nsym_used=',nsym_used,ch10,& 242 & ' This usually indicates that',& 243 & ' the input symmetries do not form a group.',ch10,& 244 & ' Action : check the input symmetries carefully do they',& 245 & ' form a group ? If they do, there is a code bug.' 246 MSG_ERROR(message) 247 end if 248 249 ! Compute phases for any nonsymmorphic symmetries 250 ! exp(-2*Pi*I*G dot tau(j)) for each symmetry j with 251 ! (possibly zero) nonsymmorphic translation tau(j) 252 do jj=1,nsym_used 253 ! First get nonsymmorphic translation and see if nonzero 254 ! (iperm grabs the symmetries in the new order after sorting) 255 isym=iperm(jj) 256 tau1=tnons_used(1,isym) 257 tau2=tnons_used(2,isym) 258 tau3=tnons_used(3,isym) 259 if (abs(tau1)>tol12.or.abs(tau2)>tol12& 260 & .or.abs(tau3)>tol12) then 261 ! compute exp(-2*Pi*I*G dot tau) using original G 262 arg=two_pi*(dble(l1)*tau1+dble(l2)*tau2+dble(l3)*tau3) 263 work2(1,jj)=cos(arg) 264 work2(2,jj)=-sin(arg) 265 else 266 work2(1,jj)=one 267 work2(2,jj)=zero 268 end if 269 end do 270 271 ! All phases arising from symmetries which map to the same 272 ! G vector must actually be the same because 273 ! rho(Strans*G)=exp(2*Pi*I*(G) dot tau_S) rho(G) 274 ! must be satisfied; if exp(2*Pi*I*(G) dot tau_S) can be different 275 ! for two different symmetries S which both take G to the same St*G, 276 ! then the related Fourier components rho(St*G) must VANISH. 277 ! Thus: set "phase" to ZERO here in that case. 278 ! The G mappings occur in sets of irep members; if irep=1 then 279 ! all the G are unique. 280 ! MT 090212: 281 ! In the case of antiferromagn. symetries, one can have 282 ! rho(Strans*G)= -exp(2*Pi*I*(G) dot tau_S) rho(G) 283 ! (look at the minus !) 284 ! A special treatment is then operated on phons. 285 ! The later must be consistent with the use of phnons array 286 ! in symrhg.F90 routine. 287 ! XG 001108 : 288 ! Note that there is a tolerance on the 289 ! accuracy of tnons, especially when they are found from 290 ! the symmetry finder (with xred that might be a bit inaccurate) 291 if (irep > 1) then 292 do jj=1,nsym_used,irep 293 setzer=0 294 ph1r=work2(1,jj);ph1i=work2(2,jj) 295 do j1=jj,jj+irep-1 296 ph2r=work2(1,j1);ph2i=work2(2,j1) 297 if (((ph2r+ph1r)**2+(ph2i+ph1i)**2) <= tol14) then 298 if (setzer/=1) setzer=-1 299 else if (((ph2r-ph1r)**2+(ph2i-ph1i)**2) > tol14) then 300 setzer=1 301 end if 302 end do 303 ! Setzer= 0: phnons are all equal 304 ! Setzer=-1: phnons are equal in absolute value 305 ! Setzer= 1: some phnons are different 306 if (setzer/=0) then 307 if (setzer==-1) then 308 if (afm_noncoll.and.nspden==4) then 309 arg=symafm_used(iperm(jj)) 310 if (all(symafm_used(iperm(jj:jj+irep-1))==arg)) then 311 setzer=1 312 else 313 do j1=jj,jj+irep-1 314 work2(:,j1)=work2(:,j1)*dble(symafm_used(iperm(j1))) 315 end do 316 end if 317 else 318 setzer=1 319 end if 320 end if 321 if (setzer==1) work2(:,jj:jj+irep-1)=zero 322 end if 323 end do 324 ! Compress data if irep>1: 325 jj=0 326 do isym=1,nsym_used,irep 327 jj=jj+1 328 class(jj)=class(isym) 329 work2(1,jj)=work2(1,isym) 330 work2(2,jj)=work2(2,isym) 331 end do 332 end if 333 334 ! Put new unique points into irrzon array: 335 irrzon(1+npt:ipt+npt,1,imagn)=class(1:ipt) 336 337 ! Put repetition number into irrzon array: 338 irrzon(1+nzone,2,imagn)=irep 339 340 ! DEBUG 341 ! write(std_out,'(a,6i7)' )' irrzg : izone,i1,i2,i3,imagn,irrzon(859,2,1)=',& 342 ! & 1+nzone,i1,i2,i3,imagn,irrzon(859,2,1) 343 ! ENDDEBUG 344 345 ! Put phases (or 0) in phnons array: 346 phnons(:,1+npt:ipt+npt,imagn)=work2(:,1:ipt) 347 348 ! Update number of points in irrzon array: 349 ! (irep must divide evenly into nsym_used !) 350 npt=npt+ipt 351 352 ! Update number of classes: 353 nzone=nzone+1 354 355 end if 356 ! 357 ! End of loop on reciprocal space points, with indices i1, i2, i3 358 end do 359 end do 360 end do 361 362 if (allocated(symafm_used)) then 363 ABI_DEALLOCATE(symafm_used) 364 end if 365 if (allocated(symrel_used)) then 366 ABI_DEALLOCATE(symrel_used) 367 end if 368 if (allocated(tnons_used)) then 369 ABI_DEALLOCATE(tnons_used) 370 end if 371 372 end do ! imagn 373 374 !Make sure number of real space points accounted for equals actual number of grid points 375 if (npt/=n1*n2*n3) then 376 write(message, '(a,a,a,a,i10,a,i10,a,a,a,a,a,a,a,a,a)' ) ch10,& 377 & ' irrzg : ERROR -',ch10,& 378 & ' npt=',npt,' and n1*n2*n3=',n1*n2*n3,' are not equal',ch10,& 379 & ' This says that the total of all points in the irreducible',& 380 & ' sector in real space',ch10,& 381 & ' and all symmetrically equivalent',& 382 & ' points, npt, does not equal the actual number',ch10,& 383 & ' of real space grid points.' 384 call wrtout(std_out,message,'COLL') 385 write(message,'(3a)') & 386 & ' This may mean that the input symmetries do not form a group',ch10,& 387 & ' Action : check input symmetries carefully for errors.' 388 MSG_ERROR(message) 389 end if 390 391 !Perform some checks 392 do imagn=1,sppoldbl 393 394 do ifft=1,nfftot 395 if (irrzon(ifft,1,imagn)<1.or.irrzon(ifft,1,imagn)>nfftot) then 396 write(message,'(a,4i0,a,a)')& 397 & ' ifft,irrzon(ifft,1,imagn),nfftot,imagn=',ifft,irrzon(ifft,1,imagn),nfftot,imagn,ch10,& 398 & ' =>irrzon goes outside acceptable bounds.' 399 MSG_BUG(message) 400 end if 401 end do 402 403 izonemax=0 404 do izone=1,nfftot 405 ! Global bounds 406 if (irrzon(izone,2,imagn)<0.or.irrzon(izone,2,imagn)>(nsym/imagn)) then 407 write(message, '(a,5i7,a,a)' )& 408 & ' izone,nzone,irrzon(izone,2,imagn),nsym,imagn =',izone,nzone,irrzon(izone,2,imagn),nsym,imagn,ch10,& 409 & ' =>irrzon goes outside acceptable bounds.' 410 MSG_BUG(message) 411 end if 412 ! Second index only goes up to nzone 413 if(izonemax==0)then 414 if (irrzon(izone,2,imagn)==0)izonemax=izone-1 415 end if 416 if(izonemax/=0)then 417 if (irrzon(izone,2,imagn)/=0) then 418 message = ' beyond izonemax, irrzon(izone,2,imagn) should be zero' 419 MSG_BUG(message) 420 end if 421 end if 422 end do 423 424 end do ! imagn 425 426 ABI_DEALLOCATE(class) 427 ABI_DEALLOCATE(iperm) 428 ABI_DEALLOCATE(work1) 429 ABI_DEALLOCATE(work2) 430 431 end subroutine irrzg