scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Source "normalization" using Norbert's script dev-tools/src-normal. 

Changes should not affect functionalities (please verify!)


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6825 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2010-06-11 09:00:04 +00:00
parent 687bd89938
commit 13e61c2110
15 changed files with 2378 additions and 2378 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
upftools/any2upf.f90
View File

@ -7,72 +7,72 @@
!
!
!---------------------------------------------------------------------
program mypp2upf
PROGRAM mypp2upf
!---------------------------------------------------------------------
!
! Convert a pseudopotential written in a user-supplied format
! to unified pseudopotential format - sample program
!
implicit none
character(len=256) filein, fileout
IMPLICIT NONE
CHARACTER(len=256) filein, fileout
!
!
call get_file ( filein )
open (unit = 1, file = filein, status = 'old', form = 'formatted')
call read_mypp(1)
close (1)
CALL get_file ( filein )
OPEN (unit = 1, file = filein, status = 'old', form = 'formatted')
CALL read_mypp(1)
CLOSE (1)
! convert variables read from user-supplied format into those needed
! by the upf format - add missing quantities
call convert_mypp
CALL convert_mypp
fileout=trim(filein)//'.UPF'
print '(''Output PP file in UPF format : '',a)', fileout
PRINT '(''Output PP file in UPF format : '',a)', fileout
open(unit=2,file=fileout,status='unknown',form='formatted')
call write_upf(2)
close (unit=2)
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
CALL write_upf(2)
CLOSE (unit=2)
stop
20 write (6,'("mypp2upf: error reading pseudopotential file name")')
stop
end program mypp2upf
STOP
20 WRITE (6,'("mypp2upf: error reading pseudopotential file name")')
STOP
END PROGRAM mypp2upf
module mypp
MODULE mypp
!
! All variables read from user-supplied file format
! Must have distinct names from variables in the "upf" module
!
end module mypp
END MODULE mypp
!
! ----------------------------------------------------------
subroutine read_mypp(iunps)
SUBROUTINE read_mypp(iunps)
! ----------------------------------------------------------
!
use mypp
implicit none
integer :: iunps
USE mypp
IMPLICIT NONE
INTEGER :: iunps
!
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential successfully read'
WRITE (6,'(a)') 'Pseudopotential successfully read'
! ----------------------------------------------------------
!
return
100 write (6,'("read_mypp: error reading pseudopotential file")')
stop
end subroutine read_mypp
RETURN
100 WRITE (6,'("read_mypp: error reading pseudopotential file")')
STOP
END SUBROUTINE read_mypp
! ----------------------------------------------------------
subroutine convert_mypp
SUBROUTINE convert_mypp
! ----------------------------------------------------------
!
use mypp
use upf
implicit none
USE mypp
USE upf
IMPLICIT NONE
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential successfully converted'
WRITE (6,'(a)') 'Pseudopotential successfully converted'
! ----------------------------------------------------------
return
end subroutine convert_mypp
RETURN
END SUBROUTINE convert_mypp

124
upftools/casino2upf.f90
View File

@ -44,7 +44,7 @@ PROGRAM casino2upf
CALL convert_casino
fileout=TRIM(filein)//'.UPF'
fileout=trim(filein)//'.UPF'
PRINT '(''Output PP file in US format : '',a)', fileout
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
@ -112,7 +112,7 @@ SUBROUTINE read_casino(iunps,nofiles)
LOGICAL :: groundstate, found
CHARACTER(len=1), DIMENSION(0:3) :: convel=(/'s','p','d','f'/)
CHARACTER(len=2) :: label, rellab
REAL(DP), parameter :: r_exp=20._dp/1500._dp
REAL(DP), PARAMETER :: r_exp=20._dp/1500._dp
INTEGER :: l, i, ir, nb, gsorbs, j,k,m,tmp, lquant, orbs, nquant
INTEGER, ALLOCATABLE :: gs(:,:)
@ -122,14 +122,14 @@ SUBROUTINE read_casino(iunps,nofiles)
nlc = 0 !These two values are always 0 for numeric pps
nnl = 0 !so lets just hard code them
nlcc_ = .FALSE. !Again these two are alwas false for CASINO pps
bhstype = .FALSE.
nlcc_ = .false. !Again these two are alwas false for CASINO pps
bhstype = .false.
READ(iunps,'(a2,35x,a2)') rellab, psd_
READ(iunps,*)
IF ( rellab .EQ. 'DF' ) THEN
IF ( rellab == 'DF' ) THEN
rel_=1
ELSE
rel_=0
@ -140,9 +140,9 @@ SUBROUTINE read_casino(iunps,nofiles)
READ(iunps,*)
ENDDO
READ(iunps,*) lmax_ !reading in lmax
IF ( zp_.LE.0d0 ) &
IF ( zp_<=0d0 ) &
CALL errore( 'read_casino','Wrong zp ',1 )
IF ( lmax_.GT.3.OR.lmax_.LT.0 ) &
IF ( lmax_>3.or.lmax_<0 ) &
CALL errore( 'read_casino','Wrong lmax ',1 )
lloc=lmax_ !think lloc shoudl always = lmax for this case yes/no ??
@ -165,7 +165,7 @@ SUBROUTINE read_casino(iunps,nofiles)
ENDDO
DO ir = 1, mesh_
rab_(ir) = r_exp * r_(ir) !hardcoded at the moment
END DO
ENDDO
ALLOCATE(vnl(mesh_,0:lmax_))
@ -192,7 +192,7 @@ SUBROUTINE read_casino(iunps,nofiles)
!
DO l = 0, lmax_
IF ( l.NE.lloc ) vnl(:,l) = vnl(:,l) - vnl(:,lloc)
IF ( l/=lloc ) vnl(:,l) = vnl(:,l) - vnl(:,lloc)
ENDDO
@ -202,7 +202,7 @@ SUBROUTINE read_casino(iunps,nofiles)
mptr => mhead
NULLIFY(mtail%p)
groundstate=.TRUE.
groundstate=.true.
DO j=1,nofiles
DO i=1,4
@ -218,7 +218,7 @@ SUBROUTINE read_casino(iunps,nofiles)
gs = 0
gsorbs = orbs
END IF
ENDIF
DO i=1,2
READ(j,*)
@ -231,20 +231,20 @@ SUBROUTINE read_casino(iunps,nofiles)
READ(j,*) tmp, nquant, lquant
IF ( groundstate ) THEN
found = .TRUE.
found = .true.
DO m=1,orbs
IF ( (nquant==gs(m,1) .AND. lquant==gs(m,2)) ) THEN
IF ( (nquant==gs(m,1) .and. lquant==gs(m,2)) ) THEN
gs(m,3) = gs(m,3) + 1
EXIT
END IF
exit
ENDIF
found = .FALSE.
found = .false.
ENDDO
IF (.NOT. found ) THEN
IF (.not. found ) THEN
DO m=1,orbs
@ -253,7 +253,7 @@ SUBROUTINE read_casino(iunps,nofiles)
gs(m,2) = lquant
gs(m,3) = 1
EXIT
exit
ENDIF
ENDDO
@ -275,34 +275,34 @@ SUBROUTINE read_casino(iunps,nofiles)
READ(j,'(13x,a2)', err=300) label
READ(j,*) tmp, nquant, lquant
IF ( .NOT. groundstate ) THEN
found = .FALSE.
IF ( .not. groundstate ) THEN
found = .false.
DO m = 1,gsorbs
IF ( nquant == gs(m,1) .AND. lquant == gs(m,2) ) THEN
found = .TRUE.
EXIT
END IF
END DO
IF ( nquant == gs(m,1) .and. lquant == gs(m,2) ) THEN
found = .true.
exit
ENDIF
ENDDO
mptr => mhead
DO
IF ( .NOT. ASSOCIATED(mptr) )EXIT
IF ( nquant == mptr%nquant .AND. lquant == mptr%lquant ) found = .TRUE.
IF ( .not. associated(mptr) )exit
IF ( nquant == mptr%nquant .and. lquant == mptr%lquant ) found = .true.
mptr =>mptr%p
END DO
ENDDO
IF ( found ) THEN
DO i=1,mesh_
READ(j,*)
ENDDO
CYCLE
END IF
END IF
ENDIF
ENDIF
#ifdef __STD_F95
IF ( ASSOCIATED(mtail%wavefunc) ) THEN
IF ( associated(mtail%wavefunc) ) THEN
#else
IF ( ALLOCATED(mtail%wavefunc) ) THEN
IF ( allocated(mtail%wavefunc) ) THEN
#endif
ALLOCATE(mtail%p)
mtail=>mtail%p
@ -310,7 +310,7 @@ SUBROUTINE read_casino(iunps,nofiles)
ALLOCATE( mtail%wavefunc(mesh_) )
ELSE
ALLOCATE( mtail%wavefunc(mesh_) )
END IF
ENDIF
mtail%label = label
mtail%nquant = nquant
mtail%lquant = lquant
@ -318,17 +318,17 @@ SUBROUTINE read_casino(iunps,nofiles)
READ(j, *, err=300) (mtail%wavefunc(ir),ir=1,mesh_)
ENDDO
groundstate = .FALSE.
groundstate = .false.
ENDDO
nchi =0
mptr => mhead
DO
IF ( .NOT. ASSOCIATED(mptr) )EXIT
IF ( .not. associated(mptr) )exit
nchi=nchi+1
mptr =>mptr%p
END DO
ENDDO
ALLOCATE(lchi_(nchi), els(nchi), nns_(nchi))
ALLOCATE(oc_(nchi))
@ -339,12 +339,12 @@ SUBROUTINE read_casino(iunps,nofiles)
DO i=1,gsorbs
oc_(i)=gs(i,3)
ENDDO
deallocate( gs )
DEALLOCATE( gs )
i=1
mptr => mhead
DO
IF ( .NOT. ASSOCIATED(mptr) )EXIT
IF ( .not. associated(mptr) )exit
nns_(i) = mptr%nquant
lchi_(i) = mptr%lquant
els(i) = mptr%label
@ -353,18 +353,18 @@ SUBROUTINE read_casino(iunps,nofiles)
chi_(ir:,i) = mptr%wavefunc(ir)
ENDDO
deallocate( mptr%wavefunc )
DEALLOCATE( mptr%wavefunc )
mptr =>mptr%p
i=i+1
END DO
ENDDO
!Clean up the linked list (deallocate it)
DO
IF ( .NOT. ASSOCIATED(mhead) )EXIT
IF ( .not. associated(mhead) )exit
mptr => mhead
mhead => mhead%p
deallocate( mptr )
END DO
DEALLOCATE( mptr )
ENDDO
!
@ -373,9 +373,9 @@ SUBROUTINE read_casino(iunps,nofiles)
ALLOCATE(rho_at_(mesh_))
rho_at_(:)=0.d0
DO nb = 1, nchi
IF( oc_(nb).NE.0.d0) &
IF( oc_(nb)/=0.d0) &
& rho_at_(:) = rho_at_(:) + oc_(nb)*chi_(:,nb)**2
END DO
ENDDO
! ----------------------------------------------------------
WRITE (6,'(a)') 'Pseudopotential successfully read'
! ----------------------------------------------------------
@ -416,7 +416,7 @@ SUBROUTINE convert_casino
rcutus(i)= 0.0d0
oc (i) = oc_(i)
epseu(i) = 0.0d0
END DO
ENDDO
DEALLOCATE (lchi_, oc_, nns_)
psd = psd_
@ -430,7 +430,7 @@ SUBROUTINE convert_casino
lmax = lmax_-1
ELSE
lmax = lmax_
END IF
ENDIF
nbeta= lmax_
mesh = mesh_
ntwfc= nchi
@ -439,7 +439,7 @@ SUBROUTINE convert_casino
lchiw(i) = lchi(i)
ocw(i) = oc(i)
elsw(i) = els(i)
END DO
ENDDO
CALL set_dft_from_name(dft_)
iexch = get_iexch()
icorr = get_icorr()
@ -465,12 +465,12 @@ SUBROUTINE convert_casino
DO ir = 1,mesh
IF ( r(ir) > rmax ) THEN
kkbeta=ir
EXIT
END IF
END DO
exit
ENDIF
ENDDO
! make sure kkbeta is odd as required for simpson
IF(MOD(kkbeta,2) == 0) kkbeta=kkbeta-1
IF(mod(kkbeta,2) == 0) kkbeta=kkbeta-1
ikk2(:) = kkbeta
ALLOCATE(aux(kkbeta))
ALLOCATE(betar(mesh,nbeta))
@ -483,18 +483,18 @@ SUBROUTINE convert_casino
iv=0
DO i=1,nchi
l=lchi(i)
IF (l.NE.lloc) THEN
IF (l/=lloc) THEN
iv=iv+1
lll(iv)=l
DO ir=1,kkbeta
betar(ir,iv)=chi_(ir,i)*vnl(ir,l)
aux(ir) = chi_(ir,i)**2*vnl(ir,l)
END DO
ENDDO
CALL simpson(kkbeta,aux,rab,vll)
dion(iv,iv) = 1.0d0/vll
END IF
IF(iv >= nbeta) EXIT ! skip additional pseudo wfns
ENDIF
IF(iv >= nbeta) exit ! skip additional pseudo wfns
ENDDO
@ -506,13 +506,13 @@ SUBROUTINE convert_casino
DO iv=1,nbeta
ikk2(iv)=kkbeta
DO ir = kkbeta,1,-1
IF ( ABS(betar(ir,iv)) > 1.d-12 ) THEN
IF ( abs(betar(ir,iv)) > 1.d-12 ) THEN
ikk2(iv)=ir
EXIT
END IF
END DO
END DO
END IF
exit
ENDIF
ENDDO
ENDDO
ENDIF
ALLOCATE (rho_at(mesh))
rho_at = rho_at_
DEALLOCATE (rho_at_)

418
upftools/cpmd2upf.f90
View File

