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

There can be conflicts between the erf and erfc in QE and those provided 

by external libraries (e.g. recent ESSL). In order to prevent such problems,
erf has been renamed qe_erf and erfc qe_erfc


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@5644 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2009-07-08 10:29:32 +00:00
parent b3daf98bf0
commit fd7a11d0dc
29 changed files with 112 additions and 128 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CPV/adjef.f90
View File

@ -118,9 +118,9 @@
USE kinds
IMPLICIT NONE
REAL(DP) :: x
REAL(DP), EXTERNAL :: erfc
REAL(DP), PARAMETER :: c=0.5641895835D0
! stepf=erfc(x)
! REAL(DP), EXTERNAL :: qe_erfc
! stepf=qe_erfc(x)
stepf=1.d0/(exp(min(x,100.d0))+1.d0)
END FUNCTION stepf

5
CPV/cp_fpmd.f90
View File

@ -970,7 +970,7 @@ END SUBROUTINE gshcount
implicit none
!
REAL(DP) :: alat, b1_(3),b2_(3),b3_(3), gmax
REAL(DP), external :: erf
REAL(DP), external :: qe_erf
REAL(DP) :: b1(3),b2(3),b3(3), tpiba2, gcutz
!
integer i1,i2,i3,ig
@ -998,7 +998,8 @@ END SUBROUTINE gshcount
!
IF( gcutz > 0.0d0 ) THEN
do ig=1,ngw
ggp(ig) = g(ig) + gcutz * ( 1.0d0 + erf( ( tpiba2 * g(ig) - ecfix ) / ecsig ) )
ggp(ig) = g(ig) + gcutz * &
( 1.0d0 + qe_erf( ( tpiba2 * g(ig) - ecfix ) / ecsig ) )
enddo
ELSE
ggp( 1 : ngw ) = g( 1 : ngw )

24
CPV/efermi.f90
View File

