[skip-CI] V.-1 of PSML reading

This commit is contained in:
Paolo Giannozzi 2023-04-13 18:11:50 +02:00
parent 2d0190377d
commit 889ecb4a9a
5 changed files with 469 additions and 38 deletions

View File

@ -35,6 +35,7 @@ read_cpmd.o \
read_fhi.o \
read_ncpp.o \
read_ps.o \
read_psml.o \
read_upf_new.o \
read_upf_v1.o \
read_uspp.o \

View File

@ -50,8 +50,6 @@ SUBROUTINE read_ps ( filein, upf_in )
END IF
!
ELSE
!
OPEN ( UNIT=iunps, FILE=filein, STATUS = 'old', FORM = 'formatted' )
!
! The type of the pseudopotential is determined by the file name:
! *.xml or *.XML UPF format with schema pp_format=0
@ -61,8 +59,21 @@ SUBROUTINE read_ps ( filein, upf_in )
! *.RRKJ3 PWSCF US PP format ("atomic" code) pp_format=4
! *.fhi or *cpi FHI/Abinit numerical NC PP pp_format=6
! *.cpmd CPMD NC pseudopot. file format pp_format=7
! *.psml Siesta/Abinit psml NCPP file format pp_format=8
! none of the above: PWSCF norm-conserving format pp_format=5
!
IF ( upf_get_pp_format( filein ) == 8 ) THEN
!
! Unlike all other cases, file must be opened with xml tools
!
WRITE( stdout, "('file type is PSML (experimental)')")
CALL read_psml (filein, upf_in)
RETURN
!
END IF
!
OPEN ( UNIT=iunps, FILE=filein, STATUS = 'old', FORM = 'formatted' )
!
IF ( upf_get_pp_format( filein ) == 2 ) THEN
!
WRITE( stdout, "('file type is Vanderbilt US PP')")
@ -93,6 +104,7 @@ SUBROUTINE read_ps ( filein, upf_in )
WRITE( stdout, "('file type is CPMD NC format')")
CALL read_cpmd (iunps, upf_in)
!
!
ELSE
!
CALL upf_error('readpp', 'file '//TRIM(filein)//' not readable',1)

409
upflib/read_psml.f90 Normal file
View File

@ -0,0 +1,409 @@
!
! Copyright (C) 2023 Quantum ESPRESSO Foundation
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------
subroutine read_psml ( filename, upf )
!-----------------------------------------------------
!! Read pseudopotential files in PSML format using "xmltools"
!! stores data into the "upf" structure
!! Written by P. Giannozzi, April 2023
!
USE xmltools
USE pseudo_types, ONLY: pseudo_upf
!
IMPLICIT NONE
INTEGER, PARAMETER :: dp = selected_real_kind(14,200)
REAL, PARAMETER :: e2 = 2.0_dp
!
CHARACTER(len=*), INTENT(IN) :: filename
!! input : name of file in psml format
TYPE(pseudo_upf), INTENT(OUT) :: upf
!! the derived type storing the pseudo data
INTEGER :: ierr
!! error code (0 if correctly read)
CHARACTER(len=30) :: tag
!! tag where error (ierr != 0) was detected
INTEGER :: iun
!! unit for reading data
!
ierr = 0
iun = xml_open_file ( filename )
IF ( iun == -1 ) THEN
ierr = 1
tag = 'file'
GO TO 10
END IF
!
tag = 'psml'
call xmlr_opentag ( trim(tag), IERR = ierr )
IF ( ierr /= 0 ) RETURN
call get_attr ( 'version', upf%nv ); print *, 'version=',upf%nv
call get_attr ( 'uuid', upf%author )
!
tag = 'provenance'
call xmlr_opentag ( trim(tag), IERR = ierr )
IF ( ierr /= 0 ) GO TO 10
call get_attr ( 'creator', upf%generated )
upf%author='UUID: '//trim(upf%author)
call get_attr ( 'date', upf%date )
upf%comment='PSML file v. '//trim(upf%nv)
call xmlr_closetag ( ) ! provenance'
!
tag = 'pseudo-atom-spec'
call read_psml_pseudo_atom_spec ( ierr )
IF ( ierr /= 0 ) GO TO 10
!
tag = 'grid'
call read_psml_grid ( ierr )
IF ( ierr /= 0 ) GO TO 10
!
tag = 'valence-charge'
call read_psml_radialfunc ( tag, upf%rho_at,ierr )
IF ( ierr /= 0 ) GO TO 10
!
IF ( upf%nlcc ) THEN
tag = 'pseudocore-charge'
call read_psml_radialfunc ( tag, upf%rho_atc,ierr )
IF ( ierr /= 0 ) GO TO 10
END IF
!
tag = 'local-potential'
call read_psml_radialfunc ( tag, upf%vloc, ierr )
IF ( ierr /= 0 ) GO TO 10
!
tag = 'nonlocal-projectors'
call read_psml_nonlocal_projectors ( ierr )
IF ( ierr /= 0 ) GO TO 10
!
tag = 'pseudo-wave-functions'
call read_psml_pseudo_wave_functions ( ierr )
IF ( ierr /= 0 ) THEN
print *, 'read_psml: tag ',trim(tag),' not present'
upf%nwfc = 0
END IF
!
call xmlr_closetag ( ) ! psml
call xml_closefile( )
!
! Convert from Hartree (PSML) to Ry
!
upf%vloc(:) = e2*upf%vloc(:)
upf%beta(:,:) = e2*upf%beta(:,:)
upf%dion(:,:) = upf%dion(:,:)/e2
!
RETURN
10 print *, 'read_psml: error reading tag ',trim(tag)
stop
!
CONTAINS
!
SUBROUTINE read_psml_pseudo_atom_spec ( ierr )
!
INTEGER :: ierr
INTEGER :: n, nxc, ndum
INTEGER :: xc(6)
CHARACTER(len=3) :: cc
CHARACTER(len=25), external :: libxc_to_qe
!
upf%tvanp = .false.
upf%tpawp = .false.
upf%has_so = .false.
upf%has_gipaw = .false.
upf%tcoulombp = .false.
upf%is_gth = .false.
upf%is_multiproj = .false.
call xmlr_opentag ( 'pseudo-atom-spec', IERR = ierr )
if (ierr /= 0) return
call get_attr ( 'atomic-label', upf%psd )
call get_attr ( 'atomic-number', upf%zmesh )
call get_attr ( 'relativity', upf%rel )
if ( upf%rel(1:5) == 'dirac' ) upf%rel='full'
call get_attr ( 'core-corrections', cc )
upf%nlcc = (cc == 'yes')
call xmlr_opentag ( 'exchange-correlation', IERR = ierr )
if (ierr /= 0) return
call xmlr_opentag ( 'libxc-info', IERR = ierr )
if (ierr /= 0) return
call get_attr ( 'number-of-functionals', nxc )
do n=1,nxc
call xmlr_readtag ( 'functional', xc(n), IERR = ierr )
if (ierr > 0) return
call get_attr ( 'id', xc(n) )
end do
call xmlr_closetag ( ) ! libxc-info
call xmlr_closetag ( ) ! exchange-correlation
!
upf%dft = libxc_to_qe (nxc, xc)
!
call xmlr_opentag ( 'valence-configuration', IERR = ierr )
if (ierr /= 0) return
call get_attr ( 'total-valence-charge', upf%zp )
! here just count the number of valence wavefunctions
n=0
do
call xmlr_readtag ( 'shell', cc, IERR = ierr )
if (ierr /= -1) exit
n = n+1
!call get_attr ( 'n', ndum )
!call get_attr ( 'l', cc )
!call get_attr ( 'occupation', ndum )
end do
upf%nwfc = n
call xmlr_closetag ( ) ! valence-configuration
call xmlr_closetag ( ) ! pseudo-atom-spec
!
ierr = 0
!
END SUBROUTINE read_psml_pseudo_atom_spec
!
SUBROUTINE read_psml_grid ( ierr )
!
INTEGER :: ierr, i, n
CHARACTER(LEN=1) :: dum
REAL(dp) :: step, delta, r0, r1, r2
!
call xmlr_opentag ( 'grid', IERR = ierr )
if (ierr /= 0) return
call get_attr ( 'npts', upf%mesh )
allocate (upf%r(upf%mesh))
call xmlr_readtag ( 'annotation', dum, IERR = ierr )
if (ierr > 0) return
call get_attr ( 'step', step )
call get_attr ( 'scale', r0 )
call get_attr ( 'delta', delta )
call xmlr_readtag ( 'grid-data', upf%r, IERR = ierr )
if (ierr /= 0) return
call xmlr_closetag ( ) ! grid
! compute, or more exactly guess, suitable dr/dx for integration
allocate (upf%rab(upf%mesh))
upf%rab(:) = 0.0_dp
! i is the index of the original logarithmic grid
! n is the index of the PSML grid
i = 0
n = 0
r1= 0.0_dp
do
i = i + 1
r2 = r0*exp((i-1)*step)
if ( r2-r1 > delta) then
! points at distance > delta are those selected
n = n+1
! check: is this the correct grid?
! write (100,*) n, i, upf%r(n+1), r2
if ( abs(r2-upf%r(n+1)) > 1.e-8 ) then
print *, 'Unknown grid !'
return
end if
upf%rab(n+1) = (r2-r1)
r1 = r2
end if
if ( n == upf%mesh -1 ) exit
end do
upf%rab(1) = delta
!
END SUBROUTINE read_psml_grid
!
SUBROUTINE read_psml_radialfunc ( tag, rho, ierr )
!
CHARACTER(len=*) :: tag
REAL(dp), allocatable :: rho(:)
INTEGER :: ierr
INTEGER :: npt
!
call xmlr_opentag ( trim(tag), IERR = ierr )
if (ierr /= 0) return
call xmlr_opentag ( 'radfunc', IERR = ierr )
if (ierr /= 0) return
call xmlr_opentag ( 'data', IERR = ierr )
if (ierr /= 0) return
call get_attr ( 'npts', npt )
! may differ
if ( npt > upf%mesh ) then
ierr = 1
return
end if
allocate ( rho(upf%mesh) )
read (iun,*) rho(1:npt)
if ( npt < upf%mesh ) rho(npt+1:) = 0.0_dp
call xmlr_closetag ( ) ! data
call xmlr_closetag ( ) ! radfunc
call xmlr_closetag ( ) ! tag
!
END SUBROUTINE read_psml_radialfunc
!
SUBROUTINE read_psml_nonlocal_projectors ( ierr )
!
INTEGER :: ierr
INTEGER :: nb, npt, ndum
real(dp) :: ekb
CHARACTER(len=1) :: spdf
!
call xmlr_opentag('nonlocal-projectors', IERR = ierr )
if (ierr /= 0) return
upf%nbeta = 0
upf%kkbeta= 0
nb=0
do
call xmlr_opentag ( 'proj', IERR = ierr )
call get_attr ( 'l', spdf )
! call get_attr ( 'seq',ndum )
call get_attr ( 'ekb',ekb )
call xmlr_opentag ( 'radfunc' )
if ( ierr == -10 .and. upf%nbeta == 0 ) then
! first scan of the file completed, file has been rewound:
! set number of projectors, allocate and read arrays
upf%nbeta = nb
ALLOCATE (upf%els_beta(nb), &
upf%lll(nb), &
upf%kbeta(nb), &
upf%rcut(nb), &
upf%rcutus(nb), &
upf%dion(nb,nb), &
upf%qqq(nb,nb) )
allocate (upf%beta(upf%mesh,nb))
upf%rcut(:) = 0.0_dp
upf%rcutus(:)= 0.0_dp
upf%dion(:,:)= 0.0_dp
upf%qqq(:,:) = 0.0_dp
! reset counter
nb = 0
end if
nb = nb+1
call xmlr_opentag ( 'data', IERR = ierr )
call get_attr ( 'npts', npt )
if ( upf%nbeta > 0 ) then
! actual read is done here during the second scan
read (iun,*) upf%beta(1:npt,nb)
if ( npt < upf%mesh) upf%beta(npt+1:,nb) = 0.0_dp
upf%dion(nb,nb) = 1.0_dp/ekb
upf%els_beta(nb) = '*'//spdf
if ( spdf == 's' .or. spdf == 'S' ) then
upf%lll(nb) = 0
else if ( spdf == 'p' .or. spdf == 'P' ) then
upf%lll(nb) = 1
else if ( spdf == 'd' .or. spdf == 'D' ) then
upf%lll(nb) = 2
else if ( spdf == 'f' .or. spdf == 'F' ) then
upf%lll(nb) = 3
end if
upf%kbeta(nb) = npt
else
! set max length of projectors during the first scan
upf%kkbeta = max ( upf%kkbeta, npt )
end if
call xmlr_closetag () ! data
call xmlr_closetag () ! radfun
call xmlr_closetag () ! proj
if ( nb == upf%nbeta ) exit
end do
call xmlr_closetag ( ) ! nonlocal-projectors
!
END SUBROUTINE read_psml_nonlocal_projectors
!
SUBROUTINE read_psml_pseudo_wave_functions ( ierr )
!
INTEGER :: ierr
INTEGER :: n, npt, ndum
CHARACTER(len=1) :: spdf
!
call xmlr_opentag('pseudo-wave-functions', IERR = ierr )
if (ierr /= 0) return
allocate ( upf%chi(upf%mesh,upf%nwfc) )
allocate ( upf%els(upf%nwfc), &
upf%oc(upf%nwfc), &
upf%lchi(upf%nwfc), &
upf%nchi(upf%nwfc), &
upf%rcut_chi(upf%nwfc), &
upf%rcutus_chi(upf%nwfc), &
upf%epseu(upf%nwfc) )
upf%rcut_chi(:) = 0.0_dp
upf%rcutus_chi(:)= 0.0_dp
upf%oc(upf%nwfc) = 0.0_dp
!IF ( upf%has_so ) THEN
! allocate ( upf%nn(upf%nwfc) )
! allocate ( upf%jchi(upf%nwfc) )
!END IF
do n=1,upf%nwfc
call xmlr_opentag ( 'pswf', IERR = ierr )
if ( ierr /= 0 ) return
call get_attr ( 'l', spdf )
call get_attr ( 'n', upf%nchi(n) )
call get_attr ( 'energy_level', upf%epseu(n) )
if ( spdf == 's' .or. spdf == 'S' ) then
upf%lchi(n) = 0
else if ( spdf == 'p' .or. spdf == 'P' ) then
upf%lchi(n) = 1
else if ( spdf == 'd' .or. spdf == 'D' ) then
upf%lchi(n) = 2
else if ( spdf == 'f' .or. spdf == 'F' ) then
upf%lchi(n) = 3
end if
write(upf%els(n),'(i1,a1)') upf%nchi(n), spdf
call xmlr_opentag ( 'radfunc' )
call xmlr_opentag ( 'data', IERR = ierr )
call get_attr ( 'npts', npt )
if ( npt > upf%mesh ) then
ierr = 1
return
end if
read (iun,*) upf%chi(1:npt,n)
if ( npt < upf%mesh ) upf%chi(npt+1:upf%mesh,n) = 0.0_dp
call xmlr_closetag () ! data
call xmlr_closetag () ! radfun
call xmlr_closetag () ! pwscf
end do
call xmlr_closetag ( ) ! nonlocal-projectors
!
END SUBROUTINE read_psml_pseudo_wave_functions
!
END subroutine read_psml
function libxc_to_qe (nxc, xc)
integer :: nxc
integer :: xc(nxc)
character(len=25) :: libxc_to_qe
!
if ( nxc < 1 ) go to 10
if ( xc(1) == 1 ) then
libxc_to_qe = 'SLA'
else
go to 10
end if
if ( nxc < 2 ) return
if ( xc(2) == 9 ) then
libxc_to_qe = trim(libxc_to_qe) // '-PZ'
else if ( xc(2) == 12 ) then
libxc_to_qe = trim(libxc_to_qe) // '-PW'
else
go to 10
end if
if ( nxc < 3 ) return
if ( xc(3) == 101 ) then
libxc_to_qe = trim(libxc_to_qe) // '-PBX'
else if ( xc(3) == 106 ) then
libxc_to_qe = trim(libxc_to_qe) // '-B88'
else if ( xc(3) == 116 ) then
libxc_to_qe = trim(libxc_to_qe) // '-PSX'
else
go to 10
end if
if ( nxc < 4 ) return
if ( xc(4) == 130 ) then
libxc_to_qe = trim(libxc_to_qe) // '-PBC'
else if ( xc(4) == 131 ) then
libxc_to_qe = trim(libxc_to_qe) // '-BLYP'
else if ( xc(4) == 133 ) then
libxc_to_qe = trim(libxc_to_qe) // '-PSC'
else
go to 10
end if
return
10 libxc_to_qe = 'Not Recognized'
return
!
end function libxc_to_qe

View File

@ -26,26 +26,30 @@ FUNCTION upf_get_pp_format (psfile) result(pp_format)
IMPLICIT NONE
INTEGER :: pp_format
CHARACTER (LEN=*) :: psfile
INTEGER :: l
INTEGER :: l, lm3, lm4, lm5
!
l = LEN_TRIM (psfile)
pp_format = 5
IF (l > 3) THEN
IF (psfile (l-3:l) =='.xml' .OR. psfile (l-3:l) =='.XML') THEN
pp_format = 0
ELSE IF (psfile (l-3:l) =='.upf' .OR. psfile (l-3:l) =='.UPF') THEN
pp_format = 1
ELSE IF (psfile (l-3:l) =='.vdb' .OR. psfile (l-3:l) =='.van') THEN
pp_format = 2
ELSE IF (psfile (l-3:l) =='.gth') THEN
pp_format = 3
ELSE IF (psfile (l-3:l) =='.cpi' .OR. psfile (l-3:l) =='.fhi') THEN
pp_format = 6
ELSE IF (psfile (l-4:l) =='.cpmd') THEN
pp_format = 7
ELSE IF (l > 5) THEN
If (psfile (l-5:l) =='.RRKJ3') pp_format = 4
END IF
lm3 = max(l-3,1)
lm4 = max(l-4,1)
lm5 = max(l-5,1)
IF (psfile (lm3:l) =='.xml' .OR. psfile (lm3:l) =='.XML') THEN
pp_format = 0
ELSE IF (psfile (lm3:l) =='.upf' .OR. psfile (lm3:l) =='.UPF') THEN
pp_format = 1
ELSE IF (psfile (lm3:l) =='.vdb' .OR. psfile (lm3:l) =='.van') THEN
pp_format = 2
ELSE IF (psfile (lm3:l) =='.gth') THEN
pp_format = 3
ELSE IF (psfile (lm3:l) =='.cpi' .OR. psfile (l-3:l) =='.fhi') THEN
pp_format = 6
ELSE IF (psfile (lm4:l) =='.cpmd') THEN
pp_format = 7
ELSE IF (psfile (lm4:l) =='.psml' ) THEN
pp_format = 8
ELSE IF (psfile (lm5:l) =='.RRKJ3') THEN
pp_format = 4
ELSE
pp_format = 5
END IF
!
END FUNCTION upf_get_pp_format

View File

@ -13,11 +13,12 @@ PROGRAM upfconv
! Pseudopotential conversion utility, can
! - convert from:
! UPF v.1 or v.2 containing "&" characters,
! old PWSCF norm-conserving and Ultrasoft formats
! Vanderbilt ultrasoft PP generation code format
! norm-conserving PPs in .psml format (experimental),
! old PWSCF norm-conserving and Ultrasoft formats,
! Vanderbilt ultrasoft PP generation code format,
! CPMD format (TYPE=NUMERIC, LOGARITHMIC, CAR, GOEDECKER)
! to:
! UPF v.2 clean, or UPF with schema (xml)
! UPF v.2, or new UPF with schema (xml)
! - extract and write to separate files:
! wavefunctions
! projectors ("beta" functions)
@ -74,6 +75,8 @@ PROGRAM upfconv
prefix_len = INDEX(TRIM(filein),'.UPF') - 1
ELSE IF (INDEX(TRIM(filein),'.upf') > 0 ) THEN
prefix_len = INDEX(TRIM(filein),'.upf') - 1
ELSE IF (INDEX(TRIM(filein),'.psml') > 0 ) THEN
prefix_len = INDEX(TRIM(filein),'.psml') - 1
ELSE
prefix_len = LEN_TRIM(filein)
ENDIF
@ -228,21 +231,23 @@ SUBROUTINE conv_upf2xml( upf )
IF ( version_compare(upf%nv,"2.0.1") == 'equal') RETURN
upf%nv="2.0.1"
!
IF ( .NOT. ALLOCATED(upf%nchi) ) THEN
ALLOCATE(upf%nchi(upf%nwfc))
upf%nchi(:) = 0
END IF
IF ( .NOT. ALLOCATED(upf%rcut_chi) ) THEN
ALLOCATE(upf%rcut_chi(upf%nwfc))
upf%rcut_chi(:) = upf%rcut(:)
END IF
IF ( .NOT. ALLOCATED(upf%rcutus_chi) ) THEN
ALLOCATE(upf%rcutus_chi(upf%nwfc))
upf%rcutus_chi(:) = upf%rcutus(:)
END IF
IF ( .NOT. ALLOCATED(upf%epseu) ) THEN
ALLOCATE(upf%epseu(upf%nwfc))
upf%epseu(:) = 0.0
IF ( upf%nwfc > 0 ) THEN
IF ( .NOT. ALLOCATED(upf%nchi) ) THEN
ALLOCATE(upf%nchi(upf%nwfc))
upf%nchi(:) = 0
END IF
IF ( .NOT. ALLOCATED(upf%rcut_chi) ) THEN
ALLOCATE(upf%rcut_chi(upf%nwfc))
upf%rcut_chi(:) = upf%rcut(:)
END IF
IF ( .NOT. ALLOCATED(upf%rcutus_chi) ) THEN
ALLOCATE(upf%rcutus_chi(upf%nwfc))
upf%rcutus_chi(:) = upf%rcutus(:)
END IF
IF ( .NOT. ALLOCATED(upf%epseu) ) THEN
ALLOCATE(upf%epseu(upf%nwfc))
upf%epseu(:) = 0.0
END IF
END IF
IF ( TRIM(upf%rel) == '' ) THEN
IF (upf%has_so) THEN