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
2010-01-27 18:50:07 +00:00 · 2010-01-27 18:50:07 +00:00 · c122ffed47
parent 65b4006959
commit c122ffed47
13 changed files with 150 additions and 106 deletions
--- 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
--- 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
  !
--- 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
--- 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
  !
--- 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, &
--- 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 \
--- 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()
--- 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
--- 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
  !
--- 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
  !
--- 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
  !
--- 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
--- 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
  !