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

387 lines
9.8 KiB
Fortran
Raw Permalink Blame History

!
! Copyright (C) 2001-2009 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file 'License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!---------------------------------------------------------------------
PROGRAM ncpp2upf
!---------------------------------------------------------------------
!
! Convert a pseudopotential written in PWSCF format
! (norm-conserving) to unified pseudopotential format
IMPLICIT NONE
CHARACTER(len=256) filein, fileout
!
!
CALL get_file ( filein )
OPEN(unit=1,file=filein,status='old',form='formatted')
! underscore added to prevent conflict with read_ncpp in Modules/
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
fileout=trim(filein)//'.UPF'
PRINT '(''Output PP file in US format : '',a)', fileout
OPEN(unit=2,file=fileout,status='unknown',form='formatted')
CALL write_upf_v1(2)
CLOSE (unit=2)
STOP
20 CALL errore ('ncpp2upf', 'Reading pseudo file name ', 1)
END PROGRAM ncpp2upf
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_
real(8) :: zp_
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:: vnl(:,:), rho_atc_(:), rho_at_(:)
INTEGER, ALLOCATABLE:: lchi_(:)
real(8), ALLOCATABLE:: chi_(:,:), oc_(:)
END MODULE ncpp
!
! ----------------------------------------------------------
SUBROUTINE read_ncpp_(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
real (8) :: x, qe_erf
INTEGER :: l, i, ir, nb, n
CHARACTER (len=255) line
EXTERNAL qe_erf
READ(iunps, *, end=300, err=300 ) dft_
IF (dft_(1:2)=='**') dft_ = 'PZ'
READ (iunps, *, err=300) psd_, zp_, lmax_, nlc, nnl, nlcc_, &
lloc, bhstype
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==-1000) lloc=lmax_
!
! In numeric pseudopotentials both nlc and nnl are zero.
!
numeric = nlc<=0 .and. nnl<=0
IF (.not.numeric) THEN
!
! read pseudopotentials in analytic form
!
READ(iunps, *, err=300) &
( alpc(i), i=1, 2 ), ( cc(i), i=1,2 )
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
IF (nlcc_) THEN
READ(iunps, *, err=300) &
a_nlcc, b_nlcc, alpha_nlcc
IF (alpha_nlcc<=0.d0) &
CALL errore('read_ncpp','nlcc but alpha=0',1)
ENDIF
IF (bhstype) CALL bachel(alps,aps,1,lmax_)
ENDIF
READ(iunps, *, err=300) zmesh, xmin, dx, mesh_, nchi
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_))
DO ir = 1, mesh_
x = xmin + dble(ir-1) * dx
r_ (ir) = exp(x) / zmesh
rab_(ir) = dx * r_(ir)
ENDDO
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
ALLOCATE(rho_atc_(mesh_))
IF(nlcc_) THEN
READ(iunps, *, err=300) ( rho_atc_(ir), ir=1,mesh_ )
ENDIF
ELSE
!
! convert analytic to numeric form
!
DO l=0,lmax_
!
! DO NOT USE f90 ARRAY SYNTAX: qe_erf IS NOT AN INTRINSIC FUNCTION!!!
!
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)
ENDDO
DO n=1,nnl
vnl(:,l)= vnl(:,l)+ (aps(n,l)+ aps(n+3,l)*r_(:)**2 )* &
exp(-alps(n,l)*r_(:)**2)
ENDDO
!
! convert to Rydberg
!
vnl(:,l) = vnl(:,l)*2.0d0
ENDDO
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)
rho_atc_ = 0
END WHERE
ENDIF
ENDIF
!
! subtract the local part
!
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
! read wavefunction label and store for later
READ(iunps, '(a)', err=300) line
READ(iunps, *, err=300) lchi_( nb), oc_( nb )
!
! Test lchi and occupation numbers
!
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)/=convel(lchi_(nb))) GOTO 222
DO l = 0, lmax_
IF (label(2:2)==convel(l)) THEN
els(nb) = label(1:2)
GOTO 223
ENDIF
ENDDO
222 CONTINUE
els(nb) = '*'//convel(lchi_(nb))
223 CONTINUE
!
! finally read the wavefunction
READ(iunps, *, err=300) (chi_(ir,nb),ir=1,mesh_)
ENDDO
!
! compute the atomic charges
!
ALLOCATE(rho_at_(mesh_))
rho_at_(:)=0.d0
DO nb = 1, nchi
IF( oc_(nb)/=0.d0) &
rho_at_(:) = rho_at_(:) + oc_(nb)*chi_(:,nb)**2
ENDDO
! ----------------------------------------------------------
WRITE (6,'(a)') 'Pseudopotential successfully read'
! ----------------------------------------------------------
RETURN
300 CALL errore('read_ncpp','pseudo file is empty or wrong',1)
END SUBROUTINE read_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(:)
real(8) :: vll
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")')
comment = 'Info: automatically converted from PWSCF format'
! reasonable assumption
IF (zmesh > 18) THEN
rel = 1
ELSE
rel = 0
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
nns (i) = 0
lchi(i) = lchi_(i)
rcut(i) = 0.0d0
rcutus(i)= 0.0d0
oc (i) = oc_(i)
epseu(i) = 0.0d0
ENDDO
DEALLOCATE (lchi_, oc_)
psd = psd_
pseudotype = 'NC'
nlcc = nlcc_
zp = zp_
etotps = 0.0d0
ecutrho=0.0d0
ecutwfc=0.0d0
IF ( lmax_ == lloc) THEN
lmax = lmax_-1
ELSE
lmax = lmax_
ENDIF
nbeta= lmax_
mesh = mesh_
ntwfc= nchi
ALLOCATE( elsw(ntwfc), ocw(ntwfc), lchiw(ntwfc) )
DO i=1, nchi
lchiw(i) = lchi(i)
ocw(i) = oc(i)
elsw(i) = els(i)
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))
rab = rab_
r = r_
ALLOCATE (rho_atc(mesh))
rho_atc = rho_atc_
DEALLOCATE (rho_atc_)
ALLOCATE (vloc0(mesh))
vloc0(:) = vnl(:,lloc)
IF (nbeta > 0) THEN
ALLOCATE(ikk2(nbeta), lll(nbeta))
kkbeta=mesh
DO ir = 1,mesh
IF ( r(ir) > rmax ) THEN
kkbeta=ir
exit
ENDIF
ENDDO
! make sure kkbeta is odd as required for simpson
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))
qfunc(:,:,:)=0.0d0
dion(:,:) =0.d0
qqq(:,:) =0.d0
iv=0
DO i=1,nchi
l=lchi(i)
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)
ENDDO
CALL simpson(kkbeta,aux,rab,vll)
dion(iv,iv) = 1.0d0/vll
ENDIF
IF(iv >= nbeta) exit ! skip additional pseudo wfns
ENDDO
DEALLOCATE (vnl, aux)
!
! redetermine ikk2
!
DO iv=1,nbeta
ikk2(iv)=kkbeta
DO ir = kkbeta,1,-1
IF ( abs(betar(ir,iv)) > 1.d-12 ) THEN
ikk2(iv)=ir
exit
ENDIF
ENDDO
ENDDO
ENDIF
ALLOCATE (rho_at(mesh))
rho_at = rho_at_
DEALLOCATE (rho_at_)
ALLOCATE (chi(mesh,ntwfc))
chi = chi_
DEALLOCATE (chi_)
RETURN
END SUBROUTINE convert_ncpp
Reference in New Issue View Git Blame Copy Permalink