TABLE OF CONTENTS


ABINIT/invars1 [ Functions ]

[ Top ] [ Functions ]

NAME

 invars1

FUNCTION

 Initialize the dimensions needed to allocate the input arrays
 for one dataset characterized by jdtset, by
 taking from string the necessary data.
 Perform some preliminary checks and echo these dimensions.

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, AR, MKV)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt.

INPUTS

  iout=unit number of output file
  jdtset=number of the dataset looked for
  lenstr=actual length of string
  msym=default maximal number of symmetries
  npsp1= number of pseudopotential files
  zionpsp(npsp1)= valence charge over all psps

OUTPUT

  mband_upper=estimation of the maximum number of bands for any k-point

SIDE EFFECTS

 Input/Output (the default value is given in the calling routine)
  dtset=<type datafiles_type>contains all input variables,
   some of which are initialized here, while other were already
   initialized, while some others will still be initialized later.
   The list of records of dtset initialized in the present routine is:
   acell_orig,densty,iatfix,kptopt,kptrlatt,
   mkmem,mkqmem,mk1mem,natsph,natvshift,nconeq,nkpt,nkptgw,nkpthf,
   nqptdm,nshiftk,nucdipmom,nzchempot,optdriver,
   rprim_orig,rprimd_orig,shiftk,
   spgroup,spinat,typat,vel_orig,vel_cell_orig,xred_orig
  bravais(11)=characteristics of Bravais lattice (see symlatt.F90)
  symafm(1:msym)=(anti)ferromagnetic part of symmetry operations
  symrel(3,3,1:msym)=symmetry operations in real space in terms
   of primitive translations
  tnons(3,1:msym)=nonsymmorphic translations for symmetry operations
  string*(*)=string of characters containing all input variables and data

NOTES

 Must set up the geometry of the system, needed to compute
 k point grids in an automatic fashion.
 Treat separately mband_upper, since
 fband, charge and zionpsp must be known for being able to initialize it.

 Defaults are provided in the calling routine.
 Defaults are also provided here for the following variables :
 mband_upper, occopt, fband, charge
 They should be kept consistent with defaults of the same variables
 provided to the invars routines.

PARENTS

      invars1m

CHILDREN

      atomdata_from_znucl,chkint_ge,ingeo,inkpts,inqpt,intagm,inupper
      invacuum,mkrdim,wrtout

SOURCE

 69 #if defined HAVE_CONFIG_H
 70 #include "config.h"
 71 #endif
 72 
 73 #include "abi_common.h"
 74 
 75 subroutine invars1(bravais,dtset,iout,jdtset,lenstr,mband_upper,msym,npsp1,&
 76 & string,symafm,symrel,tnons,zionpsp)
 77 
 78  use defs_basis
 79  use defs_abitypes
 80  use m_errors
 81  use m_profiling_abi
 82  use m_xmpi
 83  use m_atomdata
 84 
 85  use m_fstrings, only : inupper
 86  use m_geometry, only : mkrdim
 87  use m_parser,   only : intagm, chkint_ge
 88 
 89 !This section has been created automatically by the script Abilint (TD).
 90 !Do not modify the following lines by hand.
 91 #undef ABI_FUNC
 92 #define ABI_FUNC 'invars1'
 93  use interfaces_14_hidewrite
 94  use interfaces_57_iovars, except_this_one => invars1
 95 !End of the abilint section
 96 
 97  implicit none
 98 
 99 !Arguments ------------------------------------
