Changeset 6353

Timestamp:

03/16/10 15:56:54 (6 months ago)

Author:

nitsche

Message:

Rewrote much of the transport code to simplify the modularization (as the next step).
Due to some algorithmic changes most of the complicated restart procedure is not neccessary anymore and is replaced by standard routines.

* removed the symmetric version of the interface matrix vector multiplication
* the central region is no more enlarged
* the Lippmann-Schwinger calculation now shows a progress bar
* now a magnetic groundstate is also supported
* removed some push/pop in the left hand side Lippmann-Schwinger subroutine as these are called many times during a magnetic gs calculation
* The projected part of the lead wave function is now written to a file and then read by the td run.
* the self energy is now directly calculated and stored instead of the surface Green's function
* improved debugging output for the surface Green's function
* Most of the reading part for the extended central region is now obsolete.
* The subroutine restart_read now correctly reads the weights of the kpoints.
* The INNER/OUTER logic of the interface wave function now is obsolete.
* The progress of the convergence of the Crank-Nicholson solver is shown in the debug mode.
* The zeroth source term now uses the lead wave function calculated via the self energy.
* removed an apostroph which caused a compiler warning
* The transport td run does not need a special open boundary gs calculation anymore, a normal gs would also be sufficient (although not recommended for the transport mode).
* one test of the open boundary testsuite failed because energy was the matching criteria. Since the central region is not artificially enlarged anymore the resulting energy for this region is also different.
* adjusted some tests as the values coming from the new algorithm slightly changed

Location:

trunk

Files:

• src/grid/ob_interface.F90 (modified) (2 diffs)
• src/ions/simul_box.F90 (modified) (2 diffs)
• src/scf/ground_state.F90 (modified) (1 diff)
• src/scf/ob_lippmann_schwinger.F90 (modified) (13 diffs)
• src/scf/scf.F90 (modified) (1 diff)
• src/states/ob_green.F90 (modified) (13 diffs)
• src/states/restart.F90 (modified) (8 diffs)
• src/states/states.F90 (modified) (11 diffs)
• src/td/ob_rti.F90 (modified) (4 diffs)
• src/td/ob_src.F90 (modified) (4 diffs)
• src/td/pes_rc_inc.F90 (modified) (1 diff)
• src/td/td.F90 (modified) (2 diffs)
• testsuite/open_systems/01-wavepacket.test (modified) (1 diff)
• testsuite/open_systems/03-extended_eigenstate.test (modified) (1 diff)
• testsuite/open_systems/04-propagate_eigenstate.test (modified) (1 diff)

trunk/src/grid/ob_interface.F90

@@ -38 +38 @@
   public ::                 &
     lead_t,                 &
-    interface_apply_sym_op, &
     interface_apply_op,     &
     interface_t,            &
@@ -261 +260 @@
 
   ! ---------------------------------------------------------
-  ! Apply an np x np symmetric operator op to the interface region of the wavefunction.
-  ! np is the number of points in the interface: res <- res + alpha op wf
-  subroutine interface_apply_sym_op(intf, alpha, op, wf, res)
-    type(interface_t), intent(in)    :: intf
-    CMPLX,             intent(in)    :: alpha
-    CMPLX,             intent(in)    :: op(:, :)
-    CMPLX,             intent(in)    :: wf(:)
-    CMPLX,             intent(inout) :: res(:)
-    
-    CMPLX, allocatable :: intf_wf(:), op_intf_wf(:)
-
-    call push_sub('ob_interface.interface_apply_sym_op')
-
-    SAFE_ALLOCATE(intf_wf(1:intf%np_uc))
-    SAFE_ALLOCATE(op_intf_wf(1:intf%np_uc))
-
-    call get_intf_wf(intf, wf, intf_wf)
-    call get_intf_wf(intf, res, op_intf_wf)
-    call lalg_symv(intf%np_uc, alpha, op, intf_wf, M_z1, op_intf_wf)
-    call put_intf_wf(intf, op_intf_wf, res)
-
-    SAFE_DEALLOCATE_A(intf_wf)
-    SAFE_DEALLOCATE_A(op_intf_wf)
-
-    call pop_sub()
-  end subroutine interface_apply_sym_op
-
-
-  ! ---------------------------------------------------------
   ! Write number of interface points.
   subroutine interface_write_info(intf, iunit)

trunk/src/ions/simul_box.F90

@@ -77 +77 @@
     BEFORE    = 7,              & ! for 4D open system
     AFTER     = 8,              & !
-    OUTER     = 1,              & ! Block indices of interface wavefunctions
-    INNER     = 2,              & ! for source term.
     TRANS_DIR = 1                 ! Transport is in x-direction.
 
@@ -361 +359 @@
               call write_fatal(1)
             end if
-            ! If we are doing a ground-state calculation add one more unit
-            ! cell at both ends to calculate the extended eigenstates.
-            if(calc_mode_is(CM_GS)) then
-              sb%add_unit_cells(1:NLEADS) = sb%add_unit_cells(1:NLEADS) + 1
-            end if
           case(TD_POT_FORMULA)
             call parse_block_string(blk, nr, 1, sb%lead_td_pot_formula(LEFT))

trunk/src/scf/ground_state.F90

@@ -70 +70 @@
       call states_allocate_free_states(sys%st, sys%gr)
       call read_free_states(sys%st, sys%gr)
