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-2024 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 wtk = one 216 call intagm(dprarr,intarr,jdtset,marr,nkpt,string(1:lenstr),'wtk',tread,'DPR') 217 if(tread==1) wtk(1:nkpt)=dprarr(1:nkpt) 218 219 wtksum=sum(wtk(:)) 220 write(msg,'(a,i0,a,f12.6)')' inkpts: Sum of ',nkpt,' k point weights is',wtksum 221 call wrtout(std_out,msg) 222 223 if (wtksum < tol6) then 224 write(msg, '(3a)' )& 225 'This sum is too close to zero. ',ch10,& 226 'Action: correct the array wtk in the input file.' 227 ABI_ERROR(msg) 228 end if 229 if (abs(wtksum - one) > tol6) then 230 if (occopt==2) then 231 write(msg, '(a,1p,e18.8,a,a,a)' )& 232 'wtksum= ',wtksum,' /= 1.0 means wts do not add to 1 , while occopt=2.',ch10,& 233 'Action: correct the array wtk in input file.' 234 ABI_ERROR(msg) 235 else 236 write(msg,'(a,i0,a)')' With present occopt= ',occopt,', renormalize it to one' 237 call wrtout(std_out,msg) 238 norm=one/wtksum 239 wtk(1:nkpt)=wtk(1:nkpt)*norm 240 end if 241 end if 242 end if 243 244 else if (kptopt == 0 .and. getkerange_filepath /= ABI_NOFILE) then 245 ! Initialize k-points from kerange_path file. 246 use_kerange = .True. 247 ABI_MALLOC(krange2ibz, (nkpt)) 248 if (my_rank == master) then 249 #ifdef HAVE_NETCDF 250 NCF_CHECK(nctk_open_read(ncid, getkerange_filepath, xmpi_comm_self)) 251 call hdr_ncread(hdr, ncid, fform) 252 ABI_CHECK(fform == fform_from_ext("KERANGE.nc"), sjoin("Error while reading:", getkerange_filepath, ", fform:", itoa(fform))) 253 ! TODO Add code for consistency check 254 !kptopt, nsym, occopt 255 !ABI_CHECK(nkpt == hdr%nkpt, "nkpt from kerange != nkpt") 256 NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "krange2ibz"), krange2ibz)) 257 NCF_CHECK(nf90_close(ncid)) 258 #else 259 ABI_ERROR("getkerange_filepath requires NETCDF support") 260 #endif 261 end if 262 263 call xmpi_bcast(krange2ibz, master, comm, ierr) 264 call hdr%bcast(master, my_rank, comm) 265 ! Hdr contains the kpts in the IBZ. Extract k-points in the pockets via krange2ibz. 266 nshiftk = hdr%nshiftk; nshiftk_orig = hdr%nshiftk_orig 267 istwfk = hdr%istwfk(krange2ibz(:)) 268 kptrlatt = hdr%kptrlatt; kptrlatt_orig = hdr%kptrlatt_orig 269 ABI_CHECK(isdiagmat(hdr%kptrlatt), "kptrlatt is not diagonal!") 270 ngkpt(1) = hdr%kptrlatt(1, 1); ngkpt(2) = hdr%kptrlatt(2, 2); ngkpt(3) = hdr%kptrlatt(3, 3) 271 kpt = hdr%kptns(:, krange2ibz(:)) !; kpthf(3,nkpthf) 272 shiftk(:,1:nshiftk) = hdr%shiftk; shiftk_orig(:, 1:nshiftk_orig) = hdr%shiftk_orig 273 wtk = hdr%wtk(krange2ibz(:)) 274 call hdr%free() 275 kptnrm = one 276 ABI_FREE(krange2ibz) 277 278 else if (kptopt < 0) then 279 ! Band structure calculation 280 nsegment=abs(kptopt) 281 282 if (iscf /= -2)then 283 write(msg, '(3a,i0,3a)' ) & 284 'For a negative kptopt, iscf must be -2,',ch10,& 285 'while it is found to be ',iscf,'.',ch10,& 286 'Action: change the value of iscf in your input file, or change kptopt.' 287 ABI_ERROR(msg) 288 end if 289 290 if(marr<3*nsegment+3)then 291 marr=3*nsegment+3 292 ABI_FREE(dprarr) 293 ABI_FREE(intarr) 294 ABI_MALLOC(dprarr,(marr)) 295 ABI_MALLOC(intarr,(marr)) 296 end if 297 298 ABI_MALLOC(kptbounds,(3,nsegment+1)) 299 ABI_MALLOC(ndivk,(nsegment)) 300 301 call intagm(dprarr,intarr,jdtset,marr,3*nsegment+3,string(1:lenstr),'kptbounds',tread,'DPR') 302 303 if(tread==1)then 304 kptbounds(:,:)=reshape( dprarr(1:3*nsegment+3), [3,nsegment+1]) 305 else 306 write(msg,'(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 ABI_ERROR(msg) 311 end if 312 313 call intagm(dprarr,intarr,jdtset,marr,nsegment,string(1:lenstr),'ndivk',tread,'INT') 314 if(tread==1)then 315 ndivk(1:nsegment)=intarr(1:nsegment) 316 ! The 1 stand for the first point 317 nkpt_computed=1+sum(ndivk(1:nsegment)) 318 319 ! ndivk and ndivsm are mutually exclusive 320 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ndivsm',tread,'INT') 321 if (tread == 1) then 322 ABI_ERROR("ndivk and ndivsm are mutually exclusive. Choose only one variable") 323 end if 324 325 else 326 ! Calculate ndivk such as the path is normalized 327 ! Note that if both ndivk and ndivsm are defined in in input file, only ndivk is used ! 328 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ndivsm',tread,'INT') 329 if(tread==1)then 330 ndiv_small=intarr(1) 331 call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol) 332 call mknormpath(nsegment+1,kptbounds,gmet,ndiv_small,ndivk,nkpt_computed) 333 else 334 write(msg,'(5a)') & 335 'When kptopt is negative, ndivsm or ndivk must be initialized ',ch10,& 336 'in the input file, which is not the case.',ch10,& 337 'Action: initialize ndivsm or ndivk in your input file, or change kptopt.' 338 ABI_ERROR(msg) 339 end if 340 end if 341 342 ! Check that the argument nkpt is coherent with nkpt_computed 343 if (nkpt/=0 .and. nkpt /= nkpt_computed) then 344 write(msg, '(a,i0,5a,i0,7a)' ) & 345 'The argument nkpt = ',nkpt,', does not match',ch10,& 346 'the number of k points generated by kptopt, ndivk, kptbounds,',ch10,& 347 'and the eventual symmetries, that is, nkpt= ',nkpt_computed,'.',ch10,& 348 'However, note that it might due to the user,',ch10,& 349 'if nkpt is explicitely defined in the input file.',ch10,& 350 'In this case, please check your input file.' 351 ABI_ERROR(msg) 352 end if 353 354 if (nkpt/=0) then 355 ! The array kpt has the right dimension and we can generate the k-path 356 call intagm(dprarr,intarr,jdtset,marr,3*nsegment+3,string(1:lenstr),'kptbounds',tread,'DPR') 357 if(tread==1)then 358 kptbounds(:,:)=reshape( dprarr(1:3*nsegment+3), [3,nsegment+1]) 359 else 360 write(msg, '(5a)') & 361 'When kptopt is negative, kptbounds must be initialized ',ch10,& 362 'in the input file, which is not the case.',ch10,& 363 'Action: initialize kptbounds in your input file, or change kptopt.' 364 ABI_ERROR(msg) 365 end if 366 367 ! First k point 368 jkpt=1 369 kpt(:,1)=kptbounds(:,1) 370 do ii=1,nsegment 371 dkpt=ndivk(ii) 372 do ikpt=1,dkpt 373 fraction=dble(ikpt)/dble(dkpt) 374 kpt(:,ikpt+jkpt)=fraction *kptbounds(:,ii+1)+(one-fraction)*kptbounds(:,ii) 375 end do 376 jkpt=jkpt+dkpt 377 end do 378 379 end if 380 381 kptnrm=one 382 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptnrm',tread,'DPR') 383 if(tread==1) kptnrm=dprarr(1) 384 385 ABI_FREE(kptbounds) 386 ABI_FREE(ndivk) 387 388 else if (kptopt>=1 .and. kptopt<=4) then 389 ! Read ngkpt 390 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'ngkpt',tread_ngkpt,'INT') 391 if(tread_ngkpt==1)then 392 ngkpt(1:3)=intarr(1:3) 393 do ii=1,3 394 if(ngkpt(ii)<1)then 395 write(msg,'(a,i0,3a,i0,3a)') & 396 'The input variable ngkpt(',ii,') must be strictly positive,',ch10,& 397 'while it is found to be ',ngkpt(ii),'.',ch10,& 398 'Action: change it in your input file, or change kptopt.' 399 ABI_ERROR(msg) 400 end if 401 end do 402 end if 403 404 call intagm(dprarr,intarr,jdtset,marr,9,string(1:lenstr),'kptrlatt',tread_kptrlatt,'INT') 405 if(tread_kptrlatt==1) kptrlatt(:,:)=reshape(intarr(1:9), [3,3]) 406 407 if(tread_ngkpt==1 .and. tread_kptrlatt==1)then 408 write(msg, '(5a)') & 409 'The input variables ngkpt and kptrlatt cannot both ',ch10,& 410 'be defined in the input file.',ch10,& 411 'Action: change one of ngkpt or kptrlatt in your input file.' 412 ABI_ERROR(msg) 413 else if(tread_ngkpt==1)then 414 kptrlatt(:,:)=0 415 kptrlatt(1,1)=ngkpt(1) 416 kptrlatt(2,2)=ngkpt(2) 417 kptrlatt(3,3)=ngkpt(3) 418 ! Save kptrlatt for reference. 419 kptrlatt_orig = kptrlatt 420 end if 421 422 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nshiftk',tread,'INT') 423 if(tread==1)nshiftk=intarr(1) 424 425 if (nshiftk < 1 .or. nshiftk > MAX_NSHIFTK )then 426 write(msg, '(a,i0,2a,i0,3a)' )& 427 'The only allowed values of nshiftk are between 1 and ',MAX_NSHIFTK,ch10,& 428 'while it is found to be',nshiftk,'.',ch10,& 429 'Action: change the value of nshiftk in your input file, or change kptopt.' 430 ABI_ERROR(msg) 431 end if 432 433 call intagm(dprarr,intarr,jdtset,marr,3*nshiftk,string(1:lenstr),'shiftk',tread,'DPR') 434 if(tread==1)then 435 shiftk(:,1:nshiftk)=reshape( dprarr(1:3*nshiftk), [3,nshiftk]) 436 ! Save input shifts as they will be changes in getkgrid. 437 nshiftk_orig = nshiftk 438 shiftk_orig(:,1:nshiftk) = shiftk(:,1:nshiftk) 439 else 440 if(nshiftk/=1)then 441 write(msg, '(3a,i0,2a)' )& 442 'When nshiftk is not equal to 1, shiftk must be defined in the input file.',ch10,& 443 'However, shiftk is not defined, while nshiftk=',nshiftk,ch10,& 444 'Action: change the value of nshiftk in your input file, or define shiftk.' 445 ABI_ERROR(msg) 446 end if 447 ! Default values used in indefo 448 nshiftk_orig = 1 449 shiftk_orig(:,1:nshiftk) = half 450 end if 451 452 prtkpt=0 453 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'prtkpt',tread,'INT') 454 if(tread==1)prtkpt=intarr(1) 455 456 if(sum(abs(kptrlatt(:,:)))==0)then 457 kptrlatt_orig = 0 458 do ii=1,3 459 kptrlatt_orig(ii,ii) = ngkpt(ii) 460 end do 461 ! The parameters of the k lattice are not known, compute kptrlatt, nshiftk, shiftk. 462 call testkgrid(bravais,iout,kptrlatt,kptrlen, msym,nshiftk,nsym,prtkpt,rprimd,shiftk,symafm,symrel,vacuum) 463 end if 464 465 ! TODO: Avoid call to getkgrid if eph 466 !eph_task = -1 467 !call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'eph_task',tread,'INT') 468 !if(tread==1) eph_task=intarr(1) 469 470 fockdownsampling(:)=1 471 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'fockdownsampling',tread,'INT') 472 if(tread==1)fockdownsampling=intarr(1:3) 473 474 call getkgrid(chksymbreak,0,iscf,kpt,kptopt,kptrlatt,kptrlen,& 475 msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,& 476 shiftk,symafm,symrel,vacuum,wtk,nkpthf=nkpthf,kpthf=kpthf,downsampling=fockdownsampling) 477 478 kptnrm=one 479 480 else 481 write(msg,'(3a,i0,3a)' ) & 482 'The only values of kptopt allowed are smaller than 4.',ch10,& 483 'The input value of kptopt is: ',kptopt,'.',ch10,& 484 'Action: change kptopt in your input file.' 485 ABI_ERROR(msg) 486 end if 487 488 if (kptnrm < tol10) then 489 write(msg, '(5a)' )& 490 'The input variable kptnrm is lower than 1.0d-10,',ch10,& 491 'while it must be a positive, non-zero number. ',ch10,& 492 'Action: correct the kptnrm in the input file.' 493 ABI_ERROR(msg) 494 end if 495 496 ! The k point number has been computed, and, if nkpt/=0, also the list of k points. 497 ! Also nkpthf has been computed, and, if nkpt/=0, also the list kpthf. 498 ! Now, determine istwfk, and eventually shift the k points by the value of qptn. 499 if (nkpt /= 0) then 500 istwfk(1:nkpt)=0 501 call intagm(dprarr,intarr,jdtset,marr,nkpt,string(1:lenstr),'istwfk',tread,'INT') 502 if(tread==1) istwfk(1:nkpt)=intarr(1:nkpt) 503 504 ! Impose istwfk=1 for RF calculations or NSCF calculation with kpts from kerange. 505 if (response == 1 .or. use_kerange) istwfk(1:nkpt)=1 506 507 ! Also impose istwfk=1 for spinor calculations 508 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nspinor',tread,'INT') 509 if(tread/=0 .and. intarr(1)/=1)istwfk(1:nkpt)=1 510 511 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'pawspnorb',tread,'INT') 512 if(tread/=0 .and. intarr(1)/=0)istwfk(1:nkpt)=1 513 514 do ikpt=1,nkpt 515 if(istwfk(ikpt)==0)then 516 kpoint=kpt(:,ikpt)/kptnrm; if (nqpt/=0.and.response==0) kpoint=kpoint+qptn 517 istwfk(ikpt) = set_istwfk(kpoint) 518 end if 519 if (present(impose_istwf_1)) then 520 if (impose_istwf_1==1) then 521 istwfk(ikpt)=1 522 else if (impose_istwf_1==2.and.any(kpt(:,ikpt)>tol10)) then 523 istwfk(ikpt)=1 524 end if 525 end if 526 end do 527 end if 528 529 ! If nkpt was to be computed, transfer it from nkpt_computed 530 if (nkpt == 0) nkpt = nkpt_computed 531 532 ABI_FREE(intarr) 533 ABI_FREE(dprarr) 534 535 call timab(192,2,tsec) 536 537 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
568 subroutine inqpt(chksymbreak,iout,jdtset,lenstr,msym,natom,qptn,wtqc,rprimd,spinat,string,typat,vacuum,xred,qptrlatt) 569 570 !Arguments ------------------------------------ 571 !scalars 572 integer,intent(in) :: chksymbreak,iout,jdtset,lenstr,msym,natom 573 real(dp),intent(inout) :: wtqc 574 character(len=*),intent(in) :: string 575 !arrays 576 integer,intent(in) :: typat(natom),vacuum(3) 577 real(dp),intent(out) :: qptn(3) 578 integer,intent(inout) :: qptrlatt(3,3) !vz_i 579 real(dp),intent(in) :: rprimd(3,3) 580 real(dp),intent(in) :: spinat(3,natom) 581 real(dp),intent(in) :: xred(3,natom) 582 583 !Local variables------------------------------- 584 !scalars 585 integer :: ii,iqpt,iscf_fake,marr,nptsym,nqpt_max,nqpt_computed,nshiftq,nsym_new,qptopt 586 integer :: tread,tread_q_sum,tread_qptrlatt,tread_ngqpt,use_inversion 587 real(dp) :: qptnrm,qptrlen,tolsym,ucvol 588 character(len=500) :: msg 589 !arrays 590 integer :: bravais(11), ngqpt(3) 591 integer, allocatable :: symafm_new(:), ptsymrel(:,:,:),symrel_new(:,:,:), intarr(:) 592 real(dp) :: gmet(3,3),gprimd(3,3),qpt(3),rmet(3,3),shiftq(3,MAX_NSHIFTK) 593 real(dp),allocatable :: qpts(:,:),tnons_new(:,:),wtq(:), dprarr(:) 594 595 ! ************************************************************************* 596 597 ! Compute the maximum size of arrays intarr and dprarr (nshiftq is MAX_NSHIFTK at maximum) 598 marr=630 599 ABI_MALLOC(intarr,(marr)) 600 ABI_MALLOC(dprarr,(marr)) 601 tread_q_sum=0 602 603 ! Find the method to generate the q-points 604 qptopt=0 605 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'qptopt',tread,'INT') 606 tread_q_sum=tread_q_sum+tread 607 if(tread==1)qptopt=intarr(1) 608 609 if(qptopt==0)then 610 ! Read qpt and qptnrm 611 qpt=zero 612 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'qpt',tread,'DPR') 613 tread_q_sum=tread_q_sum+tread 614 if(tread==1) qpt(1:3)=dprarr(1:3) 615 616 qptnrm=1.0_dp 617 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'qptnrm',tread,'DPR') 618 tread_q_sum=tread_q_sum+tread 619 620 if(tread==1) qptnrm=dprarr(1) 621 if(qptnrm<tol10)then 622 write(msg, '(5a)' )& 623 'The input variable qptnrm is lower than 1.0d-10,',ch10,& 624 'while it must be a positive, non-zero number. ',ch10,& 625 'Action: correct the qptnrm in the input file.' 626 ABI_ERROR(msg) 627 end if 628 629 qptn(:)=qpt(:)/qptnrm 630 631 ! DBSP: one could want ot define wtq in order to reproduce what is obtained 632 ! with ngqpt but without having to do initialize the qgrid (extremly slow in case of large grid > 50x50x50 633 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'wtq',tread,'DPR') 634 tread_q_sum=tread_q_sum+tread 635 if(tread==1) wtqc=dprarr(1) 636 637 else if (qptopt>=1 .and. qptopt<=4) then 638 ngqpt(:)=0 639 call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'ngqpt',tread_ngqpt,'INT') 640 tread_q_sum=tread_q_sum+tread_ngqpt 641 642 if(tread_ngqpt==1)then 643 ngqpt(1:3)=intarr(1:3) 644 do ii=1,3 645 if(ngqpt(ii)<1)then 646 write(msg, '(a,i0,a,a,a,i0,a,a,a)' ) & 647 'The input variable ngqpt(',ii,') must be strictly positive,',ch10,& 648 'while it is found to be',ngqpt(ii),'.',ch10,& 649 'Action: change it in your input file, or change qptopt.' 650 ABI_ERROR(msg) 651 end if 652 end do 653 end if 654 655 call intagm(dprarr,intarr,jdtset,marr,9,string(1:lenstr),'qptrlatt',tread_qptrlatt,'INT') 656 if(tread_qptrlatt==1) qptrlatt(:,:)=reshape(intarr(1:9), (/3,3/) ) 657 tread_q_sum=tread_q_sum+tread_qptrlatt 658 659 if(tread_ngqpt==1 .and. tread_qptrlatt==1)then 660 write(msg, '(5a)' ) & 661 'The input variables ngqpt and qptrlatt cannot both ',ch10,& 662 'be defined in the input file.',ch10,& 663 'Action: change one of ngqpt or qptrlatt in your input file.' 664 ABI_ERROR(msg) 665 else if(tread_ngqpt==1)then 666 qptrlatt(:,:)=0 667 qptrlatt(1,1)=ngqpt(1) 668 qptrlatt(2,2)=ngqpt(2) 669 qptrlatt(3,3)=ngqpt(3) 670 end if 671 672 nshiftq=1 673 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nshiftq',tread,'INT') 674 tread_q_sum=tread_q_sum+tread 675 if(tread==1)nshiftq=intarr(1) 676 677 if (nshiftq<1 .or. nshiftq>MAX_NSHIFTK) then 678 write(msg, '(a,i0,2a,i0,3a)' )& 679 'The only allowed values of nshiftq are between 1 and,',MAX_NSHIFTK,ch10,& 680 'while it is found to be',nshiftq,'.',ch10,& 681 'Action: change the value of nshiftq in your input file, or change qptopt.' 682 ABI_ERROR(msg) 683 end if 684 685 shiftq=zero 686 call intagm(dprarr,intarr,jdtset,marr,3*nshiftq,string(1:lenstr),'shiftq',tread,'DPR') 687 tread_q_sum=tread_q_sum+tread 688 689 if(tread==1)then 690 shiftq(:,1:nshiftq)=reshape( dprarr(1:3*nshiftq), (/3,nshiftq/) ) 691 else 692 if(nshiftq/=1)then 693 write(msg, '(3a,i0,2a)' )& 694 'When nshiftq is not equal to 1, shiftq must be defined in the input file.',ch10,& 695 'However, shiftq is not defined, while nshiftq=',nshiftq,ch10,& 696 'Action: change the value of nshiftq in your input file, or define shiftq.' 697 ABI_ERROR(msg) 698 end if 699 end if 700 701 !DEBUG 702 !write(std_out,'(a)')' m_inkpts%inqpt : before symlatt ' 703 !call flush(std_out) 704 !ENDDEBUG 705 706 ! Re-generate symmetry operations from the lattice and atomic coordinates 707 ! This is a fundamental difference with respect to the k point generation. 708 tolsym=tol8 709 ABI_MALLOC(ptsymrel,(3,3,msym)) 710 ABI_MALLOC(symafm_new,(msym)) 711 ABI_MALLOC(symrel_new,(3,3,msym)) 712 ABI_MALLOC(tnons_new,(3,msym)) 713 call symlatt(bravais,dev_null,msym,nptsym,ptsymrel,rprimd,tolsym) 714 use_inversion=1 715 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 716 717 !DEBUG 718 !write(std_out,'(a)')' m_inkpts%inqpt : before symfind ' 719 !call flush(std_out) 720 !ENDDEBUG 721 722 call symfind(gprimd,msym,natom,nptsym,1,nsym_new,0,& 723 ptsymrel,spinat,symafm_new,symrel_new,tnons_new,tolsym,typat,use_inversion,xred) 724 725 !DEBUG 726 !write(std_out,'(a)')' m_inkpts%inqpt : after symfind ' 727 !call flush(std_out) 728 !ENDDEBUG 729 730 ! Prepare to compute the q-point grid in the ZB or IZB 731 iscf_fake=0 ! Do not need the weights 732 733 ! Compute the maximum number of q points 734 nqpt_max=0 735 ABI_MALLOC(qpts,(3,nqpt_max)) 736 ABI_MALLOC(wtq,(nqpt_max)) 737 call getkgrid(chksymbreak,0,iscf_fake,qpts,qptopt,qptrlatt,qptrlen,& 738 msym,nqpt_max,nqpt_computed,nshiftq,nsym_new,rprimd,& 739 shiftq,symafm_new,symrel_new,vacuum,wtq) 740 741 nqpt_max=nqpt_computed 742 ABI_FREE(qpts) 743 ABI_FREE(wtq) 744 745 ! Find the index of the q point within the set of q points that will be generated 746 iqpt=0 747 call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'iqpt',tread,'INT') 748 tread_q_sum=tread_q_sum+tread 749 if(tread==1)iqpt=intarr(1) 750 751 ! Checks that iqpt is among the computed q points 752 if(iqpt<0)then 753 write(msg, '(a,i0,3a)' )& 754 'The input variable iqpt,',iqpt,' is negative, while it should be 0 or positive.',ch10,& 755 'Action: correct iqpt in the input file.' 756 ABI_ERROR(msg) 757 end if 758 759 if (iqpt > nqpt_computed) then 760 write(msg, '(a,i0,3a,i0,7a)' )& 761 'The input variable iqpt,',iqpt,' is bigger than the computed number of q-points in the grid,',ch10,& 762 'which is ',nqpt_max,'.',ch10,& 763 'The latter has been computed from the input variables qptrlatt, ngqpt, nshiftq,',ch10,& 764 'shiftq, as well as qptopt, the symmetries of the lattice, and spinat.',ch10,& 765 'Action: correct iqpt in the input file, or correct the computed q-point grid.' 766 ABI_ERROR(msg) 767 end if 768 769 ! Compute the q-point grid in the BZ or the IBZ 770 ABI_MALLOC(qpts,(3,nqpt_max)) 771 ABI_MALLOC(wtq,(nqpt_max)) 772 773 call getkgrid(chksymbreak,iout,iscf_fake,qpts,qptopt,qptrlatt,qptrlen,& 774 msym,nqpt_max,nqpt_computed,nshiftq,nsym_new,rprimd,& 775 shiftq,symafm_new,symrel_new,vacuum,wtq) 776 777 ! Transfer to qptn, and deallocate 778 qptn(:)=zero 779 if(iqpt/=0)then 780 qptn(:)=qpts(:,iqpt) 781 wtqc = wtq(iqpt) 782 end if 783 784 ABI_FREE(ptsymrel) 785 ABI_FREE(symafm_new) 786 ABI_FREE(symrel_new) 787 ABI_FREE(tnons_new) 788 ABI_FREE(qpts) 789 ABI_FREE(wtq) 790 791 else 792 write(msg, '(3a,i0,3a)' ) & 793 'The only values of qptopt allowed are smaller than 4.',ch10,& 794 'The input value of qptopt is',qptopt,'.',ch10,& 795 'Action: change qptopt in your input file.' 796 ABI_ERROR(msg) 797 end if 798 799 ! See issue #31 on gitlab. Not really a good idea. 800 !if(nqpt==0 .and. tread_q_sum/=0)then 801 ! write(msg, '(5a)' ) & 802 ! 'When nqpt is zero, the following input variables cannot be defined :',ch10, & 803 ! ' iqpt, ngqpt, nshiftq, qptopt, qpt, qptnrm, qptrlatt, shiftq, wtq . ',ch10, & 804 ! 'Action: change nqpt to 1, or un-define all the variables above.' 805 ! ABI_ERROR(msg) 806 !endif 807 808 ABI_FREE(intarr) 809 ABI_FREE(dprarr) 810 811 end subroutine inqpt