...
 
Commits (18)
Problems fixed in development branch :
* Unphysical fractional translations (tau/n with n/=2,3,4,6) were not
explicitly discarded, thus leading in unfortunate cases to strange
values for FFT factors and grids. Also: if "nosym" is true, inversion
symmetry flag (invsym) and info on FFT factors (fft_fact) must also
be reset (problem spotted by Thomas Brumme, Leipzig)
* PPACF wasn't working any longer in v.6.4 and 6.4.1 for nspin=2 and
for hybrid functionals (fixed by Yang Jiao, Chalmers)
* option "write_unkg" of pw2wannier90.f90 wasn't working as expected
* The input parameters (for restarting DFPT+U calculations) read_dns_bare
and d2ns_type were missing in the PH input namelist, and moreover they
were not broadcasted.
* Input parameters (for restarting DFPT+U calculations) read_dns_bare
and d2ns_type were missing in the PH input namelist, and moreover
they were not broadcasted.
Incompatible changes in development branch :
* fractional translations "ftau" in FFT grid units no longer existing as
......
......@@ -54,32 +54,32 @@
!
!
CALL start_clock('a2F')
IF (mpime .eq. ionode_id) THEN
IF (mpime == ionode_id) THEN
!
DO isig = 1, nsmear
!
IF ( isig .lt. 10 ) THEN
IF ( isig < 10 ) THEN
WRITE(fila2f_suffix,'(a,a6,i1)') TRIM(prefix),'.a2f.0', isig
ELSE
WRITE(fila2f_suffix,'(a,a5,i2)') TRIM(prefix),'.a2f.', isig
ENDIF
OPEN (unit = iua2ffil, file = fila2f_suffix, form = 'formatted')
!
IF ( isig .lt. 10 ) THEN
IF ( isig < 10 ) THEN
WRITE(fila2ftr,'(a,a9,i1)') TRIM(prefix),'.a2f_tr.0', isig
ELSE
WRITE(fila2ftr,'(a,a8,i2)') TRIM(prefix),'.a2f_tr.', isig
ENDIF
OPEN (unit = iua2ftrfil, file = fila2ftr, form = 'formatted')
!
IF ( isig .lt. 10 ) THEN
IF ( isig < 10 ) THEN
WRITE(filres,'(a,a6,i1)') TRIM(prefix),'.res.0', isig
ELSE
WRITE(filres,'(a,a5,i2)') TRIM(prefix),'.res.', isig
ENDIF
OPEN (unit = iures, file = filres, form = 'formatted')
!
IF ( isig .lt. 10 ) THEN
IF ( isig < 10 ) THEN
WRITE(fildos,'(a,a8,i1)') TRIM(prefix),'.phdos.0', isig
ELSE
WRITE(fildos,'(a,a7,i2)') TRIM(prefix),'.phdos.', isig
......@@ -90,16 +90,16 @@
WRITE(stdout,'(5x,"Eliashberg Spectral Function in the Migdal Approximation")')
WRITE(stdout,'(5x,a/)') REPEAT('=',67)
!
IF ( .not. ALLOCATED(a2F) ) ALLOCATE( a2F(nqstep, nqsmear) )
IF ( .not. ALLOCATED(a2F_tr) ) ALLOCATE( a2F_tr(nqstep, nqsmear) )
IF ( .not. ALLOCATED(dosph) ) ALLOCATE( dosph(nqstep, nqsmear) )
IF ( .not. ALLOCATED(l_a2F) ) ALLOCATE( l_a2F(nqsmear) )
IF ( .not. ALLOCATED(l_a2F_tr) ) ALLOCATE( l_a2F_tr(nqsmear) )
IF ( .not. ALLOCATED(logavg) ) ALLOCATE( logavg(nqsmear) )
IF ( .NOT. ALLOCATED(a2F) ) ALLOCATE ( a2F(nqstep, nqsmear) )
IF ( .NOT. ALLOCATED(a2F_tr) ) ALLOCATE ( a2F_tr(nqstep, nqsmear) )
IF ( .NOT. ALLOCATED(dosph) ) ALLOCATE ( dosph(nqstep, nqsmear) )
IF ( .NOT. ALLOCATED(l_a2F) ) ALLOCATE ( l_a2F(nqsmear) )
IF ( .NOT. ALLOCATED(l_a2F_tr) ) ALLOCATE ( l_a2F_tr(nqsmear) )
IF ( .NOT. ALLOCATED(logavg) ) ALLOCATE ( logavg(nqsmear) )
!
! The resitivity is computed for temperature between 0K-1000K by step of 10
! This is hardcoded and needs to be changed here if one wants to modify it
IF ( .not. ALLOCATED(rho) ) ALLOCATE( rho(100, nqsmear) )
IF ( .NOT. ALLOCATED(rho) ) ALLOCATE ( rho(100, nqsmear) )
!
!om_max = ( MAXVAL( wf(:,:) ) - MINVAL( wf(:,:) ) ) + 5.d0/ryd2mev
!om_max = MAXVAL( wf(:,:) ) + 1.d0 / ryd2mev
......@@ -131,10 +131,10 @@
!
w0 = wf(imode,iq)
!
IF ( w0 .gt. eps_acustic ) THEN
IF ( w0 > eps_acustic ) THEN
!
l = lambda_all(imode,iq,isig)
IF (lambda_all(imode,iq,isig) .lt. 0.d0) l = 0.d0 ! sanity check
IF (lambda_all(imode,iq,isig) < 0.d0) l = 0.d0 ! sanity check
!
a2F_tmp = wqf(iq) * w0 * l / two
!
......@@ -143,7 +143,7 @@
dosph(iw,ismear) = dosph(iw,ismear) + wqf(iq) * weight
!
l_tr = lambda_v_all(imode,iq,isig)
IF (lambda_v_all(imode,iq,isig) .lt. 0.d0) l_tr = 0.d0 !sanity check
IF (lambda_v_all(imode,iq,isig) < 0.d0) l_tr = 0.d0 !sanity check
!
a2F_tr_tmp = wqf(iq) * w0 * l_tr / two
!
......@@ -157,13 +157,13 @@
!
! output a2F
!
IF (ismear .eq. nqsmear) WRITE (iua2ffil, '(f12.7, 15f12.7)') iomega*ryd2mev, a2F(iw,:)
IF (ismear .eq. nqsmear) WRITE (iua2ftrfil, '(f12.7, 15f12.7)') iomega*ryd2mev, a2F_tr(iw,:)
IF (ismear .eq. nqsmear) WRITE (iudosfil, '(f12.7, 15f12.7)') iomega*ryd2mev, dosph(iw,:)/ryd2mev
IF (ismear == nqsmear) WRITE (iua2ffil, '(f12.7, 15f12.7)') iomega*ryd2mev, a2F(iw,:)
IF (ismear == nqsmear) WRITE (iua2ftrfil, '(f12.7, 15f12.7)') iomega*ryd2mev, a2F_tr(iw,:)
IF (ismear == nqsmear) WRITE (iudosfil, '(f12.7, 15f12.7)') iomega*ryd2mev, dosph(iw,:)/ryd2mev
!
! do the integral 2 int (a2F(w)/w dw)
!
!IF (iomega .gt. eps_acustic) &
!IF (iomega > eps_acustic) &
l_a2F(ismear) = l_a2F(ismear) + two * a2F(iw,ismear) / iomega * dw
l_a2F_tr(ismear) = l_a2F_tr(ismear) + two * a2F_tr(iw,ismear) / iomega * dw
logavg(ismear) = logavg(ismear) + two * a2F(iw,ismear) * log(iomega) / iomega * dw
......@@ -176,9 +176,9 @@
!
DO iq = 1, nqtotf ! loop over q-points
DO imode = 1, nmodes ! loop over modes
IF (lambda_all(imode,iq,isig) .gt. 0.d0 .and. wf(imode,iq) .gt. eps_acustic ) &
IF (lambda_all(imode,iq,isig) > 0.d0 .and. wf(imode,iq) > eps_acustic ) &
lambda_tot = lambda_tot + wqf(iq) * lambda_all(imode,iq,isig)
IF (lambda_v_all(imode,iq,isig) .gt. 0.d0 .and. wf(imode,iq) .gt. eps_acustic ) &
IF (lambda_v_all(imode,iq,isig) > 0.d0 .and. wf(imode,iq) > eps_acustic ) &
lambda_tr_tot = lambda_tr_tot + wqf(iq) * lambda_v_all(imode,iq,isig)
ENDDO
ENDDO
......@@ -236,7 +236,7 @@
! Conductivity 1 a.u. = 2.2999241E6 S/m
! Now to go from Ohm*m to micro Ohm cm we need to multiply by 1E8
rho(itemp,ismear) = rho(itemp,ismear) * 1E8 / 2.2999241E6
IF (ismear .eq. nqsmear) WRITE (iures, '(i8, 15f12.7)') itemp * 10, rho(itemp,:)
IF (ismear == nqsmear) WRITE (iures, '(i8, 15f12.7)') itemp * 10, rho(itemp,:)
ENDDO
ENDDO
CLOSE(iures)
......@@ -261,13 +261,13 @@
!
CLOSE(iudosfil)
!
IF ( ALLOCATED(l_a2F) ) DEALLOCATE(l_a2F)
IF ( ALLOCATED(l_a2F_tr) ) DEALLOCATE(l_a2F_tr)
IF ( ALLOCATED(a2F) ) DEALLOCATE(a2F)
IF ( ALLOCATED(a2F_tr) ) DEALLOCATE(a2F_tr)
IF ( ALLOCATED(rho) ) DEALLOCATE(rho)
IF ( ALLOCATED(dosph) ) DEALLOCATE(dosph)
IF ( ALLOCATED(logavg) ) DEALLOCATE(logavg)
IF ( ALLOCATED(l_a2F) ) DEALLOCATE (l_a2F)
IF ( ALLOCATED(l_a2F_tr) ) DEALLOCATE (l_a2F_tr)
IF ( ALLOCATED(a2F) ) DEALLOCATE (a2F)
IF ( ALLOCATED(a2F_tr) ) DEALLOCATE (a2F_tr)
IF ( ALLOCATED(rho) ) DEALLOCATE (rho)
IF ( ALLOCATED(dosph) ) DEALLOCATE (dosph)
IF ( ALLOCATED(logavg) ) DEALLOCATE (logavg)
!
ENDDO ! isig
!
......
......@@ -73,7 +73,7 @@
COMPLEX(kind=DP) :: sum_nc(npol)
!! auxiliary sum variable non-collinear case
!
IF (.not.okvan) RETURN
IF ( .NOT. okvan) RETURN
!
CALL start_clock('adddvscf2')
!
......@@ -83,7 +83,7 @@
DO nt = 1, ntyp
IF ( upf(nt)%tvanp ) THEN
DO na = 1, nat
IF (ityp(na) .eq. nt) THEN
IF (ityp(na) == nt) THEN
!
! we multiply the integral for the becp term and the beta_n
!
......@@ -124,7 +124,7 @@
ENDDO
ELSE
DO na = 1, nat
IF (ityp(na) .eq. nt) ijkb0 = ijkb0 + nh(nt)
IF (ityp(na) == nt) ijkb0 = ijkb0 + nh(nt)
ENDDO
ENDIF
ENDDO
......
......@@ -240,7 +240,7 @@ SUBROUTINE bcast_epw_input1
! integers
!
CALL mp_bcast (nat_todo, meta_ionode_id, world_comm)
IF (nat_todo.gt.0) THEN
IF (nat_todo > 0) THEN
CALL mp_bcast (atomo, meta_ionode_id, world_comm)
ENDIF
#endif
......
......@@ -120,7 +120,7 @@
!
et_tmp = zero
et_opt = zero
IF (nexband_tmp .gt. 0) THEN
IF (nexband_tmp > 0) THEN
DO ik = 1, nks
ibnd = 0
DO i = 1, nbnd
......@@ -215,9 +215,9 @@
! [mind when comparing with wannier code (eV and angstrom units) with
! write_hr=.true.]
!
IF (mpime .eq. ionode_id) THEN
IF (mpime == ionode_id) THEN
!
OPEN(unit=iudecayH,file='decay.H')
OPEN(UNIT=iudecayH,FILE='decay.H')
WRITE(iudecayH, '(/3x,a/)') '#Spatial decay of Hamiltonian in Wannier basis'
DO ir = 1, nrr
!
......@@ -351,7 +351,7 @@
!
dmec_opt = czero
dmec_tmp = czero
IF (nexband_tmp .gt. 0) THEN
IF (nexband_tmp > 0) THEN
DO ik = 1,nks
jbnd = 0
DO j = 1, nbnd
......@@ -414,8 +414,8 @@
DO ik = 1, nks
DO ipol = 1, 3
!
! dmec_utmp(:,:) = matmul( dmec_opt(ipol,:,:,ik), cu(:,:,ik) )
! cps(ipol,:,:,ik) = matmul( conjg(transpose( cu(:,:,ik))), dmec_utmp(:,:) )
! dmec_utmp(:,:) = MATMUL( dmec_opt(ipol,:,:,ik), cu(:,:,ik) )
! cps(ipol,:,:,ik) = MATMUL( conjg(transpose( cu(:,:,ik))), dmec_utmp(:,:) )
!
CALL zgemm ('n', 'n', nbnd, nbndsub, nbnd, cone, dmec_opt(ipol,:,:,ik), &
nbnd, cu(:,:,ik), nbnd, czero, dmec_utmp(:,:), nbnd)
......@@ -459,8 +459,8 @@
! Check spatial decay of Dipole in Wannier basis
! the unit in r-space is angstrom
!
IF (mpime.eq.ionode_id) THEN
OPEN(unit=iudecayP,file='decay.P')
IF (mpime == ionode_id) THEN
OPEN(UNIT=iudecayP,FILE='decay.P')
WRITE(iudecayP, '(/3x,a/)') '#Spatial decay of dipole in Wannier basis'
DO ir = 1, nrr
!
......@@ -599,8 +599,8 @@
! the unit in r-space is angstrom, and I am plotting
! the matrix for the first mode only
!
IF (mpime.eq.ionode_id) THEN
OPEN(unit=iudecaydyn,file='decay.dynmat')
IF (mpime == ionode_id) THEN
OPEN(UNIT=iudecaydyn,FILE='decay.dynmat')
WRITE(iudecaydyn, '(/3x,a/)') '#Spatial decay of Dynamical matrix in Wannier basis'
DO ir = 1, nrr
!
......@@ -730,7 +730,7 @@
REAL(kind=DP) :: zero_vect(3)
!! temporary zero vector
REAL(kind=DP) :: delta
!! \delta_nm = 1 if n .eq. m and 0 if n .neq. m
!! \delta_nm = 1 if n == m and 0 if n /= m
!
COMPLEX(kind=DP) :: Apos(3,nbndsub,nbndsub,nks)
!! A^W_{mn,\alpha}(k)
......@@ -764,8 +764,8 @@
! RM - bvec can be writen on file by making a small change in
! W90/hamiltonian.F90/hamilotonian_write_rmn
!
tempfile=trim(prefix)//'.bvec'
OPEN(iubvec, file=tempfile, action='read', iostat=ios)
tempFILE=trim(prefix)//'.bvec'
OPEN(iubvec, FILE=tempfile, action='read', iostat=ios)
IF (ios /= 0) THEN
!
! if it doesn't exist, then we just set the bvec and wb to zero
......@@ -800,9 +800,9 @@
ALLOCATE ( M_mn(nbnd, nbnd, nnb, nkstot) )
M_mn = czero
!
IF (mpime.eq.ionode_id) THEN
tempfile=trim(prefix)//'.mmn'
OPEN(iummn, file=tempfile, status = 'old', form = 'formatted', iostat=ios)
IF (mpime == ionode_id) THEN
tempFILE=trim(prefix)//'.mmn'
OPEN(iummn, FILE=tempfile, status = 'old', form = 'formatted', iostat=ios)
!
IF (ios /= 0) THEN
! if it doesn't exist, then we just set the mmn to zero
......@@ -841,7 +841,7 @@
m_mat_tmp(:,:,:,:) = czero
zero_vect(:) = zero
!
IF (nexband_tmp .gt. 0) THEN
IF (nexband_tmp > 0) THEN
DO ik = 1, nks
CALL ktokpmq ( xk(:,ik), zero_vect, +1, ipool, nkk, nkk_abs)
!
......@@ -921,8 +921,8 @@
b_tmp(:) = alat / (twopi) * bvec(:,ib,nkk_abs)
CALL ktokpmq ( xk(:,ik), b_tmp(:), +1, ipool, nkb, nkb_abs)
!
! M_mn_utmp(:,:) = matmul( m_mat_opt(:,:,ib,ik), cu_big(:,:,nkb_abs) )
! cvs(:,:,ib,ik) = matmul( conjg(transpose(cu(:,:,ik))), M_mn_utmp(:,:) )
! M_mn_utmp(:,:) = MATMUL( m_mat_opt(:,:,ib,ik), cu_big(:,:,nkb_abs) )
! cvs(:,:,ib,ik) = MATMUL( conjg(transpose(cu(:,:,ik))), M_mn_utmp(:,:) )
!
CALL zgemm ('n', 'n', nbnd, nbndsub, nbnd, cone, m_mat_opt(:,:,ib,ik), &
nbnd, cu_big(:,:,nkb_abs), nbnd, czero, M_mn_utmp(:,:), nbnd)
......@@ -1004,8 +1004,8 @@
! position matrix cvmew and spatial dimensions are in units of bohr
! [mind when comparing with wannier code (angstrom units) with write_rmn=.true.]
!
IF (mpime.eq.ionode_id) then
OPEN(unit=iudecayv,file='decay.v')
IF (mpime == ionode_id) then
OPEN(UNIT=iudecayv,FILE='decay.v')
WRITE(iudecayv, '(/3x,a/)') '#Spatial decay of Velocity matrix element in Wannier basis'
DO ir = 1, nrr
!
......@@ -1171,8 +1171,8 @@
! the unit in r-space is angstrom, and I am plotting
! the matrix for the first mode only
!
IF (mpime.eq.ionode_id) THEN
OPEN(unit=iuwane,file='decay.epwane')
IF (mpime == ionode_id) THEN
OPEN(UNIT=iuwane,FILE='decay.epwane')
WRITE(iuwane, '(a)') '# Spatial decay of e-p matrix elements in Wannier basis'
DO ir = 1, nrr
!
......@@ -1291,9 +1291,9 @@
!
! we plot: R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)|
!
IF (mpime.eq.ionode_id) THEN
IF (ir.eq.1) open(unit=iuwanep,file='decay.epmat_wanep',status='unknown')
IF (ir.eq.1) WRITE(iuwanep, '(a)') '# R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)| '
IF (mpime == ionode_id) THEN
IF (ir == 1) open(UNIT=iuwanep,FILE='decay.epmat_wanep',status='unknown')
IF (ir == 1) WRITE(iuwanep, '(a)') '# R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)| '
DO ire = 1, nrr_k
!
rvec1 = dble(irvec_k(1,ire))*at(:,1) + &
......@@ -1312,7 +1312,7 @@
WRITE(iuwanep, '(5f15.10)') len1 * alat * bohr2ang, &
len2 * alat * bohr2ang, tmp
ENDDO
IF (ir.eq.nrr_g) CLOSE(iuwanep)
IF (ir == nrr_g) CLOSE(iuwanep)
ENDIF
!
ENDDO
......@@ -1393,7 +1393,7 @@
COMPLEX(KIND=DP), ALLOCATABLE :: epmatwp_mem(:,:,:,:)
!! e-p matrix in Wannier basis
!
ALLOCATE (epmatwp_mem( nbnd, nbnd, nrr_k, nmodes))
ALLOCATE (epmatwp_mem(nbnd, nbnd, nrr_k, nmodes))
!
!----------------------------------------------------------
! Fourier transform to go into Wannier basis
......@@ -1429,7 +1429,7 @@
! we plot: R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)|
!
IF (mpime == ionode_id) THEN
IF (ir == 1) OPEN(unit=iuwanep, file='decay.epmat_wanep', status='unknown')
IF (ir == 1) OPEN(UNIT=iuwanep, FILE='decay.epmat_wanep', status='unknown')
IF (ir == 1) WRITE(iuwanep, '(a)') '# R_e, R_p, max_{m,n,nu} |g(m,n,nu;R_e,R_p)| '
DO ire = 1, nrr_k
!
......@@ -1458,7 +1458,7 @@
!
CALL cryst_to_cart (nq, xk, bg, 1)
!
IF ( ALLOCATED (epmatwp_mem) ) DEALLOCATE (epmatwp_mem)
DEALLOCATE (epmatwp_mem)
!
END SUBROUTINE ephbloch2wanp_mem
!
......
......@@ -39,7 +39,7 @@
!
! Here the local variables
!
! max number of iterations used in mixing: n_iter must be .le. maxter
! max number of iterations used in mixing: n_iter must be <= maxter
INTEGER, PARAMETER :: maxter = 8
!
INTEGER :: n, i, j, iwork(maxter), info, iter_used, ipos, inext
......@@ -54,20 +54,20 @@
REAL(DP) wg(maxter), wg0
DATA wg0 / 0.01d0 /, wg / maxter * 1.d0 /
!
IF ( iter .lt. 1 ) CALL errore('mix_broyden','n_iter is smaller than 1',1)
IF ( n_iter .gt. maxter ) CALL errore('mix_broyden','n_iter is too big',1)
IF ( ndim .le. 0 ) CALL errore('mix_broyden','ndim .le. 0',1)
IF ( iter < 1 ) CALL errore('mix_broyden','n_iter is smaller than 1',1)
IF ( n_iter > maxter ) CALL errore('mix_broyden','n_iter is too big',1)
IF ( ndim <= 0 ) CALL errore('mix_broyden','ndim <= 0',1)
!
IF ( iter .eq. 1 ) THEN
IF ( .not. ALLOCATED(df) ) ALLOCATE( df(ndim,n_iter) )
IF ( .not. ALLOCATED(dv) ) ALLOCATE( dv(ndim,n_iter) )
IF ( iter == 1 ) THEN
IF ( .NOT. ALLOCATED(df) ) ALLOCATE ( df(ndim,n_iter) )
IF ( .NOT. ALLOCATED(dv) ) ALLOCATE ( dv(ndim,n_iter) )
ENDIF
IF ( conv .OR. iter .eq. nsiter ) THEN
IF ( ALLOCATED(df) ) DEALLOCATE(df)
IF ( ALLOCATED(dv) ) DEALLOCATE(dv)
IF ( conv .OR. iter == nsiter ) THEN
IF ( ALLOCATED(df) ) DEALLOCATE (df)
IF ( ALLOCATED(dv) ) DEALLOCATE (dv)
RETURN
ENDIF
IF ( .not. ALLOCATED(deltainsave) ) ALLOCATE( deltainsave(ndim) )
IF ( .NOT. ALLOCATED(deltainsave) ) ALLOCATE ( deltainsave(ndim) )
deltainsave(:) = deltain(:)
!
! iter_used = iter-1 IF iter <= n_iter
......@@ -84,7 +84,7 @@
deltaout(n) = deltaout(n) - deltain(n)
ENDDO
!
IF ( iter .gt. 1 ) THEN
IF ( iter > 1 ) THEN
DO n = 1, ndim
df(n,ipos) = deltaout(n) - df(n,ipos)
dv(n,ipos) = deltain(n) - dv(n,ipos)
......@@ -141,7 +141,7 @@
df(:,inext) = deltaout(:)
dv(:,inext) = deltainsave(:)
!
IF ( ALLOCATED(deltainsave) ) DEALLOCATE(deltainsave)
IF ( ALLOCATED(deltainsave) ) DEALLOCATE (deltainsave)
!
RETURN
!
......@@ -179,7 +179,7 @@
!
! Here the local variables
!
! max number of iterations used in mixing: n_iter must be .le. maxter
! max number of iterations used in mixing: n_iter must be <= maxter
INTEGER, PARAMETER :: maxter = 8
!
INTEGER :: n, i, j, iwork(maxter), info, iter_used, ipos, inext
......@@ -194,20 +194,20 @@
REAL(DP) wg(maxter), wg0
DATA wg0 / 0.01d0 /, wg / maxter * 1.d0 /
!
IF ( iter .lt. 1 ) CALL errore('mix_broyden2','n_iter is smaller than 1',1)
IF ( n_iter .gt. maxter ) CALL errore('mix_broyden2','n_iter is too big',1)
IF ( ndim .le. 0 ) CALL errore('mix_broyden2','ndim .le. 0',1)
IF ( iter < 1 ) CALL errore('mix_broyden2','n_iter is smaller than 1',1)
IF ( n_iter > maxter ) CALL errore('mix_broyden2','n_iter is too big',1)
IF ( ndim <= 0 ) CALL errore('mix_broyden2','ndim <= 0',1)
!
IF ( iter .eq. 1 ) THEN
IF ( .not. ALLOCATED(df2) ) ALLOCATE( df2(ndim,n_iter) )
IF ( .not. ALLOCATED(dv2) ) ALLOCATE( dv2(ndim,n_iter) )
IF ( iter == 1 ) THEN
IF ( .NOT. ALLOCATED(df2) ) ALLOCATE ( df2(ndim,n_iter) )
IF ( .NOT. ALLOCATED(dv2) ) ALLOCATE ( dv2(ndim,n_iter) )
ENDIF
IF ( conv .OR. iter .eq. nsiter ) THEN
IF ( ALLOCATED(df2) ) DEALLOCATE(df2)
IF ( ALLOCATED(dv2) ) DEALLOCATE(dv2)
IF ( conv .OR. iter == nsiter ) THEN
IF ( ALLOCATED(df2) ) DEALLOCATE (df2)
IF ( ALLOCATED(dv2) ) DEALLOCATE (dv2)
RETURN
ENDIF
IF ( .not. ALLOCATED(deltainsave) ) ALLOCATE( deltainsave(ndim) )
IF ( .NOT. ALLOCATED(deltainsave) ) ALLOCATE ( deltainsave(ndim) )
deltainsave(:) = deltain(:)
!
! iter_used = iter-1 IF iter <= n_iter
......@@ -224,7 +224,7 @@
deltaout(n) = deltaout(n) - deltain(n)
ENDDO
!
IF ( iter .gt. 1 ) THEN
IF ( iter > 1 ) THEN
DO n = 1, ndim
df2(n,ipos) = deltaout(n) - df2(n,ipos)
dv2(n,ipos) = deltain(n) - dv2(n,ipos)
......@@ -281,7 +281,7 @@
df2(:,inext) = deltaout(:)
dv2(:,inext) = deltainsave(:)
!
IF ( ALLOCATED(deltainsave) ) DEALLOCATE(deltainsave)
IF ( ALLOCATED(deltainsave) ) DEALLOCATE (deltainsave)
!
RETURN
!
......@@ -324,7 +324,7 @@
!
! Here the local variables
!
! max number of iterations used in mixing: n_iter must be .le. maxter
! max number of iterations used in mixing: n_iter must be <= maxter
INTEGER, PARAMETER :: maxter = 8
!
INTEGER :: n, i, j, iwork(maxter), info, iter_used, ipos, inext
......@@ -340,20 +340,20 @@
DATA wg0 / 0.01d0 /, wg / maxter * 1.d0 /
REAL(DP) :: df_(ndim,n_iter), dv_(ndim,n_iter)
!
IF ( iter .lt. 1 ) CALL errore('mix_broyden','n_iter is smaller than 1',1)
IF ( n_iter .gt. maxter ) CALL errore('mix_broyden','n_iter is too big',1)
IF ( ndim .le. 0 ) CALL errore('mix_broyden','ndim .le. 0',1)
IF ( iter < 1 ) CALL errore('mix_broyden','n_iter is smaller than 1',1)
IF ( n_iter > maxter ) CALL errore('mix_broyden','n_iter is too big',1)
IF ( ndim <= 0 ) CALL errore('mix_broyden','ndim <= 0',1)
!
IF ( iter .eq. 1 ) THEN
IF ( .not. ALLOCATED(df) ) ALLOCATE( df(nbndfs,nkfs,ndim,n_iter) )
IF ( .not. ALLOCATED(dv) ) ALLOCATE( dv(nbndfs,nkfs,ndim,n_iter) )
IF ( iter == 1 ) THEN
IF ( .NOT. ALLOCATED(df) ) ALLOCATE ( df(nbndfs,nkfs,ndim,n_iter) )
IF ( .NOT. ALLOCATED(dv) ) ALLOCATE ( dv(nbndfs,nkfs,ndim,n_iter) )
ENDIF
IF ( conv .OR. iter .eq. nsiter ) THEN
IF (ALLOCATED(df)) DEALLOCATE(df)
IF (ALLOCATED(dv)) DEALLOCATE(dv)
IF ( conv .OR. iter == nsiter ) THEN
IF (ALLOCATED(df)) DEALLOCATE (df)
IF (ALLOCATED(dv)) DEALLOCATE (dv)
RETURN
ENDIF
IF ( .not. ALLOCATED(deltainsave) ) ALLOCATE( deltainsave(ndim) )
IF ( .NOT. ALLOCATED(deltainsave) ) ALLOCATE ( deltainsave(ndim) )
deltainsave(:) = deltain(:)
!
! iter_used = iter-1 IF iter <= n_iter
......@@ -370,7 +370,7 @@
deltaout(n) = deltaout(n) - deltain(n)
ENDDO
!
IF ( iter .gt. 1 ) THEN
IF ( iter > 1 ) THEN
DO n = 1, ndim
df(ibnd,ik,n,ipos) = deltaout(n) - df(ibnd,ik,n,ipos)
dv(ibnd,ik,n,ipos) = deltain(n) - dv(ibnd,ik,n,ipos)
......@@ -429,7 +429,7 @@
df(ibnd,ik,:,inext) = deltaout(:)
dv(ibnd,ik,:,inext) = deltainsave(:)
!
IF ( ALLOCATED(deltainsave) ) DEALLOCATE(deltainsave)
IF ( ALLOCATED(deltainsave) ) DEALLOCATE (deltainsave)
!
RETURN
!
......@@ -471,7 +471,7 @@
!
! Here the local variables
!
! max number of iterations used in mixing: n_iter must be .le. maxter
! max number of iterations used in mixing: n_iter must be <= maxter
INTEGER, PARAMETER :: maxter = 8
!
INTEGER :: n, i, j, iwork(maxter), info, iter_used, ipos, inext
......@@ -487,20 +487,20 @@
DATA wg0 / 0.01d0 /, wg / maxter * 1.d0 /
REAL(DP) :: df_(ndim,n_iter), dv_(ndim,n_iter)
!
IF ( iter .lt. 1 ) CALL errore('mix_broyden','n_iter is smaller than 1',1)
IF ( n_iter .gt. maxter ) CALL errore('mix_broyden','n_iter is too big',1)
IF ( ndim .le. 0 ) CALL errore('mix_broyden','ndim .le. 0',1)
IF ( iter < 1 ) CALL errore('mix_broyden','n_iter is smaller than 1',1)
IF ( n_iter > maxter ) CALL errore('mix_broyden','n_iter is too big',1)
IF ( ndim <= 0 ) CALL errore('mix_broyden','ndim <= 0',1)
!
IF ( iter .eq. 1 ) THEN
IF ( .not. ALLOCATED(df2) ) ALLOCATE( df2(nbndfs,nkfs,ndim,n_iter) )
IF ( .not. ALLOCATED(dv2) ) ALLOCATE( dv2(nbndfs,nkfs,ndim,n_iter) )
IF ( iter == 1 ) THEN
IF ( .NOT. ALLOCATED(df2) ) ALLOCATE ( df2(nbndfs,nkfs,ndim,n_iter) )
IF ( .NOT. ALLOCATED(dv2) ) ALLOCATE ( dv2(nbndfs,nkfs,ndim,n_iter) )
ENDIF
IF ( conv .OR. iter .eq. nsiter ) THEN
IF (ALLOCATED(df2)) DEALLOCATE(df2)
IF (ALLOCATED(dv2)) DEALLOCATE(dv2)
IF ( conv .OR. iter == nsiter ) THEN
IF (ALLOCATED(df2)) DEALLOCATE (df2)
IF (ALLOCATED(dv2)) DEALLOCATE (dv2)
RETURN
ENDIF
IF ( .not. ALLOCATED(deltainsave) ) ALLOCATE( deltainsave(ndim) )
IF ( .NOT. ALLOCATED(deltainsave) ) ALLOCATE ( deltainsave(ndim) )
deltainsave(:) = deltain(:)
!
! iter_used = iter-1 IF iter <= n_iter
......@@ -517,7 +517,7 @@
deltaout(n) = deltaout(n) - deltain(n)
ENDDO
!
IF ( iter .gt. 1 ) THEN
IF ( iter > 1 ) THEN
DO n = 1, ndim
df2(ibnd,ik,n,ipos) = deltaout(n) - df2(ibnd,ik,n,ipos)
dv2(ibnd,ik,n,ipos) = deltain(n) - dv2(ibnd,ik,n,ipos)
......@@ -576,7 +576,7 @@
df2(ibnd,ik,:,inext) = deltaout(:)
dv2(ibnd,ik,:,inext) = deltainsave(:)
!
IF ( ALLOCATED(deltainsave) ) DEALLOCATE(deltainsave)
IF ( ALLOCATED(deltainsave) ) DEALLOCATE (deltainsave)
!
RETURN
!
......
......@@ -45,13 +45,6 @@
INTEGER :: ierr
!! Error status
!
DEALLOCATE (inv_tau_all)
DEALLOCATE (zi_allvb)
IF (mp_mesh_k .AND. iterative_bte .AND. epmatkqread) DEALLOCATE (s_BZtoIBZ_full)
IF (mp_mesh_k .AND. iterative_bte .AND. epmatkqread) DEALLOCATE (ixkqf_tr)
IF (int_mob .AND. carrier) DEALLOCATE (inv_tau_allcb)
IF (int_mob .AND. carrier) DEALLOCATE (zi_allcb)
!
#if defined(__MPI)
IF (etf_mem == 1) then
CALL MPI_FILE_CLOSE(iunepmatwp2,ierr)
......@@ -104,7 +97,7 @@
!! Imported the noncolinear case implemented by xlzhang
!!
!----------------------------------------------------------------------
USE phcom, ONLY : alphap, dmuxc, drc, dyn, evq, dvpsi, &
USE phcom, ONLY : alphap, dmuxc, drc, dyn, dvpsi, &
int5, vlocq, int2_so, int5_so
USE lrus, ONLY : becp1, int3, int3_nc
USE phus, ONLY : int1, int1_nc, int2, int4, int4_nc
......@@ -113,7 +106,7 @@
USE control_lr, ONLY : nbnd_occ
USE becmod, ONLY : becp, deallocate_bec_type
USE elph2, ONLY : el_ph_mat, epf17, epsi, etf,&
etq, wf, wkf, wqf, &
etq, wkf, wqf, &
xkq, zstar, xkf, xqf, epmatwp, eps_rpa
USE klist_epw, ONLY : xk_all, xk_loc, xk_cryst, et_all, et_loc, &
isk_loc, isk_all
......@@ -129,54 +122,69 @@
INTEGER :: ipol
!! Polarization number
!
IF ( epwread .and. .not. epbread ) THEN
IF ( epwread .and. .NOT. epbread ) THEN
! EPW variables ONLY
!
IF(ALLOCATED(el_ph_mat)) DEALLOCATE (el_ph_mat)
IF(ALLOCATED(epmatwp)) DEALLOCATE (epmatwp)
IF(ALLOCATED(epf17)) DEALLOCATE (epf17)
IF(ALLOCATED(etq)) DEALLOCATE (etq)
IF(ALLOCATED(etf)) DEALLOCATE (etf)
IF(ALLOCATED(wf)) DEALLOCATE (wf)
IF(ALLOCATED(xkq)) DEALLOCATE (xkq)
IF(ALLOCATED(xkf)) DEALLOCATE (xkf)
IF(ALLOCATED(wkf)) DEALLOCATE (wkf)
IF(ALLOCATED(xqf)) DEALLOCATE (xqf)
IF(ALLOCATED(wqf)) DEALLOCATE (wqf)
IF(ALLOCATED(et_all)) DEALLOCATE (et_all)
IF(ALLOCATED(eps_rpa)) DEALLOCATE (eps_rpa)
DEALLOCATE (vlocq)
DEALLOCATE (dmuxc)
DEALLOCATE (eigqts)
DEALLOCATE (rtau)
DEALLOCATE (u)
DEALLOCATE (name_rap_mode)
DEALLOCATE (num_rap_mode)
DEALLOCATE (npert)
IF (ALLOCATED(alphap)) THEN
DO ik=1, nks
DO ipol=1, 3
CALL deallocate_bec_type(alphap(ipol, ik))
ENDDO
ENDDO
DEALLOCATE (alphap)
ENDIF
IF (ALLOCATED(becp1)) THEN
DO ik=1, size(becp1)
CALL deallocate_bec_type(becp1(ik))
ENDDO
DEALLOCATE (becp1)
ENDIF
CALL deallocate_bec_type(becp)
!
ELSE
!
IF(ASSOCIATED(evq)) DEALLOCATE(evq)
IF(ASSOCIATED(igkq)) DEALLOCATE(igkq)
!
IF(ALLOCATED(dvpsi)) DEALLOCATE (dvpsi)
!
IF(ALLOCATED(vlocq)) DEALLOCATE (vlocq)
IF(ALLOCATED(dmuxc)) DEALLOCATE (dmuxc)
IF(ASSOCIATED(igkq)) DEALLOCATE (igkq)
!
IF(ALLOCATED(eigqts)) DEALLOCATE (eigqts)
IF(ALLOCATED(rtau)) DEALLOCATE (rtau)
IF(ASSOCIATED(u)) DEALLOCATE (u)
if(allocated(name_rap_mode)) deallocate (name_rap_mode)
if(allocated(num_rap_mode)) deallocate (num_rap_mode)
DEALLOCATE (vlocq)
DEALLOCATE (dmuxc)
DEALLOCATE (eigqts)
DEALLOCATE (rtau)
DEALLOCATE (u)
DEALLOCATE (name_rap_mode)
DEALLOCATE (num_rap_mode)
IF(ALLOCATED(dyn)) DEALLOCATE (dyn)
IF(ALLOCATED(epsi)) DEALLOCATE (epsi)
IF(ALLOCATED(zstar)) DEALLOCATE (zstar)
!
IF(ALLOCATED(npert)) DEALLOCATE (npert)
DEALLOCATE (npert)
!
IF(ALLOCATED(int1)) DEALLOCATE (int1)
IF(ALLOCATED(int2)) DEALLOCATE (int2)
IF(ALLOCATED(int3)) DEALLOCATE (int3)
IF(ALLOCATED(int4)) DEALLOCATE (int4)
IF(ALLOCATED(int5)) DEALLOCATE (int5)
IF(ALLOCATED(int1_nc)) DEALLOCATE(int1_nc)
IF(ALLOCATED(int3_nc)) DEALLOCATE(int3_nc)
IF(ALLOCATED(int4_nc)) DEALLOCATE(int4_nc)
IF(ALLOCATED(int2_so)) DEALLOCATE(int2_so)
IF(ALLOCATED(int5_so)) DEALLOCATE(int5_so)
IF(ALLOCATED(int1_nc)) DEALLOCATE (int1_nc)
IF(ALLOCATED(int3_nc)) DEALLOCATE (int3_nc)
IF(ALLOCATED(int4_nc)) DEALLOCATE (int4_nc)
IF(ALLOCATED(int2_so)) DEALLOCATE (int2_so)
IF(ALLOCATED(int5_so)) DEALLOCATE (int5_so)
!
IF (allocated(alphap)) THEN
DO ik = 1, nks
......@@ -184,29 +192,26 @@
CALL deallocate_bec_type( alphap(ipol,ik) )
ENDDO
ENDDO
DEALLOCATE(alphap)
DEALLOCATE (alphap)
ENDIF
IF (allocated(becp1)) THEN
DO ik = 1, size(becp1)
CALL deallocate_bec_type( becp1(ik) )
ENDDO
DEALLOCATE(becp1)
DEALLOCATE (becp1)
ENDIF
CALL deallocate_bec_type ( becp )
!
IF(ALLOCATED(nbnd_occ)) DEALLOCATE(nbnd_occ)
IF(ALLOCATED(m_loc)) DEALLOCATE(m_loc)
IF(ALLOCATED(nbnd_occ)) DEALLOCATE (nbnd_occ)
IF(ALLOCATED(m_loc)) DEALLOCATE (m_loc)
!
IF(ALLOCATED(drc)) DEALLOCATE(drc)
IF(ALLOCATED(drc)) DEALLOCATE (drc)
!
! EPW variables
!
IF(ALLOCATED(el_ph_mat)) DEALLOCATE (el_ph_mat)
IF(ALLOCATED(epmatwp)) DEALLOCATE (epmatwp)
IF(ALLOCATED(epf17)) DEALLOCATE (epf17)
IF(ALLOCATED(etq)) DEALLOCATE (etq)
IF(ALLOCATED(etf)) DEALLOCATE (etf)
IF(ALLOCATED(wf)) DEALLOCATE (wf)
IF(ALLOCATED(xkq)) DEALLOCATE (xkq)
IF(ALLOCATED(xkf)) DEALLOCATE (xkf)
IF(ALLOCATED(wkf)) DEALLOCATE (wkf)
......@@ -219,8 +224,7 @@
IF(ALLOCATED(et_loc)) DEALLOCATE (et_loc)
IF(ALLOCATED(isk_loc)) DEALLOCATE (isk_loc)
IF(ALLOCATED(isk_all)) DEALLOCATE (isk_all)
IF(ALLOCATED(eps_rpa)) DEALLOCATE (eps_rpa)
ENDIF ! epwread .and. .not. epbread
ENDIF ! epwread .and. .NOT. epbread
!
END SUBROUTINE deallocate_epw
! ---------------------------------------------------------------
......
......@@ -34,6 +34,7 @@
!
REAL(DP), PARAMETER :: ang2cm = 1.0d-8
REAL(DP), PARAMETER :: ang2m = 1.0d-10
REAL(DP), PARAMETER :: cm2m = 1.0d-2
REAL(DP), PARAMETER :: bohr = 0.52917721092d0
REAL(DP), PARAMETER :: ryd2mev = 13605.6981d0
REAL(DP), PARAMETER :: ryd2ev = 13.6056981d0
......@@ -62,6 +63,7 @@
REAL(DP), PARAMETER :: eps12 = 1.0E-12_DP
REAL(DP), PARAMETER :: eps14 = 1.0E-14_DP
REAL(DP), PARAMETER :: eps16 = 1.0E-16_DP
REAL(DP), PARAMETER :: eps20 = 1.0E-20_DP
REAL(DP), PARAMETER :: eps24 = 1.0E-24_DP
REAL(DP), PARAMETER :: eps32 = 1.0E-32_DP
REAL(DP), PARAMETER :: eps80 = 1.0E-80_DP
......
This diff is collapsed.
......@@ -105,7 +105,7 @@
COMPLEX(kind=DP), POINTER :: qgmq(:)
!! the augmentation function at q+G
!
IF (.not.okvan) RETURN
IF ( .NOT. okvan) RETURN
!
CALL start_clock('dvanqq2')
!
......@@ -113,17 +113,17 @@
int2(:,:,:,:,:) = czero
int4(:,:,:,:,:) = czero
int5(:,:,:,:,:) = czero
ALLOCATE( sk(ngm) )
ALLOCATE( aux1(ngm) )
ALLOCATE( aux2(ngm) )
ALLOCATE( aux3(ngm) )
ALLOCATE( aux5(ngm) )
ALLOCATE( qmodg(ngm) )
ALLOCATE( qmod(ngm) )
ALLOCATE( qgmq(ngm) )
ALLOCATE( qgm(ngm))
ALLOCATE( ylmk0(ngm, lmaxq * lmaxq) )
ALLOCATE( ylmkq(ngm, lmaxq * lmaxq) )
ALLOCATE ( sk(ngm) )
ALLOCATE ( aux1(ngm) )
ALLOCATE ( aux2(ngm) )
ALLOCATE ( aux3(ngm) )
ALLOCATE ( aux5(ngm) )
ALLOCATE ( qmodg(ngm) )
ALLOCATE ( qmod(ngm) )
ALLOCATE ( qgmq(ngm) )
ALLOCATE ( qgm(ngm))
ALLOCATE ( ylmk0(ngm, lmaxq * lmaxq) )
ALLOCATE ( ylmkq(ngm, lmaxq * lmaxq) )
!
! compute spherical harmonics
!
......@@ -132,10 +132,10 @@
qmodg(ig) = sqrt( gg(ig) )
ENDDO
!
ALLOCATE( qpg(3, ngm) )
ALLOCATE ( qpg(3, ngm) )
CALL setqmod( ngm, xq, g, qmod, qpg )
CALL ylmr2(lmaxq * lmaxq, ngm, qpg, qmod, ylmkq)
DEALLOCATE(qpg)
DEALLOCATE (qpg)
DO ig = 1, ngm
qmod(ig) = sqrt( qmod(ig) )
ENDDO
......@@ -314,18 +314,18 @@
!SUM((REAL(REAL(int5(:,:,:,:,:))))**2)+SUM((REAL(AIMAG(int5(:,:,:,:,:))))**2)
!END
!
DEALLOCATE(sk)
DEALLOCATE(aux1)
DEALLOCATE(aux2)
DEALLOCATE(aux3)
DEALLOCATE(aux5)
DEALLOCATE(qmodg)
DEALLOCATE(qmod)
DEALLOCATE(qgmq)
DEALLOCATE(qgm)
DEALLOCATE(ylmk0)
DEALLOCATE(ylmkq)
DEALLOCATE(veff)
DEALLOCATE (sk)
DEALLOCATE (aux1)
DEALLOCATE (aux2)
DEALLOCATE (aux3)
DEALLOCATE (aux5)
DEALLOCATE (qmodg)
DEALLOCATE (qmod)
DEALLOCATE (qgmq)
DEALLOCATE (qgm)
DEALLOCATE (ylmk0)
DEALLOCATE (ylmkq)
DEALLOCATE (veff)
!
CALL stop_clock ('dvanqq2')
RETURN
......
......@@ -10,7 +10,7 @@
! adapted from PH/dvqpsi_us (QE)
!
!----------------------------------------------------------------------
SUBROUTINE dvqpsi_us3( ik, uact, addnlcc, xxkq, xq0 )
SUBROUTINE dvqpsi_us3( ik, uact, addnlcc, xxkq, xq0, igk, igkq, npw, npwq )
!----------------------------------------------------------------------
!!
!! This routine calculates dV_bare/dtau * psi for one perturbation
......@@ -42,12 +42,12 @@
USE nlcc_ph, ONLY : drc
USE uspp, ONLY : nlcc_any
USE eqv, ONLY : dvpsi, dmuxc, vlocq
USE qpoint, ONLY : eigqts, npwq
USE qpoint, ONLY : eigqts
USE klist, ONLY : ngk
USE klist_epw, ONLY : isk_loc
USE gc_lr, ONLY : grho, dvxc_rr, dvxc_sr, dvxc_ss, dvxc_s
USE funct, ONLY : dft_is_gradient, dft_is_nonlocc
USE elph2, ONLY : igkq, igk, lower_band, upper_band
USE elph2, ONLY : lower_band, upper_band
USE constants_epw, ONLY : czero, eps12
!
IMPLICIT NONE
......@@ -57,6 +57,14 @@
!
INTEGER, INTENT(in) :: ik
!! Counter on k-point
INTEGER, INTENT(in) :: npw
!! Number of k+G-vectors inside 'ecut sphere'
INTEGER, INTENT(in) :: npwq
!! Number of k+G-vectors inside 'ecut sphere'
INTEGER, INTENT(in) :: igk(npw)
!! k+G mapping
INTEGER, INTENT(in) :: igkq(npwq)
!! k+G+q mapping
!
REAL(kind=DP), INTENT (in) :: xq0(3)
!! Current coarse q-point coordinate
......@@ -84,8 +92,6 @@
!! counter on spin
INTEGER :: ip
!! counter on polarizations
INTEGER :: npw
!! Number of k+G-vectors inside 'ecut sphere'
!
REAL(kind=DP) :: fac
!! spin degeneracy factor
......@@ -108,15 +114,13 @@
!
CALL start_clock('dvqpsi_us3')
!
npw = ngk(ik)
!
IF (nlcc_any .AND. addnlcc) THEN
ALLOCATE( drhoc(dfftp%nnr) )
ALLOCATE( aux (dfftp%nnr) )
ALLOCATE( auxs(dffts%nnr) )
ALLOCATE (drhoc(dfftp%nnr))
ALLOCATE (aux(dfftp%nnr))
ALLOCATE (auxs(dffts%nnr))
ENDIF
ALLOCATE( aux1(dffts%nnr) )
ALLOCATE( aux2(dffts%nnr) )
ALLOCATE (aux1(dffts%nnr))
ALLOCATE (aux2(dffts%nnr))
!
! We start by computing the contribution of the local potential.
! The computation of the derivative of the local potential is done in
......@@ -131,7 +135,7 @@
u1 = uact(mu+1)
u2 = uact(mu+2)
u3 = uact(mu+3)
IF (abs(u1) + abs(u2) + abs(u3) .gt. eps12) THEN
IF (abs(u1) + abs(u2) + abs(u3) > eps12) THEN
nt = ityp(na)
gu0 = xq0(1) * u1 + xq0(2) * u2 + xq0(3) * u3
DO ig = 1, ngms
......@@ -154,7 +158,7 @@
u1 = uact(mu+1)
u2 = uact(mu+2)
u3 = uact(mu+3)
IF (abs(u1) + abs(u2) + abs(u3) .gt. eps12) THEN
IF (abs(u1) + abs(u2) + abs(u3) > eps12) THEN
nt = ityp(na)
gu0 = xq0(1) * u1 + xq0(2) * u2 + xq0(3) * u3
IF (upf(nt)%nlcc) THEN
......@@ -171,7 +175,7 @@
!
CALL invfft('Rho', drhoc, dfftp)
!
IF (.not.lsda) THEN
IF ( .NOT. lsda) THEN
DO ir = 1, dfftp%nnr
aux(ir) = drhoc(ir) * dmuxc(ir,1,1)
ENDDO
......@@ -250,18 +254,18 @@
ENDDO
!
IF (nlcc_any .AND. addnlcc) THEN
DEALLOCATE(drhoc)
DEALLOCATE(aux)
DEALLOCATE(auxs)
DEALLOCATE (drhoc)
DEALLOCATE (aux)
DEALLOCATE (auxs)
ENDIF
DEALLOCATE(aux1)
DEALLOCATE(aux2)
DEALLOCATE (aux1)
DEALLOCATE (aux2)
!
! We add the contribution of the nonlocal potential in the US form
! First a term similar to the KB case.
! Then a term due to the change of the D coefficients in the perturbat
!
CALL dvqpsi_us_only3( ik, uact, xxkq )
CALL dvqpsi_us_only3(ik, uact, xxkq, igkq, npwq)
!
CALL stop_clock('dvqpsi_us3')
!
......
......@@ -10,7 +10,7 @@
! adapted from PH/dvqpsi_us_only (QE)
!
!----------------------------------------------------------------------
subroutine dvqpsi_us_only3( ik, uact, xxkq )
subroutine dvqpsi_us_only3 (ik, uact, xxkq, igkq, npwq)
!----------------------------------------------------------------------
!!
!! This routine calculates dV_bare/dtau * psi for one perturbation
......@@ -30,11 +30,10 @@
USE wvfct, ONLY : npwx, et
USE uspp, ONLY : okvan, nkb, vkb
USE uspp_param, ONLY : nh, nhm
USE qpoint, ONLY : npwq
USE phus, ONLY : int1, int1_nc, int2, int2_so, alphap
USE lrus, ONLY : becp1
USE eqv, ONLY : dvpsi
USE elph2, ONLY : igkq, lower_band, upper_band
USE elph2, ONLY : lower_band, upper_band
USE noncollin_module, ONLY : noncolin, npol
USE constants_epw, ONLY : czero, cone, eps12
USE klist_epw, ONLY : isk_loc
......@@ -43,6 +42,10 @@
!
INTEGER, INTENT(in) :: ik
!! the k point
INTEGER, INTENT(in) :: npwq
!! Number of k+G-vectors inside 'ecut sphere'
INTEGER, INTENT(in) :: igkq(npwq)
!! k+G+q mapping
REAL(kind=DP), INTENT(in) :: xxkq(3)
!! the k+q point (cartesian coordinates)
COMPLEX(kind=DP), INTENT(in) :: uact(3 * nat)
......@@ -94,15 +97,15 @@
!
CALL start_clock('dvqpsi_us_on')
IF (noncolin) THEN
ALLOCATE( ps1_nc(nkb, npol, lower_band:upper_band) )
ALLOCATE( ps2_nc(nkb, npol, lower_band:upper_band, 3) )
ALLOCATE( deff_nc(nhm, nhm, nat, nspin) )
ALLOCATE ( ps1_nc(nkb, npol, lower_band:upper_band) )
ALLOCATE ( ps2_nc(nkb, npol, lower_band:upper_band, 3) )
ALLOCATE ( deff_nc(nhm, nhm, nat, nspin) )
ELSE
ALLOCATE( ps1(nkb, lower_band:upper_band) )
ALLOCATE( ps2(nkb, lower_band:upper_band, 3) )
ALLOCATE( deff(nhm, nhm, nat) )
ALLOCATE ( ps1(nkb, lower_band:upper_band) )
ALLOCATE ( ps2(nkb, lower_band:upper_band, 3) )
ALLOCATE ( deff(nhm, nhm, nat) )
ENDIF
ALLOCATE( aux(npwx) )
ALLOCATE ( aux(npwx) )
IF (lsda) current_spin = isk_loc(ik)
!
! we first compute the coefficients of the vectors
......@@ -125,7 +128,7 @@
ijkb0 = 0
DO nt = 1, ntyp
DO na = 1, nat
IF (ityp(na) .eq. nt) THEN
IF (ityp(na) == nt) THEN
mu = 3 * (na - 1)
DO ih = 1, nh(nt)
ikb = ijkb0 + ih
......@@ -211,29 +214,29 @@
!
! This term is proportional to beta(k+q+G)
!
IF (nkb.gt.0) THEN
IF (nkb > 0) THEN
IF (noncolin) THEN
CALL zgemm( 'n', 'n', npwq, (upper_band-lower_band+1)*npol, nkb, &
cone, vkb, npwx, ps1_nc, nkb, cone, dvpsi, npwx )
CALL zgemm('n', 'n', npwq, (upper_band-lower_band+1)*npol, nkb, &
cone, vkb, npwx, ps1_nc, nkb, cone, dvpsi, npwx)
ELSE
CALL zgemm( 'n', 'n', npwq, (upper_band-lower_band+1), nkb, &
cone, vkb, npwx, ps1, nkb, cone, dvpsi, npwx )
CALL zgemm('n', 'n', npwq, (upper_band-lower_band+1), nkb, &
cone, vkb, npwx, ps1, nkb, cone, dvpsi, npwx)
ENDIF
ENDIF
!
! This term is proportional to (k+q+G)_\alpha*beta(k+q+G)
!
DO ikb = 1, nkb
DO ipol = 1, 3
DO ikb=1, nkb
DO ipol=1, 3
ok = .false.
IF (noncolin) THEN
DO ibnd = lower_band, upper_band
ok = ok .OR. ( abs( ps2_nc(ikb,1,ibnd,ipol) ) .gt. eps12 ) .OR. &
( abs( ps2_nc(ikb,2,ibnd,ipol) ) .gt. eps12 )
ok = ok .OR. (ABS(ps2_nc(ikb,1,ibnd,ipol) ) > eps12 ) .OR. &
(ABS(ps2_nc(ikb,2,ibnd,ipol) ) > eps12 )
ENDDO
ELSE
DO ibnd = lower_band, upper_band
ok = ok .OR. ( abs( ps2(ikb,ibnd,ipol) ) .gt. eps12 )
DO ibnd=lower_band, upper_band
ok = ok .OR. (ABS(ps2(ikb,ibnd,ipol) ) > eps12)
ENDDO
ENDIF
IF (ok) THEN
......@@ -254,15 +257,15 @@
ENDDO
ENDDO
!
DEALLOCATE(aux)
DEALLOCATE (aux)
IF (noncolin) THEN
DEALLOCATE(ps2_nc)
DEALLOCATE(ps1_nc)
DEALLOCATE(deff_nc)
DEALLOCATE (ps2_nc)
DEALLOCATE (ps1_nc)
DEALLOCATE (deff_nc)
ELSE
DEALLOCATE(ps2)
DEALLOCATE(ps1)
DEALLOCATE(deff)
DEALLOCATE (ps2)
DEALLOCATE (ps1)
DEALLOCATE (deff)
ENDIF
!
CALL stop_clock('dvqpsi_us_on')
......
......@@ -43,7 +43,7 @@
ENDIF
!
CALL estimate_tc_gap
IF ( gap_edge .gt. 0.d0 ) THEN
IF ( gap_edge > 0.d0 ) THEN
gap0 = gap_edge
ENDIF
IF ( lreal ) CALL eliashberg_iso_raxis
......@@ -62,13 +62,13 @@
CALL eliashberg_init
CALL evaluate_a2f_lambda
CALL estimate_tc_gap
IF ( gap_edge .gt. 0.d0 ) THEN
IF ( gap_edge > 0.d0 ) THEN
gap0 = gap_edge
ENDIF
IF ( limag ) CALL eliashberg_aniso_iaxis
ENDIF
!
IF ( .not. liso .AND. .not. laniso ) THEN
IF ( .NOT. liso .AND. .NOT. laniso ) THEN
WRITE(stdout,'(/5x,a)') REPEAT('=',67)
WRITE(stdout,'(5x,"Calculate Eliashberg spectral function")')
WRITE(stdout,'(5x,a/)') REPEAT('=',67)
......