@ -7,236 +7,236 @@
!
!
!---------------------------------------------------------------------
program cpmd2upf
PROGRAM cpmd2upf
!---------------------------------------------------------------------
!
! Convert a pseudopotential written in the CPMD format
! (TYPE=NORMCONSERVING NUMERIC only, single radial grid)
! to unified pseudopotential format
!
implicit none
character(len=256) filein, fileout
IMPLICIT NONE
CHARACTER(len=256) filein, fileout
!
!
call get_file ( filein )
open (unit = 1, file = filein, status = 'old', form = 'formatted')
call read_cpmd(1)
close (1)
CALL get_file ( filein )
OPEN (unit = 1, file = filein, status = 'old', form = 'formatted')
CALL read_cpmd(1)
CLOSE (1)
! convert variables read from CPMD format into those needed
! by the upf format - add missing quantities
call convert_cpmd
CALL convert_cpmd
fileout=trim(filein)//'.UPF'
print '(''Output PP file in UPF format : '',a)', fileout
PRINT '(''Output PP file in UPF format : '',a)', fileout
open(unit=2,file=fileout,status='unknown',form='formatted')
call write_upf(2)
close (unit=2)
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
CALL write_upf(2)
CLOSE (unit=2)
stop
20 call errore ('cpmd2upf', 'Reading pseudo file name ', 1)
STOP
20 CALL errore ('cpmd2upf', 'Reading pseudo file name ', 1)
end program cpmd2upf
END PROGRAM cpmd2upf
module cpmd
MODULE cpmd
!
! All variables read from CPMD file format
!
character (len=80) title
CHARACTER (len=80) title
!
integer :: ixc
INTEGER :: ixc
real(8) :: alphaxc
integer :: z, zv
INTEGER :: z, zv
!
integer :: mesh_
INTEGER :: mesh_
real(8) :: amesh, amesh_
real(8), allocatable :: r_(:)
real(8), ALLOCATABLE :: r_(:)
!
integer ::lmax_
real(8), allocatable :: vnl(:,:)
real(8), allocatable :: chi_(:,:)
INTEGER ::lmax_
real(8), ALLOCATABLE :: vnl(:,:)
real(8), ALLOCATABLE :: chi_(:,:)
!
logical :: nlcc_
real(8), allocatable :: rho_atc_(:)
LOGICAL :: nlcc_
real(8), ALLOCATABLE :: rho_atc_(:)
!
integer :: maxinfo_, info_lines_
parameter (maxinfo_ = 100)
character (len=80), allocatable :: info_sect_(:)
INTEGER :: maxinfo_, info_lines_
PARAMETER (maxinfo_ = 100)
CHARACTER (len=80), ALLOCATABLE :: info_sect_(:)
!------------------------------
end module cpmd
END MODULE cpmd
!
! ----------------------------------------------------------
subroutine read_cpmd(iunps)
SUBROUTINE read_cpmd(iunps)
! ----------------------------------------------------------
!
use cpmd
implicit none
integer :: iunps
USE cpmd
IMPLICIT NONE
INTEGER :: iunps
!
integer :: found = 0, closed = 0, unknown = 0
integer :: i, l, ios
character (len=80) line
character (len=4) token
INTEGER :: found = 0, closed = 0, unknown = 0
INTEGER :: i, l, ios
CHARACTER (len=80) line
CHARACTER (len=4) token
real (8) :: vnl0(0:3)
logical, external :: matches
integer, external :: locate
LOGICAL, EXTERNAL :: matches
INTEGER, EXTERNAL :: locate
!
nlcc_ = .false.
info_lines_ = 0
10 read (iunps,'(A)',end=20,err=20) line
if (matches ("&ATOM", trim(line)) ) then
10 READ (iunps,'(A)',end=20,err=20) line
IF (matches ("&ATOM", trim(line)) ) THEN
found = found + 1
! Z
read (iunps,'(a)',end=200,err=200) line
READ (iunps,'(a)',end=200,err=200) line
l = len_trim(line)
i = locate('=',line)
read (line(i+1:l),*) z
READ (line(i+1:l),*) z
! ZV
read (iunps,'(a)',end=200,err=200) line
READ (iunps,'(a)',end=200,err=200) line
l = len_trim(line)
i = locate('=',line)
read (line(i+1:l),*) zv
READ (line(i+1:l),*) zv
! XC
read (iunps,'(a)',end=200,err=200) line
READ (iunps,'(a)',end=200,err=200) line
l = len_trim(line)
i = locate('=',line)
read (line(i+1:l),*) ixc, alphaxc
READ (line(i+1:l),*) ixc, alphaxc
! TYPE
read (iunps,'(a)',end=200,err=200) line
if (.not. matches("NORMCONSERVING",line) .or. &
READ (iunps,'(a)',end=200,err=200) line
IF (.not. matches("NORMCONSERVING",line) .or. &
.not. matches("NUMERIC",line) ) &
call errore('read_cpmd','unknown type: '//line,1)
else if (matches ("&INFO", trim(line)) ) then
CALL errore('read_cpmd','unknown type: '//line,1)
ELSEIF (matches ("&INFO", trim(line)) ) THEN
found = found + 1
! read (iunps,'(a)') title
! store info section for later perusal (FIXME: not yet implemented. 2004/10/12, AK)
allocate (info_sect_(maxinfo_))
do i=1,maxinfo_
read (iunps,'(a)',end=20,err=20) title
if (matches ("&END", trim(title)) ) then
ALLOCATE (info_sect_(maxinfo_))
DO i=1,maxinfo_
READ (iunps,'(a)',end=20,err=20) title
IF (matches ("&END", trim(title)) ) THEN
closed = closed + 1
goto 10
else
GOTO 10
ELSE
info_sect_(i) = trim(title)
info_lines_ = i
end if
enddo
else if (matches ("&POTENTIAL", trim(line)) ) then
ENDIF
ENDDO
ELSEIF (matches ("&POTENTIAL", trim(line)) ) THEN
found = found + 1
read (iunps,'(a)') line
read (line,*,iostat=ios) mesh_, amesh
if ( ios /= 0) then
read (line,*,iostat=ios) mesh_
READ (iunps,'(a)') line
READ (line,*,iostat=ios) mesh_, amesh
IF ( ios /= 0) THEN
READ (line,*,iostat=ios) mesh_
amesh = -1.0d0
end if
allocate (r_(mesh_))
ENDIF
ALLOCATE (r_(mesh_))
!
! determine the number of angular momenta
!
read (iunps, '(a)') line
READ (iunps, '(a)') line
ios = 1
lmax_=4
do while (ios /= 0)
DO WHILE (ios /= 0)
lmax_ = lmax_ - 1
read(line,*,iostat=ios) r_(1),(vnl0(l),l=0,lmax_)
end do
allocate (vnl(mesh_,0:lmax_))
READ(line,*,iostat=ios) r_(1),(vnl0(l),l=0,lmax_)
ENDDO
ALLOCATE (vnl(mesh_,0:lmax_))
vnl(1,0:lmax_) = vnl0(0:lmax_)
do i=2,mesh_
read(iunps, *) r_(i),(vnl(i,l),l=0,lmax_)
end do
DO i=2,mesh_
READ(iunps, *) r_(i),(vnl(i,l),l=0,lmax_)
ENDDO
! get amesh if not available directly, check its value otherwise
print "('Radial grid r(i) has ',i4,' points')", mesh_
print "('Assuming log radial grid: r(i)=exp[(i-1)*amesh]*r(1), with:')"
if (amesh < 0.0d0) then
PRINT "('Radial grid r(i) has ',i4,' points')", mesh_
PRINT "('Assuming log radial grid: r(i)=exp[(i-1)*amesh]*r(1), with:')"
IF (amesh < 0.0d0) THEN
amesh = log (r_(mesh_)/r_(1))/(mesh_-1)
print "('amesh = log (r(mesh)/r(1))/(mesh-1) = ',f10.6)",amesh
else
PRINT "('amesh = log (r(mesh)/r(1))/(mesh-1) = ',f10.6)",amesh
ELSE
! not clear whether the value of amesh read from file
! matches the above definition, or if it is exp(amesh) ...
amesh_ = log (r_(mesh_)/r_(1))/(mesh_-1)
if ( abs ( amesh - amesh_ ) > 1.0d-5 ) then
if ( abs ( amesh - exp(amesh_) ) < 1.0d-5 ) then
IF ( abs ( amesh - amesh_ ) > 1.0d-5 ) THEN
IF ( abs ( amesh - exp(amesh_) ) < 1.0d-5 ) THEN
amesh = log(amesh)
print "('amesh = log (value read from file) = ',f10.6)",amesh
else
call errore ('cpmd2upf', 'unknown real-space grid',2)
end if
else
print "('amesh = value read from file = ',f10.6)",amesh
end if
end if
PRINT "('amesh = log (value read from file) = ',f10.6)",amesh
ELSE
CALL errore ('cpmd2upf', 'unknown real-space grid',2)
ENDIF
ELSE
PRINT "('amesh = value read from file = ',f10.6)",amesh
ENDIF
ENDIF
! check if the grid is what we expect
do i=2,mesh_
if ( abs(r_(i) - exp((i-1)*amesh)*r_(1)) > 1.0d-5) then
print "('grid point ',i4,': found ',f10.6,', expected ',f10.6)",&
DO i=2,mesh_
IF ( abs(r_(i) - exp((i-1)*amesh)*r_(1)) > 1.0d-5) THEN
PRINT "('grid point ',i4,': found ',f10.6,', expected ',f10.6)",&
i, r_(i), exp((i-1)*amesh)*r_(1)
call errore ('cpmd2upf', 'unknown real-space grid',1)
end if
end do
else if (matches ("&WAVEFUNCTION", trim(line)) ) then
CALL errore ('cpmd2upf', 'unknown real-space grid',1)
ENDIF
ENDDO
ELSEIF (matches ("&WAVEFUNCTION", trim(line)) ) THEN
found = found + 1
! read (iunps,*) mesh_, amesh
read (iunps,'(a)') line
read (line,*,iostat=ios) mesh_
allocate(chi_(mesh_,lmax_+1))
do i=1,mesh_
read(iunps, *) r_(i),(chi_(i,l+1),l=0,lmax_)
end do
else if (matches ("&NLCC", trim(line)) ) then
READ (iunps,'(a)') line
READ (line,*,iostat=ios) mesh_
ALLOCATE(chi_(mesh_,lmax_+1))
DO i=1,mesh_
READ(iunps, *) r_(i),(chi_(i,l+1),l=0,lmax_)
ENDDO
ELSEIF (matches ("&NLCC", trim(line)) ) THEN
found = found + 1
nlcc_ = .true.
read (iunps, '(a)') line
if (.not. matches ("NUMERIC", trim(line)) ) &
call errore('read_cpmd',' only NUMERIC core-correction supported',1)
read(iunps, *) mesh_
allocate (rho_atc_(mesh_))
read(iunps, * ) (r_(i), rho_atc_(i), i=1,mesh_)
else if (matches ("&ATDENS", trim(line)) ) then
READ (iunps, '(a)') line
IF (.not. matches ("NUMERIC", trim(line)) ) &
CALL errore('read_cpmd',' only NUMERIC core-correction supported',1)
READ(iunps, *) mesh_
ALLOCATE (rho_atc_(mesh_))
READ(iunps, * ) (r_(i), rho_atc_(i), i=1,mesh_)
ELSEIF (matches ("&ATDENS", trim(line)) ) THEN
! skip over &ATDENS section, add others here, if there are more.
do while(.not. matches("&END", trim(line)))
read (iunps,'(a)') line
end do
else if (matches ("&END", trim(line)) ) then
DO WHILE(.not. matches("&END", trim(line)))
READ (iunps,'(a)') line
ENDDO
ELSEIF (matches ("&END", trim(line)) ) THEN
closed = closed + 1
else
print*, 'line ignored: ', line
ELSE
PRINT*, 'line ignored: ', line
unknown = unknown + 1
end if
go to 10
ENDIF
GOTO 10
20 continue
if (nlcc_ .and. found /= 5 .or. .not.nlcc_ .and. found /= 4) &
call errore('read_cpmd','some &FIELD card missing',found)
if (closed /= found) &
call errore('read_cpmd','some &END card missing',closed)
if (unknown /= 0 ) print '("WARNING: ",i3," cards not read")', unknown
20 CONTINUE
IF (nlcc_ .and. found /= 5 .or. .not.nlcc_ .and. found /= 4) &
CALL errore('read_cpmd','some &FIELD card missing',found)
IF (closed /= found) &
CALL errore('read_cpmd','some &END card missing',closed)
IF (unknown /= 0 ) PRINT '("WARNING: ",i3," cards not read")', unknown
return
200 call errore('read_cpmd','error in reading file',1)
RETURN
200 CALL errore('read_cpmd','error in reading file',1)
end subroutine read_cpmd
END SUBROUTINE read_cpmd
! ----------------------------------------------------------
subroutine convert_cpmd
SUBROUTINE convert_cpmd
! ----------------------------------------------------------
!
use cpmd
use upf
implicit none
real(8), parameter :: rmax = 10.0d0
real(8), allocatable :: aux(:)
USE cpmd
USE upf
IMPLICIT NONE
real(8), PARAMETER :: rmax = 10.0d0
real(8), ALLOCATABLE :: aux(:)
real(8) :: vll
character (len=20):: dft
character (len=2), external :: atom_name
integer :: lloc, kkbeta, my_lmax
integer :: l, i, ir, iv
CHARACTER (len=20):: dft
CHARACTER (len=2), EXTERNAL :: atom_name
INTEGER :: lloc, kkbeta, my_lmax
INTEGER :: l, i, ir, iv
!
write(generated, '("Generated using unknown code")')
write(date_author,'("Author: unknown Generation date: as well")')
WRITE(generated, '("Generated using unknown code")')
WRITE(date_author,'("Author: unknown Generation date: as well")')
comment = 'Info: automatically converted from CPMD format'
! NOTE: many CPMD pseudopotentials created with the 'Hamann' code
@ -246,34 +246,34 @@ subroutine convert_cpmd
! we need to be able to ignore that part or the resulting UPF file
! will be useless. so we first print the info section and ask
! for the LMAX to really use. AK 2005/03/30.
do i=1,info_lines_
print '(A)', info_sect_(i)
enddo
print '("lmax to use. (max.",I2,") > ",$)', lmax_
read (5,*) my_lmax
if ((my_lmax <= lmax_) .and. (my_lmax >= 0)) lmax_ = my_lmax
print '("l local (max.",I2,") > ",$)', lmax_
read (5,*) lloc
DO i=1,info_lines_
PRINT '(A)', info_sect_(i)
ENDDO
PRINT '("lmax to use. (max.",I2,") > ",$)', lmax_
READ (5,*) my_lmax
IF ((my_lmax <= lmax_) .and. (my_lmax >= 0)) lmax_ = my_lmax
PRINT '("l local (max.",I2,") > ",$)', lmax_
READ (5,*) lloc
! reasonable assumption
if (z > 18) then
IF (z > 18) THEN
rel = 1
else
ELSE
rel = 0
end if
ENDIF
rcloc = 0.0d0
nwfs = lmax_+1
allocate( els(nwfs), oc(nwfs), epseu(nwfs))
allocate(lchi(nwfs), nns(nwfs) )
allocate(rcut (nwfs), rcutus (nwfs))
do i=1, nwfs
print '("Wavefunction # ",i1,": label, occupancy > ",$)', i
read (5,*) els(i), oc(i)
ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs))
ALLOCATE(lchi(nwfs), nns(nwfs) )
ALLOCATE(rcut (nwfs), rcutus (nwfs))
DO i=1, nwfs
PRINT '("Wavefunction # ",i1,": label, occupancy > ",$)', i
READ (5,*) els(i), oc(i)
nns (i) = 0
lchi(i) = i-1
rcut(i) = 0.0d0
rcutus(i)= 0.0d0
epseu(i) = 0.0d0
end do
ENDDO
psd = atom_name (z)
pseudotype = 'NC'
nlcc = nlcc_
@ -281,20 +281,20 @@ subroutine convert_cpmd
etotps =0.0d0
ecutrho=0.0d0
ecutwfc=0.0d0
if ( lmax_ == lloc) then
IF ( lmax_ == lloc) THEN
lmax = lmax_-1
else
ELSE
lmax = lmax_
end if
ENDIF
nbeta= lmax_
mesh = mesh_
ntwfc= nwfs
allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
do i=1, nwfs
ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
DO i=1, nwfs
lchiw(i) = lchi(i)
ocw(i) = oc(i)
elsw(i) = els(i)
end do
ENDDO
iexch = ixc/1000
icorr = (ixc-1000*iexch)/100
igcx = (ixc-1000*iexch-100*icorr)/10
@ -302,91 +302,91 @@ subroutine convert_cpmd
!
! We have igcc=2 (PW91) and 3 (LYP) exchanged wrt CPMD conventions
!
if (igcc.eq.3) then
IF (igcc==3) THEN
igcc=2
else if (igcc.eq.2) then
ELSEIF (igcc==2) THEN
igcc=3
end if
ENDIF
allocate(rab(mesh))
allocate( r(mesh))
ALLOCATE(rab(mesh))
ALLOCATE( r(mesh))
r = r_
rab = r * amesh
allocate (rho_atc(mesh))
if (nlcc) rho_atc = rho_atc_
ALLOCATE (rho_atc(mesh))
IF (nlcc) rho_atc = rho_atc_
allocate (vloc0(mesh))
ALLOCATE (vloc0(mesh))
! the factor 2 converts from Hartree to Rydberg
vloc0(:) = vnl(:,lloc)*2.d0
if (nbeta > 0) then
IF (nbeta > 0) THEN
allocate(ikk2(nbeta), lll(nbeta))
ALLOCATE(ikk2(nbeta), lll(nbeta))
kkbeta=mesh
do ir = 1,mesh
if ( r(ir) > rmax ) then
DO ir = 1,mesh
IF ( r(ir) > rmax ) THEN
kkbeta=ir
exit
end if
end do
ENDIF
ENDDO
ikk2(:) = kkbeta
allocate(aux(kkbeta))
allocate(betar(mesh,nbeta))
allocate(qfunc(mesh,nbeta,nbeta))
allocate(dion(nbeta,nbeta))
allocate(qqq (nbeta,nbeta))
ALLOCATE(aux(kkbeta))
ALLOCATE(betar(mesh,nbeta))
ALLOCATE(qfunc(mesh,nbeta,nbeta))
ALLOCATE(dion(nbeta,nbeta))
ALLOCATE(qqq (nbeta,nbeta))
qfunc(:,:,:)=0.0d0
dion(:,:) =0.d0
qqq(:,:) =0.d0
iv=0
do i=1,nwfs
DO i=1,nwfs
l=lchi(i)
if (l.ne.lloc) then
IF (l/=lloc) THEN
iv=iv+1
lll(iv)=l
do ir=1,kkbeta
DO ir=1,kkbeta
! the factor 2 converts from Hartree to Rydberg
betar(ir,iv) = 2.d0 * chi_(ir,l+1) * &
( vnl(ir,l) - vnl(ir,lloc) )
aux(ir) = chi_(ir,l+1) * betar(ir,iv)
end do
call simpson(kkbeta,aux,rab,vll)
ENDDO
CALL simpson(kkbeta,aux,rab,vll)
dion(iv,iv) = 1.0d0/vll
end if
enddo
ENDIF
ENDDO
end if
ENDIF
allocate (rho_at(mesh))
ALLOCATE (rho_at(mesh))
rho_at = 0.d0
do i=1,nwfs
DO i=1,nwfs
rho_at(:) = rho_at(:) + ocw(i) * chi_(:,i) ** 2
end do
ENDDO
allocate (chi(mesh,ntwfc))
ALLOCATE (chi(mesh,ntwfc))
chi = chi_
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential successfully converted'
WRITE (6,'(a)') 'Pseudopotential successfully converted'
! ----------------------------------------------------------
return
end subroutine convert_cpmd
RETURN
END SUBROUTINE convert_cpmd
!
! ------------------------------------------------------------------
integer function locate(onechar,string)
INTEGER FUNCTION locate(onechar,string)
! ------------------------------------------------------------------
!
character(len=1) :: onechar
character(len=*) :: string
CHARACTER(len=1) :: onechar
CHARACTER(len=*) :: string
!
integer:: i
INTEGER:: i
!
do i=1,len_trim(string)
if (string(i:i) .eq. "=") then
DO i=1,len_trim(string)
IF (string(i:i) == "=") THEN
locate = i
return
end if
end do
RETURN
ENDIF
ENDDO
locate = 0
return
end function locate
RETURN
END FUNCTION locate

378
upftools/fhi2upf.f90
View File

@ -7,7 +7,7 @@
!
!
!---------------------------------------------------------------------
program fhi2upf
PROGRAM fhi2upf
!---------------------------------------------------------------------
!
! Convert a pseudopotential file in Fritz-Haber numerical format
@ -16,237 +16,237 @@ program fhi2upf
! May or may not work: carefully check what you get
! Adapted from the converter written by Andrea Ferretti
!
implicit none
character(len=256) filein, fileout
IMPLICIT NONE
CHARACTER(len=256) filein, fileout
!
!
call get_file ( filein )
open (unit = 1, file = filein, status = 'old', form = 'formatted')
call read_fhi(1)
close (1)
CALL get_file ( filein )
OPEN (unit = 1, file = filein, status = 'old', form = 'formatted')
CALL read_fhi(1)
CLOSE (1)
! convert variables read from FHI format into those needed
! by the upf format - add missing quantities
call convert_fhi
CALL convert_fhi
fileout=trim(filein)//'.UPF'
print '(''Output PP file in UPF format : '',a)', fileout
PRINT '(''Output PP file in UPF format : '',a)', fileout
open(unit=2,file=fileout,status='unknown',form='formatted')
call write_upf(2)
close (unit=2)
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
CALL write_upf(2)
CLOSE (unit=2)
stop
20 write (6,'("fhi2upf: error reading pseudopotential file name")')
stop
end program fhi2upf
STOP
20 WRITE (6,'("fhi2upf: error reading pseudopotential file name")')
STOP
END PROGRAM fhi2upf
module fhi
MODULE fhi
!
! All variables read from FHI file format
!
type angular_comp
real(8), pointer :: pot(:)
real(8), pointer :: wfc(:)
real(8), pointer :: grid(:)
TYPE angular_comp
real(8), POINTER :: pot(:)
real(8), POINTER :: wfc(:)
real(8), POINTER :: grid(:)
real(8) :: amesh
integer :: nmesh
integer :: lcomp
end type angular_comp
INTEGER :: nmesh
INTEGER :: lcomp
END TYPE angular_comp
!------------------------------
real(8) :: Zval ! valence charge
integer :: lmax_ ! max l-component used
INTEGER :: lmax_ ! max l-component used
logical :: nlcc_
real(8), allocatable :: rho_atc_(:) ! core charge
LOGICAL :: nlcc_
real(8), ALLOCATABLE :: rho_atc_(:) ! core charge
type (angular_comp), pointer :: comp(:) ! PP numerical info
TYPE (angular_comp), POINTER :: comp(:) ! PP numerical info
! (wfc, grid, potentials...)
!------------------------------
! variables for the abinit header
real(8) :: Zatom, Zion, r2well, rchrg, fchrg, qchrg
integer :: pspdat = 0, pspcod = 0 , pspxc = 0, lloc_ = -1, mmax = 0
character(len=256) :: info
INTEGER :: pspdat = 0, pspcod = 0 , pspxc = 0, lloc_ = -1, mmax = 0
CHARACTER(len=256) :: info
end module fhi
END MODULE fhi
!
! ----------------------------------------------------------
subroutine read_fhi(iunps)
SUBROUTINE read_fhi(iunps)
! ----------------------------------------------------------
!
use fhi
implicit none
integer, parameter :: Nl=7 ! max number of l-components
integer :: iunps
USE fhi
IMPLICIT NONE
INTEGER, PARAMETER :: Nl=7 ! max number of l-components
INTEGER :: iunps
real(8) :: r, rhoc, drhoc, d2rhoc
!
integer :: l, i, idum, mesh
INTEGER :: l, i, idum, mesh
! Start reading file
read(iunps,'(a)') info
read(info,*,iostat=i) Zval, l
if ( i /= 0 .or. zval <= 0.0 .or. zval > 100.0 ) then
write (6,'("read_fhi: assuming abinit format")')
read(iunps,*) Zatom, Zion, pspdat
read(iunps,*) pspcod, pspxc, lmax_,lloc_, mmax, r2well
if (pspcod /= 6) then
write (6,'("read_fhi: unknown PP type ",i1,"...stopping")') pspcod
stop
end if
read(iunps,*) rchrg, fchrg, qchrg
READ(iunps,'(a)') info
READ(info,*,iostat=i) Zval, l
IF ( i /= 0 .or. zval <= 0.0 .or. zval > 100.0 ) THEN
WRITE (6,'("read_fhi: assuming abinit format")')
READ(iunps,*) Zatom, Zion, pspdat
READ(iunps,*) pspcod, pspxc, lmax_,lloc_, mmax, r2well
IF (pspcod /= 6) THEN
WRITE (6,'("read_fhi: unknown PP type ",i1,"...stopping")') pspcod
STOP
ENDIF
READ(iunps,*) rchrg, fchrg, qchrg
!
read(iunps,*)
read(iunps,*)
read(iunps,*)
READ(iunps,*)
READ(iunps,*)
READ(iunps,*)
!
read(iunps,*) Zval, l
if (abs(Zion-Zval) > 1.0d-8) then
write (6,'("read_fhi: Zval/Zion mismatch...stopping")')
stop
end if
if (l-1 /= lmax_) then
write (6,'("read_fhi: lmax mismatch...stopping")')
stop
end if
else
READ(iunps,*) Zval, l
IF (abs(Zion-Zval) > 1.0d-8) THEN
WRITE (6,'("read_fhi: Zval/Zion mismatch...stopping")')
STOP
ENDIF
IF (l-1 /= lmax_) THEN
WRITE (6,'("read_fhi: lmax mismatch...stopping")')
STOP
ENDIF
ELSE
info = ' '
end if
ENDIF
lmax_ = l - 1
if (lmax_+1 > Nl) then
write (6,'("read_fhi: too many l-components...stopping")')
stop
end if
IF (lmax_+1 > Nl) THEN
WRITE (6,'("read_fhi: too many l-components...stopping")')
STOP
ENDIF
do i=1,10
read(iunps,*) ! skipping 11 lines
end do
DO i=1,10
READ(iunps,*) ! skipping 11 lines
ENDDO
allocate( comp(0:lmax_) )
ALLOCATE( comp(0:lmax_) )
do l=0,lmax_
DO l=0,lmax_
comp(l)%lcomp = l
read(iunps,*) comp(l)%nmesh, comp(l)%amesh
if (mmax > 0 .and. mmax /= comp(l)%nmesh) then
write (6,'("read_fhi: mismatched number of grid points...stopping")')
stop
end if
if ( l > 0) then
if (comp(l)%nmesh /= comp(0)%nmesh .or. &
comp(l)%amesh /= comp(0)%amesh ) then
write(6,'("read_fhi: different radial grids not allowed...stopping")')
stop
end if
end if
READ(iunps,*) comp(l)%nmesh, comp(l)%amesh
IF (mmax > 0 .and. mmax /= comp(l)%nmesh) THEN
WRITE (6,'("read_fhi: mismatched number of grid points...stopping")')
STOP
ENDIF
IF ( l > 0) THEN
IF (comp(l)%nmesh /= comp(0)%nmesh .or. &
comp(l)%amesh /= comp(0)%amesh ) THEN
WRITE(6,'("read_fhi: different radial grids not allowed...stopping")')
STOP
ENDIF
ENDIF
mesh = comp(l)%nmesh
allocate( comp(l)%wfc(mesh), & ! wave-functions
ALLOCATE( comp(l)%wfc(mesh), & ! wave-functions
comp(l)%pot(mesh), & ! potentials
comp(l)%grid(mesh) ) ! real space radial grid
! read the above quantities
do i=1,mesh
read(iunps,*) idum, comp(l)%grid(i), &
DO i=1,mesh
READ(iunps,*) idum, comp(l)%grid(i), &
comp(l)%wfc(i), &
comp(l)%pot(i)
end do
end do
ENDDO
ENDDO
nlcc_ =.false.
allocate(rho_atc_(comp(0)%nmesh))
ALLOCATE(rho_atc_(comp(0)%nmesh))
mesh = comp(0)%nmesh
do i=1,mesh
read(iunps,*,end=10, err=20) r, rho_atc_(i), drhoc, d2rhoc
if ( abs( r - comp(0)%grid(i) ) > 1.d-6 ) then
write(6,'("read_fhi: radial grid for core charge? stopping")')
stop
end if
end do
DO i=1,mesh
READ(iunps,*,end=10, err=20) r, rho_atc_(i), drhoc, d2rhoc
IF ( abs( r - comp(0)%grid(i) ) > 1.d-6 ) THEN
WRITE(6,'("read_fhi: radial grid for core charge? stopping")')
STOP
ENDIF
ENDDO
nlcc_ = .true.
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential with NLCC successfully read'
WRITE (6,'(a)') 'Pseudopotential with NLCC successfully read'
! ----------------------------------------------------------
return
10 continue
RETURN
10 CONTINUE
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential without NLCC successfully read'
WRITE (6,'(a)') 'Pseudopotential without NLCC successfully read'
! ----------------------------------------------------------
return
RETURN
!
20 write(6,'("read_fhi: error reading core charge")')
stop
20 WRITE(6,'("read_fhi: error reading core charge")')
STOP
!
100 write(6,'("read_fhi: error reading pseudopotential file")')
stop
100 WRITE(6,'("read_fhi: error reading pseudopotential file")')
STOP
end subroutine read_fhi
END SUBROUTINE read_fhi
! ----------------------------------------------------------
subroutine convert_fhi
SUBROUTINE convert_fhi
! ----------------------------------------------------------
!
use fhi
use upf
use funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc
use constants, ONLY : fpi
implicit none
real(8), parameter :: rmax = 10.0d0
real(8), allocatable :: aux(:)
USE fhi
USE upf
USE funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc
USE constants, ONLY : fpi
IMPLICIT NONE
real(8), PARAMETER :: rmax = 10.0d0
real(8), ALLOCATABLE :: aux(:)
real(8) :: vll
character (len=20):: dft
character (len=2), external:: atom_name
integer :: lloc, kkbeta
integer :: l, i, ir, iv
CHARACTER (len=20):: dft
CHARACTER (len=2), EXTERNAL:: atom_name
INTEGER :: lloc, kkbeta
INTEGER :: l, i, ir, iv
!
if (nint(Zatom) > 0) then
IF (nint(Zatom) > 0) THEN
psd = atom_name(nint(Zatom))
else
print '("Atom name > ",$)'
read (5,'(a)') psd
end if
if ( lloc_ < 0 ) then
print '("l local (max: ",i1,") > ",$)', lmax_
read (5,*) lloc
else
ELSE
PRINT '("Atom name > ",$)'
READ (5,'(a)') psd
ENDIF
IF ( lloc_ < 0 ) THEN
PRINT '("l local (max: ",i1,") > ",$)', lmax_
READ (5,*) lloc
ELSE
lloc = lloc_
end if
if (pspxc == 7) then
ENDIF
IF (pspxc == 7) THEN
dft = 'PW'
else
if (pspxc > 0) then
print '("DFT read from abinit file: ",i1)', pspxc
end if
print '("DFT > ",$)'
read (5,'(a)') dft
end if
write(generated, '("Generated using Fritz-Haber code")')
write(date_author,'("Author: unknown Generation date: as well")')
if (trim(info) /= ' ') then
ELSE
IF (pspxc > 0) THEN
PRINT '("DFT read from abinit file: ",i1)', pspxc
ENDIF
PRINT '("DFT > ",$)'
READ (5,'(a)') dft
ENDIF
WRITE(generated, '("Generated using Fritz-Haber code")')
WRITE(date_author,'("Author: unknown Generation date: as well")')
IF (trim(info) /= ' ') THEN
comment = trim(info)
else
ELSE
comment = 'Info: automatically converted from FHI format'
end if
ENDIF
! reasonable assumption
rel = 1
rcloc = 0.0d0
nwfs = lmax_+1
allocate( els(nwfs), oc(nwfs), epseu(nwfs))
allocate(lchi(nwfs), nns(nwfs) )
allocate(rcut (nwfs), rcutus (nwfs))
do i=1, nwfs
print '("Wavefunction # ",i1,": label, occupancy > ",$)', i
read (5,*) els(i), oc(i)
ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs))
ALLOCATE(lchi(nwfs), nns(nwfs) )
ALLOCATE(rcut (nwfs), rcutus (nwfs))
DO i=1, nwfs
PRINT '("Wavefunction # ",i1,": label, occupancy > ",$)', i
READ (5,*) els(i), oc(i)
nns (i) = 0
lchi(i) = i-1
rcut(i) = 0.0d0
rcutus(i)= 0.0d0
epseu(i) = 0.0d0
end do
ENDDO
pseudotype = 'NC'
nlcc = nlcc_
@ -254,91 +254,91 @@ subroutine convert_fhi
etotps = 0.0d0
ecutrho=0.0d0
ecutwfc=0.0d0
if ( lmax_ == lloc) then
IF ( lmax_ == lloc) THEN
lmax = lmax_-1
else
ELSE
lmax = lmax_
end if
ENDIF
nbeta= lmax_
mesh = comp(0)%nmesh
ntwfc= nwfs
allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
do i=1, nwfs
ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
DO i=1, nwfs
lchiw(i) = lchi(i)
ocw(i) = oc(i)
elsw(i) = els(i)
end do
call set_dft_from_name(dft)
ENDDO
CALL set_dft_from_name(dft)
iexch = get_iexch()
icorr = get_icorr()
igcx = get_igcx()
igcc = get_igcc()
allocate(rab(mesh))
allocate( r(mesh))
ALLOCATE(rab(mesh))
ALLOCATE( r(mesh))
r = comp(0)%grid
rab = r * log( comp(0)%amesh )
if (nlcc) then
allocate (rho_atc(mesh))
IF (nlcc) THEN
ALLOCATE (rho_atc(mesh))
rho_atc(:) = rho_atc_(:) / fpi
end if
ENDIF
allocate (vloc0(mesh))
ALLOCATE (vloc0(mesh))
! the factor 2 converts from Hartree to Rydberg
vloc0 = 2.d0*comp(lloc)%pot
if (nbeta > 0) then
IF (nbeta > 0) THEN
allocate(ikk2(nbeta), lll(nbeta))
ALLOCATE(ikk2(nbeta), lll(nbeta))
kkbeta=mesh
do ir = 1,mesh
if ( r(ir) > rmax ) then
DO ir = 1,mesh
IF ( r(ir) > rmax ) THEN
kkbeta=ir
exit
end if
end do
ENDIF
ENDDO
ikk2(:) = kkbeta
allocate(aux(kkbeta))
allocate(betar(mesh,nbeta))
allocate(qfunc(mesh,nbeta,nbeta))
allocate(dion(nbeta,nbeta))
allocate(qqq (nbeta,nbeta))
ALLOCATE(aux(kkbeta))
ALLOCATE(betar(mesh,nbeta))
ALLOCATE(qfunc(mesh,nbeta,nbeta))
ALLOCATE(dion(nbeta,nbeta))
ALLOCATE(qqq (nbeta,nbeta))
qfunc(:,:,:)=0.0d0
dion(:,:) =0.d0
qqq(:,:) =0.d0
iv=0
do i=1,nwfs
DO i=1,nwfs
l=lchi(i)
if (l.ne.lloc) then
IF (l/=lloc) THEN
iv=iv+1
lll(iv)=l
do ir=1,kkbeta
DO ir=1,kkbeta
! FHI potentials are in Hartree
betar(ir,iv) = 2.d0 * comp(l)%wfc(ir) * &
( comp(l)%pot(ir) - comp(lloc)%pot(ir) )
aux(ir) = comp(l)%wfc(ir) * betar(ir,iv)
end do
call simpson(kkbeta,aux,rab,vll)
ENDDO
CALL simpson(kkbeta,aux,rab,vll)
dion(iv,iv) = 1.0d0/vll
end if
enddo
ENDIF
ENDDO
end if
ENDIF
allocate (rho_at(mesh))
ALLOCATE (rho_at(mesh))
rho_at = 0.d0
do i=1,nwfs
DO i=1,nwfs
l=lchi(i)
rho_at = rho_at + ocw(i) * comp(l)%wfc ** 2
end do
ENDDO
allocate (chi(mesh,ntwfc))
do i=1,ntwfc
ALLOCATE (chi(mesh,ntwfc))
DO i=1,ntwfc
chi(:,i) = comp(i-1)%wfc(:)
end do
ENDDO
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential successfully converted'
WRITE (6,'(a)') 'Pseudopotential successfully converted'
! ----------------------------------------------------------
return
end subroutine convert_fhi
RETURN
END SUBROUTINE convert_fhi

494
upftools/fpmd2upf.f90
View File

