TABLE OF CONTENTS


ABINIT/m_inkpts [ Modules ]

[ Top ] [ 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