From c122ffed47a80153f0cb6e76f60145039e690aa5 Mon Sep 17 00:00:00 2001 From: obm Date: Wed, 27 Jan 2010 18:50:07 +0000 Subject: [PATCH] 1) For TDDFPT purposes, subroutine newd does too many things at once, and accepts too few parameters. I have splitted them and collected all in a module. Usual calls to newd is not changed, apart from necessity to include the module, newq is the part I use for calculating response charge density. 2) Some gamma only additions to PH/dv_of_drho.f90 proved to be unnecessary, removing. I am still trying to find an efficient/minimal impact way to cast this subroutine to use real instead of complex input. As usual, I have tested before posting, however be sure to check before in your applications. git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6322 c92efa57-630b-4861-b058-cf58834340f0 --- GIPAW/gipaw_module.f90 | 1 + Gamma/cg_setup.f90 | 1 + PH/dv_of_drho.f90 | 21 +--- PP/do_initial_state.f90 | 1 + PP/pw2casino.f90 | 1 + PW/Makefile | 2 +- PW/compute_fes_grads.f90 | 1 + PW/electrons.f90 | 1 + PW/hinit1.f90 | 1 + PW/init_run.f90 | 1 + PW/move_ions.f90 | 1 + PW/newd.f90 | 223 +++++++++++++++++++++++---------------- PW/read_file.f90 | 1 + 13 files changed, 150 insertions(+), 106 deletions(-) diff --git a/GIPAW/gipaw_module.f90 b/GIPAW/gipaw_module.f90 index f4ab43d72..0727542a1 100644 --- a/GIPAW/gipaw_module.f90 +++ b/GIPAW/gipaw_module.f90 @@ -341,6 +341,7 @@ CONTAINS USE uspp_param, ONLY : upf USE mp_global, ONLY : inter_pool_comm USE mp, ONLY : mp_max, mp_min + USE dfunct, only : newd IMPLICIT none integer :: ik, nt, ibnd diff --git a/Gamma/cg_setup.f90 b/Gamma/cg_setup.f90 index 8d7eb18b9..a34a7c221 100644 --- a/Gamma/cg_setup.f90 +++ b/Gamma/cg_setup.f90 @@ -21,6 +21,7 @@ SUBROUTINE cg_setup USE io_files, ONLY: prefix, iunpun, iunres USE cgcom USE funct, only : dft_is_gradient, dmxc + USE dfunct, only : newd ! IMPLICIT NONE ! diff --git a/PH/dv_of_drho.f90 b/PH/dv_of_drho.f90 index 74a507b61..a49b7e129 100644 --- a/PH/dv_of_drho.f90 +++ b/PH/dv_of_drho.f90 @@ -55,14 +55,12 @@ subroutine dv_of_drho (mode, dvscf, flag) complex(DP), allocatable :: dvhart (:,:) !required in gamma_only call start_clock ('dv_of_drho') - allocate (dvaux( nrxx, nspin_mag)) - + dvaux (:,:) = (0.d0, 0.d0) if (flag) allocate (drhoc( nrxx)) ! ! the exchange-correlation contribution is computed in real space ! - dvaux (:,:) = (0.d0, 0.d0) if (lrpa) goto 111 fac = 1.d0 / DBLE (nspin_lsda) if (nlcc_any.and.flag) then @@ -72,16 +70,6 @@ subroutine dv_of_drho (mode, dvscf, flag) dvscf(:, is) = dvscf(:, is) + fac * drhoc (:) enddo endif - !OBM:not totally convinced of the necessity of the following change for gamma point, inquire. - if (gamma_only) then - do is = 1, nspin_mag - do is1 = 1, nspin_mag - do ir = 1, nrxx - dvaux(ir,is) = dvaux(ir,is) + cmplx(DBLE(dmuxc(ir,is,is1) * dvscf(ir,is1)),0.d0,dp) - enddo - enddo - enddo - else do is = 1, nspin_mag do is1 = 1, nspin_mag do ir = 1, nrxx @@ -90,7 +78,6 @@ subroutine dv_of_drho (mode, dvscf, flag) enddo enddo - endif ! ! add gradient correction to xc, NB: if nlcc is true we need to add here @@ -106,6 +93,8 @@ subroutine dv_of_drho (mode, dvscf, flag) dvscf(:, is) = dvscf(:, is) - fac * drhoc (:) enddo endif + + 111 continue ! ! copy the total (up+down) delta rho in dvscf(*,1) and go to G-space @@ -136,9 +125,9 @@ subroutine dv_of_drho (mode, dvscf, flag) enddo ! ! at the end the two contributes are added - dvscf (:,:) = dvaux (:,:) + dvhart(:,:) + dvscf = dvaux + dvhart !OBM : Again not totally convinced about this trimming. - dvscf (:,:) = cmplx(DBLE(dvscf(:,:)),0.0d0,dp) + !dvscf (:,:) = cmplx(DBLE(dvscf(:,:)),0.0d0,dp) deallocate(dvhart) else do is = 1, nspin_lsda diff --git a/PP/do_initial_state.f90 b/PP/do_initial_state.f90 index 24dce2b8b..849febdbf 100644 --- a/PP/do_initial_state.f90 +++ b/PP/do_initial_state.f90 @@ -39,6 +39,7 @@ SUBROUTINE do_initial_state (excite) USE ener, ONLY : ef USE parameters, ONLY : ntypx USE control_flags, ONLY: gamma_only + USE DFUNCT, ONLY : newd ! IMPLICIT NONE ! diff --git a/PP/pw2casino.f90 b/PP/pw2casino.f90 index 4d3892f9a..c1ffacceb 100644 --- a/PP/pw2casino.f90 +++ b/PP/pw2casino.f90 @@ -109,6 +109,7 @@ SUBROUTINE compute_casino USE funct, ONLY : dft_is_meta USE mp_global, ONLY: inter_pool_comm, intra_pool_comm USE mp, ONLY: mp_sum + USE dfunct, only : newd IMPLICIT NONE INTEGER :: ig, ibnd, ik, io, na, j, ispin, nbndup, nbnddown, & diff --git a/PW/Makefile b/PW/Makefile index 121acb8c3..fe75a4ae6 100644 --- a/PW/Makefile +++ b/PW/Makefile @@ -23,6 +23,7 @@ atomic_rho.o \ atomic_wfc.o \ average_pp.o \ becmod.o \ +newd.o \ bp_c_phase.o \ bp_calc_btq.o \ bp_qvan3.o \ @@ -130,7 +131,6 @@ n_plane_waves.o \ new_ns.o \ new_occ.o \ ns_adj.o \ -newd.o \ noncol.o \ non_scf.o \ offset_atom_wfc.o \ diff --git a/PW/compute_fes_grads.f90 b/PW/compute_fes_grads.f90 index 29cbcb99f..669e757b6 100644 --- a/PW/compute_fes_grads.f90 +++ b/PW/compute_fes_grads.f90 @@ -698,6 +698,7 @@ END SUBROUTINE metadyn SUBROUTINE reset_init_mag() !---------------------------------------------------------------------------- ! + USE dfunct, only : newd IMPLICIT NONE ! CALL hinit0() diff --git a/PW/electrons.f90 b/PW/electrons.f90 index b6d95a559..d45c8207e 100644 --- a/PW/electrons.f90 +++ b/PW/electrons.f90 @@ -78,6 +78,7 @@ SUBROUTINE electrons() n_charge_compensation, & vloc_of_g_zero, mr1, mr2, mr3, & vcomp, comp_thr, icomp + USE dfunct, only : newd ! ! IMPLICIT NONE diff --git a/PW/hinit1.f90 b/PW/hinit1.f90 index e2282317f..a7a5ce5d6 100644 --- a/PW/hinit1.f90 +++ b/PW/hinit1.f90 @@ -22,6 +22,7 @@ SUBROUTINE hinit1() USE realus, ONLY : qpointlist USE wannier_new, ONLY : use_wannier USE martyna_tuckerman, ONLY : tag_wg_corr_as_obsolete + USE dfunct, only : newd ! IMPLICIT NONE ! diff --git a/PW/init_run.f90 b/PW/init_run.f90 index 1fde2d957..08cdd49b7 100644 --- a/PW/init_run.f90 +++ b/PW/init_run.f90 @@ -25,6 +25,7 @@ SUBROUTINE init_run() USE ee_mod, ONLY : do_comp, do_coarse ! Wannier_ac USE wannier_new, ONLY : use_wannier + USE dfunct, only : newd ! IMPLICIT NONE ! diff --git a/PW/move_ions.f90 b/PW/move_ions.f90 index e554ee94d..a880eed30 100644 --- a/PW/move_ions.f90 +++ b/PW/move_ions.f90 @@ -41,6 +41,7 @@ SUBROUTINE move_ions() USE bfgs_module, ONLY : bfgs, terminate_bfgs USE basic_algebra_routines, ONLY : norm USE dynamics_module, ONLY : verlet, langevin_md, proj_verlet + USE dfunct, only : newd ! IMPLICIT NONE ! diff --git a/PW/newd.f90 b/PW/newd.f90 index 073fca32a..8be07c0ee 100644 --- a/PW/newd.f90 +++ b/PW/newd.f90 @@ -6,71 +6,44 @@ ! or http://www.gnu.org/copyleft/gpl.txt . ! ! -SUBROUTINE newd() - USE uspp, ONLY : deeq - USE uspp_param, ONLY : upf, nh - USE lsda_mod, ONLY : nspin - USE ions_base, ONLY : nat, ntyp => nsp, ityp - USE realus, ONLY : newd_r - USE control_flags, ONLY : tqr - USE paw_variables, ONLY : okpaw, ddd_paw - IMPLICIT NONE - integer :: na, nt, ih, jh, ijh - if (tqr) then - call newd_r() - else - call newd_g() - end if - if (okpaw) then - ! Add paw contributions to deeq (computed in paw_potential) - do na=1,nat - nt = ityp(na) - IF (.not.upf(nt)%tpawp) cycle - ijh=0 - do ih=1,nh(nt) - do jh=ih,nh(nt) - ijh=ijh+1 - deeq(ih,jh,na,1:nspin) = deeq(ih,jh,na,1:nspin) & - + ddd_paw(ijh,na,1:nspin) - deeq(jh,ih,na,1:nspin) = deeq(ih,jh,na,1:nspin) - end do - end do - end do - end IF +MODULE DFUNCT - return -END SUBROUTINE newd -!---------------------------------------------------------------------------- -SUBROUTINE newd_g() - !---------------------------------------------------------------------------- +CONTAINS +!--------------------------------------- +SUBROUTINE newq(vr,deeq,skip_vltot) ! - ! ... This routine computes the integral of the effective potential with - ! ... the Q function and adds it to the bare ionic D term which is used - ! ... to compute the non-local term in the US scheme. + ! This routine computes the integral of the perturbed potential with + ! the Q function ! + USE kinds, ONLY : DP USE ions_base, ONLY : nat, ntyp => nsp, ityp USE cell_base, ONLY : omega USE gvect, ONLY : nr1, nr2, nr3, nrx1, nrx2, nrx3, & g, gg, ngm, gstart, ig1, ig2, ig3, & - eigts1, eigts2, eigts3, nl + eigts1, eigts2, eigts3, nl, nrxx USE lsda_mod, ONLY : nspin - USE scf, ONLY : v, vltot - USE uspp, ONLY : deeq, dvan, deeq_nc, dvan_so, okvan, indv + USE scf, ONLY : vltot + USE uspp, ONLY : okvan, indv USE uspp_param, ONLY : upf, lmaxq, nh, nhm USE control_flags, ONLY : gamma_only USE wavefunctions_module, ONLY : psic USE spin_orb, ONLY : lspinorb, domag - USE noncollin_module, ONLY : noncolin, nspin_mag + USE noncollin_module, ONLY : nspin_mag USE mp_global, ONLY : intra_pool_comm USE mp, ONLY : mp_sum - USE uspp, ONLY : nhtol, nhtolm ! IMPLICIT NONE ! + ! + ! Input: potential , output: contribution to integral + REAL(kind=dp), intent(in) :: vr(nrxx,nspin) + REAL(kind=dp), intent(out) :: deeq( nhm, nhm, nat, nspin ) + LOGICAL, intent(in) :: skip_vltot !If .false. vltot is added to vr when necessary + ! INTERNAL INTEGER :: ig, nt, ih, jh, na, is, nht, nb, mb - ! counters on g vectors, atom type, beta functions x 2, - ! atoms, spin, aux, aux, beta func x2 (again) + ! counters on g vectors, atom type, beta functions x 2, + ! atoms, spin, aux, aux, beta func x2 (again) #ifdef __OPENMP INTEGER :: mytid, ntids, omp_get_thread_num, omp_get_num_threads #endif @@ -80,46 +53,7 @@ SUBROUTINE newd_g() REAL(DP), ALLOCATABLE :: ylmk0(:,:), qmod(:) ! spherical harmonics, modulus of G REAL(DP) :: fact, ddot - ! - ! - IF ( .NOT. okvan ) THEN - ! - ! ... no ultrasoft potentials: use bare coefficients for projectors - ! - DO na = 1, nat - ! - nt = ityp(na) - nht = nh(nt) - ! - IF ( lspinorb ) THEN - ! - deeq_nc(1:nht,1:nht,na,1:nspin) = dvan_so(1:nht,1:nht,1:nspin,nt) - ! - ELSE IF ( noncolin ) THEN - ! - deeq_nc(1:nht,1:nht,na,1) = dvan(1:nht,1:nht,nt) - deeq_nc(1:nht,1:nht,na,2) = ( 0.D0, 0.D0 ) - deeq_nc(1:nht,1:nht,na,3) = ( 0.D0, 0.D0 ) - deeq_nc(1:nht,1:nht,na,4) = dvan(1:nht,1:nht,nt) - ! - ELSE - ! - DO is = 1, nspin - ! - deeq(1:nht,1:nht,na,is) = dvan(1:nht,1:nht,nt) - ! - END DO - ! - END IF - ! - END DO - ! - ! ... early return - ! - RETURN - ! - END IF - ! + IF ( gamma_only ) THEN ! fact = 2.D0 @@ -145,13 +79,13 @@ SUBROUTINE newd_g() ! DO is = 1, nspin_mag ! - IF ( nspin_mag == 4 .AND. is /= 1 ) THEN + IF ( (nspin_mag == 4 .AND. is /= 1) .or. skip_vltot ) THEN ! - psic(:) = v%of_r(:,is) + psic(:) = vr(:,is) ! ELSE ! - psic(:) = vltot(:) + v%of_r(:,is) + psic(:) = vltot(:) + vr(:,is) ! END IF ! @@ -239,6 +173,115 @@ SUBROUTINE newd_g() CALL mp_sum( deeq( :, :, :, 1:nspin_mag ), intra_pool_comm ) ! DEALLOCATE( aux, qgm, qmod, ylmk0 ) + +END SUBROUTINE newq +!--------------------------------------- +SUBROUTINE add_paw_to_deeq(deeq) + ! Add paw contributions to deeq (computed in paw_potential) + USE kinds, ONLY : DP + USE ions_base, ONLY : nat, ntyp => nsp, ityp + USE uspp_param, ONLY : upf, nh, nhm + USE paw_variables, ONLY : okpaw, ddd_paw + USE lsda_mod, ONLY : nspin + IMPLICIT NONE + integer :: na, nt, ih, jh, ijh + REAL(kind=dp), intent(inout) :: deeq( nhm, nhm, nat, nspin ) + + if (okpaw) then + do na=1,nat + nt = ityp(na) + IF (.not.upf(nt)%tpawp) cycle + ijh=0 + do ih=1,nh(nt) + do jh=ih,nh(nt) + ijh=ijh+1 + deeq(ih,jh,na,1:nspin) = deeq(ih,jh,na,1:nspin) & + + ddd_paw(ijh,na,1:nspin) + deeq(jh,ih,na,1:nspin) = deeq(ih,jh,na,1:nspin) + end do + end do + end do + end IF + +END SUBROUTINE add_paw_to_deeq +!--------------------------------------- +SUBROUTINE newd() + USE uspp, ONLY : deeq + USE realus, ONLY : newd_r + USE control_flags, ONLY : tqr + IMPLICIT NONE + if (tqr) then + call newd_r() + else + call newd_g() + end if + call add_paw_to_deeq(deeq) + return +END SUBROUTINE newd +!---------------------------------------------------------------------------- +SUBROUTINE newd_g() + !---------------------------------------------------------------------------- + ! + ! ... This routine computes the integral of the effective potential with + ! ... the Q function and adds it to the bare ionic D term which is used + ! ... to compute the non-local term in the US scheme. + ! + USE kinds, ONLY : DP + USE ions_base, ONLY : nat, ntyp => nsp, ityp + USE lsda_mod, ONLY : nspin + USE uspp, ONLY : deeq, dvan, deeq_nc, dvan_so, okvan, indv + USE uspp_param, ONLY : upf, lmaxq, nh, nhm + USE spin_orb, ONLY : lspinorb, domag + USE noncollin_module, ONLY : noncolin, nspin_mag + USE uspp, ONLY : nhtol, nhtolm + USE scf, ONLY : v + ! + IMPLICIT NONE + ! + INTEGER :: ig, nt, ih, jh, na, is, nht, nb, mb + ! counters on g vectors, atom type, beta functions x 2, + ! atoms, spin, aux, aux, beta func x2 (again) + ! + ! + IF ( .NOT. okvan ) THEN + ! + ! ... no ultrasoft potentials: use bare coefficients for projectors + ! + DO na = 1, nat + ! + nt = ityp(na) + nht = nh(nt) + ! + IF ( lspinorb ) THEN + ! + deeq_nc(1:nht,1:nht,na,1:nspin) = dvan_so(1:nht,1:nht,1:nspin,nt) + ! + ELSE IF ( noncolin ) THEN + ! + deeq_nc(1:nht,1:nht,na,1) = dvan(1:nht,1:nht,nt) + deeq_nc(1:nht,1:nht,na,2) = ( 0.D0, 0.D0 ) + deeq_nc(1:nht,1:nht,na,3) = ( 0.D0, 0.D0 ) + deeq_nc(1:nht,1:nht,na,4) = dvan(1:nht,1:nht,nt) + ! + ELSE + ! + DO is = 1, nspin + ! + deeq(1:nht,1:nht,na,is) = dvan(1:nht,1:nht,nt) + ! + END DO + ! + END IF + ! + END DO + ! + ! ... early return + ! + RETURN + ! + END IF + ! + call newq(v%of_r,deeq,.false.) ! atoms : & DO na = 1, nat @@ -413,3 +456,5 @@ SUBROUTINE newd_g() END SUBROUTINE newd_nc ! END SUBROUTINE newd_g + +END MODULE DFUNCT diff --git a/PW/read_file.f90 b/PW/read_file.f90 index bee7dfdd0..9abb7580a 100644 --- a/PW/read_file.f90 +++ b/PW/read_file.f90 @@ -48,6 +48,7 @@ SUBROUTINE read_file() USE ldaU, ONLY : lda_plus_u, eth, oatwfc USE realus, ONLY : qpointlist,betapointlist,init_realspace_vars,real_space USE io_global, ONLY : stdout + USE dfunct, only : newd ! IMPLICIT NONE !