@ -29,20 +29,20 @@
! Example:
module fpmd2upf_module
MODULE fpmd2upf_module
USE kinds, ONLY: DP
USE parameters
use radial_grids, ONLY: ndmx
USE radial_grids, ONLY: ndmx
implicit none
save
IMPLICIT NONE
SAVE
REAL(DP), PRIVATE :: TOLMESH = 1.d-5
TYPE pseudo_ncpp
CHARACTER(LEN=4) :: psd ! Element label
CHARACTER(LEN=20) :: pottyp ! Potential type
CHARACTER(len=4) :: psd ! Element label
CHARACTER(len=20) :: pottyp ! Potential type
LOGICAL :: tmix
LOGICAL :: tnlcc
INTEGER :: igau
@ -72,53 +72,53 @@ module fpmd2upf_module
END TYPE pseudo_ncpp
contains
CONTAINS
subroutine read_pseudo_fpmd( ap, psfile )
type(pseudo_ncpp) :: ap
character(len=256) :: psfile
character(len=80) :: error_msg
integer :: info, iunit
SUBROUTINE read_pseudo_fpmd( ap, psfile )
TYPE(pseudo_ncpp) :: ap
CHARACTER(len=256) :: psfile
CHARACTER(len=80) :: error_msg
INTEGER :: info, iunit
iunit = 11
OPEN(UNIT=iunit,FILE=psfile,STATUS='OLD')
REWIND( iunit )
CALL read_head_pp( iunit, ap, error_msg, info)
IF( info /= 0 ) GO TO 200
IF( info /= 0 ) GOTO 200
IF( ap%pottyp == 'GIANNOZ' ) THEN
CALL read_giannoz(iunit, ap, info)
IF( info /= 0 ) GO TO 200
IF( info /= 0 ) GOTO 200
ELSE IF( ap%pottyp == 'NUMERIC' ) THEN
ELSEIF( ap%pottyp == 'NUMERIC' ) THEN
CALL read_numeric_pp( iunit, ap, error_msg, info)
IF( info /= 0 ) GO TO 200
IF( info /= 0 ) GOTO 200
ELSE IF( ap%pottyp == 'ANALYTIC' ) THEN
ELSEIF( ap%pottyp == 'ANALYTIC' ) THEN
CALL read_analytic_pp( iunit, ap, error_msg, info)
IF( info /= 0 ) GO TO 200
IF( info /= 0 ) GOTO 200
ELSE
info = 1
error_msg = ' Pseudopotential type '//TRIM(ap%pottyp)//' not implemented '
GO TO 200
error_msg = ' Pseudopotential type '//trim(ap%pottyp)//' not implemented '
GOTO 200
END IF
ENDIF
200 CONTINUE
IF( info /= 0 ) THEN
CALL errore(' readpseudo ', error_msg, ABS(info) )
END IF
CALL errore(' readpseudo ', error_msg, abs(info) )
ENDIF
CLOSE(iunit)
return
end subroutine read_pseudo_fpmd
RETURN
END SUBROUTINE read_pseudo_fpmd
!=----------------------------------------------------------------------------=!
@ -127,28 +127,28 @@ contains
! ... This sub. check if a given fortran unit 'iunit' contains a UPF pseudopot.
INTEGER, INTENT(IN) :: iunit
INTEGER, INTENT(OUT) :: info
CHARACTER(LEN=80) :: dummy
INTEGER, INTENT(in) :: iunit
INTEGER, INTENT(out) :: info
CHARACTER(len=80) :: dummy
INTEGER :: ios
LOGICAL, EXTERNAL :: matches
info = 0
ios = 0
header_loop: do while (ios == 0)
read (iunit, *, iostat = ios, err = 200) dummy
if (matches ("<PP_HEADER>", dummy) ) then
header_loop: DO WHILE (ios == 0)
READ (iunit, *, iostat = ios, err = 200) dummy
IF (matches ("<PP_HEADER>", dummy) ) THEN
info = 1
exit header_loop
endif
enddo header_loop
200 continue
ENDIF
ENDDO header_loop
200 CONTINUE
RETURN
END SUBROUTINE check_file_type
!=----------------------------------------------------------------------------=!
SUBROUTINE analytic_to_numeric(ap)
TYPE (pseudo_ncpp), INTENT(INOUT) :: ap
TYPE (pseudo_ncpp), INTENT(inout) :: ap
INTEGER :: ir, mesh, lmax, l, n, il, ib, ll
REAL(DP) :: xmin, zmesh, dx, x
! REAL(DP) :: pi = 3.14159265358979323846_DP
@ -158,36 +158,36 @@ contains
IF( ap%mesh == 0 ) THEN
! ... Local pseudopotential, define a logaritmic grid
mesh = SIZE( ap%rw )
mesh = size( ap%rw )
xmin = -5.0d0
zmesh = 6.0d0
dx = 0.025d0
DO ir = 1, mesh
x = xmin + DBLE(ir-1) * dx
ap%rw(ir) = EXP(x) / zmesh
IF( ap%rw(ir) > 1000.0d0 ) EXIT
END DO
x = xmin + dble(ir-1) * dx
ap%rw(ir) = exp(x) / zmesh
IF( ap%rw(ir) > 1000.0d0 ) exit
ENDDO
ap%mesh = mesh
ap%dx = dx
ap%rab = ap%dx * ap%rw
END IF
ENDIF
ap%vnl = 0.0d0
ap%vloc = 0.0d0
ap%vrps = 0.0d0
do l = 1, 3
do ir = 1, ap%mesh
ap%vnl(ir,l)= - ( ap%wrc(1) * qe_erf(SQRT(ap%rc(1))*ap%rw(ir)) + &
ap%wrc(2) * qe_erf(SQRT(ap%rc(2))*ap%rw(ir)) ) *&
DO l = 1, 3
DO ir = 1, ap%mesh
ap%vnl(ir,l)= - ( ap%wrc(1) * qe_erf(sqrt(ap%rc(1))*ap%rw(ir)) + &
ap%wrc(2) * qe_erf(sqrt(ap%rc(2))*ap%rw(ir)) ) *&
ap%zv / ap%rw(ir)
end do
do ir = 1, ap%mesh
do n = 1, ap%igau
ENDDO
DO ir = 1, ap%mesh
DO n = 1, ap%igau
ap%vnl(ir,l)= ap%vnl(ir,l)+ (ap%al(n,l)+ ap%bl(n,l)*ap%rw(ir)**2 )* &
EXP(-ap%rcl(n,l)*ap%rw(ir)**2)
end do
end do
end do
exp(-ap%rcl(n,l)*ap%rw(ir)**2)
ENDDO
ENDDO
ENDDO
! ... Copy local component to a separate array
ap%vloc(:) = ap%vnl(:,ap%lloc)
@ -195,7 +195,7 @@ contains
ll=ap%lll(l) + 1 ! find out the angular momentum (ll-1) of the component stored
! in position l
ap%vrps(:,l) = ( ap%vnl(:,ll) - ap%vloc(:) ) * ap%rps(:,ll)
END DO
ENDDO
RETURN
END SUBROUTINE analytic_to_numeric
@ -205,24 +205,24 @@ contains
SUBROUTINE read_giannoz(uni, ap, ierr)
! USE constants, ONLY : fpi
IMPLICIT NONE
TYPE (pseudo_ncpp), INTENT(INOUT) :: ap
INTEGER, INTENT(IN) :: uni
INTEGER, INTENT(OUT) :: ierr
REAL(DP) :: chi( SIZE(ap%rps, 1), SIZE(ap%rps, 2) )
REAL(DP) :: vnl( SIZE(ap%vnl, 1), SIZE(ap%vnl, 2) )
REAL(DP) :: rho_core( SIZE(ap%rhoc, 1) )
TYPE (pseudo_ncpp), INTENT(inout) :: ap
INTEGER, INTENT(in) :: uni
INTEGER, INTENT(out) :: ierr
REAL(DP) :: chi( size(ap%rps, 1), size(ap%rps, 2) )
REAL(DP) :: vnl( size(ap%vnl, 1), size(ap%vnl, 2) )
REAL(DP) :: rho_core( size(ap%rhoc, 1) )
REAL(DP) :: r, ra, rb, fac
REAL(DP) :: oc( SIZE(ap%rps, 2) )
REAL(DP) :: enl( SIZE(ap%rps, 2) )
REAL(DP) :: oc( size(ap%rps, 2) )
REAL(DP) :: enl( size(ap%rps, 2) )
REAL(DP) :: zmesh, xmin, dx, etot
REAL(DP) :: zval
INTEGER :: nn(SIZE(ap%rps, 2)), ll(SIZE(ap%rps, 2))
INTEGER :: nn(size(ap%rps, 2)), ll(size(ap%rps, 2))
INTEGER :: nwf, mesh, i, j, in1, in2, in3, in4, m
INTEGER :: lmax, nlc, nnl, lloc, l, il
LOGICAL :: nlcc
CHARACTER(len=80) :: dft
CHARACTER(len=4) :: atom
CHARACTER(len=2) :: el( SIZE(ap%rps, 2) )
CHARACTER(len=2) :: el( size(ap%rps, 2) )
CHARACTER(len=80) :: ppinfo
CHARACTER(len=80) :: strdum
CHARACTER(len=2) :: sdum1, sdum2
@ -236,38 +236,38 @@ contains
! WRITE(6,*) ' DEBUG ', atom, zval,lmax, nlc, nnl, nlcc, lloc, ppinfo
IF( (lmax+1) > SIZE(ap%vnl, 2) ) THEN
IF( (lmax+1) > size(ap%vnl, 2) ) THEN
ierr = 1
RETURN
END IF
IF( (nlcc .AND. .NOT.ap%tnlcc) .OR. (.NOT.nlcc .AND. ap%tnlcc) ) THEN
ENDIF
IF( (nlcc .and. .not.ap%tnlcc) .or. (.not.nlcc .and. ap%tnlcc) ) THEN
ierr = 2
RETURN
END IF
ENDIF
READ(uni,fmt='(f8.2,f8.4,f10.6,2i6)') zmesh, xmin, dx, mesh, nwf
IF( mesh > SIZE(ap%rps, 1) ) THEN
IF( mesh > size(ap%rps, 1) ) THEN
ierr = 3
RETURN
END IF
IF( nwf > SIZE(ap%rps, 2) ) THEN
ENDIF
IF( nwf > size(ap%rps, 2) ) THEN
ierr = 4
RETURN
END IF
ENDIF
DO j = 0, lmax
READ(uni,fmt="(A16,i1)") strdum, l
READ(uni,'(4e16.8)') (vnl(i,j+1), i=1,mesh)
END DO
ENDDO
IF (nlcc) THEN
READ(uni,fmt='(4e16.8)') (rho_core(i), i=1,mesh)
END IF
ENDIF
DO j = 1, nwf
READ(uni,fmt="(A16,a2)") strdum,el(j)
READ(uni,fmt='(i5,f6.2)') ll(j),oc(j)
READ(uni,fmt='(4e16.8)') (chi(i,j), i=1,mesh)
END DO
ENDDO
ap%zv = zval
ap%nchan = lmax+1
@ -280,16 +280,16 @@ contains
! WRITE(6,*) ' DEBUG ', ap%lloc, ap%numeric, ap%nbeta, ap%raggio, ap%zv
DO i = 1, mesh
r = EXP(xmin+DBLE(i-1)*dx)/zmesh
r = exp(xmin+dble(i-1)*dx)/zmesh
ap%rw(i) = r
DO j = 1, lmax+1
ap%vnl(i,j) = vnl(i,j) * fac
END DO
END DO
IF( MINVAL( ap%rw(1:mesh) ) <= 0.0d0 ) THEN
ENDDO
ENDDO
IF( minval( ap%rw(1:mesh) ) <= 0.0d0 ) THEN
ierr = 5
RETURN
END IF
ENDIF
ap%dx = dx
ap%rab = ap%dx * ap%rw
ap%vloc(:) = ap%vnl(:,ap%lloc)
@ -302,27 +302,27 @@ contains
! fac = 1.0d0/SQRT(fpi)
fac = 1.0d0
DO i = 1, mesh
r = EXP(xmin+DBLE(i-1)*dx)/zmesh
r = exp(xmin+dble(i-1)*dx)/zmesh
DO j = 1, nwf
ap%rps(i,j) = chi(i,j) * fac
END DO
END DO
ENDDO
ENDDO
DO l = 1, ap%nbeta
il=ap%lll(l) + 1 ! find out the angular momentum (il-1) of the component stored
! in position l
DO i = 1, mesh
ap%vrps(i,l) = ( ap%vnl(i,il) - ap%vloc(i) ) * ap%rps(i,il)
END DO
END DO
ENDDO
ENDDO
IF( nlcc ) THEN
ap%rhoc = 0.0d0
DO i = 1, mesh
r = EXP(xmin+DBLE(i-1)*dx)/zmesh
r = exp(xmin+dble(i-1)*dx)/zmesh
ap%rhoc(i) = rho_core(i)
END DO
END IF
ENDDO
ENDIF
RETURN
END SUBROUTINE read_giannoz
@ -331,7 +331,7 @@ contains
SUBROUTINE ap_info( ap )
TYPE (pseudo_ncpp), INTENT(IN) :: ap
TYPE (pseudo_ncpp), INTENT(in) :: ap
INTEGER :: in1, in2, in3, in4, m, il, ib, l, i
IF (ap%nbeta > 0) THEN
@ -342,13 +342,13 @@ contains
WRITE(6,105) (ap%wgv(l),l=1,ap%nbeta)
ELSE
WRITE(6,50) ap%lloc
END IF
ENDIF
WRITE(6,60) (ap%lll(l),l=1,ap%nbeta)
ELSE
! ... A local pseudopotential has been read.
WRITE(6,11) ap%pottyp
WRITE(6,50) ap%lloc
END IF
ENDIF
10 FORMAT( 3X,'Type is ',A10,' and NONLOCAL. ')
107 FORMAT( 3X,'Mixed reference potential:')
@ -391,15 +391,15 @@ contains
DO il=1,3
DO ib=1,ap%igau
WRITE(6,103) ap%rcl(ib,il),ap%al(ib,il),ap%bl(ib,il)
END DO
END DO
ENDDO
ENDDO
40 FORMAT( 3X,'Hsc radii and coeff. A and B :')
103 FORMAT(3X,F8.4,2(3X,F15.7))
END IF
ENDIF
IF( ap%nrps > 0 .AND. ap%mesh > 0 ) THEN
IF( ap%nrps > 0 .and. ap%mesh > 0 ) THEN
WRITE(6,141) ap%nrps, ap%mesh, ap%dx
in1=1
in2=ap%mesh/4
@ -411,7 +411,7 @@ contains
WRITE(6,120) in2,ap%rw(in2),(ap%rps(in2,m),m=1,ap%nrps)
WRITE(6,120) in3,ap%rw(in3),(ap%rps(in3,m),m=1,ap%nrps)
WRITE(6,120) in4,ap%rw(in4),(ap%rps(in4,m),m=1,ap%nrps)
END IF
ENDIF
141 FORMAT(/, 3X,'Atomic wavefunction Grid : Channels = ',I2,&
', Mesh = ',I5,/,30X,'dx = ',F16.14)
@ -429,7 +429,7 @@ contains
WRITE(6,120) in2,ap%rw(in2),ap%rhoc(in2)
WRITE(6,120) in3,ap%rw(in3),ap%rhoc(in3)
WRITE(6,120) in4,ap%rw(in4),ap%rhoc(in4)
END IF
ENDIF
151 FORMAT(/, 3X,'Core correction Grid : Mesh = ',I5, &
', dx = ',F16.14)
@ -441,15 +441,15 @@ contains
!=----------------------------------------------------------------------------=!
REAL(DP) FUNCTION calculate_dx( a, m )
REAL(DP), INTENT(IN) :: a(:)
INTEGER, INTENT(IN) :: m
REAL(DP), INTENT(in) :: a(:)
INTEGER, INTENT(in) :: m
INTEGER :: n
REAL(DP) :: ra, rb
n = MIN( SIZE( a ), m )
n = min( size( a ), m )
ra = a(1)
rb = a(n)
calculate_dx = LOG( rb / ra ) / DBLE( n - 1 )
write(6,*) 'amesh (dx) = ', calculate_dx
calculate_dx = log( rb / ra ) / dble( n - 1 )
WRITE(6,*) 'amesh (dx) = ', calculate_dx
RETURN
END FUNCTION calculate_dx
@ -457,12 +457,12 @@ contains
SUBROUTINE read_atomic_wf( iunit, ap, err_msg, ierr)
USE parser, ONLY: field_count
IMPLICIT NONE
INTEGER, INTENT(IN) :: iunit
TYPE (pseudo_ncpp), INTENT(INOUT) :: ap
CHARACTER(LEN=*) :: err_msg
INTEGER, INTENT(OUT) :: ierr
INTEGER, INTENT(in) :: iunit
TYPE (pseudo_ncpp), INTENT(inout) :: ap
CHARACTER(len=*) :: err_msg
INTEGER, INTENT(out) :: ierr
!
CHARACTER(LEN=80) :: input_line
CHARACTER(len=80) :: input_line
INTEGER :: i, j, m, strlen, info, nf, mesh
REAL(DP) :: rdum
@ -488,35 +488,35 @@ SUBROUTINE read_atomic_wf( iunit, ap, err_msg, ierr)
IF( nf == 2 ) THEN
READ(input_line(1:strlen),*,IOSTAT=ierr) mesh, ap%nrps
ELSE
READ(input_line(1:strlen),*,IOSTAT=ierr) mesh, ap%nrps, ( ap%oc(j), j=1, MIN(ap%nrps,SIZE(ap%oc)) )
END IF
IF( ap%nrps > SIZE(ap%rps,2) ) THEN
READ(input_line(1:strlen),*,IOSTAT=ierr) mesh, ap%nrps, ( ap%oc(j), j=1, min(ap%nrps,size(ap%oc)) )
ENDIF
IF( ap%nrps > size(ap%rps,2) ) THEN
ierr = 2
WRITE( 6, * ) ' nchan = (wf) ', ap%nrps
err_msg = ' NCHAN NOT PROGRAMMED '
GO TO 110
END IF
IF( mesh > SIZE(ap%rw) .OR. mesh < 0) THEN
GOTO 110
ENDIF
IF( mesh > size(ap%rw) .or. mesh < 0) THEN
ierr = 4
err_msg = ' WAVMESH OUT OF RANGE '
GO TO 110
END IF
GOTO 110
ENDIF
DO j = 1, mesh
READ(iunit,*,IOSTAT=ierr) rdum, (ap%rps(j,m),m=1,ap%nrps)
IF( ap%mesh == 0 ) ap%rw(j) = rdum
IF( ABS(rdum - ap%rw(j))/(rdum+ap%rw(j)) > TOLMESH ) THEN
IF( abs(rdum - ap%rw(j))/(rdum+ap%rw(j)) > TOLMESH ) THEN
ierr = 5
err_msg = ' radial meshes do not match '
GO TO 110
END IF
END DO
GOTO 110
ENDIF
ENDDO
IF( ap%mesh == 0 ) THEN
ap%mesh = mesh
ap%dx = calculate_dx( ap%rw, ap%mesh )
ap%rab = ap%dx * ap%rw
END IF
ENDIF
GOTO 110
100 ierr = 1
@ -529,12 +529,12 @@ END SUBROUTINE read_atomic_wf
SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr)
IMPLICIT NONE
INTEGER, INTENT(IN) :: iunit
TYPE (pseudo_ncpp), INTENT(INOUT) :: ap
CHARACTER(LEN=*) :: err_msg
INTEGER, INTENT(OUT) :: ierr
INTEGER, INTENT(in) :: iunit
TYPE (pseudo_ncpp), INTENT(inout) :: ap
CHARACTER(len=*) :: err_msg
INTEGER, INTENT(out) :: ierr
!
CHARACTER(LEN=80) :: input_line
CHARACTER(len=80) :: input_line
INTEGER :: i, j, m, strlen, info, nf, l, ll
! ... read numeric atomic pseudopotential
@ -545,22 +545,22 @@ SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr)
IF(ap%tmix) THEN
READ(iunit,*) (ap%wgv(l),l=1,ap%nbeta)
END IF
ENDIF
READ(iunit,*,IOSTAT=ierr) ap%zv
READ(iunit,*,IOSTAT=ierr) ap%mesh, ap%nchan
IF((ap%nchan > SIZE(ap%vnl,2) ) .OR. (ap%nchan < 1)) THEN
IF((ap%nchan > size(ap%vnl,2) ) .or. (ap%nchan < 1)) THEN
ierr = 1
WRITE( 6, * ) ' nchan (pp) = ', ap%nchan
err_msg = ' NCHAN NOT PROGRAMMED '
GO TO 110
END IF
IF((ap%mesh > SIZE(ap%rw) ) .OR. (ap%mesh < 0)) THEN
GOTO 110
ENDIF
IF((ap%mesh > size(ap%rw) ) .or. (ap%mesh < 0)) THEN
info = 2
err_msg = ' NPOTMESH OUT OF RANGE '
GO TO 110
END IF
GOTO 110
ENDIF
ap%rw = 0.0d0
ap%vnl = 0.0d0
@ -568,13 +568,13 @@ SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr)
ap%vrps = 0.0d0
DO j = 1, ap%mesh
READ(iunit,*,IOSTAT=ierr) ap%rw(j), (ap%vnl(j,l),l=1,ap%nchan)
END DO
ENDDO
IF( MINVAL( ap%rw(1:ap%mesh) ) <= 0.0d0 ) THEN
IF( minval( ap%rw(1:ap%mesh) ) <= 0.0d0 ) THEN
info = 30
err_msg = ' ap rw too small '
GO TO 110
END IF
GOTO 110
ENDIF
! ... mixed reference potential is in vr(lloc)
IF(ap%tmix) THEN
@ -583,26 +583,26 @@ SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr)
DO l=1,ap%nchan
IF(l /= ap%lloc) THEN
ap%vnl(j,ap%lloc)= ap%vnl(j,ap%lloc) + ap%wgv(l) * ap%vnl(j,l)
END IF
END DO
END DO
END IF
ENDIF
ENDDO
ENDDO
ENDIF
ap%vloc(:) = ap%vnl(:,ap%lloc)
ap%dx = calculate_dx( ap%rw, ap%mesh )
ap%rab = ap%dx * ap%rw
CALL read_atomic_wf( iunit, ap, err_msg, ierr)
IF( ierr /= 0 ) GO TO 110
IF( ierr /= 0 ) GOTO 110
DO l = 1, ap%nbeta
ll=ap%lll(l) + 1
ap%vrps(:,l) = ( ap%vnl(:,ll) - ap%vloc(:) ) * ap%rps(:,ll)
END DO
ENDDO
IF(ap%tnlcc) THEN
CALL read_atomic_cc( iunit, ap, err_msg, ierr)
IF( ierr /= 0 ) GO TO 110
END IF
IF( ierr /= 0 ) GOTO 110
ENDIF
GOTO 110
100 ierr = 1
@ -615,10 +615,10 @@ END SUBROUTINE read_numeric_pp
SUBROUTINE read_head_pp( iunit, ap, err_msg, ierr)
IMPLICIT NONE
INTEGER, INTENT(IN) :: iunit
TYPE (pseudo_ncpp), INTENT(INOUT) :: ap
CHARACTER(LEN=*) :: err_msg
INTEGER, INTENT(OUT) :: ierr
INTEGER, INTENT(in) :: iunit
TYPE (pseudo_ncpp), INTENT(inout) :: ap
CHARACTER(len=*) :: err_msg
INTEGER, INTENT(out) :: ierr
!
INTEGER :: i, l
@ -629,39 +629,39 @@ SUBROUTINE read_head_pp( iunit, ap, err_msg, ierr)
ap%lll = 0
READ(iunit, *) ap%tnlcc, ap%tmix
READ(iunit, *) ap%pottyp, ap%lloc, ap%nbeta, (ap%lll(l), l = 1, MIN(ap%nbeta, SIZE(ap%lll)) )
READ(iunit, *) ap%pottyp, ap%lloc, ap%nbeta, (ap%lll(l), l = 1, min(ap%nbeta, size(ap%lll)) )
ap%lll = ap%lll - 1
IF( ap%nbeta > SIZE(ap%lll) .OR. ap%nbeta < 0 ) THEN
IF( ap%nbeta > size(ap%lll) .or. ap%nbeta < 0 ) THEN
ierr = 1
err_msg = 'LNL out of range'
GO TO 110
END IF
IF( ap%lloc < 0 .OR. ap%lloc > SIZE(ap%vnl,2) ) THEN
GOTO 110
ENDIF
IF( ap%lloc < 0 .or. ap%lloc > size(ap%vnl,2) ) THEN
ierr = 3
err_msg = 'LLOC out of range'
GO TO 110
END IF
IF( ap%tmix .AND. ap%pottyp /= 'NUMERIC' ) THEN
GOTO 110
ENDIF
IF( ap%tmix .and. ap%pottyp /= 'NUMERIC' ) THEN
ierr = 4
err_msg = 'tmix not implemented for pseudo ' // ap%pottyp
GO TO 110
END IF
GOTO 110
ENDIF
DO l = 2, ap%nbeta
IF( ap%lll(l) <= ap%lll(l-1)) THEN
ierr = 5
err_msg =' NONLOCAL COMPONENTS MUST BE GIVEN IN ASCENDING ORDER'
GO TO 110
END IF
END DO
GOTO 110
ENDIF
ENDDO
DO l = 1, ap%nbeta
IF( ap%lll(l)+1 == ap%lloc) THEN
ierr = 6
err_msg = ' LLOC.EQ.L NON LOCAL!!'
GO TO 110
END IF
END DO
GOTO 110
ENDIF
ENDDO
GOTO 110
100 ierr = 1
@ -674,10 +674,10 @@ END SUBROUTINE read_head_pp
SUBROUTINE read_analytic_pp( iunit, ap, err_msg, ierr)
IMPLICIT NONE
INTEGER, INTENT(IN) :: iunit
TYPE (pseudo_ncpp), INTENT(INOUT) :: ap
CHARACTER(LEN=*) :: err_msg
INTEGER, INTENT(OUT) :: ierr
INTEGER, INTENT(in) :: iunit
TYPE (pseudo_ncpp), INTENT(inout) :: ap
CHARACTER(len=*) :: err_msg
INTEGER, INTENT(out) :: ierr
!
INTEGER :: i, l
@ -708,22 +708,22 @@ SUBROUTINE read_analytic_pp( iunit, ap, err_msg, ierr)
CASE DEFAULT
ierr = 1
err_msg = ' IGAU NOT PROGRAMMED '
GO TO 110
GOTO 110
END SELECT
DO l=1,3
DO i=1,ap%igau
READ(iunit,*,IOSTAT=ierr) ap%rcl(i,l), ap%al(i,l), ap%bl(i,l)
END DO
END DO
ENDDO
ENDDO
CALL read_atomic_wf( iunit, ap, err_msg, ierr)
IF( ierr /= 0 ) GO TO 110
IF( ierr /= 0 ) GOTO 110
IF(ap%tnlcc) THEN
CALL read_atomic_cc( iunit, ap, err_msg, ierr)
IF( ierr /= 0 ) GO TO 110
END IF
IF( ierr /= 0 ) GOTO 110
ENDIF
! ... Analytic pseudo are not supported anymore, conversion
! ... to numeric form is forced
@ -741,12 +741,12 @@ END SUBROUTINE read_analytic_pp
SUBROUTINE read_atomic_cc( iunit, ap, err_msg, ierr)
IMPLICIT NONE
INTEGER, INTENT(IN) :: iunit
TYPE (pseudo_ncpp), INTENT(INOUT) :: ap
CHARACTER(LEN=*) :: err_msg
INTEGER, INTENT(OUT) :: ierr
INTEGER, INTENT(in) :: iunit
TYPE (pseudo_ncpp), INTENT(inout) :: ap
CHARACTER(len=*) :: err_msg
INTEGER, INTENT(out) :: ierr
!
CHARACTER(LEN=80) :: input_line
CHARACTER(len=80) :: input_line
INTEGER :: j, mesh
REAL(DP) :: rdum
@ -758,26 +758,26 @@ SUBROUTINE read_atomic_cc( iunit, ap, err_msg, ierr)
ap%rhoc = 0.0d0
READ(iunit,*,IOSTAT=ierr) mesh
IF(mesh > SIZE(ap%rw) .OR. mesh < 0 ) THEN
IF(mesh > size(ap%rw) .or. mesh < 0 ) THEN
ierr = 17
err_msg = ' CORE CORRECTION MESH OUT OF RANGE '
GO TO 110
END IF
GOTO 110
ENDIF
DO j = 1, mesh
READ(iunit,*,IOSTAT=ierr) rdum, ap%rhoc(j)
IF( ap%mesh == 0 ) ap%rw(j) = rdum
IF( ABS(rdum - ap%rw(j))/(rdum+ap%rw(j)) > TOLMESH ) THEN
IF( abs(rdum - ap%rw(j))/(rdum+ap%rw(j)) > TOLMESH ) THEN
ierr = 5
err_msg = ' core cor. radial mesh does not match '
GO TO 110
END IF
END DO
GOTO 110
ENDIF
ENDDO
IF( ap%mesh == 0 ) THEN
ap%mesh = mesh
ap%dx = calculate_dx( ap%rw, ap%mesh )
ap%rab = ap%dx * ap%rw
END IF
ENDIF
GOTO 110
100 ierr = 1
@ -787,12 +787,12 @@ SUBROUTINE read_atomic_cc( iunit, ap, err_msg, ierr)
END SUBROUTINE read_atomic_cc
end module fpmd2upf_module
END MODULE fpmd2upf_module
program fpmd2upf
PROGRAM fpmd2upf
!
! Convert a pseudopotential written in the FPMD format
@ -807,39 +807,39 @@ program fpmd2upf
IMPLICIT NONE
TYPE (pseudo_ncpp) :: ap
CHARACTER(LEN=256) :: psfile
CHARACTER(LEN=2) :: wfl( 10 )
CHARACTER(len=256) :: psfile
CHARACTER(len=2) :: wfl( 10 )
REAL(8) :: wfoc( 10 )
INTEGER :: nsp, nspnl, i, lloc, l, ir, iv, kkbeta
REAL(8) :: rmax = 10
REAL(8) :: vll
REAL(8), allocatable :: aux(:)
REAL(8), ALLOCATABLE :: aux(:)
namelist / fpmd_pseudo / psfile, nwfs, wfl, wfoc, psd, &
NAMELIST / fpmd_pseudo / psfile, nwfs, wfl, wfoc, psd, &
iexch, icorr, igcx, igcc, zp
! ... end of declarations
call input_from_file()
CALL input_from_file()
read( 5, fpmd_pseudo )
READ( 5, fpmd_pseudo )
nsp = 1
CALL read_pseudo_fpmd(ap, psfile)
write(generated, '("Generated using unknown code")')
write(date_author,'("Author: unknown Generation date: as well")')
WRITE(generated, '("Generated using unknown code")')
WRITE(date_author,'("Author: unknown Generation date: as well")')
comment = 'Info: automatically converted from CPMD format'
rcloc = 0.0d0
allocate( els(nwfs), oc(nwfs), epseu(nwfs) )
allocate( lchi(nwfs), nns(nwfs) )
allocate( rcut (nwfs), rcutus (nwfs) )
ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs) )
ALLOCATE( lchi(nwfs), nns(nwfs) )
ALLOCATE( rcut (nwfs), rcutus (nwfs) )
els = '?'
oc = 0.0d0
do i = 1, nwfs
DO i = 1, nwfs
els(i) = wfl(i)
oc(i) = wfoc(i)
lchi(i) = i - 1
@ -847,110 +847,110 @@ program fpmd2upf
rcut(i) = 0.0d0
rcutus(i)= 0.0d0
epseu(i) = 0.0d0
end do
ENDDO
pseudotype = 'NC'
nlcc = ap%tnlcc
if( ap%zv > 0.0d0 ) zp = ap%zv
IF( ap%zv > 0.0d0 ) zp = ap%zv
etotps = 0.0d0
lloc = ap%lloc
lmax = MAX( MAXVAL( ap%lll( 1:ap%nbeta ) ), ap%lloc - 1 )
lmax = max( maxval( ap%lll( 1:ap%nbeta ) ), ap%lloc - 1 )
nbeta = ap%nbeta
mesh = ap%mesh
ntwfc = nwfs
allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
do i = 1, nwfs
ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
DO i = 1, nwfs
lchiw(i) = lchi(i)
ocw(i) = oc(i)
elsw(i) = els(i)
end do
ENDDO
allocate(rab(mesh))
allocate( r(mesh))
ALLOCATE(rab(mesh))
ALLOCATE( r(mesh))
r = ap%rw
ap%dx = calculate_dx( ap%rw, ap%mesh )
rab = ap%rw * ap%dx
write(6,*) ap%lloc, ap%lll( 1:ap%nbeta ) , ap%nbeta, ap%dx
WRITE(6,*) ap%lloc, ap%lll( 1:ap%nbeta ) , ap%nbeta, ap%dx
allocate (rho_atc(mesh))
if (nlcc) rho_atc = ap%rhoc
ALLOCATE (rho_atc(mesh))
IF (nlcc) rho_atc = ap%rhoc
allocate (vloc0(mesh))
ALLOCATE (vloc0(mesh))
! the factor 2 converts from Hartree to Rydberg
vloc0(:) = ap%vloc * 2.0d0
if (nbeta > 0) then
IF (nbeta > 0) THEN
allocate(ikk2(nbeta), lll(nbeta))
ALLOCATE(ikk2(nbeta), lll(nbeta))
kkbeta = mesh
do ir = 1,mesh
if ( r(ir) > rmax ) then
DO ir = 1,mesh
IF ( r(ir) > rmax ) THEN
kkbeta=ir
exit
end if
end do
ENDIF
ENDDO
ikk2(:) = kkbeta
allocate(aux(kkbeta))
allocate(betar(mesh,nbeta))
allocate(qfunc(mesh,nbeta,nbeta))
allocate(dion(nbeta,nbeta))
allocate(qqq (nbeta,nbeta))
ALLOCATE(aux(kkbeta))
ALLOCATE(betar(mesh,nbeta))
ALLOCATE(qfunc(mesh,nbeta,nbeta))
ALLOCATE(dion(nbeta,nbeta))
ALLOCATE(qqq (nbeta,nbeta))
qfunc(:,:,:)=0.0d0
dion(:,:) =0.d0
qqq(:,:) =0.d0
iv = 0
do i = 1, nwfs
DO i = 1, nwfs
l = lchi(i)
if ( l .ne. (lloc-1) ) then
IF ( l /= (lloc-1) ) THEN
iv = iv + 1
lll( iv ) = l
do ir = 1, kkbeta
DO ir = 1, kkbeta
! the factor 2 converts from Hartree to Rydberg
betar(ir, iv) = 2.d0 * ap%vrps( ir, iv )
aux(ir) = ap%rps(ir, (l+1) ) * betar(ir, iv)
end do
call simpson2(kkbeta, aux(1), rab(1), vll)
ENDDO
CALL simpson2(kkbeta, aux(1), rab(1), vll)
dion(iv,iv) = 1.0d0/vll
write(6,*) aux(2), rab(2), kkbeta, vll
end if
enddo
WRITE(6,*) aux(2), rab(2), kkbeta, vll
ENDIF
ENDDO
end if
ENDIF
allocate (rho_at(mesh))
ALLOCATE (rho_at(mesh))
rho_at = 0.d0
do i = 1, nwfs
DO i = 1, nwfs
rho_at(:) = rho_at(:) + ocw(i) * ap%rps(:, i) ** 2
end do
ENDDO
allocate (chi(mesh,ntwfc))
ALLOCATE (chi(mesh,ntwfc))
chi = ap%rps
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential successfully converted'
WRITE (6,'(a)') 'Pseudopotential successfully converted'
! ----------------------------------------------------------
call write_upf( 10 )
CALL write_upf( 10 )
100 continue
100 CONTINUE
end program fpmd2upf
END PROGRAM fpmd2upf
!----------------------------------------------------------------------
subroutine simpson2(mesh,func,rab,asum)
SUBROUTINE simpson2(mesh,func,rab,asum)
!-----------------------------------------------------------------------
!
! simpson's rule integrator for function stored on the
! radial logarithmic mesh
!
implicit none
IMPLICIT NONE
integer :: i, mesh
INTEGER :: i, mesh
real(8) :: rab(mesh), func(mesh), f1, f2, f3, r12, asum
! routine assumes that mesh is an odd number so run check
@ -964,15 +964,15 @@ subroutine simpson2(mesh,func,rab,asum)
r12 = 1.0d0 / 12.0d0
f3 = func(1) * rab(1) * r12
do i = 2,mesh-1,2
DO i = 2,mesh-1,2
f1 = f3
f2 = func(i) * rab(i) * r12
f3 = func(i+1) * rab(i+1) * r12
asum = asum + 4.0d0*f1 + 16.0d0*f2 + 4.0d0*f3
enddo
ENDDO
return
end subroutine simpson2
RETURN
END SUBROUTINE simpson2