-      ! allocate Green`s function and calculate
-      call states_init_green(sys%st, sys%gr, hm%d%nspin, hm%d%ispin, hm%lead)
+      ! allocate self_energy and calculate
+      call states_init_self_energy(sys%st, sys%gr, hm%d%nspin, hm%d%ispin, hm%lead)
     end if
 

trunk/src/scf/ob_lippmann_schwinger.F90

@@ -32 +32 @@
   use io_function_m
   use lalg_basic_m
+  use loct_m
   use math_m
   use messages_m
@@ -50 +51 @@
     lippmann_schwinger
 
-  type pg
-    CMPLX, pointer           :: green(:, :, :)
-  end type pg
+  type p_se_t
+    CMPLX, pointer           :: self_energy(:, :, :)
+  end type p_se_t
 
   ! Pointers to communicate with iterative linear solver.
   integer, pointer             :: ist_p, ik_p
   FLOAT, pointer               :: energy_p
-  type(pg), pointer            :: lead_p(:)
+  type(p_se_t), pointer        :: lead_p(:)
   type(grid_t), pointer        :: gr_p
   type(hamiltonian_t), pointer :: hm_p
@@ -78 +79 @@
     FLOAT                      :: tol, res
     CMPLX, allocatable         :: rhs(:, :)
-    type(pg), target           :: lead(2*MAX_DIM)
+    type(p_se_t), target       :: lead(2*MAX_DIM)
     logical                    :: conv
 #ifdef HAVE_MPI
@@ -94 +95 @@
         np = gr%intf(il)%np_uc
       end if
-      SAFE_ALLOCATE(lead(il)%green(1:np, 1:np, 1:st%d%dim))
+      SAFE_ALLOCATE(lead(il)%self_energy(1:np, 1:np, 1:st%d%dim))
     end do
 
@@ -108 +109 @@
     energy_p => energy
 
+    ASSERT(ubound(st%zphi, dim = 1) >= gr%mesh%np_part)
+
+    ! We have many k-points, so show the progress
+    if(mpi_grp_is_root(mpi_world)) call loct_progress_bar(-1, st%nst*st%d%nik)
+
     do ik = st%d%kpt%start, st%d%kpt%end
       do ist = st%st_start, st%st_end
-        ! Solve Schroedinger equation for this energy.
+        ! Solve Lippmann-Schwinger equation for this energy.
         energy = st%ob_eigenval(ist, ik)
         st%eigenval(ist, ik) = energy
-
-        ASSERT(ubound(st%zphi, dim = 1) >= gr%mesh%np_part)
+        do il = 1, NLEADS
+          lead(il)%self_energy(:, :, :) = st%ob_lead(il)%self_energy(:, :, :, ist, ik)
+        end do
 
         ! Calculate right hand side e-T-V0-sum(a)[H_ca*g_a*H_ac].
-        rhs(:, :) = M_z0
-        call zhamiltonian_apply(hm, gr%der, st%zphi(:, :, ist, ik), rhs(:, :), ist, ik, kinetic_only=.true.)
-
-        ! Apply lead potential.
-        ! FIXME: this wrong when non-symmetric leads are used, here the periodicity is assumed
-        ! for the Lippmann-Schwinger solution
-        do idim = 1, st%d%dim
-          do ip = 1, gr%mesh%np
-            ip_lead = mod(ip-1, gr%intf(LEFT)%np_intf) + 1
-            rhs(ip, idim) = rhs(ip, idim) + hm%lead(LEFT)%vks(ip_lead, idim)*st%zphi(ip, idim, ist, ik)
-          end do
-        end do
-
-        ! Add energy.
-        do idim = 1, st%d%dim
-          rhs(1:gr%mesh%np, idim) = energy*st%zphi(1:gr%mesh%np, idim, ist, ik) - rhs(1:gr%mesh%np, idim)
-        end do
-
-        ! Apply term with lead Green functions.
-        do il = 1, NLEADS
-          do idim = 1, st%d%dim
-            np_intf = gr%intf(il)%np_intf
-            call apply_coupling(st%ob_lead(il)%green(1:np_intf, 1:np_intf, idim, ist, ik), &
-                            hm%lead(il)%h_offdiag(:, :), lead(il)%green(:, :, idim), np_intf, il)
-          end do
-        end do
-        do il = 1, NLEADS
-          do idim = 1, st%d%dim
-            if (associated(hm%ep%A_static)) then ! magnetic gs
-              call interface_apply_op(gr%intf(il), -M_z1, lead(il)%green(:, :, idim), &
-                st%zphi(:, idim, ist, ik), rhs(:, idim))
-            else
-              call interface_apply_sym_op(gr%intf(il), -M_z1, lead(il)%green(:, :, idim), &
-                st%zphi(:, idim, ist, ik), rhs(:, idim))
-            end if
-          end do
-        end do
-
-        if (associated(hm%ep%A_static)) then ! magnetic gs
-        ! FIXME: multiply rhs with (e-h-g)^T
-        else
-        end if
+        call calc_rhs(st%zphi(:, :, ist, ik), rhs)
+
+        if (associated(hm%ep%A_static)) call calc_rhs(rhs, rhs, .true.) ! multiply transposed version
+
         ! Solve linear system lhs psi = rhs.
         iter = eigens%es_maxiter
@@ -164 +133 @@
 
         conv = .false.
-        call zqmr_sym(gr%mesh%np_part*st%d%dim, st%zpsi(:, 1, ist, ik), rhs(:, 1), lhs, dotu, nrm2, precond, &
-          iter, residue=res, threshold=tol, converged=conv)
-        ! call zqmr(gr%mesh%np_part*st%d%dim, st%zpsi(:, 1, ist, ik), rhs(:, 1), lhs, lhs_t, &
-        !   precond, precond, iter, residue=dres, threshold=tol, converged=conv, showprogress=.true.)
+        if (associated(hm%ep%A_static)) then ! magnetic gs
+          call zqmr_sym(gr%mesh%np_part*st%d%dim, st%zpsi(:, 1, ist, ik), rhs(:, 1), lhs_symmetrized, dotu, &
+            nrm2, precond, iter, residue=res, threshold=tol, converged=conv, showprogress=in_debug_mode)
+        else
+          call zqmr_sym(gr%mesh%np_part*st%d%dim, st%zpsi(:, 1, ist, ik), rhs(:, 1), lhs, dotu, nrm2, precond, &
+            iter, residue=res, threshold=tol, converged=conv, showprogress=in_debug_mode)
+        end if
+        if(in_debug_mode) then ! write info
+          write(message(1), '(a,i8,e10.3)') 'Iterations, Residual: ', iter, res
+          call write_info(1)
+        end if
 
         eigens%matvec = eigens%matvec + iter + 1 + 2
-        if(conv) then
-          eigens%converged = eigens%converged + 1
-        end if
+        if(conv) eigens%converged = eigens%converged + 1
         eigens%diff(ist, ik) = res
+        if(mpi_grp_is_root(mpi_world)) call loct_progress_bar(st%nst*(ik - 1), st%nst*st%d%nik)
       end do
     end do
@@ -199 +174 @@
     SAFE_DEALLOCATE_A(rhs)
     do il = 1, NLEADS
-      SAFE_DEALLOCATE_P(lead(il)%green)
+      SAFE_DEALLOCATE_P(lead(il)%self_energy)
     end do
 
     call pop_sub()
   end subroutine lippmann_schwinger
+
+
+  ! ---------------------------------------------------------
+  ! The right hand side of the Lippmann-Schwinger equation
+  ! e-T-V0-sum(a)[H_ca*g_a*H_ac].
+  subroutine calc_rhs(zphi, rhs, transposed)
+    CMPLX, intent(in)             :: zphi(:, :)
+    CMPLX, intent(out)            :: rhs(:, :)
+    logical, optional, intent(in) :: transposed ! needed only for the non hermitian part
+
+    integer :: ip, idim, il
+    logical :: transposed_
+    CMPLX, allocatable :: tmp(:, :)
+
+    call push_sub('ob_lippmann_schwinger.calc_rhs')
+
+    SAFE_ALLOCATE(tmp(1:gr_p%mesh%np_part, 1:st_p%d%dim))
+
+    transposed_ = .false.
+    if(present(transposed)) transposed_ = transposed
+
+    if(transposed_) then ! the usual conjugate trick for the hermitian part
+      tmp = conjg(zphi)
+    else
+      tmp = zphi
+    end if
+    ! Calculate right hand side e-T-V0-sum(a)[H_ca*g_a*H_ac].
+    rhs(:, :) = M_z0
+
+    call zhamiltonian_apply(hm_p, gr_p%der, tmp, rhs, ist_p, ik_p, kinetic_only=.true.)
+
+    ! Apply lead potential. Left and right lead potential are assumed to be equal.
+    forall(ip = 1:gr_p%mesh%np )
+      rhs(ip, :) = rhs(ip, :) + hm_p%lead(LEFT)%vks(mod(ip-1, gr_p%intf(LEFT)%np_intf) + 1, :) * tmp(ip, :)
+    end forall
+
+    ! Add energy.
+    forall(ip = 1:gr_p%mesh%np ) rhs(ip, :) = energy_p * tmp(ip, :) - rhs(ip, :)
+
+    if(transposed_) then
+      rhs = conjg(rhs)
+      tmp = conjg(tmp) ! restore original
+    end if
+
+    do il = 1, NLEADS
+      do idim = 1, st_p%d%dim
+        if(transposed_) then
+          call interface_apply_op(gr_p%intf(il), -M_z1, transpose(lead_p(il)%self_energy(:, :, idim)), &
+                                tmp(:, idim), rhs(:, idim))
+        else
+          call interface_apply_op(gr_p%intf(il), -M_z1, lead_p(il)%self_energy(:, :, idim), &
+                                tmp(:, idim), rhs(:, idim))
+        end if
+      end do
+    end do
+    
+    SAFE_DEALLOCATE_A(tmp)
+    call pop_sub()
+  end subroutine calc_rhs
 
 
@@ -212 +246 @@
     CMPLX, intent(in) :: y(:)
 
-    integer :: np_part, np, dim, idim
+    integer :: np_part, np, idim
     CMPLX   :: dot
 
-    call push_sub('ob_lippmann_schwinger.dotu')
+!    call push_sub('ob_lippmann_schwinger.dotu')
 
     np_part = gr_p%mesh%np_part
     np      = gr_p%mesh%np
-    dim     = st_p%d%dim
 
     dot = M_ZERO
-    do idim = 1, dim
+    do idim = 1, st_p%d%dim
       dot = dot + blas_dotu(np, x(l(idim)), 1, y(l(idim)), 1)
     end do
     dotu = dot
 
-    call pop_sub()
+!    call pop_sub()
 
   contains
@@ -245 +278 @@
     CMPLX, intent(in) :: x(:)
 
-    integer :: np_part, np, dim, idim
+    integer :: np_part, np, idim
     FLOAT   :: nrm
 
-    call push_sub('ob_lippmann_schwinger.nrm2')
+!    call push_sub('ob_lippmann_schwinger.nrm2')
 
     np_part = gr_p%mesh%np_part
     np      = gr_p%mesh%np
-    dim     = st_p%d%dim
 
     nrm = M_ZERO
-    do idim = 1, dim
+    do idim = 1, st_p%d%dim
       nrm = nrm + lalg_nrm2(np, x(l(idim):u(idim)))
     end do
     nrm2 = nrm
 
-    call pop_sub()
+!    call pop_sub()
 
   contains
@@ -290 +322 @@
     integer            :: np_part, idim, il, dim
 
-    call push_sub('ob_lippmann_schwinger.lhs')
+!    call push_sub('ob_lippmann_schwinger.lhs')
 
     np_part = gr_p%mesh%np_part
@@ -305 +337 @@
     ! y <- e x - tmp_y
     do idim = 1, dim
-      y(l(idim):u(idim)) = energy_p*x(l(idim):u(idim)) - tmp_y(1:np_part, idim)
+      tmp_y(:, idim) = energy_p * x(l(idim):u(idim)) - tmp_y(1:np_part, idim)
     end do
 
     do il = 1, NLEADS
       do idim = 1, dim
-        call interface_apply_sym_op(gr_p%intf(il), &
-          -M_z1, lead_p(il)%green(:, :, idim), tmp_x(:, idim), y(l(idim):u(idim)))
+        call interface_apply_op(gr_p%intf(il), -M_z1, &
+          lead_p(il)%self_energy(:, :, idim), tmp_x(:, idim), tmp_y(:, idim))
       end do
+    end do
+
+    do idim = 1, dim
+      y(l(idim):u(idim)) = tmp_y(1:np_part, idim)
     end do
 
     SAFE_DEALLOCATE_A(tmp_x)
     SAFE_DEALLOCATE_A(tmp_y)
-    call pop_sub()
+!    call pop_sub()
 
   contains
@@ -333 +369 @@
     end function u
   end subroutine lhs
+
+
+  ! ---------------------------------------------------------
+  ! The left hand side of the Lippmann-Schwinger equation
+  ! (e-H-sum(a)[H_ca*g_a*H_ac])^T.
+  ! Used by the iterative linear solver.
+  subroutine lhs_t(x, y)
+    CMPLX, intent(in)  :: x(:)
+    CMPLX, intent(out) :: y(:)
+
+    CMPLX, allocatable :: tmp_x(:, :)
+    CMPLX, allocatable :: tmp_y(:, :)
+    integer            :: np_part, idim, il, dim
+
+!    call push_sub('ob_lippmann_schwinger.lhs_t')
+
+    np_part = gr_p%mesh%np_part
+    dim     = st_p%d%dim
+
+    SAFE_ALLOCATE(tmp_x(1:np_part, 1:dim))
+    SAFE_ALLOCATE(tmp_y(1:np_part, 1:dim))
+
+    do idim = 1, dim
+      tmp_x(1:np_part, idim) = conjg(x(l(idim):u(idim)))
+    end do
+    call zhamiltonian_apply(hm_p, gr_p%der, tmp_x, tmp_y, ist_p, ik_p)
+
+    ! y <- e x - tmp_y
+    do idim = 1, dim
+      tmp_y(1:np_part, idim) = energy_p * tmp_x(1:np_part, idim) - tmp_y(1:np_part, idim)
+    end do
+    tmp_y = conjg(tmp_y)
+    tmp_x = conjg(tmp_x) ! restore for the non-hermitian part
+
+    do il = 1, NLEADS
+      do idim = 1, dim
+        call interface_apply_op(gr_p%intf(il), -M_z1, transpose(lead_p(il)%self_energy(:, :, idim)), &
+                                tmp_x(:, idim), tmp_y(:, idim))
+      end do
+    end do
+
+    do idim = 1, dim
+      y(l(idim):u(idim)) = tmp_y(1:np_part, idim)
+    end do
+
+    SAFE_DEALLOCATE_A(tmp_x)
+    SAFE_DEALLOCATE_A(tmp_y)
+!    call pop_sub()
+
+  contains
+
+    integer function l(idim)
+      integer, intent(in) :: idim
+
+      l = (idim-1)*np_part+1
+    end function l
+
+    integer function u(idim)
+      integer, intent(in) :: idim
+
+      u = l(idim)+np_part-1
+    end function u
+  end subroutine lhs_t
+
+
+  ! ---------------------------------------------------------
+  ! The left hand side of the Lippmann-Schwinger equation
+  ! (e-H-sum(a)[H_ca*g_a*H_ac])^T*(e-H-sum(a)[H_ca*g_a*H_ac]).
+  ! Used by the iterative linear solver.
+  subroutine lhs_symmetrized(x, y)
+    CMPLX, intent(in)  :: x(:)
+    CMPLX, intent(out) :: y(:)
+
+!    call push_sub('ob_lippmann_schwinger.lhs_symmetrized')
+
+    call lhs(x, y)
+    call lhs_t(y, y)
+
+!    call pop_sub()
+  end subroutine lhs_symmetrized
 
 
@@ -343 +459 @@
     CMPLX, intent(out) :: y(:)
 
-    call push_sub('ob_lippmann_schwinger.precond')
+!    call push_sub('ob_lippmann_schwinger.precond')
 
     y(:) = x(:)
 
-    call pop_sub()
+!    call pop_sub()
   end subroutine precond
 end module ob_lippmann_schwinger_m

trunk/src/scf/scf.F90

@@ -576 +576 @@
       call scf_write_static(STATIC_DIR, "info")
       call h_sys_output_all(outp, gr, geo, st, hm, STATIC_DIR)
+
+      if(gr%sb%open_boundaries) then ! write part of the source term s(0)
+        call states_write_proj_lead_wf('open_boundaries/', gr%intf, st)
+      end if
+
     end if
 

trunk/src/states/ob_green.F90

@@ -24 +24 @@
 module ob_green_m
   use global_m
-  use grid_m
   use lalg_adv_m
   use lalg_basic_m
-  use parser_m
   use math_m
   use messages_m
-  use nl_operator_m
   use ob_interface_m
   use profiling_m
   use simul_box_m
-  use string_m
 
   implicit none
@@ -40 +36 @@
 
   public ::         &
-    lead_green
+    lead_self_energy
 
 contains
 
-  ! calculate the surface Green`s function
+  ! calculate the lead self energy
   ! prefer the fast method if possible, but if not converged check also other method