@ -167,7 +167,7 @@ SUBROUTINE EFERMI(NEL,NBANDS,DEL,NKPTS,OCC,EF,EIGVAL, &
REAL(kind=DP) :: weight(nkpts), sort(nbands*nkpts)
REAL(kind=DP), EXTERNAL :: ERFC,FERMID,DELTHM,POSHM,POSHM2, SPLINE
REAL(kind=DP), EXTERNAL :: qe_erfc,FERMID,DELTHM,POSHM,POSHM2, SPLINE
INTEGER, PARAMETER :: JMAX =300
REAL(kind=DP), PARAMETER :: XACC=1.0D-17
@ -387,7 +387,7 @@ SUBROUTINE EFERMI(NEL,NBANDS,DEL,NKPTS,OCC,EF,EIGVAL, &
DO J = 1,NBANDS
X = (XE2 - EIGVAL(J,ISPPT))/DEL
IF(ISMEAR.EQ.1) THEN
Z1 = Z1 + WEIGHT(ISPPT)*( 2.d0 - ERFC(X) )/fspin
Z1 = Z1 + WEIGHT(ISPPT)*( 2.d0 - qe_erfc(X) )/fspin
ELSEIF(ISMEAR.EQ.2) THEN
Z1 = Z1 + WEIGHT(ISPPT)*FERMID(-X)/fspin
ELSEIF(ISMEAR.EQ.3) THEN
@ -413,7 +413,7 @@ SUBROUTINE EFERMI(NEL,NBANDS,DEL,NKPTS,OCC,EF,EIGVAL, &
DO J = 1,NBANDS
X = (XE1 - EIGVAL(J,ISPPT))/DEL
IF(ISMEAR.EQ.1) THEN
Z1 = Z1 + WEIGHT(ISPPT)*( 2.d0 - ERFC(X) )/fspin
Z1 = Z1 + WEIGHT(ISPPT)*( 2.d0 - qe_erfc(X) )/fspin
ELSEIF(ISMEAR.EQ.2) THEN
Z1 = Z1 + WEIGHT(ISPPT)*FERMID(-X)/fspin
ELSEIF(ISMEAR.EQ.3) THEN
@ -454,7 +454,7 @@ SUBROUTINE EFERMI(NEL,NBANDS,DEL,NKPTS,OCC,EF,EIGVAL, &
DO J2 = 1,NBANDS
X = (XMID - EIGVAL(J2,ISPPT))/DEL
IF(ISMEAR.EQ.1) THEN
Z1 = Z1 + WEIGHT(ISPPT)*( 2.d0 - ERFC(X) )/fspin
Z1 = Z1 + WEIGHT(ISPPT)*( 2.d0 - qe_erfc(X) )/fspin
ELSEIF(ISMEAR.EQ.2) THEN
Z1 = Z1 + WEIGHT(ISPPT)*FERMID(-X)/fspin
ELSEIF(ISMEAR.EQ.3) THEN
@ -488,7 +488,7 @@ SUBROUTINE EFERMI(NEL,NBANDS,DEL,NKPTS,OCC,EF,EIGVAL, &
DO J = 1,NBANDS
X = ( EF-EIGVAL(J,ISPPT))/DEL
IF(ISMEAR.EQ.1) THEN
OCC(J,ISPPT) = 2.d0 - ERFC(X)
OCC(J,ISPPT) = 2.d0 - qe_erfc(X)
ELSEIF(ISMEAR.EQ.2) THEN
OCC(J,ISPPT) = FERMID(-X)
ELSEIF(ISMEAR.EQ.3) THEN
@ -530,7 +530,7 @@ SUBROUTINE EFERMI(NEL,NBANDS,DEL,NKPTS,OCC,EF,EIGVAL, &
& *(2.0d0*x*x-1.d0)*exp(-x*x)/(2.0d0*sqrt(pi))
ELSEIF(ISMEAR.EQ.4) THEN
x=abs(x)
zeta=eesh*abs(x)*exp(-(x+sq2i)**2)+piesqq*erfc(x+sq2i)
zeta=eesh*abs(x)*exp(-(x+sq2i)**2)+piesqq*qe_erfc(x+sq2i)
delcor=delcor-del*WEIGHT(ISPPT)*zeta
ELSEIF(ISMEAR.EQ.5) THEN
a=-0.5634d0
@ -629,7 +629,7 @@ FUNCTION delthm(xx)
REAL(kind=DP) :: delthm
REAL(kind=DP), INTENT(in) :: xx
REAL(kind=DP), EXTERNAL :: ERFC
REAL(kind=DP), EXTERNAL :: qe_erfc
REAL(kind=DP) :: pi
@ -639,7 +639,7 @@ FUNCTION delthm(xx)
ELSEIF(XX .LT. -10.D0) THEN
DELTHM=0.D0
ELSE
DELTHM=(2.D0-ERFC(XX))+XX*EXP(-XX*XX)/SQRT(PI)
DELTHM=(2.D0-qe_erfc(XX))+XX*EXP(-XX*XX)/SQRT(PI)
ENDIF
!
RETURN
@ -680,7 +680,7 @@ FUNCTION poshm(x)
REAL(kind=DP) :: poshm
REAL(kind=DP), INTENT(in) :: x
REAL(kind=DP), EXTERNAL :: ERFC
REAL(kind=DP), EXTERNAL :: qe_erfc
REAL(kind=DP) :: pi,a
@ -692,7 +692,7 @@ FUNCTION poshm(x)
ELSEIF(X .LT. -10.D0) THEN
POSHM=0.D0
ELSE
POSHM=(2.D0-ERFC(X))+(-2.d0*a*x*x+2*x+a)*EXP(-X*X)/SQRT(PI)/2.d0
POSHM=(2.D0-qe_erfc(X))+(-2.d0*a*x*x+2*x+a)*EXP(-X*X)/SQRT(PI)/2.d0
ENDIF
!
RETURN
@ -710,7 +710,7 @@ FUNCTION poshm2(x)
REAL(kind=DP) :: poshm2
REAL(kind=DP), INTENT(in) :: x
REAL(kind=DP), EXTERNAL :: ERFC
REAL(kind=DP), EXTERNAL :: qe_erfc
REAL(kind=DP) :: pi
@ -720,7 +720,7 @@ FUNCTION poshm2(x)
ELSEIF(X .LT. -10.D0) THEN
POSHM2=0.D0
ELSE
POSHM2=(2.D0-ERFC(X-1.d0/sqrt(2.d0)))+ &
POSHM2=(2.D0-qe_erfc(X-1.d0/sqrt(2.d0)))+ &
& sqrt(2.d0)*exp(-x*x+sqrt(2.d0)*x-0.5d0)/sqrt(pi)
ENDIF
!

8
CPV/potentials.f90
View File

@ -674,7 +674,7 @@
REAL(DP), INTENT(IN) :: omega, hmat( 3, 3 )
COMPLEX(DP) :: screen_coul( ngm )
REAL(DP), EXTERNAL :: erf
REAL(DP), EXTERNAL :: qe_erf
! ... Locals
!
@ -715,7 +715,7 @@
IF( rmod < gsmall ) THEN
grr( ir ) = fact * 2.0d0 * rc / SQRT( pi )
ELSE
grr( ir ) = fact * erf( rc * rmod ) / rmod
grr( ir ) = fact * qe_erf( rc * rmod ) / rmod
END IF
END DO
END DO
@ -1055,7 +1055,7 @@
LOGICAL, INTENT(IN) :: TSTRESS
REAL(DP), INTENT(in) :: hmat( 3, 3 )
REAL(DP), EXTERNAL :: erfc
REAL(DP), EXTERNAL :: qe_erfc
INTEGER :: ldim_block, gind_block
EXTERNAL ldim_block, gind_block
@ -1208,7 +1208,7 @@
ELSE
ESRTZERO=1.D0
END IF
ADDESR = ZV2_KJ * erfc(ARG) / RLM
ADDESR = ZV2_KJ * qe_erfc(ARG) / RLM
ESR = ESR + ESRTZERO*ADDESR
ADDPRE = FACT_PRE * EXP(-ARG*ARG)
REPAND = ESRTZERO*(ADDESR + ADDPRE)/ERRE2

6
CPV/pseudo_base.f90
View File

@ -209,7 +209,7 @@
integer :: ig, ir, irmax
real(DP), allocatable:: f(:),vscr(:), figl(:)
real(DP), allocatable:: df(:), dfigl(:)
real(DP), external :: erf, erfc
real(DP), external :: qe_erf, qe_erfc
!
allocate( figl(ngs), f(mesh), vscr(mesh) )
if (tpre) then
@ -224,7 +224,7 @@
end do
!
do ir = 1, irmax
vscr(ir) = 0.5d0 * r(ir) * vloc_at(ir) + zv * erf( r(ir) / rcmax )
vscr(ir) = 0.5d0 * r(ir) * vloc_at(ir) + zv * qe_erf( r(ir) / rcmax )
end do
do ir = irmax + 1, mesh
vscr(ir)=0.0d0
@ -240,7 +240,7 @@
! ... reproduce the results from other PW codes
!
DO ir = 1, irmax
f(ir) = fpi * ( zv * erfc( r(ir)/rcmax ) ) * r(ir)
f(ir) = fpi * ( zv * qe_erfc( r(ir)/rcmax ) ) * r(ir)
END DO
DO ir = irmax + 1, mesh
f(ir)=0.0d0

8
D3/d3ionq.f90
View File

@ -72,7 +72,7 @@ subroutine d3ionq (nat, ntyp, ityp, zv, tau, alat, omega, q, at, &
upperbound, charge, fac, gtq2, gt2, facq, d2f, d3f, rmax, &
r (3, mxr), r2 (mxr), dtau (3), rr, ar, qrg
! auxilary variables
real (DP), external :: erfc
real (DP), external :: qe_erfc
complex (DP), allocatable :: d3dy1 (:,:,:), d3dy2 (:,:,:), d3dy3 (:,:,:)
! first term dynamical matrix
@ -106,7 +106,7 @@ subroutine d3ionq (nat, ntyp, ityp, zv, tau, alat, omega, q, at, &
if (alpha == 0.d0) call errore ('d3ion', 'optimal alpha not found', 1)
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * alpha / tpi) * &
erfc (sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
qe_erfc (sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
if (upperbound > 1.d-9) goto 11
@ -248,10 +248,10 @@ subroutine d3ionq (nat, ntyp, ityp, zv, tau, alat, omega, q, at, &
qrg = tpi * (q (1) * (r (1, nr) + dtau (1) ) + &
q (2) * (r (2, nr) + dtau (2) ) + &
q (3) * (r (3, nr) + dtau (3) ) )
d2f = (3.d0 * erfc (ar) + sqrt (8.d0 / tpi) * exp ( - ar**2) &
d2f = (3.d0 * qe_erfc (ar) + sqrt (8.d0 / tpi) * exp ( - ar**2) &
* ar * (3.d0 + 2.d0 * ar**2) ) / rr**5
d3f = ( - 15.d0 * erfc (ar) - sqrt (8.d0 / tpi) * exp ( - ar**2) &
d3f = ( - 15.d0 * qe_erfc (ar) - sqrt (8.d0 / tpi) * exp ( - ar**2) &
* ar * (15.d0 + 10.d0 * ar**2 + 4.d0 * ar**4) ) / rr**7
do jcart = 1, 3
nu_j = (nb - 1) * 3 + jcart

2
EE/setlocalcoul.f90
View File

@ -69,8 +69,6 @@
REAL( DP ) :: adist
REAL( DP ) :: temp
!
REAL( DP ), EXTERNAL :: erf
!
CHARACTER( LEN = 6 ) :: nts
!
INTEGER, EXTERNAL :: COMPINDEX

8
Gamma/d2ion.f90
View File

@ -29,7 +29,7 @@ subroutine d2ion (nat,ntyp,ityp,zv,tau,alat,omega, &
parameter(mxr=50)
real(DP) :: facg(nat), arg, tpiba2, alpha, r(3,mxr), r2(mxr), dtau(3), &
rmax, rr, upperbound, charge, gt2, fac, fnat, df, d2f, ar
real(DP), external:: erfc
real(DP), external:: qe_erfc
!
!
tpiba2 = (tpi/alat)**2
@ -56,7 +56,7 @@ subroutine d2ion (nat,ntyp,ityp,zv,tau,alat,omega, &
! upperbound is a safe upper bound for the error ON THE ENERGY
!
upperbound=e2*charge**2*sqrt(2.0d0*alpha/tpi)* &
& erfc(sqrt(tpiba2*gg(ng)/4.d0/alpha))
& qe_erfc(sqrt(tpiba2*gg(ng)/4.d0/alpha))
if(upperbound.lt.1.0d-6) go to 20
!
gt2 = gg(ng)*tpiba2
@ -131,9 +131,9 @@ subroutine d2ion (nat,ntyp,ityp,zv,tau,alat,omega, &
do nr=1, nrm
rr=sqrt(r2(nr))*alat
ar = sqrt(alpha)*rr
d2f = ( 3.d0*erfc(ar) + sqrt(8.d0/tpi)*ar* &
d2f = ( 3.d0*qe_erfc(ar) + sqrt(8.d0/tpi)*ar* &
(3.d0+2.d0*ar**2)*exp(-ar**2) ) / rr**5
df = ( -erfc(ar) - sqrt(8.d0/tpi)*ar*exp(-ar**2) ) / rr**3
df = ( -qe_erfc(ar) - sqrt(8.d0/tpi)*ar*exp(-ar**2) ) / rr**3
do nu_i = 1,nmodes
if (has_equivalent( (nu_i-1)/3+1 ).eq.1 ) go to 25
arg = r(1,nr)*u(mu_i+1,nu_i) + &

4
Modules/read_ncpp.f90
View File

@ -25,7 +25,7 @@ subroutine read_ncpp (iunps, np, upf)
real(DP) :: x, vll
real(DP), allocatable:: vnl(:,:)
real(DP), parameter :: rcut = 10.d0, e2 = 2.d0
real(DP), external :: erf
real(DP), external :: qe_erf
integer :: nlc, nnl, lmax, lloc
integer :: nb, ios, i, l, ir
logical :: bhstype, numeric
@ -179,7 +179,7 @@ subroutine read_ncpp (iunps, np, upf)
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)
qe_erf ( sqrt (alpc(i)) * upf%r(ir) ) / upf%r(ir)
end do
end do
do l = 0, lmax

8
PH/d2ionq.f90
View File

@ -71,7 +71,7 @@ subroutine d2ionq (nat, ntyp, ityp, zv, tau, alat, omega, q, at, &
complex(DP) :: dy1 (3 * nat, nmodes), dy2 (3 * nat, nmodes), &
dy3 (3 * nat, nmodes), facg, fnat, work
! work spaces, factors
real(DP), external :: erfc
real(DP), external :: qe_erfc
call start_clock ('d2ionq')
@ -92,7 +92,7 @@ subroutine d2ionq (nat, ntyp, ityp, zv, tau, alat, omega, q, at, &
if (alpha == 0.d0) call errore ('d2ion', 'optimal alpha not found',1)
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * alpha / tpi) * &
erfc ( sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
qe_erfc ( sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
if (upperbound > 1.d-9) goto 11
@ -199,9 +199,9 @@ subroutine d2ionq (nat, ntyp, ityp, zv, tau, alat, omega, q, at, &
qrg = tpi * (q (1) * (r (1, nr) + dtau (1) ) + &
q (2) * (r (2, nr) + dtau (2) ) + &
q (3) * (r (3, nr) + dtau (3) ) )
d2f = (3.d0 * erfc (ar) + sqrt (8.d0 / tpi) * ar * &
d2f = (3.d0 * qe_erfc (ar) + sqrt (8.d0 / tpi) * ar * &
(3.d0 + 2.d0 * ar**2) * exp ( - ar**2) ) / rr**5
df = ( - erfc (ar) - sqrt (8.d0 / tpi) * ar * exp ( - ar**2) ) &
df = ( - qe_erfc (ar) - sqrt (8.d0 / tpi) * ar * exp ( - ar**2) ) &
/ rr**3
do icart = 1, 3
na_icart = 3 * (na - 1) + icart

10
PH/hdiag.f90
View File

@ -37,8 +37,6 @@ subroutine hdiag( max_iter, avg_iter, xk_, et_ )
!
REAL(DP) :: cg_iter
! number of iteration in CG
REAL(DP), EXTERNAL :: erf
! error function
INTEGER :: ig, ntry, notconv
! counter on G vectors
! number or repeated call to diagonalization in case of non convergence
@ -64,14 +62,6 @@ subroutine hdiag( max_iter, avg_iter, xk_, et_ )
(xk_ (2) + g (2, igk (ig) ) ) **2 + &
(xk_ (3) + g (3, igk (ig) ) ) **2 ) * tpiba2
enddo
!
! if (qcutz > 0.d0) then
! do ig = 1, npw
! g2kin (ig) = g2kin (ig) + qcutz * (1.d0 + erf ( (g2kin (ig) &
! - ecfixed) / q2sigma) )
! enddo
! endif
!
! Conjugate-Gradient diagonalization
!

4
PH/setlocq.f90
View File

@ -49,7 +49,7 @@ subroutine setlocq (xq, mesh, msh, rab, r, vloc_at, zp, tpiba2, ngm, &
aux1 (mesh), gx
! auxiliary variables
! gx = modulus of g vectors
real(DP), external :: erf
real(DP), external :: qe_erf
! the erf function
integer :: ig, ir
! counters
@ -69,7 +69,7 @@ subroutine setlocq (xq, mesh, msh, rab, r, vloc_at, zp, tpiba2, ngm, &
! indipendent of |G| in real space
!
do ir = 1, msh
aux1 (ir) = r (ir) * vloc_at (ir) + zp * e2 * erf (r (ir) )
aux1 (ir) = r (ir) * vloc_at (ir) + zp * e2 * qe_erf (r (ir) )
enddo
fac = zp * e2 / tpiba2
!

8
PP/do_shift_ew.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! 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,
@ -76,7 +76,7 @@ subroutine do_shift_ew (alat, nat, ntyp, ityp, zv, delta_zv, at, bg, tau, &
! buffer variable
! used to optimize alpha
complex(DP), allocatable :: rhon(:)
real(DP), external :: erfc
real(DP), external :: qe_erfc
allocate (rhon(ngm))
@ -94,7 +94,7 @@ subroutine do_shift_ew (alat, nat, ntyp, ityp, zv, delta_zv, at, bg, tau, &
! upperbound is a safe upper bound for the error in the sum over G
!
if (alpha.le.0.d0) call errore ('do_shift_ew', 'optimal alpha not found', 1)
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * alpha / tpi) * erfc ( &
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * alpha / tpi) * qe_erfc ( &
sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
if (upperbound.gt.1.0d-7) goto 100
!
@ -164,7 +164,7 @@ subroutine do_shift_ew (alat, nat, ntyp, ityp, zv, delta_zv, at, bg, tau, &
rr = sqrt (r2 (nr) ) * alat
shift_ion(na) = shift_ion(na) + &
delta_zv(ityp(na)) * zv (ityp (nb) ) * &
erfc ( sqrt (alpha) * rr) / rr
qe_erfc ( sqrt (alpha) * rr) / rr
enddo
enddo
enddo

4
PW/dvloc_of_g.f90
View File

@ -40,7 +40,7 @@ subroutine dvloc_of_g (mesh, msh, rab, r, vloc_at, zp, tpiba2, ngl, gl, &
!
real(DP) :: vlcp, g2a, gx
real(DP), allocatable :: aux (:), aux1 (:)
real(DP), external :: erf
real(DP), external :: qe_erf
integer :: i, igl, igl0
! counter on erf functions or gaussians
@ -66,7 +66,7 @@ subroutine dvloc_of_g (mesh, msh, rab, r, vloc_at, zp, tpiba2, ngl, gl, &
! indipendent of |G| in real space
!
do i = 1, msh
aux1 (i) = r (i) * vloc_at (i) + zp * e2 * erf (r (i) )
aux1 (i) = r (i) * vloc_at (i) + zp * e2 * qe_erf (r (i) )
enddo
do igl = igl0, ngl
gx = sqrt (gl (igl) * tpiba2)

6
PW/ewald.f90
View File

@ -76,7 +76,7 @@ function ewald (alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, &
! buffer variable
! used to optimize alpha
complex(DP) :: rhon
real(DP), external :: erfc
real(DP), external :: qe_erfc
tpiba2 = (tpi / alat) **2
charge = 0.d0
@ -90,7 +90,7 @@ function ewald (alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, &
! upperbound is a safe upper bound for the error in the sum over G
!
if (alpha.le.0.d0) call errore ('ewald', 'optimal alpha not found', 1)
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * alpha / tpi) * erfc ( &
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * alpha / tpi) * qe_erfc ( &
sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
if (upperbound.gt.1.0d-7) goto 100
!
@ -146,7 +146,7 @@ function ewald (alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, &
!
do nr = 1, nrm
rr = sqrt (r2 (nr) ) * alat
ewaldr = ewaldr + zv (ityp (na) ) * zv (ityp (nb) ) * erfc ( &
ewaldr = ewaldr + zv (ityp (na) ) * zv (ityp (nb) ) * qe_erfc ( &
sqrt (alpha) * rr) / rr
enddo
enddo

8
PW/ewald_dipole.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! 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,
@ -31,7 +31,7 @@ subroutine ewald_dipole (tens,dipole)
real(DP) :: dipole(ntyp),charge, eta, arg, upperbound, temp
complex(DP) :: tens(nat,3,3)
complex(DP) :: rhon
real(DP), external :: erfc
real(DP), external :: qe_erfc
complex(DP), allocatable:: ewaldg(:,:,:), ewaldr(:,:,:)
integer :: alpha, beta, na, ng, nt, ipol, nb, nrm, nr
@ -57,7 +57,7 @@ subroutine ewald_dipole (tens,dipole)
! upperbound is a safe upper bound for the error in the sum over G
!
if (eta.le.0.d0) call errore ('ewald_dipole', 'optimal eta not found', 1)
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * eta / tpi) * erfc ( &
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * eta / tpi) * qe_erfc ( &
sqrt (tpiba2 * gcutm / 4.d0 / eta) )
if (upperbound.gt.1.0d-7) goto 100
!
@ -118,7 +118,7 @@ subroutine ewald_dipole (tens,dipole)
do beta=1,3
rr = sqrt (r2 (nr) ) * alat
temp= dipole (ityp (na)) * ( 3.d0 / rr**3 * &
erfc ( sqrt (eta) * rr) + (6.d0 * sqrt (eta/pi) * &
qe_erfc ( sqrt (eta) * rr) + (6.d0 * sqrt (eta/pi) * &
1.d0 / rr*2 + 4.d0 * sqrt (eta**3/pi))* exp(-eta* rr**2))
ewaldr(na, alpha,beta) = ewaldr(na, alpha,beta)+ temp *&
r(alpha,nr) * r(beta,nr) / rr**2

8
PW/force_ew.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! 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,
@ -78,7 +78,7 @@ subroutine force_ew (alat, nat, ntyp, ityp, zv, at, bg, tau, &
complex(DP), allocatable :: aux (:)
! auxiliary space
real(DP), external :: erfc
real(DP), external :: qe_erfc
!
forceion(:,:) = 0.d0
tpiba2 = (tpi / alat) **2
@ -94,7 +94,7 @@ subroutine force_ew (alat, nat, ntyp, ityp, zv, at, bg, tau, &
10 alpha = alpha - 0.1d0
if (alpha.eq.0.d0) call errore ('force_ew', 'optimal alpha not found', 1)
upperbound = e2 * charge**2 * sqrt (2.d0 * alpha / tpi) * &
erfc ( sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
qe_erfc ( sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
if (upperbound > 1.0d-6) goto 10
!
! G-space sum here
@ -150,7 +150,7 @@ subroutine force_ew (alat, nat, ntyp, ityp, zv, at, bg, tau, &
do n = 1, nrm
rr = sqrt (r2 (n) ) * alat
factor = zv (ityp (na) ) * zv (ityp (nb) ) * e2 / rr**2 * &
(erfc (sqrt (alpha) * rr) / rr + &
(qe_erfc (sqrt (alpha) * rr) / rr + &
sqrt (8.0d0 * alpha / tpi) * exp ( - alpha * rr**2) ) * alat
do ipol = 1, 3
forceion (ipol, na) = forceion (ipol, na) - factor * r (ipol, n)

4
PW/g2_kin.f90
View File

@ -27,7 +27,7 @@ SUBROUTINE g2_kin ( ik )
! ... local variables
!
INTEGER :: ig
REAL(DP), EXTERNAL :: erf
REAL(DP), EXTERNAL :: qe_erf
!
!
g2kin(1:npw) = ( ( xk(1,ik) + g(1,igk(1:npw)) )**2 + &
@ -39,7 +39,7 @@ SUBROUTINE g2_kin ( ik )
DO ig = 1, npw
!
g2kin(ig) = g2kin(ig) + qcutz * &
( 1.D0 + erf( ( g2kin(ig) - ecfixed ) / q2sigma ) )
( 1.D0 + qe_erf( ( g2kin(ig) - ecfixed ) / q2sigma ) )
!
END DO
!

10
PW/stres_ewa.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! 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,
@ -74,7 +74,7 @@ subroutine stres_ewa (alat, nat, ntyp, ityp, zv, at, bg, tau, &
! diagonal term
! nondiagonal term
complex(DP) :: rhostar
real(DP), external :: erfc
real(DP), external :: qe_erfc
! the erfc function
!
tpiba2 = (tpi / alat) **2
@ -91,8 +91,8 @@ subroutine stres_ewa (alat, nat, ntyp, ityp, zv, at, bg, tau, &
12 alpha = alpha - 0.1d0
if (alpha.eq.0.0) call errore ('stres_ew', 'optimal alpha not found &
&', 1)
upperbound = e2 * charge**2 * sqrt (2 * alpha / tpi) * erfc (sqrt &
(tpiba2 * gcutm / 4.0d0 / alpha) )
upperbound = e2 * charge**2 * sqrt (2 * alpha / tpi) * &
qe_erfc ( sqrt (tpiba2 * gcutm / 4.0d0 / alpha) )
if (upperbound.gt.eps6) goto 12
!
! G-space sum here
@ -150,7 +150,7 @@ subroutine stres_ewa (alat, nat, ntyp, ityp, zv, at, bg, tau, &
do nr = 1, nrm
rr = sqrt (r2 (nr) ) * alat
fac = - e2 / 2.0d0 / omega * alat**2 * zv (ityp (na) ) * &
zv ( ityp (nb) ) / rr**3 * (erfc (sqrt (alpha) * rr) + &
zv ( ityp (nb) ) / rr**3 * (qe_erfc (sqrt (alpha) * rr) + &
rr * sqrt (8 * alpha / tpi) * exp ( - alpha * rr**2) )
do l = 1, 3
do m = 1, l

4
PW/vloc_of_g.f90
View File

@ -54,7 +54,7 @@ subroutine vloc_of_g (mesh, msh, rab, r, vloc_at, zp, tpiba2, ngl, &
! igl0:first shell with g != 0
! ir :counter on mesh points
!
real(DP), external :: erf
real(DP), external :: qe_erf
!
allocate ( aux(msh), aux1(msh) )
if (gl (1) < eps8) then
@ -75,7 +75,7 @@ subroutine vloc_of_g (mesh, msh, rab, r, vloc_at, zp, tpiba2, ngl, &
! function independent of |G| in real space
!
do ir = 1, msh
aux1 (ir) = r (ir) * vloc_at (ir) + zp * e2 * erf (r (ir) )
aux1 (ir) = r (ir) * vloc_at (ir) + zp * e2 * qe_erf (r (ir) )
enddo
fac = zp * e2 / tpiba2
!

4
PW/wgauss.f90
View File

@ -41,7 +41,7 @@ function wgauss (x, n)
integer :: i, ni
! counter on the n indices
! counter on 2n
real(DP), external :: gauss_freq, erf
real(DP), external :: gauss_freq, qe_erf
real(DP), parameter :: maxarg = 200.d0
! maximum value for the argument of the exponential
@ -61,7 +61,7 @@ function wgauss (x, n)
if (n.eq. - 1) then
xp = x - 1.0d0 / sqrt (2.0d0)
arg = min (maxarg, xp**2)
wgauss = 0.5d0 * erf (xp) + 1.0d0 / sqrt (2.0d0 * pi) * exp ( - &
wgauss = 0.5d0 * qe_erf (xp) + 1.0d0 / sqrt (2.0d0 * pi) * exp ( - &
arg) + 0.5d0
return

4
VdW/check_v_eff.f90
View File

@ -73,7 +73,7 @@ SUBROUTINE check_v_eff ( veff, charge )
!
! ... external functions
!
REAL(KIND=DP), EXTERNAL :: erf
REAL(KIND=DP), EXTERNAL :: qe_erf
! error function
!
!
@ -173,7 +173,7 @@ SUBROUTINE check_v_eff ( veff, charge )
IF ( qcutz > 0.D0 ) THEN
DO ig = 1, npw
g2kin (ig) = g2kin(ig) + qcutz * &
( 1.D0 + erf( ( g2kin(ig) - ecfixed ) / q2sigma ) )
( 1.D0 + qe_erf( ( g2kin(ig) - ecfixed ) / q2sigma ) )
END DO
END IF
!

6
VdW/eff_pot.f90
View File

@ -50,7 +50,7 @@ subroutine eff_pot (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
integer :: nnn, nite, ite = 0
real (kind=DP) :: vstart, thresh_veff
character (len=10):: str_ite
real(kind=DP), external :: erf
real(kind=DP), external :: qe_erf
!
call start_clock('eff_pot')
!
@ -179,8 +179,8 @@ subroutine eff_pot (rho, nr1, nr2, nr3, nrx1, nrx2, nrx3, nrxx, nl, &
vstart=20000
do ir = 1, nrxx
vrs(ir,1) = v%of_r(ir,1) + vltot(ir)
vv(ir,is) = erf(abs(psi_smooth(1,ir)*dble(vstart)))*vv(ir,is) + &
(1.d0-erf( abs( psi_smooth(1,ir)*dble(vstart) ) ))*&
vv(ir,is) = qe_erf(abs(psi_smooth(1,ir)*dble(vstart)))*vv(ir,is) + &
(1.d0-qe_erf( abs( psi_smooth(1,ir)*dble(vstart) ) ))*&
vrs(ir,is)
enddo

7
atomic/read_pseudo_ncpp.f90
View File

@ -53,7 +53,8 @@ subroutine read_pseudo_ncpp (file_pseudo,zed,grid,ndmx,&
real(DP) :: &
vnloc, a_core, b_core, alfa_core, &
cc(2),alpc(2),alc(6,0:3),alps(3,0:3),erf
cc(2),alpc(2),alc(6,0:3),alps(3,0:3)
real(DP), external :: qe_erf
logical :: &
bhstype, numeric
@ -153,8 +154,8 @@ subroutine read_pseudo_ncpp (file_pseudo,zed,grid,ndmx,&
enddo
do ir=1,mesh
vpsloc(ir) = -2.0_dp*zval/grid%r(ir)* &
( cc(1)*erf(grid%r(ir)*sqrt(alpc(1))) &
+ cc(2)*erf(grid%r(ir)*sqrt(alpc(2))) )
( cc(1)*qe_erf(grid%r(ir)*sqrt(alpc(1))) &
+ cc(2)*qe_erf(grid%r(ir)*sqrt(alpc(2))) )
end do
endif

39
flib/erf.f90
View File

@ -1,13 +1,12 @@
!
! Copyright (C) 2002-2007 Quantum-Espresso group
! Copyright (C) 2002-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 .
!
#include "f_defs.h"
!---------------------------------------------------------------------
function erf (x)
function qe_erf (x)
!---------------------------------------------------------------------
!
! Error function - computed from the rational approximations of
@ -20,8 +19,8 @@ function erf (x)
implicit none
real(DP), intent(in) :: x
real(DP) :: x2, p1 (4), q1 (4)
real(DP), external :: erfc
real(DP) :: erf
real(DP), external :: qe_erfc
real(DP) :: qe_erf
data p1 / 2.426679552305318E2_DP, 2.197926161829415E1_DP, &
6.996383488619136_DP, -3.560984370181538E-2_DP /
data q1 / 2.150588758698612E2_DP, 9.116490540451490E1_DP, &
@ -31,22 +30,22 @@ function erf (x)
!
! erf(6)=1-10^(-17) cannot be distinguished from 1
!
erf = sign (1.0_DP, x)
qe_erf = sign (1.0_DP, x)
else
if (abs (x) <= 0.47_DP) then
x2 = x**2
erf = x * (p1 (1) + x2 * (p1 (2) + x2 * (p1 (3) + x2 * p1 (4) ) ) ) &
qe_erf=x *(p1 (1) + x2 * (p1 (2) + x2 * (p1 (3) + x2 * p1 (4) ) ) ) &
/ (q1 (1) + x2 * (q1 (2) + x2 * (q1 (3) + x2 * q1 (4) ) ) )
else
erf = 1.0_DP - erfc (x)
qe_erf = 1.0_DP - qe_erfc (x)
endif
endif
!
return
end function erf
end function qe_erf
!
!---------------------------------------------------------------------
function erfc (x)
function qe_erfc (x)
!---------------------------------------------------------------------
!
! erfc(x) = 1-erf(x) - See comments in erf
@ -54,9 +53,9 @@ function erfc (x)
use kinds, only : DP
implicit none
real(DP),intent(in) :: x
real(DP) :: erfc
real(DP) :: qe_erfc
real(DP) :: ax, x2, xm2, p2 (8), q2 (8), p3 (5), q3 (5), pim1
real(DP), external :: erf
real(DP), external :: qe_erf
data p2 / 3.004592610201616E2_DP, 4.519189537118719E2_DP, &
3.393208167343437E2_DP, 1.529892850469404E2_DP, &
4.316222722205674E1_DP, 7.211758250883094_DP, &
@ -79,31 +78,31 @@ function erfc (x)
!
! erfc(26.0)=10^(-296); erfc( 9.0)=10^(-37);
!
erfc = 0.0_DP
qe_erfc = 0.0_DP
elseif (ax > 4.0_DP) then
x2 = x**2
xm2 = (1.0_DP / ax) **2
erfc = (1.0_DP / ax) * exp ( - x2) * (pim1 + xm2 * (p3 (1) &
qe_erfc = (1.0_DP / ax) * exp ( - x2) * (pim1 + xm2 * (p3 (1) &
+ xm2 * (p3 (2) + xm2 * (p3 (3) + xm2 * (p3 (4) + xm2 * p3 (5) &
) ) ) ) / (q3 (1) + xm2 * (q3 (2) + xm2 * (q3 (3) + xm2 * &
(q3 (4) + xm2 * q3 (5) ) ) ) ) )
elseif (ax > 0.47_DP) then
x2 = x**2
erfc = exp ( - x2) * (p2 (1) + ax * (p2 (2) + ax * (p2 (3) &
qe_erfc = exp ( - x2) * (p2 (1) + ax * (p2 (2) + ax * (p2 (3) &
+ ax * (p2 (4) + ax * (p2 (5) + ax * (p2 (6) + ax * (p2 (7) &
+ ax * p2 (8) ) ) ) ) ) ) ) / (q2 (1) + ax * (q2 (2) + ax * &
(q2 (3) + ax * (q2 (4) + ax * (q2 (5) + ax * (q2 (6) + ax * &
(q2 (7) + ax * q2 (8) ) ) ) ) ) ) )
else
erfc = 1.0_DP - erf (ax)
qe_erfc = 1.0_DP - qe_erf (ax)
endif
!
! erf(-x)=-erf(x) => erfc(-x) = 2-erfc(x)
!
if (x < 0.0_DP) erfc = 2.0_DP - erfc
if (x < 0.0_DP) qe_erfc = 2.0_DP - qe_erfc
!
return
end function erfc
end function qe_erfc
!
!---------------------------------------------------------------------
function gauss_freq (x)
@ -118,9 +117,9 @@ function gauss_freq (x)
real(DP) :: gauss_freq
real(DP), parameter :: c = 0.7071067811865475_DP
! ( c= sqrt(1/2) )
real(DP), external :: erfc
real(DP), external :: qe_erfc
!
gauss_freq = 0.5_DP * erfc ( - x * c)
gauss_freq = 0.5_DP * qe_erfc ( - x * c)
!
return
end function gauss_freq

5
upftools/casino2upf.f90
View File

@ -109,14 +109,9 @@ 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 (8) :: erf
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(:,:)
EXTERNAL erf
NULLIFY ( mhead, mptr, mtail )
dft_ = 'HF' !Hardcoded at the moment should eventually be HF anyway

10
upftools/fpmd2upf.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! 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,
@ -154,8 +154,7 @@ contains
! REAL(DP) :: pi = 3.14159265358979323846_DP
! ... declare external function
REAL(DP) :: erf, erfc
EXTERNAL erf, erfc
REAL(DP), EXTERNAL :: qe_erf
IF( ap%mesh == 0 ) THEN
! ... Local pseudopotential, define a logaritmic grid
@ -178,8 +177,9 @@ contains
ap%vrps = 0.0d0
do l = 1, 3
do ir = 1, ap%mesh
ap%vnl(ir,l)= - ( ap%wrc(1) * erf(SQRT(ap%rc(1))*ap%rw(ir)) + &
ap%wrc(2) * erf(SQRT(ap%rc(2))*ap%rw(ir)) ) * ap%zv / ap%rw(ir)
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

12
upftools/ncpp2upf.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! 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,
@ -75,10 +75,10 @@ subroutine read_ncpp(iunps)
!
character(len=1), dimension(0:3) :: convel=(/'S','P','D','F'/)
character(len=2) :: label
real (8) :: x, erf
real (8) :: x, qe_erf
integer :: l, i, ir, nb, n
character (len=255) line
external erf
external qe_erf
read(iunps, '(a)', end=300, err=300 ) dft_
if (dft_(1:2).eq.'**') dft_ = 'PZ'
@ -162,11 +162,11 @@ subroutine read_ncpp(iunps)
!
do l=0,lmax_
!
! DO NOT USE f90 ARRAY SYNTAX: erf IS NOT AN INTRINSIC FUNCTION!!!
! DO NOT USE f90 ARRAY SYNTAX: qe_erf IS NOT AN INTRINSIC FUNCTION!!!
!
do ir=1,mesh_
vnl(ir,l)= - ( cc(1)*erf(sqrt(alpc(1))*r_(ir)) + &
cc(2)*erf(sqrt(alpc(2))*r_(ir)) ) * zp_/r_(ir)
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
do n=1,nnl

10
upftools/oldcp2upf.f90
View File

@ -1,5 +1,5 @@
!
! Copyright (C) 2001 PWSCF group
! 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,
@ -61,7 +61,7 @@ subroutine read_oldcp(iunps)
implicit none
integer :: iunps
!
real(8), external :: erf
real(8), external :: qe_erf
integer :: i, l, j, jj
!
read(iunps,*, end=10, err=10) z, zv, nbeta_, lloc, exfact
@ -85,11 +85,11 @@ subroutine read_oldcp(iunps)
allocate (vnl(mesh_,0:nbeta_))
do l=0,nbeta_
!
! DO NOT USE f90 ARRAY SYNTAX: erf IS NOT AN INTRINSIC FUNCTION!!!
! DO NOT USE f90 ARRAY SYNTAX: qe_erf IS NOT AN INTRINSIC FUNCTION!!!
!
do j=1, mesh_
vnl(j,l)= - (wrc1*erf(sqrt(rc1)*r_(j)) + &
wrc2*erf(sqrt(rc2)*r_(j)) ) * zv/r_(j)
vnl(j,l)= - (wrc1*qe_erf(sqrt(rc1)*r_(j)) + &
wrc2*qe_erf(sqrt(rc2)*r_(j)) ) * zv/r_(j)
end do
!
do i=1,3