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-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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_inkpts
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32 
33  use m_time,      only : timab
34  use m_geometry,  only : metric
35  use m_symfind,   only : symfind, symlatt
36  use m_cgtools,   only : set_istwfk
37  use m_parser,    only : intagm
38  use m_kpts,      only : getkgrid, testkgrid, mknormpath
39 
40  implicit none
41 
42  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
 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

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

SIDE EFFECTS

PARENTS

      invars1

CHILDREN

      getkgrid,intagm,metric,symfind,symlatt

SOURCE

549 subroutine inqpt(chksymbreak,iout,jdtset,lenstr,msym,natom,qptn,wtqc,rprimd,spinat,string,typat,vacuum,xred,qptrlatt)
550 
551 
552 !This section has been created automatically by the script Abilint (TD).
553 !Do not modify the following lines by hand.
554 #undef ABI_FUNC
555 #define ABI_FUNC 'inqpt'
556 !End of the abilint section
557 
558  implicit none
559 
560 !Arguments ------------------------------------
561 !scalars
562  integer,intent(in)   :: chksymbreak,iout,jdtset,lenstr,msym,natom
563  real(dp),intent(inout) :: wtqc
564  character(len=*),intent(in) :: string
565 !arrays
566  integer,intent(in) :: typat(natom),vacuum(3)
567  real(dp),intent(out) :: qptn(3)
568  integer,intent(inout) :: qptrlatt(3,3) !vz_i
569  real(dp),intent(in) :: rprimd(3,3)
570  real(dp),intent(in) :: spinat(3,natom)
571  real(dp),intent(in) :: xred(3,natom)
572 
573 !Local variables-------------------------------
574 !scalars
575  integer :: ii,iqpt,iscf_fake,marr,nptsym,nqpt,nqpt_computed,nshiftq,nsym_new,qptopt
576  integer :: tread,tread_qptrlatt,tread_ngqpt,use_inversion
577  real(dp) :: qptnrm,qptrlen,tolsym,ucvol
578  character(len=500) :: message
579 !arrays
580  integer :: bravais(11)
581  integer :: ngqpt(3)
582  integer, allocatable :: symafm_new(:)
583  integer, allocatable :: ptsymrel(:,:,:),symrel_new(:,:,:)
584  integer,allocatable :: intarr(:)
585  real(dp) :: gmet(3,3),gprimd(3,3),qpt(3),rmet(3,3),shiftq(3,210)
586  real(dp),allocatable :: qpts(:,:),tnons_new(:,:),wtq(:)
587  real(dp),allocatable :: dprarr(:)
588 
589 ! *************************************************************************
590 
591 !Compute the maximum size of arrays intarr and dprarr (nshiftq is 210 at maximum)
592  marr=630
593  ABI_ALLOCATE(intarr,(marr))
594  ABI_ALLOCATE(dprarr,(marr))
595 
596 !Find the method to generate the k points
597  qptopt=0
598  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'qptopt',tread,'INT')
599  if(tread==1)qptopt=intarr(1)
600 
601  if(qptopt==0)then
602 
603 !  Read qpt and qptnrm
604    qpt=zero
605    call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'qpt',tread,'DPR')
606    if(tread==1) qpt(1:3)=dprarr(1:3)
607 
608    qptnrm=1.0_dp
609    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'qptnrm',tread,'DPR')
610 
611    if(tread==1) qptnrm=dprarr(1)
612    if(qptnrm<tol10)then
613      write(message, '(5a)' )&
614 &     'The input variable qptnrm is lower than 1.0d-10,',ch10,&
615 &     'while it must be a positive, non-zero number.   ',ch10,&
616 &     'Action: correct the qptnrm in the input file.'
617      MSG_ERROR(message)
618    end if
619 
620    qptn(:)=qpt(:)/qptnrm
621 
622 !  DBSP: one could want ot define wtq in order to reproduce what is obtained
623 !  with ngqpt but without having to do initialize the qgrid (extremly slow in
624 !  case of large grid > 50x50x50
625    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'wtq',tread,'DPR')
626    if(tread==1) wtqc=dprarr(1)
627 
628  else if(qptopt>=1 .and. qptopt<=4)then
629 
630    ngqpt(:)=0
631    call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'ngqpt',tread_ngqpt,'INT')
632 
633    if(tread_ngqpt==1)then
634      ngqpt(1:3)=intarr(1:3)
635      do ii=1,3
636        if(ngqpt(ii)<1)then
637          write(message,  '(a,i0,a,a,a,i0,a,a,a)' ) &
638 &         'The input variable ngqpt(',ii,') must be strictly positive,',ch10,&
639 &         'while it is found to be',ngqpt(ii),'.',ch10,&
640 &         'Action: change it in your input file, or change qptopt.'
641          MSG_ERROR(message)
642        end if
643      end do
644    end if
645 
646    call intagm(dprarr,intarr,jdtset,marr,9,string(1:lenstr),'qptrlatt',tread_qptrlatt,'INT')
647    if(tread_qptrlatt==1) qptrlatt(:,:)=reshape(intarr(1:9), (/3,3/) )
648 
649    if(tread_ngqpt==1 .and. tread_qptrlatt==1)then
650      write(message,  '(a,a,a,a,a)' ) &
651 &     'The input variables ngqpt and qptrlatt cannot both ',ch10,&
652 &     'be defined in the input file.',ch10,&
653 &     'Action: change one of ngqpt or qptrlatt in your input file.'
654      MSG_ERROR(message)
655    else if(tread_ngqpt==1)then
656      qptrlatt(:,:)=0
657      qptrlatt(1,1)=ngqpt(1)
658      qptrlatt(2,2)=ngqpt(2)
659      qptrlatt(3,3)=ngqpt(3)
660    end if
661 
662    nshiftq=1
663    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nshiftq',tread,'INT')
664    if(tread==1)nshiftq=intarr(1)
665 
666    if(nshiftq<1 .or. nshiftq>210 )then
667      write(message,  '(3a,i4,a,a,a)' )&
668 &     'The only allowed values of nshiftq are between 1 and 210,',ch10,&
669 &     'while it is found to be',nshiftq,'.',ch10,&
670 &     'Action: change the value of nshiftq in your input file, or change qptopt.'
671      MSG_ERROR(message)
672    end if
673 
674    shiftq=half
675    call intagm(dprarr,intarr,jdtset,marr,3*nshiftq,string(1:lenstr),'shiftq',tread,'DPR')
676 
677    if(tread==1)then
678      shiftq(:,1:nshiftq)=reshape( dprarr(1:3*nshiftq), (/3,nshiftq/) )
679    else
680      if(nshiftq/=1)then
681        write(message,  '(3a,i4,a,a)' )&
682 &       'When nshiftq is not equal to 1, shiftq must be defined in the input file.',ch10,&
683 &       'However, shiftq is not defined, while nshiftq=',nshiftq,ch10,&
684 &       'Action: change the value of nshiftq in your input file, or define shiftq.'
685        MSG_ERROR(message)
686      end if
687    end if
688 
689 !  Re-generate symmetry operations from the lattice and atomic coordinates
690 !  This is a fundamental difference with respect to the k point generation.
691    tolsym=tol8
692    ABI_ALLOCATE(ptsymrel,(3,3,msym))
693    ABI_ALLOCATE(symafm_new,(msym))
694    ABI_ALLOCATE(symrel_new,(3,3,msym))
695    ABI_ALLOCATE(tnons_new,(3,msym))
696    call symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym)
697    use_inversion=1
698    call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
699    call symfind(0,(/zero,zero,zero/),gprimd,0,msym,natom,0,nptsym,nsym_new,0,0,&
700 &   ptsymrel,spinat,symafm_new,symrel_new,tnons_new,tolsym,typat,use_inversion,xred)
701 
702 !  Prepare to compute the q-point grid in the ZB or IZB
703    iscf_fake=0 ! Do not need the weights
704 
705 !  Compute the maximum number of q points
706    nqpt=0
707    ABI_ALLOCATE(qpts,(3,nqpt))
708    ABI_ALLOCATE(wtq,(nqpt))
709    call getkgrid(chksymbreak,0,iscf_fake,qpts,qptopt,qptrlatt,qptrlen,&
710 &   msym,nqpt,nqpt_computed,nshiftq,nsym_new,rprimd,&
711 &   shiftq,symafm_new,symrel_new,vacuum,wtq)
712    nqpt=nqpt_computed
713    ABI_DEALLOCATE(qpts)
714    ABI_DEALLOCATE(wtq)
715 
716 !  Find the index of the q point within the set of q points that will be generated
717    iqpt=0
718    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'iqpt',tread,'INT')
719    if(tread==1)iqpt=intarr(1)
720 
721 !  Checks that iqpt is among the computed q points
722    if(iqpt<0)then
723      write(message, '(a,i5,3a)' )&
724 &     'The input variable iqpt,',iqpt,' is negative, while it should be 0 or positive.',ch10,&
725 &     'Action: correct iqpt in the input file.'
726      MSG_ERROR(message)
727    end if
728 
729    if(iqpt>nqpt_computed)then
730      write(message, '(a,i5,3a,i5,7a)' )&
731 &     'The input variable iqpt,',iqpt,' is bigger than the computed number of q-points in the grid,',ch10,&
732 &     'which is ',nqpt,'.',ch10,&
733 &     'The latter has been computed from the input variables qptrlatt, ngqpt, nshiftq,',ch10,&
734 &     'shiftq, as well as qptopt, the symmetries of the lattice, and spinat.',ch10,&
735 &     'Action: correct iqpt in the input file, or correct the computed q-point grid.'
736      MSG_ERROR(message)
737    end if
738 
739 !  Compute the q-point grid in the BZ or the IBZ
740    ABI_ALLOCATE(qpts,(3,nqpt))
741    ABI_ALLOCATE(wtq,(nqpt))
742 
743    call getkgrid(chksymbreak,iout,iscf_fake,qpts,qptopt,qptrlatt,qptrlen,&
744 &   msym,nqpt,nqpt_computed,nshiftq,nsym_new,rprimd,&
745 &   shiftq,symafm_new,symrel_new,vacuum,wtq)
746 
747 !  Transfer to qptn, and deallocate
748    qptn(:)=zero
749    if(iqpt/=0)then
750      qptn(:)=qpts(:,iqpt)
751      wtqc = wtq(iqpt)
752    end if
753 
754    ABI_DEALLOCATE(ptsymrel)
755    ABI_DEALLOCATE(symafm_new)
756    ABI_DEALLOCATE(symrel_new)
757    ABI_DEALLOCATE(tnons_new)
758    ABI_DEALLOCATE(qpts)
759    ABI_DEALLOCATE(wtq)
760 
761  else
762 
763    write(message,  '(a,a,a,i4,a,a,a)' ) &
764 &   'The only values of qptopt allowed are smaller than 4.',ch10,&
765 &   'The input value of qptopt is',qptopt,'.',ch10,&
766 &   'Action: change qptopt in your input file.'
767    MSG_ERROR(message)
768  end if
769 
770  ABI_DEALLOCATE(intarr)
771  ABI_DEALLOCATE(dprarr)
772 
773 end subroutine inqpt