From 676a7274c9de92178e07afe1b53f8469108a9f18 Mon Sep 17 00:00:00 2001 From: Paolo Giannozzi Date: Sun, 26 Apr 2020 13:59:11 +0000 Subject: [PATCH] Avoid divisions by zero In a material with a gap treated with tetrahedra, the density of state at Ef is exactly zero and the Ef shift gives NaN if computed with no check. Also: there are several other divisions by zero in the tetrahedron method for phonons, harmless in practice but hindering usage of compiler debug options. --- LR_Modules/dfpt_tetra_mod.f90 | 60 ++++++++++++++++++++++++++++------- PHonon/PH/ef_shift.f90 | 6 +++- PHonon/PH/phq_readin.f90 | 2 +- 3 files changed, 54 insertions(+), 14 deletions(-) diff --git a/LR_Modules/dfpt_tetra_mod.f90 b/LR_Modules/dfpt_tetra_mod.f90 index b6f58d30e..30f031002 100644 --- a/LR_Modules/dfpt_tetra_mod.f90 +++ b/LR_Modules/dfpt_tetra_mod.f90 @@ -256,7 +256,7 @@ SUBROUTINE dfpt_tetra_calc_delta(tfst,tlst,et_col,delta) REAL(dp),INTENT(IN) :: et_col(nbnd,nkstot) REAL(dp),INTENT(OUT) :: delta(nbnd,nkstot) ! - INTEGER :: nspin_lsda, ns, nt, nk, ii, ibnd, jbnd, kbnd, itetra(4), ik + INTEGER :: nspin_lsda, ns, nt, nk, i, ii, ibnd, jbnd, kbnd, itetra(4), ik REAL(dp) :: ei(4,nbnd), e(4), wdos1(4), a(4,4), V, wdos0(4,nbnd), wdos2 ! delta(1:nbnd,1:nkstot) = 0.0_dp @@ -299,7 +299,13 @@ SUBROUTINE dfpt_tetra_calc_delta(tfst,tlst,et_col,delta) call hpsort (4, e, itetra) ! DO ii = 1, 4 - a(ii,1:4) = ( 0.0_dp - e(1:4) ) / (e(ii) - e(1:4)) + DO i = 1, 4 + IF ( i == ii ) THEN + a(ii,i) = 0.0_dp + ELSE + a(ii,i) = ( 0.0_dp - e(i) ) / (e(ii) - e(i)) + END IF + END DO END DO ! IF(e(1) <= 0.0_dp .AND. 0.0_dp < e(2)) THEN @@ -396,7 +402,7 @@ SUBROUTINE dfpt_tetra_calc_beta1(iq,tfst,tlst,et_col,beta) REAL(dp),INTENT(IN) :: et_col(nbnd, nkstot) REAL(dp),INTENT(INOUT) :: beta(nbnd,nbnd,nkstot) ! - INTEGER :: nspin_lsda, ns, nt, nk, ii, jj, ibnd, itetra(4), ik + INTEGER :: nspin_lsda, ns, nt, nk, i, ii, jj, ibnd, itetra(4), ik REAL(dp) :: thr = 1e-8_dp, V, tsmall(4,4), ei0(4,nbnd), ej0(4,nbnd), e(4), & & ei1(4), ej1(4,nbnd), a(4,4), w0(nbnd,nbnd,4), w1(nbnd,4) ! @@ -439,7 +445,13 @@ SUBROUTINE dfpt_tetra_calc_beta1(iq,tfst,tlst,et_col,beta) call hpsort (4, e, itetra) ! DO ii = 1, 4 - a(ii,1:4) = ( 0.0_dp - e(1:4)) / (e(ii) - e(1:4)) + DO i = 1, 4 + IF ( i == ii ) THEN + a(ii,i) = 0.0_dp + ELSE + a(ii,i) = ( 0.0_dp - e(i) ) / (e(ii) - e(i)) + END IF + END DO END DO ! IF(e(1) <= 0.0_dp .AND. 0.0_dp < e(2)) THEN @@ -707,7 +719,7 @@ SUBROUTINE dfpt_tetra_calc_beta2(iq,tfst,tlst,et_col,beta) REAL(dp),INTENT(IN) :: et_col(nbnd,nkstot) REAL(dp),INTENT(INOUT) :: beta(nbnd,nbnd,nkstot) ! - INTEGER :: nspin_lsda, ns, nt, nk, ii, jj, ibnd, itetra(4), ik + INTEGER :: nspin_lsda, ns, nt, nk, i, ii, jj, ibnd, itetra(4), ik REAL(dp) :: thr = 1e-8_dp, V, tsmall(4,4), ei0(4,nbnd), ej0(4,nbnd), e(4), & & ei1(4), ej1(4,nbnd), a(4,4), w0(nbnd,nbnd,4), w1(nbnd,4) ! @@ -752,7 +764,13 @@ SUBROUTINE dfpt_tetra_calc_beta2(iq,tfst,tlst,et_col,beta) call hpsort (4, e, itetra) ! DO ii = 1, 4 - a(ii,1:4) = ( 0.0_dp - e(1:4)) / (e(ii) - e(1:4)) + DO i = 1, 4 + IF ( i == ii ) THEN + a(ii,i) = 0.0_dp + ELSE + a(ii,i) = ( 0.0_dp - e(i) ) / (e(ii) - e(i)) + END IF + END DO END DO ! IF(e(1) <= 0.0_dp .AND. 0.0_dp < e(2)) THEN @@ -959,7 +977,7 @@ SUBROUTINE dfpt_tetra_calc_beta3(iq,tfst,tlst,et_col,beta) REAL(dp),INTENT(IN) :: et_col(nbnd,nkstot) REAL(dp),INTENT(INOUT) :: beta(nbnd,nbnd,nkstot) ! - INTEGER :: nspin_lsda, ns, nt, nk, ii, jj, ibnd, itetra(4), ik + INTEGER :: nspin_lsda, ns, nt, nk, i, ii, jj, ibnd, itetra(4), ik REAL(dp) :: thr = 1e-8_dp, V, tsmall(4,4), ei0(4,nbnd), ej0(4,nbnd), e(4), & & ei1(4), ej1(4,nbnd), a(4,4), w0(nbnd,nbnd,4), w1(nbnd,4) ! @@ -1002,7 +1020,13 @@ SUBROUTINE dfpt_tetra_calc_beta3(iq,tfst,tlst,et_col,beta) call hpsort (4, e, itetra) ! DO ii = 1, 4 - a(ii,1:4) = ( 0.0_dp - e(1:4)) / (e(ii) - e(1:4)) + DO i = 1, 4 + IF ( i == ii ) THEN + a(ii,i) = 0.0_dp + ELSE + a(ii,i) = ( 0.0_dp - e(i) ) / (e(ii) - e(i)) + END IF + END DO END DO ! IF(e(1) <= 0.0_dp .AND. 0.0_dp < e(2)) THEN @@ -1203,7 +1227,7 @@ SUBROUTINE dfpt_tetra2_theta(ei0,ej0,w0) REAL(dp),INTENT(IN) :: ei0(4), ej0(4,nbnd) REAL(dp),INTENT(OUT) :: w0(nbnd,4) ! - INTEGER :: ii, ibnd, itetra(4) + INTEGER :: i, ii, ibnd, itetra(4) REAL(dp) :: C(3), e(4), a(4,4), thr = 1.0e-8_dp ! w0(1:nbnd,1:4) = 0.0_dp @@ -1215,7 +1239,13 @@ SUBROUTINE dfpt_tetra2_theta(ei0,ej0,w0) call hpsort (4, e, itetra) ! DO ii = 1, 4 - a(ii,1:4) = (0_dp - e(1:4)) / (e(ii) - e(1:4)) + DO i = 1, 4 + IF ( i == ii ) THEN + a(ii,i) = 0.0_dp + ELSE + a(ii,i) = ( 0.0_dp - e(i) ) / (e(ii) - e(i)) + END IF + END DO END DO ! IF(ABS(e(1)) < thr .AND. ABS(e(4)) < thr) THEN @@ -1280,7 +1310,7 @@ SUBROUTINE dfpt_tetra2_lindhard(ei0,ej0,w0) REAL(dp),INTENT(IN) :: ei0(4), ej0(4,nbnd) REAL(dp),INTENT(OUT) :: w0(nbnd,4) ! - INTEGER :: ii, ibnd, itetra(4) + INTEGER :: i, ii, ibnd, itetra(4) REAL(dp) :: V, ei1(4), ej1(4), w1(4), thr = 1e-8_dp REAL(dp) :: e(4), a(4,4), tsmall(4,4) ! @@ -1293,7 +1323,13 @@ SUBROUTINE dfpt_tetra2_lindhard(ei0,ej0,w0) call hpsort (4, e, itetra) ! DO ii = 1, 4 - a(ii,1:4) = ( 0.0_dp - e(1:4) ) / (e(ii) - e(1:4)) + DO i = 1, 4 + IF ( i == ii ) THEN + a(ii,i) = 0.0_dp + ELSE + a(ii,i) = ( 0.0_dp - e(i) ) / (e(ii) - e(i)) + END IF + END DO END DO ! IF(0_dp <= e(1)) THEN diff --git a/PHonon/PH/ef_shift.f90 b/PHonon/PH/ef_shift.f90 index bd8eca6bc..697689fc2 100644 --- a/PHonon/PH/ef_shift.f90 +++ b/PHonon/PH/ef_shift.f90 @@ -97,7 +97,11 @@ subroutine ef_shift (drhoscf, ldos, ldoss, dos_ef, irr, npe, flag) CALL invfft ('Rho', drhoscf(:,is,ipert), dfftp) enddo call mp_sum ( delta_n, intra_bgrp_comm ) - def (ipert) = - delta_n / dos_ef + IF ( ABS(dos_ef) > 1.d-18 ) THEN + def (ipert) = - delta_n / dos_ef + ELSE + def (ipert) = 0.0_dp + ENDIF enddo ! ! symmetrizes the Fermi energy shift diff --git a/PHonon/PH/phq_readin.f90 b/PHonon/PH/phq_readin.f90 index c2c636f21..7671fcca0 100644 --- a/PHonon/PH/phq_readin.f90 +++ b/PHonon/PH/phq_readin.f90 @@ -225,7 +225,7 @@ SUBROUTINE phq_readin() ! Rewind the input if the title is actually the beginning of inputph namelist ! IF( imatches("&inputph", title) ) THEN - WRITE(*, '(6x,a)') "Title line not specified: using 'default'." + WRITE(stdout,'(6x,a)') "Title line not specified: using 'default'." title='default' IF (meta_ionode) REWIND(5, iostat=ios) CALL mp_bcast(ios, meta_ionode_id, world_comm )