TABLE OF CONTENTS
ABINIT/testkgrid [ Functions ]
NAME
testkgrid
FUNCTION
Test different grids of k points. The algorithm used is based on the idea of testing different one-dimensional sets of possible k point grids. It is not exhaustive (other families could be included), but should do a respectable job in all cases. The Monkhorst-Pack set of grids (defined with respect to symmetry axes, and not primitive axes) is always tested.
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
bravais(11): bravais(1)=iholohedry
bravais(2)=center
bravais(3:11)=coordinates of rprim in the axes of the conventional bravais lattice (*2 if center/=0)
iout=unit number for echoed output
msym=default maximal number of symmetries
nsym=number of symetries
prtkpt=if non-zero, will write the characteristics of k grids, then stop
rprimd(3,3)=dimensional real space primitive translations (bohr)
symafm(nsym)=(anti)ferromagnetic part of symmetry operations
symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations
vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum
OUTPUT
kptrlatt(3,3)=k-point lattice specification nshiftk=number of k-point shifts in shiftk (always 1 from this routine) shiftk(3,210)=shift vectors for k point generation
SIDE EFFECTS
kptrlen=length of the smallest real space supercell vector associated with the lattice of k points.
NOTES
Note that nkpt can be computed by calling this routine with input value nkpt=0 Note that kptopt is always =1 in this routine.
PARENTS
inkpts,m_ab7_kpoints
CHILDREN
getkgrid,leave_new,matr3inv,metric,smallprim,wrtout,xmpi_abort
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 subroutine testkgrid(bravais,iout,kptrlatt,kptrlen,& 60 & msym,nshiftk,nsym,prtkpt,rprimd,shiftk,symafm,symrel,vacuum) 61 62 use defs_basis 63 use m_errors 64 use m_profiling_abi 65 use m_xmpi 66 67 use m_geometry, only : metric 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'testkgrid' 73 use interfaces_14_hidewrite 74 use interfaces_16_hideleave 75 use interfaces_32_util 76 use interfaces_41_geometry 77 use interfaces_56_recipspace, except_this_one => testkgrid 78 !End of the abilint section 79 80 implicit none 81 82 !Arguments ------------------------------------ 83 !scalars 84 integer,intent(in) :: iout,msym,nsym,prtkpt 85 integer,intent(out) :: nshiftk 86 real(dp),intent(inout) :: kptrlen 87 !arrays 88 integer,intent(in) :: bravais(11),symafm(msym),symrel(3,3,msym),vacuum(3) 89 integer,intent(out) :: kptrlatt(3,3) 90 real(dp),intent(in) :: rprimd(3,3) 91 real(dp),intent(inout) :: shiftk(3,210) !vz_i 92 93 !Local variables------------------------------- 94 !scalars 95 integer,parameter :: kptopt=1,mkpt_list=100000 96 integer :: ang90,center,dirvacuum,equal,igrid,igrid_current,iholohedry,ii,init_mult,iscale,iscf 97 integer :: iset,mult1,mult2,mult3,ndims,nkpt,nkpt_current,nkpt_trial,nset 98 real(dp) :: buffer_scale,determinant,fact,factor,kptrlen_current,kptrlen_max,kptrlen_target 99 real(dp) :: kptrlen_trial,length1,length2,length3,length_axis1,length_axis2 100 real(dp) :: length_axis3,merit_factor,mult1h,mult2h,mult3h,reduceda,reducedb 101 real(dp) :: sca,scb,scc,surface,ucvol 102 character(len=500) :: message 103 !arrays 104 integer :: kptrlatt_current(3,3),kptrlatt_trial(3,3) 105 integer,allocatable :: grid_list(:) 106 real(dp) :: axes(3,3),gmet(3,3),gprimd(3,3),matrix1(3,3),matrix2(3,3) 107 real(dp) :: metmin(3,3),minim(3,3),r2d(3,3),rmet(3,3),rsuper(3,3) 108 real(dp) :: shiftk_current(3,210),shiftk_trial(3,210) 109 real(dp),allocatable :: kpt(:,:),kptrlen_list(:),wtk(:) 110 111 ! ************************************************************************* 112 113 kptrlen_target=kptrlen 114 115 !The vacuum array must be made of 0 or 1 116 do ii=1,3 117 if(vacuum(ii)/=0 .and. vacuum(ii)/=1)then 118 write(message,'(a,a,a,i1,a,i3,a,a)')& 119 & 'The values of vacuum must be 0 or 1.',ch10,& 120 & 'However, the input vacuum(',ii,') is',vacuum(ii),ch10,& 121 & 'Action : correct vacuum in your input file.' 122 MSG_ERROR(message) 123 end if 124 end do 125 126 !Specific preparation for 2-dimensional system 127 if(sum(vacuum(:))==1)then 128 129 ! Make the non-active vector orthogonal to the active vectors, 130 ! and take it along the z direction 131 if(vacuum(1)==1)then 132 r2d(1,3)=rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3) 133 r2d(2,3)=rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3) 134 r2d(3,3)=rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3) 135 r2d(:,1)=rprimd(:,2) 136 r2d(:,2)=rprimd(:,3) 137 dirvacuum=1 138 else if(vacuum(2)==1)then 139 r2d(1,3)=rprimd(2,3)*rprimd(3,1)-rprimd(3,3)*rprimd(2,1) 140 r2d(2,3)=rprimd(3,3)*rprimd(1,1)-rprimd(1,3)*rprimd(3,1) 141 r2d(3,3)=rprimd(1,3)*rprimd(2,1)-rprimd(2,3)*rprimd(1,1) 142 r2d(:,1)=rprimd(:,3) 143 r2d(:,2)=rprimd(:,1) 144 dirvacuum=2 145 else if(vacuum(3)==1)then 146 r2d(1,3)=rprimd(2,1)*rprimd(3,2)-rprimd(3,1)*rprimd(2,2) 147 r2d(2,3)=rprimd(3,1)*rprimd(1,2)-rprimd(1,1)*rprimd(3,2) 148 r2d(3,3)=rprimd(1,1)*rprimd(2,2)-rprimd(2,1)*rprimd(1,2) 149 r2d(:,1)=rprimd(:,1) 150 r2d(:,2)=rprimd(:,2) 151 dirvacuum=3 152 end if 153 surface=sqrt(sum(r2d(:,3)**2)) 154 ! Identify the 2-D Bravais lattice 155 ! DEBUG 156 ! write(std_out,*)' r2d=',r2d(:,:) 157 ! ENDDEBUG 158 call metric(gmet,gprimd,-1,rmet,r2d,ucvol) 159 call smallprim(metmin,minim,r2d) 160 ! DEBUG 161 ! write(std_out,*)' minim=',minim(:,:) 162 ! ENDDEBUG 163 ang90=0 ; equal=0 ; center=0 164 axes(:,:)=minim(:,:) 165 if(abs(metmin(1,2))<tol8)ang90=1 166 if(abs(metmin(1,1)-metmin(2,2))<tol8)equal=1 167 if(ang90==1)then 168 if(equal==1)iholohedry=4 169 if(equal==0)iholohedry=2 170 else if(equal==1)then 171 reduceda=metmin(1,2)/metmin(1,1) 172 if(abs(reduceda+0.5_dp)<tol8)then 173 iholohedry=3 174 else if(abs(reduceda-0.5_dp)<tol8)then 175 iholohedry=3 176 ! Use conventional axes 177 axes(:,2)=minim(:,2)-minim(:,1) 178 else 179 iholohedry=2 ; center=1 180 axes(:,1)=minim(:,1)+minim(:,2) 181 axes(:,2)=minim(:,2)-minim(:,1) 182 end if 183 else 184 reduceda=metmin(1,2)/metmin(1,1) 185 reducedb=metmin(1,2)/metmin(2,2) 186 if(abs(reduceda+0.5_dp)<tol8)then 187 iholohedry=2 ; center=1 188 axes(:,2)=2.0_dp*minim(:,2)+minim(:,1) 189 else if(abs(reduceda-0.5_dp)<tol8)then 190 iholohedry=2 ; center=1 191 axes(:,2)=2.0_dp*minim(:,2)-minim(:,1) 192 else if(abs(reducedb+0.5_dp)<tol8)then 193 iholohedry=2 ; center=1 194 axes(:,1)=2.0_dp*minim(:,1)+minim(:,2) 195 else if(abs(reducedb-0.5_dp)<tol8)then 196 iholohedry=2 ; center=1 197 axes(:,1)=2.0_dp*minim(:,1)-minim(:,2) 198 else 199 iholohedry=1 200 end if 201 end if 202 ! Make sure that axes form a right-handed coordinate system 203 determinant=axes(1,1)*axes(2,2)*axes(3,3) & 204 & +axes(1,2)*axes(2,3)*axes(3,1) & 205 & +axes(1,3)*axes(3,2)*axes(2,1) & 206 & -axes(1,1)*axes(3,2)*axes(2,3) & 207 & -axes(1,3)*axes(2,2)*axes(3,1) & 208 & -axes(1,2)*axes(2,1)*axes(3,3) 209 if(determinant<zero)then 210 axes(:,1)=-axes(:,1) 211 end if 212 ! Prefer symmetry axes on the same side as the primitive axes 213 sca=axes(1,1)*r2d(1,1)+axes(2,1)*r2d(2,1)+axes(3,1)*r2d(3,1) 214 scb=axes(1,2)*r2d(1,2)+axes(2,2)*r2d(2,2)+axes(3,2)*r2d(3,2) 215 scc=axes(1,3)*rprimd(1,dirvacuum)& 216 & +axes(2,3)*rprimd(2,dirvacuum)& 217 & +axes(3,3)*rprimd(3,dirvacuum) 218 if(sca<-tol8 .and. scb<-tol8)then 219 axes(:,1)=-axes(:,1) ; sca=-sca 220 axes(:,2)=-axes(:,2) ; scb=-scb 221 end if 222 ! Doing this might change the angle between vectors, so that 223 ! the cell is not conventional anymore 224 ! if(sca<-tol8 .and. scc<-tol8)then 225 ! axes(:,1)=-axes(:,1) ; sca=-sca 226 ! axes(:,3)=-axes(:,3) ; scc=-scc 227 ! end if 228 ! if(scb<-tol8 .and. scc<-tol8)then 229 ! axes(:,2)=-axes(:,2) ; scb=-scb 230 ! axes(:,3)=-axes(:,3) ; scc=-scc 231 ! end if 232 length_axis1=sqrt(axes(1,1)**2+axes(2,1)**2+axes(3,1)**2) 233 length_axis2=sqrt(axes(1,2)**2+axes(2,2)**2+axes(3,2)**2) 234 235 ! DEBUG 236 ! write(std_out,*)' testkgrid : iholohedry, center =',iholohedry,center 237 ! write(std_out,*)' testkgrid : axis 1=',axes(:,1) 238 ! write(std_out,*)' testkgrid : axis 2=',axes(:,2) 239 ! write(std_out,*)' testkgrid : axis 3=',axes(:,3) 240 ! write(std_out,*)' testkgrid : length_axis=',length_axis1,length_axis2 241 ! ENDDEBUG 242 243 ! End special treatment of 2-D case 244 end if 245 246 !3-dimensional system 247 if(sum(vacuum(:))==0)then 248 iholohedry=bravais(1) 249 center=bravais(2) 250 fact=1.0_dp 251 if(center/=0)fact=0.5_dp 252 matrix1(:,1)=bravais(3:5)*fact 253 matrix1(:,2)=bravais(6:8)*fact 254 matrix1(:,3)=bravais(9:11)*fact 255 call matr3inv(matrix1,matrix2) 256 do ii=1,3 257 axes(:,ii)=rprimd(:,1)*matrix2(ii,1)+rprimd(:,2)*matrix2(ii,2)+rprimd(:,3)*matrix2(ii,3) 258 end do 259 length_axis1=sqrt(axes(1,1)**2+axes(2,1)**2+axes(3,1)**2) 260 length_axis2=sqrt(axes(1,2)**2+axes(2,2)**2+axes(3,2)**2) 261 length_axis3=sqrt(axes(1,3)**2+axes(2,3)**2+axes(3,3)**2) 262 ! DEBUG 263 ! write(std_out,*)' testkgrid : axes=',axes(:,:) 264 ! write(std_out,*)' length_axis=',length_axis1,length_axis2,length_axis3 265 ! ENDDEBUG 266 end if 267 268 !This routine examine only primitive k lattices. 269 nshiftk=1 270 271 !If prtkpt/=0, will examine more grids than strictly needed 272 buffer_scale=one 273 if(prtkpt/=0)buffer_scale=two 274 275 if(prtkpt/=0)then 276 write(message,'(a,a,a,a,a,a,a,a)' )ch10,& 277 & ' testkgrid : will perform the analysis of a series of k-grids.',ch10,& 278 & ' Note that kptopt=1 in this analysis, irrespective of its input value.',ch10,ch10,& 279 & ' Grid# kptrlatt shiftk kptrlen nkpt iset',ch10 280 call wrtout(std_out,message,'COLL') 281 call wrtout(iout,message,'COLL') 282 ABI_ALLOCATE(grid_list,(mkpt_list)) 283 ABI_ALLOCATE(kptrlen_list,(mkpt_list)) 284 grid_list(:)=0 285 kptrlen_list(:)=0.0_dp 286 end if 287 288 if(sum(vacuum(:))==3)then 289 290 kptrlatt(:,:)=0 291 kptrlatt(1,1)=1 292 kptrlatt(2,2)=1 293 kptrlatt(3,3)=1 294 shiftk(:,1)=0.0_dp 295 kptrlen=1000.0_dp 296 nkpt_current=1 297 igrid_current=1 298 299 if(prtkpt/=0)then 300 write(message,& 301 & '(a,3i4,a,es14.4,a,es14.4,i8,i6,a,a,3i4,a,es14.4,a,a,3i4,a,es14.4,a)' )& 302 & ' 1 ',kptrlatt(:,1),' ',shiftk(1,1),' ',kptrlen,1,1,ch10,& 303 & ' ',kptrlatt(:,2),' ',shiftk(2,1),ch10,& 304 & ' ',kptrlatt(:,3),' ',shiftk(3,1),ch10 305 call wrtout(std_out,message,'COLL') 306 call wrtout(iout,message,'COLL') 307 ! The unit cell volume is fake 308 ucvol=kptrlen**3 309 end if 310 311 else 312 313 nkpt=0 ; nkpt_current=0 ; iscf=1 ; iset=1 314 kptrlen_current=0.0_dp 315 mult1=0 ; mult2=0 ; mult3=0 ; init_mult=1 316 ABI_ALLOCATE(kpt,(3,nkpt)) 317 ABI_ALLOCATE(wtk,(nkpt)) 318 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 319 320 ! Loop on different grids, the upper limit is only to avoid an infinite loop 321 do igrid=1,1000 322 323 kptrlatt_trial(:,:)=0 324 kptrlatt_trial(1,1)=1 325 kptrlatt_trial(2,2)=1 326 kptrlatt_trial(3,3)=1 327 shiftk_trial(:,1)=0.0_dp 328 329 ! 1-dimensional system 330 if(sum(vacuum(:))==2)then 331 if(vacuum(1)==0)then 332 kptrlatt_trial(1,1)=2*igrid ; shiftk_trial(1,1)=0.5_dp 333 else if(vacuum(2)==0)then 334 kptrlatt_trial(2,2)=2*igrid ; shiftk_trial(2,1)=0.5_dp 335 else if(vacuum(3)==0)then 336 kptrlatt_trial(3,3)=2*igrid ; shiftk_trial(3,1)=0.5_dp 337 end if 338 end if 339 340 ! 2-dimensional system 341 if(sum(vacuum(:))==1)then 342 343 ! Treat hexagonal holohedries separately 344 if(iholohedry==3)then 345 346 ! write(std_out,*)' testkgrid : 2D, hexagonal' 347 348 mult1=mult1+1 349 nset=4 350 if(iset==1)then 351 rsuper(:,1)=axes(:,1)*mult1 352 rsuper(:,2)=axes(:,2)*mult1 353 shiftk_trial(:,1)=0.0_dp 354 else if(iset==2)then 355 rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1 356 rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1 357 shiftk_trial(1,1)=1.0_dp/3.0_dp 358 shiftk_trial(2,1)=1.0_dp/3.0_dp 359 else if(iset==3)then 360 rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1 361 rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1 362 shiftk_trial(:,1)=0.0_dp 363 else if(iset==4)then 364 rsuper(:,1)=axes(:,1)*mult1 365 rsuper(:,2)=axes(:,2)*mult1 366 shiftk_trial(1,1)=0.5_dp 367 shiftk_trial(2,1)=0.5_dp 368 end if 369 370 else 371 ! Now treat all other holohedries 372 length1=length_axis1*mult1 373 length2=length_axis2*mult2 374 ! DEBUG 375 ! write(std_out,*)' testkgrid : (2d) length=',length1,length2 376 ! ENDDEBUG 377 if(abs(length1-length2)<tol8)then 378 mult1=mult1+1 379 mult2=mult2+1 380 else if(length1>length2)then 381 mult2=mult2+1 382 else if(length2>length1)then 383 mult1=mult1+1 384 end if 385 nset=4 386 ! iset==5 and 6 are allowed only for centered lattice 387 if(center==1)nset=6 388 if(iset==1 .or. iset==2)then 389 rsuper(:,1)=axes(:,1)*mult1 390 rsuper(:,2)=axes(:,2)*mult2 391 else if(iset==3 .or. iset==4)then 392 rsuper(:,1)=axes(:,1)*mult1-axes(:,2)*mult2 393 rsuper(:,2)=axes(:,1)*mult1+axes(:,2)*mult2 394 else if(iset==5 .or. iset==6)then 395 rsuper(:,1)=axes(:,1)*(mult1-0.5_dp)-axes(:,2)*(mult2-0.5_dp) 396 rsuper(:,2)=axes(:,1)*(mult1-0.5_dp)+axes(:,2)*(mult2-0.5_dp) 397 end if 398 ! This was the easiest way to code all even mult1 and mult2 pairs : 399 ! make separate series for this possibility. 400 if(iset==2 .or. iset==4 .or. iset==6)then 401 rsuper(:,1)=2.0_dp*rsuper(:,1) 402 rsuper(:,2)=2.0_dp*rsuper(:,2) 403 end if 404 shiftk_trial(1,1)=0.5_dp 405 shiftk_trial(2,1)=0.5_dp 406 407 end if 408 409 ! Put back the inactive direction 410 if(dirvacuum==1)then 411 rsuper(:,3)=rsuper(:,1) 412 shiftk_trial(3,1)=shiftk_trial(1,1) 413 rsuper(:,1)=rprimd(:,1) 414 shiftk_trial(1,1)=0.0_dp 415 else if(dirvacuum==2)then 416 rsuper(:,3)=rsuper(:,1) 417 shiftk_trial(3,1)=shiftk_trial(1,1) 418 rsuper(:,1)=rsuper(:,2) 419 shiftk_trial(1,1)=shiftk_trial(2,1) 420 rsuper(:,2)=rprimd(:,2) 421 shiftk_trial(2,1)=0.0_dp 422 else if(dirvacuum==3)then 423 rsuper(:,3)=rprimd(:,3) 424 shiftk_trial(3,1)=0.0_dp 425 end if 426 427 ! The supercell and the corresponding shift have been generated ! 428 ! Convert cartesian coordinates into kptrlatt_trial 429 do ii=1,3 430 kptrlatt_trial(:,ii)=nint( gprimd(1,:)*rsuper(1,ii)+& 431 & gprimd(2,:)*rsuper(2,ii)+& 432 & gprimd(3,:)*rsuper(3,ii) ) 433 end do 434 435 ! End of 2-dimensional system 436 end if 437 438 ! 3-dimensional system 439 if(sum(vacuum(:))==0)then 440 ! Treat hexagonal holohedries separately 441 if(iholohedry==6)then 442 length1=length_axis1*mult1 443 length3=length_axis3*mult3 444 ! DEBUG 445 ! write(std_out,*)' testkgrid : (hex) lengths=',length1,length2 446 ! ENDDEBUG 447 if(abs(length1-length3)<tol8)then 448 mult1=mult1+1 449 mult3=mult3+1 450 else if(length1>length3)then 451 mult3=mult3+1 452 else if(length3>length1)then 453 mult1=mult1+1 454 end if 455 nset=4 456 if(iset==1)then 457 rsuper(:,1)=axes(:,1)*mult1 458 rsuper(:,2)=axes(:,2)*mult1 459 rsuper(:,3)=axes(:,3)*mult3 460 shiftk_trial(:,1)=0.0_dp 461 shiftk_trial(3,1)=0.5_dp 462 else if(iset==2)then 463 rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1 464 rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1 465 rsuper(:,3)=axes(:,3)*mult3 466 shiftk_trial(1,1)=1.0_dp/3.0_dp 467 shiftk_trial(2,1)=1.0_dp/3.0_dp 468 shiftk_trial(3,1)=0.5_dp 469 else if(iset==3)then 470 rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1 471 rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1 472 rsuper(:,3)=axes(:,3)*mult3 473 shiftk_trial(:,1)=0.0_dp 474 shiftk_trial(3,1)=0.5_dp 475 else if(iset==4)then 476 rsuper(:,1)=axes(:,1)*mult1 477 rsuper(:,2)=axes(:,2)*mult1 478 rsuper(:,3)=axes(:,3)*mult3 479 shiftk_trial(:,1)=0.5_dp 480 end if 481 482 else 483 ! Now treat all other holohedries 484 length1=length_axis1*mult1 485 length2=length_axis2*mult2 486 length3=length_axis3*mult3 487 ! DEBUG 488 ! write(std_out,*)' testkgrid : length=',length1,length2,length3 489 ! ENDDEBUG 490 if(length2>length1+tol8 .and. length3>length1+tol8)then 491 mult1=mult1+1 492 else if(length1>length2+tol8 .and. length3>length2+tol8)then 493 mult2=mult2+1 494 else if(length1>length3+tol8 .and. length2>length3+tol8)then 495 mult3=mult3+1 496 else if(abs(length2-length3)<tol8 .and. & 497 & abs(length1-length3)<tol8 .and. & 498 & abs(length1-length2)<tol8 )then 499 mult1=mult1+1 ; mult2=mult2+1 ; mult3=mult3+1 500 else if(abs(length1-length2)<tol8)then 501 mult1=mult1+1 ; mult2=mult2+1 502 else if(abs(length1-length3)<tol8)then 503 mult1=mult1+1 ; mult3=mult3+1 504 else if(abs(length2-length3)<tol8)then 505 mult2=mult2+1 ; mult3=mult3+1 506 end if 507 nset=6 508 if(center==-1 .or. center==-3)nset=8 509 if(iset==1 .or. iset==2)then 510 ! Simple lattice of k points 511 rsuper(:,1)=axes(:,1)*mult1 512 rsuper(:,2)=axes(:,2)*mult2 513 rsuper(:,3)=axes(:,3)*mult3 514 shiftk_trial(:,1)=0.5_dp 515 else if(iset==3 .or. iset==4)then 516 ! FCC lattice of k points = BCC lattice in real space 517 rsuper(:,1)=-axes(:,1)*mult1+axes(:,2)*mult2+axes(:,3)*mult3 518 rsuper(:,2)= axes(:,1)*mult1-axes(:,2)*mult2+axes(:,3)*mult3 519 rsuper(:,3)= axes(:,1)*mult1+axes(:,2)*mult2-axes(:,3)*mult3 520 shiftk_trial(:,1)=0.5_dp 521 else if(iset==5 .or. iset==6)then 522 ! BCC lattice of k points = FCC lattice in real space 523 rsuper(:,1)= axes(:,2)*mult2+axes(:,3)*mult3 524 rsuper(:,2)= axes(:,1)*mult1 +axes(:,3)*mult3 525 rsuper(:,3)= axes(:,1)*mult1+axes(:,2)*mult2 526 ! The BCC lattice has no empty site with full symmetry 527 shiftk_trial(:,1)=0.0_dp 528 else if(iset==7 .or. iset==8)then 529 ! iset==7 and 8 are allowed only for centered lattice 530 mult1h=mult1-0.5_dp 531 mult2h=mult2-0.5_dp 532 mult3h=mult3-0.5_dp 533 if(center==-1)then 534 ! FCC lattice of k points = BCC lattice in real space 535 rsuper(:,1)=-axes(:,1)*mult1h+axes(:,2)*mult2h+axes(:,3)*mult3h 536 rsuper(:,2)= axes(:,1)*mult1h-axes(:,2)*mult2h+axes(:,3)*mult3h 537 rsuper(:,3)= axes(:,1)*mult1h+axes(:,2)*mult2h-axes(:,3)*mult3h 538 shiftk_trial(:,1)=0.5_dp 539 else if(center==-3)then 540 ! BCC lattice of k points = FCC lattice in real space 541 rsuper(:,1)= axes(:,2)*mult2h+axes(:,3)*mult3h 542 rsuper(:,2)= axes(:,1)*mult1h +axes(:,3)*mult3h 543 rsuper(:,3)= axes(:,1)*mult1h+axes(:,2)*mult2h 544 ! The BCC lattice has no empty site with full symmetry 545 shiftk_trial(:,1)=0.0_dp 546 end if 547 end if 548 ! This was the easiest way to code all even mult1, mult2, mult3 triplets : 549 ! make separate series for this possibility. 550 if(2*(iset/2)==iset)then 551 rsuper(:,1)=2.0_dp*rsuper(:,1) 552 rsuper(:,2)=2.0_dp*rsuper(:,2) 553 rsuper(:,3)=2.0_dp*rsuper(:,3) 554 end if 555 end if 556 557 ! DEBUG 558 ! write(std_out,*)' testkgrid : gprimd=',gprimd(:,:) 559 ! write(std_out,*)' testkgrid : rsuper=',rsuper(:,:) 560 ! write(std_out,*)' testkgrid : iset =',iset 561 ! ENDDEBUG 562 563 564 ! The supercell and the corresponding shift have been generated ! 565 ! Convert cartesian coordinates into kptrlatt_trial 566 do ii=1,3 567 kptrlatt_trial(:,ii)=nint( gprimd(1,:)*rsuper(1,ii)+& 568 & gprimd(2,:)*rsuper(2,ii)+& 569 & gprimd(3,:)*rsuper(3,ii) ) 570 end do 571 572 ! End of 3-dimensional system 573 end if 574 575 ! DEBUG 576 ! write(std_out,*)' testkgrid : before getkgrid' 577 ! write(std_out,*)' testkgrid : rprimd=',rprimd(:,:) 578 ! write(std_out,*)' testkgrid : kptrlatt_trial=',kptrlatt_trial(:,:) 579 ! ENDDEBUG 580 581 call getkgrid(0,0,iscf,kpt,& 582 & kptopt,kptrlatt_trial,kptrlen_trial,& 583 & msym,nkpt,nkpt_trial,nshiftk,nsym,rprimd,& 584 & shiftk_trial,symafm,symrel,vacuum,wtk) 585 586 ! DEBUG 587 ! write(std_out,*)' testkgrid : after getkgrid' 588 ! ENDDEBUG 589 590 ! In case one does not need the full list of grids, will take a shortcut, and go to one of the last grids of the series, 591 ! that generates a kptrlen_trial that is just below kptrlen. 592 if(prtkpt==0 .and. init_mult==1 .and. kptrlen_trial<(half-tol8)*kptrlen )then 593 iscale=int((one-tol8)*kptrlen/kptrlen_trial) 594 mult1=mult1*iscale 595 mult2=mult2*iscale 596 mult3=mult3*iscale 597 init_mult=0 598 ! DEBUG 599 ! write(std_out,*)' testkgrid : iscale=',iscale 600 ! ENDDEBUG 601 kptrlatt_trial(:,:)=kptrlatt_trial(:,:)*iscale 602 call getkgrid(0,0,iscf,kpt,& 603 & kptopt,kptrlatt_trial,kptrlen_trial,& 604 & msym,nkpt,nkpt_trial,nshiftk,nsym,rprimd,& 605 & shiftk_trial,symafm,symrel,vacuum,wtk) 606 end if 607 608 if( (kptrlen_trial+tol8>kptrlen*(1.0_dp+tol8) .and. nkpt_current==0) .or. & 609 & (kptrlen_trial+tol8>kptrlen*(1.0_dp+tol8) .and. nkpt_trial<nkpt_current) .or. & 610 & (nkpt_trial==nkpt_current .and. kptrlen_trial>kptrlen_current*(1.0_dp+tol8)))then 611 612 kptrlatt_current(:,:)=kptrlatt_trial(:,:) 613 nkpt_current=nkpt_trial 614 shiftk_current(:,:)=shiftk_trial(:,:) 615 kptrlen_current=kptrlen_trial 616 igrid_current=igrid 617 end if 618 619 if(prtkpt/=0)then 620 write(message,'(i5,a,3i4,a,es14.4,a,es14.4,i8,i6,a,a,3i4,a,es14.4,a,a,3i4,a,es14.4,a)' )& 621 & igrid,' ',kptrlatt_trial(:,1),' ',shiftk_trial(1,1),& 622 & ' ',kptrlen_trial,nkpt_trial,iset,ch10,& 623 & ' ',kptrlatt_trial(:,2),' ',shiftk_trial(2,1),ch10,& 624 & ' ',kptrlatt_trial(:,3),' ',shiftk_trial(3,1),ch10 625 call wrtout(std_out,message,'COLL') 626 call wrtout(iout,message,'COLL') 627 628 ! Keep track of this grid, if it is worth 629 if(kptrlen_trial > kptrlen_list(nkpt_trial)*(1.0_dp+tol8))then 630 grid_list(nkpt_trial)=igrid 631 kptrlen_list(nkpt_trial)=kptrlen_trial 632 end if 633 end if 634 635 ! Treat 1-D case 636 if( sum(vacuum(:))==2 .and. kptrlen_trial>buffer_scale*(1.0_dp+tol8)*kptrlen )exit 637 638 ! Treat 2-D case or 3-D case 639 if( sum(vacuum(:))<=1 .and. kptrlen_trial>buffer_scale*(1.0_dp+tol8)*kptrlen )then 640 ! The present set of sets of k points is finished : 641 ! either it was the last, or one has to go to the next one 642 if(iset==nset)exit 643 iset=iset+1 644 mult1=0 ; mult2=0 ; mult3=0 ; init_mult=1 645 end if 646 647 end do ! igrid=1,1000 648 649 ABI_DEALLOCATE(kpt) 650 ABI_DEALLOCATE(wtk) 651 652 kptrlatt(:,:)=kptrlatt_current(:,:) 653 shiftk(:,:)=shiftk_current(:,:) 654 kptrlen=kptrlen_current 655 656 end if ! test on the number of dimensions 657 658 if(prtkpt/=0)then 659 660 ! sqrt(1/2) comes from the FCC packing, the best one 661 factor=sqrt(0.5_dp)/ucvol/dble(nsym) 662 ndims=3 663 if(sum(vacuum(:))/=0)then 664 if(sum(vacuum(:))==1)then 665 ! sqrt(3/4) comes from the hex packing, the best one 666 ! one multiplies by 2 because nsym is likely twice the number 667 ! of symmetries that can be effectively used in 2D 668 ndims=2 ; factor=sqrt(0.75_dp)/surface/dble(nsym)*2 669 write(message,'(2a)' )ch10,' Note that the system is bi-dimensional.' 670 else if(sum(vacuum(:))==2)then 671 ndims=1 ; factor=1/ucvol 672 write(message,'(2a)' )ch10,' Note that the system is uni-dimensional.' 673 else if(sum(vacuum(:))==3)then 674 ndims=0 675 write(message,'(2a)' )ch10,' Note that the system is zero-dimensional.' 676 end if 677 call wrtout(std_out,message,'COLL') 678 call wrtout(iout,message,'COLL') 679 end if 680 681 ! The asymptotic value of the merit factor is determined 682 ! by the set of symmetries : in 3D, if it includes the 683 ! inversion symmetry, the limit will be 1, if not, it 684 ! will be two. In 2D, if it includes the inversion symmetry 685 ! and an operation that maps z on -z, it will tend to one, 686 ! while if only one of these operations is present, 687 ! it will tend to two, and if none is present, it will tend to four. 688 write(message,'(11a)' )ch10,& 689 & ' List of best grids, ordered by nkpt.',ch10,& 690 & ' (stop at a value of kptrlen 20% larger than the target value).',ch10,& 691 & ' (the merit factor will tend to one or two in 3 dimensions)',ch10,& 692 & ' (and to one, two or four in 2 dimensions)',ch10,ch10,& 693 & ' nkpt kptrlen grid# merit_factor' 694 call wrtout(std_out,message,'COLL') 695 call wrtout(iout,message,'COLL') 696 697 kptrlen_max=0.0_dp 698 do ii=1,mkpt_list 699 if(kptrlen_list(ii)>kptrlen_max*(1.0_dp+tol8))then 700 kptrlen_max=kptrlen_list(ii) 701 merit_factor=kptrlen_max**ndims/dble(ii)*factor 702 write(message, '(i6,es14.4,i6,f12.4)' )ii,kptrlen_max,grid_list(ii),merit_factor 703 call wrtout(std_out,message,'COLL') 704 call wrtout(iout,message,'COLL') 705 end if 706 if(kptrlen_max>1.2_dp*(1.0_dp-tol8)*kptrlen_target)exit 707 end do 708 709 write(message,'(a,a,es14.4,a,a,i6,a,a,a,es14.4,a,i6)' )ch10,& 710 & ' For target kptrlen=',kptrlen_target,',',& 711 & ' the selected grid is number',igrid_current,',',ch10,& 712 & ' giving kptrlen=',kptrlen_current,' with nkpt=',nkpt_current 713 call wrtout(std_out,message,'COLL') 714 call wrtout(iout,message,'COLL') 715 716 write(message,'(a,a,a,a)' )ch10,& 717 & ' testkgrid : stop after analysis of a series of k-grids.',ch10,& 718 & ' For usual production runs, set prtkpt back to 0 (the default).' 719 call wrtout(std_out,message,'COLL',do_flush=.True.) 720 call wrtout(iout,message,'COLL',do_flush=.True.) 721 722 call xmpi_abort() 723 call leave_new('PERS',exit_status=0,print_config=.false.) 724 725 end if 726 727 end subroutine testkgrid