392
upftools/ncpp2upf.f90
View File

@ -7,281 +7,281 @@
!
!
!---------------------------------------------------------------------
program ncpp2upf
PROGRAM ncpp2upf
!---------------------------------------------------------------------
!
! Convert a pseudopotential written in PWSCF format
! (norm-conserving) to unified pseudopotential format
implicit none
character(len=256) filein, fileout
IMPLICIT NONE
CHARACTER(len=256) filein, fileout
!
!
call get_file ( filein )
open(unit=1,file=filein,status='old',form='formatted')
call read_ncpp(1)
close (unit=1)
CALL get_file ( filein )
OPEN(unit=1,file=filein,status='old',form='formatted')
CALL read_ncpp(1)
CLOSE (unit=1)
! convert variables read from NCPP format into those needed
! by the upf format - add missing quantities
call convert_ncpp
CALL convert_ncpp
fileout=trim(filein)//'.UPF'
print '(''Output PP file in US format : '',a)', fileout
PRINT '(''Output PP file in US format : '',a)', fileout
open(unit=2,file=fileout,status='unknown',form='formatted')
call write_upf(2)
close (unit=2)
stop
20 call errore ('ncpp2upf', 'Reading pseudo file name ', 1)
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
CALL write_upf(2)
CLOSE (unit=2)
STOP
20 CALL errore ('ncpp2upf', 'Reading pseudo file name ', 1)
end program ncpp2upf
END PROGRAM ncpp2upf
module ncpp
MODULE ncpp
!
! All variables read from NCPP file format
!
! trailing underscore means that a variable with the same name
! is used in module 'upf' containing variables to be written
!
character(len=20) :: dft_
character(len=2) :: psd_
CHARACTER(len=20) :: dft_
CHARACTER(len=2) :: psd_
real(8) :: zp_
integer nlc, nnl, lmax_, lloc, nchi
logical :: numeric, bhstype, nlcc_
INTEGER nlc, nnl, lmax_, lloc, nchi
LOGICAL :: numeric, bhstype, nlcc_
real(8) :: alpc(2), cc(2), alps(3,0:3), aps(6,0:3)
real(8) :: a_nlcc, b_nlcc, alpha_nlcc
real(8) :: zmesh, xmin, dx
real(8), allocatable:: r_(:), rab_(:)
integer :: mesh_
real(8), ALLOCATABLE:: r_(:), rab_(:)
INTEGER :: mesh_
real(8), allocatable:: vnl(:,:), rho_atc_(:), rho_at_(:)
integer, allocatable:: lchi_(:)
real(8), allocatable:: chi_(:,:), oc_(:)
real(8), ALLOCATABLE:: vnl(:,:), rho_atc_(:), rho_at_(:)
INTEGER, ALLOCATABLE:: lchi_(:)
real(8), ALLOCATABLE:: chi_(:,:), oc_(:)
end module ncpp
END MODULE ncpp
!
! ----------------------------------------------------------
subroutine read_ncpp(iunps)
SUBROUTINE read_ncpp(iunps)
! ----------------------------------------------------------
!
use ncpp
use upf , only : els
implicit none
integer :: iunps
USE ncpp
USE upf , ONLY : els
IMPLICIT NONE
INTEGER :: iunps
!
character(len=1), dimension(0:3) :: convel=(/'S','P','D','F'/)
character(len=2) :: label
CHARACTER(len=1), DIMENSION(0:3) :: convel=(/'S','P','D','F'/)
CHARACTER(len=2) :: label
real (8) :: x, qe_erf
integer :: l, i, ir, nb, n
character (len=255) line
external qe_erf
INTEGER :: l, i, ir, nb, n
CHARACTER (len=255) line
EXTERNAL qe_erf
read(iunps, '(a)', end=300, err=300 ) dft_
if (dft_(1:2).eq.'**') dft_ = 'PZ'
READ(iunps, '(a)', end=300, err=300 ) dft_
IF (dft_(1:2)=='**') dft_ = 'PZ'
read (iunps, *, err=300) psd_, zp_, lmax_, nlc, nnl, nlcc_, &
READ (iunps, *, err=300) psd_, zp_, lmax_, nlc, nnl, nlcc_, &
lloc, bhstype
if ( nlc.gt.2 .or. nnl.gt.3) &
call errore( 'read_ncpp','Wrong nlc or nnl',1 )
if ( nlc* nnl .lt. 0 ) &
call errore( 'read_ncpp','nlc*nnl < 0 ? ',1 )
if ( zp_.le.0d0 ) &
call errore( 'read_ncpp','Wrong zp ',1 )
if ( lmax_.gt.3.or.lmax_.lt.0 ) &
call errore( 'read_ncpp','Wrong lmax ',1 )
IF ( nlc>2 .or. nnl>3) &
CALL errore( 'read_ncpp','Wrong nlc or nnl',1 )
IF ( nlc* nnl < 0 ) &
CALL errore( 'read_ncpp','nlc*nnl < 0 ? ',1 )
IF ( zp_<=0d0 ) &
CALL errore( 'read_ncpp','Wrong zp ',1 )
IF ( lmax_>3.or.lmax_<0 ) &
CALL errore( 'read_ncpp','Wrong lmax ',1 )
if (lloc.eq.-1000) lloc=lmax_
IF (lloc==-1000) lloc=lmax_
!
! In numeric pseudopotentials both nlc and nnl are zero.
!
numeric = nlc.le.0 .and. nnl.le.0
numeric = nlc<=0 .and. nnl<=0
if (.not.numeric) then
IF (.not.numeric) THEN
!
! read pseudopotentials in analytic form
!
read(iunps, *, err=300) &
READ(iunps, *, err=300) &
( alpc(i), i=1, 2 ), ( cc(i), i=1,2 )
if ( abs(cc(1)+cc(2)-1.d0).gt.1.0d-6) &
call errore ('read_ncpp','wrong pseudopotential coefficients',1)
do l = 0, lmax_
read (iunps, *, err=300) &
IF ( abs(cc(1)+cc(2)-1.d0)>1.0d-6) &
CALL errore ('read_ncpp','wrong pseudopotential coefficients',1)
DO l = 0, lmax_
READ (iunps, *, err=300) &
( alps(i,l),i=1,3 ), (aps(i,l),i=1,6)
enddo
ENDDO
if (nlcc_) then
read(iunps, *, err=300) &
IF (nlcc_) THEN
READ(iunps, *, err=300) &
a_nlcc, b_nlcc, alpha_nlcc
if (alpha_nlcc.le.0.d0) &
call errore('read_ncpp','nlcc but alpha=0',1)
end if
IF (alpha_nlcc<=0.d0) &
CALL errore('read_ncpp','nlcc but alpha=0',1)
ENDIF
if (bhstype) call bachel(alps,aps,1,lmax_)
end if
IF (bhstype) CALL bachel(alps,aps,1,lmax_)
ENDIF
read(iunps, *, err=300) zmesh, xmin, dx, mesh_, nchi
READ(iunps, *, err=300) zmesh, xmin, dx, mesh_, nchi
if ( mesh_.le.0) call errore( 'read_ncpp', 'mesh too small', 1)
if ( (nchi.lt.lmax_ .and. lloc.eq.lmax_).or. &
(nchi.lt.lmax_+1 .and. lloc.ne.lmax_) ) &
call errore( 'read_ncpp', 'wrong no. of wfcts', 1 )
IF ( mesh_<=0) CALL errore( 'read_ncpp', 'mesh too small', 1)
IF ( (nchi<lmax_ .and. lloc==lmax_).or. &
(nchi<lmax_+1 .and. lloc/=lmax_) ) &
CALL errore( 'read_ncpp', 'wrong no. of wfcts', 1 )
!
! compute the radial mesh
!
allocate( r_(mesh_))
allocate(rab_(mesh_))
ALLOCATE( r_(mesh_))
ALLOCATE(rab_(mesh_))
do ir = 1, mesh_
x = xmin + DBLE(ir-1) * dx
DO ir = 1, mesh_
x = xmin + dble(ir-1) * dx
r_ (ir) = exp(x) / zmesh
rab_(ir) = dx * r_(ir)
end do
ENDDO
allocate(vnl(mesh_,0:lmax_))
if (numeric) then
ALLOCATE(vnl(mesh_,0:lmax_))
IF (numeric) THEN
!
! read pseudopotentials in numeric form
!
do l = 0, lmax_
read(iunps, '(a)', err=300)
read(iunps, *, err=300) (vnl(ir,l),ir=1,mesh_)
enddo
DO l = 0, lmax_
READ(iunps, '(a)', err=300)
READ(iunps, *, err=300) (vnl(ir,l),ir=1,mesh_)
ENDDO
allocate(rho_atc_(mesh_))
if(nlcc_) then
read(iunps, *, err=300) ( rho_atc_(ir), ir=1,mesh_ )
endif
ALLOCATE(rho_atc_(mesh_))
IF(nlcc_) THEN
READ(iunps, *, err=300) ( rho_atc_(ir), ir=1,mesh_ )
ENDIF
else
ELSE
!
! convert analytic to numeric form
!
do l=0,lmax_
DO l=0,lmax_
!
! DO NOT USE f90 ARRAY SYNTAX: qe_erf IS NOT AN INTRINSIC FUNCTION!!!
!
do ir=1,mesh_
DO ir=1,mesh_
vnl(ir,l)= - ( cc(1)*qe_erf(sqrt(alpc(1))*r_(ir)) + &
cc(2)*qe_erf(sqrt(alpc(2))*r_(ir)) ) * zp_/r_(ir)
end do
ENDDO
do n=1,nnl
DO n=1,nnl
vnl(:,l)= vnl(:,l)+ (aps(n,l)+ aps(n+3,l)*r_(:)**2 )* &
exp(-alps(n,l)*r_(:)**2)
end do
ENDDO
!
! convert to Rydberg
!
vnl(:,l) = vnl(:,l)*2.0d0
end do
ENDDO
allocate(rho_atc_(mesh_))
if (nlcc_) then
ALLOCATE(rho_atc_(mesh_))
IF (nlcc_) THEN
rho_atc_(:) =(a_nlcc+b_nlcc*(r_(:)**2))*exp(-alpha_nlcc*r_(:)**2)
where(abs(rho_atc_) < 1.0d-15)
WHERE(abs(rho_atc_) < 1.0d-15)
rho_atc_ = 0
end where
end if
endif
END WHERE
ENDIF
ENDIF
!
! subtract the local part
!
do l = 0, lmax_
if ( l.ne.lloc ) vnl(:,l) = vnl(:,l) - vnl(:,lloc)
enddo
DO l = 0, lmax_
IF ( l/=lloc ) vnl(:,l) = vnl(:,l) - vnl(:,lloc)
ENDDO
!
! read pseudowavefunctions
!
allocate(lchi_(nchi), els(nchi))
allocate(oc_(nchi))
allocate(chi_(mesh_,nchi))
do nb = 1, nchi
ALLOCATE(lchi_(nchi), els(nchi))
ALLOCATE(oc_(nchi))
ALLOCATE(chi_(mesh_,nchi))
DO nb = 1, nchi
! read wavefunction label and store for later
read(iunps, '(a)', err=300) line
read(iunps, *, err=300) lchi_( nb), oc_( nb )
READ(iunps, '(a)', err=300) line
READ(iunps, *, err=300) lchi_( nb), oc_( nb )
!
! Test lchi and occupation numbers
!
if ( nb.le.lmax_.and.lchi_(nb)+1.ne.nb) &
call errore('read_ncpp','order of wavefunctions',nb)
if (lchi_(nb).gt.lmax_ .or. lchi_(nb).lt.0) &
call errore('read_ncpp','wrong lchi',nb)
if ( oc_(nb).lt.0.d0 .or. &
oc_(nb).gt.2.d0*(2*lchi_(nb)+1)) &
call errore('read_ncpp','wrong oc',nb)
IF ( nb<=lmax_.and.lchi_(nb)+1/=nb) &
CALL errore('read_ncpp','order of wavefunctions',nb)
IF (lchi_(nb)>lmax_ .or. lchi_(nb)<0) &
CALL errore('read_ncpp','wrong lchi',nb)
IF ( oc_(nb)<0.d0 .or. &
oc_(nb)>2.d0*(2*lchi_(nb)+1)) &
CALL errore('read_ncpp','wrong oc',nb)
!
! parse and check wavefunction label
read(line,'(14x,a2)', err=222, end=222) label
if (label(2:2).ne.convel(lchi_(nb))) goto 222
do l = 0, lmax_
if (label(2:2).eq.convel(l)) then
READ(line,'(14x,a2)', err=222, end=222) label
IF (label(2:2)/=convel(lchi_(nb))) GOTO 222
DO l = 0, lmax_
IF (label(2:2)==convel(l)) THEN
els(nb) = label(1:2)
goto 223
endif
end do
222 continue
GOTO 223
ENDIF
ENDDO
222 CONTINUE
els(nb) = '*'//convel(lchi_(nb))
223 continue
223 CONTINUE
!
! finally read the wavefunction
read(iunps, *, err=300) (chi_(ir,nb),ir=1,mesh_)
enddo
READ(iunps, *, err=300) (chi_(ir,nb),ir=1,mesh_)
ENDDO
!
! compute the atomic charges
!
allocate(rho_at_(mesh_))
ALLOCATE(rho_at_(mesh_))
rho_at_(:)=0.d0
do nb = 1, nchi
if( oc_(nb).ne.0.d0) &
DO nb = 1, nchi
IF( oc_(nb)/=0.d0) &
rho_at_(:) = rho_at_(:) + oc_(nb)*chi_(:,nb)**2
end do
ENDDO
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential successfully read'
WRITE (6,'(a)') 'Pseudopotential successfully read'
! ----------------------------------------------------------
return
RETURN
300 call errore('read_ncpp','pseudo file is empty or wrong',1)
300 CALL errore('read_ncpp','pseudo file is empty or wrong',1)
end subroutine read_ncpp
END SUBROUTINE read_ncpp
! ----------------------------------------------------------
subroutine convert_ncpp
SUBROUTINE convert_ncpp
! ----------------------------------------------------------
use ncpp
use upf
use funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc
implicit none
real(8), parameter :: rmax = 10.0d0
real(8), allocatable :: aux(:)
USE ncpp
USE upf
USE funct, ONLY : set_dft_from_name, get_iexch, get_icorr, get_igcx, get_igcc
IMPLICIT NONE
real(8), PARAMETER :: rmax = 10.0d0
real(8), ALLOCATABLE :: aux(:)
real(8) :: vll
integer :: kkbeta, l, iv, ir, i
INTEGER :: kkbeta, l, iv, ir, i
write(generated, '("Generated using ld1 code (maybe, or maybe not)")')
write(date_author,'("Author: unknown Generation date: as well")')
WRITE(generated, '("Generated using ld1 code (maybe, or maybe not)")')
WRITE(date_author,'("Author: unknown Generation date: as well")')
comment = 'Info: automatically converted from PWSCF format'
! reasonable assumption
if (zmesh > 18) then
IF (zmesh > 18) THEN
rel = 1
else
ELSE
rel = 0
end if
ENDIF
rcloc = 0.0d0
nwfs = nchi
allocate( oc(nwfs), epseu(nwfs))
allocate(lchi(nwfs), nns(nwfs) )
allocate(rcut (nwfs), rcutus (nwfs))
do i=1, nwfs
ALLOCATE( oc(nwfs), epseu(nwfs))
ALLOCATE(lchi(nwfs), nns(nwfs) )
ALLOCATE(rcut (nwfs), rcutus (nwfs))
DO i=1, nwfs
nns (i) = 0
lchi(i) = lchi_(i)
rcut(i) = 0.0d0
rcutus(i)= 0.0d0
oc (i) = oc_(i)
epseu(i) = 0.0d0
end do
deallocate (lchi_, oc_)
ENDDO
DEALLOCATE (lchi_, oc_)
psd = psd_
pseudotype = 'NC'
@ -290,96 +290,96 @@ subroutine convert_ncpp
etotps = 0.0d0
ecutrho=0.0d0
ecutwfc=0.0d0
if ( lmax_ == lloc) then
IF ( lmax_ == lloc) THEN
lmax = lmax_-1
else
ELSE
lmax = lmax_
end if
ENDIF
nbeta= lmax_
mesh = mesh_
ntwfc= nchi
allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
do i=1, nchi
ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
DO i=1, nchi
lchiw(i) = lchi(i)
ocw(i) = oc(i)
elsw(i) = els(i)
end do
call set_dft_from_name(dft_)
ENDDO
CALL set_dft_from_name(dft_)
iexch = get_iexch()
icorr = get_icorr()
igcx = get_igcx()
igcc = get_igcc()
allocate(rab(mesh))
allocate( r(mesh))
ALLOCATE(rab(mesh))
ALLOCATE( r(mesh))
rab = rab_
r = r_
allocate (rho_atc(mesh))
ALLOCATE (rho_atc(mesh))
rho_atc = rho_atc_
deallocate (rho_atc_)
DEALLOCATE (rho_atc_)
allocate (vloc0(mesh))
ALLOCATE (vloc0(mesh))
vloc0(:) = vnl(:,lloc)
if (nbeta > 0) then
IF (nbeta > 0) THEN
allocate(ikk2(nbeta), lll(nbeta))
ALLOCATE(ikk2(nbeta), lll(nbeta))
kkbeta=mesh
do ir = 1,mesh
if ( r(ir) > rmax ) then
DO ir = 1,mesh
IF ( r(ir) > rmax ) THEN
kkbeta=ir
exit
end if
end do
ENDIF
ENDDO
! make sure kkbeta is odd as required for simpson
if(mod(kkbeta,2) == 0) kkbeta=kkbeta-1
IF(mod(kkbeta,2) == 0) kkbeta=kkbeta-1
ikk2(:) = kkbeta
allocate(aux(kkbeta))
allocate(betar(mesh,nbeta))
allocate(qfunc(mesh,nbeta,nbeta))
allocate(dion(nbeta,nbeta))
allocate(qqq (nbeta,nbeta))
ALLOCATE(aux(kkbeta))
ALLOCATE(betar(mesh,nbeta))
ALLOCATE(qfunc(mesh,nbeta,nbeta))
ALLOCATE(dion(nbeta,nbeta))
ALLOCATE(qqq (nbeta,nbeta))
qfunc(:,:,:)=0.0d0
dion(:,:) =0.d0
qqq(:,:) =0.d0
iv=0
do i=1,nchi
DO i=1,nchi
l=lchi(i)
if (l.ne.lloc) then
IF (l/=lloc) THEN
iv=iv+1
lll(iv)=l
do ir=1,kkbeta
DO ir=1,kkbeta
betar(ir,iv)=chi_(ir,i)*vnl(ir,l)
aux(ir) = chi_(ir,i)**2*vnl(ir,l)
end do
call simpson(kkbeta,aux,rab,vll)
ENDDO
CALL simpson(kkbeta,aux,rab,vll)
dion(iv,iv) = 1.0d0/vll
end if
if(iv >= nbeta) exit ! skip additional pseudo wfns
enddo
deallocate (vnl, aux)
ENDIF
IF(iv >= nbeta) exit ! skip additional pseudo wfns
ENDDO
DEALLOCATE (vnl, aux)
!
! redetermine ikk2
!
do iv=1,nbeta
DO iv=1,nbeta
ikk2(iv)=kkbeta
do ir = kkbeta,1,-1
if ( abs(betar(ir,iv)) > 1.d-12 ) then
DO ir = kkbeta,1,-1
IF ( abs(betar(ir,iv)) > 1.d-12 ) THEN
ikk2(iv)=ir
exit
end if
end do
end do
end if
allocate (rho_at(mesh))
ENDIF
ENDDO
ENDDO
ENDIF
ALLOCATE (rho_at(mesh))
rho_at = rho_at_
deallocate (rho_at_)
DEALLOCATE (rho_at_)
allocate (chi(mesh,ntwfc))
ALLOCATE (chi(mesh,ntwfc))
chi = chi_
deallocate (chi_)
DEALLOCATE (chi_)
return
end subroutine convert_ncpp
RETURN
END SUBROUTINE convert_ncpp

