1!! NAME
2!!  input_dataset_mod
3!!
4!! FUNCTION
5!!  This module contains objects and routines used to parse the atompaw
6!!  input file. The complete input dataset is stored in the `input_dataset`
7!!  Fortran data-structure.
8
9#if defined HAVE_CONFIG_H
10#include "config.h"
11#endif
12
13MODULE input_dataset_mod
14
15 USE io_tools
16 USE tools
17
18 IMPLICIT NONE
19
20 PRIVATE
21
22!Public functions
23 PUBLIC :: input_dataset_read        ! Initialize input_dataset by reading it from file
24 PUBLIC :: input_dataset_free        ! Destroy input_dataset data-structure
25 PUBLIC :: input_dataset_copy        ! Copy a input data-structure into another
26 PUBLIC :: input_dataset_read_occ    ! Read and change modified occupations in input_dataset
27 PUBLIC :: input_dataset_read_abinit ! Read ABINIT options and store them in input_dataset
28 PUBLIC :: input_dataset_read_xml    ! Read XML options and store them in input_dataset
29 PUBLIC :: input_dataset_read_upf    ! Read UPF options and store them in input_dataset
30
31!Public structured datatype
32!All data from input file
33 TYPE,PUBLIC :: input_dataset_t
34   CHARACTER(2) :: atomic_symbol    ! Atomic symbol
35   INTEGER      :: atomic_charge    ! Atomic charge
36   LOGICAL :: scalarrelativistic    ! Flag activating the scalar relativistic scheme
37   LOGICAL :: diracrelativistic     ! Flag activating the Dirac relativistic scheme
38   LOGICAL :: finitenucleus         ! Flag activating finite nucleus model
39   INTEGER :: finitenucleusmodel    ! Option for the finite nucleus model
40   LOGICAL :: HFpostprocess         ! Option for the post-processing of Hartree-Fock
41   LOGICAL :: localizedcoreexchange ! Flag activating the localization of the core exchange (HF)
42   CHARACTER(10) :: gridkey         ! Type of radial grid
43   REAL(8) :: gridrange             ! Range of the radial grid
44   REAL(8) :: gridmatch             ! A matching radius in the radial grid
45   INTEGER :: gridpoints            ! Number of points of the radial grid
46   INTEGER :: nlogderiv             ! Number of points for the logarithmic derivative calculation
47   REAL(8) :: minlogderiv           ! Minimum energy value for the logarithmic derivative calculation
48   REAL(8) :: maxlogderiv           ! Maximum energy value for the logarithmic derivative calculation
49   LOGICAL :: BDsolve               ! Flag activating the Block-Davidson solver
50   LOGICAL :: fixed_zero            ! Flag activating the "fixed zero" exact exchange potential calculation
51   INTEGER :: fixed_zero_index      ! Option for "fixed zero" calculation in the exact exchange potential
52   CHARACTER(132) :: exctype        ! Exchange-correlation type (string)
53   INTEGER :: np(5)                 ! Electronic configuration: number of s,p,d,f,g shells
54   INTEGER :: norbit                ! Electronic configuration: number of orbitals
55   LOGICAL,ALLOCATABLE :: orbit_iscore(:)  ! Electronic configuration: TRUE for the core orbitals
56   INTEGER :: norbit_mod            ! Electronic configuration: number of orbitals with modified occupations
57   INTEGER,ALLOCATABLE :: orbit_mod_n(:)   ! Electronic config.: n number of the modified orbital
58   INTEGER,ALLOCATABLE :: orbit_mod_l(:)   ! Electronic config.: l number of the modified orbital
59   INTEGER,ALLOCATABLE :: orbit_mod_k(:)   ! Electronic config.: kappa number of the modified orbital
60   REAL(8),ALLOCATABLE :: orbit_mod_occ(:) ! Electronic config.: occupation of the modified orbital
61   INTEGER :: norbit_val            ! Electronic configuration: number of valence orbitals
62   INTEGER,ALLOCATABLE :: orbit_val_n(:)   ! Electronic config.: n number of the valence orbital
63   INTEGER,ALLOCATABLE :: orbit_val_l(:)   ! Electronic config.: l number of the valence orbital
64   INTEGER,ALLOCATABLE :: orbit_val_k(:)   ! Electronic config.: kappa number of the valence orbital
65   INTEGER :: lmax=-1               ! PAW Basis: maximum l value
66   REAL(8) :: rc=0.d0               ! PAW basis: cut-off radius for the augmentation regions
67   REAL(8) :: rc_shap=0.d0          ! PAW basis: cut-off radius of the compensation charge shape function
68   REAL(8) :: rc_vloc=0.d0          ! PAW basis: matching radius for the local potential
69   REAL(8) :: rc_core=0.d0          ! PAW basis: matching radius for the pseudo-core density
70   INTEGER :: nbasis                ! PAW basis : number of basis functions
71   INTEGER :: nbasis_add            ! PAW basis: number of additional basis functions (unbound states)
72   INTEGER,ALLOCATABLE :: basis_add_l(:)      ! PAW basis: l number for the additional basis func.
73   INTEGER,ALLOCATABLE :: basis_add_k(:)      ! PAW basis: kappa number for the additional basis func.
74   REAL(8),ALLOCATABLE :: basis_add_energy(:) ! PAW basis: ref. energy for the additional basis func.
75   REAL(8),ALLOCATABLE :: basis_func_rc(:)    ! PAW basis: matching radii for the pseudo partial waves
76   INTEGER :: projector_type        ! Type of projectors (Bloechl, Vanderbilt, ...)
77   INTEGER :: pseudo_type           ! Type of pseudization scheme (Bessel, polynom, ...)
78   INTEGER :: ortho_type            ! Type of orthogonalization scheme (Gram-Schmidt, ...)
79   INTEGER :: pseudo_polynom2_pdeg  ! Polynom2 projectors: degree of the polynom
80   REAL(8) :: pseudo_polynom2_qcut  ! Polynom2 projectors: q-value for Fourier filtering
81   INTEGER :: shapefunc_type           ! Compensation shape function type (sinc2, gaussian, ...)
82   REAL(8) :: shapefunc_gaussian_param ! Compensation shape function: parameter for gaussian type
83   REAL(8) :: hf_coretol            ! Tolerance for core density (Hartree-Fock only)
84   INTEGER :: vloc_type             ! Type of local potential pseudization
85   INTEGER :: vloc_l                ! Local potential: l quantum number (MTrouillier, Ultrasoft)
86   REAL(8) :: vloc_ene              ! Local potential: reference energy (MTrouillier, Ultrasoft)
87   REAL(8) :: vloc_setvloc_coef     ! "SetVloc" local potential: coefficient
88   REAL(8) :: vloc_setvloc_rad      ! "SetVloc" local potential: radius
89   INTEGER :: vloc_kerker_power(4)  ! "Kerker" locazl potential: polynomial powers
90   LOGICAL :: abinit_usexcnhat      ! ABINIT option: TRUE if nhat density (compensation) is included in XC
91   LOGICAL :: abinit_prtcorewf      ! ABINIT option: TRUE is printing of core WF is required
92   CHARACTER(132) :: abinit_author  ! ABINIT option: author to be printed in the PAW dataset file
93   LOGICAL :: abinit_uselog         ! ABINIT option: TRUE if data are transfered on a log. grid
94   INTEGER :: abinit_log_meshsz     ! ABINIT option: mesh size for the logarithmic grid
95   REAL(8) :: abinit_log_step       ! ABINIT option: logarithmic step for the logarithmic grid
96   LOGICAL :: abinit_userso         ! ABINIT option: TRUE if REAL Space Optimization is required
97   REAL(8) :: abinit_rso_ecut       ! ABINIT option: Real Space Optimization parameter - plane wave cutoff
98   REAL(8) :: abinit_rso_gfact      ! ABINIT option: Real Space Optimization parameter - Gamma/Gmax ratio
99   REAL(8) :: abinit_rso_werror     ! ABINIT option: Real Space Optimization parameter: max. error
100   LOGICAL :: xml_usexcnhat         ! XML option: TRUE if nhat density (compensation) is included in XC
101   LOGICAL :: xml_prtcorewf         ! XML option: TRUE is printing of core WF is required
102   CHARACTER(132) :: xml_author     ! XML option: author to be printed in the PAW dataset file
103   CHARACTER(132) :: xml_comment    ! XML option: comment to be printed in the header of PAW dataset file
104   LOGICAL :: xml_usespl            ! XML option: TRUE if data are interpolated on a log. grid
105   INTEGER :: xml_spl_meshsz        ! XML option: mesh size for the reduced grid
106   LOGICAL :: xml_userso            ! XML option: TRUE if REAL Space Optimization is required
107   REAL(8) :: xml_rso_ecut          ! XML option: Real Space Optimization parameter - plane wave cutoff
108   REAL(8) :: xml_rso_gfact         ! XML option: Real Space Optimization parameter - Gamma/Gmax ratio
109   REAL(8) :: xml_rso_werror        ! XML option: Real Space Optimization parameter: max. error
110   LOGICAL :: xml_uselda12          ! XML option: TRUE if LDA-1/2 potential calculation is required
111   INTEGER :: xml_lda12_orb_l       ! XML option: LDA-1/2 parameter - l quantum number of the ionized orbital
112   INTEGER :: xml_lda12_orb_n       ! XML option: LDA-1/2 parameter - n quantum number of the ionized orbital
113   REAL(8) :: xml_lda12_ion         ! XML option: LDA-1/2 parameter - amount of charge removed from the ionized orbital
114   REAL(8) :: xml_lda12_rcut        ! XML option: LDA-1/2 parameter - cut-off radius (in bohr)
115   CHARACTER(15) :: xml_lda12_logfile ! XML option: LDA-1/2 parameter - name of the log file
116   REAL(8) :: upf_grid_xmin         ! UPF option: minimum radius given by the grid
117   REAL(8) :: upf_grid_zmesh        ! UPF option: inverse of the rdial step of the grid
118   REAL(8) :: upf_grid_dx           ! UPF option: logarithmic step of the grid
119   REAL(8) :: upf_grid_range        ! UPF option: range of the grid
120 END TYPE input_dataset_t
121
122!Public variable containing the complete input dataset
123 TYPE(input_dataset_t),PUBLIC,SAVE,TARGET :: input_dataset
124
125!Public parameters
126 INTEGER,PARAMETER,PUBLIC :: PROJECTOR_TYPE_BLOECHL   = 1
127 INTEGER,PARAMETER,PUBLIC :: PROJECTOR_TYPE_VANDERBILT= 2
128 INTEGER,PARAMETER,PUBLIC :: PROJECTOR_TYPE_CUSTOM    = 3
129 INTEGER,PARAMETER,PUBLIC :: PROJECTOR_TYPE_MODRRKJ   = 4
130 INTEGER,PARAMETER,PUBLIC :: PROJECTOR_TYPE_HF        = 5
131 INTEGER,PARAMETER,PUBLIC :: PSEUDO_TYPE_BLOECHL      = 1
132 INTEGER,PARAMETER,PUBLIC :: PSEUDO_TYPE_POLYNOM      = 2
133 INTEGER,PARAMETER,PUBLIC :: PSEUDO_TYPE_POLYNOM2     = 3
134 INTEGER,PARAMETER,PUBLIC :: PSEUDO_TYPE_RRKJ         = 4
135 INTEGER,PARAMETER,PUBLIC :: PSEUDO_TYPE_BLOECHL_K    = 5
136 INTEGER,PARAMETER,PUBLIC :: PSEUDO_TYPE_HF           = 6
137 INTEGER,PARAMETER,PUBLIC :: ORTHO_TYPE_GRAMSCHMIDT   = 1
138 INTEGER,PARAMETER,PUBLIC :: ORTHO_TYPE_VANDERBILT    = 2
139 INTEGER,PARAMETER,PUBLIC :: ORTHO_TYPE_SVD           = 3
140 INTEGER,PARAMETER,PUBLIC :: ORTHO_TYPE_HF            = 4
141 INTEGER,PARAMETER,PUBLIC :: SHAPEFUNC_TYPE_GAUSSIAN  = 1
142 INTEGER,PARAMETER,PUBLIC :: SHAPEFUNC_TYPE_SINC      = 2
143 INTEGER,PARAMETER,PUBLIC :: SHAPEFUNC_TYPE_BESSEL    = 3
144 INTEGER,PARAMETER,PUBLIC :: VLOC_TYPE_MTROULLIER     = 1
145 INTEGER,PARAMETER,PUBLIC :: VLOC_TYPE_ULTRASOFT      = 2
146 INTEGER,PARAMETER,PUBLIC :: VLOC_TYPE_BESSEL         = 3
147 INTEGER,PARAMETER,PUBLIC :: VLOC_TYPE_SETVLOC        = 4
148 INTEGER,PARAMETER,PUBLIC :: VLOC_TYPE_KERKER_EXPF    = 5
149 INTEGER,PARAMETER,PUBLIC :: VLOC_TYPE_KERKER_POLY    = 6
150 INTEGER,PARAMETER,PUBLIC :: UNKNOWN_TYPE             =-1
151
152 !Default values
153 LOGICAL,PARAMETER,PRIVATE :: PRTCOREWF_DEF      =.false.
154 LOGICAL,PARAMETER,PRIVATE :: USEXCNHAT_DEF      =.false.
155 LOGICAL,PARAMETER,PRIVATE :: USELOG_DEF         =.false.
156 LOGICAL,PARAMETER,PRIVATE :: USERSO_DEF         =.false.
157 LOGICAL,PARAMETER,PRIVATE :: USELDA12_DEF       =.false.
158 INTEGER,PARAMETER,PRIVATE :: LOGGRD_SIZE_DEF    =350
159 REAL(8),PARAMETER,PRIVATE :: LOGGRD_STEP_DEF    =0.035d0
160 REAL(8),PARAMETER,PRIVATE :: RSO_ECUT_DEF       =10.0d0
161 REAL(8),PARAMETER,PRIVATE :: RSO_GFACT_DEF      =2.d0
162 REAL(8),PARAMETER,PRIVATE :: RSO_WERROR_DEF     =0.0001d0
163 INTEGER,PARAMETER,PRIVATE :: LDA12_ORB_L_DEF    =-1
164 INTEGER,PARAMETER,PRIVATE :: LDA12_ORB_N_DEF    =-1
165 REAL(8),PARAMETER,PRIVATE :: LDA12_ION_DEF      =0.d0
166 REAL(8),PARAMETER,PRIVATE :: LDA12_RCUT_DEF     =0.d0
167 CHARACTER(15),PARAMETER,PRIVATE :: lda12_logfile='lda-12.log'
168 REAL(8),PARAMETER,PRIVATE :: UPF_DX_DEF         =0.005d0
169 REAL(8),PARAMETER,PRIVATE :: UPF_XMIN_DEF       =-9.d0
170 REAL(8),PARAMETER,PRIVATE :: UPF_ZMESH_DEF      =1.d0
171 REAL(8),PARAMETER,PRIVATE :: UPF_RANGE_DEF      =15.d0
172 CHARACTER(132),PARAMETER,PRIVATE :: AUTHOR_DEF  ='', COMMENT_DEF=''
173
174!Private parameters
175 INTEGER,PARAMETER,PRIVATE :: norbit_max=50,nbasis_add_max=25
176 REAL(8),PARAMETER,PRIVATE :: linrange=50.d0,mxgridlin=20001
177 REAL(8),PARAMETER,PRIVATE :: logrange=80.d0,v4logrange=100.d0,mxgridlog=2001
178 REAL(8),PARAMETER,PRIVATE :: logder_min=-5.d0,logder_max=4.95d0,logder_pts=200
179 REAL(8),PARAMETER,PRIVATE :: polynom2_pdeg_def=4
180 REAL(8),PARAMETER,PRIVATE :: polynom2_qcut_def=10.d0
181 REAL(8),PARAMETER,PRIVATE :: gausstol_def=1.d-4
182 REAL(8),PARAMETER,PRIVATE :: hf_coretol_def=1.d-4
183
184!Private variables
185 LOGICAL,PRIVATE,SAVE :: has_to_print=.true.
186 LOGICAL,PRIVATE,SAVE :: has_to_ask=.true.
187
188CONTAINS
189
190!!=================================================================
191!! NAME
192!!  input_dataset_read
193!!
194!! FUNCTION
195!!  Initialize an input_dataset datastructure by reading it from
196!!  a file. If file is omitted, then read from standard input.
197!!  Note: we only read here data used to generate the PAW dataset,
198!!    not data used for the post-processing (output, explore, scfpaw, ...)
199!!
200!! INPUTS (all optionals)
201!!  [inputfile]= name of input file to be read
202!!  [echofile]= name of a file to echo input file content
203!!  [read_global_data]= if TRUE, read global data (atom, XC, grid, ...) - Default TRUE
204!!  [read_elec_data]= if TRUE, read electronic configuration (orbital & occupations) - Default TRUE
205!!  [read_coreval_data]= if TRUE, read electronic config (core and valence) - Default TRUE
206!!  [read_basis_data]= if TRUE, read basis data (radii, pseudo scheme, ...) - Default TRUE
207!!
208!! OUTPUT
209!!  [input_dt]= datastructure containing the complete input file.
210!!              If omitted, then the global public `input_dataset`
211!!              is used.
212!!
213!!=================================================================
214 SUBROUTINE input_dataset_read(input_dt,inputfile,echofile,&
215&           read_global_data,read_elec_data,read_coreval_data,read_basis_data)
216
217!---- Arguments
218 CHARACTER*(*),INTENT(IN),OPTIONAL :: inputfile,echofile
219 LOGICAL,INTENT(IN),OPTIONAL :: read_global_data,read_elec_data,&
220&                               read_coreval_data,read_basis_data
221 TYPE(input_dataset_t),INTENT(INOUT),OPTIONAL,TARGET :: input_dt
222
223!---- Local variables
224 INTEGER,PARAMETER :: ifunit=111,ecunit=222
225 INTEGER,PARAMETER :: nkappa(5)=(/1,2,2,2,2/)
226 INTEGER :: input_unit
227 INTEGER :: ii,io,nadd,norb,nval,nbl,nn,ik,kk
228 INTEGER :: ilin,ilog,inrl,iscl,ipnt,ifin,iend,ihfpp,ilcex
229 INTEGER :: igrid,irelat,ilogder,ilogv4,ibd,idirac,ifixz,ll,nstart
230 LOGICAL :: has_to_echo
231 LOGICAL :: read_global_data_,read_elec_data_,read_coreval_data_,read_basis_data_
232 CHARACTER(200) :: inputline,inputword
233 CHARACTER(128) :: exchangecorrelationandgridline
234 CHARACTER(1) :: CHR
235 REAL(8) :: x1,x2
236 INTEGER :: basis_add_l(nbasis_add_max)
237 INTEGER :: basis_add_k(nbasis_add_max)
238 REAL(8) :: basis_add_energy(nbasis_add_max)
239 TYPE(input_dataset_t),POINTER :: dataset
240
241!------------------------------------------------------------------
242
243!Select datastruture to be read
244 IF (PRESENT(input_dt)) THEN
245   dataset => input_dt
246 ELSE
247   dataset => input_dataset
248 ENDIF
249
250!Select input file logical unit
251 IF (PRESENT(inputfile)) THEN
252   input_unit=ifunit
253   OPEN(ifunit,file=trim(inputfile),form='formatted')
254 ELSE
255   input_unit=STD_IN
256 END IF
257
258!Check if we read a TTY or a FILE
259 has_to_ask=unit_isatty(input_unit)
260
261!Do we echo input file content?
262 has_to_echo=PRESENT(echofile)
263 IF (has_to_echo) THEN
264   OPEN(ecunit,file=trim(echofile),form='formatted')
265 END IF
266
267!Select which components have to be read
268 read_global_data_=.true.;if (PRESENT(read_global_data)) read_global_data_=read_global_data
269 read_elec_data_=.true.;if (PRESENT(read_elec_data)) read_elec_data_=read_elec_data
270 read_coreval_data_=.true.;if (PRESENT(read_coreval_data)) read_coreval_data_=read_coreval_data
271 read_basis_data_=.true.;if (PRESENT(read_basis_data)) read_basis_data_=read_basis_data
272
273!Print a title
274 IF (read_global_data_.OR.read_elec_data_.OR.read_coreval_data_.OR.read_basis_data_) THEN
275   IF (has_to_ask) THEN
276     WRITE(STD_OUT,'(/,1x,a)') "===== READING OF INPUT DATA ====="
277   ELSE
278     WRITE(STD_OUT,'(/,3x,a)') "===== READING OF INPUT FILE ====="
279   END IF
280 END IF
281
282
283!------------------------------------------------------------------
284!Start reading of AE data
285 IF (read_global_data_) THEN
286
287!------------------------------------------------------------------
288!=== 1st line: read atomic symbol, atomic number
289
290 IF (has_to_ask) THEN
291   WRITE(STD_OUT,*) 'Enter atomic symbol and atomic number'
292 END IF
293
294 READ(input_unit,'(a)') inputline
295 IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
296 CALL eliminate_comment(inputline)
297
298 READ(inputline,*) dataset%atomic_symbol,dataset%atomic_charge
299
300!Print read data
301 IF (has_to_print) THEN
302   WRITE(STD_OUT,'(3x,a,a2)') "Atomic symbol : ",dataset%atomic_symbol
303   WRITE(STD_OUT,'(3x,a,i0)') "Atomic charge : ",dataset%atomic_charge
304 END IF
305
306!------------------------------------------------------------------
307!=== 2nd line: read XC type, grid data, relativistic,point-nucleus,
308!              logderiv data, HF data, Block-Davidson keyword
309
310 IF (has_to_ask) THEN
311   WRITE(STD_OUT,*) 'Enter exchange-correlation type, among the following:'
312   WRITE(STD_OUT,*) '    * LDA-PW (default), GGA-PBE, GGA-PBESOL'
313   WRITE(STD_OUT,*) '    * LibXC keyword beginning with XC_'
314   WRITE(STD_OUT,*) '    * EXX, EXXOCC, EXXKLI, EXXCS'
315   WRITE(STD_OUT,*) '    * HF, HFV'
316   WRITE(STD_OUT,*) ' further optionally (space) "nonrelativistic/scalarrelativistic/diracrelativistic" keyword'
317   WRITE(STD_OUT,*) ' further optionally (space) "point-nucleus/finite-nucleus" keyword'
318   WRITE(STD_OUT,*) ' further optionally (space) "loggrid/lineargrid" keyword if appropriate'
319   WRITE(STD_OUT,*) '   Note: "loggridv4" puts more points near origin'
320   WRITE(STD_OUT,*) ' further optionally n (number of grid points)'
321   WRITE(STD_OUT,*) '                       r_max (max. grid radius)'
322   WRITE(STD_OUT,*) '                       r_match (exact value of r(n))'
323   WRITE(STD_OUT,*) ' further optionally (space) "logderivrange" keyword'
324   WRITE(STD_OUT,*) '  further optionally emin (minimum energy for log. deriv. plot)'
325   WRITE(STD_OUT,*) '                     emax (maximum energy for log. deriv. plot)'
326   WRITE(STD_OUT,*) '                     ne   (#  of energies for log. deriv. plot)'
327   WRITE(STD_OUT,*) ' and optionally "BDsolve" keyword for Block-Davidson solver'
328 END IF
329
330!Read full line
331 READ(input_unit,'(a)') exchangecorrelationandgridline
332 IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(exchangecorrelationandgridline)
333 CALL eliminate_comment(inputline)
334
335 CALL Uppercase(exchangecorrelationandgridline)
336 exchangecorrelationandgridline=trim(exchangecorrelationandgridline)
337
338!Retrieve keyword indexes
339 ilin=0;ilin=0;ilog=0;ilogv4=0;inrl=0;iscl=0;ipnt=0;ifin=0
340 ihfpp=0;ilcex=0;igrid=0;irelat=0;ilogder=0;ibd=0;idirac=0
341 ilin=INDEX(exchangecorrelationandgridline,'LINEARGRID')
342 ilog=INDEX(exchangecorrelationandgridline,'LOGGRID')
343 ilogv4=INDEX(exchangecorrelationandgridline,'LOGGRIDV4')
344 ibd=INDEX(exchangecorrelationandgridline,'BDSOLVE')
345 inrl=INDEX(exchangecorrelationandgridline,'NONRELATIVISTIC')
346 iscl=INDEX(exchangecorrelationandgridline,'SCALARRELATIVISTIC')
347 idirac=INDEX(exchangecorrelationandgridline,'DIRACRELATIVISTIC')
348 ipnt=INDEX(exchangecorrelationandgridline,'POINT-NUCLEUS')
349 ifin=INDEX(exchangecorrelationandgridline,'FINITE-NUCLEUS')
350 ilogder=INDEX(exchangecorrelationandgridline,'LOGDERIVRANGE')
351 ihfpp=INDEX(exchangecorrelationandgridline,'HFPOSTPROCESS')
352 ilcex=INDEX(exchangecorrelationandgridline,'LOCALIZEDCOREEXCHANGE')
353 ifixz=INDEX(exchangecorrelationandgridline,'FIXED_ZERO')
354 igrid=max(ilin,ilog)  !This line may need attention....
355 irelat=max(inrl,iscl) !This line may need attention....
356
357!Treat simple logical variables
358 dataset%scalarrelativistic=(iscl>0.and.inrl==0)
359 dataset%diracrelativistic=(idirac>0.and.inrl==0)
360 dataset%finitenucleus=(ifin>0.and.ipnt==0)
361 dataset%BDsolve=(ibd>0)
362 dataset%HFpostprocess=(ihfpp>0)
363
364!Treat finite nucleus option
365 dataset%finitenucleusmodel=-1
366 IF (dataset%finitenucleus) THEN
367   READ(exchangecorrelationandgridline(ifin+14:ifin+14),'(a)') CHR
368   IF (CHR=="2") dataset%finitenucleusmodel=2
369   IF (CHR=="3") dataset%finitenucleusmodel=3
370   IF (CHR=="4") dataset%finitenucleusmodel=4
371   IF (CHR=="5") dataset%finitenucleusmodel=5
372 END IF
373
374!Treat grid data
375 dataset%gridkey='LINEAR'
376 dataset%gridpoints=mxgridlin
377 dataset%gridrange=linrange
378 dataset%gridmatch=linrange
379 IF (ilog>0.and.ilin==0.and.ilogv4==0) THEN
380   dataset%gridkey='LOGGRID'
381   dataset%gridpoints=mxgridlog;
382   dataset%gridrange=logrange
383   dataset%gridmatch=logrange
384 END IF
385 IF (ilog>0.and.ilin==0.and.ilogv4>0) THEN
386   dataset%gridkey='LOGGRID4'
387   dataset%gridpoints=mxgridlog;
388   dataset%gridrange=v4logrange
389   dataset%gridmatch=v4logrange
390 END IF
391 IF (igrid>0) THEN
392   iend=128
393   IF (irelat >igrid.and.irelat-1 <iend) iend=irelat -1
394   IF (ilogder>igrid.and.ilogder-1<iend) iend=ilogder-1
395   IF (ibd>igrid.and.ibd-1<iend) iend=ibd-1
396   inputline=""
397   IF (ilog>0.and.ilogv4==0.and.iend>igrid+7) &
398&    inputline=TRIM(exchangecorrelationandgridline(igrid+7:iend))
399   IF (ilog>0.and.ilogv4>0.and.iend>igrid+9) &
400&    inputline=TRIM(exchangecorrelationandgridline(igrid+9:iend))
401   IF (ilin>0.and.iend>igrid+10) &
402&    inputline=TRIM(exchangecorrelationandgridline(igrid+10:iend))
403   IF (inputline/="") THEN
404     CALL extractword(1,inputline,inputword);inputword=trim(inputword)
405     IF (inputword/="") THEN
406       READ(inputword,*) dataset%gridpoints
407       CALL extractword(2,inputline,inputword);inputword=trim(inputword)
408       IF (inputword/="") THEN
409         READ(inputword,*) dataset%gridrange
410         dataset%gridmatch=dataset%gridrange
411         CALL extractword(3,inputline,inputword);inputword=trim(inputword)
412         IF (inputword/="") read(inputword,*) dataset%gridmatch
413       END IF
414     END IF
415   END IF
416   IF (dataset%gridpoints<=0) STOP "input_dataset: error -- number of grid points should be >0!"
417 END IF
418
419!Treat logderiv data
420 dataset%minlogderiv=logder_min
421 dataset%maxlogderiv=logder_max
422 dataset%nlogderiv=logder_pts
423 IF (ilogder>0) THEN
424   iend=128
425   IF (igrid >ilogder.and.igrid-1 <iend) iend=igrid -1
426   IF (irelat>ilogder.and.irelat-1<iend) iend=irelat-1
427   inputline=""
428   IF (iend>ilogder+13) inputline=trim(exchangecorrelationandgridline(ilogder+13:iend))
429   IF (inputline/="") THEN
430     CALL extractword(1,inputline,inputword);inputword=trim(inputword)
431     IF (inputword/="") THEN
432       READ(inputword,*) dataset%minlogderiv
433       CALL extractword(2,inputline,inputword);inputword=trim(inputword)
434       IF (inputword/="") THEN
435         READ(inputword,*) dataset%maxlogderiv
436         CALL extractword(3,inputline,inputword);inputword=trim(inputword)
437         IF (inputword/="") READ(inputword,*) dataset%nlogderiv
438       END IF
439     END IF
440   END IF
441 END IF
442
443!Treat XC/HF
444 READ(unit=exchangecorrelationandgridline,fmt=*) dataset%exctype
445 dataset%localizedcoreexchange=(ilcex>0)
446 dataset%fixed_zero=(ifixz>0) ; dataset%fixed_zero_index=-1
447 IF (dataset%fixed_zero) &
448&   READ(unit=exchangecorrelationandgridline(ifixz+10:),fmt=*) dataset%fixed_zero_index
449
450!Print read data
451 IF (has_to_print) THEN
452   WRITE(STD_OUT,'(3x,2a)')     "Scalar-relativistic calculation : ",MERGE("YES"," NO",dataset%scalarrelativistic)
453   WRITE(STD_OUT,'(3x,2a)')     "Dirac-relativistic calculation  : ",MERGE("YES"," NO",dataset%diracrelativistic)
454   WRITE(STD_OUT,'(3x,2a)')     "Finite-nucleus calculation      : ",MERGE("YES"," NO",dataset%finitenucleus)
455   IF (dataset%finitenucleus) THEN
456     WRITE(STD_OUT,'(3x,a,i0)') "    - Finite-nucleus model      : ",dataset%finitenucleusmodel
457   END IF
458   WRITE(STD_OUT,'(3x,2a)')     "Block-Davidson calculation      : ",MERGE("YES"," NO",dataset%BDsolve)
459   WRITE(STD_OUT,'(3x,2a)')     "Grid type                       : ",TRIM(dataset%gridkey)
460   WRITE(STD_OUT,'(3x,a,i0)')   "Grid size                       : ",dataset%gridpoints
461   WRITE(STD_OUT,'(3x,a,f7.3)') "Grid maximum value              : ",dataset%gridrange
462   WRITE(STD_OUT,'(3x,a,f7.3)') "Grid imposed value              : ",dataset%gridmatch
463   WRITE(STD_OUT,'(3x,a,i0)')   "Log. derivative, number of pts  : ",dataset%nlogderiv
464   WRITE(STD_OUT,'(3x,a,f7.3)') "Log. derivative, min. energy    : ",dataset%minlogderiv
465   WRITE(STD_OUT,'(3x,a,f7.3)') "Log. derivative, max. energy    : ",dataset%maxlogderiv
466   WRITE(STD_OUT,'(3x,2a)')     "Hartree-Fock, post-processing   : ",MERGE("YES"," NO",dataset%HFpostprocess)
467   WRITE(STD_OUT,'(3x,2a)')     "Hartree-Fock, localized core ex.: ",MERGE("YES"," NO",dataset%localizedcoreexchange)
468   WRITE(STD_OUT,'(3x,2a)')     "Hartree-Fock, fixed zero        : ",MERGE("YES"," NO",dataset%fixed_zero)
469   IF (dataset%fixed_zero) THEN
470     WRITE(STD_OUT,'(3x,a,i0)') "    - HF fixed zero index       : ",dataset%fixed_zero_index
471   END IF
472   IF (dataset%BDsolve.and.dataset%gridkey=='LINEAR') THEN
473     WRITE(STD_OUT,'(/,3x,a)') "WARNING: BlockDavidson solver works very slowly with linear grid!"
474   END IF
475END IF
476
477
478!------------------------------------------------------------------
479!End reading of global data. Start reading of electronic configuration data
480 ENDIF
481 IF (read_elec_data_) THEN
482
483
484!------------------------------------------------------------------
485!=== 3rd line and following: electronic configuration of atom
486
487 IF (has_to_ask) THEN
488   WRITE(STD_OUT,*) 'Enter maximum principal quantum numbers for s,p,d,f,g'
489 END IF
490
491 READ(input_unit,'(a)') inputline
492 IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
493 CALL eliminate_comment(inputline)
494
495 READ(inputline,*) dataset%np(1:5)
496
497 DO ll=1,5
498   IF(dataset%np(ll)<0) dataset%np(ll)=0
499 END DO
500 dataset%norbit=dataset%np(1)+max(dataset%np(2)-1,0)+max(dataset%np(3)-2,0) &
501&                            +max(dataset%np(4)-3,0)+max(dataset%np(5)-4,0)
502 IF (dataset%diracrelativistic) dataset%norbit=2*dataset%norbit-dataset%np(1)
503
504!Print read data
505 IF (has_to_print) THEN
506   WRITE(STD_OUT,'(3x,a,5(1x,i0))') "Max. quantum numbers (s,p,d,f,g): ",dataset%np(1:5)
507   WRITE(STD_OUT,'(3x,a,i0)') "Total number of orbitals: ",dataset%norbit
508 END IF
509
510 CALL input_dataset_read_occ(dataset%norbit_mod,dataset%orbit_mod_l,&
511&                   dataset%orbit_mod_n,dataset%orbit_mod_k,dataset%orbit_mod_occ,&
512&                   dataset%np,dataset%diracrelativistic,&
513&                   inputfile_unit=input_unit,echofile_unit=ecunit)
514
515
516!------------------------------------------------------------------
517!End reading of electronic data. Start reading of core/valence data
518 ENDIF
519 IF (read_coreval_data_) THEN
520
521
522!------------------------------------------------------------------
523!=== Core and valence states
524
525 IF (has_to_ask) THEN
526   WRITE(STD_OUT,*) 'For each state enter c for core or v for valence'
527 END IF
528
529!Read core and valence states
530 IF (ALLOCATED(dataset%orbit_iscore)) DEALLOCATE(dataset%orbit_iscore)
531 ALLOCATE(dataset%orbit_iscore(dataset%norbit))
532 DO io=1,dataset%norbit
533   DO
534     READ(input_unit,'(a)') inputline
535     CALL eliminate_comment(inputline)
536     READ(inputline,*) CHR
537     IF (CHR=='c'.OR.CHR=='C'.OR.&
538&        CHR=='v'.OR.CHR=='V') THEN
539       IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
540       EXIT
541     ELSE
542       WRITE(STD_OUT,*) ' >> Please input c or v!'
543     END IF
544   END DO
545   dataset%orbit_iscore(io)=(CHR=='c'.OR.CHR=='C')
546 END DO
547
548!Store valence states
549 dataset%norbit_val=dataset%norbit-COUNT(dataset%orbit_iscore(:))
550 IF (ALLOCATED(dataset%orbit_val_n)) DEALLOCATE(dataset%orbit_val_n)
551 IF (ALLOCATED(dataset%orbit_val_l)) DEALLOCATE(dataset%orbit_val_l)
552 IF (ALLOCATED(dataset%orbit_val_k)) DEALLOCATE(dataset%orbit_val_k)
553 ALLOCATE(dataset%orbit_val_n(dataset%norbit_val))
554 ALLOCATE(dataset%orbit_val_l(dataset%norbit_val))
555 ALLOCATE(dataset%orbit_val_k(dataset%norbit_val))
556 kk=0
557 io=0;nval=0
558 DO ll=0,4
559   nn=dataset%np(ll+1)
560   IF (nn>0) THEN
561     DO ik=1,MERGE(nkappa(ll+1),1,dataset%diracrelativistic)
562       kk=MERGE(ll,-(ll+1),ik==1);IF (ll==0) kk=-1
563       IF (.NOT.dataset%diracrelativistic) kk=0
564       DO ii=1+ll,nn
565         io=io+1
566         IF (.NOT.dataset%orbit_iscore(io)) THEN
567           nval=nval+1
568           dataset%orbit_val_n(nval)=ii
569           dataset%orbit_val_l(nval)=ll
570           dataset%orbit_val_k(nval)=kk
571         END IF
572       END DO
573     END DO
574   END IF
575 END DO
576 IF (dataset%norbit_val/=nval) STOP 'input_dataset: bug -- wrong nval!'
577
578!Print read data
579 IF (has_to_print) THEN
580   WRITE(STD_OUT,'(3x,a)') "Core and valence orbitals:"
581   IF (.NOT.dataset%diracrelativistic) WRITE(STD_OUT,'(7x,a)') "n l : type"
582   IF (dataset%diracrelativistic)      WRITE(STD_OUT,'(7x,a)') "n l kappa : type"
583   io=0
584   DO ll=0,4
585     nn=dataset%np(ll+1)
586     IF (nn>0) THEN
587       IF (.NOT.dataset%diracrelativistic) THEN
588         DO ii=1+ll,nn
589           io=io+1
590           WRITE(STD_OUT,'(7x,i1,1x,i1,2a)') ii,ll," : ", &
591&            MERGE("CORE   ","VALENCE",dataset%orbit_iscore(io))
592         END DO
593       ELSE
594         DO ik=1,nkappa(ll+1)
595           kk=MERGE(ll,-(ll+1),ik==1);IF (ll==0) kk=-1
596           DO ii=1+ll,nn
597             io=io+1
598             WRITE(STD_OUT,'(7x,i1,1x,i1,2x,i2,2x,2a)') ii,ll,kk," : ", &
599&              MERGE("CORE   ","VALENCE",dataset%orbit_iscore(io))
600           END DO
601         END DO
602       END IF
603     END IF
604   END DO
605 END IF
606
607
608!------------------------------------------------------------------
609!End reading of AE data. Start reading of basis data
610 ENDIF
611 IF (read_basis_data_) THEN
612
613
614!------------------------------------------------------------------
615!=== Maximum L for basis functions
616
617 IF (has_to_ask) THEN
618   WRITE(STD_OUT,*) 'Enter maximum L for basis and projector functions'
619 END IF
620
621 READ(STD_IN,'(a)') inputline
622 IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
623 CALL eliminate_comment(inputline)
624
625 READ(inputline,*) dataset%lmax
626
627!Print read data
628 IF (has_to_print) THEN
629   WRITE(STD_OUT,'(3x,a,i0)') "Basis, maximum L : ",dataset%lmax
630 END IF
631
632
633!------------------------------------------------------------------
634!=== Cut-off radii
635
636 IF (has_to_ask) THEN
637   WRITE(STD_OUT,*) 'Enter rc [and possibly: rc_shape, rc_vloc, rc_core]'
638 END IF
639
640 READ(STD_IN,'(a)') inputline
641 IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
642 CALL eliminate_comment(inputline)
643
644 CALL extractword(1,inputline,inputword);inputword=trim(inputword)
645 IF (inputword/="") READ(inputword,*) dataset%rc
646 IF (dataset%rc<=1.d-12) THEN
647   WRITE(STD_OUT,*) 'input_dataset: error -- rc too small ',dataset%rc,'!'
648   STOP
649 END IF
650
651 dataset%rc_shap=dataset%rc
652 dataset%rc_vloc=dataset%rc
653 dataset%rc_core=dataset%rc
654
655 CALL extractword(2,inputline,inputword);inputword=trim(inputword)
656 IF (inputword/="") THEN
657   READ(inputword,*) dataset%rc_shap
658   CALL extractword(3,inputline,inputword);inputword=trim(inputword)
659   IF (inputword/="") THEN
660     READ(inputword,*) dataset%rc_vloc
661     CALL extractword(4,inputline,inputword);inputword=trim(inputword)
662     IF (inputword/="") THEN
663       READ(inputword,*) dataset%rc_core
664     ELSE
665       WRITE(STD_OUT,*) 'input_dataset: error -- rc(core) is missing!'
666       STOP
667     END IF
668   ELSE
669     WRITE(STD_OUT,*) 'input_dataset: error -- rc(Vloc) is missing!'
670     STOP
671   END IF
672   IF (dataset%rc_shap<=1.d-12.OR.dataset%rc_vloc<=1.d-12.OR.&
673&      dataset%rc_core<=1.d-12) THEN
674     WRITE(STD_OUT,*) 'input_dataset: error -- one rc is too small!'
675     STOP
676   END IF
677   IF (dataset%rc_shap>dataset%rc.OR.dataset%rc_vloc>dataset%rc.OR.&
678&      dataset%rc_core>dataset%rc) THEN
679     WRITE(STD_OUT,*) 'input_dataset: error -- rc_shape, rc_vloc and rc_core must be <rc!'
680     STOP
681   END IF
682 ENDIF
683
684!Print read data
685 IF (has_to_print) THEN
686   WRITE(STD_OUT,'(3x,a,f7.4)') "Augmentation region radius : ",dataset%rc
687   WRITE(STD_OUT,'(3x,a,f7.4)') "Core dens. matching radius : ",dataset%rc_core
688   WRITE(STD_OUT,'(3x,a,f7.4)') "Local pot. matching radius : ",dataset%rc_vloc
689   WRITE(STD_OUT,'(3x,a,f7.4)') "Compens. shape func radius : ",dataset%rc_shap
690 END IF
691
692
693!------------------------------------------------------------------
694!=== Additional basis functions
695
696 nstart=0 ; dataset%nbasis_add=0 ; basis_add_k(:)=0
697 DO ll=0,dataset%lmax
698   nbl=0
699   nadd = MERGE(nkappa(ll+1),1,dataset%diracrelativistic)
700   IF (dataset%np(ll+1)>0) THEN
701     nbl=COUNT(.NOT.dataset%orbit_iscore(nstart+1:nstart+dataset%np(ll+1)-ll))
702     nstart=nstart+dataset%np(ll+1)-ll
703   END IF
704   DO
705     IF (has_to_ask) THEN
706       WRITE(STD_OUT,*) 'For l = ',ll,' there are currently ',nbl,' basis functions'
707       IF (.NOT.dataset%diracrelativistic) THEN
708         WRITE(STD_OUT,*) 'Enter y to add an additional function or n to go to next l'
709       ELSE
710         WRITE(STD_OUT,*) 'Enter y to add additional functions (one per kappa) or n to go to next l'
711       END IF
712     END IF
713     READ(STD_IN,'(a)') inputline
714     IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
715     CALL eliminate_comment(inputline)
716     READ(inputline,*) CHR
717     IF (CHR/='y'.AND.CHR/='Y') THEN
718       IF (CHR/='n'.AND.CHR/='N') STOP 'input_dataset: error -- Please enter Y or N!'
719       EXIT
720     END IF
721     dataset%nbasis_add=dataset%nbasis_add+nadd
722     IF (dataset%nbasis_add>nbasis_add_max) STOP 'Too many additional basis functions!'
723     basis_add_l(dataset%nbasis_add-nadd+1:dataset%nbasis_add)=ll
724     IF (dataset%diracrelativistic) THEN
725       basis_add_k(dataset%nbasis_add)=-1
726       IF (ll/=0) THEN
727         basis_add_k(dataset%nbasis_add-1)=ll
728         basis_add_k(dataset%nbasis_add)=-(ll+1)
729       END IF
730     END IF
731     IF (has_to_ask) THEN
732       IF (.NOT.dataset%diracrelativistic) THEN
733         WRITE(STD_OUT,*) 'Enter energy for generalized function'
734       ELSE
735         IF (ll==0) WRITE(STD_OUT,*) 'Enter 1 energy for generalized function (kappa=-1)'
736         IF (ll/=0) WRITE(STD_OUT,*) 'Enter 2 energies for generalized functions (one energy per kappa)'
737       END IF
738     END IF
739     READ(STD_IN,'(a)') inputline
740     IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
741     CALL eliminate_comment(inputline)
742     READ(inputline,*) basis_add_energy(dataset%nbasis_add-nadd+1:dataset%nbasis_add)
743   END DO
744 END DO
745
746 IF (ALLOCATED(dataset%basis_add_l)) DEALLOCATE(dataset%basis_add_l)
747 IF (ALLOCATED(dataset%basis_add_k)) DEALLOCATE(dataset%basis_add_k)
748 IF (ALLOCATED(dataset%basis_add_energy)) DEALLOCATE(dataset%basis_add_energy)
749 ALLOCATE(dataset%basis_add_l(dataset%nbasis_add))
750 ALLOCATE(dataset%basis_add_k(dataset%nbasis_add))
751 ALLOCATE(dataset%basis_add_energy(dataset%nbasis_add))
752 IF (dataset%nbasis_add>0) THEN
753   dataset%basis_add_l(1:dataset%nbasis_add)=basis_add_l(1:dataset%nbasis_add)
754   dataset%basis_add_k(1:dataset%nbasis_add)=basis_add_k(1:dataset%nbasis_add)
755   dataset%basis_add_energy(1:dataset%nbasis_add)=basis_add_energy(1:dataset%nbasis_add)
756 END IF
757
758 dataset%nbasis=COUNT(.NOT.dataset%orbit_iscore(:))+dataset%nbasis_add
759
760!Print read data
761 IF (has_to_print) THEN
762   WRITE(STD_OUT,'(3x,a,i0)') "Initial number of basis functions    : ",dataset%nbasis-dataset%nbasis_add
763   WRITE(STD_OUT,'(3x,a,i0)') "Number of additional basis functions : ",dataset%nbasis_add
764   WRITE(STD_OUT,'(3x,a,i0)') "Total number of basis functions      : ",dataset%nbasis
765   WRITE(STD_OUT,'(3x,a)') "Additional basis functions:"
766   IF (.NOT.dataset%diracrelativistic) THEN
767     WRITE(STD_OUT,'(7x,a)') "l : energy"
768     DO io=1,dataset%nbasis_add
769       WRITE(STD_OUT,'(7x,i1,a,f7.4)') dataset%basis_add_l(io)," : " ,dataset%basis_add_energy(io)
770     END DO
771   ELSE
772     WRITE(STD_OUT,'(7x,a)') "l kappa : energy"
773     DO io=1,dataset%nbasis_add
774       WRITE(STD_OUT,'(7x,i1,2x,i2,2x,a,f7.4)') dataset%basis_add_l(io), &
775&          dataset%basis_add_k(io)," : " ,dataset%basis_add_energy(io)
776     END DO
777   END IF
778 END IF
779
780
781!------------------------------------------------------------------
782!=== Projectors, compensation charge shape function, core tolerance
783
784 IF (has_to_ask) THEN
785   WRITE(STD_OUT,*) 'Enter "Bloechl", "Vanderbilt", "modrrkj" or "custom" keywords',&
786&             ' for projector generation method.'
787   WRITE(STD_OUT,*) ' In case of "custom" choice, enter additional (optional) keywords:'
788   WRITE(STD_OUT,*) ' - For partial waves pseudization scheme:'
789   WRITE(STD_OUT,*) '     "bloechlps", "polynom", "polynom2 p qcut" or "RRKJ"'
790   WRITE(STD_OUT,*) ' - For orthogonalization scheme:'
791   WRITE(STD_OUT,*) '     "gramschmidtortho" or "vanderbiltortho" or "svdortho"'
792   WRITE(STD_OUT,*) ' - Sinc^2 shape is set by default,'
793   WRITE(STD_OUT,*) '   Gaussian shape can be specified by adding "Gaussian" keyword and tol,'
794   WRITE(STD_OUT,*) '   Bessel shape can be specified by adding "Besselshape" keyword'
795 END IF
796
797 READ(STD_IN,'(a)') inputline
798 IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
799 CALL eliminate_comment(inputline)
800
801 CALL Uppercase(inputline)
802 inputline=TRIM(inputline)
803
804 dataset%pseudo_type=PSEUDO_TYPE_BLOECHL
805 dataset%ortho_type=ORTHO_TYPE_GRAMSCHMIDT
806 dataset%pseudo_polynom2_pdeg=polynom2_pdeg_def
807 dataset%pseudo_polynom2_qcut=polynom2_qcut_def
808 dataset%shapefunc_type=SHAPEFUNC_TYPE_SINC
809 dataset%shapefunc_gaussian_param=gausstol_def
810 dataset%hf_coretol=hf_coretol_def
811
812 READ(unit=inputline,fmt=*) inputword
813
814 IF (TRIM(inputword)=='BLOECHL'.OR.TRIM(inputword)=='VNCT') THEN
815   dataset%projector_type=PROJECTOR_TYPE_BLOECHL
816   dataset%pseudo_type=PSEUDO_TYPE_BLOECHL
817   dataset%ortho_type=ORTHO_TYPE_GRAMSCHMIDT
818 ELSE IF (TRIM(inputword)=='VNCK') THEN
819   dataset%projector_type=PROJECTOR_TYPE_BLOECHL
820   dataset%pseudo_type=PSEUDO_TYPE_BLOECHL_K
821   dataset%ortho_type=ORTHO_TYPE_GRAMSCHMIDT
822 ELSE IF (TRIM(inputword)=='VANDERBILT'.OR.TRIM(inputword)=='VNCTV') THEN
823   dataset%projector_type=PROJECTOR_TYPE_VANDERBILT
824   dataset%pseudo_type=PSEUDO_TYPE_POLYNOM
825   dataset%ortho_type=ORTHO_TYPE_VANDERBILT
826 ELSE IF(TRIM(inputword)=='MODRRKJ') THEN
827   dataset%projector_type=PROJECTOR_TYPE_MODRRKJ
828   dataset%pseudo_type=PSEUDO_TYPE_RRKJ
829   dataset%ortho_type=ORTHO_TYPE_VANDERBILT
830   IF (INDEX(inputline,'VANDERBILTORTHO')>0) dataset%ortho_type=ORTHO_TYPE_VANDERBILT
831   IF (INDEX(inputline,'GRAMSCHMIDTORTHO')>0) dataset%ortho_type=ORTHO_TYPE_GRAMSCHMIDT
832   IF (INDEX(inputline,'SVDORTHO')>0) dataset%ortho_type=ORTHO_TYPE_SVD
833 ELSE IF (TRIM(inputword)=='CUSTOM') THEN
834   dataset%projector_type=PROJECTOR_TYPE_CUSTOM
835   IF (INDEX(inputline,'BLOECHLPS')>0) THEN
836     dataset%pseudo_type=PSEUDO_TYPE_BLOECHL
837     dataset%ortho_type=ORTHO_TYPE_GRAMSCHMIDT
838   ELSE IF (INDEX(inputline,'POLYNOM2')>0) THEN
839     dataset%pseudo_type=PSEUDO_TYPE_POLYNOM2
840     nstart=INDEX(inputline,'POLYNOM2')
841     READ(unit=inputline(nstart+8:),fmt=*,err=111,end=111,iostat=io) &
842&         dataset%pseudo_polynom2_pdeg,dataset%pseudo_polynom2_qcut
843111  CONTINUE
844   ELSE IF (INDEX(inputline,'POLYNOM')>0) THEN
845     dataset%pseudo_type=PSEUDO_TYPE_POLYNOM
846   ELSE IF (INDEX(inputline,'RRKJ')>0) THEN
847     dataset%pseudo_type=PSEUDO_TYPE_RRKJ
848   END IF
849   IF (INDEX(inputline,'VANDERBILTORTHO')>0) dataset%ortho_type=ORTHO_TYPE_VANDERBILT
850   IF (INDEX(inputline,'GRAMSCHMIDTORTHO')>0) dataset%ortho_type=ORTHO_TYPE_GRAMSCHMIDT
851 END IF
852
853 IF (TRIM(dataset%exctype)=='HF') THEN
854   dataset%projector_type=PROJECTOR_TYPE_HF
855   dataset%pseudo_type=PSEUDO_TYPE_HF
856   dataset%ortho_type=ORTHO_TYPE_HF
857   WRITE(STD_OUT,'(3x,a)') '>> You are using HF XC type: pseudo and orthogonalization line will be ignored!'
858 END IF
859
860 IF ((dataset%pseudo_type==PSEUDO_TYPE_BLOECHL.OR.dataset%pseudo_type==PSEUDO_TYPE_BLOECHL_K)&
861&   .AND.dataset%ortho_type==ORTHO_TYPE_VANDERBILT) STOP &
862&  'input_dataset: error -- Vanderbilt orthogonalization not compatible with Bloechls projector scheme!'
863
864 IF (INDEX(inputline,'SINC2')>0) THEN
865   dataset%shapefunc_type=SHAPEFUNC_TYPE_SINC
866 ELSE IF (INDEX(inputline,'GAUSSIAN')>0) THEN
867   dataset%shapefunc_type=SHAPEFUNC_TYPE_GAUSSIAN
868   nstart=INDEX(inputline,'GAUSSIAN')
869   READ(unit=inputline(nstart+8:),fmt=*,err=222,end=222,iostat=io) &
870&       dataset%shapefunc_gaussian_param
871222 CONTINUE
872 ELSE IF (INDEX(inputline,'BESSELSHAPE')>0) THEN
873   dataset%shapefunc_type=SHAPEFUNC_TYPE_BESSEL
874 END IF
875
876 nstart=INDEX(inputline,'CORETOL')
877 IF (nstart>0) THEN
878   READ(unit=inputline(nstart+7:),fmt=*) dataset%hf_coretol
879 END IF
880
881!Print read data
882 IF (has_to_print) THEN
883   WRITE(STD_OUT,'(3x,a)') "Projectors description:"
884   IF (dataset%projector_type==PROJECTOR_TYPE_BLOECHL) &
885&    WRITE(STD_OUT,'(7x,a)') "Type              : BLOECHL"
886   IF (dataset%projector_type==PROJECTOR_TYPE_VANDERBILT) &
887&    WRITE(STD_OUT,'(7x,a)') "Type              : VANDERBILT"
888   IF (dataset%projector_type==PROJECTOR_TYPE_MODRRKJ) &
889&    WRITE(STD_OUT,'(7x,a)') "Type              : MODRRKJ"
890   IF (dataset%projector_type==PROJECTOR_TYPE_CUSTOM) &
891&    WRITE(STD_OUT,'(7x,a)') "Type              : CUSTOM"
892   IF (dataset%projector_type==PROJECTOR_TYPE_HF) &
893&    WRITE(STD_OUT,'(7x,a)') "Type : HARTREE-FOCK"
894   IF (dataset%projector_type/=PROJECTOR_TYPE_HF) THEN
895     IF (dataset%pseudo_type==PSEUDO_TYPE_BLOECHL) &
896&      WRITE(STD_OUT,'(7x,a)') "Pseudization      : BLOECHL"
897     IF (dataset%pseudo_type==PSEUDO_TYPE_POLYNOM) &
898&      WRITE(STD_OUT,'(7x,a)') "Pseudization      : POLYNOM"
899     IF (dataset%pseudo_type==PSEUDO_TYPE_RRKJ) &
900&      WRITE(STD_OUT,'(7x,a)') "Pseudization      : RRKJ"
901     IF (dataset%pseudo_type==PSEUDO_TYPE_BLOECHL_K) &
902&      WRITE(STD_OUT,'(7x,a)') "Pseudization      : BLOECHL KERKER"
903     IF (dataset%pseudo_type==PSEUDO_TYPE_POLYNOM2) &
904&      WRITE(STD_OUT,'(7x,a,i0,a,es9.3)') "Pseudization      : POLYNOM2, pdeg=", &
905&       dataset%pseudo_polynom2_pdeg,", qcut=",dataset%pseudo_polynom2_qcut
906     IF (dataset%ortho_type==ORTHO_TYPE_GRAMSCHMIDT) &
907&      WRITE(STD_OUT,'(7x,a)') "Orthogonalisation : GRAM-SCHMIDT"
908     IF (dataset%ortho_type==ORTHO_TYPE_VANDERBILT) &
909&      WRITE(STD_OUT,'(7x,a)') "Orthogonalisation : VANDERBILT"
910     IF (dataset%ortho_type==ORTHO_TYPE_SVD) &
911&      WRITE(STD_OUT,'(7x,a)') "Orthogonalisation : SVD"
912   END IF
913   IF (dataset%shapefunc_type==SHAPEFUNC_TYPE_GAUSSIAN) &
914&    WRITE(STD_OUT,'(3x,a,es9.3)') "Compensation charge shape function : GAUSSIAN, tol=",&
915&    dataset%shapefunc_gaussian_param
916   IF (dataset%shapefunc_type==SHAPEFUNC_TYPE_SINC) &
917&    WRITE(STD_OUT,'(3x,a)') "Compensation charge shape function : SINC2"
918   IF (dataset%shapefunc_type==SHAPEFUNC_TYPE_BESSEL) &
919&    WRITE(STD_OUT,'(3x,a)') "Compensation charge shape function : BESSEL"
920   IF (INDEX(inputline,'CORETOL')>0) &
921&    WRITE(STD_OUT,'(3x,a,es9.3)') "Core tolerance for Hartree-Fock : ",dataset%hf_coretol
922 END IF
923
924
925!------------------------------------------------------------------
926!=== Local pseudopotential
927
928 IF (has_to_ask) THEN
929   WRITE(STD_OUT,*) 'To generate the local pseudopotential, this code can use:'
930   WRITE(STD_OUT,*) '  1- a Troullier-Martins scheme for specified l value and energy'
931   WRITE(STD_OUT,*) '  2- a non norm-conserving pseudopotential scheme for specified l value and energy'
932   WRITE(STD_OUT,*) '  3- a simple pseudization scheme using a single spherical Bessel function'
933   WRITE(STD_OUT,*) '  4- Vloc ==  VlocCoef*Shapefunc'
934   WRITE(STD_OUT,*) 'For choice 1, enter (high) l value and energy e'
935   WRITE(STD_OUT,*) 'For choice 2, enter (high) l value, energy e and "ultrasoft"'
936   WRITE(STD_OUT,*) 'For choice 3, enter "bessel"'
937   WRITE(STD_OUT,*) 'For choice 4, enter "setvloc x y" - x is VlocCoef y is VlocRad'
938 END IF
939
940 READ(STD_IN,'(a)') inputline
941 IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
942 CALL eliminate_comment(inputline)
943
944 call Uppercase(inputline)
945 inputline=TRIM(inputline)
946
947 dataset%vloc_type=VLOC_TYPE_MTROULLIER
948 dataset%vloc_l=-1
949 dataset%vloc_ene=0.d0
950 dataset%vloc_setvloc_coef=0.d0
951 dataset%vloc_setvloc_rad=dataset%rc
952 dataset%vloc_kerker_power(:)=0
953
954 IF (INDEX(inputline,'MTROULLIER')>0) THEN
955   dataset%vloc_type=VLOC_TYPE_MTROULLIER
956 ELSE IF (INDEX(inputline,'ULTRASOFT')>0) THEN
957   dataset%vloc_type=VLOC_TYPE_ULTRASOFT
958 ELSE IF (INDEX(inputline,'BESSEL')>0) THEN
959   dataset%vloc_type=VLOC_TYPE_BESSEL
960 ELSE IF (INDEX(inputline,'SETVLOC')>0) THEN
961   dataset%vloc_type=VLOC_TYPE_SETVLOC
962   nstart=INDEX(inputline,'SETVLOC')
963   READ(unit=inputline(nstart+8:),fmt=*,err=333,end=333,iostat=io) x1,x2
964   IF (x1<1.d3.AND.x1>-1.d3) dataset%vloc_setvloc_coef=x1
965   IF (x2>1.d-8.AND.x2<dataset%rc) dataset%vloc_setvloc_rad=x2
966333  CONTINUE
967 ELSE IF (INDEX(inputline,'KERKER')>0.OR.dataset%pseudo_type==PSEUDO_TYPE_BLOECHL_K) THEN
968   IF (INDEX(inputline,'EXPF')>0) THEN
969     dataset%vloc_type=VLOC_TYPE_KERKER_EXPF
970     nstart=INDEX(inputline,'EXPF')
971   ELSE IF (INDEX(inputline,'POLY')>0) THEN
972     dataset%vloc_type=VLOC_TYPE_KERKER_POLY
973     nstart=INDEX(inputline,'POLY')
974   ELSE
975     STOP "EXPF or POLY keyword missing!"
976   END IF
977   READ(unit=inputline(nstart+5:),fmt=*,err=334,end=334,iostat=io) &
978&    dataset%vloc_kerker_power(1:4)
979334  CONTINUE
980 END IF
981
982 IF (dataset%vloc_type==VLOC_TYPE_MTROULLIER.OR.dataset%vloc_type==VLOC_TYPE_ULTRASOFT) THEN
983   READ(unit=inputline,fmt=*,err=444,end=444,iostat=io) dataset%vloc_l,dataset%vloc_ene
984444  CONTINUE
985   IF (dataset%vloc_l<0.or.dataset%vloc_l>10) STOP 'input_dataset: error while reading Vloc parameters!'
986 END IF
987
988!Print read data
989 IF (has_to_print) THEN
990   IF (dataset%vloc_type==VLOC_TYPE_MTROULLIER) &
991&    WRITE(STD_OUT,'(7x,a,i0,a,f7.4)') "Local pseudopotential type : MTROULLIER, l=",&
992&          dataset%vloc_l,", energy=",dataset%vloc_ene
993   IF (dataset%vloc_type==VLOC_TYPE_ULTRASOFT) &
994&    WRITE(STD_OUT,'(7x,a,i0,a,f7.4)') "Local pseudopotential type : ULTRASOFT, l=",&
995&          dataset%vloc_l,", energy=",dataset%vloc_ene
996   IF (dataset%vloc_type==VLOC_TYPE_BESSEL) &
997&    WRITE(STD_OUT,'(7x,a)') "Local pseudopotential type : BESSEL"
998   IF (dataset%vloc_type==VLOC_TYPE_SETVLOC) &
999&    WRITE(STD_OUT,'(7x,a,es9.4,a,es9.4)') "Local pseudopotential type : SETVLOC, coef=",&
1000&          dataset%vloc_setvloc_coef,", rad=",dataset%vloc_setvloc_rad
1001   IF (dataset%vloc_type==VLOC_TYPE_KERKER_EXPF) &
1002&    WRITE(STD_OUT,'(7x,a,4(1x,i0))') "Local pseudopotential type : KERKER EXPF, powers=",&
1003&          dataset%vloc_kerker_power(1:4)
1004   IF (dataset%vloc_type==VLOC_TYPE_KERKER_POLY) &
1005&    WRITE(STD_OUT,'(7x,a,4(1x,i0))') "Local pseudopotential type : KERKER POLY, powers=",&
1006&          dataset%vloc_kerker_power(1:4)
1007 END IF
1008
1009
1010!------------------------------------------------------------------
1011!=== Matching radii for the basis functions
1012
1013!Not for all choice of projectors
1014 IF (dataset%projector_type==PROJECTOR_TYPE_CUSTOM.OR.&
1015&    dataset%projector_type==PROJECTOR_TYPE_VANDERBILT.OR.&
1016&    dataset%projector_type==PROJECTOR_TYPE_MODRRKJ.OR.&
1017&    dataset%projector_type==PROJECTOR_TYPE_HF.AND.dataset%vloc_type==VLOC_TYPE_MTROULLIER) THEN
1018
1019   IF (ALLOCATED(dataset%basis_func_rc)) DEALLOCATE(dataset%basis_func_rc)
1020   ALLOCATE(dataset%basis_func_rc(dataset%nbasis))
1021
1022   IF (has_to_ask) WRITE(STD_OUT,*) 'For each of the following basis functions enter rc'
1023
1024   norb=0
1025   DO ll=0,dataset%lmax
1026     DO ik=1,MERGE(nkappa(ll+1),1,dataset%diracrelativistic)
1027       kk=MERGE(ll,-(ll+1),ik==1);IF (ll==0) kk=-1
1028       IF (.NOT.dataset%diracrelativistic) kk=0
1029       DO io=1,dataset%norbit_val
1030         IF (dataset%orbit_val_l(io)==ll.AND. &
1031&           ((.NOT.dataset%diracrelativistic).OR.dataset%orbit_val_k(io)==kk)) THEN
1032           norb=norb+1
1033           IF (has_to_ask.AND.(.NOT.dataset%diracrelativistic)) &
1034&            WRITE(STD_OUT,'(a,i2,a,2i4)') '  rc for basis function ',norb,' - n,l= ',norb,ll
1035           IF (has_to_ask.AND.dataset%diracrelativistic) &
1036&            WRITE(STD_OUT,'(a,i2,a,3i4)') '  rc for basis function ',norb,' - n,l,kappa= ',norb,ll,kk
1037           READ(STD_IN,'(a)') inputline
1038           IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
1039           CALL eliminate_comment(inputline)
1040           READ(inputline,*) dataset%basis_func_rc(norb)
1041         END IF
1042       END DO
1043       IF (dataset%nbasis_add>0) THEN
1044         DO io=1,dataset%nbasis_add
1045           IF (dataset%basis_add_l(io)==ll.AND. &
1046&             ((.NOT.dataset%diracrelativistic).OR.dataset%basis_add_k(io)==kk)) THEN
1047             norb=norb+1
1048             IF (has_to_ask.AND.(.NOT.dataset%diracrelativistic)) &
1049&              WRITE(STD_OUT,'(a,i2,a,2i4)') '  rc for basis function ',norb,' - n,l= ',999,ll
1050             IF (has_to_ask.AND.dataset%diracrelativistic) &
1051&              WRITE(STD_OUT,'(a,i2,a,3i4)') '  rc for basis function ',norb,' - n,l,kappa= ',999,ll,kk
1052             READ(STD_IN,'(a)') inputline
1053             IF (has_to_echo) WRITE(ecunit,'(a)') TRIM(inputline)
1054             CALL eliminate_comment(inputline)
1055             READ(inputline,*) dataset%basis_func_rc(norb)
1056           END IF
1057         END DO
1058       END IF
1059     END DO
1060   END DO
1061
1062   IF (dataset%nbasis/=norb) STOP 'input_dataset: error -- inconsistency in the number of basis functions!'
1063
1064!  Print read data
1065   IF (has_to_print) THEN
1066     WRITE(STD_OUT,'(3x,a)') "Matching radius for basis functions:"
1067     IF (.NOT.dataset%diracrelativistic) WRITE(STD_OUT,'(7x,a)') " # - n l : radius"
1068     IF (dataset%diracrelativistic) WRITE(STD_OUT,'(7x,a)') " # - n l kappa : radius"
1069     norb=0
1070     DO ll=0,dataset%lmax
1071       DO ik=1,MERGE(nkappa(ll+1),1,dataset%diracrelativistic)
1072         kk=MERGE(ll,-(ll+1),ik==1);IF (ll==0) kk=-1
1073         IF (.NOT.dataset%diracrelativistic) kk=0
1074         DO io=1,dataset%norbit_val
1075           IF (dataset%orbit_val_l(io)==ll.AND. &
1076&             ((.NOT.dataset%diracrelativistic).OR.dataset%orbit_val_k(io)==kk)) THEN
1077             norb=norb+1
1078             IF (.NOT.dataset%diracrelativistic) &
1079&              WRITE(STD_OUT,'(7x,i2,a,i1,1x,i1,a,f7.4)') &
1080&              norb," - ",dataset%orbit_val_n(io),ll," : ",dataset%basis_func_rc(norb)
1081             IF (dataset%diracrelativistic) &
1082&              WRITE(STD_OUT,'(7x,i2,a,i1,1x,i1,2x,i2,2x,a,f7.4)') &
1083&              norb," - ",dataset%orbit_val_n(io),ll,kk," : ",dataset%basis_func_rc(norb)
1084           END IF
1085         END DO
1086         IF (dataset%nbasis_add>0) THEN
1087           DO io=1,dataset%nbasis_add
1088             IF (dataset%basis_add_l(io)==ll.AND. &
1089&             ((.NOT.dataset%diracrelativistic).OR.dataset%basis_add_k(io)==kk)) THEN
1090               norb=norb+1
1091               IF (.NOT.dataset%diracrelativistic) &
1092&                WRITE(STD_OUT,'(7x,i2,a,a1,1x,i1,a,f7.4)') &
1093&                norb," - ",".",ll," : ",dataset%basis_func_rc(norb)
1094               IF (dataset%diracrelativistic) &
1095&                WRITE(STD_OUT,'(7x,i2,a,a1,1x,i1,2x,i2,2x,a,f7.4)') &
1096&                norb," - ",".",ll,kk," : ",dataset%basis_func_rc(norb)
1097             END IF
1098           END DO
1099         END IF
1100       END DO
1101     END DO
1102   END IF
1103
1104 ELSE ! Other projectors
1105   IF (ALLOCATED(dataset%basis_func_rc)) DEALLOCATE(dataset%basis_func_rc)
1106   ALLOCATE(dataset%basis_func_rc(0))
1107 END IF
1108
1109
1110!------------------------------------------------------------------
1111!End reading of basis data
1112 ENDIF
1113
1114!Final message
1115 IF (read_global_data_.OR.read_elec_data_.OR.read_coreval_data_.OR.read_basis_data_) THEN
1116   IF (has_to_ask) THEN
1117     WRITE(STD_OUT,'(1x,a)') "===== END READING OF INPUT DATA ====="
1118   ELSE
1119     WRITE(STD_OUT,'(3x,a)') "===== END READING OF INPUT FILE ====="
1120   END IF
1121 END IF
1122 WRITE(STD_OUT,'(2/)')
1123
1124
1125!------------------------------------------------------------------
1126!Close files
1127 IF (PRESENT(inputfile)) THEN
1128   CLOSE(ifunit)
1129 END IF
1130 IF (has_to_echo) THEN
1131   CLOSE(ecunit)
1132 END IF
1133
1134 END SUBROUTINE input_dataset_read
1135
1136
1137!!=================================================================
1138!! NAME
1139!!  input_dataset_free
1140!!
1141!! FUNCTION
1142!!  Free and destroy a input_dataset datastructure
1143!!
1144!! SIDE EFFECT
1145!!  [input_dt]= datastructure to be destroyed (optional)
1146!!              If omitted, then the global public `input_dataset`
1147!!              is used.
1148!!
1149!!=================================================================
1150 SUBROUTINE input_dataset_free(input_dt)
1151
1152!---- Arguments
1153 TYPE(input_dataset_t),INTENT(inout),OPTIONAL,TARGET :: input_dt
1154
1155!---- Local variables
1156 TYPE(input_dataset_t),POINTER :: dataset
1157
1158!------------------------------------------------------------------
1159
1160!Select datastruture to be read
1161 IF (PRESENT(input_dt)) THEN
1162   dataset => input_dt
1163 ELSE
1164   dataset => input_dataset
1165 ENDIF
1166
1167!Integers
1168 dataset%atomic_charge            =-1
1169 dataset%finitenucleusmodel       =-1
1170 dataset%gridpoints               = 0
1171 dataset%nlogderiv                = 0
1172 dataset%Fixed_zero_index         = 0
1173 dataset%np(1:5)                  = 0
1174 dataset%norbit                   = 0
1175 dataset%norbit_mod               = 0
1176 dataset%norbit_val               = 0
1177 dataset%lmax                     =-1
1178 dataset%nbasis                   = 0
1179 dataset%nbasis_add               = 0
1180 dataset%projector_type           = UNKNOWN_TYPE
1181 dataset%pseudo_type              = UNKNOWN_TYPE
1182 dataset%ortho_type               = UNKNOWN_TYPE
1183 dataset%vloc_type                = UNKNOWN_TYPE
1184 dataset%vloc_l                   = 0
1185 dataset%pseudo_polynom2_pdeg     = 0
1186 dataset%shapefunc_type           = UNKNOWN_TYPE
1187 dataset%abinit_log_meshsz        = 0
1188 dataset%xml_spl_meshsz           = 0
1189 dataset%xml_lda12_orb_l          = -1
1190 dataset%xml_lda12_orb_n          = -1
1191 dataset%vloc_kerker_power(1:4)   = 0
1192
1193!Logicals
1194 dataset%scalarrelativistic       = .false.
1195 dataset%diracrelativistic        = .false.
1196 dataset%finitenucleus            = .false.
1197 dataset%HFpostprocess            = .false.
1198 dataset%localizedcoreexchange    = .false.
1199 dataset%BDsolve                  = .false.
1200 dataset%Fixed_zero               = .false.
1201 dataset%abinit_usexcnhat         = .false.
1202 dataset%abinit_prtcorewf         = .false.
1203 dataset%abinit_uselog            = .false.
1204 dataset%abinit_userso            = .false.
1205 dataset%xml_usexcnhat            = .false.
1206 dataset%xml_prtcorewf            = .false.
1207 dataset%xml_usespl               = .false.
1208 dataset%xml_userso               = .false.
1209 dataset%xml_uselda12             = .false.
1210
1211!Reals
1212 dataset%gridrange                = 0.d0
1213 dataset%gridmatch                = 0.d0
1214 dataset%minlogderiv              = 0.d0
1215 dataset%maxlogderiv              = 0.d0
1216 dataset%rc                       = 0.d0
1217 dataset%rc_shap                  = 0.d0
1218 dataset%rc_vloc                  = 0.d0
1219 dataset%rc_core                  = 0.d0
1220 dataset%pseudo_polynom2_qcut     = 0.d0
1221 dataset%shapefunc_gaussian_param = 0.d0
1222 dataset%hf_coretol               = 0.d0
1223 dataset%vloc_ene                 = 0.d0
1224 dataset%vloc_setvloc_coef        = 0.d0
1225 dataset%vloc_setvloc_rad         = 0.d0
1226 dataset%abinit_log_step          = 0.d0
1227 dataset%abinit_rso_ecut          = 0.d0
1228 dataset%abinit_rso_gfact         = 0.d0
1229 dataset%abinit_rso_werror        = 0.d0
1230 dataset%xml_rso_ecut             = 0.d0
1231 dataset%xml_rso_gfact            = 0.d0
1232 dataset%xml_rso_werror           = 0.d0
1233 dataset%xml_lda12_ion            = 0.d0
1234 dataset%xml_lda12_rcut           = 0.d0
1235 dataset%upf_grid_xmin            = 0.d0
1236 dataset%upf_grid_zmesh           = 0.d0
1237 dataset%upf_grid_dx              = 0.d0
1238 dataset%upf_grid_range           = 0.d0
1239
1240!Character strings
1241 dataset%atomic_symbol            = ''
1242 dataset%gridkey                  = ''
1243 dataset%exctype                  = ''
1244 dataset%abinit_author            = ''
1245 dataset%xml_author               = ''
1246 dataset%xml_comment              = ''
1247 dataset%xml_lda12_logfile        = ''
1248
1249!Integer allocatable arrays
1250 IF (ALLOCATED(dataset%orbit_mod_l)) DEALLOCATE(dataset%orbit_mod_l)
1251 IF (ALLOCATED(dataset%orbit_mod_n)) DEALLOCATE(dataset%orbit_mod_n)
1252 IF (ALLOCATED(dataset%orbit_mod_k)) DEALLOCATE(dataset%orbit_mod_k)
1253 IF (ALLOCATED(dataset%orbit_val_l)) DEALLOCATE(dataset%orbit_val_l)
1254 IF (ALLOCATED(dataset%orbit_val_n)) DEALLOCATE(dataset%orbit_val_n)
1255 IF (ALLOCATED(dataset%orbit_val_k)) DEALLOCATE(dataset%orbit_val_k)
1256 IF (ALLOCATED(dataset%basis_add_l)) DEALLOCATE(dataset%basis_add_l)
1257 IF (ALLOCATED(dataset%basis_add_k)) DEALLOCATE(dataset%basis_add_k)
1258
1259!Logical allocatable arrays
1260 IF (ALLOCATED(dataset%orbit_iscore)) DEALLOCATE(dataset%orbit_iscore)
1261
1262!Real allocatable arrays
1263 IF (ALLOCATED(dataset%orbit_mod_occ))    DEALLOCATE(dataset%orbit_mod_occ)
1264 IF (ALLOCATED(dataset%basis_add_energy)) DEALLOCATE(dataset%basis_add_energy)
1265 IF (ALLOCATED(dataset%basis_func_rc))    DEALLOCATE(dataset%basis_func_rc)
1266
1267 END SUBROUTINE input_dataset_free
1268
1269
1270!!=================================================================
1271!! NAME
1272!!  input_dataset_copy
1273!!
1274!! FUNCTION
1275!!  Copy a input_dataset datastructure into another
1276!!
1277!! INPUTS
1278!!  input_dt= datastructure to copy
1279!!
1280!! OUTPUT
1281!!  output_dt= copied datastructure
1282!!
1283!!=================================================================
1284 SUBROUTINE input_dataset_copy(input_dt,output_dt)
1285
1286!---- Arguments
1287 TYPE(input_dataset_t),INTENT(in) :: input_dt
1288 TYPE(input_dataset_t),INTENT(inout) :: output_dt
1289
1290!------------------------------------------------------------------
1291
1292!Integers
1293 output_dt%atomic_charge           =input_dt%atomic_charge
1294 output_dt%finitenucleusmodel      =input_dt%finitenucleusmodel
1295 output_dt%gridpoints              =input_dt%gridpoints
1296 output_dt%nlogderiv               =input_dt%nlogderiv
1297 output_dt%Fixed_zero_index        =input_dt%Fixed_zero_index
1298 output_dt%np(1:5)                 =input_dt%np(1:5)
1299 output_dt%norbit                  =input_dt%norbit
1300 output_dt%norbit_mod              =input_dt%norbit_mod
1301 output_dt%norbit_val              =input_dt%norbit_val
1302 output_dt%lmax                    =input_dt%lmax
1303 output_dt%nbasis                  =input_dt%nbasis
1304 output_dt%nbasis_add              =input_dt%nbasis_add
1305 output_dt%projector_type          =input_dt%projector_type
1306 output_dt%pseudo_type             =input_dt%pseudo_type
1307 output_dt%ortho_type              =input_dt%ortho_type
1308 output_dt%pseudo_polynom2_pdeg    =input_dt%pseudo_polynom2_pdeg
1309 output_dt%shapefunc_type          =input_dt%shapefunc_type
1310 output_dt%vloc_type               =input_dt%vloc_type
1311 output_dt%vloc_l                  =input_dt%vloc_l
1312 output_dt%abinit_log_meshsz       =input_dt%abinit_log_meshsz
1313 output_dt%xml_spl_meshsz          =input_dt%xml_spl_meshsz
1314 output_dt%xml_lda12_orb_l         =input_dt%xml_lda12_orb_l
1315 output_dt%xml_lda12_orb_n         =input_dt%xml_lda12_orb_n
1316 output_dt%vloc_kerker_power(1:4)  =input_dt%vloc_kerker_power(1:4)
1317
1318!Logicals
1319 output_dt%scalarrelativistic      =input_dt%scalarrelativistic
1320 output_dt%diracrelativistic       =input_dt%diracrelativistic
1321 output_dt%finitenucleus           =input_dt%finitenucleus
1322 output_dt%HFpostprocess           =input_dt%HFpostprocess
1323 output_dt%localizedcoreexchange   =input_dt%localizedcoreexchange
1324 output_dt%BDsolve                 =input_dt%BDsolve
1325 output_dt%Fixed_zero              =input_dt%Fixed_zero
1326 output_dt%abinit_usexcnhat        =input_dt%abinit_usexcnhat
1327 output_dt%abinit_prtcorewf        =input_dt%abinit_prtcorewf
1328 output_dt%abinit_uselog           =input_dt%abinit_uselog
1329 output_dt%abinit_userso           =input_dt%abinit_userso
1330 output_dt%xml_usexcnhat           =input_dt%xml_usexcnhat
1331 output_dt%xml_prtcorewf           =input_dt%xml_prtcorewf
1332 output_dt%xml_usespl              =input_dt%xml_usespl
1333 output_dt%xml_userso              =input_dt%xml_userso
1334 output_dt%xml_uselda12            =input_dt%xml_uselda12
1335
1336!Reals
1337 output_dt%gridrange               =input_dt%gridrange
1338 output_dt%gridmatch               =input_dt%gridmatch
1339 output_dt%minlogderiv             =input_dt%minlogderiv
1340 output_dt%maxlogderiv             =input_dt%maxlogderiv
1341 output_dt%rc                      =input_dt%rc
1342 output_dt%rc_shap                 =input_dt%rc_shap
1343 output_dt%rc_vloc                 =input_dt%rc_vloc
1344 output_dt%rc_core                 =input_dt%rc_core
1345 output_dt%pseudo_polynom2_qcut    =input_dt%pseudo_polynom2_qcut
1346 output_dt%shapefunc_gaussian_param=input_dt%shapefunc_gaussian_param
1347 output_dt%hf_coretol              =input_dt%hf_coretol
1348 output_dt%vloc_ene                =input_dt%vloc_ene
1349 output_dt%vloc_setvloc_coef       =input_dt%vloc_setvloc_coef
1350 output_dt%vloc_setvloc_rad        =input_dt%vloc_setvloc_rad
1351 output_dt%abinit_log_step         =input_dt%abinit_log_step
1352 output_dt%abinit_rso_ecut         =input_dt%abinit_rso_ecut
1353 output_dt%abinit_rso_gfact        =input_dt%abinit_rso_gfact
1354 output_dt%abinit_rso_werror       =input_dt%abinit_rso_werror
1355 output_dt%xml_rso_ecut            =input_dt%xml_rso_ecut
1356 output_dt%xml_rso_gfact           =input_dt%xml_rso_gfact
1357 output_dt%xml_rso_werror          =input_dt%xml_rso_werror
1358 output_dt%xml_lda12_ion           =input_dt%xml_lda12_ion
1359 output_dt%xml_lda12_rcut          =input_dt%xml_lda12_rcut
1360 output_dt%upf_grid_xmin           =input_dt%upf_grid_xmin
1361 output_dt%upf_grid_zmesh          =input_dt%upf_grid_zmesh
1362 output_dt%upf_grid_dx             =input_dt%upf_grid_dx
1363 output_dt%upf_grid_range          =input_dt%upf_grid_range
1364
1365!Character strings
1366 output_dt%atomic_symbol           =TRIM(input_dt%atomic_symbol)
1367 output_dt%gridkey                 =TRIM(input_dt%gridkey)
1368 output_dt%exctype                 =TRIM(input_dt%exctype)
1369 output_dt%abinit_author           =TRIM(input_dt%abinit_author)
1370 output_dt%xml_author              =TRIM(input_dt%xml_author)
1371 output_dt%xml_comment             =TRIM(input_dt%xml_comment)
1372 output_dt%xml_lda12_logfile       =TRIM(input_dt%xml_lda12_logfile)
1373
1374!Integer allocatable arrays
1375 IF (ALLOCATED(output_dt%orbit_mod_l)) DEALLOCATE(output_dt%orbit_mod_l)
1376 IF (ALLOCATED(output_dt%orbit_mod_n)) DEALLOCATE(output_dt%orbit_mod_n)
1377 IF (ALLOCATED(output_dt%orbit_mod_k)) DEALLOCATE(output_dt%orbit_mod_k)
1378 IF (ALLOCATED(output_dt%orbit_val_l)) DEALLOCATE(output_dt%orbit_val_l)
1379 IF (ALLOCATED(output_dt%orbit_val_n)) DEALLOCATE(output_dt%orbit_val_n)
1380 IF (ALLOCATED(output_dt%orbit_val_k)) DEALLOCATE(output_dt%orbit_val_k)
1381 IF (ALLOCATED(output_dt%basis_add_l)) DEALLOCATE(output_dt%basis_add_l)
1382 IF (ALLOCATED(output_dt%basis_add_k)) DEALLOCATE(output_dt%basis_add_k)
1383 ALLOCATE(output_dt%orbit_mod_l(input_dt%norbit_mod))
1384 ALLOCATE(output_dt%orbit_mod_n(input_dt%norbit_mod))
1385 ALLOCATE(output_dt%orbit_mod_k(input_dt%norbit_mod))
1386 ALLOCATE(output_dt%orbit_val_l(input_dt%norbit_val))
1387 ALLOCATE(output_dt%orbit_val_n(input_dt%norbit_val))
1388 ALLOCATE(output_dt%orbit_val_k(input_dt%norbit_val))
1389 ALLOCATE(output_dt%basis_add_l(input_dt%nbasis_add))
1390 ALLOCATE(output_dt%basis_add_k(input_dt%nbasis_add))
1391 output_dt%orbit_mod_l(1:input_dt%norbit_mod)=input_dt%orbit_mod_l(1:input_dt%norbit_mod)
1392 output_dt%orbit_mod_n(1:input_dt%norbit_mod)=input_dt%orbit_mod_n(1:input_dt%norbit_mod)
1393 output_dt%orbit_mod_k(1:input_dt%norbit_mod)=input_dt%orbit_mod_k(1:input_dt%norbit_mod)
1394 output_dt%orbit_val_l(1:input_dt%norbit_val)=input_dt%orbit_val_l(1:input_dt%norbit_val)
1395 output_dt%orbit_val_n(1:input_dt%norbit_val)=input_dt%orbit_val_n(1:input_dt%norbit_val)
1396 output_dt%orbit_val_k(1:input_dt%norbit_val)=input_dt%orbit_val_k(1:input_dt%norbit_val)
1397 output_dt%basis_add_l(1:input_dt%nbasis_add)=input_dt%basis_add_l(1:input_dt%nbasis_add)
1398 output_dt%basis_add_k(1:input_dt%nbasis_add)=input_dt%basis_add_k(1:input_dt%nbasis_add)
1399
1400!Logical allocatable arrays
1401 IF (ALLOCATED(output_dt%orbit_iscore)) DEALLOCATE(output_dt%orbit_iscore)
1402 ALLOCATE(output_dt%orbit_iscore(input_dt%norbit))
1403 output_dt%orbit_iscore(1:input_dt%norbit)=input_dt%orbit_iscore(1:input_dt%norbit)
1404
1405!Real allocatable arrays
1406 IF (ALLOCATED(output_dt%orbit_mod_occ))    DEALLOCATE(output_dt%orbit_mod_occ)
1407 IF (ALLOCATED(output_dt%basis_add_energy)) DEALLOCATE(output_dt%basis_add_energy)
1408 IF (ALLOCATED(output_dt%basis_func_rc))    DEALLOCATE(output_dt%basis_func_rc)
1409 ALLOCATE(output_dt%orbit_mod_occ(input_dt%norbit_mod))
1410 ALLOCATE(output_dt%basis_add_energy(input_dt%nbasis_add))
1411 ALLOCATE(output_dt%basis_func_rc(input_dt%nbasis))
1412 output_dt%orbit_mod_occ(1:input_dt%norbit_mod)=input_dt%orbit_mod_occ(1:input_dt%norbit_mod)
1413 output_dt%basis_add_energy(1:input_dt%nbasis_add)=input_dt%basis_add_energy(1:input_dt%nbasis_add)
1414 output_dt%basis_func_rc(1:input_dt%nbasis)=input_dt%basis_func_rc(1:input_dt%nbasis)
1415
1416 END SUBROUTINE input_dataset_copy
1417
1418
1419!!=================================================================
1420!! NAME
1421!!  input_dataset_read_occ
1422!!
1423!! FUNCTION
1424!!  Initialize only modified occupations in a data-structure by reading them
1425!!  from a file. If file is omitted, then read from standard input.
1426!!  Note: file is supposed to be opened.
1427!!
1428!! INPUTS
1429!!  np(5)= electronic configuration: number of s,p,d,f,g shells
1430!!  diracrel= TRUE if we perform a Dirac relativistic calculation
1431!!  [inputfile_unit]= logical unit of input file to be read (optional)
1432!!  [echofile_unit]= logical unit of a file to echo input file content (optional)
1433!!
1434!! OUTPUT
1435!!  norbit_mod= number of orbital with modified occupations
1436!!  orbit_mod_l(norbit_mod)= l quantum numbers for the modifed orbitals
1437!!  orbit_mod_n(norbit_mod)= n quantum numbers for the modifed orbitals
1438!!  orbit_mod_k(norbit_mod)= kappa quantum numbers for the modifed orbitals
1439!!  orbit_mod_occ(norbit_mod)= occupations for the modifed orbitals
1440!!  [input_dt]= data-structure containing the complete input file (optional).
1441!!              If omitted, then the global public `input_dataset`
1442!!              is used.
1443!!
1444!!=================================================================
1445 SUBROUTINE input_dataset_read_occ(norbit_mod,orbit_mod_l,orbit_mod_n,orbit_mod_k,&
1446 &                          orbit_mod_occ,np,diracrel,inputfile_unit,echofile_unit)
1447
1448!---- Arguments
1449 INTEGER,INTENT(IN) :: np(5)
1450 LOGICAL,INTENT(IN) :: diracrel
1451 INTEGER,INTENT(IN),OPTIONAL :: inputfile_unit,echofile_unit
1452 INTEGER,INTENT(OUT) :: norbit_mod
1453 INTEGER,ALLOCATABLE,INTENT(INOUT) :: orbit_mod_l(:),orbit_mod_n(:),orbit_mod_k(:)
1454 REAL(8),ALLOCATABLE,INTENT(INOUT) :: orbit_mod_occ(:)
1455
1456!---- Local variables
1457 INTEGER,PARAMETER :: nkappa(5)=(/1,2,2,2,2/)
1458 INTEGER :: input_unit,echo_unit,ll,nn,io,ii,kk,ik
1459 LOGICAL :: has_to_echo
1460 CHARACTER(200) :: inputline
1461 REAL(8) :: xocc
1462 INTEGER :: tmp_n(norbit_max),tmp_l(norbit_max),tmp_k(norbit_max)
1463 REAL(8) :: tmp_occ(norbit_max),basis_add_energy(nbasis_add_max)
1464
1465!------------------------------------------------------------------
1466
1467!Select input file logical unit
1468 input_unit=STD_IN;IF (PRESENT(inputfile_unit)) input_unit=inputfile_unit
1469
1470!Do we echo input file content?
1471 has_to_echo=PRESENT(echofile_unit)
1472 echo_unit=-1;IF (has_to_echo) echo_unit=echofile_unit
1473
1474 IF (has_to_ask) THEN
1475   IF (.NOT.diracrel) THEN
1476     WRITE(STD_OUT,*)'Enter np l occ for all occupations for all revisions'
1477     WRITE(STD_OUT,*)'Enter 0 0 0 to end'
1478   ELSE
1479     WRITE(STD_OUT,*)'Enter np l kappa occ for all occupations for all revisions'
1480     WRITE(STD_OUT,*)'Enter 0 0 0 0 to end'
1481   END IF
1482 END IF
1483
1484 norbit_mod=0 ; kk=0
1485 DO
1486   READ(input_unit,'(a)') inputline
1487   IF (has_to_echo) WRITE(echo_unit,'(a)') TRIM(inputline)
1488   if (.not.diracrel) READ(inputline,*) nn,ll,xocc
1489   if (     diracrel) READ(inputline,*) nn,ll,kk,xocc
1490   IF (nn<=0) EXIT
1491   IF (xocc<0.d0.OR.&
1492&    ((.NOT.diracrel).AND.(xocc>2.d0*(2*ll+1))).OR.&
1493&    ((     diracrel).AND.(xocc>2.d0*ABS(kk)))) THEN
1494     WRITE(STD_OUT,*) 'input_dataset: error in occupations -- ip,l,kap,xocc: ',nn,ll,kk,xocc,'!'
1495     STOP
1496   END IF
1497   norbit_mod=norbit_mod+1
1498   if (norbit_mod>norbit_max) stop 'input_dataset_occ: error -- to many occupation lines!'
1499   tmp_l(norbit_mod)=ll
1500   tmp_n(norbit_mod)=nn
1501   tmp_k(norbit_mod)=kk
1502   tmp_occ(norbit_mod)=xocc
1503 END DO
1504
1505 IF (ALLOCATED(orbit_mod_l)) DEALLOCATE(orbit_mod_l)
1506 IF (ALLOCATED(orbit_mod_n)) DEALLOCATE(orbit_mod_n)
1507 IF (ALLOCATED(orbit_mod_n)) DEALLOCATE(orbit_mod_k)
1508 IF (ALLOCATED(orbit_mod_occ)) DEALLOCATE(orbit_mod_occ)
1509 ALLOCATE(orbit_mod_l(norbit_mod))
1510 ALLOCATE(orbit_mod_n(norbit_mod))
1511 ALLOCATE(orbit_mod_k(norbit_mod))
1512 ALLOCATE(orbit_mod_occ(norbit_mod))
1513 orbit_mod_l(1:norbit_mod)=tmp_l(1:norbit_mod)
1514 orbit_mod_n(1:norbit_mod)=tmp_n(1:norbit_mod)
1515 orbit_mod_k(1:norbit_mod)=tmp_k(1:norbit_mod)
1516 orbit_mod_occ(1:norbit_mod)=tmp_occ(1:norbit_mod)
1517
1518!Print read data
1519 IF (has_to_print) THEN
1520   WRITE(STD_OUT,'(3x,a)') "Occupations of orbitals:"
1521   IF (.NOT.diracrel) THEN
1522     WRITE(STD_OUT,'(7x,a)') "n l :  occ"
1523     DO ll=0,4
1524       IF (np(ll+1)>0) THEN
1525         DO ii=1+ll,np(ll+1)
1526           nn=-1
1527           DO io=1,norbit_mod
1528             IF (orbit_mod_l(io)==ll.AND.orbit_mod_n(io)==ii) THEN
1529               nn=io ; EXIT
1530             END IF
1531           END DO
1532           IF (nn<=0) WRITE(STD_OUT,'(7x,i1,1x,i1,a,f4.1)') ii,ll," : ",real(2*(2*ll+1))
1533           IF (nn >0) WRITE(STD_OUT,'(7x,i1,1x,i1,a,f4.1)') ii,ll," : ",orbit_mod_occ(nn)
1534         END DO
1535       END IF
1536     END DO
1537   ELSE
1538     WRITE(STD_OUT,'(7x,a)') "n l kappa :  occ"
1539     DO ll=0,4
1540       IF (np(ll+1)>0) THEN
1541         DO ik=1,nkappa(ll+1)
1542           kk=MERGE(ll,-(ll+1),ik==1);IF (ll==0) kk=-1
1543           DO ii=1+ll,np(ll+1)
1544             nn=-1
1545             DO io=1,norbit_mod
1546               IF (orbit_mod_l(io)==ll.and.orbit_mod_n(io)==ii.AND.orbit_mod_k(io)==kk) THEN
1547                 nn=io ; EXIT
1548               END IF
1549             END DO
1550             IF (nn<=0) WRITE(STD_OUT,'(7x,i1,1x,i1,2x,i2,2x,a,f4.1)') ii,ll,kk," : ",real(2*abs(kk))
1551             IF (nn >0) WRITE(STD_OUT,'(7x,i1,1x,i1,2x,i2,2x,a,f4.1)') ii,ll,kk," : ",orbit_mod_occ(nn)
1552           END DO
1553         END DO
1554       END IF
1555     END DO
1556   END IF
1557 END IF
1558
1559 END SUBROUTINE input_dataset_read_occ
1560
1561
1562!!=================================================================
1563!! NAME
1564!! input_dataset_read_abinit
1565!!
1566!! FUNCTION
1567!! Read the input file (or standard input) in order to get ABINIT options.
1568!! Put them in a input_dataset datastructure.
1569!! If file is omitted, then read from standard input.
1570!! File is supposed to be opened.
1571!!
1572!! INPUTS
1573!!  [inputfile_unit]= logical unit of input file to be read (optional)
1574!!  [echofile_unit]= logical unit of a file to echo input file content (optional)
1575!!
1576!! OUTPUT
1577!!  [input_dt]= data-structure containing the complete input file (optional).
1578!!              If omitted, then the global public `input_dataset`
1579!!    %abinit_usexcnhat=TRUE if nhat density (compensation) is included in XC potential
1580!!    %abinit_prtcorewf= TRUE is printing of core WF is required
1581!!    %abinit_author= author to be printed in the PAW dataset file
1582!!    %abinit_uselog=TRUE if data are transfered on a log. grid before being written
1583!!    %abinit_log_meshsz=mesh size for the logarithmic grid
1584!!    %abinit_log_step=logarithmic step for the logarithmic grid
1585!!    %abinit_userso=TRUE if REAL Space Optimization is required
1586!!    %abinit_rso_ecut=Real Space Optimization parameter: plane wave cutoff = 1/2 Gmax**2
1587!!    %abinit_rso_gfact=Real Space Optimization parameter: Gamma/Gmax ratio
1588!!    %abinit_rso_werror=Real Space Optimization parameter: max. error W_l allowed
1589!!  [abinit_string]= character string containing the ABINIT options line from input file
1590!!
1591!!=================================================================
1592 SUBROUTINE input_dataset_read_abinit(input_dt,inputfile_unit,echofile_unit,abinit_string)
1593
1594!---- Arguments
1595 INTEGER,INTENT(IN),OPTIONAL :: inputfile_unit,echofile_unit
1596 CHARACTER(*),INTENT(OUT),OPTIONAL :: abinit_string
1597 TYPE(input_dataset_t),INTENT(INOUT),OPTIONAL,TARGET :: input_dt
1598
1599!---- Local variables
1600 INTEGER :: input_unit,echo_unit
1601 INTEGER :: i_author,i_usexcnhat,i_prtcorewf,i_logspline,i_rsoptim,iend
1602 LOGICAL :: has_to_echo
1603 CHARACTER(200) :: inputline,inputstring,inputword
1604 TYPE(input_dataset_t),POINTER :: dataset
1605
1606!------------------------------------------------------------------
1607
1608!Select datastruture to be read
1609 IF (PRESENT(input_dt)) THEN
1610   dataset => input_dt
1611 ELSE
1612   dataset => input_dataset
1613 ENDIF
1614
1615!Select input file logical unit
1616 input_unit=STD_IN;IF (PRESENT(inputfile_unit)) input_unit=inputfile_unit
1617
1618!Do we echo input file content?
1619 has_to_echo=PRESENT(echofile_unit)
1620 echo_unit=-1;IF (has_to_echo) echo_unit=echofile_unit
1621
1622!Reading of keywords
1623 READ(input_unit,'(a)') inputline
1624 IF (has_to_echo) WRITE(echo_unit,'(a)') TRIM(inputline)
1625 IF (PRESENT(abinit_string)) abinit_string=TRIM(inputline)
1626 CALL Uppercase(inputline)
1627 i_usexcnhat=INDEX(inputline,'USEXCNHAT')
1628 i_prtcorewf=INDEX(inputline,'PRTCOREWF')
1629 i_logspline=INDEX(inputline,'LOGSPLINE')
1630 i_rsoptim  =INDEX(inputline,'RSOPTIM')
1631 i_author   =INDEX(inputline,'AUTHOR')
1632
1633!Option for core WF printing
1634 dataset%abinit_prtcorewf=MERGE(.true.,PRTCOREWF_DEF,i_prtcorewf>0)
1635
1636!Option for use of NHAT in XC
1637 dataset%abinit_usexcnhat=MERGE(.true.,USEXCNHAT_DEF,i_usexcnhat>0)
1638
1639!To be activated later:
1640!dataset%abinit_vbare=merge(.true.,vbare_def,i_vbare>0)
1641
1642!Options related to the use of REAL SPACE OPTIMIZATION
1643 dataset%abinit_userso=MERGE(.true.,USERSO_DEF,i_rsoptim>0)
1644 dataset%abinit_rso_ecut=RSO_ECUT_DEF
1645 dataset%abinit_rso_gfact=RSO_GFACT_DEF
1646 dataset%abinit_rso_werror=RSO_WERROR_DEF
1647 IF (dataset%abinit_userso) THEN
1648   iend=200
1649   IF (i_usexcnhat>i_rsoptim.AND.i_usexcnhat-1<iend) iend=i_usexcnhat-1
1650   IF (i_prtcorewf>i_rsoptim.AND.i_prtcorewf-1<iend) iend=i_prtcorewf-1
1651   IF (i_logspline>i_rsoptim.AND.i_logspline-1<iend) iend=i_logspline-1
1652   IF (i_author   >i_rsoptim.AND.i_author   -1<iend) iend=i_author   -1
1653   inputstring="";IF (iend>i_rsoptim+7) inputstring=TRIM(inputline(i_rsoptim+7:iend))
1654   IF (inputstring/="") THEN
1655     CALL extractword(1,inputstring,inputword);inputword=TRIM(inputword)
1656     IF (inputword/="") THEN
1657       READ(inputword,*) dataset%abinit_rso_ecut
1658       CALL extractword(2,inputstring,inputword);inputword=TRIM(inputword)
1659       IF (inputword/="") THEN
1660         READ(inputword,*) dataset%abinit_rso_gfact
1661         CALL extractword(3,inputstring,inputword);inputword=TRIM(inputword)
1662         IF (inputword/="") READ(inputword,*) dataset%abinit_rso_werror
1663       END IF
1664     END IF
1665   END IF
1666 END IF
1667
1668!Options for related to the transfer to a reduced logarihmic grid
1669 dataset%abinit_uselog=MERGE(.true.,USELOG_DEF,i_logspline>0)
1670 dataset%abinit_log_meshsz=LOGGRD_SIZE_DEF
1671 dataset%abinit_log_step=LOGGRD_STEP_DEF
1672 IF (dataset%abinit_uselog) THEN
1673   iend=200
1674   IF (i_usexcnhat>i_logspline.AND.i_usexcnhat-1<iend) iend=i_usexcnhat-1
1675   IF (i_prtcorewf>i_logspline.AND.i_prtcorewf-1<iend) iend=i_prtcorewf-1
1676   IF (i_rsoptim  >i_logspline.AND.i_rsoptim  -1<iend) iend=i_rsoptim  -1
1677   IF (i_author   >i_logspline.AND.i_author   -1<iend) iend=i_author   -1
1678   inputstring="";IF (iend>i_logspline+9) inputstring=TRIM(inputline(i_logspline+9:iend))
1679   IF (inputstring/="") THEN
1680     CALL extractword(1,inputstring,inputword);inputword=TRIM(inputword)
1681     IF (inputword/="") THEN
1682       READ(inputword,*) dataset%abinit_log_meshsz
1683       CALL extractword(2,inputstring,inputword);inputword=TRIM(inputword)
1684       IF (inputword/="") READ(inputword,*) dataset%abinit_log_step
1685     END IF
1686   END IF
1687 END IF
1688
1689!Author to be mentioned in the ABINIT file
1690 dataset%abinit_author=TRIM(AUTHOR_DEF)
1691 IF (i_author>0) then
1692   iend=200
1693   IF (i_usexcnhat>i_author.AND.i_usexcnhat-1<iend) iend=i_usexcnhat-1
1694   IF (i_prtcorewf>i_author.AND.i_prtcorewf-1<iend) iend=i_prtcorewf-1
1695   IF (i_logspline>i_author.AND.i_logspline-1<iend) iend=i_logspline-1
1696   IF (i_rsoptim  >i_author.AND.i_rsoptim  -1<iend) iend=i_rsoptim  -1
1697   inputstring=TRIM(inputline(i_author+6:iend))
1698   READ(unit=inputstring,fmt=*) dataset%abinit_author
1699   dataset%abinit_author=TRIM(dataset%abinit_author)
1700 END IF
1701
1702!Print read data
1703 IF (has_to_print) THEN
1704   WRITE(STD_OUT,'(/,2x,a)')'Options for ABINIT file:'
1705   WRITE(STD_OUT,'(2x,2a)') 'ABINIT option : output of core wave funcs =',MERGE("YES"," NO",dataset%abinit_prtcorewf)
1706   WRITE(STD_OUT,'(2x,2a)') 'ABINIT option : transfer to a log grid    =',MERGE("YES"," NO",dataset%abinit_uselog)
1707   IF (dataset%abinit_uselog) THEN
1708     WRITE(STD_OUT,'(5x,a,i5)') '- mesh size of the log grid =',dataset%abinit_log_meshsz
1709     WRITE(STD_OUT,'(5x,a,g8.2)') '- step of the log grid      = ',dataset%abinit_log_step
1710   END IF
1711   WRITE(STD_OUT,'(2x,2a)') 'ABINIT option : Real Space Optimization   =',MERGE("YES"," NO",dataset%abinit_userso)
1712   IF (dataset%abinit_userso) THEN
1713     WRITE(STD_OUT,'(5x,a,f7.2)') '- RSO plane wave cutoff =',dataset%abinit_rso_ecut
1714     WRITE(STD_OUT,'(5x,a,1x,f6.2)') '- RSO Gamma/Gmax ratio  =',dataset%abinit_rso_gfact
1715     WRITE(STD_OUT,'(5x,a,g11.3)') '- RSO maximum error     = ',dataset%abinit_rso_werror
1716   END IF
1717   WRITE(STD_OUT,'(2x,2a)') 'ABINIT option : use of compensation in XC =',MERGE("YES"," NO",dataset%abinit_usexcnhat)
1718   IF (dataset%abinit_usexcnhat) THEN
1719     WRITE(STD_OUT,'(5x,a)') '- Kresse local ionic potential output in XML file'
1720   ELSE
1721     WRITE(STD_OUT,'(5x,a)') '- Blochl local ionic potential output in XML file'
1722   END IF
1723 END IF
1724
1725 END SUBROUTINE input_dataset_read_abinit
1726
1727
1728!!=================================================================
1729!! NAME
1730!! input_dataset_read_xml
1731!!
1732!! FUNCTION
1733!! Read the input file (or standard input) in order to get UPF options.
1734!! Put them in a input_dataset datastructure.
1735!! If file is omitted, then read from standard input.
1736!! File is supposed to be opened.
1737!!
1738!! INPUTS
1739!!  [inputfile_unit]= logical unit of input file to be read (optional)
1740!!  [echofile_unit]= logical unit of a file to echo input file content (optional)
1741!!
1742!! OUTPUT
1743!!  [input_dt]= data-structure containing the complete input file (optional).
1744!!              If omitted, then the global public `input_dataset`
1745!!    %_usexcnhat=TRUE if nhat density (compensation) is included in XC potential
1746!!    %xml_prtcorewf= TRUE is printing of core WF is required
1747!!    %xml_author= author to be printed in the PAW dataset file
1748!!    %xml_comment= comment line to be printed in the heder of the PAW dataset file
1749!!    %xml_uselog=TRUE if data are transfered on a log. grid before being written
1750!!    %xml_log_meshsz=mesh size for the logarithmic grid
1751!!    %xml_log_step=logarithmic step for the logarithmic grid
1752!!    %xml_userso=TRUE if REAL Space Optimization is required
1753!!    %xml_rso_ecut=Real Space Optimization parameter: plane wave cutoff = 1/2 Gmax**2
1754!!    %xml_rso_gfact=Real Space Optimization parameter: Gamma/Gmax ratio
1755!!    %xml_rso_werror=Real Space Optimization parameter: max. error W_l allowed
1756!!    %xml_uselda12= TRUE if LDA-1/2 potential calculation is required
1757!!    %xml_lda12_orb_l=LDA-1/2 parameter: l quantum number of the orbital to be ionized
1758!!    %xml_lda12_orb_n=LDA-1/2 parameter: n quantum number of the orbital to be ionized
1759!!    %xml_lda12_ion=LDA-1/2 parameter: amount of charge to be removed from the ionized orbital
1760!!    %xml_lda12_rcut=LDA-1/2 parameter: cut-off radius (in bohr)
1761!!    %xml_lda12_logfile=LDA-1/2 parameter: name of the log file
1762!!  [xml_string]= character string containing the ABINIT options line from input file
1763!!
1764!!=================================================================
1765 SUBROUTINE input_dataset_read_xml(input_dt,inputfile_unit,echofile_unit,xml_string)
1766
1767!---- Arguments
1768 INTEGER,INTENT(IN),OPTIONAL :: inputfile_unit,echofile_unit
1769 CHARACTER(*),INTENT(OUT),OPTIONAL :: xml_string
1770 TYPE(input_dataset_t),INTENT(INOUT),OPTIONAL,TARGET :: input_dt
1771
1772!---- Local variables
1773 INTEGER :: input_unit,echo_unit
1774 INTEGER :: i_author,i_comment,i_usexcnhat,i_prtcorewf,i_logspline,i_rsoptim,i_lda12,iend
1775 LOGICAL :: has_to_echo
1776 CHARACTER(200) :: inputline,inputstring,inputword
1777 TYPE(input_dataset_t),POINTER :: dataset
1778
1779!------------------------------------------------------------------
1780
1781!Select datastruture to be read
1782 IF (PRESENT(input_dt)) THEN
1783   dataset => input_dt
1784 ELSE
1785   dataset => input_dataset
1786 ENDIF
1787
1788!Select input file logical unit
1789 input_unit=STD_IN;IF (PRESENT(inputfile_unit)) input_unit=inputfile_unit
1790
1791!Do we echo input file content?
1792 has_to_echo=PRESENT(echofile_unit)
1793 echo_unit=-1;IF (has_to_echo) echo_unit=echofile_unit
1794
1795!Reading of keywords
1796 READ(input_unit,'(a)') inputline
1797 IF (has_to_echo) WRITE(echo_unit,'(a)') TRIM(inputline)
1798 IF (PRESENT(xml_string)) xml_string=TRIM(inputline)
1799 CALL Uppercase(inputline)
1800 i_usexcnhat=INDEX(inputline,'USEXCNHAT')
1801 i_prtcorewf=INDEX(inputline,'PRTCOREWF')
1802 i_logspline=INDEX(inputline,'WITHSPLGRID')
1803 i_rsoptim  =INDEX(inputline,'RSOPTIM')
1804 i_lda12    =INDEX(inputline,'LDA12')
1805 i_author   =INDEX(inputline,'AUTHOR')
1806 i_comment  =INDEX(inputline,'COMMENT')
1807
1808!Option for core WF printing
1809 dataset%xml_prtcorewf=MERGE(.true.,PRTCOREWF_DEF,i_prtcorewf>0)
1810
1811!Option for use of NHAT in XC
1812 dataset%xml_usexcnhat=MERGE(.true.,USEXCNHAT_DEF,i_usexcnhat>0)
1813
1814!To be activated later:
1815!dataset%xml_vbare=merge(.true.,vbare_def,i_vbare>0)
1816
1817!Options related to the use of REAL SPACE OPTIMIZATION
1818 dataset%xml_userso=MERGE(.true.,USERSO_DEF,i_rsoptim>0)
1819 dataset%xml_rso_ecut=RSO_ECUT_DEF
1820 dataset%xml_rso_gfact=RSO_GFACT_DEF
1821 dataset%xml_rso_werror=RSO_WERROR_DEF
1822 IF (dataset%xml_userso) THEN
1823   iend=200
1824   IF (i_usexcnhat>i_rsoptim.AND.i_usexcnhat-1<iend) iend=i_usexcnhat-1
1825   IF (i_prtcorewf>i_rsoptim.AND.i_prtcorewf-1<iend) iend=i_prtcorewf-1
1826   IF (i_logspline>i_rsoptim.AND.i_logspline-1<iend) iend=i_logspline-1
1827   IF (i_lda12    >i_rsoptim.AND.i_lda12    -1<iend) iend=i_lda12    -1
1828   IF (i_author   >i_rsoptim.AND.i_author   -1<iend) iend=i_author   -1
1829   IF (i_comment  >i_rsoptim.AND.i_comment  -1<iend) iend=i_comment  -1
1830   inputstring="";IF (iend>i_rsoptim+7) inputstring=TRIM(inputline(i_rsoptim+7:iend))
1831   IF (inputstring/="") THEN
1832     CALL extractword(1,inputstring,inputword);inputword=TRIM(inputword)
1833     IF (inputword/="") THEN
1834       READ(inputword,*) dataset%xml_rso_ecut
1835       CALL extractword(2,inputstring,inputword);inputword=TRIM(inputword)
1836       IF (inputword/="") THEN
1837         READ(inputword,*) dataset%xml_rso_gfact
1838         CALL extractword(3,inputstring,inputword);inputword=TRIM(inputword)
1839         IF (inputword/="") READ(inputword,*) dataset%xml_rso_werror
1840       END IF
1841     END IF
1842   END IF
1843 END IF
1844
1845!Options related to the use of LDA-1/2 technique
1846 dataset%xml_uselda12=MERGE(.true.,USELDA12_DEF,i_lda12>0)
1847 dataset%xml_lda12_orb_l=LDA12_ORB_L_DEF
1848 dataset%xml_lda12_orb_n=LDA12_ORB_N_DEF
1849 dataset%xml_lda12_ion=LDA12_ION_DEF
1850 dataset%xml_lda12_rcut=LDA12_RCUT_DEF
1851 dataset%xml_lda12_logfile=TRIM(LDA12_LOGFILE)
1852 IF (dataset%xml_uselda12) THEN
1853   iend=200
1854   IF (i_usexcnhat>i_lda12.AND.i_usexcnhat-1<iend) iend=i_usexcnhat-1
1855   IF (i_prtcorewf>i_lda12.AND.i_prtcorewf-1<iend) iend=i_prtcorewf-1
1856   IF (i_logspline>i_lda12.AND.i_logspline-1<iend) iend=i_logspline-1
1857   IF (i_rsoptim  >i_lda12.AND.i_rsoptim  -1<iend) iend=i_rsoptim  -1
1858   IF (i_author   >i_lda12.AND.i_author   -1<iend) iend=i_author   -1
1859   IF (i_comment  >i_lda12.AND.i_comment  -1<iend) iend=i_comment  -1
1860   inputstring="";IF (iend>i_lda12+5) inputstring=TRIM(inputline(i_lda12+5:iend))
1861   IF (inputstring/="") THEN
1862     CALL extractword(1,inputstring,inputword);inputword=TRIM(inputword)
1863     IF (inputword/="") THEN
1864       READ(inputword,*) dataset%xml_lda12_orb_n
1865       CALL extractword(2,inputstring,inputword);inputword=TRIM(inputword)
1866       IF (inputword/="") THEN
1867         READ(inputword,*) dataset%xml_lda12_orb_l
1868         CALL extractword(3,inputstring,inputword);inputword=TRIM(inputword)
1869         IF (inputword/="") THEN
1870           READ(inputword,*) dataset%xml_lda12_ion
1871           CALL extractword(4,inputstring,inputword);inputword=TRIM(inputword)
1872           IF (inputword/="") READ(inputword,*) dataset%xml_lda12_rcut
1873         END IF
1874       END IF
1875     END IF
1876   END IF
1877 END IF
1878
1879!Options for related to the transfer to a reduced logarihmic grid
1880 dataset%xml_usespl=MERGE(.true.,USELOG_DEF,i_logspline>0)
1881 dataset%xml_spl_meshsz=LOGGRD_SIZE_DEF
1882 IF (dataset%xml_usespl) THEN
1883   iend=200
1884   IF (i_usexcnhat>i_logspline.AND.i_usexcnhat-1<iend) iend=i_usexcnhat-1
1885   IF (i_prtcorewf>i_logspline.AND.i_prtcorewf-1<iend) iend=i_prtcorewf-1
1886   IF (i_rsoptim  >i_logspline.AND.i_rsoptim  -1<iend) iend=i_rsoptim  -1
1887   IF (i_lda12    >i_logspline.AND.i_lda12    -1<iend) iend=i_lda12    -1
1888   IF (i_author   >i_logspline.AND.i_author   -1<iend) iend=i_author   -1
1889   IF (i_comment  >i_logspline.AND.i_comment  -1<iend) iend=i_comment  -1
1890   inputstring="";IF (iend>i_logspline+11) inputstring=TRIM(inputline(i_logspline+11:iend))
1891   IF (inputstring/="") THEN
1892     CALL extractword(1,inputstring,inputword);inputword=TRIM(inputword)
1893     IF (inputword/="") READ(inputword,*) dataset%xml_spl_meshsz
1894   END IF
1895 END IF
1896
1897!Author to be mentioned in the XML file
1898 dataset%xml_author=TRIM(AUTHOR_DEF)
1899 IF (i_author>0) then
1900   iend=200
1901   IF (i_usexcnhat>i_author.AND.i_usexcnhat-1<iend) iend=i_usexcnhat-1
1902   IF (i_prtcorewf>i_author.AND.i_prtcorewf-1<iend) iend=i_prtcorewf-1
1903   IF (i_logspline>i_author.AND.i_logspline-1<iend) iend=i_logspline-1
1904   IF (i_rsoptim  >i_author.AND.i_rsoptim  -1<iend) iend=i_rsoptim  -1
1905   IF (i_lda12    >i_author.AND.i_lda12    -1<iend) iend=i_lda12    -1
1906   IF (i_comment  >i_author.AND.i_comment  -1<iend) iend=i_comment  -1
1907   inputstring=TRIM(inputline(i_author+6:))
1908   READ(unit=inputstring,fmt=*) dataset%xml_author
1909   dataset%xml_author=TRIM(dataset%xml_author)
1910 END IF
1911
1912!Comment to be added in the XML file
1913 dataset%xml_comment=TRIM(COMMENT_DEF)
1914 IF (i_comment>0) then
1915   iend=200
1916   IF (i_usexcnhat>i_comment.AND.i_usexcnhat-1<iend) iend=i_usexcnhat-1
1917   IF (i_prtcorewf>i_comment.AND.i_prtcorewf-1<iend) iend=i_prtcorewf-1
1918   IF (i_logspline>i_comment.AND.i_logspline-1<iend) iend=i_logspline-1
1919   IF (i_rsoptim  >i_comment.AND.i_rsoptim  -1<iend) iend=i_rsoptim  -1
1920   IF (i_lda12    >i_comment.AND.i_lda12    -1<iend) iend=i_lda12    -1
1921   IF (i_author   >i_comment.AND.i_author   -1<iend) iend=i_author   -1
1922   inputstring=TRIM(inputline(i_comment+7:))
1923   READ(unit=inputstring,fmt=*) dataset%xml_comment
1924   dataset%xml_comment=TRIM(dataset%xml_comment)
1925 END IF
1926
1927!Print read data
1928 IF (has_to_print) THEN
1929   WRITE(STD_OUT,'(/,2x,a)')'Options for XML file:'
1930   WRITE(STD_OUT,'(2x,2a)') 'XML option : output of core wave funcs =',MERGE("YES"," NO",dataset%xml_prtcorewf)
1931   WRITE(STD_OUT,'(2x,2a)') 'XML option : spline to a log grid      =',MERGE("YES"," NO",dataset%xml_usespl)
1932   IF (dataset%xml_usespl) THEN
1933     WRITE(STD_OUT,'(5x,a,i5)') '- mesh size of the log grid =',dataset%xml_spl_meshsz
1934   END IF
1935   WRITE(STD_OUT,'(2x,2a)') 'XML option : Real Space Optimization   =',MERGE("YES"," NO",dataset%xml_userso)
1936   IF (dataset%xml_userso) THEN
1937     WRITE(STD_OUT,'(5x,a,f7.2)') '- RSO plane wave cutoff =',dataset%xml_rso_ecut
1938     WRITE(STD_OUT,'(5x,a,1x,f6.2)') '- RSO Gamma/Gmax ratio  =',dataset%xml_rso_gfact
1939     WRITE(STD_OUT,'(5x,a,g11.3)') '- RSO maximum error     = ',dataset%xml_rso_werror
1940   END IF
1941   WRITE(STD_OUT,'(2x,2a)') 'XML option : output of LDA-1/2 pot.    =',MERGE("YES"," NO",dataset%xml_uselda12)
1942   IF (dataset%xml_uselda12) THEN
1943     WRITE(STD_OUT,'(5x,a,i2)') '- LDA-1/2 ionized orbital n =',dataset%xml_lda12_orb_n
1944     WRITE(STD_OUT,'(5x,a,i2)') '- LDA-1/2 ionized orbital l =',dataset%xml_lda12_orb_l
1945     WRITE(STD_OUT,'(5x,a,f5.2)') '- LDA-1/2 ionization        =',dataset%xml_lda12_ion
1946     WRITE(STD_OUT,'(5x,a,f7.4)') '- LDA-1/2 cutoff radius (au)=',dataset%xml_lda12_rcut
1947     WRITE(STD_OUT,'(5x,3a)') '- LDA-1/2: see ''',trim(dataset%xml_lda12_logfile),&
1948&                             ''' file to check convergence'
1949   END IF
1950   WRITE(STD_OUT,'(2x,2a)') 'XML option : use of compensation in XC =',MERGE("YES"," NO",dataset%xml_usexcnhat)
1951   IF (dataset%xml_usexcnhat) THEN
1952     WRITE(STD_OUT,'(5x,a)') '- Zero potential and Kresse local ionic potential output in XML file'
1953   ELSE
1954     WRITE(STD_OUT,'(5x,a)') '- Zero potential and Blochl local ionic potential output in XML file'
1955   END IF
1956 END IF
1957
1958 END SUBROUTINE input_dataset_read_xml
1959
1960
1961!!=================================================================
1962!! NAME
1963!! input_dataset_read_upf
1964!!
1965!! FUNCTION
1966!! Read the input file (or standard input) in order to get XML options.
1967!! Put them in a input_dataset datastructure.
1968!! If file is omitted, then read from standard input.
1969!! File is supposed to be opened.
1970!!
1971!! INPUTS
1972!!  [inputfile_unit]= logical unit of input file to be read (optional)
1973!!  [echofile_unit]= logical unit of a file to echo input file content (optional)
1974!!
1975!! OUTPUT
1976!!  [input_dt]= data-structure containing the complete input file (optional).
1977!!              If omitted, then the global public `input_dataset`
1978!!    %upf_grid_xmin= minimum radius given by the grid
1979!!    %upf_grid_zmesh= inverse of the radial step of the grid
1980!!    %upf_grid_dx= logarithmic step of the grid
1981!!    %upf_grid_range= range of the grid
1982!!  [upf_string]= character string containing the ABINIT options line from input file
1983!!
1984!!=================================================================
1985 SUBROUTINE input_dataset_read_upf(input_dt,inputfile_unit,echofile_unit,upf_string)
1986
1987!---- Arguments
1988 INTEGER,INTENT(IN),OPTIONAL :: inputfile_unit,echofile_unit
1989 CHARACTER(*),INTENT(OUT),OPTIONAL :: upf_string
1990 TYPE(input_dataset_t),INTENT(INOUT),OPTIONAL,TARGET :: input_dt
1991
1992!---- Local variables
1993 INTEGER :: input_unit,echo_unit
1994 INTEGER :: i_all,i_dx,i_xmin,i_zmesh,i_range
1995 LOGICAL :: has_to_echo
1996 CHARACTER(200) :: inputline
1997 TYPE(input_dataset_t),POINTER :: dataset
1998
1999!------------------------------------------------------------------
2000
2001!Select datastruture to be read
2002 IF (PRESENT(input_dt)) THEN
2003   dataset => input_dt
2004 ELSE
2005   dataset => input_dataset
2006 ENDIF
2007
2008!Select input file logical unit
2009 input_unit=STD_IN;IF (PRESENT(inputfile_unit)) input_unit=inputfile_unit
2010
2011!Do we echo input file content?
2012 has_to_echo=PRESENT(echofile_unit)
2013 echo_unit=-1;IF (has_to_echo) echo_unit=echofile_unit
2014
2015!Reading of keywords
2016 READ(input_unit,'(a)') inputline
2017 IF (has_to_echo) WRITE(echo_unit,'(a)') TRIM(inputline)
2018 IF (PRESENT(upf_string)) upf_string=TRIM(inputline)
2019 CALL Uppercase(inputline)
2020
2021!Default values
2022 dataset%upf_grid_dx=UPF_DX_DEF
2023 dataset%upf_grid_xmin=UPF_XMIN_DEF
2024 dataset%upf_grid_zmesh=UPF_ZMESH_DEF
2025 dataset%upf_grid_range=UPF_RANGE_DEF
2026
2027 i_all=INDEX(inputline,'UPFDX')+INDEX(inputline,'UPFXMIN')+&
2028&      INDEX(inputline,'UPFZMESH')+INDEX(inputline,'UPFRANGE')
2029 IF (i_all>0) then
2030
2031!  Logarithmic step
2032   i_dx=INDEX(inputline,'UPFDX')
2033   IF (i_dx>0) READ(inputline(i_dx+5:256),*) dataset%upf_grid_dx
2034
2035!  Minimum radius
2036   i_xmin=INDEX(inputline,'UPFXMIN')
2037   IF (i_xmin>0) READ(inputline(i_xmin+7:256),*) dataset%upf_grid_xmin
2038
2039!  Inverse of radial step
2040   i_zmesh=INDEX(inputline,'UPFZMESH')
2041   IF (i_zmesh>0) READ(inputline(i_zmesh+8:256),*) dataset%upf_grid_zmesh
2042
2043!  Range
2044   i_range=INDEX(inputline,'UPFRANGE')
2045   IF (i_range>0) READ(inputline(i_range+8:256),*) dataset%upf_grid_range
2046
2047 END IF
2048
2049!Print read data
2050 IF (has_to_print) THEN
2051   WRITE(STD_OUT,'(/,2x,a)')      'Options for UPF file:'
2052   WRITE(STD_OUT,'(2x,a,es10.3)') 'UPF grid option : logarithmic step (upfdx)      =',dataset%upf_grid_dx
2053   WRITE(STD_OUT,'(2x,a,es10.3)') 'UPF grid option : minimum grid x (upfxmin)      =',dataset%upf_grid_xmin
2054   WRITE(STD_OUT,'(2x,a,es10.3)') 'UPF grid option : inv. of radial step (upfzmesh)=',dataset%upf_grid_zmesh
2055   WRITE(STD_OUT,'(2x,a,es10.3)') 'UPF grid option : grid range (upf_gridrange)    =',dataset%upf_grid_range
2056 END IF
2057
2058 END SUBROUTINE input_dataset_read_upf
2059
2060!!=================================================================
2061
2062END MODULE input_dataset_mod
2063