Index: /trunk/src/math.F90 =================================================================== --- /trunk/src/math.F90 (revision 2903) +++ /trunk/src/math.F90 (revision 2910) @@ -59,4 +59,5 @@ set_app_threshold, & operator(.app.), & + math_xor, & ddelta @@ -750,4 +751,13 @@ end function ddelta + logical function math_xor(a, b) + logical, intent(in) :: a + logical, intent(in) :: b + + math_xor = ( a .or. b ) + if ( a .and. b ) math_xor = .false. + + end function math_xor + #include "undef.F90" #include "complex.F90" Index: /trunk/src/em_resp_calc_inc.F90 =================================================================== --- /trunk/src/em_resp_calc_inc.F90 (revision 2806) +++ /trunk/src/em_resp_calc_inc.F90 (revision 2910) @@ -329,7 +329,12 @@ do idir = 1, ndim do idim = 1, sys%st%d%dim - - call X(sternheimer_calc_hvar)(sh, sys, h, lr(idir, :, ifreq), 2, sys%gr%m%x(:, idir), & - hvar(:, :, :, idim, idir, ifreq)) + + do isigma = 1, 2 + do ispin = 1, sys%st%d%nspin + hvar(1:np, ispin, isigma, idim, idir, ifreq) = sys%gr%m%x(1:np, idir) + end do + end do + + call X(sternheimer_calc_hvar)(sh, sys, h, lr(idir, :, ifreq), 2, hvar(:, :, :, idim, idir, ifreq)) end do !idim Index: /trunk/src/sternheimer_inc.F90 =================================================================== --- /trunk/src/sternheimer_inc.F90 (revision 2781) +++ /trunk/src/sternheimer_inc.F90 (revision 2910) @@ -26,7 +26,7 @@ subroutine X(sternheimer_solve)(& - this, sys, h, lr, nsigma, omega, vext, & + this, sys, h, lr, nsigma, omega, & restart_dir, rho_tag, wfs_tag, & - have_restart_rho) + have_restart_rho, vext, vext_psi) type(sternheimer_t), intent(inout) :: this type(system_t), target, intent(inout) :: sys @@ -35,8 +35,9 @@ integer, intent(in) :: nsigma R_TYPE, intent(in) :: omega - FLOAT, intent(in) :: vext(:) character(len=*), intent(in) :: restart_dir character(len=*), intent(in) :: rho_tag character(len=*), intent(in) :: wfs_tag + FLOAT, optional, intent(in) :: vext(:) + R_TYPE, optional, intent(in) :: vext_psi(:,:,:,:) logical, optional, intent(in) :: have_restart_rho @@ -61,5 +62,6 @@ ASSERT( nsigma==1 .or. nsigma ==2 ) - + ASSERT( math_xor(present(vext), present(vext_psi)) ) + m => sys%gr%m st => sys%st @@ -106,5 +108,15 @@ dl_rhoin(1:m%np, 1:st%d%nspin, 1) = lr(1)%X(dl_rho)(1:m%np, 1:st%d%nspin) - call X(sternheimer_calc_hvar)(this, sys, h, lr, nsigma, vext, hvar) + do sigma = 1, nsigma + do ik = 1, st%d%nspin + if(present(vext)) then + hvar(1:m%np, ik, sigma) = vext(1:m%np) + else + hvar(1:m%np, ik, sigma) = M_ZERO + end if + end do + end do + + call X(sternheimer_calc_hvar)(this, sys, h, lr, nsigma, hvar) do ik = 1, st%d%nspin @@ -116,4 +128,8 @@ !calculate the RHS of the Sternheimer eq Y(1:m%np, 1, sigma) = -hvar(1:m%np, ik, sigma)*st%X(psi)(1:m%np, 1, ist, ik) + + if(present(vext_psi)) then + Y(1:m%np, 1, sigma) = Y(1:m%np, 1, sigma) - vext_psi(1:m%np, 1, ist, ik) + end if !and project it into the unoccupied states @@ -236,5 +252,5 @@ end subroutine X(sternheimer_solve) -subroutine X(sternheimer_calc_hvar)(this, sys, h, lr, nsigma, vext, hvar) +subroutine X(sternheimer_calc_hvar)(this, sys, h, lr, nsigma, hvar) type(sternheimer_t), intent(inout) :: this type(system_t), target, intent(inout) :: sys @@ -242,5 +258,4 @@ type(lr_t), intent(inout) :: lr(:) integer, intent(in) :: nsigma - FLOAT, intent(in) :: vext(:) R_TYPE, intent(out) :: hvar(:,:,:) @@ -268,10 +283,6 @@ do ik = 1, sys%st%d%nspin - !* Vext - hvar(1:np, ik, 1) = vext(1:np) - !* hartree - if (this%add_hartree) hvar(1:np, ik, 1) = hvar(1:np, ik, 1) & - + hartree(1:np) + if (this%add_hartree) hvar(1:np, ik, 1) = hvar(1:np, ik, 1) + hartree(1:np) !* fxc @@ -279,6 +290,5 @@ if(.not. this%oep_kernel) then do ik2 = 1, sys%st%d%nspin - hvar(1:np, ik, 1) = hvar(1:np, ik, 1) + & - this%fxc(1:np, ik, ik2)*lr(1)%X(dl_rho)(1:np, ik2) + hvar(1:np, ik, 1) = hvar(1:np, ik, 1) + this%fxc(1:np, ik, ik2)*lr(1)%X(dl_rho)(1:np, ik2) end do else Index: /trunk/src/phonons_lr.F90 =================================================================== --- /trunk/src/phonons_lr.F90 (revision 2908) +++ /trunk/src/phonons_lr.F90 (revision 2910) @@ -120,6 +120,6 @@ do idir = 1, sys%NDIM - call dsternheimer_solve(sh, sys, h, lr, 1, M_ZERO, dv(1:sys%NP, idir, iatom), & - RESTART_DIR, phn_rho_tag(iatom, idir), phn_wfs_tag(iatom, idir)) + call dsternheimer_solve(sh, sys, h, lr, 1, M_ZERO, RESTART_DIR, phn_rho_tag(iatom, idir), phn_wfs_tag(iatom, idir), & + vext= dv(1:sys%NP, idir, iatom)) do jatom = 1, sys%geo%natoms Index: /trunk/src/vdW.F90 =================================================================== --- /trunk/src/vdW.F90 (revision 2799) +++ /trunk/src/vdW.F90 (revision 2910) @@ -243,6 +243,6 @@ call write_info(1) - call zsternheimer_solve(sh, sys, h, lr(dir, :), 1, omega, sys%gr%m%x(:,dir), & - RESTART_DIR, em_rho_tag(real(omega),dir), em_wfs_tag(dir,1)) + call zsternheimer_solve(sh, sys, h, lr(dir, :), 1, omega, & + RESTART_DIR, em_rho_tag(real(omega),dir), em_wfs_tag(dir,1), vext=sys%gr%m%x(:,dir)) end do Index: /trunk/src/em_resp.F90 =================================================================== --- /trunk/src/em_resp.F90 (revision 2908) +++ /trunk/src/em_resp.F90 (revision 2910) @@ -246,13 +246,15 @@ call zsternheimer_solve(sh, sys, h, em_vars%lr(dir, :, ifactor), 2 , & em_vars%freq_factor(ifactor)*em_vars%omega(iomega) + M_zI * em_vars%eta, & - sys%gr%m%x(:, dir), RESTART_DIR,& + RESTART_DIR,& em_rho_tag(em_vars%freq_factor(ifactor)*em_vars%omega(iomega), dir),& - em_wfs_tag(dir, ifactor), have_restart_rho=(ierr==0)) + em_wfs_tag(dir, ifactor), have_restart_rho=(ierr==0),& + vext = sys%gr%m%x(:, dir)) else call dsternheimer_solve(sh, sys, h, em_vars%lr(dir, :, ifactor), 2 , & em_vars%freq_factor(ifactor)*em_vars%omega(iomega), & - sys%gr%m%x(:, dir), RESTART_DIR,& + RESTART_DIR,& em_rho_tag(em_vars%freq_factor(ifactor)*em_vars%omega(iomega), dir),& - em_wfs_tag(dir, ifactor), have_restart_rho=(ierr==0)) + em_wfs_tag(dir, ifactor), have_restart_rho=(ierr==0), & + vext = sys%gr%m%x(:, dir)) end if