-  subroutine lead_green(energy, diag, offdiag, intf, green, h_is_real)
+  subroutine lead_self_energy(energy, diag, offdiag, intf, self_energy, h_is_real)
     FLOAT,             intent(in)  :: energy        ! Energy to calculate Green`s function for.
     CMPLX,             intent(in)  :: diag(:, :)    ! Diagonal block of lead Hamiltonian.
     CMPLX,             intent(in)  :: offdiag(:, :) ! Off-diagonal block of lead Hamiltonian.
     type(interface_t), intent(in)  :: intf          ! => gr%intf(il)
-    CMPLX,             intent(out) :: green(:, :)   ! The calculated Green`s function.
+    CMPLX,             intent(inout) :: self_energy(:, :) ! The calculated self_energy.
     logical,           intent(in)  :: h_is_real     ! Is the Hamiltonian real? (no vector potential)
 
@@ -59 +55 @@
     CMPLX, allocatable :: green2(:, :)
 
-    call push_sub('ob_green.lead_green')
+    call push_sub('ob_green.lead_self_energy')
 
     np = intf%np_intf
@@ -72 +68 @@
     ! we have only one way of calculating the Green`s function: with the Sancho method
     if(.not.intf%reducible) then ! use large matrices (rank = np_uc)
-      call lead_green_sancho(c_energy, diag, offdiag, np_uc, green, threshold, h_is_real)
-      residual = calc_residual_green(energy, green, diag, offdiag, intf)
+      call lead_green_sancho(c_energy, diag, offdiag, np_uc, self_energy, threshold, h_is_real)
+      residual = calc_residual_green(energy, self_energy, diag, offdiag, intf)
       if(in_debug_mode) then ! write info
         write(message(1), '(a)') 'Non-reducible unit cell'
@@ -81 +77 @@
     else
       ! 1. calculate with the fastest method
-      call lead_green_umerski(c_energy, diag, offdiag, intf, green)
+      call lead_green_umerski(c_energy, diag, offdiag, intf, self_energy)
       ! 2. check if converged
-      residual = calc_residual_green(energy, green, diag, offdiag, intf)
+      residual = calc_residual_green(energy, self_energy, diag, offdiag, intf)
       if(in_debug_mode) then ! write info
         write(message(1), '(a,e10.3)') 'Umerski-Residual = ', residual
@@ -99 +95 @@
         ! 4. if also not converged choose the best one and give warning
         if(residual2.lt.eps) then ! finally converged
-          green = green2
+          self_energy = green2
         else ! write warning
           if(residual2.lt.residual) then
-            green = green2
+            self_energy = green2
           end if
           message(1) = 'The surface Green`s function did not converge properly'
@@ -115 +111 @@
       end if
     end if