224
upftools/oldcp2upf.f90
View File

@ -7,139 +7,139 @@
!
!
!---------------------------------------------------------------------
program oldcp2upf
PROGRAM oldcp2upf
!---------------------------------------------------------------------
!
! Convert a pseudopotential written in the old CP90 format
! (without core correction) to unified pseudopotential format
!
implicit none
character(len=256) filein, fileout
IMPLICIT NONE
CHARACTER(len=256) filein, fileout
!
!
call get_file ( filein )
open (unit = 1, file = filein, status = 'old', form = 'formatted')
call read_oldcp(1)
close (1)
CALL get_file ( filein )
OPEN (unit = 1, file = filein, status = 'old', form = 'formatted')
CALL read_oldcp(1)
CLOSE (1)
! convert variables read from old CP90 format into those needed
! by the upf format - add missing quantities
call convert_oldcp
CALL convert_oldcp
fileout=trim(filein)//'.UPF'
print '(''Output PP file in UPF format : '',a)', fileout
PRINT '(''Output PP file in UPF format : '',a)', fileout
open(unit=2,file=fileout,status='unknown',form='formatted')
call write_upf(2)
close (unit=2)
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
CALL write_upf(2)
CLOSE (unit=2)
stop
20 call errore ('oldcp2upf', 'Reading pseudo file name ', 1)
STOP
20 CALL errore ('oldcp2upf', 'Reading pseudo file name ', 1)
end program oldcp2upf
END PROGRAM oldcp2upf
module oldcp
MODULE oldcp
!
! All variables read from old CP90 file format
!
real(8) :: amesh, z, zv
integer :: exfact, lloc, nbeta_, mesh_
INTEGER :: exfact, lloc, nbeta_, mesh_
real(8) :: wrc1, rc1, wrc2, rc2, rcl(3,3), al(3,3), bl(3,3)
real(8), allocatable :: r_(:), vnl(:,:), chi_(:,:)
real(8), ALLOCATABLE :: r_(:), vnl(:,:), chi_(:,:)
!
!------------------------------
end module oldcp
END MODULE oldcp
!
! ----------------------------------------------------------
subroutine read_oldcp(iunps)
SUBROUTINE read_oldcp(iunps)
! ----------------------------------------------------------
!
use oldcp
implicit none
integer :: iunps
USE oldcp
IMPLICIT NONE
INTEGER :: iunps
!
real(8), external :: qe_erf
integer :: i, l, j, jj
real(8), EXTERNAL :: qe_erf
INTEGER :: i, l, j, jj
!
read(iunps,*, end=10, err=10) z, zv, nbeta_, lloc, exfact
if (z < 1 .or. z > 100 .or. zv < 1 .or. zv > 25 ) &
call errore ('read_oldcp','wrong potential read',1)
read(iunps,*, end=10, err=10) wrc1, rc1, wrc2, rc2
read(iunps,*, end=10, err=10) ( ( rcl(i,l), al(i,l), &
READ(iunps,*, end=10, err=10) z, zv, nbeta_, lloc, exfact
IF (z < 1 .or. z > 100 .or. zv < 1 .or. zv > 25 ) &
CALL errore ('read_oldcp','wrong potential read',1)
READ(iunps,*, end=10, err=10) wrc1, rc1, wrc2, rc2
READ(iunps,*, end=10, err=10) ( ( rcl(i,l), al(i,l), &
bl(i,l), i = 1, 3), l = 1, 3)
read(iunps,*, end=10, err=10) mesh_, amesh
allocate(r_(mesh_))
allocate (chi_(mesh_,nbeta_))
do l = 1, nbeta_
if (l > 1) read(iunps,*, end=10, err=10) mesh_, amesh
do j = 1, mesh_
read(iunps,*, end=10, err=10) jj, r_(j), chi_(j,l)
end do
end do
READ(iunps,*, end=10, err=10) mesh_, amesh
ALLOCATE(r_(mesh_))
ALLOCATE (chi_(mesh_,nbeta_))
DO l = 1, nbeta_
IF (l > 1) READ(iunps,*, end=10, err=10) mesh_, amesh
DO j = 1, mesh_
READ(iunps,*, end=10, err=10) jj, r_(j), chi_(j,l)
ENDDO
ENDDO
!
! convert analytic to numeric form
!
allocate (vnl(mesh_,0:nbeta_))
do l=0,nbeta_
ALLOCATE (vnl(mesh_,0:nbeta_))
DO l=0,nbeta_
!
! DO NOT USE f90 ARRAY SYNTAX: qe_erf IS NOT AN INTRINSIC FUNCTION!!!
!
do j=1, mesh_
DO j=1, mesh_
vnl(j,l)= - (wrc1*qe_erf(sqrt(rc1)*r_(j)) + &
wrc2*qe_erf(sqrt(rc2)*r_(j)) ) * zv/r_(j)
end do
ENDDO
!
do i=1,3
DO i=1,3
vnl(:,l)= vnl(:,l)+ (al(i,l+1)+ bl(i,l+1)*r_(:)**2) * &
exp(-rcl(i,l+1)*r_(:)**2)
end do
end do
ENDDO
ENDDO
return
10 call errore('read_oldcp','error in reading file',1)
RETURN
10 CALL errore('read_oldcp','error in reading file',1)
end subroutine read_oldcp
END SUBROUTINE read_oldcp
! ----------------------------------------------------------
subroutine convert_oldcp
SUBROUTINE convert_oldcp
! ----------------------------------------------------------
!
use oldcp
use upf
implicit none
real(8), parameter :: rmax = 10.0d0
real(8), allocatable :: aux(:)
USE oldcp
USE upf
IMPLICIT NONE
real(8), PARAMETER :: rmax = 10.0d0
real(8), ALLOCATABLE :: aux(:)
real(8) :: vll
character (len=20):: dft
character (len=2), external :: atom_name
integer :: kkbeta
integer :: l, i, ir, iv
CHARACTER (len=20):: dft
CHARACTER (len=2), EXTERNAL :: atom_name
INTEGER :: kkbeta
INTEGER :: l, i, ir, iv
!
write(generated, '("Generated using unknown code")')
write(date_author,'("Author: unknown Generation date: as well")')
WRITE(generated, '("Generated using unknown code")')
WRITE(date_author,'("Author: unknown Generation date: as well")')
comment = 'Info: automatically converted from old CP90 format'
! reasonable assumption
if (z > 18) then
IF (z > 18) THEN
rel = 1
else
ELSE
rel = 0
end if
ENDIF
rcloc = 0.0d0
nwfs = nbeta_
allocate( els(nwfs), oc(nwfs), epseu(nwfs))
allocate(lchi(nwfs), nns(nwfs) )
allocate(rcut (nwfs), rcutus (nwfs))
do i=1, nwfs
print '("Wavefunction # ",i1,": label, occupancy > ",$)', i
read (5,*) els(i), oc(i)
ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs))
ALLOCATE(lchi(nwfs), nns(nwfs) )
ALLOCATE(rcut (nwfs), rcutus (nwfs))
DO i=1, nwfs
PRINT '("Wavefunction # ",i1,": label, occupancy > ",$)', i
READ (5,*) els(i), oc(i)
nns (i) = 0
lchi(i) = i-1
rcut(i) = 0.0d0
rcutus(i)= 0.0d0
epseu(i) = 0.0d0
end do
ENDDO
psd = atom_name (nint(z))
pseudotype = 'NC'
nlcc = .false.
@ -151,90 +151,90 @@ subroutine convert_oldcp
nbeta = nbeta_
mesh = mesh_
ntwfc = nwfs
allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
do i=1, nwfs
ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
DO i=1, nwfs
lchiw(i) = lchi(i)
ocw(i) = oc(i)
elsw(i) = els(i)
end do
ENDDO
!
if ( exfact.eq.0) then
IF ( exfact==0) THEN
iexch=1; icorr=1; igcx=0; igcc=0 ! Perdew-Zunger
else if ( exfact.eq.1) then
ELSEIF ( exfact==1) THEN
iexch=1; icorr=3; igcx=1; igcc=3 ! Becke-Lee-Yang-Parr
else if ( exfact.eq.2) then
ELSEIF ( exfact==2) THEN
iexch=1; icorr=1; igcx=1; igcc=0 ! Becke88 exchange
else if (exfact.eq.-5.or.exfact.eq.3) then
ELSEIF (exfact==-5.or.exfact==3) THEN
iexch=1; icorr=1; igcx=1; igcc=1 ! Becke88-Perdew 86
else if (exfact.eq.-6.or.exfact.eq.4) then
ELSEIF (exfact==-6.or.exfact==4) THEN
iexch=1; icorr=4; igcx=2; igcc=2 ! Perdew-Wang 91
else if (exfact.eq. 5) then
ELSEIF (exfact== 5) THEN
iexch=1; icorr=4; igcx=3; igcc=4 ! Perdew-Becke-Erkerhof
else
call errore('convert','Wrong xc in pseudopotential',1)
end if
ELSE
CALL errore('convert','Wrong xc in pseudopotential',1)
ENDIF
allocate(rab(mesh))
allocate( r(mesh))
ALLOCATE(rab(mesh))
ALLOCATE( r(mesh))
r = r_
rab = r * log( amesh )
!
! convert analytic to numeric form
!
!
allocate (vloc0(mesh))
ALLOCATE (vloc0(mesh))
! the factor 2 converts from Hartree to Rydberg
vloc0(:) = vnl(:,lloc)*2.d0
if (nbeta > 0) then
IF (nbeta > 0) THEN
allocate(ikk2(nbeta), lll(nbeta))
ALLOCATE(ikk2(nbeta), lll(nbeta))
kkbeta=mesh
do ir = 1,mesh
if ( r(ir) > rmax ) then
DO ir = 1,mesh
IF ( r(ir) > rmax ) THEN
kkbeta=ir
exit
end if
end do
ENDIF
ENDDO
ikk2(:) = kkbeta
allocate(aux(kkbeta))
allocate(betar(mesh,nbeta))
allocate(qfunc(mesh,nbeta,nbeta))
allocate(dion(nbeta,nbeta))
allocate(qqq (nbeta,nbeta))
ALLOCATE(aux(kkbeta))
ALLOCATE(betar(mesh,nbeta))
ALLOCATE(qfunc(mesh,nbeta,nbeta))
ALLOCATE(dion(nbeta,nbeta))
ALLOCATE(qqq (nbeta,nbeta))
qfunc(:,:,:)=0.0d0
dion(:,:) =0.d0
qqq(:,:) =0.d0
iv=0
do i=1,nwfs
DO i=1,nwfs
l=lchi(i)
if (l.ne.lloc) then
IF (l/=lloc) THEN
iv=iv+1
lll(iv)=l
do ir=1,kkbeta
DO ir=1,kkbeta
! the factor 2 converts from Hartree to Rydberg
betar(ir,iv) = 2.d0 * chi_(ir,l+1) * &
( vnl(ir,l) - vnl(ir,lloc) )
aux(ir) = chi_(ir,l+1) * betar(ir,iv)
end do
call simpson(kkbeta,aux,rab,vll)
ENDDO
CALL simpson(kkbeta,aux,rab,vll)
dion(iv,iv) = 1.0d0/vll
end if
enddo
ENDIF
ENDDO
end if
ENDIF
allocate (rho_at(mesh))
ALLOCATE (rho_at(mesh))
rho_at = 0.d0
do i=1,nwfs
DO i=1,nwfs
rho_at(:) = rho_at(:) + ocw(i) * chi_(:,i) ** 2
end do
ENDDO
allocate (chi(mesh,ntwfc))
ALLOCATE (chi(mesh,ntwfc))
chi = chi_
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential successfully converted'
WRITE (6,'(a)') 'Pseudopotential successfully converted'
! ----------------------------------------------------------
return
end subroutine convert_oldcp
RETURN
END SUBROUTINE convert_oldcp

434
upftools/read_upf.f90
View File

