mirror of https://gitlab.com/QEF/q-e.git
375 lines
11 KiB
Fortran
375 lines
11 KiB
Fortran
!
|
|
! Copyright (C) 2001-2007 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 .
|
|
!
|
|
!
|
|
!-----------------------------------------------------------------------
|
|
subroutine read_ncpp (iunps, upf, ierr)
|
|
!-----------------------------------------------------------------------
|
|
!
|
|
USE upf_kinds, only: dp
|
|
USE upf_params, ONLY: lmaxx
|
|
USE upf_io , ONLY: stdout
|
|
USE pseudo_types
|
|
|
|
implicit none
|
|
!
|
|
TYPE (pseudo_upf) :: upf
|
|
integer, intent(in) :: iunps
|
|
integer, intent(out):: ierr
|
|
!
|
|
real(DP) :: cc(2), alpc(2), aps(6,0:3), alps(3,0:3), &
|
|
a_nlcc, b_nlcc, alpha_nlcc
|
|
real(DP) :: x, vll
|
|
real(DP), allocatable:: vnl(:,:)
|
|
real(DP), parameter :: rcut = 10.d0, e2 = 2.d0
|
|
integer :: nbeta, nlc, nnl, lmax, lloc, ll(1)
|
|
integer :: nb, i, l, ir, ios = 0
|
|
logical :: bhstype, numeric
|
|
!
|
|
!====================================================================
|
|
! read norm-conserving PPs
|
|
!
|
|
ierr = 1
|
|
!
|
|
read (iunps, *, end=300, err=300, iostat=ios) upf%dft
|
|
if (upf%dft(1:2) .eq.'**') upf%dft = 'PZ'
|
|
read (iunps, *, err=300, iostat=ios) upf%psd, upf%zp, lmax, nlc, &
|
|
nnl, upf%nlcc, lloc, bhstype
|
|
if (nlc > 2 .or. nnl > 3 .or. nlc*nnl < 0 ) then
|
|
write(stdout,'(5x,"read_ncpp: Wrong nlc and/or nnl")')
|
|
return
|
|
end if
|
|
if (upf%zp <= 0d0 .or. upf%zp > 100 ) then
|
|
write(stdout,'(5x,"read_ncpp: Wrong zp")')
|
|
return
|
|
end if
|
|
!
|
|
! In numeric pseudopotentials both nlc and nnl are zero.
|
|
!
|
|
numeric = (nlc <= 0) .and. (nnl <= 0)
|
|
!
|
|
if (lloc == -1000) lloc = lmax
|
|
if (lloc < 0 .or. lmax < 0 .or. &
|
|
.not.numeric .and. (lloc > min(lmax+1,lmaxx+1) .or. &
|
|
lmax > max(lmaxx,lloc)) .or. &
|
|
numeric .and. (lloc > lmax .or. lmax > lmaxx) ) then
|
|
write(stdout,'(5x,"read_ncpp: Wrong lmax and/or lloc")')
|
|
return
|
|
end if
|
|
if (.not.numeric ) then
|
|
!
|
|
! read here pseudopotentials in analytic form
|
|
!
|
|
read (iunps, *, err=300, iostat=ios) &
|
|
(alpc(i), i=1,2), (cc(i), i=1,2)
|
|
if (abs (cc(1)+cc(2)-1.d0) > 1.0d-6) then
|
|
write(stdout,'(5x,"read_ncpp: Wrong PP coefficients")')
|
|
return
|
|
end if
|
|
do l = 0, lmax
|
|
read (iunps, *, err=300, iostat=ios) (alps(i,l), i=1,3), &
|
|
(aps(i,l), i=1,6)
|
|
enddo
|
|
if (upf%nlcc ) then
|
|
read (iunps, *, err=300, iostat=ios) &
|
|
a_nlcc, b_nlcc, alpha_nlcc
|
|
if (alpha_nlcc <= 0.d0) then
|
|
write(stdout,'(5x,"read_ncpp: wrong alpha_nlcc")')
|
|
return
|
|
end if
|
|
endif
|
|
endif
|
|
read (iunps, *, err=300, iostat=ios) upf%zmesh, upf%xmin, upf%dx, &
|
|
upf%mesh, upf%nwfc
|
|
if ( upf%mesh <= 0) then
|
|
write(stdout,'(5x,"read_ncpp: wrong number of mesh points")')
|
|
return
|
|
end if
|
|
if ( upf%nwfc < 0 .or. &
|
|
(upf%nwfc < lmax .and. lloc == lmax) .or. &
|
|
(upf%nwfc < lmax+1 .and. lloc /= lmax) ) then
|
|
write(stdout,'(5x,"read_ncpp: wrong number of wavefunctions")')
|
|
return
|
|
end if
|
|
!
|
|
! Here pseudopotentials in numeric form are read
|
|
!
|
|
ALLOCATE ( upf%chi(upf%mesh,upf%nwfc), upf%rho_atc(upf%mesh) )
|
|
upf%rho_atc(:) = 0.d0
|
|
ALLOCATE ( upf%lchi(upf%nwfc), upf%oc(upf%nwfc) )
|
|
allocate (vnl(upf%mesh, 0:lmax))
|
|
if (numeric ) then
|
|
do l = 0, lmax
|
|
read (iunps, '(a)', err=300, iostat=ios)
|
|
read (iunps, *, err=300, iostat=ios) (vnl(ir,l), ir=1,upf%mesh )
|
|
enddo
|
|
if ( upf%nlcc ) then
|
|
read (iunps, *, err=300, iostat=ios) (upf%rho_atc(ir), ir=1,upf%mesh)
|
|
endif
|
|
endif
|
|
!
|
|
! Here pseudowavefunctions (in numeric form) are read
|
|
!
|
|
do nb = 1, upf%nwfc
|
|
read (iunps, '(a)', err=300, iostat=ios)
|
|
read (iunps, *, err=300, iostat=ios) upf%lchi(nb), upf%oc(nb)
|
|
!
|
|
! Test lchi and occupation numbers
|
|
!
|
|
if (nb <= lmax .and. upf%lchi(nb)+1 /= nb) then
|
|
write(stdout,'(5x,"read_ncpp: wrong order of wavefunctions")')
|
|
return
|
|
end if
|
|
if (upf%lchi(nb) > lmaxx .or. upf%lchi(nb) < 0) then
|
|
write(stdout,'(5x,"read_ncpp: wrong lchi")')
|
|
return
|
|
end if
|
|
if (upf%oc(nb) < 0.d0 .or. upf%oc(nb) > 2.d0*(2*upf%lchi(nb)+1)) then
|
|
write(stdout,'(5x,"read_ncpp: wrong oc")')
|
|
return
|
|
end if
|
|
read (iunps, *, err=300, iostat=ios) ( upf%chi(ir,nb), ir=1,upf%mesh )
|
|
enddo
|
|
!
|
|
ierr = 0
|
|
!====================================================================
|
|
! PP succesfully read: now setup
|
|
!
|
|
IF ( numeric ) THEN
|
|
upf%generated='Generated by old ld1 code (numerical format)'
|
|
ELSE
|
|
upf%generated='From published tables, or generated by old fitcar code (analytical format)'
|
|
END IF
|
|
!
|
|
! calculate the number of beta functions
|
|
!
|
|
upf%nbeta = 0
|
|
do l = 0, lmax
|
|
if (l /= lloc ) upf%nbeta = upf%nbeta + 1
|
|
enddo
|
|
!
|
|
! FIXME: to avoid problems when broadcasting PPs with 0 nonlocal betas,
|
|
! use "nbeta" instead of "upf%nbeta" for allocations
|
|
!
|
|
nbeta = MAX(1,upf%nbeta)
|
|
ALLOCATE ( upf%lll( nbeta) )
|
|
nb = 0
|
|
do l = 0, lmax
|
|
if (l /= lloc ) then
|
|
nb = nb + 1
|
|
upf%lll (nb) = l
|
|
end if
|
|
enddo
|
|
!
|
|
! compute the radial mesh
|
|
!
|
|
ALLOCATE ( upf%r(upf%mesh), upf%rab(upf%mesh) )
|
|
do ir = 1, upf%mesh
|
|
x = upf%xmin + DBLE (ir - 1) * upf%dx
|
|
upf%r(ir) = exp (x) / upf%zmesh
|
|
upf%rab(ir) = upf%dx * upf%r(ir)
|
|
enddo
|
|
do ir = 1, upf%mesh
|
|
if ( upf%r(ir) > rcut) then
|
|
upf%kkbeta = ir
|
|
go to 5
|
|
end if
|
|
end do
|
|
upf%kkbeta = upf%mesh
|
|
!
|
|
! ... force kkbeta to be odd for simpson integration (obsolete?)
|
|
!
|
|
5 upf%kkbeta = 2 * ( ( upf%kkbeta + 1 ) / 2) - 1
|
|
!
|
|
ALLOCATE ( upf%kbeta(nbeta) )
|
|
upf%kbeta(:) = upf%kkbeta
|
|
ALLOCATE ( upf%vloc(upf%mesh) )
|
|
upf%vloc (:) = 0.d0
|
|
!
|
|
if (.not. numeric) then
|
|
!
|
|
! bring analytic potentials into numerical form
|
|
!
|
|
IF ( nlc == 2 .AND. nnl == 3 .AND. bhstype ) THEN
|
|
ll(1) = lmax ! workaround for NAG compiler
|
|
CALL bachel( alps(1,0), aps(1,0), 1, ll )
|
|
END IF
|
|
!
|
|
do i = 1, nlc
|
|
do ir = 1, upf%kkbeta
|
|
upf%vloc (ir) = upf%vloc (ir) - upf%zp * e2 * cc (i) * &
|
|
erf ( sqrt (alpc(i)) * upf%r(ir) ) / upf%r(ir)
|
|
end do
|
|
end do
|
|
do l = 0, lmax
|
|
vnl (:, l) = upf%vloc (1:upf%mesh)
|
|
do i = 1, nnl
|
|
vnl (:, l) = vnl (:, l) + e2 * (aps (i, l) + &
|
|
aps (i + 3, l) * upf%r (:) **2) * &
|
|
exp ( - upf%r(:) **2 * alps (i, l) )
|
|
enddo
|
|
enddo
|
|
if ( upf%nlcc ) then
|
|
upf%rho_atc(:) = ( a_nlcc + b_nlcc*upf%r(:)**2 ) * &
|
|
exp ( -upf%r(:)**2 * alpha_nlcc )
|
|
end if
|
|
!
|
|
end if
|
|
!
|
|
! assume l=lloc as local part and subtract from the other channels
|
|
!
|
|
if (lloc <= lmax ) &
|
|
upf%vloc (:) = vnl (:, lloc)
|
|
! lloc > lmax is allowed for PP in analytical form only
|
|
! it means that only the erf part is taken as local part
|
|
do l = 0, lmax
|
|
if (l /= lloc) vnl (:, l) = vnl(:, l) - upf%vloc(:)
|
|
enddo
|
|
!
|
|
! compute the atomic charges
|
|
!
|
|
ALLOCATE ( upf%rho_at (upf%mesh) )
|
|
upf%rho_at(:) = 0.d0
|
|
do nb = 1, upf%nwfc
|
|
if ( upf%oc(nb) > 0.d0) then
|
|
do ir = 1, upf%mesh
|
|
upf%rho_at(ir) = upf%rho_at(ir) + upf%oc(nb) * upf%chi(ir,nb)**2
|
|
enddo
|
|
endif
|
|
enddo
|
|
!====================================================================
|
|
! convert to separable (KB) form
|
|
!
|
|
ALLOCATE ( upf%beta (upf%mesh, nbeta) )
|
|
ALLOCATE ( upf%dion (nbeta,nbeta) )
|
|
upf%dion (:,:) = 0.d0
|
|
nb = 0
|
|
do l = 0, lmax
|
|
if (l /= lloc ) then
|
|
nb = nb + 1
|
|
! upf%beta is used here as work space
|
|
do ir = 1, upf%kkbeta
|
|
upf%beta (ir, nb) = upf%chi(ir, l+1) **2 * vnl(ir, l)
|
|
end do
|
|
call simpson (upf%kkbeta, upf%beta (1, nb), upf%rab, vll )
|
|
upf%dion (nb, nb) = 1.d0 / vll
|
|
! upf%beta stores projectors |beta(r)> = |V_nl(r)phi(r)>
|
|
do ir = 1, upf%kkbeta
|
|
upf%beta (ir, nb) = vnl (ir, l) * upf%chi (ir, l + 1)
|
|
enddo
|
|
upf%lll (nb) = l
|
|
endif
|
|
enddo
|
|
deallocate (vnl)
|
|
!
|
|
! for compatibility with USPP and other formats
|
|
!
|
|
upf%nqf = 0
|
|
upf%nqlc= 0
|
|
upf%tvanp =.false.
|
|
upf%tpawp =.false.
|
|
upf%has_so=.false.
|
|
upf%has_wfc=.false.
|
|
upf%has_gipaw=.false.
|
|
upf%tcoulombp=.false.
|
|
upf%is_gth=.false.
|
|
upf%is_multiproj=.false.
|
|
!
|
|
! Set additional, not present, variables to dummy values
|
|
allocate(upf%els(upf%nwfc), upf%nchi(upf%nwfc), upf%epseu(upf%nwfc))
|
|
upf%els(:) = 'nX'
|
|
upf%nchi(:) = 0
|
|
upf%epseu(:)= 0._dp
|
|
allocate(upf%rcut_chi(upf%nwfc), upf%rcutus_chi(upf%nwfc))
|
|
upf%rcut_chi(:) = 0._dp
|
|
upf%rcutus_chi(:) = 0._dp
|
|
allocate(upf%els_beta( nbeta) )
|
|
upf%els_beta(:) = 'nX'
|
|
allocate(upf%rcut(nbeta), upf%rcutus(nbeta))
|
|
upf%rcut(:) = 0._dp
|
|
upf%rcutus(:) = 0._dp
|
|
!
|
|
return
|
|
|
|
300 write(stdout,'(5x,"read_ncpp: PP file is empty or wrong")')
|
|
end subroutine read_ncpp
|
|
!
|
|
!----------------------------------------------------------------------
|
|
subroutine bachel (alps, aps, npseu, lmax)
|
|
!----------------------------------------------------------------------
|
|
!
|
|
USE upf_kinds, ONLY : DP
|
|
USE upf_const, ONLY : pi
|
|
implicit none
|
|
!
|
|
! First I/O variables
|
|
!
|
|
integer :: npseu, lmax (npseu)
|
|
! input: number of pseudopotential
|
|
! input: max. angul. momentum of the ps
|
|
real(DP) :: alps (3, 0:3, npseu), aps (6, 0:3, npseu)
|
|
! input: the b_l coefficient
|
|
! in/out: the a_l coefficient
|
|
!
|
|
! Here local variables
|
|
!
|
|
integer :: np, lmx, l, i, j, k, ia, ka, nik
|
|
! counter on number of pseudopot.
|
|
! aux. var. (max. ang. mom. of a fix. ps
|
|
! counter on angular momentum
|
|
|
|
real(DP) :: s (6, 6), alpl, alpi, ail
|
|
! auxiliary array
|
|
! first real aux. var. (fix. value of al
|
|
! second real aux. var. (fix. value of a
|
|
! third real aux. var.
|
|
!
|
|
do np = 1, npseu
|
|
lmx = lmax (np)
|
|
do l = 0, lmx
|
|
do k = 1, 6
|
|
ka = mod (k - 1, 3) + 1
|
|
alpl = alps (ka, l, np)
|
|
do i = 1, k
|
|
ia = mod (i - 1, 3) + 1
|
|
alpi = alps (ia, l, np)
|
|
ail = alpi + alpl
|
|
s (i, k) = sqrt (pi / ail) / 4.d0 / ail
|
|
nik = int ( (k - 1) / 3) + int ( (i - 1) / 3) + 1
|
|
do j = 2, nik
|
|
s (i, k) = s (i, k) / 2.d0 / ail * (2 * j - 1)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
!
|
|
do i = 1, 6
|
|
do j = i, 6
|
|
do k = 1, i - 1
|
|
s (i, j) = s (i, j) - s (k, i) * s (k, j)
|
|
enddo
|
|
if (i.eq.j) then
|
|
s (i, i) = sqrt (s (i, i) )
|
|
else
|
|
s (i, j) = s (i, j) / s (i, i)
|
|
endif
|
|
enddo
|
|
enddo
|
|
!
|
|
aps (6, l, np) = - aps (6, l, np) / s (6, 6)
|
|
do i = 5, 1, - 1
|
|
aps (i, l, np) = - aps (i, l, np)
|
|
do k = i + 1, 6
|
|
aps (i, l, np) = aps (i, l, np) - aps (k, l, np) * s (i, k)
|
|
enddo
|
|
aps (i, l, np) = aps (i, l, np) / s (i, i)
|
|
enddo
|
|
enddo
|
|
|
|
enddo
|
|
return
|
|
end subroutine bachel
|