TABLE OF CONTENTS
ABINIT/m_inkpts [ Modules ]
NAME
m_inkpts
FUNCTION
Routines to initialize k-point and q-point sampling from input file.
COPYRIGHT
Copyright (C) 1998-2022 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_inkpts 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 use m_nctk 29 #ifdef HAVE_NETCDF 30 use netcdf 31 #endif 32 use m_hdr 33 34 use m_time, only : timab 35 use m_fstrings, only : sjoin, itoa 36 use m_numeric_tools, only : isdiagmat 37 use m_geometry, only : metric 38 use m_symfind, only : symfind, symlatt 39 use m_cgtools, only : set_istwfk 40 use m_parser, only : intagm 41 use m_kpts, only : getkgrid, testkgrid, mknormpath 42 43 implicit none 44 45 private
m_inkpts/inkpts [ Functions ]
[ Top ] [ m_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.
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 getkerange_filepath= Path of KERANGE.nc file used to initialize k-point sampling if kptopt == 0 and string != ABI_NOFILE 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 comm= MPI communicator
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,MAX_NSHIFTK)=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,MAX_NSHIFTK)=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!
SOURCE
125 subroutine inkpts(bravais,chksymbreak,fockdownsampling,iout,iscf,istwfk,jdtset,& 126 & kpt,kpthf,kptopt,kptnrm,kptrlatt_orig,kptrlatt,kptrlen,lenstr,msym, getkerange_filepath, & 127 & nkpt,nkpthf,nqpt,ngkpt,nshiftk,nshiftk_orig,shiftk_orig,nsym,& 128 & occopt,qptn,response,rprimd,shiftk,string,symafm,symrel,vacuum,wtk,comm,& 129 & impose_istwf_1) ! Optional argument 130 131 !Arguments ------------------------------------ 132 !scalars 133 integer,intent(in) :: chksymbreak,iout,iscf,jdtset,kptopt,lenstr,msym,nqpt,nsym,occopt 134 integer,intent(in) :: response, comm 135 integer,intent(in),optional :: impose_istwf_1 136 integer,intent(inout) :: nkpt,nkpthf 137 integer,intent(out) :: nshiftk,nshiftk_orig 138 integer,intent(out) :: fockdownsampling(3) 139 real(dp),intent(out) :: kptnrm,kptrlen 140 character(len=*),intent(in) :: string 141 character(len=*),intent(in) :: getkerange_filepath 142 !arrays 143 integer,intent(in) :: bravais(11),symafm(msym),symrel(3,3,msym),vacuum(3) 144 integer,intent(out) :: istwfk(nkpt),kptrlatt(3,3),kptrlatt_orig(3,3),ngkpt(3) 145 real(dp),intent(in) :: rprimd(3,3),qptn(3) 146 real(dp),intent(out) :: kpt(3,nkpt),kpthf(3,nkpthf),shiftk(3,MAX_NSHIFTK),wtk(nkpt),shiftk_orig(3,MAX_NSHIFTK) 147 148 !Local variables------------------------------- 149 !scalars 150 integer,parameter :: master = 0 151 integer :: dkpt,ii,ikpt,jkpt,marr,ndiv_small,nkpt_computed,my_rank,nprocs 152 integer :: nsegment,prtkpt,tread,tread_kptrlatt,tread_ngkpt, ncid, fform, ierr 153 logical :: use_kerange 154 real(dp) :: fraction,norm,ucvol,wtksum 155 character(len=500) :: msg 156 type(hdr_type) :: hdr 157 !arrays 158 integer,allocatable :: ndivk(:),intarr(:), krange2ibz(:) 159 real(dp) :: gmet(3,3),gprimd(3,3),kpoint(3),rmet(3,3),tsec(2) 160 real(dp),allocatable :: kptbounds(:,:),dprarr(:) 161 162 ! ************************************************************************* 163 164 call timab(192,1,tsec) 165 166 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 167 168 ! Compute the maximum size of arrays intarr and dprarr 169 marr = max(3*nkpt,3*MAX_NSHIFTK) 170 ABI_MALLOC(intarr,(marr)) 171 ABI_MALLOC(dprarr,(marr)) 172 173 ! Use zero to signal that these values have not been read. 174 ngkpt = 0 175 shiftk_orig = zero 176 kptrlatt_orig = 0; kptrlatt = 0 177 nshiftk_orig = 1; nshiftk = 1 178 use_kerange = .False. 179 180 !fockdownsampling(:)=1 181 !kptnrm = one 182 !kpthf = zero 183 184 ! MG: FIXME These values should be initialized because they are intent(out) 185 ! but several tests fails. So we keep this bug to avoid problems somewhere else 186 ! The initialization of the kpoints should be rewritten in a cleaner way 187 ! without all these side effects! 188 !shiftk = zero 189 ! Initializing these three variables is OK but we keep the bug to preserve the old behavior 190 !wtk = one 191 !kpt = zero 192 !istwfk = 1 193 194 ! Initialize kptrlen 195 kptrlen=30.0_dp 196 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptrlen',tread,'DPR') 197 if(tread==1)kptrlen=dprarr(1) 198 199 ! Initialize kpt, kptnrm and wtk according to kptopt. 200 if (kptopt == 0 .and. getkerange_filepath == ABI_NOFILE) then 201 ! For kptopt==0, one must have nkpt defined. 202 kpt(:,:)=zero 203 call intagm(dprarr,intarr,jdtset,marr,3*nkpt,string(1:lenstr),'kpt',tread,'DPR') 204 if(tread==1) kpt(:,:)=reshape( dprarr(1:3*nkpt), [3,nkpt]) 205 206 kptnrm=one 207 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptnrm',tread,'DPR') 208 if(tread==1) kptnrm=dprarr(1) 209 210 ! Only read wtk when iscf >0 or iscf=-1 or iscf=-3 or (iscf=-2 and response=1) 211 ! (this last option is for Zach Levine) 212 ! Normalize the k-point weights when occopt/=2 213 ! Check that k point weights add to 1 when occopt==2 214 if (iscf>0.or.iscf==-1.or.iscf==-3.or.(iscf==-2.and.response==1)) then 215 call intagm(dprarr,intarr,jdtset,marr,nkpt,string(1:lenstr),'wtk',tread,'DPR') 216 if(tread==1) wtk(1:nkpt)=dprarr(1:nkpt) 217 218 wtksum=sum(wtk(:)) 219 write(msg,'(a,i0,a,f12.6)')' inkpts: Sum of ',nkpt,' k point weights is',wtksum 220 call wrtout(std_out,msg) 221 222 if (wtksum < tol6) then 223 write(msg, '(3a)' )& 224 'This sum is too close to zero. ',ch10,& 225 'Action: correct the array wtk in the input file.' 226 ABI_ERROR(msg) 227 end if 228 if (abs(wtksum - one) > tol6) then 229 if (occopt==2) then 230 write(msg, '(a,1p,e18.8,a,a,a)' )& 231 'wtksum= ',wtksum,' /= 1.0 means wts do not add to 1 , while occopt=2.',ch10,& 232 'Action: correct the array wtk in input file.' 233 ABI_ERROR(msg) 234 else 235 write(msg,'(a,i0,a)')' With present occopt= ',occopt,', renormalize it to one' 236 call wrtout(std_out,msg) 237 norm=one/wtksum 238 wtk(1:nkpt)=wtk(1:nkpt)*norm 239 end if 240 end if 241 end if 242 243 else if (kptopt == 0 .and. getkerange_filepath /= ABI_NOFILE) then 244 ! Initialize k-points from kerange_path file. 245 use_kerange = .True. 246 ABI_MALLOC(krange2ibz, (nkpt)) 247 if (my_rank == master) then 248 #ifdef HAVE_NETCDF 249 NCF_CHECK(nctk_open_read(ncid, getkerange_filepath, xmpi_comm_self)) 250 call hdr_ncread(hdr, ncid, fform) 251 ABI_CHECK(fform == fform_from_ext("KERANGE.nc"), sjoin("Error while reading:", getkerange_filepath, ", fform:", itoa(fform))) 252 ! TODO Add code for consistency check 253 !kptopt, nsym, occopt 254 !ABI_CHECK(nkpt == hdr%nkpt, "nkpt from kerange != nkpt") 255 NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "krange2ibz"), krange2ibz)) 256 NCF_CHECK(nf90_close(ncid)) 257 #else 258 ABI_ERROR("getkerange_filepath requires NETCDF support") 259 #endif 260 end if 261 262 call xmpi_bcast(krange2ibz, master, comm, ierr) 263 call hdr%bcast(master, my_rank, comm) 264 ! Hdr contains the kpts in the IBZ. Extract k-points in the pockets via krange2ibz. 265 nshiftk = hdr%nshiftk; nshiftk_orig = hdr%nshiftk_orig 266 istwfk = hdr%istwfk(krange2ibz(:)) 267 kptrlatt = hdr%kptrlatt; kptrlatt_orig = hdr%kptrlatt_orig 268 ABI_CHECK(isdiagmat(hdr%kptrlatt), "kptrlatt is not diagonal!") 269 ngkpt(1) = hdr%kptrlatt(1, 1); ngkpt(2) = hdr%kptrlatt(2, 2); ngkpt(3) = hdr%kptrlatt(3, 3) 270 kpt = hdr%kptns(:, krange2ibz(:)) !; kpthf(3,nkpthf) 271 shiftk(:,1:nshiftk) = hdr%shiftk; shiftk_orig(:, 1:nshiftk_orig) = hdr%shiftk_orig 272 wtk = hdr%wtk(krange2ibz(:)) 273 call hdr%free() 274 kptnrm = one 275 ABI_FREE(krange2ibz) 276 277 else if (kptopt < 0) then 278 ! Band structure calculation 279 nsegment=abs(kptopt) 280 281 if (iscf /= -2)then 282 write(msg, '(3a,i0,3a)' ) & 283 'For a negative kptopt, iscf must be -2,',ch10,& 284 'while it is found to be ',iscf,'.',ch10,& 285 'Action: change the value of iscf in your input file, or change kptopt.' 286 ABI_ERROR(msg) 287 end if 288 289 if(marr<3*nsegment+3)then 290 marr=3*nsegment+3 291 ABI_FREE(dprarr) 292 ABI_FREE(intarr) 293 ABI_MALLOC(dprarr,(marr)) 294 ABI_MALLOC(intarr,(marr)) 295 end if 296 297 ABI_MALLOC(kptbounds,(3,nsegment+1)) 298 ABI_MALLOC(ndivk,(nsegment)) 299 300 call intagm(dprarr,intarr,jdtset,marr,3*nsegment+3,string(1:lenstr),'kptbounds',tread,'DPR') 301 302 if(tread==1)then 303 kptbounds(:,:)=reshape( dprarr(1:3*nsegment+3), [3,nsegment+1]) 304 else 305 write(msg,'(5a)') & 306 'When kptopt is negative, kptbounds must be initialized ',ch10,& 307 'in the input file, which is not the case.',ch10,& 308 'Action: initialize kptbounds in your input file, or change kptopt.' 309 ABI_ERROR(msg) 310 end if 311 312 call intagm(dprarr,intarr,jdtset,marr,nsegment,string(1:lenstr),'ndivk',tread,'INT') 313 if(tread==1)then 314 ndivk(1:nsegment)=intarr(1:nsegment) 315 ! The 1 stand for the first point 316 nkpt_computed=1+sum(ndivk(1:nsegment)) 317 318 ! ndivk and ndivsm are mutually exclusive 319 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ndivsm',tread,'INT') 320 if (tread == 1) then 321 ABI_ERROR("ndivk and ndivsm are mutually exclusive. Choose only one variable") 322 end if 323 324 else 325 ! Calculate ndivk such as the path is normalized 326 ! Note that if both ndivk and ndivsm are defined in in input file, only ndivk is used ! 327 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ndivsm',tread,'INT') 328 if(tread==1)then 329 ndiv_small=intarr(1) 330 call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol) 331 call mknormpath(nsegment+1,kptbounds,gmet,ndiv_small,ndivk,nkpt_computed) 332 else 333 write(msg,'(5a)') & 334 'When kptopt is negative, ndivsm or ndivk must be initialized ',ch10,& 335 'in the input file, which is not the case.',ch10,& 336 'Action: initialize ndivsm or ndivk in your input file, or change kptopt.' 337 ABI_ERROR(msg) 338 end if 339 end if 340 341 ! Check that the argument nkpt is coherent with nkpt_computed 342 if (nkpt/=0 .and. nkpt /= nkpt_computed) then 343 write(msg, '(a,i0,5a,i0,7a)' ) & 344 'The argument nkpt = ',nkpt,', does not match',ch10,& 345 'the number of k points generated by kptopt, ndivk, kptbounds,',ch10,& 346 'and the eventual symmetries, that is, nkpt= ',nkpt_computed,'.',ch10,& 347 'However, note that it might due to the user,',ch10,& 348 'if nkpt is explicitely defined in the input file.',ch10,& 349 'In this case, please check your input file.' 350 ABI_ERROR(msg) 351 end if 352 353 if (nkpt/=0) then 354 ! The array kpt has the right dimension and we can generate the k-path 355 call intagm(dprarr,intarr,jdtset,marr,3*nsegment+3,string(1:lenstr),'kptbounds',tread,'DPR') 356 if(tread==1)then 357 kptbounds(:,:)=reshape( dprarr(1:3*nsegment+3), [3,nsegment+1]) 358 else 359 write(msg, '(5a)') & 360 'When kptopt is negative, kptbounds must be initialized ',ch10,& 361 'in the input file, which is not the case.',ch10,& 362 'Action: initialize kptbounds in your input file, or change kptopt.' 363 ABI_ERROR(msg) 364 end if 365 366 ! First k point 367 jkpt=1 368 kpt(:,1)=kptbounds(:,1) 369 do ii=1,nsegment 370 dkpt=ndivk(ii) 371 do ikpt=1,dkpt 372 fraction=dble(ikpt)/dble(dkpt) 373 kpt(:,ikpt+jkpt)=fraction *kptbounds(:,ii+1)+(one-fraction)*kptbounds(:,ii) 374 end do 375 jkpt=jkpt+dkpt 376 end do 377 378 end if 379 380 kptnrm=one 381 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptnrm',tread,'DPR') 382 if(tread==1) kptnrm=dprarr(1) 383 384 ABI_FREE(kptbounds) 385 ABI_FREE(ndivk) 386 387 else if (kptopt>=1 .and. kptopt<=4) then 388 ! Read ngkpt 389 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'ngkpt',tread_ngkpt,'INT') 390 if(tread_ngkpt==1)then 391 ngkpt(1:3)=intarr(1:3) 392 do ii=1,3 393 if(ngkpt(ii)<1)then 394 write(msg,'(a,i0,3a,i0,3a)') & 395 'The input variable ngkpt(',ii,') must be strictly positive,',ch10,& 396 'while it is found to be ',ngkpt(ii),'.',ch10,& 397 'Action: change it in your input file, or change kptopt.' 398 ABI_ERROR(msg) 399 end if 400 end do 401 end if 402 403 call intagm(dprarr,intarr,jdtset,marr,9,string(1:lenstr),'kptrlatt',tread_kptrlatt,'INT') 404 if(tread_kptrlatt==1) kptrlatt(:,:)=reshape(intarr(1:9), [3,3]) 405 406 if(tread_ngkpt==1 .and. tread_kptrlatt==1)then 407 write(msg, '(5a)') & 408 'The input variables ngkpt and kptrlatt cannot both ',ch10,& 409 'be defined in the input file.',ch10,& 410 'Action: change one of ngkpt or kptrlatt in your input file.' 411 ABI_ERROR(msg) 412 else if(tread_ngkpt==1)then 413 kptrlatt(:,:)=0 414 kptrlatt(1,1)=ngkpt(1) 415 kptrlatt(2,2)=ngkpt(2) 416 kptrlatt(3,3)=ngkpt(3) 417 ! Save kptrlatt for reference. 418 kptrlatt_orig = kptrlatt 419 end if 420 421 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nshiftk',tread,'INT') 422 if(tread==1)nshiftk=intarr(1) 423 424 if (nshiftk < 1 .or. nshiftk > MAX_NSHIFTK )then 425 write(msg, '(a,i0,2a,i0,3a)' )& 426 'The only allowed values of nshiftk are between 1 and ',MAX_NSHIFTK,ch10,& 427 'while it is found to be',nshiftk,'.',ch10,& 428 'Action: change the value of nshiftk in your input file, or change kptopt.' 429 ABI_ERROR(msg) 430 end if 431 432 call intagm(dprarr,intarr,jdtset,marr,3*nshiftk,string(1:lenstr),'shiftk',tread,'DPR') 433 if(tread==1)then 434 shiftk(:,1:nshiftk)=reshape( dprarr(1:3*nshiftk), [3,nshiftk]) 435 ! Save input shifts as they will be changes in getkgrid. 436 nshiftk_orig = nshiftk 437 shiftk_orig(:,1:nshiftk) = shiftk(:,1:nshiftk) 438 else 439 if(nshiftk/=1)then 440 write(msg, '(3a,i0,2a)' )& 441 'When nshiftk is not equal to 1, shiftk must be defined in the input file.',ch10,& 442 'However, shiftk is not defined, while nshiftk=',nshiftk,ch10,& 443 'Action: change the value of nshiftk in your input file, or define shiftk.' 444 ABI_ERROR(msg) 445 end if 446 ! Default values used in indefo 447 nshiftk_orig = 1 448 shiftk_orig(:,1:nshiftk) = half 449 end if 450 451 prtkpt=0 452 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'prtkpt',tread,'INT') 453 if(tread==1)prtkpt=intarr(1) 454 455 if(sum(abs(kptrlatt(:,:)))==0)then 456 kptrlatt_orig = 0 457 do ii=1,3 458 kptrlatt_orig(ii,ii) = ngkpt(ii) 459 end do 460 ! The parameters of the k lattice are not known, compute kptrlatt, nshiftk, shiftk. 461 call testkgrid(bravais,iout,kptrlatt,kptrlen, msym,nshiftk,nsym,prtkpt,rprimd,shiftk,symafm,symrel,vacuum) 462 end if 463 464 ! TODO: Avoid call to getkgrid if eph 465 !eph_task = -1 466 !call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'eph_task',tread,'INT') 467 !if(tread==1) eph_task=intarr(1) 468 469 fockdownsampling(:)=1 470 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'fockdownsampling',tread,'INT') 471 if(tread==1)fockdownsampling=intarr(1:3) 472 473 call getkgrid(chksymbreak,0,iscf,kpt,kptopt,kptrlatt,kptrlen,& 474 msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,& 475 shiftk,symafm,symrel,vacuum,wtk,nkpthf=nkpthf,kpthf=kpthf,downsampling=fockdownsampling) 476 477 kptnrm=one 478 479 else 480 write(msg,'(3a,i0,3a)' ) & 481 'The only values of kptopt allowed are smaller than 4.',ch10,& 482 'The input value of kptopt is: ',kptopt,'.',ch10,& 483 'Action: change kptopt in your input file.' 484 ABI_ERROR(msg) 485 end if 486 487 if (kptnrm < tol10) then 488 write(msg, '(5a)' )& 489 'The input variable kptnrm is lower than 1.0d-10,',ch10,& 490 'while it must be a positive, non-zero number. ',ch10,& 491 'Action: correct the kptnrm in the input file.' 492 ABI_ERROR(msg) 493 end if 494 495 ! The k point number has been computed, and, if nkpt/=0, also the list of k points. 496 ! Also nkpthf has been computed, and, if nkpt/=0, also the list kpthf. 497 ! Now, determine istwfk, and eventually shift the k points by the value of qptn. 498 if (nkpt /= 0) then 499 istwfk(1:nkpt)=0 500 call intagm(dprarr,intarr,jdtset,marr,nkpt,string(1:lenstr),'istwfk',tread,'INT') 501 if(tread==1) istwfk(1:nkpt)=intarr(1:nkpt) 502 503 ! Impose istwfk=1 for RF calculations or NSCF calculation with kpts from kerange. 504 if (response == 1 .or. use_kerange) istwfk(1:nkpt)=1 505 506 ! Also impose istwfk=1 for spinor calculations 507 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nspinor',tread,'INT') 508 if(tread/=0 .and. intarr(1)/=1)istwfk(1:nkpt)=1 509 510 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'pawspnorb',tread,'INT') 511 if(tread/=0 .and. intarr(1)/=0)istwfk(1:nkpt)=1 512 513 do ikpt=1,nkpt 514 if(istwfk(ikpt)==0)then 515 kpoint=kpt(:,ikpt)/kptnrm; if (nqpt/=0.and.response==0) kpoint=kpoint+qptn 516 istwfk(ikpt) = set_istwfk(kpoint) 517 end if 518 if (present(impose_istwf_1)) then 519 if (impose_istwf_1==1) then 520 istwfk(ikpt)=1 521 else if (impose_istwf_1==2.and.any(kpt(:,ikpt)>tol10)) then 522 istwfk(ikpt)=1 523 end if 524 end if 525 end do 526 end if 527 528 ! If nkpt was to be computed, transfer it from nkpt_computed 529 if (nkpt == 0) nkpt = nkpt_computed 530 531 ABI_FREE(intarr) 532 ABI_FREE(dprarr) 533 534 call timab(192,2,tsec) 535 536 end subroutine inkpts
m_inkpts/inqpt [ Functions ]
[ Top ] [ m_inkpts ] [ Functions ]
NAME
inqpt
FUNCTION
Initialize the q point for one particular dataset, characterized by jdtset.
INPUTS
chksymbreak= if 1, will check whether the k point grid is symmetric, and stop if not. iout=unit number for echoed output jdtset=number of the dataset looked for lenstr=actual length of the string msym=default maximal number of symmetries natom=number of atoms rprimd(3,3)=dimensional real space primitive translations (bohr) spinat(3,1:natom)=spin-magnetization of the atoms string*(*)=character string containing all the input data. Initialized previously in instrng. typat(natom)=type for each atom vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum xred(3,natom,nimage) =reduced coordinates of atoms
OUTPUT
qptn(3)=reduced coordinates of eventual q point (normalisation is already included) kptrlatt(3,3)=q-point lattice specification (if kptopt/=0) wtqc=weigth of the eventual current q point
SOURCE
567 subroutine inqpt(chksymbreak,iout,jdtset,lenstr,msym,natom,qptn,wtqc,rprimd,spinat,string,typat,vacuum,xred,qptrlatt) 568 569 !Arguments ------------------------------------ 570 !scalars 571 integer,intent(in) :: chksymbreak,iout,jdtset,lenstr,msym,natom 572 real(dp),intent(inout) :: wtqc 573 character(len=*),intent(in) :: string 574 !arrays 575 integer,intent(in) :: typat(natom),vacuum(3) 576 real(dp),intent(out) :: qptn(3) 577 integer,intent(inout) :: qptrlatt(3,3) !vz_i 578 real(dp),intent(in) :: rprimd(3,3) 579 real(dp),intent(in) :: spinat(3,natom) 580 real(dp),intent(in) :: xred(3,natom) 581 582 !Local variables------------------------------- 583 !scalars 584 integer :: ii,iqpt,iscf_fake,marr,nptsym,nqpt_max,nqpt_computed,nshiftq,nsym_new,qptopt 585 integer :: tread,tread_q_sum,tread_qptrlatt,tread_ngqpt,use_inversion 586 real(dp) :: qptnrm,qptrlen,tolsym,ucvol 587 character(len=500) :: msg 588 !arrays 589 integer :: bravais(11), ngqpt(3) 590 integer, allocatable :: symafm_new(:), ptsymrel(:,:,:),symrel_new(:,:,:), intarr(:) 591 real(dp) :: gmet(3,3),gprimd(3,3),qpt(3),rmet(3,3),shiftq(3,MAX_NSHIFTK) 592 real(dp),allocatable :: qpts(:,:),tnons_new(:,:),wtq(:), dprarr(:) 593 594 ! ************************************************************************* 595 596 ! Compute the maximum size of arrays intarr and dprarr (nshiftq is MAX_NSHIFTK at maximum) 597 marr=630 598 ABI_MALLOC(intarr,(marr)) 599 ABI_MALLOC(dprarr,(marr)) 600 tread_q_sum=0 601 602 ! Find the method to generate the q-points 603 qptopt=0 604 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'qptopt',tread,'INT') 605 tread_q_sum=tread_q_sum+tread 606 if(tread==1)qptopt=intarr(1) 607 608 if(qptopt==0)then 609 ! Read qpt and qptnrm 610 qpt=zero 611 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'qpt',tread,'DPR') 612 tread_q_sum=tread_q_sum+tread 613 if(tread==1) qpt(1:3)=dprarr(1:3) 614 615 qptnrm=1.0_dp 616 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'qptnrm',tread,'DPR') 617 tread_q_sum=tread_q_sum+tread 618 619 if(tread==1) qptnrm=dprarr(1) 620 if(qptnrm<tol10)then 621 write(msg, '(5a)' )& 622 'The input variable qptnrm is lower than 1.0d-10,',ch10,& 623 'while it must be a positive, non-zero number. ',ch10,& 624 'Action: correct the qptnrm in the input file.' 625 ABI_ERROR(msg) 626 end if 627 628 qptn(:)=qpt(:)/qptnrm 629 630 ! DBSP: one could want ot define wtq in order to reproduce what is obtained 631 ! with ngqpt but without having to do initialize the qgrid (extremly slow in case of large grid > 50x50x50 632 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'wtq',tread,'DPR') 633 tread_q_sum=tread_q_sum+tread 634 if(tread==1) wtqc=dprarr(1) 635 636 else if (qptopt>=1 .and. qptopt<=4) then 637 ngqpt(:)=0 638 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'ngqpt',tread_ngqpt,'INT') 639 tread_q_sum=tread_q_sum+tread_ngqpt 640 641 if(tread_ngqpt==1)then 642 ngqpt(1:3)=intarr(1:3) 643 do ii=1,3 644 if(ngqpt(ii)<1)then 645 write(msg, '(a,i0,a,a,a,i0,a,a,a)' ) & 646 'The input variable ngqpt(',ii,') must be strictly positive,',ch10,& 647 'while it is found to be',ngqpt(ii),'.',ch10,& 648 'Action: change it in your input file, or change qptopt.' 649 ABI_ERROR(msg) 650 end if 651 end do 652 end if 653 654 call intagm(dprarr,intarr,jdtset,marr,9,string(1:lenstr),'qptrlatt',tread_qptrlatt,'INT') 655 if(tread_qptrlatt==1) qptrlatt(:,:)=reshape(intarr(1:9), (/3,3/) ) 656 tread_q_sum=tread_q_sum+tread_qptrlatt 657 658 if(tread_ngqpt==1 .and. tread_qptrlatt==1)then 659 write(msg, '(5a)' ) & 660 'The input variables ngqpt and qptrlatt cannot both ',ch10,& 661 'be defined in the input file.',ch10,& 662 'Action: change one of ngqpt or qptrlatt in your input file.' 663 ABI_ERROR(msg) 664 else if(tread_ngqpt==1)then 665 qptrlatt(:,:)=0 666 qptrlatt(1,1)=ngqpt(1) 667 qptrlatt(2,2)=ngqpt(2) 668 qptrlatt(3,3)=ngqpt(3) 669 end if 670 671 nshiftq=1 672 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nshiftq',tread,'INT') 673 tread_q_sum=tread_q_sum+tread 674 if(tread==1)nshiftq=intarr(1) 675 676 if (nshiftq<1 .or. nshiftq>MAX_NSHIFTK) then 677 write(msg, '(a,i0,2a,i0,3a)' )& 678 'The only allowed values of nshiftq are between 1 and,',MAX_NSHIFTK,ch10,& 679 'while it is found to be',nshiftq,'.',ch10,& 680 'Action: change the value of nshiftq in your input file, or change qptopt.' 681 ABI_ERROR(msg) 682 end if 683 684 shiftq=zero 685 call intagm(dprarr,intarr,jdtset,marr,3*nshiftq,string(1:lenstr),'shiftq',tread,'DPR') 686 tread_q_sum=tread_q_sum+tread 687 688 if(tread==1)then 689 shiftq(:,1:nshiftq)=reshape( dprarr(1:3*nshiftq), (/3,nshiftq/) ) 690 else 691 if(nshiftq/=1)then 692 write(msg, '(3a,i0,2a)' )& 693 'When nshiftq is not equal to 1, shiftq must be defined in the input file.',ch10,& 694 'However, shiftq is not defined, while nshiftq=',nshiftq,ch10,& 695 'Action: change the value of nshiftq in your input file, or define shiftq.' 696 ABI_ERROR(msg) 697 end if 698 end if 699 700 ! Re-generate symmetry operations from the lattice and atomic coordinates 701 ! This is a fundamental difference with respect to the k point generation. 702 tolsym=tol8 703 ABI_MALLOC(ptsymrel,(3,3,msym)) 704 ABI_MALLOC(symafm_new,(msym)) 705 ABI_MALLOC(symrel_new,(3,3,msym)) 706 ABI_MALLOC(tnons_new,(3,msym)) 707 call symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym) 708 use_inversion=1 709 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 710 call symfind(0,(/zero,zero,zero/),gprimd,0,msym,natom,0,nptsym,nsym_new,0,0,& 711 ptsymrel,spinat,symafm_new,symrel_new,tnons_new,tolsym,typat,use_inversion,xred) 712 713 ! Prepare to compute the q-point grid in the ZB or IZB 714 iscf_fake=0 ! Do not need the weights 715 716 ! Compute the maximum number of q points 717 nqpt_max=0 718 ABI_MALLOC(qpts,(3,nqpt_max)) 719 ABI_MALLOC(wtq,(nqpt_max)) 720 call getkgrid(chksymbreak,0,iscf_fake,qpts,qptopt,qptrlatt,qptrlen,& 721 msym,nqpt_max,nqpt_computed,nshiftq,nsym_new,rprimd,& 722 shiftq,symafm_new,symrel_new,vacuum,wtq) 723 724 nqpt_max=nqpt_computed 725 ABI_FREE(qpts) 726 ABI_FREE(wtq) 727 728 ! Find the index of the q point within the set of q points that will be generated 729 iqpt=0 730 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'iqpt',tread,'INT') 731 tread_q_sum=tread_q_sum+tread 732 if(tread==1)iqpt=intarr(1) 733 734 ! Checks that iqpt is among the computed q points 735 if(iqpt<0)then 736 write(msg, '(a,i0,3a)' )& 737 'The input variable iqpt,',iqpt,' is negative, while it should be 0 or positive.',ch10,& 738 'Action: correct iqpt in the input file.' 739 ABI_ERROR(msg) 740 end if 741 742 if (iqpt > nqpt_computed) then 743 write(msg, '(a,i0,3a,i0,7a)' )& 744 'The input variable iqpt,',iqpt,' is bigger than the computed number of q-points in the grid,',ch10,& 745 'which is ',nqpt_max,'.',ch10,& 746 'The latter has been computed from the input variables qptrlatt, ngqpt, nshiftq,',ch10,& 747 'shiftq, as well as qptopt, the symmetries of the lattice, and spinat.',ch10,& 748 'Action: correct iqpt in the input file, or correct the computed q-point grid.' 749 ABI_ERROR(msg) 750 end if 751 752 ! Compute the q-point grid in the BZ or the IBZ 753 ABI_MALLOC(qpts,(3,nqpt_max)) 754 ABI_MALLOC(wtq,(nqpt_max)) 755 756 call getkgrid(chksymbreak,iout,iscf_fake,qpts,qptopt,qptrlatt,qptrlen,& 757 msym,nqpt_max,nqpt_computed,nshiftq,nsym_new,rprimd,& 758 shiftq,symafm_new,symrel_new,vacuum,wtq) 759 760 ! Transfer to qptn, and deallocate 761 qptn(:)=zero 762 if(iqpt/=0)then 763 qptn(:)=qpts(:,iqpt) 764 wtqc = wtq(iqpt) 765 end if 766 767 ABI_FREE(ptsymrel) 768 ABI_FREE(symafm_new) 769 ABI_FREE(symrel_new) 770 ABI_FREE(tnons_new) 771 ABI_FREE(qpts) 772 ABI_FREE(wtq) 773 774 else 775 write(msg, '(3a,i0,3a)' ) & 776 'The only values of qptopt allowed are smaller than 4.',ch10,& 777 'The input value of qptopt is',qptopt,'.',ch10,& 778 'Action: change qptopt in your input file.' 779 ABI_ERROR(msg) 780 end if 781 782 ! See issue #31 on gitlab. Not really a good idea. 783 !if(nqpt==0 .and. tread_q_sum/=0)then 784 ! write(msg, '(5a)' ) & 785 ! 'When nqpt is zero, the following input variables cannot be defined :',ch10, & 786 ! ' iqpt, ngqpt, nshiftq, qptopt, qpt, qptnrm, qptrlatt, shiftq, wtq . ',ch10, & 787 ! 'Action: change nqpt to 1, or un-define all the variables above.' 788 ! ABI_ERROR(msg) 789 !endif 790 791 ABI_FREE(intarr) 792 ABI_FREE(dprarr) 793 794 end subroutine inqpt