-
-    call pop_sub()
-  end subroutine lead_green
+    ! now calculate the self energy
+    if(intf%il.eq.LEFT) then
+      call lalg_trmm(np, np, 'U', 'N', 'L', M_z1, offdiag, self_energy)
+      call lalg_trmm(np, np, 'U', 'T', 'R', M_z1, offdiag, self_energy)
+    else
+      call lalg_trmm(np, np, 'L', 'N', 'L', M_z1, offdiag, self_energy)
+      call lalg_trmm(np, np, 'L', 'T', 'R', M_z1, offdiag, self_energy)
+    end if
+
+    call pop_sub()
+  end subroutine lead_self_energy
 
 
   ! check if the Green`s function gives the correct density of states
   ! if not compute its Hermitian conjugate
-  subroutine fix_green(np, green, dos)
+  subroutine fix_green(np, green)
     integer, intent(in)    :: np
     CMPLX,   intent(inout) :: green(:, :)
-    FLOAT,   intent(out)   :: dos
 
     integer  :: j
+    FLOAT    :: dos
 
     call push_sub('ob_green.fix_green')
@@ -136 +140 @@
       dos = dos - aimag(green(j, j))
     end do
+
+    if(in_debug_mode) then ! write info
+      write(message(1), '(a,e10.3)') 'density of states = ', abs(dos)
+      call write_info(1)
+    end if
+
     if(dos.lt.M_ZERO) then
       green(:, :) = transpose(conjg(green(:, :)))
+      if(in_debug_mode) then ! write info
+        message(1) = "surface Green's function changed to its hermitian conjugate"
+        call write_info(1)
+      end if
     end if
 
@@ -207 +221 @@
     CMPLX, allocatable :: tmp1(:, :), tmp2(:, :), tmp3(:, :)
     integer            :: i, j
-    FLOAT              :: dos, old_norm, norm, res
+    FLOAT              :: det, old_norm, norm, res
 
     call push_sub('ob_green.lead_green_sancho')
@@ -236 +250 @@
 
       forall (j = 1:np) inv(j, j) = inv(j, j) + energy
-      dos = lalg_inverter(np, inv, invert=.true.)
+      det = lalg_inverter(np, inv, invert=.true.)
 
       call lalg_gemm(np, np, np, M_z1, a, inv, M_z0, tmp2)
@@ -262 +276 @@
 
     forall (j = 1:np) green(j, j) = green(j, j) + energy
-    dos = lalg_inverter(np, green, invert = .true.)
+    det = lalg_inverter(np, green, invert = .true.)
     if (h_is_real) then ! the Green`s function is complex symmetric
       call matrix_symmetric_average(green, np)
     end if ! otherwise it is general complex
 
-    call fix_green(np, green, dos)
-
     if(in_debug_mode) then ! write some info
-      write(message(1), '(a,e10.3)') 'Sancho-DOS = ', dos
+      message(1) = "surface Green's function: iterative algorithm"
       call write_info(1)
     end if
+
+    call fix_green(np, green)
 
     SAFE_DEALLOCATE_A(e)
@@ -393 +407 @@
     integer              :: ib, np, np2
     CMPLX, allocatable   :: x(:,:), x2(:,:), o2(:,:), o4(:,:), d(:), h(:, :, :), v(:, :, :)
-    FLOAT                :: dos
+    FLOAT                :: det
 
     call push_sub('ob_green.lead_green_umerski')
@@ -435 +449 @@
     o4(1:np, 1:np) = x(np+1:np2, np+1:np2)
     ! 4. calculate green = o2*o4^(-1)
-    dos = lalg_inverter(np, o4, invert = .true.)
+    det = lalg_inverter(np, o4, invert = .true.)
     call zgemm('N', 'N', np, np, np, M_z1, o2, np, o4, np, M_z0, green, np)
 
-    call fix_green(np, green, dos)
     if(in_debug_mode) then ! write some info
-      write(message(1), '(a,e10.5)') 'Umerski: DOS = ', dos
+      message(1) = "surface Green's function: direct algorithm (Moebius transformation)"
       call write_info(1)
     end if
+
+    call fix_green(np, green)
 
     SAFE_DEALLOCATE_A(x)

trunk/src/states/restart.F90

@@ -28 +28 @@
   use io_m
   use io_function_m
+  use lalg_basic_m
   use linear_response_m
   use loct_m
@@ -36 +37 @@
   use messages_m
   use mpi_m
+  use ob_green_m
   use parser_m
   use profiling_m
@@ -56 +58 @@
     restart_write,           &
     restart_read,            &
-    restart_read_ob_intf,    &
-    restart_read_ob_central, &
+    restart_get_ob_intf,     &
     restart_format,          &
     restart_dir,             &
@@ -365 +366 @@
 
   ! ---------------------------------------------------------
-  ! Reads the interface regions of the ground-state wavefunctions of
-  ! an open-boundaries calculation.
-  ! Can also read the inner interface wf of a non-open boundary gs calculation
-  ! Returns (in ierr)
-  ! <0 => Fatal error,
-  ! =0 => read all wavefunctions,
-  ! >0 => could only read x wavefunctions.
-  subroutine restart_read_ob_intf(dir, st, gr, ierr)
-    character(len=*), intent(in)    :: dir
+  ! Reads the interface regions of the wavefunctions
+  subroutine restart_get_ob_intf(st, gr)
     type(states_t),   intent(inout) :: st
-    type(grid_t),     intent(inout) :: gr
-    integer,          intent(out)   :: ierr
-
-    integer              :: io_wfns, io_occs, io_mesh, np, lead_np, np_uc, int, err, il, ip
-    integer              :: ik, ist, idim, tnp
-    FLOAT                :: flt
-    character            :: char
-    character(len=256)   :: line, filename
-    CMPLX, allocatable   :: tmp(:)
-    type(mesh_t)         :: gs_mesh
-    type(mpi_grp_t)      :: mpi_grp
-
-    call push_sub('restart.restart_read_ob_intf')
+    type(grid_t),     intent(in) :: gr
+
+    integer              :: ik, ist, idim, il, ip
+
+    call push_sub('restart.restart_get_ob_intf')
 
     write(message(1), '(a,i5)') 'Info: Reading ground-state interface wavefunctions.'
     call write_info(1)
 
-    mpi_grp = mpi_world
     ! Sanity check.
     do il = 1, NLEADS
@@ -399 +384 @@
     end do
 