@ -5,32 +5,32 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
module pseudo
MODULE pseudo
!
! All variables to be read from the UPF file
! (UPF = unified pseudopotential format)
!
integer ,parameter :: npsx = 6
INTEGER ,PARAMETER :: npsx = 6
! npsx : maximum number of different pseudopotentials
integer, parameter :: lmaxx = 3, nchix = 6, ndm = 2000
INTEGER, PARAMETER :: lmaxx = 3, nchix = 6, ndm = 2000
! lmaxx : maximum non local angular momentum in PP
! nchix : maximum number of atomic wavefunctions per PP
! ndm : maximum number of points in the radial mesh
integer, parameter :: nbrx = 8, lqmax = 5, nqfx = 8
INTEGER, PARAMETER :: nbrx = 8, lqmax = 5, nqfx = 8
! nbrx : maximum number of beta functions
! lqmax : maximum number of angular momentum of Q
! nqfx : maximum number of coefficients in Q smoothing
!
! pp_header
character (len=80):: generated, date_author, comment
character (len=2) :: psd(npsx), pseudotype
character (len=20):: dft(npsx)
integer :: lmax(npsx), mesh(npsx), nbeta(npsx), ntwfc(npsx)
logical :: nlcc(npsx), isus(npsx)
CHARACTER (len=80):: generated, date_author, comment
CHARACTER (len=2) :: psd(npsx), pseudotype
CHARACTER (len=20):: dft(npsx)
INTEGER :: lmax(npsx), mesh(npsx), nbeta(npsx), ntwfc(npsx)
LOGICAL :: nlcc(npsx), isus(npsx)
real(8) :: zp(npsx), ecutrho, ecutwfc, etotps
real(8) :: oc(nchix,npsx)
character(len=2) :: els(nchix,npsx)
integer :: lchi(nchix,npsx)
CHARACTER(len=2) :: els(nchix,npsx)
INTEGER :: lchi(nchix,npsx)
!
! pp_mesh
real(8) :: r(ndm,npsx), rab(ndm,npsx)
@ -43,11 +43,11 @@ module pseudo
! pp_nonlocal
! pp_beta
real(8) :: betar(ndm, nbrx, npsx)
integer :: lll(nbrx,npsx), ikk2(nbrx,npsx)
INTEGER :: lll(nbrx,npsx), ikk2(nbrx,npsx)
! pp_dij
real(8) :: dion(nbrx,nbrx,npsx)
! pp_qij
integer :: nqf(npsx), nqlc(npsx)
INTEGER :: nqf(npsx), nqlc(npsx)
real(8) :: rinner(lqmax,npsx), qqq(nbrx,nbrx,npsx), &
qfunc(ndm,nbrx,nbrx,npsx)
! pp_qfcoef
@ -58,339 +58,339 @@ module pseudo
!
! pp_rhoatom
real(8) :: rho_at(ndm,npsx)
end module pseudo
END MODULE pseudo
!
!---------------------------------------------------------------------
program read_ps
PROGRAM read_ps
!---------------------------------------------------------------------
!
! Read pseudopotentials in the Unified Pseudopotential Format (UPF)
!
implicit none
integer :: is, ios, iunps = 4
character (len=256) :: filein
IMPLICIT NONE
INTEGER :: is, ios, iunps = 4
CHARACTER (len=256) :: filein
!
is = 0
10 print '('' Input PP file # '',i2,'' in UPF format > '',$)', is+1
read (5, '(a)', end = 20, err = 20) filein
open(unit=iunps,file=filein,status='old',form='formatted',iostat=ios)
if (ios.ne.0) stop
10 PRINT '('' Input PP file # '',i2,'' in UPF format > '',$)', is+1
READ (5, '(a)', end = 20, err = 20) filein
OPEN(unit=iunps,file=filein,status='old',form='formatted',iostat=ios)
IF (ios/=0) STOP
is = is + 1
call read_pseudo(is, iunps)
close (unit=iunps)
go to 10
20 stop
end program read_ps
CALL read_pseudo(is, iunps)
CLOSE (unit=iunps)
GOTO 10
20 STOP
END PROGRAM read_ps
!
!---------------------------------------------------------------------
subroutine read_pseudo (is, iunps)
SUBROUTINE read_pseudo (is, iunps)
!---------------------------------------------------------------------
!
! Read pseudopotential in the Unified Pseudopotential Format (UPF)
!
use pseudo
implicit none
USE pseudo
IMPLICIT NONE
!
integer :: is, iunps
INTEGER :: is, iunps
! is : index of this pseudopotential
! iunps: unit connected with pseudopotential file
!
if (is < 0 .or. is > npsx ) call errore ('read_pseudo', 'Wrong is number', 1)
write ( *, * ) " Reading pseudopotential file in UPF format..."
IF (is < 0 .or. is > npsx ) CALL errore ('read_pseudo', 'Wrong is number', 1)
WRITE ( *, * ) " Reading pseudopotential file in UPF format..."
!------->Search for Header
call scan_begin (iunps, "HEADER", .true.)
call read_pseudo_header (is, iunps)
call scan_end (iunps, "HEADER")
CALL scan_begin (iunps, "HEADER", .true.)
CALL read_pseudo_header (is, iunps)
CALL scan_end (iunps, "HEADER")
!-------->Search for mesh information
call scan_begin (iunps, "MESH", .true.)
call read_pseudo_mesh (is, iunps)
call scan_end (iunps, "MESH")
CALL scan_begin (iunps, "MESH", .true.)
CALL read_pseudo_mesh (is, iunps)
CALL scan_end (iunps, "MESH")
!-------->If present, search for nlcc
if (nlcc (is) ) then
call scan_begin (iunps, "NLCC", .true.)
call read_pseudo_nlcc (is, iunps)
call scan_end (iunps, "NLCC")
endif
IF (nlcc (is) ) THEN
CALL scan_begin (iunps, "NLCC", .true.)
CALL read_pseudo_nlcc (is, iunps)
CALL scan_end (iunps, "NLCC")
ENDIF
!-------->Search for Local potential
call scan_begin (iunps, "LOCAL", .true.)
call read_pseudo_local (is, iunps)
call scan_end (iunps, "LOCAL")
CALL scan_begin (iunps, "LOCAL", .true.)
CALL read_pseudo_local (is, iunps)
CALL scan_end (iunps, "LOCAL")
!-------->Search for Nonlocal potential
call scan_begin (iunps, "NONLOCAL", .true.)
call read_pseudo_nl (is, iunps)
call scan_end (iunps, "NONLOCAL")
CALL scan_begin (iunps, "NONLOCAL", .true.)
CALL read_pseudo_nl (is, iunps)
CALL scan_end (iunps, "NONLOCAL")
!-------->Search for atomic wavefunctions
call scan_begin (iunps, "PSWFC", .true.)
call read_pseudo_pswfc (is, iunps)
call scan_end (iunps, "PSWFC")
CALL scan_begin (iunps, "PSWFC", .true.)
CALL read_pseudo_pswfc (is, iunps)
CALL scan_end (iunps, "PSWFC")
!-------->Search for atomic charge
call scan_begin (iunps, "RHOATOM", .true.)
call read_pseudo_rhoatom (is, iunps)
call scan_end (iunps, "RHOATOM")
CALL scan_begin (iunps, "RHOATOM", .true.)
CALL read_pseudo_rhoatom (is, iunps)
CALL scan_end (iunps, "RHOATOM")
!
write ( *, * ) " ...done"
return
end subroutine read_pseudo
WRITE ( *, * ) " ...done"
RETURN
END SUBROUTINE read_pseudo
!---------------------------------------------------------------------
subroutine scan_begin (iunps, string, rew)
SUBROUTINE scan_begin (iunps, string, rew)
!---------------------------------------------------------------------
!
implicit none
IMPLICIT NONE
! Unit of the input file
integer :: iunps
INTEGER :: iunps
! Label to be matched
character (len=*) :: string
logical :: rew
CHARACTER (len=*) :: string
LOGICAL :: rew
! Flag: if .true. rewind the file
character (len=80) :: rstring
CHARACTER (len=80) :: rstring
! String read from file
integer :: ios
logical, external :: matches
INTEGER :: ios
LOGICAL, EXTERNAL :: matches
ios = 0
if (rew) rewind (iunps)
do while (ios.eq.0)
read (iunps, *, iostat = ios, err = 300) rstring
if (matches ("<PP_"//string//">", rstring) ) return
enddo
300 call errore ('scan_begin', 'No '//string//' block', abs (ios) )
IF (rew) REWIND (iunps)
DO WHILE (ios==0)
READ (iunps, *, iostat = ios, err = 300) rstring
IF (matches ("<PP_"//string//">", rstring) ) RETURN
ENDDO
300 CALL errore ('scan_begin', 'No '//string//' block', abs (ios) )
end subroutine scan_begin
END SUBROUTINE scan_begin
!---------------------------------------------------------------------
subroutine scan_end (iunps, string)
SUBROUTINE scan_end (iunps, string)
!---------------------------------------------------------------------
implicit none
IMPLICIT NONE
! Unit of the input file
integer :: iunps
INTEGER :: iunps
! Label to be matched
character (len=*) :: string
CHARACTER (len=*) :: string
! String read from file
character (len=80) :: rstring
integer :: ios
logical, external :: matches
CHARACTER (len=80) :: rstring
INTEGER :: ios
LOGICAL, EXTERNAL :: matches
read (iunps, '(a)', iostat = ios, err = 300) rstring
if (matches ("</PP_"//string//">", rstring) ) return
300 call errore ('scan_end', &
READ (iunps, '(a)', iostat = ios, err = 300) rstring
IF (matches ("</PP_"//string//">", rstring) ) RETURN
300 CALL errore ('scan_end', &
'No '//string//' block end statement, possibly corrupted file', - 1)
end subroutine scan_end
END SUBROUTINE scan_end
!
!---------------------------------------------------------------------
subroutine read_pseudo_header (is, iunps)
SUBROUTINE read_pseudo_header (is, iunps)
!---------------------------------------------------------------------
!
use pseudo
implicit none
USE pseudo
IMPLICIT NONE
!
integer :: is, iunps
INTEGER :: is, iunps
!
integer :: nv, ios, nw
character (len=75) :: dummy
logical, external :: matches
INTEGER :: nv, ios, nw
CHARACTER (len=75) :: dummy
LOGICAL, EXTERNAL :: matches
read (iunps, *, err = 100, iostat = ios) nv, dummy
read (iunps, *, err = 100, iostat = ios) psd (is), dummy
read (iunps, *, err = 100, iostat = ios) pseudotype
if (matches (pseudotype, "US") ) isus (is) = .true.
read (iunps, *, err = 100, iostat = ios) nlcc (is), dummy
read (iunps, '(a20,t24,a)', err = 100, iostat = ios) dft(is), dummy
read (iunps, * ) zp (is), dummy
read (iunps, * ) etotps, dummy
read (iunps, * ) ecutwfc, ecutrho
read (iunps, * ) lmax (is), dummy
read (iunps, *, err = 100, iostat = ios) mesh (is), dummy
read (iunps, *, err = 100, iostat = ios) ntwfc(is), nbeta (is), dummy
read (iunps, '(a)', err = 100, iostat = ios) dummy
do nw = 1, ntwfc(is)
read (iunps, * ) els (nw,is), lchi (nw, is), oc (nw, is)
enddo
return
100 call errore ('read_pseudo_header', 'Reading pseudo file', abs (ios))
end subroutine read_pseudo_header
READ (iunps, *, err = 100, iostat = ios) nv, dummy
READ (iunps, *, err = 100, iostat = ios) psd (is), dummy
READ (iunps, *, err = 100, iostat = ios) pseudotype
IF (matches (pseudotype, "US") ) isus (is) = .true.
READ (iunps, *, err = 100, iostat = ios) nlcc (is), dummy
READ (iunps, '(a20,t24,a)', err = 100, iostat = ios) dft(is), dummy
READ (iunps, * ) zp (is), dummy
READ (iunps, * ) etotps, dummy
READ (iunps, * ) ecutwfc, ecutrho
READ (iunps, * ) lmax (is), dummy
READ (iunps, *, err = 100, iostat = ios) mesh (is), dummy
READ (iunps, *, err = 100, iostat = ios) ntwfc(is), nbeta (is), dummy
READ (iunps, '(a)', err = 100, iostat = ios) dummy
DO nw = 1, ntwfc(is)
READ (iunps, * ) els (nw,is), lchi (nw, is), oc (nw, is)
ENDDO
RETURN
100 CALL errore ('read_pseudo_header', 'Reading pseudo file', abs (ios))
END SUBROUTINE read_pseudo_header
!
!---------------------------------------------------------------------
subroutine read_pseudo_local (is, iunps)
SUBROUTINE read_pseudo_local (is, iunps)
!---------------------------------------------------------------------
!
use pseudo
implicit none
USE pseudo
IMPLICIT NONE
!
integer :: is, iunps
INTEGER :: is, iunps
!
integer :: ir, ios
INTEGER :: ir, ios
!
read (iunps, *, err=100, iostat=ios) (vloc0(ir,is) , ir=1,mesh(is))
READ (iunps, *, err=100, iostat=ios) (vloc0(ir,is) , ir=1,mesh(is))
100 call errore ('read_pseudo_local','Reading pseudo file', abs(ios) )
100 CALL errore ('read_pseudo_local','Reading pseudo file', abs(ios) )
return
end subroutine read_pseudo_local
RETURN
END SUBROUTINE read_pseudo_local
!
!---------------------------------------------------------------------
subroutine read_pseudo_mesh (is, iunps)
SUBROUTINE read_pseudo_mesh (is, iunps)
!---------------------------------------------------------------------
!
use pseudo
implicit none
USE pseudo
IMPLICIT NONE
!
integer :: is, iunps
INTEGER :: is, iunps
!
integer :: ir, ios
INTEGER :: ir, ios
!
call scan_begin (iunps, "R", .false.)
read (iunps, *, err = 100, iostat = ios) (r(ir,is), ir=1,mesh(is) )
call scan_end (iunps, "R")
call scan_begin (iunps, "RAB", .false.)
read (iunps, *, err = 100, iostat = ios) (rab(ir,is), ir=1,mesh(is) )
call scan_end (iunps, "RAB")
CALL scan_begin (iunps, "R", .false.)
READ (iunps, *, err = 100, iostat = ios) (r(ir,is), ir=1,mesh(is) )
CALL scan_end (iunps, "R")
CALL scan_begin (iunps, "RAB", .false.)
READ (iunps, *, err = 100, iostat = ios) (rab(ir,is), ir=1,mesh(is) )
CALL scan_end (iunps, "RAB")
return
RETURN
100 call errore ('read_pseudo_mesh', 'Reading pseudo file', abs (ios) )
end subroutine read_pseudo_mesh
100 CALL errore ('read_pseudo_mesh', 'Reading pseudo file', abs (ios) )
END SUBROUTINE read_pseudo_mesh
!
!---------------------------------------------------------------------
subroutine read_pseudo_nl (is, iunps)
SUBROUTINE read_pseudo_nl (is, iunps)
!---------------------------------------------------------------------
!
use pseudo
implicit none
USE pseudo
IMPLICIT NONE
!
integer :: is, iunps
INTEGER :: is, iunps
!
integer :: nb, mb, n, ir, nd, ios, idum, ldum, icon, lp, i
INTEGER :: nb, mb, n, ir, nd, ios, idum, ldum, icon, lp, i
! counters
character (len=75) :: dummy
CHARACTER (len=75) :: dummy
!
do nb = 1, nbeta (is)
call scan_begin (iunps, "BETA", .false.)
read (iunps, *, err = 100, iostat = ios) idum, lll(nb,is), dummy
read (iunps, '(i6)', err = 100, iostat = ios) ikk2(nb,is)
read (iunps, *, err = 100, iostat = ios) &
DO nb = 1, nbeta (is)
CALL scan_begin (iunps, "BETA", .false.)
READ (iunps, *, err = 100, iostat = ios) idum, lll(nb,is), dummy
READ (iunps, '(i6)', err = 100, iostat = ios) ikk2(nb,is)
READ (iunps, *, err = 100, iostat = ios) &
(betar(ir,nb,is), ir=1,ikk2(nb,is))
do ir = ikk2(nb,is) + 1, mesh (is)
DO ir = ikk2(nb,is) + 1, mesh (is)
betar (ir, nb, is) = 0.d0
enddo
call scan_end (iunps, "BETA")
enddo
ENDDO
CALL scan_end (iunps, "BETA")
ENDDO
call scan_begin (iunps, "DIJ", .false.)
read (iunps, *, err = 100, iostat = ios) nd, dummy
CALL scan_begin (iunps, "DIJ", .false.)
READ (iunps, *, err = 100, iostat = ios) nd, dummy
dion (:,:,is) = 0.d0
do icon = 1, nd
read (iunps, *, err = 100, iostat = ios) nb, mb, dion(nb,mb,is)
DO icon = 1, nd
READ (iunps, *, err = 100, iostat = ios) nb, mb, dion(nb,mb,is)
dion (mb,nb,is) = dion (nb,mb,is)
enddo
call scan_end (iunps, "DIJ")
ENDDO
CALL scan_end (iunps, "DIJ")
if (isus (is) ) then
call scan_begin (iunps, "QIJ", .false.)
read (iunps, *, err = 100, iostat = ios) nqf(is)
IF (isus (is) ) THEN
CALL scan_begin (iunps, "QIJ", .false.)
READ (iunps, *, err = 100, iostat = ios) nqf(is)
nqlc (is)= 2 * lmax (is) + 1
if (nqlc(is).gt.lqmax .or. nqlc(is).lt.0) &
call errore (' read_pseudo_nl', 'Wrong nqlc', nqlc (is) )
if (nqf(is).ne.0) then
call scan_begin (iunps, "RINNER", .false.)
read (iunps,*,err=100,iostat=ios) &
IF (nqlc(is)>lqmax .or. nqlc(is)<0) &
CALL errore (' read_pseudo_nl', 'Wrong nqlc', nqlc (is) )
IF (nqf(is)/=0) THEN
CALL scan_begin (iunps, "RINNER", .false.)
READ (iunps,*,err=100,iostat=ios) &
(idum,rinner(i,is),i=1,nqlc(is))
call scan_end (iunps, "RINNER")
end if
do nb = 1, nbeta(is)
do mb = nb, nbeta(is)
CALL scan_end (iunps, "RINNER")
ENDIF
DO nb = 1, nbeta(is)
DO mb = nb, nbeta(is)
read (iunps,*,err=100,iostat=ios) idum, idum, ldum, dummy
READ (iunps,*,err=100,iostat=ios) idum, idum, ldum, dummy
!" i j (l)"
if (ldum.ne.lll(mb,is) ) call errore ('read_pseudo_nl', &
IF (ldum/=lll(mb,is) ) CALL errore ('read_pseudo_nl', &
'inconsistent angular momentum for Q_ij', 1)
read (iunps,*,err=100,iostat=ios) qqq(nb,mb,is), dummy
READ (iunps,*,err=100,iostat=ios) qqq(nb,mb,is), dummy
! "Q_int"
qqq(mb,nb,is) = qqq(nb,mb,is)
read (iunps,*,err=100,iostat=ios) &
READ (iunps,*,err=100,iostat=ios) &
(qfunc(n,nb,mb,is), n=1,mesh(is))
do n = 0, mesh (is)
DO n = 0, mesh (is)
qfunc(n,mb,nb,is) = qfunc(n,nb,mb,is)
enddo
ENDDO
if (nqf(is).gt.0) then
call scan_begin (iunps, "QFCOEF", .false.)
read (iunps,*,err=100,iostat=ios) &
IF (nqf(is)>0) THEN
CALL scan_begin (iunps, "QFCOEF", .false.)
READ (iunps,*,err=100,iostat=ios) &
((qfcoef(i,lp,nb,mb,is),i=1,nqf(is)),lp=1,nqlc(is))
call scan_end (iunps, "QFCOEF")
end if
CALL scan_end (iunps, "QFCOEF")
ENDIF
enddo
enddo
call scan_end (iunps, "QIJ")
else
ENDDO
ENDDO
CALL scan_end (iunps, "QIJ")
ELSE
qqq (:,:,is) = 0.d0
qfunc(:,:,:,is) =0.d0
endif
ENDIF
100 call errore ('read_pseudo_nl', 'Reading pseudo file', abs (ios) )
return
end subroutine read_pseudo_nl
100 CALL errore ('read_pseudo_nl', 'Reading pseudo file', abs (ios) )
RETURN
END SUBROUTINE read_pseudo_nl
!
!---------------------------------------------------------------------
subroutine read_pseudo_nlcc (is, iunps)
SUBROUTINE read_pseudo_nlcc (is, iunps)
!---------------------------------------------------------------------
!
use pseudo
implicit none
USE pseudo
IMPLICIT NONE
!
integer :: is, iunps
INTEGER :: is, iunps
!
integer :: ir, ios
INTEGER :: ir, ios
read (iunps, *, err = 100, iostat = ios) (rho_atc(ir,is), ir=1,mesh(is) )
READ (iunps, *, err = 100, iostat = ios) (rho_atc(ir,is), ir=1,mesh(is) )
!
100 call errore ('read_pseudo_nlcc', 'Reading pseudo file', abs (ios) )
return
end subroutine read_pseudo_nlcc
100 CALL errore ('read_pseudo_nlcc', 'Reading pseudo file', abs (ios) )
RETURN
END SUBROUTINE read_pseudo_nlcc
!
!---------------------------------------------------------------------
subroutine read_pseudo_pswfc (is, iunps)
SUBROUTINE read_pseudo_pswfc (is, iunps)
!---------------------------------------------------------------------
!
use pseudo
implicit none
USE pseudo
IMPLICIT NONE
!
integer :: is, iunps
INTEGER :: is, iunps
!
character (len=75) :: dummy
integer :: nb, ir, ios
CHARACTER (len=75) :: dummy
INTEGER :: nb, ir, ios
!
do nb = 1, ntwfc(is)
read (iunps,*,err=100,iostat=ios) dummy !Wavefunction labels
read (iunps,*,err=100,iostat=ios) (chi(ir,nb,is), ir=1,mesh(is))
enddo
100 call errore ('read_pseudo_pswfc', 'Reading pseudo file', abs(ios))
return
DO nb = 1, ntwfc(is)
READ (iunps,*,err=100,iostat=ios) dummy !Wavefunction labels
READ (iunps,*,err=100,iostat=ios) (chi(ir,nb,is), ir=1,mesh(is))
ENDDO
100 CALL errore ('read_pseudo_pswfc', 'Reading pseudo file', abs(ios))
RETURN
end subroutine read_pseudo_pswfc
END SUBROUTINE read_pseudo_pswfc
!
!---------------------------------------------------------------------
subroutine read_pseudo_rhoatom (is, iunps)
SUBROUTINE read_pseudo_rhoatom (is, iunps)
!---------------------------------------------------------------------
!
use pseudo
implicit none
USE pseudo
IMPLICIT NONE
!
integer :: is, iunps
INTEGER :: is, iunps
!
integer :: ir, ios
INTEGER :: ir, ios
read (iunps,*,err=100,iostat=ios) (rho_at(ir,is), ir=1,mesh(is))
return
READ (iunps,*,err=100,iostat=ios) (rho_at(ir,is), ir=1,mesh(is))
RETURN
100 call errore ('read_pseudo_rhoatom','Reading pseudo file',abs(ios))
100 CALL errore ('read_pseudo_rhoatom','Reading pseudo file',abs(ios))
end subroutine read_pseudo_rhoatom
END SUBROUTINE read_pseudo_rhoatom

20
upftools/read_upf_tofile.f90
View File

@ -21,10 +21,10 @@ PROGRAM read_upf_tofile
! PWSCF modules
!
!
USE constants, only : fpi
USE constants, ONLY : fpi
USE pseudo_types
USE upf_module
USE radial_grids, only : radial_grid_type, nullify_radial_grid
USE radial_grids, ONLY : radial_grid_type, nullify_radial_grid
!
IMPLICIT NONE
!
@ -53,39 +53,39 @@ PROGRAM read_upf_tofile
CALL read_upf(upf, grid, ierr, unit=iunps)
!
IF (ierr .NE. 0) &
CALL errore('read_upf_tofile','reading pseudo upf', ABS(ierr))
IF (ierr /= 0) &
CALL errore('read_upf_tofile','reading pseudo upf', abs(ierr))
!
CLOSE(iunps)
!
OPEN(UNIT=iunps,FILE='filewfc',STATUS='unknown',FORM='formatted', &
ERR=200, IOSTAT=ios)
200 CALL errore('read_upf_tofile','open error on file filewfc',ABS(ios))
200 CALL errore('read_upf_tofile','open error on file filewfc',abs(ios))
DO n=1,upf%mesh
WRITE(iunps,'(30f12.6)') upf%r(n), (upf%chi(n,j), j=1,upf%nwfc)
END DO
ENDDO
CLOSE(iunps)
OPEN(UNIT=iunps,FILE='filebeta',STATUS='unknown',FORM='formatted', &
ERR=300, IOSTAT=ios)
300 CALL errore('read_upf_tofile','open error on file filebeta',ABS(ios))
300 CALL errore('read_upf_tofile','open error on file filebeta',abs(ios))
DO n=1,upf%mesh
WRITE(iunps,'(30f12.6)') upf%r(n), (upf%beta(n,j), j=1,upf%nbeta)
END DO
ENDDO
CLOSE(iunps)
OPEN(UNIT=iunps,FILE='filepot',STATUS='unknown',FORM='formatted', &
ERR=400, IOSTAT=ios)
400 CALL errore('read_upf_tofile','open error on file filepot',ABS(ios))
400 CALL errore('read_upf_tofile','open error on file filepot',abs(ios))
DO n=1,upf%mesh
WRITE(iunps,'(4f12.6)') upf%r(n), upf%vloc(n), &
upf%rho_at(n), upf%rho_atc(n)*fpi*upf%r(n)**2
END DO
ENDDO
CLOSE(iunps)

280
upftools/rrkj2upf.f90
View File

@ -7,178 +7,178 @@
!
!
!---------------------------------------------------------------------
program rrkj2upf
PROGRAM rrkj2upf
!---------------------------------------------------------------------
!
! Convert a pseudopotential written in "rrkj3" format
! (Rabe-Rappe-Kaxiras-Joannopoulos with 3 Bessel functions)
! to unified pseudopotential format
!
implicit none
character(len=256) filein, fileout
IMPLICIT NONE
CHARACTER(len=256) filein, fileout
!
!
call get_file ( filein )
open (unit = 1, file = filein, status = 'old', form = 'formatted')
call read_rrkj(1)
close (1)
CALL get_file ( filein )
OPEN (unit = 1, file = filein, status = 'old', form = 'formatted')
CALL read_rrkj(1)
CLOSE (1)
! convert variables read from rrkj3 format into those needed
! by the upf format - add missing quantities
call convert_rrkj
CALL convert_rrkj
fileout=trim(filein)//'.UPF'
print '(''Output PP file in UPF format : '',a)', fileout
PRINT '(''Output PP file in UPF format : '',a)', fileout
open(unit=2,file=fileout,status='unknown',form='formatted')
call write_upf(2)
close (unit=2)
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
CALL write_upf(2)
CLOSE (unit=2)
stop
20 write (6,'("rrkj2upf: error reading pseudopotential file name")')
stop
STOP
20 WRITE (6,'("rrkj2upf: error reading pseudopotential file name")')
STOP
end program rrkj2upf
END PROGRAM rrkj2upf
module rrkj3
MODULE rrkj3
!
! All variables read from RRKJ3 file format
!
! trailing underscore means that a variable with the same name
! is used in module 'upf' containing variables to be written
!
character(len=75):: titleps
character (len=2), allocatable :: els_(:)
integer :: pseudotype_, iexch_, icorr_, igcx_, igcc_, mesh_, &
CHARACTER(len=75):: titleps
CHARACTER (len=2), ALLOCATABLE :: els_(:)
INTEGER :: pseudotype_, iexch_, icorr_, igcx_, igcc_, mesh_, &
nwfs_, nbeta_, lmax_
logical :: rel_, nlcc_
LOGICAL :: rel_, nlcc_
real (8) :: zp_, etotps_, xmin, rmax, zmesh, dx, rcloc_
integer, allocatable:: lchi_(:), nns_(:), ikk2_(:)
real (8), allocatable :: rcut_(:), rcutus_(:), oc_(:), &
INTEGER, ALLOCATABLE:: lchi_(:), nns_(:), ikk2_(:)
real (8), ALLOCATABLE :: rcut_(:), rcutus_(:), oc_(:), &
beta(:,:), dion_(:,:), qqq_(:,:), ddd(:,:), qfunc_(:,:,:), &
rho_atc_(:), rho_at_(:), chi_(:,:), vloc_(:)
end module rrkj3
END MODULE rrkj3
!
! ----------------------------------------------------------
subroutine read_rrkj(iunps)
SUBROUTINE read_rrkj(iunps)
! ----------------------------------------------------------
!
use rrkj3
implicit none
integer :: iunps
integer :: nb, mb, n, ir, ios
USE rrkj3
IMPLICIT NONE
INTEGER :: iunps
INTEGER :: nb, mb, n, ir, ios
!--- > Start the header reading
read (iunps, '(a75)', err = 100) titleps
read (iunps, *, err = 100) pseudotype_
read (iunps, *, err = 100) rel_, nlcc_
read (iunps, *, err=100) iexch_, icorr_, igcx_, igcc_
read (iunps, '(2e17.11,i5)') zp_, etotps_, lmax_
read (iunps, '(4e17.11,i5)', err=100) xmin, rmax, zmesh, dx, mesh_
read (iunps, *, err=100) nwfs_, nbeta_
READ (iunps, '(a75)', err = 100) titleps
READ (iunps, *, err = 100) pseudotype_
READ (iunps, *, err = 100) rel_, nlcc_
READ (iunps, *, err=100) iexch_, icorr_, igcx_, igcc_
READ (iunps, '(2e17.11,i5)') zp_, etotps_, lmax_
READ (iunps, '(4e17.11,i5)', err=100) xmin, rmax, zmesh, dx, mesh_
READ (iunps, *, err=100) nwfs_, nbeta_
allocate(rcut_(nwfs_), rcutus_(nwfs_))
read (iunps, *, err=100) (rcut_(nb), nb=1,nwfs_)
read (iunps, *, err=100) (rcutus_(nb), nb=1,nwfs_)
ALLOCATE(rcut_(nwfs_), rcutus_(nwfs_))
READ (iunps, *, err=100) (rcut_(nb), nb=1,nwfs_)
READ (iunps, *, err=100) (rcutus_(nb), nb=1,nwfs_)
allocate(els_(nwfs_), nns_(nwfs_), lchi_(nwfs_), oc_(nwfs_))
do nb = 1, nwfs_
read (iunps, '(a2,2i3,f6.2)', err = 100) els_(nb), &
ALLOCATE(els_(nwfs_), nns_(nwfs_), lchi_(nwfs_), oc_(nwfs_))
DO nb = 1, nwfs_
READ (iunps, '(a2,2i3,f6.2)', err = 100) els_(nb), &
nns_(nb), lchi_(nb) , oc_(nb)
enddo
ENDDO
allocate(ikk2_(nbeta_))
allocate(beta( mesh_,nbeta_))
allocate(dion_(nbeta_,nbeta_))
allocate(ddd (nbeta_,nbeta_))
allocate(qqq_(nbeta_,nbeta_))
allocate(qfunc_(mesh_,nbeta_,nbeta_))
ALLOCATE(ikk2_(nbeta_))
ALLOCATE(beta( mesh_,nbeta_))
ALLOCATE(dion_(nbeta_,nbeta_))
ALLOCATE(ddd (nbeta_,nbeta_))
ALLOCATE(qqq_(nbeta_,nbeta_))
ALLOCATE(qfunc_(mesh_,nbeta_,nbeta_))
do nb = 1, nbeta_
read (iunps, *, err = 100) ikk2_(nb)
read (iunps, *, err = 100) (beta (ir, nb) , ir = 1,ikk2_(nb) )
do ir = ikk2_(nb) + 1, mesh_
DO nb = 1, nbeta_
READ (iunps, *, err = 100) ikk2_(nb)
READ (iunps, *, err = 100) (beta (ir, nb) , ir = 1,ikk2_(nb) )
DO ir = ikk2_(nb) + 1, mesh_
beta (ir, nb) = 0.d0
enddo
do mb = 1, nb
read (iunps, *, err = 100) dion_(nb, mb)
ENDDO
DO mb = 1, nb
READ (iunps, *, err = 100) dion_(nb, mb)
dion_(mb, nb) = dion_(nb, mb)
if (pseudotype_.eq.3) then
read (iunps, *, err = 100) qqq_(nb, mb)
IF (pseudotype_==3) THEN
READ (iunps, *, err = 100) qqq_(nb, mb)
qqq_(mb, nb) = qqq_(nb, mb)
read (iunps, *, err = 100) (qfunc_(n,nb, mb), n = 1, mesh_)
do n = 1, mesh_
READ (iunps, *, err = 100) (qfunc_(n,nb, mb), n = 1, mesh_)
DO n = 1, mesh_
qfunc_(n, mb, nb) = qfunc_(n, nb, mb)
enddo
else
ENDDO
ELSE
qqq_(nb, mb) = 0.d0
qqq_(mb, nb) = 0.d0
do n = 1, mesh_
DO n = 1, mesh_
qfunc_(n, nb, mb) = 0.d0
qfunc_(n, mb, nb) = 0.d0
enddo
endif
enddo
enddo
ENDDO
ENDIF
ENDDO
ENDDO
!
! read the local potential
!
allocate(vloc_(mesh_))
read (iunps, *, err = 100) rcloc_, (vloc_(ir ) , ir = 1, mesh_ )
ALLOCATE(vloc_(mesh_))
READ (iunps, *, err = 100) rcloc_, (vloc_(ir ) , ir = 1, mesh_ )
!
! read the atomic charge
!
allocate(rho_at_(mesh_))
read (iunps, *, err=100) (rho_at_(ir), ir=1,mesh_)
ALLOCATE(rho_at_(mesh_))
READ (iunps, *, err=100) (rho_at_(ir), ir=1,mesh_)
!
! if present read the core charge
!
allocate(rho_atc_(mesh_))
if (nlcc_) then
read (iunps, *, err=100) (rho_atc_(ir), ir=1, mesh_)
endif
ALLOCATE(rho_atc_(mesh_))
IF (nlcc_) THEN
READ (iunps, *, err=100) (rho_atc_(ir), ir=1, mesh_)
ENDIF
!
! read the pseudo wavefunctions of the atom
!
allocate(chi_(mesh_,nwfs_))
read (iunps, *, err=100) ( (chi_(ir,nb), ir = 1,mesh_) , nb = 1, nwfs_)
ALLOCATE(chi_(mesh_,nwfs_))
READ (iunps, *, err=100) ( (chi_(ir,nb), ir = 1,mesh_) , nb = 1, nwfs_)
!
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential successfully read'
WRITE (6,'(a)') 'Pseudopotential successfully read'
! ----------------------------------------------------------
!
return
100 write (6,'("read_rrkj: error reading pseudopotential file")')
stop
RETURN
100 WRITE (6,'("read_rrkj: error reading pseudopotential file")')
STOP
end subroutine read_rrkj
END SUBROUTINE read_rrkj
subroutine convert_rrkj
SUBROUTINE convert_rrkj
! ----------------------------------------------------------
!
use rrkj3
use upf
use constants, only : fpi
implicit none
integer i, n
USE rrkj3
USE upf
USE constants, ONLY : fpi
IMPLICIT NONE
INTEGER i, n
real(8) :: x
write(generated, '("Generated using Andrea Dal Corso code (rrkj3)")')
write(date_author,'("Author: Andrea Dal Corso Generation date: unknown")')
WRITE(generated, '("Generated using Andrea Dal Corso code (rrkj3)")')
WRITE(date_author,'("Author: Andrea Dal Corso Generation date: unknown")')
comment = 'Info:'//titleps
if (rel_) then
IF (rel_) THEN
rel = 1
else
ELSE
rel = 0
end if
ENDIF
rcloc = rcloc_
nwfs = nwfs_
allocate( els(nwfs), oc(nwfs), epseu(nwfs))
allocate(lchi(nwfs), nns(nwfs) )
allocate(rcut (nwfs), rcutus (nwfs))
do i=1, nwfs
ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs))
ALLOCATE(lchi(nwfs), nns(nwfs) )
ALLOCATE(rcut (nwfs), rcutus (nwfs))
DO i=1, nwfs
nns (i) = nns_(i)
lchi(i) = lchi_(i)
rcut(i) = rcut_(i)
@ -186,15 +186,15 @@ subroutine convert_rrkj
oc (i) = oc_(i)
els(i) = els_(i)
epseu(i) = 0.0d0
end do
deallocate (els_, oc_, rcutus_, rcut_, nns_)
ENDDO
DEALLOCATE (els_, oc_, rcutus_, rcut_, nns_)
psd = titleps (7:8)
if (pseudotype_.eq.3) then
IF (pseudotype_==3) THEN
pseudotype = 'US'
else
ELSE
pseudotype = 'NC'
endif
ENDIF
nlcc = nlcc_
zp = zp_
etotps = etotps_
@ -204,85 +204,85 @@ subroutine convert_rrkj
mesh = mesh_
nbeta = nbeta_
ntwfc = 0
do i=1, nwfs
if (oc(i) .gt. 1.0d-12) ntwfc = ntwfc + 1
end do
allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
DO i=1, nwfs
IF (oc(i) > 1.0d-12) ntwfc = ntwfc + 1
ENDDO
ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
n = 0
do i=1, nwfs
if (oc(i) .gt. 1.0d-12) then
DO i=1, nwfs
IF (oc(i) > 1.0d-12) THEN
n = n + 1
elsw(n) = els(i)
ocw (n) = oc (i)
lchiw(n)=lchi(i)
end if
end do
ENDIF
ENDDO
iexch = iexch_
icorr = icorr_
igcx = igcx_
igcc = igcc_
allocate(rab(mesh))
allocate( r(mesh))
ALLOCATE(rab(mesh))
ALLOCATE( r(mesh))
! define logarithmic mesh
do i = 1, mesh
x = xmin + DBLE(i-1) * dx
DO i = 1, mesh
x = xmin + dble(i-1) * dx
r (i) = exp(x) / zmesh
rab(i) = dx * r(i)
end do
ENDDO
allocate (rho_atc(mesh))
ALLOCATE (rho_atc(mesh))
! rrkj rho_core(r) = 4pi*r^2*rho_core(r) UPF
rho_atc (:) = rho_atc_(:) / fpi / r(:)**2
deallocate (rho_atc_)
DEALLOCATE (rho_atc_)
allocate (vloc0(mesh))
ALLOCATE (vloc0(mesh))
vloc0 = vloc_
deallocate (vloc_)
DEALLOCATE (vloc_)
allocate(ikk2(nbeta), lll(nbeta))
ALLOCATE(ikk2(nbeta), lll(nbeta))
ikk2 = ikk2_
lll = lchi_
deallocate (ikk2_, lchi_)
DEALLOCATE (ikk2_, lchi_)
! kkbeta = 0
! do nb=1,nbeta
! kkbeta = max (kkbeta , ikk2(nb) )
! end do
allocate(betar(mesh,nbeta))
ALLOCATE(betar(mesh,nbeta))
betar = 0.0d0
do i=1, nbeta
DO i=1, nbeta
betar(1:ikk2(i),i) = beta(1:ikk2(i),i)
end do
deallocate (beta)
ENDDO
DEALLOCATE (beta)
allocate(dion(nbeta,nbeta))
ALLOCATE(dion(nbeta,nbeta))
dion = dion_
deallocate (dion_)
DEALLOCATE (dion_)
allocate(qqq(nbeta,nbeta))
ALLOCATE(qqq(nbeta,nbeta))
qqq = qqq_
deallocate (qqq_)
DEALLOCATE (qqq_)
allocate(qfunc(mesh,nbeta,nbeta))
ALLOCATE(qfunc(mesh,nbeta,nbeta))
qfunc = qfunc_
nqf = 0
nqlc= 0
allocate (rho_at(mesh))
ALLOCATE (rho_at(mesh))
rho_at = rho_at_
deallocate (rho_at_)
DEALLOCATE (rho_at_)
allocate (chi(mesh,ntwfc))
ALLOCATE (chi(mesh,ntwfc))
n = 0
do i=1, nwfs
if (oc(i) .gt. 1.0d-12) then
DO i=1, nwfs
IF (oc(i) > 1.0d-12) THEN
n = n + 1
chi(:,n) = chi_(:,i)
end if
end do
deallocate (chi_)
ENDIF
ENDDO
DEALLOCATE (chi_)
return
end subroutine convert_rrkj
RETURN
END SUBROUTINE convert_rrkj

32
upftools/uspp2upf.f90
View File

@ -7,35 +7,35 @@
!
!
!---------------------------------------------------------------------
program uspp2upf
PROGRAM uspp2upf
!---------------------------------------------------------------------
!
! Convert a pseudopotential written in Vanderbilt format
! (unformatted) to unified pseudopotential format
!
implicit none
character(len=256) filein, fileout
IMPLICIT NONE
CHARACTER(len=256) filein, fileout
!
!
call get_file ( filein )
open(unit=1,file=filein,status='old',form='unformatted')
call read_uspp(1)
close (unit=1)
CALL get_file ( filein )
OPEN(unit=1,file=filein,status='old',form='unformatted')
CALL read_uspp(1)
CLOSE (unit=1)
! convert variables read from Vanderbilt format into those needed
! by the upf format - add missing quantities
call convert_uspp
CALL convert_uspp
fileout=trim(filein)//'.UPF'
print '(''Output PP file in UPF format : '',a)', fileout
PRINT '(''Output PP file in UPF format : '',a)', fileout
open(unit=2,file=fileout,status='unknown',form='formatted')
call write_upf(2)
close (unit=2)
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
CALL write_upf(2)
CLOSE (unit=2)
stop
20 write (6,'("uspp2upf: error reading pseudopotential file name")')
stop
end program uspp2upf
STOP
20 WRITE (6,'("uspp2upf: error reading pseudopotential file name")')
STOP
END PROGRAM uspp2upf

406
upftools/vanderbilt.f90
View File

@ -6,273 +6,273 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
module Vanderbilt
MODULE Vanderbilt
!
! All variables read from Vanderbilt's file format
!
! trailing underscore means that a variable with the same name
! is used in module 'upf' containing variables to be written
!
integer :: nvalps, nang, nbeta_, kkbeta, nchi, ifpcor, keyps, &
INTEGER :: nvalps, nang, nbeta_, kkbeta, nchi, ifpcor, keyps, &
mesh_, iver(3), idmy(3), nnlz, ifqopt, nqf_, irel, npf, &
nlc, lloc
real(8) :: z_, zp_, exfact, etot, eloc, rcloc_, rpcor, &
qtryc, ptryc, rinner1_
real(8), allocatable:: wwnlps(:), eeps(:), rinner_(:), rc(:), &
real(8), ALLOCATABLE:: wwnlps(:), eeps(:), rinner_(:), rc(:), &
beta(:,:), ddd0(:,:), ddd(:,:), qqq_(:,:), eee(:), rho_atc_(:), &
r_(:), rab_(:), rho_at_(:), qfunc_(:,:,:), vloc(:), vloc_(:), &
wf(:,:), qfcoef_(:,:,:,:)
integer, allocatable :: lll_(:), nnlzps(:), iptype(:)
Character(len=20):: title
end module Vanderbilt
INTEGER, ALLOCATABLE :: lll_(:), nnlzps(:), iptype(:)
CHARACTER(len=20):: title
END MODULE Vanderbilt
!
! ----------------------------------------------------------
subroutine read_uspp(iunit)
SUBROUTINE read_uspp(iunit)
! ----------------------------------------------------------
!
use Vanderbilt
implicit none
integer :: iunit
USE Vanderbilt
IMPLICIT NONE
INTEGER :: iunit
!
integer :: i, j, k, lp
INTEGER :: i, j, k, lp
real(8) :: rinner1
!
!
read (iunit) (iver(i),i=1,3),(idmy(i),i=1,3)
read (iunit) title, z_, zp_, exfact, nvalps, mesh_, etot
READ (iunit) (iver(i),i=1,3),(idmy(i),i=1,3)
READ (iunit) title, z_, zp_, exfact, nvalps, mesh_, etot
allocate(nnlzps(nvalps), wwnlps(nvalps), eeps(nvalps))
read (iunit) (nnlzps(i),wwnlps(i),eeps(i),i=1,nvalps)
ALLOCATE(nnlzps(nvalps), wwnlps(nvalps), eeps(nvalps))
READ (iunit) (nnlzps(i),wwnlps(i),eeps(i),i=1,nvalps)
read (iunit) keyps, ifpcor, rinner1
READ (iunit) keyps, ifpcor, rinner1
if ( iver(1) .eq. 1 ) then
IF ( iver(1) == 1 ) THEN
nang = nvalps
nqf_ = 3
nlc = 5
elseif ( iver(1) .eq. 2 ) then
ELSEIF ( iver(1) == 2 ) THEN
nang = nvalps
nqf_ = 3
nlc = 2 * nvalps - 1
else if ( iver(1) .ge. 3 ) then
read (iunit) nang, lloc, eloc, ifqopt, nqf_, qtryc
ELSEIF ( iver(1) >= 3 ) THEN
READ (iunit) nang, lloc, eloc, ifqopt, nqf_, qtryc
nlc = 2 * nang - 1
endif
ENDIF
allocate(rinner_(2*nang-1))
ALLOCATE(rinner_(2*nang-1))
rinner_(1) = rinner1
rinner1_ = rinner1
if (10*iver(1)+iver(2).ge.51) &
read (iunit) (rinner_(i),i=1,nang*2-1)
IF (10*iver(1)+iver(2)>=51) &
READ (iunit) (rinner_(i),i=1,nang*2-1)
if ( iver(1) .ge. 4 ) then
read (iunit) irel
else
IF ( iver(1) >= 4 ) THEN
READ (iunit) irel
ELSE
irel = 0
end if
ENDIF
allocate(rc(nang))
read (iunit) (rc(i),i=1,nang)
ALLOCATE(rc(nang))
READ (iunit) (rc(i),i=1,nang)
read (iunit) nbeta_,kkbeta
READ (iunit) nbeta_,kkbeta
!
allocate(beta(kkbeta,nbeta_))
allocate(qfunc_(kkbeta,nbeta_,nbeta_))
allocate(ddd0(nbeta_,nbeta_))
allocate(ddd (nbeta_,nbeta_))
allocate(qqq_(nbeta_,nbeta_))
allocate(lll_(nbeta_))
allocate(eee(nbeta_))
allocate(qfcoef_(nqf_,nlc,nbeta_,nbeta_))
ALLOCATE(beta(kkbeta,nbeta_))
ALLOCATE(qfunc_(kkbeta,nbeta_,nbeta_))
ALLOCATE(ddd0(nbeta_,nbeta_))
ALLOCATE(ddd (nbeta_,nbeta_))
ALLOCATE(qqq_(nbeta_,nbeta_))
ALLOCATE(lll_(nbeta_))
ALLOCATE(eee(nbeta_))
ALLOCATE(qfcoef_(nqf_,nlc,nbeta_,nbeta_))
!
do j=1,nbeta_
read (iunit) lll_(j),eee(j),(beta(i,j),i=1,kkbeta)
do k=j,nbeta_
read (iunit) ddd0(j,k),ddd(j,k),qqq_(j,k), &
DO j=1,nbeta_
READ (iunit) lll_(j),eee(j),(beta(i,j),i=1,kkbeta)
DO k=j,nbeta_
READ (iunit) ddd0(j,k),ddd(j,k),qqq_(j,k), &
(qfunc_(i,j,k),i=1,kkbeta), &
((qfcoef_(i,lp,j,k),i=1,nqf_),lp=1,2*nang-1)
end do
end do
ENDDO
ENDDO
!
allocate(iptype(nbeta_))
if (10*iver(1)+iver(2).ge.72) &
read (iunit) (iptype(j),j=1,nbeta_),npf,ptryc
ALLOCATE(iptype(nbeta_))
IF (10*iver(1)+iver(2)>=72) &
READ (iunit) (iptype(j),j=1,nbeta_),npf,ptryc
!
allocate(vloc_(mesh_))
read (iunit) rcloc_,(vloc_(i),i=1,mesh_)
ALLOCATE(vloc_(mesh_))
READ (iunit) rcloc_,(vloc_(i),i=1,mesh_)
!
allocate(rho_atc_(mesh_))
if (ifpcor.gt.0) then
read (iunit) rpcor
read (iunit) (rho_atc_(i),i=1,mesh_)
end if
ALLOCATE(rho_atc_(mesh_))
IF (ifpcor>0) THEN
READ (iunit) rpcor
READ (iunit) (rho_atc_(i),i=1,mesh_)
ENDIF
!
allocate(rho_at_(mesh_), vloc(mesh_))
read (iunit) (vloc(i),i=1,mesh_)
read (iunit) (rho_at_(i),i=1,mesh_)
ALLOCATE(rho_at_(mesh_), vloc(mesh_))
READ (iunit) (vloc(i),i=1,mesh_)
READ (iunit) (rho_at_(i),i=1,mesh_)
allocate(r_(mesh_), rab_(mesh_))
read (iunit) (r_(i),i=1,mesh_)
read (iunit) (rab_(i),i=1,mesh_)
if (iver(1) .ge. 6) then
ALLOCATE(r_(mesh_), rab_(mesh_))
READ (iunit) (r_(i),i=1,mesh_)
READ (iunit) (rab_(i),i=1,mesh_)
IF (iver(1) >= 6) THEN
nchi = nvalps
if (iver(1) .ge. 7) read (iunit) nchi
allocate(wf(mesh_,nchi))
read (iunit) ((wf(i,j), i=1,mesh_),j=1,nchi)
end if
IF (iver(1) >= 7) READ (iunit) nchi
ALLOCATE(wf(mesh_,nchi))
READ (iunit) ((wf(i,j), i=1,mesh_),j=1,nchi)
ENDIF
!
! ----------------------------------------------------------
write (6,'(a)') 'Pseudopotential successfully read'
WRITE (6,'(a)') 'Pseudopotential successfully read'
! ----------------------------------------------------------
!
end subroutine read_uspp
END SUBROUTINE read_uspp
! ----------------------------------------------------------
! ----------------------------------------------------------
subroutine read_vdb(iunit)
SUBROUTINE read_vdb(iunit)
! ----------------------------------------------------------
!
use Vanderbilt
implicit none
integer :: iunit
USE Vanderbilt
IMPLICIT NONE
INTEGER :: iunit
!
integer :: i, j, k, lp
INTEGER :: i, j, k, lp
real(8) :: rinner1
!
!
read(iunit, *) (iver(i),i=1,3),(idmy(i),i=1,3)
read(iunit,'(a20,3f15.9)' ) title, z_, zp_, exfact
read(iunit, *) nvalps, mesh_, etot
READ(iunit, *) (iver(i),i=1,3),(idmy(i),i=1,3)
READ(iunit,'(a20,3f15.9)' ) title, z_, zp_, exfact
READ(iunit, *) nvalps, mesh_, etot
allocate(nnlzps(nvalps), wwnlps(nvalps), eeps(nvalps))
do i = 1,nvalps
read(iunit, *) nnlzps(i), wwnlps(i), eeps(i)
end do
ALLOCATE(nnlzps(nvalps), wwnlps(nvalps), eeps(nvalps))
DO i = 1,nvalps
READ(iunit, *) nnlzps(i), wwnlps(i), eeps(i)
ENDDO
read(iunit, *) keyps, ifpcor, rinner1
READ(iunit, *) keyps, ifpcor, rinner1
if ( iver(1) .eq. 1 ) then
IF ( iver(1) == 1 ) THEN
nang = nvalps
nqf_ = 3
nlc = 5
elseif ( iver(1) .eq. 2 ) then
ELSEIF ( iver(1) == 2 ) THEN
nang = nvalps
nqf_ = 3
nlc = 2 * nvalps - 1
else if ( iver(1) .ge. 3 ) then
read(iunit, *) nang, lloc, eloc, ifqopt, nqf_, qtryc
ELSEIF ( iver(1) >= 3 ) THEN
READ(iunit, *) nang, lloc, eloc, ifqopt, nqf_, qtryc
nlc = 2 * nang - 1
endif
ENDIF
allocate(rinner_(2*nang-1))
ALLOCATE(rinner_(2*nang-1))
rinner_(1) = rinner1
if (10*iver(1)+iver(2).ge.51) &
read (iunit, *) (rinner_(i),i=1,nang*2-1)
if ( iver(1) .ge. 4 ) then
read (iunit, *) irel
else
IF (10*iver(1)+iver(2)>=51) &
READ (iunit, *) (rinner_(i),i=1,nang*2-1)
IF ( iver(1) >= 4 ) THEN
READ (iunit, *) irel
ELSE
irel = 0
end if
ENDIF
allocate(rc(nang))
read(iunit, *) ( rc(i), i=1,nang)
ALLOCATE(rc(nang))
READ(iunit, *) ( rc(i), i=1,nang)
read (iunit,* ) nbeta_, kkbeta
READ (iunit,* ) nbeta_, kkbeta
allocate(beta(kkbeta,nbeta_))
allocate(qfunc_(kkbeta,nbeta_,nbeta_))
allocate(ddd0(nbeta_,nbeta_))
allocate(ddd (nbeta_,nbeta_))
allocate(qqq_(nbeta_,nbeta_))
allocate(lll_(nbeta_))
allocate(eee (nbeta_))
allocate(qfcoef_(nqf_,nlc,nbeta_,nbeta_))
ALLOCATE(beta(kkbeta,nbeta_))
ALLOCATE(qfunc_(kkbeta,nbeta_,nbeta_))
ALLOCATE(ddd0(nbeta_,nbeta_))
ALLOCATE(ddd (nbeta_,nbeta_))
ALLOCATE(qqq_(nbeta_,nbeta_))
ALLOCATE(lll_(nbeta_))
ALLOCATE(eee (nbeta_))
ALLOCATE(qfcoef_(nqf_,nlc,nbeta_,nbeta_))
do j=1,nbeta_
read ( iunit, *) lll_(j)
read ( iunit, *) eee(j), ( beta(i,j), i=1,kkbeta )
do k=j,nbeta_
read( iunit, *) ddd0(j,k), ddd(j,k), qqq_(j,k), &
DO j=1,nbeta_
READ ( iunit, *) lll_(j)
READ ( iunit, *) eee(j), ( beta(i,j), i=1,kkbeta )
DO k=j,nbeta_
READ( iunit, *) ddd0(j,k), ddd(j,k), qqq_(j,k), &
(qfunc_(i,j,k),i=1,kkbeta),&
((qfcoef_(i,lp,j,k),i=1,nqf_),lp=1,2*nang-1)
enddo
enddo
ENDDO
ENDDO
allocate(iptype(nbeta_))
if (10*iver(1)+iver(2).ge.72) then
read ( iunit, * ) (iptype(i), i=1,nbeta_)
read ( iunit, * ) npf, ptryc
end if
ALLOCATE(iptype(nbeta_))
IF (10*iver(1)+iver(2)>=72) THEN
READ ( iunit, * ) (iptype(i), i=1,nbeta_)
READ ( iunit, * ) npf, ptryc
ENDIF
allocate(vloc_(mesh_))
read(iunit, *) rcloc_, ( vloc_(i), i=1,mesh_)
ALLOCATE(vloc_(mesh_))
READ(iunit, *) rcloc_, ( vloc_(i), i=1,mesh_)
allocate(rho_atc_(mesh_))
if ( ifpcor.gt.0 ) then
read(iunit, *) rpcor
read(iunit, *) ( rho_atc_(i), i=1,mesh_)
endif
ALLOCATE(rho_atc_(mesh_))
IF ( ifpcor>0 ) THEN
READ(iunit, *) rpcor
READ(iunit, *) ( rho_atc_(i), i=1,mesh_)
ENDIF
allocate(rho_at_(mesh_), vloc(mesh_))
read(iunit, *) (vloc(i), i=1,mesh_)
read(iunit, *) (rho_at_(i), i=1,mesh_)
ALLOCATE(rho_at_(mesh_), vloc(mesh_))
READ(iunit, *) (vloc(i), i=1,mesh_)
READ(iunit, *) (rho_at_(i), i=1,mesh_)
allocate(r_(mesh_),rab_(mesh_))
read(iunit, *) (r_(i), i=1,mesh_)
read(iunit, *) (rab_(i),i=1,mesh_)
ALLOCATE(r_(mesh_),rab_(mesh_))
READ(iunit, *) (r_(i), i=1,mesh_)
READ(iunit, *) (rab_(i),i=1,mesh_)
if (iver(1) .ge. 6) then
IF (iver(1) >= 6) THEN
nchi = nvalps
if (iver(1) .ge. 7) read (iunit, *) nchi
allocate(wf(mesh_,nchi))
read (iunit, *) ((wf(i,j), i=1,mesh_),j=1,nchi)
end if
IF (iver(1) >= 7) READ (iunit, *) nchi
ALLOCATE(wf(mesh_,nchi))
READ (iunit, *) ((wf(i,j), i=1,mesh_),j=1,nchi)
ENDIF
return
end subroutine read_vdb
RETURN
END SUBROUTINE read_vdb
subroutine convert_uspp
SUBROUTINE convert_uspp
! ----------------------------------------------------------
!
use Vanderbilt
use constants, only : fpi
use upf
implicit none
integer i
character(len=1), dimension(0:3) :: convel=(/'S','P','D','F'/)
USE Vanderbilt
USE constants, ONLY : fpi
USE upf
IMPLICIT NONE
INTEGER i
CHARACTER(len=1), DIMENSION(0:3) :: convel=(/'S','P','D','F'/)
write(generated, '("Generated using Vanderbilt code, version ",3i3)') iver
write(date_author,'("Author: unknown Generation date:",3i5)') idmy
write(comment,'("Automatically converted from original format")')
if (irel == 0) then
WRITE(generated, '("Generated using Vanderbilt code, version ",3i3)') iver
WRITE(date_author,'("Author: unknown Generation date:",3i5)') idmy
WRITE(comment,'("Automatically converted from original format")')
IF (irel == 0) THEN
rel = 0
else if (irel == 1) then
ELSEIF (irel == 1) THEN
rel = 2
else if (irel == 2) then
ELSEIF (irel == 2) THEN
rel = 1
end if
ENDIF
rcloc = rcloc_
nwfs = nvalps
allocate( els(nwfs), oc(nwfs), epseu(nwfs))
allocate(lchi(nwfs), nns(nwfs) )
allocate(rcut (nwfs), rcutus (nwfs))
do i=1, nwfs
ALLOCATE( els(nwfs), oc(nwfs), epseu(nwfs))
ALLOCATE(lchi(nwfs), nns(nwfs) )
ALLOCATE(rcut (nwfs), rcutus (nwfs))
DO i=1, nwfs
nns (i) = nnlzps(i)/100
lchi(i) = mod (nnlzps(i)/10,10)
rcut(i) = rinner1_
rcutus(i)= rc(lchi(i)+1)
oc (i) = wwnlps(i)
write(els(i),'(i1,a1)') nns(i), convel(lchi(i))
WRITE(els(i),'(i1,a1)') nns(i), convel(lchi(i))
epseu(i) = eeps(i)
end do
deallocate (nnlzps, rc, wwnlps, eeps)
ENDDO
DEALLOCATE (nnlzps, rc, wwnlps, eeps)
psd = title
if (keyps.le.2) then
IF (keyps<=2) THEN
pseudotype = 'NC'
else
ELSE
pseudotype = 'US'
end if
nlcc = ifpcor.gt.0
ENDIF
nlcc = ifpcor>0
zp = zp_
etotps = etot
ecutrho=0.0d0
@ -280,87 +280,87 @@ subroutine convert_uspp
lmax = nang - 1
mesh = mesh_
nbeta = nbeta_
if (nvalps .ne. nchi) then
print *, 'WARNING: verify info on atomic wavefunctions'
end if
IF (nvalps /= nchi) THEN
PRINT *, 'WARNING: verify info on atomic wavefunctions'
ENDIF
ntwfc = nchi
allocate( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
do i=1, min(ntwfc,nwfs)
ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
DO i=1, min(ntwfc,nwfs)
elsw(i) = els(i)
ocw(i) = oc (i)
lchiw(i)=lchi(i)
end do
if ( exfact.eq.0) then
ENDDO
IF ( exfact==0) THEN
iexch=1; icorr=1; igcx=0; igcc=0 ! Perdew-Zunger
else if ( exfact.eq.1) then
ELSEIF ( exfact==1) THEN
iexch=1; icorr=3; igcx=1; igcc=3 ! Becke-Lee-Yang-Parr
else if ( exfact.eq.2) then
ELSEIF ( exfact==2) THEN
iexch=1; icorr=1; igcx=1; igcc=0 ! Becke88 exchange
else if (exfact.eq.-5.or.exfact.eq.3) then
ELSEIF (exfact==-5.or.exfact==3) THEN
iexch=1; icorr=1; igcx=1; igcc=1 ! Becke88-Perdew 86
else if (exfact.eq.-6.or.exfact.eq.4) then
ELSEIF (exfact==-6.or.exfact==4) THEN
iexch=1; icorr=4; igcx=2; igcc=2 ! Perdew-Wang 91
else if (exfact.eq. 5) then
ELSEIF (exfact== 5) THEN
iexch=1; icorr=4; igcx=3; igcc=4 ! Perdew-Becke-Erkerhof
else
write (6,'("convert: wrong xc in pseudopotential ",f12.6)') exfact
stop
end if
ELSE
WRITE (6,'("convert: wrong xc in pseudopotential ",f12.6)') exfact
STOP
ENDIF
allocate (r(mesh), rab(mesh))
ALLOCATE (r(mesh), rab(mesh))
r = r_
rab=rab_
deallocate (r_, rab_)
DEALLOCATE (r_, rab_)
allocate (rho_atc(mesh))
ALLOCATE (rho_atc(mesh))
! Vanderbilt rho_core(r) = 4pi*r^2*rho_core(r) UPF
rho_atc (1) = 0.d0
rho_atc (2:mesh) = rho_atc_(2:mesh) / fpi / r(2:mesh)**2
deallocate (rho_atc_)
DEALLOCATE (rho_atc_)
allocate (vloc0(mesh))
ALLOCATE (vloc0(mesh))
vloc0(2:mesh) = vloc_(2:mesh)/r(2:mesh)
vloc0(1) = vloc0(2)
deallocate (vloc_)
DEALLOCATE (vloc_)
allocate(ikk2(nbeta), lll(nbeta))
ALLOCATE(ikk2(nbeta), lll(nbeta))
ikk2 = kkbeta
lll = lll_
deallocate (lll_)
allocate(betar(kkbeta,nbeta))
DEALLOCATE (lll_)
ALLOCATE(betar(kkbeta,nbeta))
betar = beta
deallocate (beta)
DEALLOCATE (beta)
allocate(dion(nbeta,nbeta))
ALLOCATE(dion(nbeta,nbeta))
dion = ddd0
deallocate (ddd0)
DEALLOCATE (ddd0)
allocate(qqq(nbeta,nbeta))
ALLOCATE(qqq(nbeta,nbeta))
qqq = qqq_
deallocate (qqq_)
DEALLOCATE (qqq_)
allocate(qfunc(mesh,nbeta,nbeta))
ALLOCATE(qfunc(mesh,nbeta,nbeta))
qfunc(1:kkbeta,:,:) = qfunc_(1:kkbeta,:,:)
qfunc(kkbeta+1:mesh,:,:) = 0.d0
deallocate (qfunc_)
DEALLOCATE (qfunc_)
nqf = nqf_
nqlc= nlc
allocate(rinner(nqlc))
ALLOCATE(rinner(nqlc))
rinner = rinner_
deallocate(rinner_)
allocate(qfcoef(nqf,nqlc,nbeta,nbeta))
DEALLOCATE(rinner_)
ALLOCATE(qfcoef(nqf,nqlc,nbeta,nbeta))
qfcoef = qfcoef_
deallocate (qfcoef_)
DEALLOCATE (qfcoef_)
allocate (rho_at(mesh))
ALLOCATE (rho_at(mesh))
rho_at = rho_at_
deallocate (rho_at_)
DEALLOCATE (rho_at_)
allocate (chi(mesh,ntwfc))
ALLOCATE (chi(mesh,ntwfc))
chi = wf
deallocate (wf)
DEALLOCATE (wf)
return
end subroutine convert_uspp
RETURN
END SUBROUTINE convert_uspp

28
upftools/vdb2upf.f90
View File

@ -7,32 +7,32 @@
!
!
!---------------------------------------------------------------------
program vdb2upf
PROGRAM vdb2upf
!---------------------------------------------------------------------
!
! Convert a pseudopotential written in Vanderbilt format
! (formatted) to unified pseudopotential format
!
implicit none
character(len=256) filein, fileout
IMPLICIT NONE
CHARACTER(len=256) filein, fileout
!
!
call get_file ( filein )
open(unit=1,file=filein,status='old',form='formatted')
call read_vdb(1)
close (unit=1)
CALL get_file ( filein )
OPEN(unit=1,file=filein,status='old',form='formatted')
CALL read_vdb(1)
CLOSE (unit=1)
! convert variables read from Vanderbilt format into those needed
! by the upf format - add missing quantities
call convert_uspp
CALL convert_uspp
fileout=trim(filein)//'.UPF'
print '(''Output PP file in UPF format : '',a)', fileout
PRINT '(''Output PP file in UPF format : '',a)', fileout
open(unit=2,file=fileout,status='unknown',form='formatted')
call write_upf(2)
close (unit=2)
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
CALL write_upf(2)
CLOSE (unit=2)
stop
end program vdb2upf
STOP
END PROGRAM vdb2upf

708
upftools/virtual.f90
View File

File diff suppressed because it is too large Load Diff

494
upftools/write_upf.f90
View File

@ -5,445 +5,445 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
module upf
MODULE upf
!
! All variables to be written into the UPF file
! (UPF = unified pseudopotential format)
!
! pp_info
integer :: rel
INTEGER :: rel
real(8) :: rcloc
integer :: nwfs
real(8), allocatable :: oc(:), rcut(:), rcutus(:), epseu(:)
character(len=2), allocatable :: els(:)
integer, allocatable:: lchi (:), nns (:)
INTEGER :: nwfs
real(8), ALLOCATABLE :: oc(:), rcut(:), rcutus(:), epseu(:)
CHARACTER(len=2), ALLOCATABLE :: els(:)
INTEGER, ALLOCATABLE:: lchi (:), nns (:)
!
! pp_header
character (len=80):: generated, date_author, comment
character (len=2) :: psd, pseudotype
integer :: nv = 0
integer :: iexch, icorr, igcx, igcc
integer :: lmax, mesh, nbeta, ntwfc
logical :: nlcc
CHARACTER (len=80):: generated, date_author, comment
CHARACTER (len=2) :: psd, pseudotype
INTEGER :: nv = 0
INTEGER :: iexch, icorr, igcx, igcc
INTEGER :: lmax, mesh, nbeta, ntwfc
LOGICAL :: nlcc
real(8) :: zp, ecutrho, ecutwfc, etotps
real(8), allocatable :: ocw(:)
character(len=2), allocatable :: elsw(:)
integer, allocatable:: lchiw(:)
real(8), ALLOCATABLE :: ocw(:)
CHARACTER(len=2), ALLOCATABLE :: elsw(:)
INTEGER, ALLOCATABLE:: lchiw(:)
!
! pp_mesh
real(8), allocatable :: r(:), rab(:)
real(8), ALLOCATABLE :: r(:), rab(:)
!
! pp_nlcc
real(8), allocatable :: rho_atc(:)
real(8), ALLOCATABLE :: rho_atc(:)
!
! pp_local
real(8), allocatable :: vloc0(:)
real(8), ALLOCATABLE :: vloc0(:)
!
! pp_nonlocal
! pp_beta
real(8), allocatable :: betar(:,:)
integer, allocatable:: lll(:), ikk2(:)
real(8), ALLOCATABLE :: betar(:,:)
INTEGER, ALLOCATABLE:: lll(:), ikk2(:)
! pp_dij
real(8), allocatable :: dion(:,:)
real(8), ALLOCATABLE :: dion(:,:)
! pp_qij
integer :: nqf, nqlc
real(8), allocatable :: rinner(:), qqq(:,:), qfunc(:,:,:)
INTEGER :: nqf, nqlc
real(8), ALLOCATABLE :: rinner(:), qqq(:,:), qfunc(:,:,:)
! pp_qfcoef
real(8), allocatable :: qfcoef(:,:,:,:)
real(8), ALLOCATABLE :: qfcoef(:,:,:,:)
!
! pp_pswfc
real(8), allocatable :: chi(:,:)
real(8), ALLOCATABLE :: chi(:,:)
!
! pp_rhoatom
real(8), allocatable :: rho_at(:)
end module upf
real(8), ALLOCATABLE :: rho_at(:)
END MODULE upf
!
subroutine write_upf(ounps)
SUBROUTINE write_upf(ounps)
use upf, only: nlcc
USE upf, ONLY: nlcc
integer :: ounps
INTEGER :: ounps
call write_pseudo_comment(ounps)
call write_pseudo_header(ounps)
call write_pseudo_mesh(ounps)
if (nlcc) call write_pseudo_nlcc(ounps)
call write_pseudo_local(ounps)
call write_pseudo_nl(ounps)
call write_pseudo_pswfc(ounps)
call write_pseudo_rhoatom(ounps)
CALL write_pseudo_comment(ounps)
CALL write_pseudo_header(ounps)
CALL write_pseudo_mesh(ounps)
IF (nlcc) CALL write_pseudo_nlcc(ounps)
CALL write_pseudo_local(ounps)
CALL write_pseudo_nl(ounps)
CALL write_pseudo_pswfc(ounps)
CALL write_pseudo_rhoatom(ounps)
!
print '("*** PLEASE TEST BEFORE USING!!! ***")'
print '("review the content of the PP_INFO fields")'
PRINT '("*** PLEASE TEST BEFORE USING!!! ***")'
PRINT '("review the content of the PP_INFO fields")'
!
end subroutine write_upf
END SUBROUTINE write_upf
!
!---------------------------------------------------------------------
subroutine write_pseudo_comment (ounps)
SUBROUTINE write_pseudo_comment (ounps)
!---------------------------------------------------------------------
!
!
! This routine writes the comments of the new UPF file
!
use upf
implicit none
integer :: ounps
USE upf
IMPLICIT NONE
INTEGER :: ounps
integer :: nb, ios
INTEGER :: nb, ios
write (ounps, '(a9)', err = 100, iostat = ios) "<PP_INFO>"
WRITE (ounps, '(a9)', err = 100, iostat = ios) "<PP_INFO>"
write (ounps, '(a)', err = 100, iostat = ios) generated
write (ounps, '(a)', err = 100, iostat = ios) date_author
write (ounps, '(a)', err = 100, iostat = ios) comment
if (rel==2) then
write (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel,&
WRITE (ounps, '(a)', err = 100, iostat = ios) generated
WRITE (ounps, '(a)', err = 100, iostat = ios) date_author
WRITE (ounps, '(a)', err = 100, iostat = ios) comment
IF (rel==2) THEN
WRITE (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel,&
&"The Pseudo was generated with a Full-Relativistic Calculation"
else if (rel==1) then
write (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel,&
ELSEIF (rel==1) THEN
WRITE (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel,&
&"The Pseudo was generated with a Scalar-Relativistic Calculation"
else if (rel==0) then
write (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel, &
ELSEIF (rel==0) THEN
WRITE (ounps, '(i5,t14,a)', err = 100, iostat = ios) rel, &
& "The Pseudo was generated with a Non-Relativistic Calculation"
endif
ENDIF
if (rcloc > 0.d0) &
write (ounps, '(1pe19.11,t24,a)', err = 100, iostat = ios) &
IF (rcloc > 0.d0) &
WRITE (ounps, '(1pe19.11,t24,a)', err = 100, iostat = ios) &
rcloc, "Local Potential cutoff radius"
if (nwfs>0) &
write (ounps, '(a2,2a3,a6,3a19)', err = 100, iostat = ios) "nl", &
IF (nwfs>0) &
WRITE (ounps, '(a2,2a3,a6,3a19)', err = 100, iostat = ios) "nl", &
&" pn", "l", "occ", "Rcut", "Rcut US", "E pseu"
do nb = 1, nwfs
write (ounps, '(a2,2i3,f6.2,3f19.11)') els (nb) , nns (nb) , &
DO nb = 1, nwfs
WRITE (ounps, '(a2,2i3,f6.2,3f19.11)') els (nb) , nns (nb) , &
lchi (nb) , oc (nb) , rcut (nb) , rcutus (nb) , epseu(nb)
enddo
ENDDO
write (ounps, '(a10)', err = 100, iostat = ios) "</PP_INFO>"
return
100 write(6,'("write_pseudo_comment: error writing pseudopotential file")')
stop
WRITE (ounps, '(a10)', err = 100, iostat = ios) "</PP_INFO>"
RETURN
100 WRITE(6,'("write_pseudo_comment: error writing pseudopotential file")')
STOP
end subroutine write_pseudo_comment
END SUBROUTINE write_pseudo_comment
!
!---------------------------------------------------------------------
subroutine write_pseudo_header (ounps)
SUBROUTINE write_pseudo_header (ounps)
!---------------------------------------------------------------------
!
!
! This routine writes the header of the new UPF file
!
use upf
implicit none
integer :: ounps
USE upf
IMPLICIT NONE
INTEGER :: ounps
!
character (len=4) :: shortname
character (len=20):: dft
integer :: nb, ios
CHARACTER (len=4) :: shortname
CHARACTER (len=20):: dft
INTEGER :: nb, ios
!
!
write (ounps, '(//a11)', err = 100, iostat = ios) "<PP_HEADER>"
WRITE (ounps, '(//a11)', err = 100, iostat = ios) "<PP_HEADER>"
write (ounps, '(t3,i2,t24,a)', err = 100, iostat = ios) nv, &
WRITE (ounps, '(t3,i2,t24,a)', err = 100, iostat = ios) nv, &
"Version Number"
write (ounps, '(t3,a,t24,a)', err = 100, iostat = ios) psd , &
WRITE (ounps, '(t3,a,t24,a)', err = 100, iostat = ios) psd , &
"Element"
if (pseudotype == 'NC') then
write (ounps, '(a5,t24,a)', err = 100, iostat = ios) "NC", &
IF (pseudotype == 'NC') THEN
WRITE (ounps, '(a5,t24,a)', err = 100, iostat = ios) "NC", &
"Norm - Conserving pseudopotential"
else if (pseudotype == 'US') then
write (ounps, '(a5,t24,a)', err = 100, iostat = ios) "US", &
ELSEIF (pseudotype == 'US') THEN
WRITE (ounps, '(a5,t24,a)', err = 100, iostat = ios) "US", &
"Ultrasoft pseudopotential"
else
write(6,'("write_pseudo_header: unknown PP type ",A)') pseudotype
stop
endif
write (ounps, '(l5,t24,a)', err = 100, iostat = ios) nlcc , &
ELSE
WRITE(6,'("write_pseudo_header: unknown PP type ",A)') pseudotype
STOP
ENDIF
WRITE (ounps, '(l5,t24,a)', err = 100, iostat = ios) nlcc , &
"Nonlinear Core Correction"
call dftname (iexch, icorr, igcx, igcc, dft, shortname)
write (ounps, '(a,t24,a4,a)', err = 100, iostat = ios) &
CALL dftname (iexch, icorr, igcx, igcc, dft, shortname)
WRITE (ounps, '(a,t24,a4,a)', err = 100, iostat = ios) &
dft, shortname," Exchange-Correlation functional"
write (ounps, '(f17.11,t24,a)') zp , "Z valence"
write (ounps, '(f17.11,t24,a)') etotps, "Total energy"
write (ounps, '(2f11.7,t24,a)') ecutrho, ecutwfc, &
WRITE (ounps, '(f17.11,t24,a)') zp , "Z valence"
WRITE (ounps, '(f17.11,t24,a)') etotps, "Total energy"
WRITE (ounps, '(2f11.7,t24,a)') ecutrho, ecutwfc, &
"Suggested cutoff for wfc and rho"
write (ounps, '(i5,t24,a)') lmax, "Max angular momentum component"
write (ounps, '(i5,t24,a)') mesh, "Number of points in mesh"
write (ounps, '(2i5,t24,a)', err = 100, iostat = ios) ntwfc, &
WRITE (ounps, '(i5,t24,a)') lmax, "Max angular momentum component"
WRITE (ounps, '(i5,t24,a)') mesh, "Number of points in mesh"
WRITE (ounps, '(2i5,t24,a)', err = 100, iostat = ios) ntwfc, &
nbeta , "Number of Wavefunctions, Number of Projectors"
write (ounps, '(a,t24,a2,a3,a6)', err = 100, iostat = ios) &
WRITE (ounps, '(a,t24,a2,a3,a6)', err = 100, iostat = ios) &
" Wavefunctions", "nl", "l", "occ"
do nb = 1, ntwfc
write (ounps, '(t24,a2,i3,f6.2)') elsw(nb), lchiw(nb), ocw(nb)
enddo
DO nb = 1, ntwfc
WRITE (ounps, '(t24,a2,i3,f6.2)') elsw(nb), lchiw(nb), ocw(nb)
ENDDO
!---> End header writing
write (ounps, '(a12)', err = 100, iostat = ios) "</PP_HEADER>"
return
100 write(6,'("write_pseudo_header: error writing pseudopotential file")')
stop
WRITE (ounps, '(a12)', err = 100, iostat = ios) "</PP_HEADER>"
RETURN
100 WRITE(6,'("write_pseudo_header: error writing pseudopotential file")')
STOP
end subroutine write_pseudo_header
END SUBROUTINE write_pseudo_header
!
!---------------------------------------------------------------------
subroutine write_pseudo_mesh (ounps)
SUBROUTINE write_pseudo_mesh (ounps)
!---------------------------------------------------------------------
!
!
! This routine writes the atomic charge density to the new UPF file
!
use upf
implicit none
integer :: ounps
USE upf
IMPLICIT NONE
INTEGER :: ounps
!
integer :: ir, ios
INTEGER :: ir, ios
!
write (ounps, '(//a9)', err = 100, iostat = ios) "<PP_MESH>"
WRITE (ounps, '(//a9)', err = 100, iostat = ios) "<PP_MESH>"
write (ounps, '(t3,a6)', err = 100, iostat = ios) "<PP_R>"
write (ounps, '(1p4e19.11)', err=100, iostat=ios) (r(ir), ir=1,mesh )
write (ounps, '(t3,a7)', err = 100, iostat = ios) "</PP_R>"
write (ounps, '(t3,a8)', err = 100, iostat = ios) "<PP_RAB>"
write (ounps, '(1p4e19.11)', err=100, iostat=ios) (rab(ir), ir=1,mesh )
write (ounps, '(t3,a9)', err = 100, iostat = ios) "</PP_RAB>"
WRITE (ounps, '(t3,a6)', err = 100, iostat = ios) "<PP_R>"
WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) (r(ir), ir=1,mesh )
WRITE (ounps, '(t3,a7)', err = 100, iostat = ios) "</PP_R>"
WRITE (ounps, '(t3,a8)', err = 100, iostat = ios) "<PP_RAB>"
WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) (rab(ir), ir=1,mesh )
WRITE (ounps, '(t3,a9)', err = 100, iostat = ios) "</PP_RAB>"
write (ounps, '(a10)', err = 100, iostat = ios) "</PP_MESH>"
WRITE (ounps, '(a10)', err = 100, iostat = ios) "</PP_MESH>"
return
RETURN
100 write(6,'("write_pseudo_mesh: error writing pseudopotential file")')
stop
100 WRITE(6,'("write_pseudo_mesh: error writing pseudopotential file")')
STOP
end subroutine write_pseudo_mesh
END SUBROUTINE write_pseudo_mesh
!
!---------------------------------------------------------------------
subroutine write_pseudo_nlcc (ounps)
SUBROUTINE write_pseudo_nlcc (ounps)
!---------------------------------------------------------------------
!
!
! This routine writes the core charge for the nonlinear core
! correction of the new UPF file
!
use upf
implicit none
integer :: ounps
USE upf
IMPLICIT NONE
INTEGER :: ounps
!
integer :: ir, ios
INTEGER :: ir, ios
write (ounps, '(//a9)', err = 100, iostat = ios) "<PP_NLCC>"
WRITE (ounps, '(//a9)', err = 100, iostat = ios) "<PP_NLCC>"
write (ounps, '(1p4e19.11)', err=100, iostat=ios) &
WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) &
( rho_atc(ir), ir = 1, mesh )
write (ounps, '(a10)', err = 100, iostat = ios) "</PP_NLCC>"
return
WRITE (ounps, '(a10)', err = 100, iostat = ios) "</PP_NLCC>"
RETURN
100 write(6,'("write_pseudo_nlcc: error writing pseudopotential file")')
stop
100 WRITE(6,'("write_pseudo_nlcc: error writing pseudopotential file")')
STOP
end subroutine write_pseudo_nlcc
END SUBROUTINE write_pseudo_nlcc
!
!---------------------------------------------------------------------
subroutine write_pseudo_local (ounps)
SUBROUTINE write_pseudo_local (ounps)
!---------------------------------------------------------------------
!
!
! This routine writes the local part of the new UPF file
!
use upf
implicit none
integer :: ounps
USE upf
IMPLICIT NONE
INTEGER :: ounps
!
integer :: ir, ios
INTEGER :: ir, ios
write (ounps, '(//a10)', err = 100, iostat = ios) "<PP_LOCAL>"
write (ounps, '(1p4e19.11)', err=100, iostat=ios) &
WRITE (ounps, '(//a10)', err = 100, iostat = ios) "<PP_LOCAL>"
WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) &
( vloc0(ir), ir = 1, mesh )
write (ounps, '(a11)', err = 100, iostat = ios) "</PP_LOCAL>"
return
100 write(6,'("write_pseudo_local: error writing pseudopotential file")')
stop
WRITE (ounps, '(a11)', err = 100, iostat = ios) "</PP_LOCAL>"
RETURN
100 WRITE(6,'("write_pseudo_local: error writing pseudopotential file")')
STOP
end subroutine write_pseudo_local
END SUBROUTINE write_pseudo_local
!
!---------------------------------------------------------------------
subroutine write_pseudo_nl (ounps)
SUBROUTINE write_pseudo_nl (ounps)
!---------------------------------------------------------------------
!
!
! This routine writes the non local part of the new UPF file
!
use upf
implicit none
integer :: ounps
USE upf
IMPLICIT NONE
INTEGER :: ounps
!
integer :: nb, mb, n, ir, nd, i, lp, ios
INTEGER :: nb, mb, n, ir, nd, i, lp, ios
write (ounps, '(//a13)', err = 100, iostat = ios) "<PP_NONLOCAL>"
do nb = 1, nbeta
write (ounps, '(t3,a9)', err = 100, iostat = ios) "<PP_BETA>"
write (ounps, '(2i5,t24,a)', err=100, iostat=ios) &
WRITE (ounps, '(//a13)', err = 100, iostat = ios) "<PP_NONLOCAL>"
DO nb = 1, nbeta
WRITE (ounps, '(t3,a9)', err = 100, iostat = ios) "<PP_BETA>"
WRITE (ounps, '(2i5,t24,a)', err=100, iostat=ios) &
nb, lll(nb), "Beta L"
write (ounps, '(i6)', err=100, iostat=ios) ikk2 (nb)
write (ounps, '(1p4e19.11)', err=100, iostat=ios) &
WRITE (ounps, '(i6)', err=100, iostat=ios) ikk2 (nb)
WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) &
( betar(ir,nb), ir=1,ikk2(nb) )
write (ounps, '(t3,a10)', err = 100, iostat = ios) "</PP_BETA>"
enddo
WRITE (ounps, '(t3,a10)', err = 100, iostat = ios) "</PP_BETA>"
ENDDO
write (ounps, '(t3,a8)', err = 100, iostat = ios) "<PP_DIJ>"
WRITE (ounps, '(t3,a8)', err = 100, iostat = ios) "<PP_DIJ>"
nd = 0
do nb = 1, nbeta
do mb = nb, nbeta
if ( abs(dion(nb,mb)) .gt. 1.0d-12 ) nd = nd + 1
enddo
enddo
write (ounps, '(1p,i5,t24,a)', err=100, iostat=ios) &
DO nb = 1, nbeta
DO mb = nb, nbeta
IF ( abs(dion(nb,mb)) > 1.0d-12 ) nd = nd + 1
ENDDO
ENDDO
WRITE (ounps, '(1p,i5,t24,a)', err=100, iostat=ios) &
nd, "Number of nonzero Dij"
do nb = 1, nbeta
do mb = nb, nbeta
if ( abs(dion(nb,mb)) .gt. 1.0d-12 ) &
write(ounps,'(1p,2i5,e19.11)', err=100, iostat=ios) &
DO nb = 1, nbeta
DO mb = nb, nbeta
IF ( abs(dion(nb,mb)) > 1.0d-12 ) &
WRITE(ounps,'(1p,2i5,e19.11)', err=100, iostat=ios) &
nb, mb, dion(nb,mb)
enddo
enddo
write (ounps, '(t3,a9)', err=100, iostat=ios) "</PP_DIJ>"
ENDDO
ENDDO
WRITE (ounps, '(t3,a9)', err=100, iostat=ios) "</PP_DIJ>"
if (pseudotype == 'US') then
write (ounps, '(t3,a8)', err = 100, iostat = ios) "<PP_QIJ>"
write (ounps, '(i5,a)',err=100, iostat=ios) nqf," nqf.&
IF (pseudotype == 'US') THEN
WRITE (ounps, '(t3,a8)', err = 100, iostat = ios) "<PP_QIJ>"
WRITE (ounps, '(i5,a)',err=100, iostat=ios) nqf," nqf.&
& If not zero, Qij's inside rinner are computed using qfcoef's"
if (nqf.gt.0) then
write (ounps, '(t5,a11)', err=100, iostat=ios) "<PP_RINNER>"
write (ounps,'(i5,1pe19.11)', err=100, iostat=ios) &
IF (nqf>0) THEN
WRITE (ounps, '(t5,a11)', err=100, iostat=ios) "<PP_RINNER>"
WRITE (ounps,'(i5,1pe19.11)', err=100, iostat=ios) &
(i, rinner(i), i = 1, nqlc)
write (ounps, '(t5,a12)', err=100, iostat=ios) "</PP_RINNER>"
end if
do nb = 1, nbeta
do mb = nb, nbeta
write (ounps, '(3i5,t24,a)', err=100, iostat=ios) &
WRITE (ounps, '(t5,a12)', err=100, iostat=ios) "</PP_RINNER>"
ENDIF
DO nb = 1, nbeta
DO mb = nb, nbeta
WRITE (ounps, '(3i5,t24,a)', err=100, iostat=ios) &
nb, mb, lll(mb) , "i j (l(j))"
write (ounps, '(1pe19.11,t24,a)', err=100, iostat=ios) &
WRITE (ounps, '(1pe19.11,t24,a)', err=100, iostat=ios) &
qqq(nb,mb), "Q_int"
write (ounps, '(1p4e19.11)', err=100, iostat=ios) &
WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) &
( qfunc (n,nb,mb), n=1,mesh )
if (nqf.gt.0) then
write (ounps, '(t5,a11)', err=100, iostat=ios) &
IF (nqf>0) THEN
WRITE (ounps, '(t5,a11)', err=100, iostat=ios) &
"<PP_QFCOEF>"
write(ounps,'(1p4e19.11)', err=100, iostat=ios) &
WRITE(ounps,'(1p4e19.11)', err=100, iostat=ios) &
((qfcoef(i,lp,nb,mb),i=1,nqf),lp=1,nqlc)
write (ounps, '(t5,a12)', err=100, iostat=ios) &
WRITE (ounps, '(t5,a12)', err=100, iostat=ios) &
"</PP_QFCOEF>"
end if
enddo
enddo
write (ounps, '(t3,a9)', err = 100, iostat = ios) "</PP_QIJ>"
ENDIF
ENDDO
ENDDO
WRITE (ounps, '(t3,a9)', err = 100, iostat = ios) "</PP_QIJ>"
endif
write (ounps, '(a14)', err = 100, iostat = ios) "</PP_NONLOCAL>"
return
ENDIF
WRITE (ounps, '(a14)', err = 100, iostat = ios) "</PP_NONLOCAL>"
RETURN
100 write(6,'("write_pseudo_nl: error writing pseudopotential file")')
stop
100 WRITE(6,'("write_pseudo_nl: error writing pseudopotential file")')
STOP
end subroutine write_pseudo_nl
END SUBROUTINE write_pseudo_nl
!
!---------------------------------------------------------------------
subroutine write_pseudo_pswfc (ounps)
SUBROUTINE write_pseudo_pswfc (ounps)
!---------------------------------------------------------------------
!
!
! This routine writes the pseudo atomic functions
! of the new UPF file
!
use upf
implicit none
integer :: ounps
USE upf
IMPLICIT NONE
INTEGER :: ounps
!
integer :: nb, ir, ios
INTEGER :: nb, ir, ios
write (ounps, '(//a10)', err = 100, iostat = ios) "<PP_PSWFC>"
do nb = 1, ntwfc
write (ounps,'(a2,i5,f6.2,t24,a)', err=100, iostat=ios) &
WRITE (ounps, '(//a10)', err = 100, iostat = ios) "<PP_PSWFC>"
DO nb = 1, ntwfc
WRITE (ounps,'(a2,i5,f6.2,t24,a)', err=100, iostat=ios) &
elsw(nb), lchiw(nb), ocw(nb), "Wavefunction"
write (ounps, '(1p4e19.11)', err=100, iostat=ios) &
WRITE (ounps, '(1p4e19.11)', err=100, iostat=ios) &
( chi(ir,nb), ir=1,mesh )
enddo
write (ounps, '(a11)', err = 100, iostat = ios) "</PP_PSWFC>"
return
ENDDO
WRITE (ounps, '(a11)', err = 100, iostat = ios) "</PP_PSWFC>"
RETURN
100 write(6,'("write_pseudo_pswfc: error writing pseudopotential file")')
stop
100 WRITE(6,'("write_pseudo_pswfc: error writing pseudopotential file")')
STOP
end subroutine write_pseudo_pswfc
END SUBROUTINE write_pseudo_pswfc
!
!---------------------------------------------------------------------
subroutine write_pseudo_rhoatom (ounps)
SUBROUTINE write_pseudo_rhoatom (ounps)
!---------------------------------------------------------------------
!
!
! This routine writes the atomic charge density to the new UPF file
!
use upf
implicit none
integer :: ounps
USE upf
IMPLICIT NONE
INTEGER :: ounps
!
integer :: ir, ios
INTEGER :: ir, ios
write (ounps, '(//a12)', err = 100, iostat = ios) "<PP_RHOATOM>"
write (ounps, '(1p4e19.11)', err = 100, iostat = ios) &
WRITE (ounps, '(//a12)', err = 100, iostat = ios) "<PP_RHOATOM>"
WRITE (ounps, '(1p4e19.11)', err = 100, iostat = ios) &
( rho_at(ir), ir=1,mesh )
write (ounps, '(a13)', err = 100, iostat = ios) "</PP_RHOATOM>"
return
WRITE (ounps, '(a13)', err = 100, iostat = ios) "</PP_RHOATOM>"
RETURN
100 write(6,'("write_pseudo_rhoatom: error writing pseudopotential file")')
stop
100 WRITE(6,'("write_pseudo_rhoatom: error writing pseudopotential file")')
STOP
end subroutine write_pseudo_rhoatom
END SUBROUTINE write_pseudo_rhoatom
!---------------------------------------------------------------------
subroutine dftname(iexch, icorr, igcx, igcc, longname, shortname)
SUBROUTINE dftname(iexch, icorr, igcx, igcc, longname, shortname)
!---------------------------------------------------------------------
implicit none
integer iexch, icorr, igcx, igcc
character (len=4) :: shortname
character (len=20):: longname
IMPLICIT NONE
INTEGER iexch, icorr, igcx, igcc
CHARACTER (len=4) :: shortname
CHARACTER (len=20):: longname
!
! The data used to convert iexch, icorr, igcx, igcc
! into a user-readable string
!
integer, parameter :: nxc = 6, ncc = 9, ngcx = 4, ngcc = 5
character (len=20) :: exc, corr, gradx, gradc
dimension exc (0:nxc), corr (0:ncc), gradx (0:ngcx), gradc (0:ngcc)
INTEGER, PARAMETER :: nxc = 6, ncc = 9, ngcx = 4, ngcc = 5
CHARACTER (len=20) :: exc, corr, gradx, gradc
DIMENSION exc (0:nxc), corr (0:ncc), gradx (0:ngcx), gradc (0:ngcc)
data exc / 'NOX ', 'SLA ', 'SL1 ', 'RXC ', 'OEP ', 'HF ', 'PB0X' /
data corr / 'NOC ', 'PZ ', 'VWN ', 'LYP ', 'PW ', 'WIG ', 'HL ',&
'OBZ ', 'OBW ', 'GL ' /
data gradx / 'NOGX', 'B88 ', 'GGX ', 'PBE ', 'TPSS' /
data gradc / 'NOGC', 'P86 ', 'GGC ', 'BLYP', 'PBE ', 'TPSS' /
if (iexch==1.and.igcx==0.and.igcc==0) then
IF (iexch==1.and.igcx==0.and.igcc==0) THEN
shortname = corr(icorr)
else if (iexch==1.and.icorr==3.and.igcx==1.and.igcc==3) then
ELSEIF (iexch==1.and.icorr==3.and.igcx==1.and.igcc==3) THEN
shortname = 'BLYP'
else if (iexch==1.and.icorr==1.and.igcx==1.and.igcc==0) then
ELSEIF (iexch==1.and.icorr==1.and.igcx==1.and.igcc==0) THEN
shortname = 'B88'
else if (iexch==1.and.icorr==1.and.igcx==1.and.igcc==1) then
ELSEIF (iexch==1.and.icorr==1.and.igcx==1.and.igcc==1) THEN
shortname = 'BP'
else if (iexch==1.and.icorr==4.and.igcx==2.and.igcc==2) then
ELSEIF (iexch==1.and.icorr==4.and.igcx==2.and.igcc==2) THEN
shortname = 'PW91'
else if (iexch==1.and.icorr==4.and.igcx==3.and.igcc==4) then
ELSEIF (iexch==1.and.icorr==4.and.igcx==3.and.igcc==4) THEN
shortname = 'PBE'
else if (iexch==1.and.icorr==4.and.igcx==4.and.igcc==5) then
ELSEIF (iexch==1.and.icorr==4.and.igcx==4.and.igcc==5) THEN
shortname = 'TPSS'
else
ELSE
shortname = ' '
end if
write(longname,'(4a5)') exc(iexch),corr(icorr),gradx(igcx),gradc(igcc)
ENDIF
WRITE(longname,'(4a5)') exc(iexch),corr(icorr),gradx(igcx),gradc(igcc)
return
end subroutine dftname
RETURN
END SUBROUTINE dftname