1!
2! Copyright (C) 2020 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!--------------------------------------------------------
9MODULE  read_upf_new_module
10  !-----------------------------------------------------
11  !! this module contains the simplified code for reading
12  !! pseudopotential files in either UPF v.2 or xml
13  !
14  USE xmltools
15  USE upf_kinds, ONLY: dp
16  USE pseudo_types, ONLY: pseudo_upf, pseudo_config
17  !
18  LOGICAL :: v2
19  !! true if UPF v.2 version, false if new UPF with xml schema
20  INTEGER :: iun
21  !! unit for reading data
22  !
23  PUBLIC
24  !
25CONTAINS
26  !
27  !------------------------------------------------+
28  SUBROUTINE read_upf_new (filename, upf, ierr)         !
29    !---------------------------------------------+
30    !! Reads pseudopotential in UPF format (either v.2 or upf_schema).
31    !! Derived-type variable *upf* store in output the data read from file.
32    !! File *filename* is opened and closed inside the routine
33    !
34    IMPLICIT NONE
35    CHARACTER(len=*), INTENT(IN) :: filename
36    !! i/o filename
37    TYPE(pseudo_upf),INTENT(OUT) :: upf
38    !! the derived type storing the pseudo data
39    INTEGER, INTENT(OUT) :: ierr
40    !! ierr=0  : xml schema, ierr=-2: UPF v.2
41    !! ierr=-81: error reading PP file
42    !
43    iun = xml_openfile ( filename )
44    IF ( iun == -1 ) CALL upf_error('read_upf', 'cannot open file',1)
45    call xmlr_opentag ( 'qe_pp:pseudo', IERR = ierr )
46    if ( ierr == 0 ) then
47       v2 =.false.
48    else if ( ierr == -1 ) then
49       rewind (iun)
50       call xmlr_opentag ( 'UPF', IERR = ierr )
51       if ( ierr == 0 ) then
52          v2 =.true.
53          ierr = -2
54          CALL get_attr ( 'version', upf%nv )
55       end if
56    end if
57    if ( ierr /= 0 .and. ierr /= -2 ) then
58       call xml_closefile( )
59       ierr = -81
60       return
61    end if
62    !
63    ! The header sections differ a lot between UPF v.2 and UPF with schema
64    !
65    IF ( v2 ) THEN
66       CALL read_pp_header_v2 ( upf )
67    ELSE
68       CALL read_pp_header_schema ( upf )
69    END IF
70    ! compatibility
71    upf%is_gth = .false.
72    upf%is_multiproj = .true.
73    !
74    ! From here on the format of v2 and schema do not differ much:
75    ! the most frequent difference is capitalization of tags
76    ! (see function capitalize_if_v2)
77    !
78    CALL read_pp_mesh ( upf )
79    !
80    allocate ( upf%rho_atc(upf%mesh) )
81    IF(upf%nlcc) then
82       CALL xmlr_readtag( capitalize_if_v2('pp_nlcc'), &
83            upf%rho_atc(:) )
84    else
85       upf%rho_atc(:) = 0.0_dp
86    end if
87    IF( .NOT. upf%tcoulombp) then
88       allocate ( upf%vloc(upf%mesh) )
89       CALL xmlr_readtag( capitalize_if_v2('pp_local'), &
90            upf%vloc(:), ierr )
91       !
92       ! existing PP files may have pp_nlcc first, pp_local later,
93       ! but also the other way round - check that everything was right
94       !
95       if ( ierr /= 0 ) then
96          ierr = -81
97          return
98       end if
99    end if
100    !
101    CALL read_pp_semilocal ( upf )
102    !
103    CALL read_pp_nonlocal ( upf )
104    !
105    CALL read_pp_pswfc ( upf )
106    !
107    CALL read_pp_full_wfc ( upf )
108    !
109    allocate( upf%rho_at(1:upf%mesh) )
110    CALL xmlr_readtag( capitalize_if_v2('pp_rhoatom'), &
111         upf%rho_at(1:upf%mesh) )
112    !
113    CALL read_pp_spinorb ( upf )
114    !
115    CALL read_pp_paw ( upf )
116    !
117    CALL read_pp_gipaw ( upf )
118    !
119    ! close initial tag, qe_pp:pseudo or UPF
120    !
121    CALL xmlr_closetag ( )
122    !
123    CALL xml_closefile ( )
124    !
125  END SUBROUTINE read_upf_new
126  !
127  FUNCTION capitalize_if_v2 ( strin ) RESULT ( strout )
128    !
129    ! returns a capitalized string for UPF v.2, the same string otherwise
130    ! (UPF v.2 uses capitalized tags, UPF with schema use lowercase)
131    !
132    USE upf_utils, ONLY: capital
133    IMPLICIT NONE
134    CHARACTER(LEN=*) :: strin
135    !
136    INTEGER :: n
137    CHARACTER(LEN=:), ALLOCATABLE :: strout
138    !
139    IF ( v2 ) THEN
140       strout = ''
141       DO n = 1,LEN_TRIM(strin)
142          strout = strout // capital(strin(n:n))
143       END DO
144    ELSE
145       strout = TRIM(strin)
146    END IF
147    !
148  END FUNCTION capitalize_if_v2
149  !--------------------------------------------------------
150  SUBROUTINE read_pp_header_schema ( upf )
151    !--------------------------------------------------------
152    !
153    IMPLICIT NONE
154    TYPE(pseudo_upf), INTENT(INOUT) :: upf ! the pseudo data
155    !
156    CALL xmlr_opentag( capitalize_if_v2('pp_header') )
157    !
158    CALL xmlr_readtag( 'element', upf%psd )
159    CALL xmlr_readtag( 'z_valence', upf%zp )
160    CALL xmlr_readtag( 'type', upf%typ )
161    CALL xmlr_readtag( 'functional', upf%dft )
162    CALL xmlr_readtag( 'relativistic', upf%rel )
163    CALL xmlr_readtag( 'is_ultrasoft', upf%tvanp )
164    CALL xmlr_readtag( 'is_paw', upf%tpawp )
165    CALL xmlr_readtag( 'is_coulomb', upf%tcoulombp )
166    CALL xmlr_readtag( 'has_so', upf%has_so )
167    CALL xmlr_readtag( 'has_wfc', upf%has_wfc )
168    CALL xmlr_readtag( 'has_gipaw', upf%has_gipaw )
169    CALL xmlr_readtag( 'paw_as_gipaw', upf%paw_as_gipaw)
170    CALL xmlr_readtag( 'core_correction', upf%nlcc)
171    CALL xmlr_readtag( 'total_psenergy', upf%etotps )
172    CALL xmlr_readtag( 'wfc_cutoff', upf%ecutwfc )
173    CALL xmlr_readtag( 'rho_cutoff', upf%ecutrho )
174    CALL xmlr_readtag( 'l_max', upf%lmax )
175    CALL xmlr_readtag( 'l_max_rho', upf%lmax_rho )
176    CALL xmlr_readtag( 'l_local', upf%lloc )
177    CALL xmlr_readtag( 'mesh_size', upf%mesh )
178    CALL xmlr_readtag( 'number_of_wfc', upf%nwfc )
179    CALL xmlr_readtag( 'number_of_proj', upf%nbeta )
180    !
181    CALL xmlr_closetag( )
182    !
183  END SUBROUTINE read_pp_header_schema
184  !
185  !--------------------------------------------------------
186  SUBROUTINE read_pp_header_v2 ( upf )
187    !--------------------------------------------------------
188    !
189    IMPLICIT NONE
190    TYPE(pseudo_upf), INTENT(INOUT) :: upf ! the pseudo data
191    !
192    CHARACTER(LEN=1) :: dummy
193    !
194    CALL xmlr_readtag ( capitalize_if_v2('pp_header'), dummy )
195    CALL get_attr ('generated', upf%generated)
196    CALL get_attr ('author', upf%author)
197    CALL get_attr ('date', upf%date)
198    CALL get_attr ('comment', upf%comment)
199    CALL get_attr ('element', upf%psd)
200    CALL get_attr ('pseudo_type', upf%typ)
201    CALL get_attr ('relativistic', upf%rel)
202    CALL get_attr ('is_ultrasoft', upf%tvanp)
203    CALL get_attr ('is_paw', upf%tpawp)
204    CALL get_attr ('is_coulomb', upf%tcoulombp)
205    CALL get_attr ('has_so', upf%has_so)
206    CALL get_attr ('has_wfc', upf%has_wfc)
207    CALL get_attr ('has_gipaw', upf%has_gipaw)
208    CALL get_attr ('paw_as_gipaw', upf%paw_as_gipaw)
209    CALL get_attr ('core_correction', upf%nlcc)
210    CALL get_attr ('functional', upf%dft)
211    CALL get_attr ('z_valence', upf%zp)
212    CALL get_attr ('total_psenergy', upf%etotps)
213    CALL get_attr ('wfc_cutoff', upf%ecutwfc)
214    CALL get_attr ('rho_cutoff', upf%ecutrho)
215    CALL get_attr ('l_max', upf%lmax)
216    CALL get_attr ('l_max_rho', upf%lmax_rho)
217    CALL get_attr ('l_local', upf%lloc)
218    CALL get_attr ('mesh_size', upf%mesh)
219    CALL get_attr ('number_of_wfc', upf%nwfc)
220    CALL get_attr ('number_of_proj', upf%nbeta )
221    !
222  END SUBROUTINE read_pp_header_v2
223  !
224  !--------------------------------------------------------
225  SUBROUTINE read_pp_mesh ( upf )
226    !--------------------------------------------------------
227    !
228    IMPLICIT NONE
229    TYPE(pseudo_upf),INTENT(INOUT) :: upf ! the pseudo data
230    integer :: mesh
231    !
232    CALL xmlr_opentag( capitalize_if_v2('pp_mesh') )
233    CALL get_attr ( 'mesh', mesh )
234    if ( mesh == 0 ) THEN
235       call upf_error('read_pp_mesh',&
236         'mesh size missing, using the one in header',-1)
237    else if ( mesh /= upf%mesh ) THEN
238       call upf_error('read_pp_mesh',&
239         'mismatch in mesh size, discarding the one in header',-1)
240       upf%mesh = mesh
241    end if
242    CALL get_attr ( 'dx'  , upf%dx   )
243    CALL get_attr ( 'xmin', upf%xmin )
244    CALL get_attr ( 'rmax', upf%rmax )
245    CALL get_attr ( 'zmesh', upf%zmesh )
246    allocate ( upf%r(1:upf%mesh) )
247    CALL xmlr_readtag( capitalize_if_v2('pp_r'), upf%r(1:upf%mesh) )
248    allocate ( upf%rab(1:upf%mesh) )
249    CALL xmlr_readtag( capitalize_if_v2('pp_rab'), upf%rab(1:upf%mesh) )
250    !
251    CALL xmlr_closetag( ) ! end pp_mesh
252    !
253  END SUBROUTINE read_pp_mesh
254  !
255  !--------------------------------------------------------
256  SUBROUTINE read_pp_semilocal ( upf )
257    !--------------------------------------------------------
258    !
259    IMPLICIT NONE
260    TYPE(pseudo_upf),INTENT(INOUT) :: upf ! the pseudo data
261    !
262    INTEGER :: nb, ind, l, j, ierr
263    CHARACTER(LEN=8) :: tag
264    real(dp), allocatable :: vnl(:)
265    !
266    IF ( upf%typ == "SL" ) THEN
267       !
268       IF ( upf%has_so ) then
269          ALLOCATE(upf%vnl(upf%mesh,0:upf%lmax,2))
270       else
271          ALLOCATE(upf%vnl(upf%mesh,0:upf%lmax,1))
272       end if
273       allocate ( vnl(1:upf%mesh) )
274       CALL xmlr_opentag( capitalize_if_v2('pp_semilocal') )
275       !
276       tag = 'vnl'
277       DO nb = 1,upf%nbeta
278          IF ( v2 ) THEN
279             ! NOTA BENE: v2 format follows available PP files, written
280             ! using original write_upf_v2; not FoX-based write_upf_v2
281             IF ( nb - 1 == upf%lloc ) CYCLE
282             tag = 'PP_VNL.'//i2c(nb-1)
283          END IF
284          CALL xmlr_readtag( tag, vnl, ierr )
285          if ( ierr /= 0 ) &
286               call upf_error('read_pp_semilocal','error reading SL PPs',1)
287          CALL get_attr ( 'l', l)
288          ind = 1
289          IF ( upf%has_so ) then
290             CALL get_attr ( 'j', j)
291             IF ( l > 0 .AND. ABS(j-l-0.5_dp) < 0.001_dp ) ind = 2
292             ! FIXME: what about spin-orbit case for v.2 upf?
293             if ( v2 ) &
294                  call upf_error('read_pp_semilocal','check spin-orbit',1)
295          END IF
296          upf%vnl(:,l,ind) = vnl(:)
297       END DO
298       deallocate ( vnl )
299       !
300       CALL xmlr_closetag( ) ! end pp_semilocal
301       !
302    END IF
303    !
304  END SUBROUTINE read_pp_semilocal
305  !
306  !--------------------------------------------------------
307  SUBROUTINE read_pp_nonlocal ( upf )
308    !--------------------------------------------------------
309    !
310    IMPLICIT NONE
311    TYPE(pseudo_upf),INTENT(INOUT) :: upf ! the pseudo data
312    !
313    LOGICAL :: isnull
314    INTEGER :: nb, ind, l, l_, ln, lm, mb, nmb
315    CHARACTER(LEN=15) :: tag
316    REAL(dp), ALLOCATABLE :: aux(:)
317    !
318    nb = upf%nbeta
319    IF ( nb == 0 ) nb = 1
320    ALLOCATE (upf%beta(upf%mesh,nb) )
321    ALLOCATE (upf%els_beta(nb), &
322              upf%lll(nb),      &
323              upf%kbeta(nb),    &
324              upf%rcut(nb),     &
325              upf%rcutus(nb),   &
326              upf%dion(nb,nb),  &
327              upf%qqq(nb,nb)    )
328    !
329    IF (upf%has_so) ALLOCATE( upf%jjj(upf%nbeta))
330    !
331    IF ( upf%nbeta == 0 ) THEN
332       upf%nqf = 0
333       upf%nqlc= 0
334       upf%kkbeta = 0
335       upf%qqq_eps=-1.0_dp
336       RETURN
337    END IF
338    !
339    CALL xmlr_opentag( capitalize_if_v2('pp_nonlocal') )
340    !
341    DO nb = 1,upf%nbeta
342       !
343       IF ( v2 ) THEN
344          tag = 'PP_BETA.'//i2c(nb)
345       ELSE
346          tag = 'pp_beta'
347       END IF
348       CALL xmlr_readtag( tag, upf%beta(1:upf%mesh,nb) )
349       CALL get_attr('index', mb)
350       ! not-so-strict test: index is absent or incorrect in some UPF v.2 files
351       IF ( .NOT. v2 .AND. nb /= mb ) &
352            CALL upf_error('read_pp_nonlocal','mismatch',nb)
353       CALL get_attr('label', upf%els_beta(nb))
354       CALL get_attr('angular_momentum', upf%lll(nb))
355       IF ( .NOT. v2 .AND. upf%has_so ) &
356            CALL get_attr('tot_ang_mom', upf%jjj(nb))
357       CALL get_attr('cutoff_radius_index', upf%kbeta(nb))
358       CALL get_attr('cutoff_radius', upf%rcut(nb))
359       CALL get_attr('ultrasoft_cutoff_radius', upf%rcutus(nb))
360       !
361    END DO
362    !
363    ! pp_dij (D_lm matrix)
364    !
365    CALL xmlr_opentag( capitalize_if_v2 ('pp_dij') )
366    READ(iun,*) upf%dion(1:upf%nbeta,1:upf%nbeta)
367    CALL xmlr_closetag( )
368    !
369    ! pp_augmentation
370    !
371    IF (upf%tvanp .or. upf%tpawp) THEN
372       CALL xmlr_opentag( capitalize_if_v2('pp_augmentation') )
373       !
374       IF ( v2 ) THEN
375          CALL get_attr ( 'q_with_l', upf%q_with_l )
376          CALL get_attr ( 'nqf', upf%nqf )
377          CALL get_attr ( 'nqlc', upf%nqlc )
378          IF (upf%tpawp) THEN
379             CALL get_attr ( 'shape', upf%paw%augshape )
380             CALL get_attr ( 'cutoff_r', upf%paw%raug )
381             CALL get_attr ( 'cutoff_r_index', upf%paw%iraug )
382             CALL get_attr ( 'augmentation_epsilon', upf%qqq_eps )
383             CALL get_attr ( 'l_max_aug', upf%paw%lmax_aug )
384          ENDIF
385       ELSE
386          CALL xmlr_readtag( 'q_with_l', upf%q_with_l )
387          CALL xmlr_readtag( 'nqf', upf%nqf )
388          CALL xmlr_readtag( 'nqlc', upf%nqlc )
389          IF (upf%tpawp) THEN
390             CALL xmlr_readtag( 'shape', upf%paw%augshape )
391             CALL xmlr_readtag( 'cutoff_r', upf%paw%raug )
392             CALL xmlr_readtag( 'cutoff_r_index', upf%paw%iraug )
393             CALL xmlr_readtag( 'augmentation_epsilon', upf%qqq_eps )
394             CALL xmlr_readtag( 'l_max_aug', upf%paw%lmax_aug )
395          ENDIF
396       ENDIF
397       !
398       CALL xmlr_opentag( capitalize_if_v2('pp_q') )
399       READ(iun,*) upf%qqq(1:upf%nbeta,1:upf%nbeta)
400       CALL xmlr_closetag( )
401       !
402       IF ( upf%tpawp ) THEN
403          CALL xmlr_opentag( capitalize_if_v2('pp_multipoles') )
404          ALLOCATE ( upf%paw%augmom(1:upf%nbeta,1:upf%nbeta,0:2*upf%lmax) )
405          READ(iun,*) upf%paw%augmom(1:upf%nbeta,1:upf%nbeta,0:2*upf%lmax)
406          CALL xmlr_closetag ()
407       ENDIF
408       !
409       ! read polinomial coefficients for Q_ij expansion at small radius
410       !
411       IF ( upf%nqlc == 0 ) upf%nqlc = 2*upf%lmax+1
412       ALLOCATE( upf%rinner( upf%nqlc ) )
413       IF ( v2 .AND. upf%nqf > 0) THEN
414          ALLOCATE ( upf%qfcoef(upf%nqf, upf%nqlc, upf%nbeta, upf%nbeta) )
415          CALL xmlr_opentag('PP_QFCOEF')
416          READ(iun,*) upf%qfcoef
417          CALL xmlr_closetag ()
418          CALL xmlr_readtag('PP_RINNER',upf%rinner)
419       ELSE IF ( upf%nqf == 0 ) THEN
420          ALLOCATE( upf%qfcoef(1,1,1,1) )
421          upf%qfcoef =0.0_dp
422       ENDIF
423       !
424       ! Read augmentation charge Q_ij
425       !
426       IF( upf%q_with_l ) THEN
427          ALLOCATE( upf%qfuncl(upf%mesh,upf%nbeta*(upf%nbeta+1)/2,0:2*upf%lmax) )
428          upf%qfuncl(:,:,:) = 0.0_dp
429          ! NOTE: it would be wiser to dimension qfuncl as (:,:,0:upf%lmax)
430          ! and store the q_l(r) with index l=L/2 (see loop_on_l below)
431          ! This would save some storage and avoid "holes" in the array
432          ! that may be a source of trouble if not initialized to zero
433       ELSE
434          ALLOCATE ( upf%qfunc(upf%mesh,upf%nbeta*(upf%nbeta+1)/2) )
435          upf%qfunc (:,:) = 0.0_dp
436       END IF
437       ALLOCATE ( aux(upf%mesh) )
438       loop_on_nb: DO nb = 1,upf%nbeta
439          ln = upf%lll(nb)
440          loop_on_mb: DO mb = nb,upf%nbeta
441             lm = upf%lll(mb)
442             IF( upf%q_with_l ) THEN
443                loop_on_l: DO l = abs(ln-lm),ln+lm,2 ! only even terms
444                   isnull = .FALSE.
445                   IF( upf%tpawp ) isnull = (abs(upf%paw%augmom(nb,mb,l)) < upf%qqq_eps)
446                   IF(isnull) CYCLE loop_on_l
447                   IF ( v2 ) THEN
448                      tag = 'PP_QIJL.'//i2c(nb)//'.'//i2c(mb)//'.'//i2c(l)
449                   ELSE
450                      tag = 'pp_qijl'
451                   END IF
452                   CALL xmlr_readtag( tag, aux )
453                   CALL get_attr ('composite_index', nmb)
454                   IF ( nmb /= mb*(mb-1)/2 + nb ) &
455                        CALL upf_error ('read_pp_nonlocal','mismatch',1)
456                   CALL get_attr ('angular_momentum', l_)
457                   IF ( l /= l_ ) CALL upf_error ('read_pp_nonlocal','mismatch',2)
458                   upf%qfuncl(:,nmb,l) = aux(:)
459                   IF (upf%tpawp) upf%qfuncl(upf%paw%iraug+1:,nmb,l) = 0._DP
460                ENDDO loop_on_l
461             ELSE
462                isnull = .FALSE.
463                IF  ( upf%tpawp ) isnull = ( abs(upf%qqq(nb,mb)) < upf%qqq_eps )
464                IF (isnull) CYCLE loop_on_mb
465                IF ( v2 ) THEN
466                   tag = 'PP_QIJ.'//i2c(nb)//'.'//i2c(mb)
467                ELSE
468                   tag = 'pp_qij'
469                END IF
470                CALL xmlr_readtag( tag, aux )
471                CALL get_attr ('composite_index', nmb)
472                IF ( nmb /= mb*(mb-1)/2 + nb ) &
473                     CALL upf_error ('read_pp_nonlocal','mismatch',3)
474                upf%qfunc(:,nmb) = aux(:)
475                !
476             ENDIF
477          ENDDO loop_on_mb
478       ENDDO  loop_on_nb
479       !
480       DEALLOCATE (aux)
481       CALL xmlr_closetag( ) ! end pp_augmentation
482       !
483    END IF
484    CALL xmlr_closetag( ) ! end pp_nonlocal
485    !
486    ! Maximum radius of beta projector: outer radius to integrate
487    upf%kkbeta = MAXVAL(upf%kbeta(1:upf%nbeta))
488    ! For PAW, augmentation charge may extend a bit further:
489    IF(upf%tpawp) upf%kkbeta = MAX(upf%kkbeta, upf%paw%iraug)
490    !
491  END SUBROUTINE read_pp_nonlocal
492  !
493  !--------------------------------------------------------
494  SUBROUTINE read_pp_pswfc ( upf )
495    !--------------------------------------------------------
496    !
497    IMPLICIT NONE
498    TYPE(pseudo_upf),INTENT(INOUT) :: upf ! the pseudo data
499    !
500    INTEGER :: nw, ind, l
501    CHARACTER(LEN=8) :: tag
502    !
503    allocate ( upf%chi(1:upf%mesh,upf%nwfc) )
504    allocate ( upf%els(upf%nwfc), &
505                upf%oc(upf%nwfc), &
506                upf%lchi(upf%nwfc), &
507                upf%nchi(upf%nwfc), &
508                upf%rcut_chi(upf%nwfc), &
509                upf%rcutus_chi(upf%nwfc), &
510                upf%epseu(upf%nwfc) )
511    IF ( upf%has_so ) THEN
512       allocate ( upf%nn(upf%nwfc) )
513       allocate ( upf%jchi(upf%nwfc) )
514    END IF
515    !
516    CALL xmlr_opentag( capitalize_if_v2('pp_pswfc') )
517    DO nw=1,upf%nwfc
518       IF ( v2 ) THEN
519          tag = 'PP_CHI.'//i2c(nw)
520       ELSE
521          tag = 'pp_chi'
522       END IF
523       CALL xmlr_readtag( tag, upf%chi(1:upf%mesh,nw) )
524       call get_attr('index', ind)
525       ! not-so-strict test: index is absent or incorrect in some UPF v.2 files
526       if ( .NOT. v2 .AND. ind /= nw ) &
527            call upf_error('read_pp_pswfc','mismatch reading PSWFC', nw)
528       call get_attr( 'label', upf%els(nw) )
529       call get_attr( 'l', upf%lchi(nw) )
530       IF ( .not. v2 .and. upf%has_so ) THEN
531          call get_attr( 'nn', upf%nn(nw) )
532          call get_attr( 'jchi', upf%jchi(nw) )
533       END IF
534       call get_attr( 'occupation', upf%oc(nw) )
535       call get_attr( 'n', upf%nchi(nw) )
536       call get_attr( 'pseudo_energy', upf%epseu(nw) )
537       call get_attr( 'cutoff_radius', upf%rcut_chi(nw) )
538       call get_attr( 'ultrasoft_cutoff_radius', upf%rcutus_chi(nw) )
539    END DO
540    CALL xmlr_closetag( ) ! end pp_pswfc
541    !
542  END SUBROUTINE read_pp_pswfc
543  !
544  !--------------------------------------------------------
545  SUBROUTINE read_pp_full_wfc ( upf )
546    !--------------------------------------------------------
547    !
548    IMPLICIT NONE
549    TYPE(pseudo_upf),INTENT(INOUT) :: upf ! the pseudo data
550    !
551    INTEGER :: nb, mb
552    CHARACTER(LEN=15) :: tag
553    !
554    IF ( upf%has_wfc ) THEN
555       !
556       ALLOCATE (upf%aewfc(1:upf%mesh,upf%nbeta) )
557       CALL xmlr_opentag( capitalize_if_v2('pp_full_wfc') )
558       !
559       DO nb = 1, upf%nbeta
560          IF ( v2 ) THEN
561             tag = 'PP_AEWFC.'//i2c(nb)
562          ELSE
563             tag = 'pp_aewfc'
564          END IF
565          CALL xmlr_readtag( tag, upf%aewfc(1:upf%mesh,nb) )
566          CALL get_attr ('index',mb)
567          ! not-so-strict test (and two more below):
568          ! index may be absent or incorrect in some UPF v.2 files
569          IF ( .NOT. v2 .AND. nb /= mb ) CALL upf_error('read_pp_full_wfc','mismatch',1)
570       END DO
571       !
572       IF ( upf%has_so .AND. upf%tpawp ) THEN
573          ALLOCATE (upf%paw%aewfc_rel(1:upf%mesh,upf%nbeta) )
574          DO nb = 1, upf%nbeta
575             IF ( v2 ) THEN
576                tag = 'PP_AEWFC_rel.'//i2c(nb)
577             ELSE
578                tag = 'pp_aewfc_rel'
579             END IF
580             CALL xmlr_readtag(tag, upf%paw%aewfc_rel(1:upf%mesh,nb) )
581             CALL get_attr ('index',mb)
582             IF ( .NOT. v2 .AND. nb /= mb ) CALL upf_error('read_pp_full_wfc','mismatch',2)
583          END DO
584       END IF
585       !
586       ALLOCATE (upf%pswfc(1:upf%mesh,upf%nbeta) )
587       DO nb = 1, upf%nbeta
588          IF ( v2 ) THEN
589             tag = 'PP_PSWFC.'//i2c(nb)
590          ELSE
591             tag = 'pp_pswfc'
592          END IF
593          CALL xmlr_readtag(tag, upf%pswfc(1:upf%mesh,nb) )
594          CALL get_attr ('index',mb)
595          IF ( .NOT. v2 .AND. nb /= mb ) CALL upf_error('read_pp_full_wfc','mismatch',3)
596       END DO
597       !
598       CALL xmlr_closetag( )
599       !
600    END IF
601    !
602  END SUBROUTINE read_pp_full_wfc
603  !
604  !--------------------------------------------------------
605  SUBROUTINE read_pp_spinorb ( upf )
606    !--------------------------------------------------------
607    !
608    IMPLICIT NONE
609    TYPE(pseudo_upf),INTENT(INOUT) :: upf ! the pseudo data
610    INTEGER :: nw, nb, ierr
611    CHARACTER(LEN=1) :: dummy
612    !
613    IF ( .NOT. v2 .OR. .NOT. upf%has_so ) RETURN
614    !
615    CALL xmlr_opentag( 'PP_SPIN_ORB' )
616    DO nw = 1,upf%nwfc
617       CALL xmlr_readtag( 'PP_RELWFC.'//i2c(nw), dummy )
618       CALL get_attr( 'index' , nb )
619       ! not-so-strict test: index absent or incorrect in some UPF v.2 files
620       IF ( .NOT. v2 .AND. nb /= nw ) CALL upf_error('read_pp_spinorb','mismatch',1)
621       CALL get_attr( 'els',   upf%els(nw) )
622       CALL get_attr( 'nn',    upf%nn(nw) )
623       CALL get_attr( 'lchi',  upf%lchi(nw) )
624       CALL get_attr( 'jchi',  upf%jchi(nw) )
625       CALL get_attr( 'oc',    upf%oc(nw) )
626    ENDDO
627    !
628    DO nb = 1,upf%nbeta
629       CALL xmlr_readtag( 'PP_RELBETA.'//i2c(nb), dummy, ierr )
630       !
631       ! existing PP files may have pp_relbeta first, pp_relwfc later,
632       ! but also the other way round - check that everything was right
633       !
634       if ( ierr /= 0 ) then
635          ierr = -81
636          return
637       end if
638       CALL get_attr( 'index' , nw )
639       IF ( .NOT.v2 .AND. nb /= nw ) CALL upf_error('read_pp_spinorb','mismatch',2)
640       CALL get_attr( 'lll',  upf%lll(nb) )
641       CALL get_attr( 'jjj',  upf%jjj(nb) )
642    ENDDO
643    CALL xmlr_closetag () ! end pp_spin_orb
644    !
645  END SUBROUTINE read_pp_spinorb
646  !
647  !--------------------------------------------------------
648  SUBROUTINE read_pp_paw ( upf )
649    !--------------------------------------------------------
650    !
651    IMPLICIT NONE
652    TYPE(pseudo_upf),INTENT(INOUT) :: upf ! the pseudo data
653    INTEGER :: nb, mb
654    !
655    IF ( .NOT. upf%tpawp ) RETURN
656    !
657    CALL xmlr_opentag( capitalize_if_v2('pp_paw') )
658    CALL get_attr ('paw_data_format', upf%paw_data_format)
659    CALL get_attr ('core_energy', upf%paw%core_energy)
660    ! Full occupation (not only > 0 ones)
661    ALLOCATE (upf%paw%oc(upf%nbeta) )
662    ALLOCATE (upf%paw%ae_rho_atc(upf%mesh) )
663    ALLOCATE (upf%paw%ae_vloc(upf%mesh) )
664    CALL xmlr_readtag( capitalize_if_v2('pp_occupations'), &
665         upf%paw%oc(1:upf%nbeta) )
666    ! All-electron core charge
667    CALL xmlr_readtag( capitalize_if_v2('pp_ae_nlcc'), &
668         upf%paw%ae_rho_atc(1:upf%mesh) )
669    ! All-electron local potential
670    CALL xmlr_readtag( capitalize_if_v2('pp_ae_vloc'), &
671         upf%paw%ae_vloc(1:upf%mesh) )
672    CALL xmlr_closetag () ! end pp_paw
673    !
674    ALLOCATE(upf%paw%pfunc(upf%mesh, upf%nbeta,upf%nbeta) )
675    upf%paw%pfunc(:,:,:) = 0._dp
676    IF (upf%has_so) THEN
677       ALLOCATE(upf%paw%pfunc_rel(upf%mesh, upf%nbeta,upf%nbeta) )
678       upf%paw%pfunc_rel(:,:,:) = 0._dp
679    ENDIF
680    DO nb=1,upf%nbeta
681       DO mb=1,nb
682          upf%paw%pfunc (1:upf%mesh, nb, mb) = &
683               upf%aewfc(1:upf%mesh, nb) * upf%aewfc(1:upf%mesh, mb)
684          IF (upf%has_so) THEN
685             upf%paw%pfunc_rel (1:upf%paw%iraug, nb, mb) =  &
686                  upf%paw%aewfc_rel(1:upf%paw%iraug, nb) *   &
687                  upf%paw%aewfc_rel(1:upf%paw%iraug, mb)
688!
689!    The small component is added to pfunc. pfunc_rel is useful only
690!    to add a small magnetic contribution
691!
692             upf%paw%pfunc (1:upf%paw%iraug, nb, mb) = &
693                        upf%paw%pfunc (1:upf%paw%iraug, nb, mb) + &
694                        upf%paw%pfunc_rel (1:upf%paw%iraug, nb, mb)
695          ENDIF
696          upf%paw%pfunc(upf%paw%iraug+1:,nb,mb) = 0._dp
697          !
698          upf%paw%pfunc (1:upf%mesh, mb, nb) = upf%paw%pfunc (1:upf%mesh, nb, mb)
699          IF (upf%has_so) upf%paw%pfunc_rel (1:upf%mesh, mb, nb) =  &
700               upf%paw%pfunc_rel (1:upf%mesh, nb, mb)
701       ENDDO
702    ENDDO
703    !
704    ! Pseudo wavefunctions (not only the ones for oc > 0)
705    ! All-electron wavefunctions
706    ALLOCATE(upf%paw%ptfunc(upf%mesh, upf%nbeta,upf%nbeta) )
707    upf%paw%ptfunc(:,:,:) = 0._dp
708    DO nb=1,upf%nbeta
709       DO mb=1,upf%nbeta
710          upf%paw%ptfunc (1:upf%mesh, nb, mb) = &
711               upf%pswfc(1:upf%mesh, nb) * upf%pswfc(1:upf%mesh, mb)
712          upf%paw%ptfunc(upf%paw%iraug+1:,nb,mb) = 0._dp
713          !
714          upf%paw%ptfunc (1:upf%mesh, mb, nb) = upf%paw%ptfunc (1:upf%mesh, nb, mb)
715       ENDDO
716    ENDDO
717    !
718  END SUBROUTINE read_pp_paw
719  !--------------------------------------------------------
720  SUBROUTINE read_pp_gipaw ( upf )
721    !--------------------------------------------------------
722    !
723    IMPLICIT NONE
724    TYPE(pseudo_upf),INTENT(INOUT) :: upf ! the pseudo data
725    !
726    INTEGER :: nb, mb
727    CHARACTER(LEN=24) :: tag
728    !
729    IF (.NOT. upf%has_gipaw) RETURN
730    !
731    CALL xmlr_opentag( capitalize_if_v2('pp_gipaw') )
732    CALL get_attr ('gipaw_data_format', upf%gipaw_data_format )
733    IF ( v2 ) THEN
734       CALL xmlr_opentag( 'PP_GIPAW_CORE_ORBITALS')
735       CALL get_attr ('number_of_core_orbitals', upf%gipaw_ncore_orbitals)
736    ELSE
737       CALL xmlr_readtag ('number_of_core_orbitals', upf%gipaw_ncore_orbitals)
738       IF ( .NOT. upf%paw_as_gipaw) &
739          CALL xmlr_readtag( 'number_of_valence_orbitals', upf%gipaw_wfs_nchannels)
740    END IF
741    ALLOCATE ( upf%gipaw_core_orbital(upf%mesh,upf%gipaw_ncore_orbitals) )
742    ALLOCATE ( upf%gipaw_core_orbital_n(upf%gipaw_ncore_orbitals) )
743    ALLOCATE ( upf%gipaw_core_orbital_el(upf%gipaw_ncore_orbitals) )
744    ALLOCATE ( upf%gipaw_core_orbital_l(upf%gipaw_ncore_orbitals) )
745    DO nb = 1,upf%gipaw_ncore_orbitals
746       IF ( v2 ) THEN
747          tag = "PP_GIPAW_CORE_ORBITAL."//i2c(nb)
748       ELSE
749          tag = 'pp_gipaw_core_orbital'
750       END IF
751       CALL xmlr_readtag( tag, upf%gipaw_core_orbital(1:upf%mesh,nb) )
752       CALL get_attr ('index', mb)
753       IF ( nb /= mb ) CALL upf_error('read_pp_gipaw','mismatch',1)
754       CALL get_attr ('label', upf%gipaw_core_orbital_el(nb) )
755       CALL get_attr ('n', upf%gipaw_core_orbital_n(nb) )
756       CALL get_attr ('l', upf%gipaw_core_orbital_l(nb) )
757    END DO
758    IF ( v2 ) CALL xmlr_closetag ( )
759    !
760    IF ( upf%paw_as_gipaw) THEN
761       !
762       !    PAW as GIPAW case: all-electron and pseudo-orbitals not read here
763       !
764       upf%gipaw_wfs_nchannels = upf%nbeta
765       ALLOCATE ( upf%gipaw_wfs_el(upf%gipaw_wfs_nchannels) )
766       ALLOCATE ( upf%gipaw_wfs_ll(upf%gipaw_wfs_nchannels) )
767       ALLOCATE ( upf%gipaw_wfs_rcut(upf%gipaw_wfs_nchannels) )
768       ALLOCATE ( upf%gipaw_wfs_rcutus(upf%gipaw_wfs_nchannels) )
769       ALLOCATE ( upf%gipaw_wfs_ae(upf%mesh,upf%gipaw_wfs_nchannels) )
770       ALLOCATE ( upf%gipaw_wfs_ps(upf%mesh,upf%gipaw_wfs_nchannels) )
771       DO nb = 1,upf%gipaw_wfs_nchannels
772          upf%gipaw_wfs_el(nb) = upf%els_beta(nb)
773          upf%gipaw_wfs_ll(nb) = upf%lll(nb)
774          upf%gipaw_wfs_ae(:,nb) = upf%aewfc(:,nb)
775       ENDDO
776       DO nb = 1,upf%gipaw_wfs_nchannels
777          upf%gipaw_wfs_ps(:,nb) = upf%pswfc(:,nb)
778       ENDDO
779       ALLOCATE ( upf%gipaw_vlocal_ae(upf%mesh) )
780       ALLOCATE ( upf%gipaw_vlocal_ps(upf%mesh) )
781       upf%gipaw_vlocal_ae(:)= upf%paw%ae_vloc(:)
782       upf%gipaw_vlocal_ps(:)= upf%vloc(:)
783       DO nb = 1,upf%gipaw_wfs_nchannels
784          upf%gipaw_wfs_rcut(nb)=upf%rcut(nb)
785          upf%gipaw_wfs_rcutus(nb)=upf%rcutus(nb)
786       ENDDO
787       !
788    ELSE
789       !
790       ! Read valence all-electron and pseudo orbitals
791       !
792       IF ( v2 ) THEN
793          CALL xmlr_opentag( 'PP_GIPAW_ORBITALS' )
794          CALL get_attr( 'number_of_valence_orbitals', &
795               upf%gipaw_wfs_nchannels )
796       END IF
797       ALLOCATE ( upf%gipaw_wfs_el(upf%gipaw_wfs_nchannels) )
798       ALLOCATE ( upf%gipaw_wfs_ll(upf%gipaw_wfs_nchannels) )
799       ALLOCATE ( upf%gipaw_wfs_rcut(upf%gipaw_wfs_nchannels) )
800       ALLOCATE ( upf%gipaw_wfs_rcutus(upf%gipaw_wfs_nchannels) )
801       ALLOCATE ( upf%gipaw_wfs_ae(upf%mesh,upf%gipaw_wfs_nchannels) )
802       ALLOCATE ( upf%gipaw_wfs_ps(upf%mesh,upf%gipaw_wfs_nchannels) )
803       ALLOCATE ( upf%gipaw_vlocal_ae(upf%mesh) )
804       ALLOCATE ( upf%gipaw_vlocal_ps(upf%mesh) )
805       DO nb = 1,upf%gipaw_wfs_nchannels
806          IF ( v2 ) THEN
807             tag = "PP_GIPAW_ORBITAL."//i2c(nb)
808          ELSE
809             tag = 'pp_gipaw_orbital'
810          END IF
811          CALL xmlr_opentag( tag )
812          CALL get_attr ('index', mb)
813          IF ( nb /= mb ) CALL upf_error('read_pp_gipaw','mismatch',2)
814          CALL get_attr ('label', upf%gipaw_wfs_el(nb) )
815          CALL get_attr ('l',     upf%gipaw_wfs_ll(nb) )
816          CALL get_attr ('cutoff_radius', upf%gipaw_wfs_rcut(nb) )
817          CALL get_attr ('ultrasoft_cutoff_radius', upf%gipaw_wfs_rcutus(nb) )
818          CALL xmlr_readtag( capitalize_if_v2('pp_gipaw_wfs_ae'), &
819               upf%gipaw_wfs_ae(1:upf%mesh,nb) )
820          CALL xmlr_readtag( capitalize_if_v2('pp_gipaw_wfs_ps'),&
821               upf%gipaw_wfs_ps(1:upf%mesh,nb) )
822          CALL xmlr_closetag ()
823       END DO
824       IF ( v2 ) CALL xmlr_closetag( )
825       !
826       ! Read all-electron and pseudo local potentials
827       !
828       CALL xmlr_opentag( capitalize_if_v2('pp_gipaw_vlocal') )
829       CALL xmlr_readtag( capitalize_if_v2('pp_gipaw_vlocal_ae'), &
830            upf%gipaw_vlocal_ae(1:upf%mesh) )
831       CALL xmlr_readtag( capitalize_if_v2('pp_gipaw_vlocal_ps'), &
832            upf%gipaw_vlocal_ps(1:upf%mesh) )
833       CALL xmlr_closetag ()
834    END IF
835    CALL xmlr_closetag () ! end pp_gipaw
836    !
837  END SUBROUTINE read_pp_gipaw
838  !
839END MODULE read_upf_new_module
840