-    ierr = 0
-
-    ! Read the mesh of the ground-state calculation.
-    io_mesh = io_open(trim(dir)//'/mesh', action='read', status='old', die=.false., is_tmp=.true.)
-    if(io_mesh.lt.0) then
-      ierr = -1
-      call io_close(io_mesh, grp=mpi_grp)
-      call pop_sub()
-      return
-    end if
-      
-    gs_mesh%sb => gr%mesh%sb
-    call mesh_init_from_file(gs_mesh, io_mesh)
-    call io_close(io_mesh) 
-
-    io_wfns  = io_open(trim(dir)//'/wfns', action='read', status='old', die=.false., is_tmp=.true., grp=mpi_grp)
-    if(io_wfns.lt.0) then
-      ierr = -1
-    end if
-
-    io_occs = io_open(trim(dir)//'/occs', action='read', status='old', die=.false., is_tmp=.true., grp=mpi_grp)
-    if(io_occs.lt.0) then
-      if(io_wfns.gt.0) call io_close(io_wfns, grp=mpi_grp)
-      ierr = -1
-    end if
-
-    if(ierr .ne. 0) then
-      write(message(1),'(a)') 'Could not load ground-state interface wavefunctions.'
-      call write_info(1)
-      call messages_print_stress(stdout)
-      call pop_sub()
-      return
-    end if
-
-    ! Skip two lines.
-    call iopar_read(mpi_grp, io_wfns, line, err)
-    call iopar_read(mpi_grp, io_wfns, line, err)
-
-    call iopar_read(mpi_grp, io_occs, line, err)
-    call iopar_read(mpi_grp, io_occs, line, err)
-
-    ! FIXME: make the assertion more exact!
-!    lead_np = gr%mesh%lead_unit_cell(LEFT)%np
-!    np = gr%intf(LEFT)%np_intf
-    ASSERT(gr%mesh%np.le.gs_mesh%np)
-    SAFE_ALLOCATE(tmp(1:gs_mesh%np))
-    
-    do
-      call iopar_read(mpi_grp, io_wfns, line, int)
-      read(line, '(a)') char
-      if(int .ne. 0 .or. char .eq. '%') exit
-
-      call iopar_backspace(mpi_grp, io_wfns)
-
-      call iopar_read(mpi_grp, io_wfns, line, err)
-      read(line, *) ik, char, ist, char, idim, char, filename
-
-      call iopar_read(mpi_grp, io_occs, line, err)
-      read(line, *) st%occ(ist, ik), char, st%eigenval(ist, ik), char, flt, char, flt,&
-                    char, flt, char, st%d%kweights(ik)
-
-      if(ist .ge. st%st_start .and. ist .le. st%st_end .and. &
-         ik .ge. st%d%kpt%start .and. ik .le. st%d%kpt%end) then
-        ! Open boundaries imply complex wavefunctions.
-        call zrestart_read_function(dir, filename, gs_mesh, tmp, err)
-
-        
-        if(err.le.0) then
-          ierr = ierr + 1 ! count the valid file readings
-          if(gr%mesh%np .eq. gs_mesh%np) then ! no extra unit cell present
-            do il = 1, NLEADS
-              forall(ip = 1:gr%intf(il)%np_uc)
-                st%ob_lead(il)%intf_psi(ip, INNER, idim, ist, ik) = tmp(gr%intf(il)%index(ip))
-              end forall
-              ! the outer part is zero (if the source term is still active)
-              st%ob_lead(il)%intf_psi(:, OUTER, idim, ist, ik) = M_ZERO
-            end do
-          else ! transport mode with extra unit cell
-            ! FIXME: to be generalized for more than 2 leads
-            lead_np = gr%mesh%lead_unit_cell(LEFT)%np
-            np = gr%intf(LEFT)%np_intf
-            np_uc = gr%intf(LEFT)%np_uc
-            ! Here we use the fact that transport is in x-direction (x is the index
-            ! running slowest).
-            ! Outer block.
-            tnp = gr%mesh%np + 2 * lead_np
-            ! FIXME: check for extent>der-order if this is still correct (left side)
-            st%ob_lead(LEFT)%intf_psi(1:np, OUTER, idim, ist, ik)  = tmp(lead_np - np + 1:lead_np)
-            st%ob_lead(RIGHT)%intf_psi(1:np, OUTER, idim, ist, ik) = tmp(tnp - lead_np + 1:tnp - lead_np + np)
-            ! Inner block.
-            st%ob_lead(LEFT)%intf_psi(1:np, INNER, idim, ist, ik)  = tmp(lead_np + 1:lead_np+np)
-            st%ob_lead(RIGHT)%intf_psi(1:np, INNER, idim, ist, ik) = tmp(tnp - lead_np - np + 1:tnp - lead_np)
-          end if
-        end if ! err
-      end if
+    do ik = st%d%kpt%start, st%d%kpt%end
+      do ist = st%st_start, st%st_end
+        do idim = 1, st%d%dim
+          do il = 1, NLEADS
+            forall(ip = 1:gr%intf(il)%np_uc)
+              st%ob_lead(il)%intf_psi(ip, idim, ist, ik) = st%zpsi(gr%intf(il)%index(ip), idim, ist, ik)
+            end forall
+          end do
+        end do
+      end do
     end do
-    
-#if defined(HAVE_MPI)
-    if(st%parallel_in_states) then
-      call MPI_Allreduce(ierr, err, 1, MPI_INTEGER, MPI_SUM, st%mpi_grp%comm, mpi_err)
-      ierr = err
-    end if
-    if(st%d%kpt%parallel) then
-      call MPI_Allreduce(ierr, err, 1, MPI_INTEGER, MPI_SUM, st%d%kpt%mpi_grp%comm, mpi_err)
-      ierr = err
-    end if
-#endif
-
-    if(ierr.eq.0) then
-      ierr = -1
-      call messages_print_stress(stdout, 'Loading restart information')
-      message(1) = 'No interface wavefunctions could be read.'
-      call write_info(1)
-      call messages_print_stress(stdout)
-    else
-      if(ierr .eq. st%nst * st%d%nik * st%d%dim) then
-        ierr = 0
-      else
-        call messages_print_stress(stdout, 'Loading restart information')
-        write(message(1),'(a,i4,a,i4,a)') 'Only ', ierr,' interface wavefunctions out of ', &
-          st%nst * st%d%nik * st%d%dim, ' could be read.'
-        call write_info(1)
-        call messages_print_stress(stdout)
-      end if
-    end if
-
-    call io_close(io_wfns, mpi_grp)
-    call io_close(io_occs, mpi_grp)
-
-    SAFE_DEALLOCATE_A(tmp)
 
     call pop_sub()
-  end subroutine restart_read_ob_intf
-
-
-  ! ---------------------------------------------------------
-  ! Reads the central regions of the ground-state wavefunctions of
-  ! an open-boundaries calculation.
-  ! Returns (in ierr)
-  ! <0 => Fatal error,
-  ! =0 => read all wavefunctions,
-  ! >0 => could only read x wavefunctions.
-  subroutine restart_read_ob_central(dir, st, gr, ierr)
-    character(len=*), intent(in)    :: dir
-    type(states_t),   intent(inout) :: st
-    type(grid_t),     intent(in)    :: gr
-    integer,          intent(out)   :: ierr
-
-    logical              :: psi_allocated
-    integer              :: io_wfns, io_occs, io_mesh, lead_np, int, err
-    integer              :: ik, ist, idim
-    character            :: char
-    character(len=256)   :: line, filename
-    CMPLX, allocatable   :: tmp(:)
-    type(mesh_t)         :: gs_mesh
-    type(mpi_grp_t)      :: mpi_grp
-    FLOAT                :: k_x, k_y, k_z
-
-    call push_sub('restart.restart_read_ob_central')
-
-    write(message(1), '(a,i5)') 'Info: Loading restart information'
-    call write_info(1)
-
-    mpi_grp = mpi_world
-
-    ! Sanity check.
-    psi_allocated = (associated(st%dpsi) .and. st%wfs_type.eq.M_REAL) .or. &
-      (associated(st%zpsi) .and. st%wfs_type.eq.M_CMPLX)
-    ASSERT(psi_allocated)
-
-    ierr = 0
-
-    ! Read the mesh of the ground state calculation.
-    io_mesh = io_open(trim(dir)//'/mesh', action='read', status='old', die=.false., is_tmp=.true.)
-    if(io_mesh.lt.0) then
-      ierr = -1
-      call io_close(io_mesh, grp=mpi_grp)
-      call pop_sub(); return
-    end if
-    gs_mesh%sb => gr%mesh%sb
-    call mesh_init_from_file(gs_mesh, io_mesh)
-    call io_close(io_mesh)
-
-    io_wfns  = io_open(trim(dir)//'/wfns', action='read', status='old', die=.false., is_tmp=.true., grp=mpi_grp)
-    if(io_wfns.lt.0) then
-      ierr = -1
-    end if
-
-    io_occs = io_open(trim(dir)//'/occs', action='read', status='old', die=.false., is_tmp = .true., grp = mpi_grp)
-    if(io_occs.lt.0) then
-      if(io_wfns.gt.0) call io_close(io_wfns, grp=mpi_grp)
-      ierr = -1
-    end if
-
-    if(ierr .ne. 0) then
-      write(message(1),'(a)') 'Could not load any previous restart information.'
-      call write_info(1)
-      call messages_print_stress(stdout)
-      call pop_sub()
-      return
-    end if
-
-    ! Skip two lines.
-    call iopar_read(mpi_grp, io_wfns, line, err)
-    call iopar_read(mpi_grp, io_wfns, line, err)
-
-    call iopar_read(mpi_grp, io_occs, line, err)
-    call iopar_read(mpi_grp, io_occs, line, err)
-
-    lead_np = gr%mesh%lead_unit_cell(LEFT)%np
-    SAFE_ALLOCATE(tmp(1:gr%mesh%np + 2 * lead_np))
-    
-    do
-      call iopar_read(mpi_grp, io_wfns, line, int)
-      if(int .ne. 0) exit
-      read(line, '(a)') char
-      if(char.eq.'%') exit
-
-      call iopar_backspace(mpi_grp, io_wfns)
-
-      call iopar_read(mpi_grp, io_wfns, line, err)
-      read(line, *) ik, char, ist, char, idim, char, filename
-
-      call iopar_read(mpi_grp, io_occs, line, err)
-      read(line, *) st%occ(ist, ik), char, st%eigenval(ist, ik), char, &
-                    k_x, char, k_y, char, k_z, char, st%d%kweights(ik)
-      st%d%kpoints(:, ik) = (/k_x, k_y, k_z/)
-
-      if(ist .ge. st%st_start .and. ist .le. st%st_end .and. &
-         ik .ge. st%d%kpt%start .and. ik .le. st%d%kpt%end) then
-        ! Open boundaries imply complex wavefunctions.
-        call zrestart_read_function(dir, filename, gs_mesh, tmp, err)
-        ! Here we use the fact that transport is in x-direction (x is the index
-        ! running slowest).
-        st%zpsi(1:gr%mesh%np, idim, ist, ik) = tmp(lead_np + 1:gr%mesh%np + lead_np)
-        if(err.le.0) then
-          ierr = ierr + 1
-        end if
-      end if
-    end do
-    
-#if defined(HAVE_MPI)
-    if(st%parallel_in_states) then
-      call MPI_Allreduce(ierr, err, 1, MPI_INTEGER, MPI_SUM, st%mpi_grp%comm, mpi_err)
-      ierr = err
-    end if
-    if(st%d%kpt%parallel) then
-      call MPI_Allreduce(ierr, err, 1, MPI_INTEGER, MPI_SUM, st%d%kpt%mpi_grp%comm, mpi_err)
-      ierr = err
-    end if
-#endif
-    if(ierr.eq.0) then
-      ierr = -1
-      call messages_print_stress(stdout, 'Loading restart information')
-      message(1) = 'No files could be read. No restart information can be used.'
-      call write_info(1)
-      call messages_print_stress(stdout)
-    else
-      if(ierr.eq.st%nst * st%d%nik * st%d%dim) then
-        ierr = 0
-      else
-        call messages_print_stress(stdout, 'Loading restart information')
-        write(message(1),'(a,i4,a,i4,a)') 'Only ', ierr,' files out of ', &
-          st%nst * st%d%nik * st%d%dim, ' could be read.'
-        call write_info(1)
-        call messages_print_stress(stdout)
-      end if
-    end if
-
-    call io_close(io_wfns, mpi_grp)
-    call io_close(io_occs, mpi_grp)
-
-    SAFE_DEALLOCATE_A(tmp)
-
-    call pop_sub()
-  end subroutine restart_read_ob_central
+  end subroutine restart_get_ob_intf
 
 
@@ -700 +423 @@
     character(len=50)    :: str
 
-    FLOAT                :: my_occ
+    FLOAT                :: my_occ, flt
     logical              :: read_occ_, gs_allocated, lr_allocated
 
@@ -719 +442 @@
     read_occ_ = .true.
     if(present(read_occ)) read_occ_ = .false.
-
+    
     ! sanity check
     gs_allocated = (associated(st%dpsi) .and. st%wfs_type == M_REAL) .or. &
@@ -780 +503 @@
 
       call iopar_read(gr%mesh%mpi_grp, iunit2, line, err)
-      read(line, *) my_occ, char, st%eigenval(ist, ik)
+      read(line, *) my_occ, char, st%eigenval(ist, ik), char, flt, char, flt,&
+                    char, flt, char, st%d%kweights(ik)
       if(read_occ_) st%occ(ist, ik) = my_occ
 

trunk/src/states/states.F90

@@ -95 +95 @@
     states_freeze_orbitals,           &
     states_total_density,             &
-    states_init_green
+    states_init_self_energy,          &
+    states_write_proj_lead_wf,        &
+    states_read_proj_lead_wf
 
   public ::                           &
@@ -102 +104 @@
 
   type states_lead_t
-    CMPLX, pointer     :: intf_psi(:, :, :, :, :) !< (np, 2, st%d%dim, st%nst, st%d%nik)
+    CMPLX, pointer     :: intf_psi(:, :, :, :) !< (np, st%d%dim, st%nst, st%d%nik)
     FLOAT, pointer     :: rho(:, :)   !< Density of the lead unit cells.
-    CMPLX, pointer     :: green(:, :, :, :, :) !< (np, np, nspin, ncs, nik) Green`s function of the leads.
+    CMPLX, pointer     :: self_energy(:, :, :, :, :) !< (np, np, nspin, ncs, nik) self energy of the leads.
   end type states_lead_t
 
@@ -186 +188 @@
     nullify(st%dpsi, st%zpsi, st%zphi, st%rho, st%current, st%rho_core, st%frozen_rho, st%eigenval)
     do il=1, NLEADS
-      nullify(st%ob_lead(il)%intf_psi, st%ob_lead(il)%rho, st%ob_lead(il)%green)
+      nullify(st%ob_lead(il)%intf_psi, st%ob_lead(il)%rho, st%ob_lead(il)%self_energy)
     end do
     nullify(st%ob_eigenval, st%ob_occ)
@@ -917 +919 @@
       do il = 1, NLEADS
         np = mesh%lead_unit_cell(il)%np
-        SAFE_ALLOCATE(st%ob_lead(il)%intf_psi(1:np, 1:2, 1:st%d%dim, st1:st2, k1:k2))
+        SAFE_ALLOCATE(st%ob_lead(il)%intf_psi(1:np, 1:st%d%dim, st1:st2, k1:k2))
         st%ob_lead(il)%intf_psi = M_z0
       end do
@@ -1107 +1109 @@
       SAFE_DEALLOCATE_P(st%ob_occ)
       do il = 1, NLEADS
-        SAFE_DEALLOCATE_P(st%ob_lead(il)%green)
+        SAFE_DEALLOCATE_P(st%ob_lead(il)%self_energy)
       end do
     end if
@@ -2762 +2764 @@
 
   ! ---------------------------------------------------------
-  !> initialize the surface Green`s functions of the leads
-  subroutine states_init_green(st, gr, nspin, d_ispin, lead)
+  !> initialize the self energy of the leads
+  subroutine states_init_self_energy(st, gr, nspin, d_ispin, lead)
     type(states_t),      intent(inout) :: st
     type(grid_t),        intent(in)    :: gr
@@ -2771 +2773 @@
 
     character(len=2)      :: spin
-    character(len=256)    :: fmt, fname_real, fname_imag
+    character(len=256)    :: fmt, fname, fname_real, fname_imag
     FLOAT                 :: energy
     integer  :: np, ik, ist, il, ispin, s1, s2, k1, k2
     integer  :: green_real, green_imag, irow
 
-    call push_sub('states.states_init_green')
-
-    ! Calculate Green`s function of the leads.
+    call push_sub('states.states_init_self_energy')
+
+    ! Calculate self energy of the leads.
     ! FIXME: For spinors, this calculation is almost certainly wrong.
     ASSERT(st%ob_nst == st%nst)
@@ -2786 +2788 @@
     do il = 1, NLEADS
       np = gr%intf(il)%np_intf
-      SAFE_ALLOCATE(st%ob_lead(il)%green(1:np, 1:np, 1:nspin, s1:s2, k1:k2))
+      SAFE_ALLOCATE(st%ob_lead(il)%self_energy(1:np, 1:np, 1:nspin, s1:s2, k1:k2))
     end do
-    call messages_print_stress(stdout, "Lead Green's functions")
+    call messages_print_stress(stdout, "Lead self energy")
     message(1) = ' st#     k#  Spin      Lead     Energy'
     call write_info(1)
@@ -2823 +2825 @@
               trim(spin), trim(LEAD_NAME(il)), energy
             call write_info(1)
-            ! \todo magnetic gs
-            call lead_green(energy, lead(il)%h_diag(:, :, ispin), lead(il)%h_offdiag(:, :), &
-              gr%intf(il), st%ob_lead(il)%green(:, :, ispin, ist, ik), .true.)
-
-            ! Write the entire Green`s function to a file.
+
+            call lead_self_energy(energy, lead(il)%h_diag(:, :, ispin), lead(il)%h_offdiag(:, :), &
+              gr%intf(il), st%ob_lead(il)%self_energy(:, :, ispin, ist, ik), .true.)
+
+            ! Write the entire self energy to a file.
             if(in_debug_mode) then
               call io_mkdir('debug/open_boundaries')
-              write(fname_real, '(3a,i4.4,a,i3.3,a,i1.1,a)') 'debug/open_boundaries/green-', &
+              write(fname_real, '(3a,i4.4,a,i3.3,a,i1.1,a)') 'debug/open_boundaries/self-energy-', &
                 trim(LEAD_NAME(il)), '-', ist, '-', ik, '-', ispin, '.real'
-              write(fname_imag, '(3a,i4.4,a,i3.3,a,i1.1,a)') 'debug/open_boundaries/green-', &
+              write(fname_imag, '(3a,i4.4,a,i3.3,a,i1.1,a)') 'debug/open_boundaries/self-energy-', &
                 trim(LEAD_NAME(il)), '-', ist, '-', ik, '-', ispin, '.imag'
               green_real = io_open(fname_real, action='write', grp=st%d%kpt%mpi_grp, is_tmp=.false.)
@@ -2839 +2841 @@
               write(fmt, '(a,i6,a)') '(', np, 'e24.16)'
               do irow = 1, np
-                write(green_real, fmt) real(st%ob_lead(il)%green(irow, :, ispin, ist, ik))
-                write(green_imag, fmt) aimag(st%ob_lead(il)%green(irow, :, ispin, ist, ik))
+                write(green_real, fmt) real(st%ob_lead(il)%self_energy(irow, :, ispin, ist, ik))
+                write(green_imag, fmt) aimag(st%ob_lead(il)%self_energy(irow, :, ispin, ist, ik))
               end do
               call io_close(green_real); call io_close(green_imag)
@@ -2851 +2853 @@
 
     call pop_sub()
-  end subroutine states_init_green
+  end subroutine states_init_self_energy
+
+
+  ! ---------------------------------------------------------
+  ! write H_(C,apha)*Psi_(alpha) without using Psi_(alpha)
+  subroutine states_write_proj_lead_wf(dir, intf, st)
+    character(len=*),  intent(in) :: dir   ! directory
+    type(interface_t), intent(in) :: intf(:)
+    type(states_t),    intent(in) :: st
+
+    integer :: ik, ist, idim, il, np, ip, iunit
+    CMPLX, allocatable :: psi(:, :), phi(:, :), hpsi(:, :), self_energy(:, :)
+    character(len=256) :: fname
+
+    call push_sub('states.write_proj_lead_wf')
+
+    np = maxval(intf(1:NLEADS)%np_intf)
+
+    SAFE_ALLOCATE(psi(1:np, 1:st%d%dim))
+    SAFE_ALLOCATE(phi(1:np, 1:st%d%dim))
+    SAFE_ALLOCATE(hpsi(1:np, 1:st%d%dim))
+    SAFE_ALLOCATE(self_energy(1:np, 1:np))
+
+    hpsi(:, :) = M_z0
+
+    call io_mkdir(dir, is_tmp=.true.)
+
+    do ik = st%d%kpt%start, st%d%kpt%end
+      do ist = st%st_start, st%st_end
+        do il = 1, NLEADS
+          np = intf(il)%np_intf
+          forall(ip = 1:np)
+            psi(ip, :) = st%zpsi(intf(il)%index(ip), :, ist, ik)
+            phi(ip, :) = st%zphi(intf(il)%index(ip), :, ist, ik)
+          end forall
+          do idim = 1, st%d%dim
+            self_energy(1:np, 1:np) = st%ob_lead(il)%self_energy(1:np, 1:np, idim, ist, ik)
+            call lalg_gemv(np, np, M_z1, self_energy(1:np, 1:np), psi(1:np, idim), M_z0, hpsi(1:np, idim))
+            if((il.eq.LEFT).and.(-st%d%kpoints(1, ik).gt.M_ZERO) .or. (il.eq.RIGHT).and.(-st%d%kpoints(1, ik).lt.M_ZERO)) then
+              ! add the reflecting part
+              self_energy(1:np, 1:np) = transpose(conjg(self_energy(1:np, 1:np))) - self_energy(1:np, 1:np)
+              call lalg_gemv(np, np, M_z1, self_energy(1:np, 1:np), phi(1:np, idim), M_z1, hpsi(1:np, idim))
+            end if
+            write(fname, '(3a,i4.4,a,i3.3,a,i1.1)') 'src0-', trim(LEAD_NAME(il)), '-', ist, '-', ik, '-', idim
+            iunit = io_open(trim(dir)//trim(fname), action='write', form='unformatted', is_tmp=.true.)
+            if(iunit.lt.0) then
+              message(1) = 'Cannot write term for source term to file.'
+              call write_warning(1)
+              call io_close(iunit)
+              call pop_sub(); return
+            end if
+            ! Write parameters.
+            write(iunit) np
+            ! Write matrix.
+            write(iunit) hpsi(1:np, idim)
+            call io_close(iunit)
+          end do
+        end do
+      end do
+    end do
+
+    SAFE_DEALLOCATE_A(psi)
+    SAFE_DEALLOCATE_A(phi)
+    SAFE_DEALLOCATE_A(hpsi)
+    SAFE_DEALLOCATE_A(self_energy)
+
+    call pop_sub()
+  end subroutine states_write_proj_lead_wf
+
+  ! ---------------------------------------------------------
+  subroutine states_read_proj_lead_wf(dir, intf, st, src0)
+    character(len=*),  intent(in)    :: dir   ! directory
+    type(interface_t), intent(in)    :: intf
+    type(states_t),    intent(in)    :: st
+    CMPLX,             intent(inout) :: src0(:, : ,:, :)
+
+    integer :: ik, ist, idim, np, ip, iunit
+    character(len=256) :: fname
+
+    call push_sub('states.states_read_proj_lead_wf')
+
+    do ik = st%d%kpt%start, st%d%kpt%end
+      do ist = st%st_start, st%st_end
+        do idim = 1, st%d%dim
+          ! Try to open file.
+          write(fname, '(3a,i4.4,a,i3.3,a,i1.1)') 'src0-', trim(LEAD_NAME(intf%il)), '-', ist, '-', ik, '-', idim
+          iunit = io_open(trim(dir)//trim(fname), action='read', status='old', die=.false., is_tmp=.true., form='unformatted')
+          if(iunit.lt.0) then ! no file found
+            message(1) = 'Cannot read src(0) from file.'
+            call write_fatal(1)
+          end if
+
+          ! Now read the data.
+          read(iunit) np
+
+          if(np.ne.size(src0, 1)) then
+            message(1) = 'Size mismatch! Cannot read src(0) from file.'
+            call write_fatal(1)
+          end if
+
+          read(iunit) src0(1:np, idim, ist, ik)
+
+          call io_close(iunit)
+        end do
+      end do
+    end do
+
+    call pop_sub()
+  end subroutine states_read_proj_lead_wf
+
 
 end module states_m

trunk/src/td/ob_rti.F90

@@ -220 +220 @@
       SAFE_ALLOCATE(ob%lead(il)%st_intface(1:np, s1:s2, k1:k2, 0:ob%max_mem_coeffs))
       ob%lead(il)%st_intface = M_z0
+      call states_read_proj_lead_wf('open_boundaries/', gr%intf(il), st, ob%lead(il)%src_prev)
     end do
 
@@ -324 +325 @@
             end if
             call calc_source_wf(ob%max_mem_coeffs, m, np, il, hm%lead(il)%h_offdiag, tmp_mem(1:np, 1:np), dt, &
-              st%ob_lead(il)%intf_psi(:, :, 1, ist, ik), ob%src_mem_u(:, il), f0, fac,              &
+              st%ob_lead(il)%intf_psi(:, 1, ist, ik), ob%src_mem_u(:, il), f0, fac,              &
               lambda(m, 0, max_iter, ob%src_mem_u(:, il)), ob%lead(il)%src_prev(:, 1, ist, ik))
             call apply_src(gr%intf(il), ob%lead(il)%src_prev(1:np, 1, ist, ik), st%zpsi(:, :, ist, ik))
@@ -355 +356 @@
         ! h_eff_backward must be a complex symmetric operator !
         call zqmr_sym(gr%mesh%np, st%zpsi(:, 1, ist, ik), tmp(:, 1), h_eff_backward, precond_prop, &
-          qmr_iter, residue=dres, threshold=qmr_tol, showprogress=.false., converged=conv)
+          qmr_iter, residue=dres, threshold=qmr_tol, showprogress=in_debug_mode, converged=conv)
       end do
       if(in_debug_mode) then ! write info
@@ -447 +448 @@
             if(m.gt.0) tmp_mem(1:npo) = tmp_mem(1:npo) + ob%lead(il)%q_sp(1:npo, m-1)
             call calc_source_wf_sp(max_iter, m, np, il, hm%lead(il)%h_offdiag(:, :),     &
-              tmp_mem(1:npo), dt, order, gr%sb%dim, st%ob_lead(il)%intf_psi(:, :, 1, ist, ik), &
+              tmp_mem(1:npo), dt, order, gr%sb%dim, st%ob_lead(il)%intf_psi(:, 1, ist, ik), &
               ob%lead(il)%q_s(:, :, :), ob%lead(il)%sp2full_map, ob%src_mem_u(:, il), f0, fac,   &
               lambda(m, 0, max_iter, ob%src_mem_u(:, il)), ob%lead(il)%src_prev(:, 1, ist, ik))

trunk/src/td/ob_src.F90

@@ -88 +88 @@
     CMPLX,   intent(in)    :: mem(np, np)   ! the effective memory coefficient
     FLOAT,   intent(in)    :: dt            ! Timestep.
-    CMPLX,   intent(in)    :: psi0(:, :)    ! (np, INNER/OUTER).
+    CMPLX,   intent(in)    :: psi0(:)       ! (np).
     CMPLX,   intent(in)    :: u(0:)
     CMPLX,   intent(in)    :: f0
@@ -98 +98 @@
 
     call push_sub('ob_src.calc_source_wf')
-
     if(m.eq.0) then
-      src(1:np) = psi0(1:np, OUTER)
-      if (il.eq.LEFT) then
-        call ztrmv('U', 'N', 'N', np, offdiag, np, src, 1)
-      else
-        call ztrmv('L', 'N', 'N', np, offdiag, np, src, 1)
-      end if
-      tmp = -M_zI*dt*u(0)*f0
-!      call lalg_gemv(np, np, tmp*M_zI*M_HALF*dt, mem, psi0(1:np, INNER), tmp, src)
-      call zsymv('U', np, tmp*M_zI*M_HALF*dt, mem, np, psi0(1:np, INNER), 1, tmp,  src, 1)
+      ! initial src is V*psi0(outer) (precalculated in gs mode)
+      tmp = -M_zI*dt*f0*u(0)
+      !call lalg_gemv(np, np, tmp*M_zI*M_HALF*dt, mem, psi0, tmp, src)
+      call zsymv('U', np, tmp*M_zI*M_HALF*dt, mem, np, psi0, 1, tmp, src, 1)
     else
       tmp   = factor*u(m)*u(m-1)
@@ -115 +109 @@
       else
         alpha = M_HALF*dt**2*f0*lambda/u(m)
-!        call lalg_gemv(np, np, alpha, mem, psi0(1:np, INNER), tmp, src)
-        call zsymv('U', np, alpha, mem, np, psi0(1:np, INNER), 1, tmp, src, 1)
+!        call lalg_gemv(np, np, alpha, mem, psi0, tmp, src)
+        call zsymv('U', np, alpha, mem, np, psi0, 1, tmp, src, 1)
       end if
     end if
@@ -138 +132 @@
     CMPLX,   intent(in)    :: sp_mem(:)
     FLOAT,   intent(in)    :: dt             ! Timestep.
-    CMPLX,   intent(in)    :: psi0(:, :)     ! (np, INNER/OUTER)
+    CMPLX,   intent(in)    :: psi0(:)        ! (np)
     CMPLX,   intent(in)    :: mem_s(:, :, :)
     integer, intent(in)    :: mapping(:)

trunk/src/td/pes_rc_inc.F90

@@ -38 +38 @@
   !%Section Time-Dependent::PES
   !%Description
-  !% List of points where to calculate the photoelectron spectrum by Suraud's method.
+  !% List of points where to calculate the photoelectron spectrum by Suraud`s method.
   !% The exact syntax is:
   !%

trunk/src/td/td.F90

@@ -394 +394 @@
           fromScratch = .true.
         end if
-      end if
-
-      if(.not.fromscratch .and. st%open_boundaries) then
-        call restart_read_ob_intf(trim(restart_dir)//GS_DIR, st, gr, ierr)
+        ! extract the interface wave function
+        if(st%open_boundaries) call restart_get_ob_intf(st, gr)
       end if
 
@@ -432 +430 @@
 
         if(.not. st%only_userdef_istates) then
-          ! In the open-boundary case the ground-state wavefunctions are bit too "wide".
-          ! Therefore, we need a special routine to extract the middle.
-          if(gr%sb%open_boundaries) then
-            call restart_read_ob_intf(trim(restart_dir)//GS_DIR, st, gr, ierr)
-            if(ierr.ne.0) then
-              message(1) = "Could not read interface wavefunctions from '"//trim(restart_dir)//GS_DIR//"'"
-              message(2) = "Please run an open-boundaries ground-state calculation first!"
-              call write_fatal(2)
-            end if
-            call restart_read_ob_central(trim(restart_dir)//GS_DIR, st, gr, ierr)
-          else
-            call restart_read(trim(restart_dir)//GS_DIR, st, gr, geo, ierr)
-          end if
+          call restart_read(trim(restart_dir)//GS_DIR, st, gr, geo, ierr)
           if(ierr.ne.0) then
             message(1) = "Could not read KS orbitals from '"//trim(restart_dir)//GS_DIR//"'"
             message(2) = "Please run a ground-state calculation first!"
             call write_fatal(2)
+          end if
+          
+          if(gr%sb%open_boundaries) then
+            ! extract the interface wave function
+            call restart_get_ob_intf(st, gr)
           end if
         end if

trunk/testsuite/open_systems/01-wavepacket.test

@@ -20 +20 @@
 match ; Density [step   0] ; LINE(flat_td.general/multipoles,  16, 30) ; 3.981621426799e+00
 match ; Density [step 100] ; LINE(flat_td.general/multipoles, 116, 30) ; 3.954076302565e+00
-match ; Density [step 200] ; LINE(flat_td.general/multipoles, 216, 30) ; 2.434921125600e+00
-match ; Density [step 300] ; LINE(flat_td.general/multipoles, 316, 30) ; 9.635221759786e-01
+match ; Density [step 200] ; LINE(flat_td.general/multipoles, 216, 30) ; 2.434921125591e+00
+match ; Density [step 300] ; LINE(flat_td.general/multipoles, 316, 30) ; 9.635221759715e-01

trunk/testsuite/open_systems/03-extended_eigenstate.test

@@ -9 +9 @@
 # Calculate eigenstate of 1d attractive square potential barrier.
 Input: 03-extended_eigenstate.01-square_well_1d.inp
-match ; Total energy ; GREP(well_static/info, 'Total       =', 20) ; 5.7224
+match ; Total energy ; GREP(well_static/info, 'Total       =', 20) ; 5.0117
 
 # Calculate eigenstates of 2D ring potential attached to upside-down Gaussian
 # shaped lead channels.
 Input: 03-extended_eigenstate.02-ring_leads_2d.inp
-match ; Total energy ; GREP(ring_lead_static/info, 'Total       =', 20) ; 81.6571
+match ; Total energy ; GREP(ring_lead_static/info, 'Total       =', 20) ; 58.9977

trunk/testsuite/open_systems/04-propagate_eigenstate.test

@@ -9 +9 @@
 # Calculate eigenstate of 1d attractive square-potential barrier.
 Input: 04-propagate_eigenstate.01-square_well_1d.inp
-match ; Total energy [step   0]; LINE(well_td.general/energy,   6, 30) ; 2.049018304223e+00
-match ; Total energy [step 100]; LINE(well_td.general/energy, 106, 30) ; 2.049017944732e+00
-match ; Total energy [step 200]; LINE(well_td.general/energy, 206, 30) ; 2.049017813254e+00
-match ; Total energy [step 300]; LINE(well_td.general/energy, 306, 30) ; 2.049017923834e+00
+match ; Total energy [step   0]; LINE(well_td.general/energy,   6, 30) ; 2.064297127468e+00
+match ; Total energy [step 100]; LINE(well_td.general/energy, 106, 30) ; 2.064289561608e+00
+match ; Total energy [step 200]; LINE(well_td.general/energy, 206, 30) ; 2.064286304422e+00
+match ; Total energy [step 300]; LINE(well_td.general/energy, 306, 30) ; 2.064284239477e+00
 
 # Calculate eigenstates of 2D ring potential attached to upside-down Gaussian-shaped
 # lead channels and propagate with a bias switched on.
 Input: 04-propagate_eigenstate.02-ring_leads_2d.inp
-match ; Total energy [step   0]; LINE(ring_lead_td.general/energy,   6, 30) ; 4.626426074917e+02
-match ; Total energy [step  25]; LINE(ring_lead_td.general/energy,  31, 30) ; 4.635193147887e+02
-match ; Total energy [step  50]; LINE(ring_lead_td.general/energy,  56, 30) ; 4.648048709782e+02
-match ; Total energy [step  75]; LINE(ring_lead_td.general/energy,  81, 30) ; 4.676178752149e+02
-match ; Total energy [step 100]; LINE(ring_lead_td.general/energy, 106, 30) ; 4.714202246759e+02
+match ; Total energy [step   0]; LINE(ring_lead_td.general/energy,   6, 30) ; 4.626426040993e+02
+match ; Total energy [step  25]; LINE(ring_lead_td.general/energy,  31, 30) ; 4.635193115106e+02
+match ; Total energy [step  50]; LINE(ring_lead_td.general/energy,  56, 30) ; 4.648048674987e+02
+match ; Total energy [step  75]; LINE(ring_lead_td.general/energy,  81, 30) ; 4.676178717766e+02
+match ; Total energy [step 100]; LINE(ring_lead_td.general/energy, 106, 30) ; 4.714202213930e+02