100 !scalars
101  integer,intent(in) :: iout,jdtset,lenstr,msym,npsp1
102  integer,intent(out) :: mband_upper
103  character(len=*),intent(inout) :: string
104  type(dataset_type),intent(inout) :: dtset
105 !arrays
106  integer,intent(inout) :: bravais(11),symafm(msym),symrel(3,3,msym)
107  real(dp),intent(inout) :: tnons(3,msym)
108  real(dp),intent(in) :: zionpsp(npsp1)
109 
110 !Local variables-------------------------------
111  character :: blank=' ',string1
112 !scalars
113  integer :: chksymbreak,found,ierr,iatom,ii,ikpt,iimage,index_blank,index_lower
114  integer :: index_typsymb,index_upper,ipsp,iscf,intimage,itypat,leave,marr
115  integer :: natom,nkpt,nkpthf,npsp,npspalch
116  integer :: nqpt,nspinor,nsppol,ntypat,ntypalch,ntyppure,occopt,response
117  integer :: rfddk,rfelfd,rfphon,rfstrs,rfuser,rf2_dkdk,rf2_dkde,rfmagn
118  integer :: tfband,tnband,tread,tread_alt
119  real(dp) :: charge,fband,kptnrm,kptrlen,zelect,zval
120  character(len=2) :: string2,symbol
121  character(len=500) :: message
122  type(atomdata_t) :: atom
123 !arrays
124  integer :: cond_values(4),vacuum(3)
125  integer,allocatable :: iatfix(:,:),intarr(:),istwfk(:),nband(:),typat(:)
126  real(dp) :: acell(3),rprim(3,3)
127 !real(dp) :: field(3)
128  real(dp),allocatable :: amu(:),dprarr(:),kpt(:,:),kpthf(:,:),mixalch(:,:),nucdipmom(:,:)
129  real(dp),allocatable :: ratsph(:),reaalloc(:),spinat(:,:)
130  real(dp),allocatable :: vel(:,:),vel_cell(:,:),wtk(:),xred(:,:),znucl(:)
131  character(len=32) :: cond_string(4)
132 
133 
134 !************************************************************************
135 
136 !Some initialisations
137  ierr=0
138  cond_string(1:4)=' '
139  cond_values(1:4)=(/0,0,0,0/)
140 
141 !Read parameters
142  marr=dtset%npsp;if (dtset%npsp<3) marr=3
143  marr=max(marr,dtset%nimage)
144  ABI_ALLOCATE(intarr,(marr))
145  ABI_ALLOCATE(dprarr,(marr))
146 
147 !---------------------------------------------------------------------------
148 
149  rfddk=0; rfelfd=0; rfphon=0; rfmagn=0; rfstrs=0; rfuser=0; rf2_dkdk=0; rf2_dkde=0
150  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'rfddk',tread,'INT')
151  if(tread==1) rfddk=intarr(1)
152  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'rfelfd',tread,'INT')
153  if(tread==1) rfelfd=intarr(1)
154  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'rfmagn',tread,'INT')
155  if(tread==1) rfmagn=intarr(1)
156  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'rfphon',tread,'INT')
157  if(tread==1) rfphon=intarr(1)
158  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'rfstrs',tread,'INT')
159  if(tread==1) rfstrs=intarr(1)
160  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'rfuser',tread,'INT')
161  if(tread==1) rfuser=intarr(1)
162  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'rf2_dkdk',tread,'INT')
163  if(tread==1) rf2_dkdk=intarr(1)
164  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'rf2_dkde',tread,'INT')
165  if(tread==1) rf2_dkde=intarr(1)
166 
167  response=0
168  if(rfddk/=0.or.rf2_dkdk/=0.or.rf2_dkde/=0.or.rfelfd/=0.or.rfphon/=0.or.rfstrs/=0.or.rfuser/=0.or.rfmagn/=0)response=1
169 
170  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'optdriver',tread,'INT')
171  if (tread==1) then
172    dtset%optdriver=intarr(1)
173  else
174 !  If optdriver was not read, while response=1, set optdriver to 1
175    if(response==1)dtset%optdriver=1
176  end if
177 
178 !---------------------------------------------------------------------------
179 !For now, waiting express parallelisation for recursion
180  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'tfkinfunc',tread,'INT')
181  if(tread==1) dtset%tfkinfunc=intarr(1)
182 
183 !---------------------------------------------------------------------------
184 ! wvl_bigdft_comp, done here since default values of nline, nwfshist and iscf
185 ! depend on its value (see indefo)
186  if(dtset%usewvl==1) then
187    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'wvl_bigdft_comp',tread,'INT')
188    if(tread==1) dtset%wvl_bigdft_comp=intarr(1)
189  end if
190 
191 !---------------------------------------------------------------------------
192 
193  natom=dtset%natom
194  npsp=dtset%npsp
195  ntypat=dtset%ntypat
196 
197 !No default value for znucl
198  call intagm(dprarr,intarr,jdtset,marr,dtset%npsp,string(1:lenstr),'znucl',tread,'DPR')
199  if(tread==1)then
200    dtset%znucl(1:dtset%npsp)=dprarr(1:dtset%npsp)
201  end if
202  if(tread/=1)then
203    write(message, '(3a)' )&
204 &   'The array znucl MUST be initialized in the input file while this is not done.',ch10,&
205 &   'Action: initialize znucl in your input file.'
206    MSG_ERROR(message)
207  end if
208 
209 !The default for ratsph has already been initialized
210  call intagm(dprarr,intarr,jdtset,marr,dtset%ntypat,string(1:lenstr),'ratsph',tread,'LEN')
211  if(tread==1)then
212    do ii=1,dtset%ntypat
213      dtset%ratsph(ii)=dprarr(ii)
214    end do
215  end if
216  ABI_ALLOCATE(ratsph,(dtset%ntypat))
217  do ii=1,dtset%ntypat
218    ratsph(ii)=dtset%ratsph(ii)
219  end do
220 
221 !Special treatment of _TYPAX (from a XYZ file), taking into account
222 !the fact that znucl does NOT depend on the dataset
223 !Examine all occurences of '_TYPAX'
224 
225  do
226    index_typsymb=index(string(1:lenstr),'_TYPAX')
227    if(index_typsymb==0)exit
228 !  Replace '_TYPAX' by '_TYPAT'
229    string(index_typsymb:index_typsymb+5)='_TYPAT'
230    index_upper=index_typsymb+5
231 !  Must start from the first blank after the tag (including possible dtset_char)
232    index_upper=index(string(index_upper:lenstr),blank)+index_upper-1
233    index_lower=index_upper
234 
235 !  Examine all atoms (the end of the symbol string is delimited by a XX )
236    do
237      index_blank=index(string(index_upper:lenstr),blank)+index_upper-1
238      string2=string(index_blank+1:index_blank+2)
239      if(string2=="XX")exit
240      found=0
241 !    Find the matching symbol
242      do ipsp=1,dtset%npsp
243        call atomdata_from_znucl(atom,dtset%znucl(ipsp))
244        symbol = atom%symbol
245        call inupper(symbol)
246        call inupper(string2)
247 
248 !      DEBUG
249 !      write(std_out,'(a)')' invars1 : before test, trim(adjustl(symbol)),trim(adjustl(string2))'
250 !      write(std_out,'(5a)' )'"',trim(adjustl(symbol)),'","',trim(adjustl(string2)),'"'
251 !      ENDDEBUG
252 
253        if(trim(adjustl(symbol))==trim(adjustl(string2)))then
254          found=1
255          index_upper=index_blank+1
256 !        Cannot deal properly with more that 9 psps
257          if(ipsp>=10)then
258            message = 'Need to use a pseudopotential with number larger than 9. Not allowed yet.'
259            MSG_ERROR(message)
260          end if
261 
262 !        DEBUG
263 !        write(std_out,*)' invars1 : found ipsp=',ipsp
264 !        ENDDEBUG
265 
266          write(string1,'(i1)')ipsp
267          string(index_lower:index_lower+1)=blank//string1
268          index_lower=index_lower+2
269        end if
270      end do ! ipsp
271 !    if not found ...
272      if(found==0)then
273        write(message,'(6a)' )&
274 &       'Did not find matching pseudopotential for XYZ atomic symbol,',ch10,&
275 &       'with value ',string2,ch10,&
276 &       'Action: check that the atoms required by the XYZ file correspond to one psp file.'
277        MSG_ERROR(message)
278      end if
279    end do ! Loop on atoms
280 !  One should find blanks after the last significant type value
281    string(index_lower:index_blank+2)=blank
282  end do ! loop to identify _TYPAX
283 
284 !---------------------------------------------------------------------------
285 
286 !Here, set up quantities that are related to geometrical description
287 !of the system (acell,rprim,xred), as well as
288 !initial velocity(vel), and spin of atoms (spinat), nuclear dipole moments
289 ! of atoms (nucdipmom),
290 !the symmetries (symrel,symafm, and tnons)
291 !and the list of fixed atoms (iatfix,iatfixx,iatfixy,iatfixz).
292 !Arrays have already been
293 !dimensioned thanks to the knowledge of msym and mxnatom
294 
295 !ji: We need to read the electric field before calling ingeo
296 !****** Temporary ******
297 
298  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'berryopt',tread,'INT')
299  if(tread==1) dtset%berryopt=intarr(1)
300 
301  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'berrysav',tread,'INT')
302  if(tread==1) dtset%berrysav=intarr(1)
303 
304  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'bfield',tread,'DPR')
305  if (tread==1) dtset%bfield(1:3) = dprarr(1:3)
306 
307  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'dfield',tread,'DPR')
308  if (tread==1) dtset%dfield(1:3) = dprarr(1:3)
309 
310  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'efield',tread,'DPR')
311  if (tread==1) dtset%efield(1:3) = dprarr(1:3)
312 
313  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'red_dfield',tread,'DPR')
314  if (tread==1) dtset%red_dfield(1:3) = dprarr(1:3)
315 
316  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'red_efield',tread,'DPR')
317  if (tread==1) dtset%red_efield(1:3) = dprarr(1:3)
318 
319  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'red_efieldbar',tread,'DPR')
320  if (tread==1) dtset%red_efieldbar(1:3) = dprarr(1:3)
321 
322  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'jfielddir',tread,'INT')
323  if(tread==1) dtset%jfielddir(1:3)=intarr(1:3)
324 
325 !We need to know nsppol/nspinor/nspden before calling ingeo
326  nsppol=dtset%nsppol
327  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nsppol',tread,'INT')
328  if(tread==1) nsppol=intarr(1)
329 
330 !Alternate SIESTA definition of nsppol
331  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'SpinPolarized',tread_alt,'LOG')
332  if(tread_alt==1)then
333    if(tread==1)then
334      write(message, '(a,a,a,a,a,a,a,a)' ) ch10,&
335 &     ' invars1: ERROR -',ch10,&
336 &     '  nsppol and SpinPolarized cannot be specified simultaneously',ch10,&
337 &     '  for the same dataset.',ch10,&
338 &     '  Action : check the input file.'
339      call wrtout(std_out,  message,'COLL')
340      leave=1
341    else
342 !    Note that SpinPolarized is a logical input variable
343      nsppol=1
344      if(intarr(1)==1)nsppol=2
345      tread=1
346    end if
347  end if
348  dtset%nsppol=nsppol
349 
350  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nspinor',tread,'INT')
351  if(tread==1) dtset%nspinor=intarr(1)
352 
353 !Has to read pawspnorb now, in order to adjust nspinor
354 !Also, if nspinor=1, turn on spin-orbit coupling by default, here for the PAW case. NC case is treated elsewhere.
355  if (dtset%usepaw>0)then
356 !  Change the default value
357    if(dtset%nspinor==2)dtset%pawspnorb=1
358    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'pawspnorb',tread,'INT')
359    if(tread==1)then
360      dtset%pawspnorb=intarr(1)
361      if(dtset%pawspnorb>0) dtset%nspinor=2
362    else
363      if(dtset%nspinor==2)then
364        write(message, '(4a)' ) ch10,&
365 &       ' invars1: COMMENT -',ch10,&
366 &       '  With nspinor=2 and usepaw=1, pawspnorb=1 has been switched on by default.'
367        call wrtout(iout,  message,'COLL')
368      end if
369    end if
370  end if
371  nspinor=dtset%nspinor
372 
373  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nspden',tread,'INT')
374  if(tread==1) then
375    dtset%nspden=intarr(1)
376  else
377    dtset%nspden=dtset%nsppol
378  end if
379 
380  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ntypalch',tread,'INT')
381  if(tread==1) dtset%ntypalch=intarr(1)
382 
383  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nzchempot',tread,'INT')
384  if(tread==1) dtset%nzchempot=intarr(1)
385 
386  ntypalch=dtset%ntypalch
387  if(ntypalch>ntypat)then
388    write(message, '(3a,i0,a,i0,a,a)' )&
389 &   'The input variable ntypalch must be smaller than ntypat, while it is',ch10,&
390 &   'ntypalch=',dtset%ntypalch,', and ntypat=',ntypat,ch10,&
391 &   'Action: check ntypalch vs ntypat in your input file.'
392    MSG_ERROR(message)
393  end if
394 
395  ntyppure=ntypat-ntypalch
396  dtset%ntyppure=ntyppure
397  npspalch=npsp-ntyppure
398  dtset%npspalch=npspalch
399  if(npspalch<0)then
400    write(message, '(a,i0,2a,i0,a,a)' )&
401 &   'The number of available pseudopotentials, npsp=',npsp,ch10,&
402 &   'is smaller than the requested number of types of pure atoms, ntyppure=',ntyppure,ch10,&
403 &   'Action: check ntypalch versus ntypat and npsp in your input file.'
404    MSG_ERROR(message)
405  end if
406 
407  if(ntypalch>0)then
408    call intagm(dprarr,intarr,jdtset,marr,ntypalch,string(1:lenstr),'algalch',tread,'INT')
409    if(tread==1) dtset%algalch(1:ntypalch)=intarr(1:ntypalch)
410  end if
411 
412 !Read the Zeeman field
413  call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'zeemanfield',tread,'BFI')
414  if(tread==1) then
415    if(dtset%nspden == 2)then
416      write(message,'(7a)')&
417 &     'A Zeeman field has been specified without noncollinear spins.',ch10,&
418 &     'Only the z-component of the magnetic field will be used.'
419      MSG_WARNING(message)
420    else if (dtset%nspden == 1)then
421      write(message, '(a,a,a)' )&
422 &     'A Zeeman field has been specified for a non-spin-polarized calculation.',ch10,&
423 &     'Action: check the input file.'
424      MSG_ERROR(message)
425    end if
426 
427    dtset%zeemanfield(1:3) = dprarr(1:3)
428  end if
429 
430 
431  ABI_ALLOCATE(amu,(ntypat))
432  ABI_ALLOCATE(mixalch,(npspalch,ntypalch))
433  ABI_ALLOCATE(vel,(3,natom))
434  ABI_ALLOCATE(vel_cell,(3,3))
435  ABI_ALLOCATE(xred,(3,natom))
436  intimage=2 ; if(dtset%nimage==1)intimage=1
437  do ii=1,dtset%nimage+1
438    iimage=ii
439    if(dtset%nimage==1 .and. ii==2)exit
440    if(dtset%nimage==2 .and. ii==3)exit
441    if(dtset%nimage> 2 .and. ii==intimage)cycle ! Will do the intermediate reference image at the last reading
442    if(dtset%nimage>=2 .and. ii==dtset%nimage+1)iimage=intimage
443 
444    write(message,'(a,i0)')' invars1 : treat image number: ',iimage
445    call wrtout(std_out,message,'COLL')
446 
447 !  Need to reset nsym to default value for each image
448    dtset%nsym=0
449 
450 !  Call ingeo for each image in turn, with the possible default values
451    acell=dtset%acell_orig(1:3,iimage)
452    amu=dtset%amu_orig(1:ntypat,iimage)
453    mixalch=dtset%mixalch_orig(1:npspalch,1:ntypalch,iimage)
454    rprim=dtset%rprim_orig(1:3,1:3,iimage)
455    vel=dtset%vel_orig(1:3,1:natom,iimage)
456    vel_cell=dtset%vel_cell_orig(1:3,1:3,iimage)
457    xred=dtset%xred_orig(1:3,1:natom,iimage)
458    ABI_ALLOCATE(iatfix,(3,natom))
459    ABI_ALLOCATE(nucdipmom,(3,natom))
460    ABI_ALLOCATE(spinat,(3,natom))
461    ABI_ALLOCATE(typat,(natom))
462    ABI_ALLOCATE(znucl,(dtset%npsp))
463    nucdipmom(1:3,1:natom)=dtset%nucdipmom(1:3,1:natom)
464    spinat(1:3,1:natom)=dtset%spinat(1:3,1:natom)
465    znucl(1:dtset%npsp)=dtset%znucl(1:dtset%npsp)
466 
467    call ingeo(acell,amu,dtset,bravais,dtset%genafm(1:3),iatfix,&
468 &   dtset%icoulomb,iimage,iout,jdtset,dtset%jellslab,lenstr,mixalch,&
469 &   msym,natom,dtset%nimage,dtset%npsp,npspalch,dtset%nspden,dtset%nsppol,&
470 &   dtset%nsym,ntypalch,dtset%ntypat,nucdipmom,dtset%nzchempot,&
471 &   dtset%pawspnorb,dtset%ptgroupma,ratsph,&
472 &   rprim,dtset%slabzbeg,dtset%slabzend,dtset%spgroup,spinat,&
473 &   string,symafm,dtset%symmorphi,symrel,tnons,dtset%tolsym,typat,vel,vel_cell,xred,znucl)
474    dtset%iatfix(1:3,1:natom)=iatfix(1:3,1:natom)
475    dtset%nucdipmom(1:3,1:natom)=nucdipmom(1:3,1:natom)
476    dtset%spinat(1:3,1:natom)=spinat(1:3,1:natom)
477    dtset%typat(1:natom)=typat(1:natom)
478    ABI_DEALLOCATE(iatfix)
479    ABI_DEALLOCATE(nucdipmom)
480    ABI_DEALLOCATE(spinat)
481    ABI_DEALLOCATE(typat)
482    ABI_DEALLOCATE(znucl)
483    dtset%acell_orig(1:3,iimage)=acell
484    dtset%amu_orig(1:ntypat,iimage)=amu
485    dtset%mixalch_orig(1:npspalch,1:ntypalch,iimage)=mixalch
486    dtset%rprim_orig(1:3,1:3,iimage)=rprim
487    dtset%vel_orig(1:3,1:natom,iimage)=vel
488    dtset%vel_cell_orig(1:3,1:3,iimage)=vel_cell
489    dtset%xred_orig(1:3,1:natom,iimage)=xred
490    call mkrdim(dtset%acell_orig(1:3,iimage),dtset%rprim_orig(1:3,1:3,iimage),dtset%rprimd_orig(1:3,1:3,iimage))
491  end do
492  ABI_DEALLOCATE(amu)
493  ABI_DEALLOCATE(mixalch)
494  ABI_DEALLOCATE(vel)
495  ABI_DEALLOCATE(vel_cell)
496  ABI_DEALLOCATE(xred)
497 
498 !Examine whether there is some vacuum space in the unit cell
499  call invacuum(jdtset,lenstr,natom,dtset%rprimd_orig(1:3,1:3,intimage),string,vacuum,&
500 & dtset%xred_orig(1:3,1:natom,intimage))
501 
502 !DEBUG
503 !write(std_out,*)' invars1: before inkpts, dtset%mixalch_orig(1:npspalch,1:ntypalch,:)=',&
504 !dtset%mixalch_orig(1:npspalch,1:ntypalch,1:dtset%nimage)
505 !stop
506 !ENDDEBUG
507 
508 !---------------------------------------------------------------------------
509 
510 !Set up k point grid number
511 !First, get additional information
512  dtset%kptopt=1
513  if(dtset%nspden==4)dtset%kptopt=4
514  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptopt',tread,'INT')
515  if(tread==1) dtset%kptopt=intarr(1)
516 
517  dtset%qptopt=1
518  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'qptopt',tread,'INT')
519  if(tread==1) dtset%qptopt=intarr(1)
520 
521  iscf=5
522  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'iscf',tread,'INT')
523  if(tread==1) iscf=intarr(1)
524 
525 
526  dtset%natsph=dtset%natom
527  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natsph',tread,'INT')
528  if(tread==1) dtset%natsph=intarr(1)
529 
530  dtset%natsph_extra=0
531  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natsph_extra',tread,'INT')
532  if(tread==1) dtset%natsph_extra=intarr(1)
533 
534  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'natvshift',tread,'INT')
535  if(tread==1) dtset%natvshift=intarr(1)
536 
537  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nconeq',tread,'INT')
538  if(tread==1) dtset%nconeq=intarr(1)
539 
540  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nkptgw',tread,'INT')
541  if(tread==1) dtset%nkptgw=intarr(1)
542  if (dtset%nkptgw<0) then
543    write(message, '(a,i0,a,a,a,a)' )&
544 &   'Input nkptgw must be >= 0, but was ',dtset%nkptgw,ch10,&
545 &   'This is not allowed.',ch10,&
546 &   'Action: check the input file.'
547    MSG_ERROR(message)
548  end if
549 
550 !Number of points for long wavelength limit.
551 !Default is dtset%gw_nqlwl=0
552  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'gw_nqlwl',tread,'INT')
553  if(tread==1) dtset%gw_nqlwl=intarr(1)
554  if (dtset%gw_nqlwl<0) then
555    write(message, '(a,i12,a,a,a,a)' )&
556 &   'Input gw_nqlwl must be > 0, but was ',dtset%gw_nqlwl,ch10,&
557 &   'This is not allowed.',ch10,&
558 &   'Action: check the input file.'
559    MSG_ERROR(message)
560  end if
561 
562 !Read number of k-points (if specified)
563  nkpt=0
564  if(dtset%kptopt==0)nkpt=1
565  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nkpt',tread,'INT')
566  if(tread==1) nkpt=intarr(1)
567  dtset%nkpt=nkpt
568 
569  call chkint_ge(0,0,cond_string,cond_values,ierr,'nkpt',nkpt,0,iout)
570  if(dtset%kptopt==0)then
571    cond_string(1)='kptopt' ; cond_values(1)=0
572    call chkint_ge(1,1,cond_string,cond_values,ierr,'nkpt',nkpt,1,iout)
573  end if
574 
575  nkpthf=nkpt
576  dtset%nkpthf=nkpt
577 
578 !Will compute the actual value of nkpt, if needed. Otherwise,
579 !test that the value of nkpt is OK, if kptopt/=0
580 !Set up dummy arrays istwfk, kpt, wtk
581 
582  if(nkpt/=0 .or. dtset%kptopt/=0)then
583    ABI_ALLOCATE(istwfk,(nkpt))
584    ABI_ALLOCATE(kpt,(3,nkpt))
585    ABI_ALLOCATE(kpthf,(3,nkpthf))
586    ABI_ALLOCATE(wtk,(nkpt))
587 !  Here, occopt is also a dummy argument
588    occopt=1 ; dtset%nshiftk=1 ; dtset%kptrlatt(:,:)=0
589 
590    kptrlen=20.0_dp ; wtk(:)=1.0_dp
591    dtset%shiftk(:,:)=half
592 
593    nqpt=0
594    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nqpt',tread,'INT')
595    if(tread==1) nqpt=intarr(1)
596 
597    chksymbreak=1
598    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'chksymbreak',tread,'INT')
599    if(tread==1) chksymbreak=intarr(1)
600 
601 
602 !  Use the first image to predict k and/or q points,
603 !  except if an intermediate image is available
604    intimage=1 ; if(dtset%nimage>2)intimage=(1+dtset%nimage)/2
605 
606 !  Find the q-point, if any.
607    if(nqpt==1)then
608      call inqpt(chksymbreak,std_out,jdtset,lenstr,msym,natom,dtset%qptn,dtset%wtq,&
609 &     dtset%rprimd_orig(1:3,1:3,intimage),dtset%spinat,string,dtset%typat,&
610 &     vacuum,dtset%xred_orig(1:3,1:natom,intimage),dtset%qptrlatt)
611    end if
612 
613 !  Find the k point grid
614    call inkpts(bravais,chksymbreak,dtset%fockdownsampling,iout,iscf,istwfk,jdtset,&
615 &   kpt,kpthf,dtset%kptopt,kptnrm,dtset%kptrlatt_orig,dtset%kptrlatt,kptrlen,lenstr,msym,&
616 &   nkpt,nkpthf,nqpt,dtset%ngkpt,dtset%nshiftk,dtset%nshiftk_orig,dtset%shiftk_orig,dtset%nsym,&
617 &   occopt,dtset%qptn,response,dtset%rprimd_orig(1:3,1:3,intimage),dtset%shiftk,&
618 &   string,symafm,symrel,vacuum,wtk)
619 
620    ABI_DEALLOCATE(istwfk)
621    ABI_DEALLOCATE(kpt)
622    ABI_DEALLOCATE(kpthf)
623    ABI_DEALLOCATE(wtk)
624 !  nkpt and nkpthf have been computed, as well as the k point grid, if needed
625    dtset%nkpt=nkpt
626    dtset%nkpthf=nkpthf
627  end if
628 
629  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nqptdm',tread,'INT')
630  if(tread==1) dtset%nqptdm=intarr(1)
631 
632  if (dtset%nqptdm<-1) then
633    write(message, '(a,i12,a,a,a,a)' )&
634 &   'Input nqptdm must be >= 0, but was ',dtset%nqptdm,ch10,&
635 &   'This is not allowed.',ch10,&
636 &   'Action: check the input file.'
637    MSG_ERROR(message)
638  end if
639 
640  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nzchempot',tread,'INT')
641  if(tread==1) dtset%nzchempot=intarr(1)
642 
643  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'cd_customnimfrqs',tread,'INT')
644  if(tread==1) dtset%cd_customnimfrqs=intarr(1)
645 
646  if (dtset%cd_customnimfrqs<0) then
647    write(message, '(a,i0,a,a,a,a)' )&
648 &   'Input cd_customnimfrqs must be >= 0, but was ',dtset%cd_customnimfrqs,ch10,&
649 &   'This is not allowed.',ch10,&
650 &   'Action: check the input file.'
651    MSG_ERROR(message)
652  end if
653 
654  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'gw_customnfreqsp',tread,'INT')
655  if(tread==1) dtset%gw_customnfreqsp=intarr(1)
656 
657  if (dtset%gw_customnfreqsp<0) then
658    write(message, '(a,i0,a,a,a,a)' )&
659 &   'Input gw_customnfreqsp must be >= 0, but was ',dtset%gw_customnfreqsp,ch10,&
660 &   'This is not allowed.',ch10,&
661 &   'Action: check the input file.'
662    MSG_ERROR(message)
663  end if
664 
665  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'gwls_n_proj_freq',tread,'INT')
666  if(tread==1) dtset%gwls_n_proj_freq=intarr(1)
667 
668  if (dtset%gwls_n_proj_freq<0) then
669    write(message, '(a,i0,a,a,a,a)' )&
670 &   'Input gwls_n_proj_freq must be >= 0, but was ',dtset%gwls_n_proj_freq,ch10,&
671 &   'This is not allowed.',ch10,&
672 &   'Action : check the input file.'
673    MSG_ERROR(message)
674  end if
675 
676  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'efmas_calc_dirs',tread,'INT')
677  if(tread==1) dtset%efmas_calc_dirs=intarr(1)
678 
679  if (ABS(dtset%efmas_calc_dirs)>3) then
680    write(message, '(a,i0,a,a,a,a)' )&
681 &   'Input efmas_calc_dirs must be between -3 and 3, but was ',dtset%efmas_calc_dirs,ch10,&
682 &   'This is not allowed.',ch10,&
683 &   'Action: check the input file.'
684    MSG_ERROR(message)
685  end if
686 
687  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'efmas_n_dirs',tread,'INT')
688  if(tread==1) dtset%efmas_n_dirs=intarr(1)
689 
690  if (dtset%efmas_n_dirs<0) then
691    write(message, '(a,i0,a,a,a,a)' )&
692 &   'Input efmas_n_dirs must be >= 0, but was ',dtset%efmas_n_dirs,ch10,&
693 &   'This is not allowed.',ch10,&
694 &   'Action: check the input file.'
695    MSG_ERROR(message)
696  end if
697 !---------------------------------------------------------------------------
698 
699  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nnos',tread,'INT')
700  if(tread==1) dtset%nnos=intarr(1)
701 
702  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ga_n_rules',tread,'INT')
703  if(tread==1) dtset%ga_n_rules=intarr(1)
704 
705 !Perform the first checks
706 
707  leave=0
708 
709 !Check that nkpt is greater than 0
710  if (nkpt<=0) then
711    write(message, '(a,a,a,a,i12,a,a,a,a)' ) ch10,&
712 &   ' invars1: ERROR -',ch10,&
713 &   '  After inkpts, nkpt must be > 0, but was ',nkpt,ch10,&
714 &   '  This is not allowed.',ch10,&
715 &   '  Action : check the input file.'
716    call wrtout(std_out,  message,'COLL')
717    leave=1
718  end if
719 
720 !Check that nsppol is 1 or 2
721  if (nsppol/=1 .and. nsppol/=2) then
722    write(message, '(a,a,a,a,i12,a,a,a,a)' ) ch10,&
723 &   ' invars1: ERROR -',ch10,&
724 &   '  Input nsppol must be 1 or 2, but was ',nsppol,ch10,&
725 &   '  This is not allowed.',ch10,&
726 &   '  Action : check the input file.'
727    call wrtout(std_out,message,'COLL')
728    leave=1
729  end if
730 
731 !Check that nspinor is 1 or 2
732  if (nspinor/=1 .and. nspinor/=2) then
733    write(message, '(a,a,a,a,i12,a,a,a,a)' ) ch10,&
734 &   ' invars1: ERROR -',ch10,&
735 &   '  Input nspinor must be 1 or 2, but was ',nspinor,ch10,&
736 &   '  This is not allowed.',ch10,&
737 &   '  Action : check the input file.'
738    call wrtout(std_out,message,'COLL')
739    leave=1
740  end if
741 
742 !Check that nspinor and nsppol are not 2 together
743  if (nsppol==2 .and. nspinor==2) then
744    write(message, '(8a)' ) ch10,&
745 &   ' invars1: ERROR -',ch10,&
746 &   '  nspinor and nsappol cannot be 2 together !',ch10,&
747 &   '  This is not allowed.',ch10,&
748 &   '  Action : check the input file.'
749    call wrtout(std_out,message,'COLL')
750    leave=1
751  end if
752 
753 !Here, leave if an error has been detected earlier
754  if(leave==1) then
755    message = ' Other errors might be present in the input file. '
756    MSG_ERROR(message)
757  end if
758 
759 
760 !Now, take care of mband_upper
761 
762  mband_upper=1
763  occopt=1
764  fband=0.5_dp
765 
766  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'occopt',tread,'INT')
767  if(tread==1) occopt=intarr(1)
768 
769 !Also read fband, that is an alternative to nband. The default
770 !is different for occopt==1 and for metallic occupations.
771  if(occopt==1)fband=0.125_dp
772  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'fband',tfband,'DPR')
773  if(tfband==1)fband=dprarr(1)
774 
775 !fband cannot be used when occopt==0 or occopt==2
776  if(tfband==1 .and. (occopt==0 .or. occopt==2) )then
777    write(message, '(3a)' )&
778 &   'fband cannot be used if occopt==0 or occopt==2 ',ch10,&
779 &   'Action: correct your input file, suppress fband, or change occopt.'
780    MSG_ERROR(message)
781  end if
782 
783  ABI_ALLOCATE(nband,(nkpt*nsppol))
784  tnband=0
785 
786 !Compute ziontypat
787 !When the pseudo-atom is pure, simple copy
788  if(ntyppure>0)then
789    do itypat=1,ntyppure
790      dtset%ziontypat(itypat)=zionpsp(itypat)
791    end do
792  end if
793 !When the pseudo-atom is alchemical, must make mixing
794  if(ntypalch>0)then
795    do itypat=ntyppure+1,ntypat
796      dtset%ziontypat(itypat)=zero
797      do ipsp=ntyppure+1,npsp
798        dtset%ziontypat(itypat)=dtset%ziontypat(itypat) &
799 &       +dtset%mixalch_orig(ipsp-ntyppure,itypat-ntyppure,1)*zionpsp(ipsp)
800      end do
801    end do
802  end if
803 
804 
805 
806  if (occopt==0 .or. occopt==1 .or. (occopt>=3 .and. occopt<=8) ) then
807    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nband',tnband,'INT')
808 !  Note : mband_upper is initialized, not nband
809    if(tnband==1) mband_upper=intarr(1)
810 
811    if(tfband==1 .and. tnband==1)then
812      write(message, '(3a)' )&
813 &     'fband and nband cannot be used together. ',ch10,&
814 &     'Action: correct your input file, suppress either fband or nband.'
815      MSG_ERROR(message)
816    end if
817 
818 !  In case nband was not read, use fband, either read, or the default,
819 !  to provide an upper limit for mband_upper
820    if(tnband==0)then
821 
822      charge=0.0_dp
823      call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'charge',tread,'DPR')
824      if(tread==1) charge=dprarr(1)
825 
826 !    Only take into account negative charge, to compute maximum number of bands
827      if(charge > 0.0_dp)charge=0.0_dp
828 
829 !     mband_upper=nspinor*((nint(zion_max)*natom+1)/2 - floor(charge/2.0_dp)&
830 !&     + ceiling(fband*natom-1.0d-10))
831      zval=0.0_dp
832      do iatom=1,natom
833        zval=zval+dtset%ziontypat(dtset%typat(iatom))
834      end do
835      zelect=zval-charge
836      mband_upper=nspinor * ((ceiling(zelect-1.0d-10)+1)/2 + ceiling( fband*natom - 1.0d-10 ))
837 
838      nband(:)=mband_upper
839 
840 !    DEBUG
841 !    write(std_out,*)' invars1 : zion_max,natom,fband,mband_upper '
842 !    write(std_out,*)zion_max,natom,fband,mband_upper
843 !    ENDDEBUG
844 
845    end if
846 
847    nband(:)=mband_upper
848 
849  else if (occopt==2) then
850    ABI_ALLOCATE(reaalloc,(nkpt*nsppol))
851    call intagm(reaalloc,nband,jdtset,nkpt*nsppol,nkpt*nsppol,string(1:lenstr),'nband',tnband,'INT')
852    if(tnband==1)then
853      do ikpt=1,nkpt*nsppol
854        if (nband(ikpt)>mband_upper) mband_upper=nband(ikpt)
855      end do
856    end if
857    ABI_DEALLOCATE(reaalloc)
858  else
859    write(message, '(a,i0,3a)' )&
860 &   'occopt=',occopt,' is not an allowed value.',ch10,&
861 &   'Action: correct your input file.'
862    MSG_ERROR(message)
863  end if
864 
865 !Check that mband_upper is greater than 0
866  if (mband_upper<=0) then
867    write(message, '(a,i0,4a)' )&
868 &   'Maximal nband must be > 0, but was ',mband_upper,ch10,&
869 &   'This is not allowed.',ch10,&
870 &   'Action: check the input file.'
871    MSG_ERROR(message)
872  end if
873 
874 !The following 3 values are needed to dimension the parallelism over images
875  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'imgmov',tread,'INT')
876  if(tread==1) dtset%imgmov=intarr(1)
877  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ntimimage',tread,'INT')
878  if(tread==1) dtset%ntimimage=intarr(1)
879  call intagm(dprarr,intarr,jdtset,marr,dtset%nimage,string(1:lenstr),'dynimage',tread,'INT')
880  if(tread==1)then
881    dtset%dynimage(1:dtset%nimage)=intarr(1:dtset%nimage)
882  else if (dtset%imgmov==2.or.dtset%imgmov==5) then
883    dtset%dynimage(1)=0;dtset%dynimage(dtset%nimage)=0
884  end if
885  dtset%ndynimage=count(dtset%dynimage(1:dtset%nimage)/=0)
886 
887  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'wfoptalg',tread,'INT')
888  if(tread==1) then
889    dtset%wfoptalg=intarr(1)
890  else
891    if (dtset%usepaw==0)    dtset%wfoptalg=0
892    if (dtset%usepaw/=0)    dtset%wfoptalg=10
893    if (dtset%optdriver==RUNL_GSTATE) then
894      if (dtset%paral_kgb/=0) dtset%wfoptalg=14
895    end if
896  end if
897 
898 !---------------------------------------------------------------------------
899 !Some PAW+DMFT keywords
900  dtset%usedmft=0
901  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'usedmft',tread,'INT')
902  if(tread==1) dtset%usedmft=intarr(1)
903 
904 !Some ucrpa keywords
905  dtset%ucrpa=0
906  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ucrpa',tread,'INT')
907  if(tread==1) dtset%ucrpa=intarr(1)
908 
909  if (dtset%ucrpa > 0 .and. dtset%usedmft > 0) then
910    write(message, '(7a)' )&
911 &   'usedmft and ucrpa are both activated in the input file ',ch10,&
912 &   'In the following, abinit assume you are doing a ucrpa calculation and ',ch10,&
913 &   'you define Wannier functions as in DFT+DMFT calculation',ch10,&
914 &   'If instead, you want to do a full dft+dmft calculation and not only the Wannier construction, use ucrpa=0'
915    MSG_WARNING(message)
916  end if
917 
918 !Some PAW+U keywords
919  dtset%usepawu=0
920  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'usepawu',tread,'INT')
921  if(tread==1) dtset%usepawu=intarr(1)
922  if ( dtset%usedmft > 0 .and. dtset%usepawu >= 0 ) dtset%usepawu = 1
923 
924  if (dtset%usepawu > 0 ) then
925    write(message, '(7a)' )&
926 &   'usedmft and usepawu are both activated ',ch10,&
927 &   'This is not an usual calculation:',ch10,&
928 &   'usepawu will be put to a value >= 10:',ch10,&
929 &   'LDA+U potential and energy will be put to zero'
930    MSG_WARNING(message)
931  end if
932 
933  dtset%usedmatpu=0
934  dtset%lpawu(1:dtset%ntypat)=-1
935  if (dtset%usepawu>0.or.dtset%usedmft>0) then
936 
937    call intagm(dprarr,intarr,jdtset,marr,dtset%ntypat,string(1:lenstr),'lpawu',tread,'INT')
938    if(tread==1) dtset%lpawu(1:dtset%ntypat)=intarr(1:dtset%ntypat)
939 
940    call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'usedmatpu',tread,'INT')
941    if(tread==1) dtset%usedmatpu=intarr(1)
942 
943  end if
944 
945 !Some PAW+Exact exchange keywords
946  dtset%useexexch=0
947  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'useexexch',tread,'INT')
948  if(tread==1) dtset%useexexch=intarr(1)
949 
950 
951  dtset%lexexch(1:dtset%ntypat)=-1
952 
953  if (dtset%useexexch>0) then
954    call intagm(dprarr,intarr,jdtset,marr,dtset%ntypat,string(1:lenstr),'lexexch',tread,'INT')
955    if(tread==1) dtset%lexexch(1:dtset%ntypat)=intarr(1:dtset%ntypat)
956  end if
957 ! LDA minus half keyword
958  call intagm(dprarr,intarr,jdtset,marr,dtset%ntypat,string(1:lenstr),'ldaminushalf',tread,'INT')
959  if(tread==1) dtset%ldaminushalf(1:dtset%ntypat)=intarr(1:dtset%ntypat)
960 !Some plowan data
961  dtset%plowan_natom=0
962  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'plowan_natom',tread,'INT')
963  if(tread==1) dtset%plowan_natom=intarr(1)
964 
965  dtset%plowan_nt=0
966  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'plowan_nt',tread,'INT')
967  if(tread==1) dtset%plowan_natom=intarr(1)
968 
969 
970 !PAW potential zero keyword
971  dtset%usepotzero=0
972  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'usepotzero',tread,'INT')
973  if(tread==1) dtset%usepotzero=intarr(1)
974 
975 !Macro_uj (determination of U in PAW+U), governs also allocation of atvshift
976  dtset%macro_uj = 0
977  call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'macro_uj',tread,'INT')
978  if(tread==1) dtset%macro_uj=intarr(1)
979 
980  ABI_DEALLOCATE(nband)
981  ABI_DEALLOCATE(ratsph)
982  ABI_DEALLOCATE(intarr)
983  ABI_DEALLOCATE(dprarr)
984 
985 end subroutine invars1