Changeset 2910

Timestamp:

05/14/07 15:48:34 (20 months ago)

Author:

xavier

Message:

Added the posibility of using non-local operators with the Sternheimer
equation, this is done by passing a set of ground state wavefunctions
with the operator applied to them instead of the local potential.
There is still the posibility of passing a local potential, to avoid
the memory cost of storing another set of wavefunctions.

Location:

trunk/src

Files:

• em_resp.F90 (modified) (1 diff)
• em_resp_calc_inc.F90 (modified) (1 diff)
• math.F90 (modified) (2 diffs)
• phonons_lr.F90 (modified) (1 diff)
• sternheimer_inc.F90 (modified) (9 diffs)
• vdW.F90 (modified) (1 diff)

trunk/src/em_resp.F90

@@ -246 +246 @@
               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
             

trunk/src/em_resp_calc_inc.F90

@@ -329 +329 @@
     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

trunk/src/math.F90

@@ -59 +59 @@
     set_app_threshold,          &
     operator(.app.),            &
+    math_xor,                   &
     ddelta
 
@@ -750 +751 @@
   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"

trunk/src/phonons_lr.F90

@@ -120 +120 @@
       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

trunk/src/sternheimer_inc.F90

@@ -26 +26 @@
 
 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 +35 @@
   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 +62 @@
 
   ASSERT( nsigma==1 .or. nsigma ==2 )
-  
+  ASSERT( math_xor(present(vext), present(vext_psi)) )
+
   m => sys%gr%m
   st => sys%st
@@ -106 +108 @@
     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 +128 @@
           !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 +252 @@
 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 +258 @@
   type(lr_t),             intent(inout) :: lr(:) 
   integer,                intent(in)    :: nsigma 
-  FLOAT,                  intent(in)    :: vext(:)
   R_TYPE,                 intent(out)   :: hvar(:,:,:)
 
@@ -268 +283 @@
   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 +290 @@
       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

trunk/src/vdW.F90

@@ -243 +243 @@
         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