TABLE OF CONTENTS
ABINIT/inkpts [ Functions ]
NAME
inkpts
FUNCTION
Initialize k points (list of k points, weights, storage) for one particular dataset, characterized by jdtset. Note that nkpt (and nkpthf) can be computed by calling this routine with input value of nkpt=0, provided kptopt/=0.
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)
chksymbreak= if 1, will check whether the k point grid is symmetric, and stop if not.
impose_istwf_1= (optional argument):
0: no restriction on istwfk
1: impose istwfk=1 for all k points
2: impose istwfk=1 for all k points non equal to zero
iout=unit number for echoed output
iscf= ( <= 0 =>non-SCF), >0 => SCF)
jdtset=number of the dataset looked for
lenstr=actual length of the string
kptopt=option for the generation of k points
msym=default maximal number of symmetries
nqpt=number of q points (0 or 1)
nsym=number of symetries
occopt=option for occupation numbers
qptn(3)=reduced coordinates of eventual q point shift (already normalized).
response=0 if GS case, =1 if RF case.
rprimd(3,3)=dimensional real space primitive translations (bohr)
string*(*)=character string containing all the input data.
Initialized previously in instrng.
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
fockdownsampling(3)=echo of input variable fockdownsampling(3) kptnrm=normalisation of k points kptrlatt_orig(3,3)=Original value of kptrlatt as specified in the input file (if kptopt/=0) kptrlatt(3,3)=k-point lattice specification (if kptopt/=0) kptrlen=length of the smallest real space supercell vector nshiftk_orig=Original number of k-point shifts (0 if not read) nshiftk=actual number of k-point shifts in shiftk (if kptopt/=0) shiftk(3,210)=shift vectors for k point generation (if kptopt/=0) If nkpt/=0 the following arrays are also output : istwfk(nkpt)=option parameters that describes the storage of wfs kpt(3,nkpt)=reduced coordinates of k points. kpthf(3,nkpthf)=reduced coordinates of k points for Fock operator. wtk(nkpt)=weight assigned to each k point. ngkpt(3)=Number of divisions along the three reduced directions (0 signals that this variable has not been used. shiftk_orig(3,210)=Original shifts read from the input file (0 signals that this variable has not been read).
SIDE EFFECTS
Input/output: nkpt=number of k points if non-zero at input, is only an input variable if zero at input, its actual value will be computed nkpthf=number of k points for Fock operator, computed if nkpt=0 at input
NOTES
Warning: this routine can be called with nkpt=0 (in which case it returns the true value of nkpt), which can lead to strange bugs in the debugging procedure, if one tries to print wtk or istwfk, in this case !
PARENTS
invars1,invars2
CHILDREN
getkgrid,intagm,metric,mknormpath,testkgrid,timab,wrtout
SOURCE
88 #if defined HAVE_CONFIG_H 89 #include "config.h" 90 #endif 91 92 #include "abi_common.h" 93 94 subroutine inkpts(bravais,chksymbreak,fockdownsampling,iout,iscf,istwfk,jdtset,& 95 & kpt,kpthf,kptopt,kptnrm,kptrlatt_orig,kptrlatt,kptrlen,lenstr,msym,& 96 & nkpt,nkpthf,nqpt,ngkpt,nshiftk,nshiftk_orig,shiftk_orig,nsym,& 97 & occopt,qptn,response,rprimd,shiftk,string,symafm,symrel,vacuum,wtk,& 98 & impose_istwf_1) ! Optional argument 99 100 use defs_basis 101 use m_profiling_abi 102 use m_errors 103 104 use m_geometry, only : metric 105 use m_cgtools, only : set_istwfk 106 use m_parser, only : intagm 107 108 !This section has been created automatically by the script Abilint (TD). 109 !Do not modify the following lines by hand. 110 #undef ABI_FUNC 111 #define ABI_FUNC 'inkpts' 112 use interfaces_14_hidewrite 113 use interfaces_18_timing 114 use interfaces_32_util 115 use interfaces_56_recipspace 116 !End of the abilint section 117 118 implicit none 119 120 !Arguments ------------------------------------ 121 !scalars 122 integer,intent(in) :: chksymbreak,iout,iscf,jdtset,kptopt,lenstr,msym,nqpt,nsym,occopt 123 integer,intent(in) :: response 124 integer,intent(in),optional :: impose_istwf_1 125 integer,intent(inout) :: nkpt,nkpthf 126 integer,intent(out) :: nshiftk,nshiftk_orig 127 integer,intent(out) :: fockdownsampling(3) 128 real(dp),intent(out) :: kptnrm,kptrlen 129 character(len=*),intent(in) :: string 130 !arrays 131 integer,intent(in) :: bravais(11),symafm(msym),symrel(3,3,msym),vacuum(3) 132 integer,intent(out) :: istwfk(nkpt),kptrlatt(3,3),kptrlatt_orig(3,3),ngkpt(3) 133 real(dp),intent(in) :: rprimd(3,3),qptn(3) 134 real(dp),intent(out) :: kpt(3,nkpt),kpthf(3,nkpthf),shiftk(3,210),wtk(nkpt),shiftk_orig(3,210) 135 136 !Local variables------------------------------- 137 !scalars 138 integer :: dkpt,ii,ikpt,jkpt,marr,ndiv_small,nkpt_computed 139 integer :: nsegment,prtkpt,tread,tread_kptrlatt,tread_ngkpt 140 real(dp) :: fraction,norm,ucvol,wtksum 141 character(len=500) :: message 142 !arrays 143 integer,allocatable :: ndivk(:),intarr(:) 144 real(dp) :: gmet(3,3),gprimd(3,3),kpoint(3),rmet(3,3),tsec(2) 145 real(dp),allocatable :: kptbounds(:,:),dprarr(:) 146 147 ! ************************************************************************* 148 149 call timab(192,1,tsec) 150 151 !Compute the maximum size of arrays intarr and dprarr 152 marr=max(3*nkpt,3*210) 153 ABI_ALLOCATE(intarr,(marr)) 154 ABI_ALLOCATE(dprarr,(marr)) 155 156 !Use zero to signal that these values have not been read. 157 ngkpt = 0 158 shiftk_orig = zero 159 kptrlatt_orig = 0; kptrlatt = 0 160 nshiftk_orig = 1; nshiftk = 1 161 162 ! MG: FIXME These values should be initialized because they are intent(out) 163 ! but several tests fails. So we keep this bug to avoid problems somewhere else 164 ! The initialization of the kpoints should be rewritten in a cleaner way 165 ! without all these side effects! 166 !shiftk = zero 167 ! Initializing these three variables is OK but we keep the bug to preserve the old behavior 168 !wtk = one 169 !kpt = zero 170 !istwfk = 1 171 172 !Initialize kptrlen 173 kptrlen=30.0_dp 174 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptrlen',tread,'DPR') 175 if(tread==1)kptrlen=dprarr(1) 176 177 !Initialize kpt, kptnrm and wtk according to kptopt. 178 !For kptopt==0, one must have nkpt defined. 179 if(kptopt==0)then 180 181 kpt(:,:)=zero 182 call intagm(dprarr,intarr,jdtset,marr,3*nkpt,string(1:lenstr),'kpt',tread,'DPR') 183 if(tread==1) kpt(:,:)=reshape( dprarr(1:3*nkpt), [3,nkpt]) 184 185 kptnrm=one 186 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptnrm',tread,'DPR') 187 if(tread==1) kptnrm=dprarr(1) 188 189 ! Only read wtk when iscf >0 or iscf=-1 or iscf=-3 or (iscf=-2 and response=1) 190 ! (this last option is for Zach Levine) 191 ! Normalize the k-point weights when occopt/=2 192 ! Check that k point weights add to 1 when occopt==2 193 if (iscf>0.or.iscf==-1.or.iscf==-3.or.(iscf==-2.and.response==1)) then 194 195 call intagm(dprarr,intarr,jdtset,marr,nkpt,string(1:lenstr),'wtk',tread,'DPR') 196 if(tread==1) wtk(1:nkpt)=dprarr(1:nkpt) 197 198 wtksum=sum(wtk(:)) 199 write(message,'(a,i0,a,f12.6)')' inkpts: Sum of ',nkpt,' k point weights is',wtksum 200 call wrtout(std_out,message,'COLL') 201 202 if (wtksum<1.d-6) then 203 write(message, '(6a)' )& 204 & 'This sum is too close to zero. ',ch10,& 205 & 'Action: correct the array wtk in the input file.' 206 MSG_ERROR(message) 207 end if 208 if (abs(wtksum - one)>1.d-06) then 209 if(occopt==2)then 210 write(message, '(a,1p,e18.8,a,a,a)' )& 211 & 'wtksum=',wtksum,' /=1.0 means wts do not add to 1 , while occopt=2.',ch10,& 212 & 'Action: correct the array wtk in input file.' 213 MSG_ERROR(message) 214 else 215 write(message,'(a,i0,a)')' With present occopt= ',occopt,', renormalize it to one' 216 call wrtout(std_out,message,'COLL') 217 norm=one/wtksum 218 wtk(1:nkpt)=wtk(1:nkpt)*norm 219 end if 220 end if 221 end if 222 223 else if (kptopt<0) then 224 225 ! Band structure calculation 226 nsegment=abs(kptopt) 227 228 if(iscf/=-2)then 229 write(message, '(3a,i0,3a)' ) & 230 & 'For a negative kptopt, iscf must be -2,',ch10,& 231 & 'while it is found to be ',iscf,'.',ch10,& 232 & 'Action: change the value of iscf in your input file, or change kptopt.' 233 MSG_ERROR(message) 234 end if 235 236 if(marr<3*nsegment+3)then 237 marr=3*nsegment+3 238 ABI_DEALLOCATE(dprarr) 239 ABI_DEALLOCATE(intarr) 240 ABI_ALLOCATE(dprarr,(marr)) 241 ABI_ALLOCATE(intarr,(marr)) 242 end if 243 244 ABI_ALLOCATE(kptbounds,(3,nsegment+1)) 245 ABI_ALLOCATE(ndivk,(nsegment)) 246 247 call intagm(dprarr,intarr,jdtset,marr,3*nsegment+3,string(1:lenstr),'kptbounds',tread,'DPR') 248 249 if(tread==1)then 250 kptbounds(:,:)=reshape( dprarr(1:3*nsegment+3), [3,nsegment+1]) 251 else 252 write(message,'(5a)') & 253 & 'When kptopt is negative, kptbounds must be initialized ',ch10,& 254 & 'in the input file, which is not the case.',ch10,& 255 & 'Action: initialize kptbounds in your input file, or change kptopt.' 256 MSG_ERROR(message) 257 end if 258 259 call intagm(dprarr,intarr,jdtset,marr,nsegment,string(1:lenstr),'ndivk',tread,'INT') 260 if(tread==1)then 261 ndivk(1:nsegment)=intarr(1:nsegment) 262 ! The 1 stand for the first point 263 nkpt_computed=1+sum(ndivk(1:nsegment)) 264 265 ! ndivk and ndivsm are mutually exclusive 266 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ndivsm',tread,'INT') 267 if (tread == 1) then 268 MSG_ERROR("ndivk and ndivsm are mutually exclusive. Choose only one variable") 269 end if 270 271 else 272 ! Calculate ndivk such as the path is normalized 273 ! Note that if both ndivk and ndivsm are defined in in input file, only ndivk is used ! 274 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ndivsm',tread,'INT') 275 if(tread==1)then 276 ndiv_small=intarr(1) 277 call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol) 278 call mknormpath(nsegment+1,kptbounds,gmet,ndiv_small,ndivk,nkpt_computed) 279 else 280 write(message,'(5a)') & 281 & 'When kptopt is negative, ndivsm or ndivk must be initialized ',ch10,& 282 & 'in the input file, which is not the case.',ch10,& 283 & 'Action : initialize ndivsm or ndivk in your input file, or change kptopt.' 284 MSG_ERROR(message) 285 end if 286 end if 287 288 ! Check that the argument nkpt is coherent with nkpt_computed, 289 if(nkpt/=0 .and. nkpt/=nkpt_computed)then 290 write(message, '(a,i0,5a,i0,7a)' ) & 291 & 'The argument nkpt= ',nkpt,', does not match',ch10,& 292 & 'the number of k points generated by kptopt, ndivk, kptbounds,',ch10,& 293 & 'and the eventual symmetries, that is, nkpt= ',nkpt_computed,'.',ch10,& 294 & 'However, note that it might due to the user,',ch10,& 295 & 'if nkpt is explicitely defined in the input file.',ch10,& 296 & 'In this case, please check your input file.' 297 MSG_ERROR(message) 298 end if 299 300 if (nkpt/=0) then 301 ! the array kpt has the right dimension and we can generate the k-path 302 call intagm(dprarr,intarr,jdtset,marr,3*nsegment+3,string(1:lenstr),'kptbounds',tread,'DPR') 303 if(tread==1)then 304 kptbounds(:,:)=reshape( dprarr(1:3*nsegment+3), [3,nsegment+1]) 305 else 306 write(message, '(5a)') & 307 & 'When kptopt is negative, kptbounds must be initialized ',ch10,& 308 & 'in the input file, which is not the case.',ch10,& 309 & 'Action: initialize kptbounds in your input file, or change kptopt.' 310 MSG_ERROR(message) 311 end if 312 313 ! First k point 314 jkpt=1 315 kpt(:,1)=kptbounds(:,1) 316 do ii=1,nsegment 317 dkpt=ndivk(ii) 318 do ikpt=1,dkpt 319 fraction=dble(ikpt)/dble(dkpt) 320 kpt(:,ikpt+jkpt)=fraction *kptbounds(:,ii+1)+(one-fraction)*kptbounds(:,ii) 321 end do 322 jkpt=jkpt+dkpt 323 end do 324 325 end if 326 327 kptnrm=one 328 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptnrm',tread,'DPR') 329 if(tread==1) kptnrm=dprarr(1) 330 331 ABI_DEALLOCATE(kptbounds) 332 ABI_DEALLOCATE(ndivk) 333 334 else if (kptopt>=1 .and. kptopt<=4) then 335 ! Read ngkpt 336 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'ngkpt',tread_ngkpt,'INT') 337 if(tread_ngkpt==1)then 338 ngkpt(1:3)=intarr(1:3) 339 do ii=1,3 340 if(ngkpt(ii)<1)then 341 write(message,'(a,i0,3a,i0,3a)') & 342 & 'The input variable ngkpt(',ii,') must be strictly positive,',ch10,& 343 & 'while it is found to be',ngkpt(ii),'.',ch10,& 344 & 'Action: change it in your input file, or change kptopt.' 345 MSG_ERROR(message) 346 end if 347 end do 348 end if 349 350 call intagm(dprarr,intarr,jdtset,marr,9,string(1:lenstr),'kptrlatt',tread_kptrlatt,'INT') 351 if(tread_kptrlatt==1) kptrlatt(:,:)=reshape(intarr(1:9), [3,3]) 352 353 if(tread_ngkpt==1 .and. tread_kptrlatt==1)then 354 write(message, '(5a)' ) & 355 & 'The input variables ngkpt and kptrlatt cannot both ',ch10,& 356 & 'be defined in the input file.',ch10,& 357 & 'Action : change one of ngkpt or kptrlatt in your input file.' 358 MSG_ERROR(message) 359 else if(tread_ngkpt==1)then 360 kptrlatt(:,:)=0 361 kptrlatt(1,1)=ngkpt(1) 362 kptrlatt(2,2)=ngkpt(2) 363 kptrlatt(3,3)=ngkpt(3) 364 ! Save kptrlatt for reference. 365 kptrlatt_orig = kptrlatt 366 end if 367 368 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nshiftk',tread,'INT') 369 if(tread==1)nshiftk=intarr(1) 370 371 if(nshiftk<1 .or. nshiftk>210 )then 372 write(message, '(3a,i0,3a)' )& 373 & 'The only allowed values of nshiftk are between 1 and 210,',ch10,& 374 & 'while it is found to be',nshiftk,'.',ch10,& 375 & 'Action: change the value of nshiftk in your input file, or change kptopt.' 376 MSG_ERROR(message) 377 end if 378 379 call intagm(dprarr,intarr,jdtset,marr,3*nshiftk,string(1:lenstr),'shiftk',tread,'DPR') 380 if(tread==1)then 381 shiftk(:,1:nshiftk)=reshape( dprarr(1:3*nshiftk), [3,nshiftk]) 382 ! Save input shifts as they will be changes in getkgrid. 383 nshiftk_orig = nshiftk 384 shiftk_orig(:,1:nshiftk) = shiftk(:,1:nshiftk) 385 else 386 if(nshiftk/=1)then 387 write(message, '(3a,i0,2a)' )& 388 & 'When nshiftk is not equal to 1, shiftk must be defined in the input file.',ch10,& 389 & 'However, shiftk is not defined, while nshiftk=',nshiftk,ch10,& 390 & 'Action: change the value of nshiftk in your input file, or define shiftk.' 391 MSG_ERROR(message) 392 end if 393 ! Default values used in indefo 394 nshiftk_orig = 1 395 shiftk_orig(:,1:nshiftk) = half 396 end if 397 398 prtkpt=0 399 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'prtkpt',tread,'INT') 400 if(tread==1)prtkpt=intarr(1) 401 402 if(sum(abs(kptrlatt(:,:)))==0)then 403 kptrlatt_orig = 0 404 do ii=1,3 405 kptrlatt_orig(ii,ii) = ngkpt(ii) 406 end do 407 408 ! The parameters of the k lattice are not known, compute kptrlatt, nshiftk, shiftk. 409 call testkgrid(bravais,iout,kptrlatt,kptrlen,& 410 & msym,nshiftk,nsym,prtkpt,rprimd,shiftk,symafm,symrel,vacuum) 411 412 end if 413 414 fockdownsampling(:)=1 415 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'fockdownsampling',tread,'INT') 416 if(tread==1)fockdownsampling=intarr(1:3) 417 418 call getkgrid(chksymbreak,0,iscf,kpt,kptopt,kptrlatt,kptrlen,& 419 & msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,& 420 & shiftk,symafm,symrel,vacuum,wtk,nkpthf=nkpthf,kpthf=kpthf,downsampling=fockdownsampling) 421 422 kptnrm=one 423 424 else 425 426 write(message, '(3a,i0,3a)' ) & 427 & 'The only values of kptopt allowed are smaller than 4.',ch10,& 428 & 'The input value of kptopt is',kptopt,'.',ch10,& 429 & 'Action: change kptopt in your input file.' 430 MSG_ERROR(message) 431 end if 432 433 if(kptnrm<tol10)then 434 write(message, '(5a)' )& 435 & 'The input variable kptnrm is lower than 1.0d-10,',ch10,& 436 & 'while it must be a positive, non-zero number. ',ch10,& 437 & 'Action: correct the kptnrm in the input file.' 438 MSG_ERROR(message) 439 end if 440 441 !The k point number has been computed, and, if nkpt/=0, also the list of k points. 442 !Also nkpthf has been computed, and, if nkpt/=0, also the list kpthf. 443 444 !Now, determine istwfk, and eventually shift the k points by the value of qptn. 445 446 if(nkpt/=0)then 447 448 istwfk(1:nkpt)=0 449 call intagm(dprarr,intarr,jdtset,marr,nkpt,string(1:lenstr),'istwfk',tread,'INT') 450 if(tread==1) istwfk(1:nkpt)=intarr(1:nkpt) 451 452 if(response==1)istwfk(1:nkpt)=1 ! Impose istwfk=1 for RF calculations 453 454 ! Also impose istwfk=1 for spinor calculations 455 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nspinor',tread,'INT') 456 if(tread/=0 .and. intarr(1)/=1)istwfk(1:nkpt)=1 457 458 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'pawspnorb',tread,'INT') 459 if(tread/=0 .and. intarr(1)/=0)istwfk(1:nkpt)=1 460 461 do ikpt=1,nkpt 462 if(istwfk(ikpt)==0)then 463 kpoint=kpt(:,ikpt)/kptnrm; if (nqpt/=0.and.response==0) kpoint=kpoint+qptn 464 istwfk(ikpt) = set_istwfk(kpoint) 465 end if 466 if (present(impose_istwf_1)) then 467 if (impose_istwf_1==1) then 468 istwfk(ikpt)=1 469 else if (impose_istwf_1==2.and.any(kpt(:,ikpt)>tol10)) then 470 istwfk(ikpt)=1 471 end if 472 end if 473 end do 474 end if 475 476 !If nkpt was to be computed, transfer it from nkpt_computed 477 if(nkpt==0)nkpt=nkpt_computed 478 479 ABI_DEALLOCATE(intarr) 480 ABI_DEALLOCATE(dprarr) 481 482 call timab(192,2,tsec) 483 484 end subroutine inkpts