Commit e4c171d9 authored by Matthias Redies's avatar Matthias Redies

redo it in 1e

parent 14346a90
...@@ -14,8 +14,8 @@ module m_corespec ...@@ -14,8 +14,8 @@ module m_corespec
complex, parameter :: cone = cmplx(1.0,0.0) complex, parameter :: cone = cmplx(1.0,0.0)
complex, parameter :: cimu = cmplx(0.0,1.0) complex, parameter :: cimu = cmplx(0.0,1.0)
real, parameter :: alpha = 7.29735257*10.0**-3 real, parameter :: alpha = 7.29735257e-3
real, parameter :: mec2 = 0.51099891*10.0**6 real, parameter :: mec2 = 0.51099891e6
real, parameter :: ecoredeep = 0.5 real, parameter :: ecoredeep = 0.5
integer, parameter :: edgel(11) = (/0,1,1,2,2,3,3,4,4,5,5/) integer, parameter :: edgel(11) = (/0,1,1,2,2,3,3,4,4,5,5/)
......
This diff is collapsed.
...@@ -222,7 +222,7 @@ CONTAINS ...@@ -222,7 +222,7 @@ CONTAINS
END IF END IF
IF (ANY(sym%tau(:, isym2) /= 0)) CYCLE IF (ANY(sym%tau(:, isym2) /= 0)) CYCLE
IF (ALL(ABS(MATMUL(rrot(:, :, isym), kpts%bk(:, ikpt)) - kpts%bk(:, ikpt)) < 10.0**-12)) THEN IF (ALL(ABS(MATMUL(rrot(:, :, isym), kpts%bk(:, ikpt)) - kpts%bk(:, ikpt)) < 1e-12)) THEN
isym1 = isym1 + 1 isym1 = isym1 + 1
sym1(isym1, ikpt) = isym sym1(isym1, ikpt) = isym
END IF END IF
...@@ -577,7 +577,7 @@ CONTAINS ...@@ -577,7 +577,7 @@ CONTAINS
q = MATMUL(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp), cell%bmat) q = MATMUL(kpts%bk(:, ikpt) + hybrid%gptm(:, igptp), cell%bmat)
qnorm = SQRT(SUM(q**2)) qnorm = SQRT(SUM(q**2))
iqnrm = pqnrm(igpt, ikpt) iqnrm = pqnrm(igpt, ikpt)
IF (ABS(qnrm(iqnrm) - qnorm) > 10.0**-12) then IF (ABS(qnrm(iqnrm) - qnorm) > 1e-12) then
STOP 'coulombmatrix: qnorm does not equal corresponding & element in qnrm (bug?)' ! We shouldn't stop here! STOP 'coulombmatrix: qnorm does not equal corresponding & element in qnrm (bug?)' ! We shouldn't stop here!
endif endif
...@@ -1721,7 +1721,7 @@ CONTAINS ...@@ -1721,7 +1721,7 @@ CONTAINS
! !
! Convergence parameter ! Convergence parameter
#define CONVPARAM 10.0**-18 #define CONVPARAM 1e-18
! Do some additional shells ( real-space and Fourier-space sum ) ! Do some additional shells ( real-space and Fourier-space sum )
#define ADDSHELL1 40 #define ADDSHELL1 40
#define ADDSHELL2 0 #define ADDSHELL2 0
...@@ -1886,13 +1886,13 @@ CONTAINS ...@@ -1886,13 +1886,13 @@ CONTAINS
DO i = 1, nptsh DO i = 1, nptsh
IF (ALL(conv /= HUGE(i))) EXIT IF (ALL(conv /= HUGE(i))) EXIT
IF (i /= 1) THEN IF (i /= 1) THEN
IF (ABS(radsh(i) - radsh(i - 1)) > 10.0**-10) ishell = ishell + 1 IF (ABS(radsh(i) - radsh(i - 1)) > 1e-10) ishell = ishell + 1
ENDIF ENDIF
ra = MATMUL(cell%amat, ptsh(:, i)) + rc ra = MATMUL(cell%amat, ptsh(:, i)) + rc
a = scale*SQRT(SUM(ra**2)) a = scale*SQRT(SUM(ra**2))
IF (a == 0) THEN IF (a == 0) THEN
CYCLE CYCLE
ELSE IF (ABS(a - a1) > 10.0**-10) THEN ELSE IF (ABS(a - a1) > 1e-10) THEN
a1 = a a1 = a
rexp = EXP(-a) rexp = EXP(-a)
IF (ishell <= conv(0)) g(0) = rexp/a & IF (ishell <= conv(0)) g(0) = rexp/a &
...@@ -1978,13 +1978,13 @@ CONTAINS ...@@ -1978,13 +1978,13 @@ CONTAINS
conv = HUGE(i) conv = HUGE(i)
DO i = 1, nptsh DO i = 1, nptsh
IF (i > 1) THEN IF (i > 1) THEN
IF (ABS(radsh(i) - radsh(i - 1)) > 10.0**-10) ishell = ishell + 1 IF (ABS(radsh(i) - radsh(i - 1)) > 1e-10) ishell = ishell + 1
ENDIF ENDIF
ki = ptsh(:, i) + k - NINT(k) ! -nint(...) transforms to Wigner-Seitz cell ( i.e. -0.5 <= x,y,z < 0.5 ) ki = ptsh(:, i) + k - NINT(k) ! -nint(...) transforms to Wigner-Seitz cell ( i.e. -0.5 <= x,y,z < 0.5 )
ka = MATMUL(ki, cell%bmat) ka = MATMUL(ki, cell%bmat)
a = SQRT(SUM(ka**2))/scale a = SQRT(SUM(ka**2))/scale
aa = (1 + a**2)**(-1) aa = (1 + a**2)**(-1)
IF (ABS(a - a1) > 10.0**-10) THEN IF (ABS(a - a1) > 1e-10) THEN
a1 = a a1 = a
IF (a == 0) THEN IF (a == 0) THEN
g(0) = pref*(-4) g(0) = pref*(-4)
...@@ -2055,7 +2055,7 @@ CONTAINS ...@@ -2055,7 +2055,7 @@ CONTAINS
! Calculate accuracy of Gamma-decomposition ! Calculate accuracy of Gamma-decomposition
IF (ALL(kpts%bk == 0)) GOTO 4 IF (ALL(kpts%bk == 0)) GOTO 4
a = 10.0**30 ! ikpt = index of shortest non-zero k-point a = 1e30 ! ikpt = index of shortest non-zero k-point
DO i = 2, kpts%nkpt DO i = 2, kpts%nkpt
rdum = SUM(MATMUL(kpts%bk(:, i), cell%bmat)**2) rdum = SUM(MATMUL(kpts%bk(:, i), cell%bmat)**2)
IF (rdum < a) THEN IF (rdum < a) THEN
...@@ -2165,7 +2165,7 @@ CONTAINS ...@@ -2165,7 +2165,7 @@ CONTAINS
ptsh = ptsh(:, pnt) ptsh = ptsh(:, pnt)
nshell = 1 nshell = 1
DO i = 2, nptsh DO i = 2, nptsh
IF (radsh(i) - radsh(i - 1) > 10.0**-10) nshell = nshell + 1 IF (radsh(i) - radsh(i - 1) > 1e-10) nshell = nshell + 1
END DO END DO
IF (lwrite) & IF (lwrite) &
...@@ -2200,7 +2200,7 @@ CONTAINS ...@@ -2200,7 +2200,7 @@ CONTAINS
q = MATMUL(kpts%bk(:, ikpt) + gpt(:, igptp), cell%bmat) q = MATMUL(kpts%bk(:, ikpt) + gpt(:, igptp), cell%bmat)
qnorm = SQRT(SUM(q**2)) qnorm = SQRT(SUM(q**2))
DO j = 1, i DO j = 1, i
IF (ABS(qnrm(j) - qnorm) < 10.0**-12) THEN IF (ABS(qnrm(j) - qnorm) < 1e-12) THEN
pqnrm(igpt, ikpt) = j pqnrm(igpt, ikpt) = j
CYCLE igptloop CYCLE igptloop
END IF END IF
...@@ -2289,7 +2289,7 @@ CONTAINS ...@@ -2289,7 +2289,7 @@ CONTAINS
r2 = MIN(ABS(db/b1), ABS(dc/c1)) r2 = MIN(ABS(db/b1), ABS(dc/c1))
! Ensure numerical stability. If both formulas are not sufficiently stable, the program stops. ! Ensure numerical stability. If both formulas are not sufficiently stable, the program stops.
IF (r1 > r2) THEN IF (r1 > r2) THEN
IF (r1 < 10.0**-6 .AND. l_warn) THEN IF (r1 < 1e-6 .AND. l_warn) THEN
WRITE (6, '(A,E12.5,A,E12.5,A)') 'sphbessel_integral: Warning! Formula One possibly unstable. Ratios:', & WRITE (6, '(A,E12.5,A,E12.5,A)') 'sphbessel_integral: Warning! Formula One possibly unstable. Ratios:', &
r1, '(', r2, ')' r1, '(', r2, ')'
WRITE (6, '(A,2F15.10,I4)') ' Current qnorms and atom type:', q1, q2, itype WRITE (6, '(A,2F15.10,I4)') ' Current qnorms and atom type:', q1, q2, itype
...@@ -2297,7 +2297,7 @@ CONTAINS ...@@ -2297,7 +2297,7 @@ CONTAINS
END IF END IF
sphbessel_integral = s**3/dq*da sphbessel_integral = s**3/dq*da
ELSE ELSE
IF (r2 < 10.0**-6 .AND. l_warn) THEN IF (r2 < 1e-6 .AND. l_warn) THEN
WRITE (6, '(A,E13.5,A,E13.5,A)') 'sphbessel_integral: Warning! Formula Two possibly unstable. Ratios:', & WRITE (6, '(A,E13.5,A,E13.5,A)') 'sphbessel_integral: Warning! Formula Two possibly unstable. Ratios:', &
r2, '(', r1, ')' r2, '(', r1, ')'
WRITE (6, '(A,2F15.10,I4)') ' Current qnorms and atom type:', & WRITE (6, '(A,2F15.10,I4)') ' Current qnorms and atom type:', &
......
...@@ -199,7 +199,7 @@ CONTAINS ...@@ -199,7 +199,7 @@ CONTAINS
END DO END DO
IF (l_real) THEN IF (l_real) THEN
IF (ANY(ABS(AIMAG(exchange)) > 10.0**-10)) THEN IF (ANY(ABS(AIMAG(exchange)) > 1e-10)) THEN
IF (mpi%irank == 0) WRITE (6, '(A)') 'exchangeCore: Warning! Unusually large imaginary component.' IF (mpi%irank == 0) WRITE (6, '(A)') 'exchangeCore: Warning! Unusually large imaginary component.'
WRITE (*, *) MAXVAL(ABS(AIMAG(exchange))) WRITE (*, *) MAXVAL(ABS(AIMAG(exchange)))
STOP 'exchangeCore: Unusually large imaginary component.' STOP 'exchangeCore: Unusually large imaginary component.'
...@@ -384,7 +384,7 @@ CONTAINS ...@@ -384,7 +384,7 @@ CONTAINS
END DO END DO
IF (mat_ex%l_real) THEN IF (mat_ex%l_real) THEN
IF (ANY(ABS(AIMAG(exchange)) > 10.0**-10)) THEN IF (ANY(ABS(AIMAG(exchange)) > 1e-10)) THEN
IF (mpi%irank == 0) WRITE (6, '(A)') 'exchangeCore: Warning! Unusually large imaginary component.' IF (mpi%irank == 0) WRITE (6, '(A)') 'exchangeCore: Warning! Unusually large imaginary component.'
WRITE (*, *) MAXVAL(ABS(AIMAG(exchange))) WRITE (*, *) MAXVAL(ABS(AIMAG(exchange)))
STOP 'exchangeCore: Unusually large imaginary component.' STOP 'exchangeCore: Unusually large imaginary component.'
......
...@@ -513,8 +513,8 @@ CONTAINS ...@@ -513,8 +513,8 @@ CONTAINS
REAL :: expo, rrad, k(3), kv1(3), kv2(3), kv3(3), knorm2 REAL :: expo, rrad, k(3), kv1(3), kv2(3), kv3(3), knorm2
COMPLEX :: cdum COMPLEX :: cdum
expo = 5*10.0**-3 expo = 5e-3
rrad = sqrt(-log(5*10.0**-3)/expo) rrad = sqrt(-log(5e-3)/expo)
cdum = sqrt(expo)*rrad cdum = sqrt(expo)*rrad
divergence = cell%omtil/(tpi_const**2)*sqrt(pi_const/expo)*cerf(cdum) divergence = cell%omtil/(tpi_const**2)*sqrt(pi_const/expo)*cerf(cdum)
rrad = rrad**2 rrad = rrad**2
......
...@@ -755,7 +755,7 @@ MODULE m_kp_perturbation ...@@ -755,7 +755,7 @@ MODULE m_kp_perturbation
DO iband1 = bandi1, bandf1 DO iband1 = bandi1, bandf1
DO iband2 = bandi2, bandf2 DO iband2 = bandi2, bandf2
rdum = eig_irr(iband2, nk) - eig_irr(iband1, nk) rdum = eig_irr(iband2, nk) - eig_irr(iband1, nk)
IF (abs(rdum) > 10.0**-6) THEN !10.0**-6 IF (abs(rdum) > 1e-6) THEN !10.0**-6
dcprod(iband2, iband1, :) = dcprod(iband2, iband1, :)/rdum dcprod(iband2, iband1, :) = dcprod(iband2, iband1, :)/rdum
ELSE ELSE
dcprod(iband2, iband1, :) = 0.0 dcprod(iband2, iband1, :) = 0.0
......
...@@ -644,7 +644,7 @@ CONTAINS ...@@ -644,7 +644,7 @@ CONTAINS
rdum = rdum + 2*rdum1**2 rdum = rdum + 2*rdum1**2
END IF END IF
IF (rdum1 > 10.0**-6) THEN IF (rdum1 > 1e-6) THEN
IF (mpi%irank == 0) THEN IF (mpi%irank == 0) THEN
WRITE (6, '(A)') 'mixedbasis: Bad orthonormality of ' & WRITE (6, '(A)') 'mixedbasis: Bad orthonormality of ' &
//lchar(l)//'-product basis. Increase tolerance.' //lchar(l)//'-product basis. Increase tolerance.'
......
...@@ -562,11 +562,11 @@ CONTAINS ...@@ -562,11 +562,11 @@ CONTAINS
sphi = 0 sphi = 0
cphi = 0 cphi = 0
r = SQRT(SUM(rvec**2)) r = SQRT(SUM(rvec**2))
IF (r > 10.0**-16) THEN IF (r > 1e-16) THEN
rvec1 = rvec/r rvec1 = rvec/r
ctheta = rvec1(3) ctheta = rvec1(3)
stheta = SQRT(rvec1(1)**2 + rvec1(2)**2) stheta = SQRT(rvec1(1)**2 + rvec1(2)**2)
IF (stheta > 10.0**-16) THEN IF (stheta > 1e-16) THEN
cphi = rvec1(1)/stheta cphi = rvec1(1)/stheta
sphi = rvec1(2)/stheta sphi = rvec1(2)/stheta
END IF END IF
......
...@@ -302,7 +302,7 @@ CONTAINS ...@@ -302,7 +302,7 @@ CONTAINS
if (l_real) THEN if (l_real) THEN
CALL commonphase(cdum, zhlp(:, i), dimension%nbasfcn) CALL commonphase(cdum, zhlp(:, i), dimension%nbasfcn)
IF (any(abs(aimag(zhlp(:, i)/cdum)) > 10.0**-8)) THEN IF (any(abs(aimag(zhlp(:, i)/cdum)) > 1e-8)) THEN
WRITE (*, *) maxval(abs(aimag(zhlp(:, i)/cdum))) WRITE (*, *) maxval(abs(aimag(zhlp(:, i)/cdum)))
WRITE (*, *) zhlp WRITE (*, *) zhlp
STOP 'waveftrafo1: Residual imaginary part.' STOP 'waveftrafo1: Residual imaginary part.'
...@@ -422,7 +422,7 @@ CONTAINS ...@@ -422,7 +422,7 @@ CONTAINS
! DO i = 1,dim1 ! DO i = 1,dim1
! DO j = 1,dim2 ! DO j = 1,dim2
! rdum = abs(real(mat(i,j)))+abs(aimag(mat(i,j))) ! rdum = abs(real(mat(i,j)))+abs(aimag(mat(i,j)))
! IF(rdum.gt.10.0**-6) THEN ! IF(rdum.gt.1e-6) THEN
! cfac = mat(i,j)/abs(mat(i,j)) ! cfac = mat(i,j)/abs(mat(i,j))
! GOTO 1 ! GOTO 1
! ELSE IF(rdum.gt.rmax) THEN ! ELSE IF(rdum.gt.rmax) THEN
...@@ -431,17 +431,17 @@ CONTAINS ...@@ -431,17 +431,17 @@ CONTAINS
! END IF ! END IF
! END DO ! END DO
! END DO ! END DO
! IF(1-abs(cfac) .gt.10.0**-8) THEN ; mat = 0 ; RETURN ; END IF ! IF(1-abs(cfac) .gt.1e-8) THEN ; mat = 0 ; RETURN ; END IF
! 1 IF(abs(1-cfac**2).gt.10.0**-8) mat = mat/cfac ! 1 IF(abs(1-cfac**2).gt.1e-8) mat = mat/cfac
! !
! IF(any(abs(aimag(mat)).gt.10.0**-8)) THEN ! IF(any(abs(aimag(mat)).gt.1e-8)) THEN
! WRITE(*,*) maxval(aimag(mat)) ! WRITE(*,*) maxval(aimag(mat))
! STOP 'symmetrize: Residual imaginary part. Symmetrization failed.' ! STOP 'symmetrize: Residual imaginary part. Symmetrization failed.'
! Determine common phase factor and divide by it to make the output matrix real. ! Determine common phase factor and divide by it to make the output matrix real.
CALL commonphase(cfac, mat, dim1*dim2) CALL commonphase(cfac, mat, dim1*dim2)
mat = mat/cfac mat = mat/cfac
IF (any(abs(aimag(mat)) > 10.0**-8)) & IF (any(abs(aimag(mat)) > 1e-8)) &
&STOP 'symmetrize: Residual imaginary part. Symmetrization failed.' &STOP 'symmetrize: Residual imaginary part. Symmetrization failed.'
END IF END IF
...@@ -745,7 +745,7 @@ CONTAINS ...@@ -745,7 +745,7 @@ CONTAINS
CALL commonphase(phase(j, i), vecout1(:, j, i), dim) CALL commonphase(phase(j, i), vecout1(:, j, i), dim)
vecout1(:, j, i) = vecout1(:, j, i)/phase(j, i) vecout1(:, j, i) = vecout1(:, j, i)/phase(j, i)
IF (any(abs(aimag(vecout1(:, j, i))) > 10.0**-8)) THEN IF (any(abs(aimag(vecout1(:, j, i))) > 1e-8)) THEN
WRITE (*, *) vecout1(:, j, i) WRITE (*, *) vecout1(:, j, i)
STOP 'bra_trafo2: Residual imaginary part.' STOP 'bra_trafo2: Residual imaginary part.'
END IF END IF
...@@ -918,8 +918,8 @@ CONTAINS ...@@ -918,8 +918,8 @@ CONTAINS
! rtaual = matmul(invrot,taual(:,iatom)) + tau(:,invtab(iisym)) ! rtaual = matmul(invrot,taual(:,iatom)) + tau(:,invtab(iisym))
! iatom1 = 0 ! iatom1 = 0
! DO ieq1 = 1,neq(itype) ! DO ieq1 = 1,neq(itype)
! IF( all(abs(modulo(rtaual-taual(:,iiatom+ieq1)+10.0**-12,1.0)) ! IF( all(abs(modulo(rtaual-taual(:,iiatom+ieq1)+1e-12,1.0))
! & .lt. 10.0**-10) ) THEN ! The 10.0**-12 is a dirty fix. ! & .lt. 1e-10) ) THEN ! The 1e-12 is a dirty fix.
! iatom1 = iiatom + ieq1 ! iatom1 = iiatom + ieq1
! END IF ! END IF
! END DO ! END DO
...@@ -1412,8 +1412,8 @@ CONTAINS ...@@ -1412,8 +1412,8 @@ CONTAINS
rtaual = matmul(rot, atoms%taual(:, iatom)) + sym%tau(:, iisym) rtaual = matmul(rot, atoms%taual(:, iatom)) + sym%tau(:, iisym)
iatom1 = 0 iatom1 = 0
DO ieq1 = 1, atoms%neq(itype) DO ieq1 = 1, atoms%neq(itype)
IF (all(abs(modulo(rtaual - atoms%taual(:, iiatom + ieq1) + 10.0**-12, 1.0))& IF (all(abs(modulo(rtaual - atoms%taual(:, iiatom + ieq1) + 1e-12, 1.0))&
& < 10.0**-10)) THEN ! The 10.0**-12 is a dirty fix. & < 1e-10)) THEN ! The 1e-12 is a dirty fix.
iatom1 = iiatom + ieq1 iatom1 = iiatom + ieq1
END IF END IF
END DO END DO
...@@ -1646,7 +1646,7 @@ CONTAINS ...@@ -1646,7 +1646,7 @@ CONTAINS
rmax = 0 rmax = 0
DO i = 1, n DO i = 1, n
rdum = abs(carr(i)) rdum = abs(carr(i))
IF (rdum > 10.0**-6) THEN; cfac = carr(i)/rdum; EXIT IF (rdum > 1e-6) THEN; cfac = carr(i)/rdum; EXIT
ELSE IF (rdum > rmax) THEN; cfac = carr(i)/rdum; rmax = rdum ELSE IF (rdum > rmax) THEN; cfac = carr(i)/rdum; rmax = rdum
END IF END IF
END DO END DO
......
...@@ -135,7 +135,7 @@ c -------- ...@@ -135,7 +135,7 @@ c --------
k = k + 1 k = k + 1
packmat_d(k) = mat(i,j) packmat_d(k) = mat(i,j)
# ifdef test # ifdef test
if(abs(mat(j,i)-mat(i,j)).gt.10.0**-8) if(abs(mat(j,i)-mat(i,j)).gt.1e-8)
+ STOP 'packmat_d: input matrix not symmetric' + STOP 'packmat_d: input matrix not symmetric'
# endif # endif
enddo enddo
...@@ -157,7 +157,7 @@ c -------- ...@@ -157,7 +157,7 @@ c --------
packmatcoul_d(k) = (mat(i,j)+mat(j,i))/2. packmatcoul_d(k) = (mat(i,j)+mat(j,i))/2.
! if(abs(mat(j,i)-mat(i,j)).gt.10.0**-6) then ! if(abs(mat(j,i)-mat(i,j)).gt.1e-6) then
! write(*,*) 'packmatcoul_d: input matrix not symmetric; deviation .gt. 1E-06' ! write(*,*) 'packmatcoul_d: input matrix not symmetric; deviation .gt. 1E-06'
! endif ! endif
enddo enddo
...@@ -198,7 +198,7 @@ c -------- ...@@ -198,7 +198,7 @@ c --------
k = k + 1 k = k + 1
packmat_z(k) = mat(i,j) packmat_z(k) = mat(i,j)
# ifdef test # ifdef test
if(abs(conjg(mat(j,i))-mat(i,j)).gt.10.0**-8) if(abs(conjg(mat(j,i))-mat(i,j)).gt.1e-8)
+ stop 'packmat_z: input matrix not Hermitian.' + stop 'packmat_z: input matrix not Hermitian.'
# endif # endif
enddo enddo
...@@ -219,7 +219,7 @@ c -------- ...@@ -219,7 +219,7 @@ c --------
k = k + 1 k = k + 1
packmatcoul_z(k) = (mat(i,j) + conjg(mat(j,i)))/2. packmatcoul_z(k) = (mat(i,j) + conjg(mat(j,i)))/2.
if(abs(conjg(mat(j,i))-mat(i,j)).gt.10.0**-4) then if(abs(conjg(mat(j,i))-mat(i,j)).gt.1e-4) then
stop 'packmatcoul_z: input matrix not Hermitian; & stop 'packmatcoul_z: input matrix not Hermitian; &
& deviation .gt. 1E-04.' & deviation .gt. 1E-04.'
endif endif
......
...@@ -44,7 +44,7 @@ SUBROUTINE gen_map(atoms,sym,oneD,hybrid) ...@@ -44,7 +44,7 @@ SUBROUTINE gen_map(atoms,sym,oneD,hybrid)
ratom = 0 ratom = 0
DO ieq1 = 1,atoms%neq(itype) DO ieq1 = 1,atoms%neq(itype)
IF( all(abs(modulo(rtaual-atoms%taual(:,iatom0 + ieq1)+10.0**-12,1.0)).lt. 10.0**-10) ) THEN IF( all(abs(modulo(rtaual-atoms%taual(:,iatom0 + ieq1)+1e-12,1.0)).lt. 1e-10) ) THEN
ratom = iatom0 + ieq1 ratom = iatom0 + ieq1
hybrid%map ( iatom,isym) = ratom hybrid%map ( iatom,isym) = ratom
hybrid%tvec(:,iatom,isym) = nint ( rtaual-atoms%taual(:,ratom) ) hybrid%tvec(:,iatom,isym) = nint ( rtaual-atoms%taual(:,ratom) )
......
...@@ -137,7 +137,7 @@ ...@@ -137,7 +137,7 @@
DO k = 2,nsym DO k = 2,nsym
rarr = matmul(rrot(:,:,k),bk(:,i)) * kpts%nkpt3 rarr = matmul(rrot(:,:,k),bk(:,i)) * kpts%nkpt3
iarr2 = nint(rarr) iarr2 = nint(rarr)
IF(any(abs(iarr2-rarr).gt.10.0**-10)) THEN IF(any(abs(iarr2-rarr).gt.1e-10)) THEN
WRITE(6,'(A,I3,A)') 'kptgen: Symmetry operation',k, WRITE(6,'(A,I3,A)') 'kptgen: Symmetry operation',k,
& ' incompatible with k-point set.' & ' incompatible with k-point set.'
ldum = .true. ldum = .true.
...@@ -280,7 +280,7 @@ c function modulo1 maps kpoint into first BZ ...@@ -280,7 +280,7 @@ c function modulo1 maps kpoint into first BZ
nkpt3 = (/a,b,c/) nkpt3 = (/a,b,c/)
modulo1 = kpoint*nkpt3 modulo1 = kpoint*nkpt3
help = nint(modulo1) help = nint(modulo1)
IF(any(abs(help-modulo1).gt.10.0**-8)) THEN IF(any(abs(help-modulo1).gt.1e-8)) THEN
modulo1 = kpoint*nkpt3 modulo1 = kpoint*nkpt3
WRITE(*,*) modulo1 WRITE(*,*) modulo1
help = nint(modulo1) help = nint(modulo1)
......
...@@ -138,9 +138,9 @@ C ...@@ -138,9 +138,9 @@ C
do 2200 i=1,3 do 2200 i=1,3
dist=dist+(kvc(i,j)-kvc(i,l))**2 dist=dist+(kvc(i,j)-kvc(i,l))**2
2200 continue 2200 continue
if ( dist.lt.(dm*(1.0-10.0**-6)) ) dm=dist if ( dist.lt.(dm*(1.0-1e-6)) ) dm=dist
2100 continue 2100 continue
if ( d.lt.(dm*(1.0+10.0**-6)) ) goto 2000 if ( d.lt.(dm*(1.0+1e-6)) ) goto 2000
d=dm d=dm
lmin=l lmin=l
2000 continue 2000 continue
...@@ -164,9 +164,9 @@ C ...@@ -164,9 +164,9 @@ C
do 3300 i=1,3 do 3300 i=1,3
dist=dist+(kvc(i,j)-knew(i))**2 dist=dist+(kvc(i,j)-knew(i))**2
3300 continue 3300 continue
if ( dist.lt.(dm*(1.0-10.0**-6)) ) dm=dist if ( dist.lt.(dm*(1.0-1e-6)) ) dm=dist
3200 continue 3200 continue
if ( dm.lt.(dmax*(1.0+10.0**-6)) ) goto 3000 if ( dm.lt.(dmax*(1.0+1e-6)) ) goto 3000
dmax=dm dmax=dm
do 3400 i=1,3 do 3400 i=1,3
kc(i)=knew(i) kc(i)=knew(i)
......
...@@ -161,7 +161,7 @@ c - local - ...@@ -161,7 +161,7 @@ c - local -
pure_intgrf%ierror = NO_ERROR pure_intgrf%ierror = NO_ERROR
! integral from 0 to r1 approximated by leading term in power series expansion ! integral from 0 to r1 approximated by leading term in power series expansion
IF (f(1)*f(2).gt.10.0**-10.and.h.gt.0) THEN IF (f(1)*f(2).gt.1e-10.and.h.gt.0) THEN
IF(f(2).eq.f(1)) THEN IF(f(2).eq.f(1)) THEN
pure_intgrf%value = r1*f(1) pure_intgrf%value = r1*f(1)
ELSE ELSE
...@@ -349,7 +349,7 @@ c - local arrays - ...@@ -349,7 +349,7 @@ c - local arrays -
END IF END IF
! integral from 0 to r1 approximated by leading term in power series expansion (only if h>0) ! integral from 0 to r1 approximated by leading term in power series expansion (only if h>0)
IF(h.gt.0.and.f(1)*f(2).gt.10.0**-10) THEN IF(h.gt.0.and.f(1)*f(2).gt.1e-10) THEN
IF(f(2).eq.f(1)) THEN IF(f(2).eq.f(1)) THEN
intgr = r1*f(1) intgr = r1*f(1)
ELSE ELSE
...@@ -416,7 +416,7 @@ c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ...@@ -416,7 +416,7 @@ c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
modulo1 = kpoint*nkpt3 modulo1 = kpoint*nkpt3
help = nint(modulo1) help = nint(modulo1)
IF(any(abs(help-modulo1).gt.10.0**-10)) THEN IF(any(abs(help-modulo1).gt.1e-10)) THEN
WRITE(6,'(A,F5.3,2('','',F5.3),A)')'modulo1: argument (',kpoint, WRITE(6,'(A,F5.3,2('','',F5.3),A)')'modulo1: argument (',kpoint,
+ ') is not an element of the k-point set.' + ') is not an element of the k-point set.'
CALL juDFT_error( CALL juDFT_error(
...@@ -645,7 +645,7 @@ c "Computation of Special Functions".) ...@@ -645,7 +645,7 @@ c "Computation of Special Functions".)
lsta = lsta2(x,l,15) ! lsta = lsta2(x,l,15) !
END IF ! END IF !
f0 = 0.0 ! f0 = 0.0 !
f1 = 10.0**-100 ! f1 = 1e-100 !
DO ll = lsta,0,-1 ! backward recurrence DO ll = lsta,0,-1 ! backward recurrence
f = f1 / x * (2*ll+3) - f0 ; IF(ll.le.lmax) sphbes(ll) = f ! with arbitrary start values f = f1 / x * (2*ll+3) - f0 ; IF(ll.le.lmax) sphbes(ll) = f ! with arbitrary start values
f0 = f1 ! f0 = f1 !
...@@ -779,11 +779,11 @@ c for l = 0,...,ll in Y(1,...,(ll+1)**2). ...@@ -779,11 +779,11 @@ c for l = 0,...,ll in Y(1,...,(ll+1)**2).
sphi = 0 sphi = 0
cphi = 0 cphi = 0
r = sqrt(sum(rvec**2)) r = sqrt(sum(rvec**2))
IF(r.gt.10.0**-16) THEN IF(r.gt.1e-16) THEN
rvec1 = rvec / r rvec1 = rvec / r
ctheta = rvec1(3) ctheta = rvec1(3)
stheta = sqrt(rvec1(1)**2+rvec1(2)**2) stheta = sqrt(rvec1(1)**2+rvec1(2)**2)
IF(stheta.gt.10.0**-16) THEN IF(stheta.gt.1e-16) THEN
cphi = rvec1(1) / stheta cphi = rvec1(1) / stheta
sphi = rvec1(2) / stheta sphi = rvec1(2) / stheta
END IF END IF
...@@ -848,7 +848,7 @@ c Returns the complex error function. ...@@ -848,7 +848,7 @@ c Returns the complex error function.
i = i + 1 i = i + 1
c = -c * z2 / i c = -c * z2 / i
cerf = cerf + c/(2*i+1) cerf = cerf + c/(2*i+1)
IF(abs(c/(2*i+1)).lt.10.0**-20) EXIT IF(abs(c/(2*i+1)).lt.1e-20) EXIT
END DO END DO
cerf = cerf * 2/sqrt(pi) cerf = cerf * 2/sqrt(pi)
ELSE ! continued fraction using Lentz's method ELSE ! continued fraction using Lentz's method
...@@ -862,13 +862,13 @@ c Returns the complex error function. ...@@ -862,13 +862,13 @@ c Returns the complex error function.
d = ( 2*z1 + i*d )**(-1) d = ( 2*z1 + i*d )**(-1)
delta = c*d delta = c*d
cerf = cerf * delta cerf = cerf * delta
IF(abs(1-delta).lt.10.0**-15) EXIT IF(abs(1-delta).lt.1e-15) EXIT
i = i + 1 i = i + 1
c = z1 + i/c c = z1 + i/c
d = ( z1 + i*d )**(-1) d = ( z1 + i*d )**(-1)
delta = c*d delta = c*d
cerf = cerf * delta cerf = cerf * delta
IF(abs(1-delta).lt.10.0**-15) EXIT IF(abs(1-delta).lt.1e-15) EXIT
IF(i.eq.10000) IF(i.eq.10000)
& STOP 'cerf: Lentz method not converged after 10000 steps.' & STOP 'cerf: Lentz method not converged after 10000 steps.'