diff --git a/CPV/pseudopot.f90 b/CPV/pseudopot.f90 index c63fd348d..2460792af 100644 --- a/CPV/pseudopot.f90 +++ b/CPV/pseudopot.f90 @@ -465,6 +465,7 @@ omega = ht%deth + ps%wsg = 0.0d0 do is = 1, size(ps%wsg,2) do igh = 1, nh( is ) ps%wsg( igh, is) = 4.0d0 * ( 4.0d0 * pi ) ** 2 * dvan( igh, igh, is ) / omega @@ -509,10 +510,17 @@ ALLOCATE( vloc( ps%ap(is)%mesh ) ) ! vloc = ps%ap(is)%vloc * 2.0d0 + WRITE(6,*) ' DEBUG vloc = ', SUM( vloc( : ) ) CALL formfn( vps(:,is), dvps(:,is), ps%ap(is)%rw, ps%ap(is)%rab, & vloc, ps%ap(is)%zv, rcmax(is), g, omega, & tpiba2, 0.0d0, ps%ap(is)%mesh, ngm, .false., .true. ) dvps(:,is) = -dvps(:,is) + WRITE(6,*) ' DEBUG mesh = ', ps%ap(is)%mesh + WRITE(6,*) ' DEBUG rcma = ', rcmax(is) + WRITE(6,*) ' DEBUG zv = ', ps%ap(is)%zv + WRITE(6,*) ' DEBUG rw = ', SUM( ps%ap(is)%rw( : ) ) + WRITE(6,*) ' DEBUG rab = ', SUM( ps%ap(is)%rab( : ) ) + WRITE(6,*) ' DEBUG dvps = ', SUM( dvps( :, is ) ) ! DEALLOCATE( vloc ) diff --git a/Modules/constants.f90 b/Modules/constants.f90 index dae93f019..6b64b392e 100644 --- a/Modules/constants.f90 +++ b/Modules/constants.f90 @@ -82,6 +82,7 @@ MODULE constants ! REAL(dbl), PARAMETER :: eps4 = 1.0D-4 REAL(dbl), PARAMETER :: eps8 = 1.0D-8 + REAL(dbl), PARAMETER :: eps14 = 1.0D-14 REAL(dbl), PARAMETER :: eps16 = 1.0D-16 REAL(dbl), PARAMETER :: eps32 = 1.0D-32 ! diff --git a/flib/sph_bes.f90 b/flib/sph_bes.f90 index 40942fd8e..3750222d5 100644 --- a/flib/sph_bes.f90 +++ b/flib/sph_bes.f90 @@ -19,7 +19,7 @@ subroutine sph_bes (msh, r, q, l, jl) ! ... jl(1:msh) = j_l(q*r(i)) (j_l = spherical bessel function) ! use kinds, only: DP - USE constants, ONLY : eps8 + USE constants, ONLY : eps14 ! implicit none ! @@ -32,7 +32,7 @@ subroutine sph_bes (msh, r, q, l, jl) real(kind=DP) :: qr(msh), sin_qr(msh), cos_qr(msh) #endif ! - if (abs (q) < eps8) then + if (abs (q) < eps14) then if (l == -1) then call errore ('sph_bes', 'j_{-1}(0) ?!?', 1) elseif (l == 0) then @@ -41,7 +41,7 @@ subroutine sph_bes (msh, r, q, l, jl) jl(:) = 0.d0 endif else - if (abs (q * r (1) ) > eps8) then + if (abs (q * r (1) ) > eps14) then ir0 = 1 else if (l == -1) then diff --git a/upftools/fpmd2upf.f90 b/upftools/fpmd2upf.f90 index bc71fe521..17ff8a436 100644 --- a/upftools/fpmd2upf.f90 +++ b/upftools/fpmd2upf.f90 @@ -5,6 +5,30 @@ ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! + +! This utility can be used to convert norm-conserving +! pseudopotentials from FPMD old format to UPF format +! +! Usage: +! fpmd2upf.x < input.namelist +! +! input.namelist should contain the namelist fpmd_pseudo +! fpmd_pseudo parameter are: +! +! psfile pseudopotential filename in FPMD format +! nwfs Number of wavefunction +! wfl(i) i = 1, nwfs Wavefunction label +! wfoc(i) i = 1, nwfs Wavefunction occupation +! psd element name +! zp valence charge +! iexch exchange functional +! icorr correlation functional +! igcx exchange gradient correction +! igcc correlation gradient correction +! +! Example: + + module fpmd2upf_module USE kinds, ONLY: dbl @@ -109,6 +133,7 @@ contains DO ir = 1, mesh x = xmin + REAL(ir-1) * dx ap%rw(ir) = EXP(x) / zmesh + IF( ap%rw(ir) > 1000.0d0 ) EXIT END DO ap%mesh = mesh ap%dx = dx @@ -758,19 +783,22 @@ program fpmd2upf TYPE (pseudo_ncpp) :: ap CHARACTER(LEN=256) :: psfile + CHARACTER(LEN=2) :: wfl( 10 ) + REAL(kind=8) :: wfoc( 10 ) INTEGER :: nsp, nspnl, i, lloc, l, ir, iv, kkbeta REAL(kind=8) :: rmax = 10 REAL(kind=8) :: vll REAL(kind=8), allocatable :: aux(:) -! ... end of declarations -! ---------------------------------------------- + namelist / fpmd_pseudo / psfile, nwfs, wfl, wfoc, psd, & + iexch, icorr, igcx, igcc, zp - WRITE(6,*) ' Enter the pseudopotential filename in FPMD format : ' - READ(5,'(a)') psfile + ! ... end of declarations - nsp = 1 - CALL read_pseudo_fpmd(ap, psfile) + 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")') @@ -778,16 +806,15 @@ program fpmd2upf rcloc = 0.0 - write(6, * ) 'Number of wavefunction > ' - read (5,*) 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 - write(6, * ) 'Wavefunction ',i,' label, occupation > ' - read (5,*) els(i), oc(i) + els(i) = wfl(i) + oc(i) = wfoc(i) lchi(i) = i - 1 nns (i) = 0 rcut(i) = 0.0 @@ -795,19 +822,14 @@ program fpmd2upf epseu(i) = 0.0 end do - write(6, * ) 'element > ' - read(5,*) psd - pseudotype = 'NC' nlcc = ap%tnlcc - zp = ap%zv + if( ap%zv > 0.0d0 ) zp = ap%zv etotps = 0.0 lloc = ap%lloc lmax = MAX( MAXVAL( ap%lll( 1:ap%nbeta ) ), ap%lloc - 1 ) -! go to 100 - nbeta = ap%nbeta mesh = ap%mesh ntwfc = nwfs @@ -818,9 +840,6 @@ program fpmd2upf elsw(i) = els(i) end do - write(6, * ) 'XC (iexch,icorr,igcx,igcc) > ' - read(5,*) iexch, icorr, igcx, igcc - allocate(rab(mesh)) allocate( r(mesh)) r = ap%rw @@ -829,7 +848,6 @@ program fpmd2upf write(6,*) ap%lloc, ap%lll( 1:ap%nbeta ) , ap%nbeta, ap%dx - allocate (rho_atc(mesh)) if (nlcc) rho_atc = ap%rhoc @@ -893,6 +911,8 @@ program fpmd2upf end program fpmd2upf + + !---------------------------------------------------------------------- subroutine simpson2(mesh,func,rab,asum) !-----------------------------------------------------------------------