Added accuracy control to transiesta on a per-electrode case

Through fdf-flags one may control the accuracy required for the
convergence of the self-energies. This enables finer control of
the self-energy calculations and in some cases it may be useful.

Src/m_ts_electrode.F90

@@ -21,8 +21,6 @@ module m_ts_electrode
 
   private
 
-  ! Accuracy of the surface-Green function
-  real(dp), parameter :: accur = 1.e-15_dp
   ! BLAS parameters
   complex(dp), parameter :: z_1  = dcmplx(1._dp,0._dp)
   complex(dp), parameter :: z_m1 = dcmplx(-1._dp,0._dp)
@@ -39,7 +37,7 @@ contains
   ! Calculates the surface Green function for the electrodes
   ! Handles both the left and right one
   ! this is the Sancho, Sancho and Rubio algorithm
-  subroutine SSR_sGreen_DOS(no,ZE,H00,S00,H01,S01,GS, &
+  subroutine SSR_sGreen_DOS(no,ZE,H00,S00,H01,S01,accu, GS, &
        DOS, T, &
        nwork, zwork, &
        iterations, final_invert)
@@ -51,9 +49,11 @@ contains
 ! complex(dp) S00     : Overlap matrix within the first unit cell (discarding T-direction)
 ! complex(dp) H01     : Transfer matrix from H00 to the neighbouring cell (in T-direction)
 ! complex(dp) S01     : Transfer matrix from S00 to the neighbouring cell (in T-direction)
+! real(dp) accu       : Define the accuracy needed for convergence.
 ! ***************** OUTPUT *********************************************
 ! complex(dp) GS      : Surface Green function of the electrode
 ! real(dp) DOS        : DOS of bulk electrode (additive)
+! real(dp) T          : transmission of bulk electrode (additive)
 ! **********************************************************************
     use m_pivot_array, only : ipiv
     use m_mat_invert
@@ -68,8 +68,9 @@ contains
     complex(dp), intent(in) :: ZE 
     complex(dp), intent(in) :: H00(no*no),S00(no*no)
     complex(dp), intent(in) :: H01(no*no),S01(no*no)
+    real(dp), intent(in) :: accu
 
-    integer,     intent(in) :: nwork
+    integer, intent(in) :: nwork
 
     logical, intent(in), optional :: final_invert
 
@@ -174,9 +175,9 @@ contains
 !$OMP end parallel
 
     ! Initialize loop
-    ro = accur + 1._dp
+    ro = accu + 1._dp
     as_first = .false.
-    do while ( ro > accur ) 
+    do while ( ro > accu ) 
 
        ! Increment iterations
        if ( present(iterations) ) &
@@ -507,7 +508,7 @@ contains
   ! Calculates the surface Green function for the electrodes
   ! Handles both the left and right one
   ! this is the Sancho, Sancho and Rubio algorithm
-  subroutine SSR_sGreen_NoDOS(no,ZE,H00,S00,H01,S01,GS, &
+  subroutine SSR_sGreen_NoDOS(no,ZE,H00,S00,H01,S01,accu, GS, &
        nwork, zwork, &
        iterations, final_invert)
        
@@ -532,6 +533,7 @@ contains
     complex(dp), intent(in) :: ZE 
     complex(dp), intent(in) :: H00(no*no),S00(no*no)
     complex(dp), intent(in) :: H01(no*no),S01(no*no)
+    real(dp), intent(in) :: accu
 
     integer,     intent(in) :: nwork
 
@@ -619,9 +621,9 @@ contains
 !$OMP end parallel
 
     ! Initialize loop
-    ro = accur + 1._dp
+    ro = accu + 1._dp
     as_first = .false.
-    do while ( ro > accur ) 
+    do while ( ro > accu ) 
 
        ! Increment iterations
        if ( present(iterations) ) &
@@ -1193,7 +1195,8 @@ contains
                 ! Zenergy is wrt. to the system Fermi-level
                 if ( CalcDOS ) then
                    lDOS = 0._dp
-                   call SSR_sGreen_DOS(nuo_E,ZEnergy,H00,S00,H01,S01,GS, &
+                   call SSR_sGreen_DOS(nuo_E,ZEnergy,H00,S00,H01,S01, &
+                        El%accu, GS, &
                         lDOS,i_mean,9*nS,zwork, &
                         iterations=iters(iqpt,iEn,ikpt,1), final_invert = reduce_size)
                    
@@ -1202,7 +1205,8 @@ contains
                    if ( CalcT ) T(iEn,ispin) = T(iEn,ispin) + i_mean
 
                 else
-                   call SSR_sGreen_NoDos(nuo_E,ZEnergy,H00,S00,H01,S01,GS, &
+                   call SSR_sGreen_NoDos(nuo_E,ZEnergy,H00,S00,H01,S01, &
+                        El%accu, GS, &
                         8*nS,zwork, &
                         iterations=iters(iqpt,iEn,ikpt,1), final_invert = reduce_size)
                    
@@ -1934,12 +1938,14 @@ contains
 
        ! calculate the contribution for this q-point
        if ( calc_DOS ) then
-          call SSR_sGreen_DOS(nuo_E,Z,H00,S00,H01,S01,GS, &
+          call SSR_sGreen_DOS(nuo_E,Z,H00,S00,H01,S01, &
+               El%accu, GS, &
                DOS(1:nuo_E), T, &
                nw, zwork, &
                final_invert = reduce_size)
        else
-          call SSR_sGreen_NoDOS(nuo_E,Z,H00,S00,H01,S01,GS, &
+          call SSR_sGreen_NoDOS(nuo_E,Z,H00,S00,H01,S01, &
+               El%accu, GS, &
                nw, zwork, &
                final_invert = reduce_size)
        end if

Src/m_ts_electype.F90

@@ -146,8 +146,13 @@ module m_ts_electype
      ! and NOT: i (Sigma - Sigma^\dagger)
      complex(dp), pointer :: Gamma(:), Sigma(:)
 
+     ! The accuracy required for the self-energy
+     !  == 1e-13 * eV
+     real(dp) :: accu = 7.349806700083788e-15_dp
+
      ! The imaginary part in the electrode
-     real(dp) :: Eta = 7.3498067e-6_dp ! corresponds to 0.0001 eV
+     !  == 0.0001 * eV
+     real(dp) :: Eta = 7.3498067e-6_dp
 
      ! The region of the down-folded region
      type(tRgn) :: o_inD, inDpvt
@@ -558,6 +563,14 @@ contains
           if ( fdf_bnnames(pline) < 2 ) call die('TSDE name not supplied')
           this%DEfile = trim(fdf_bnames(pline,2))
 
+#ifdef TBTRANS
+       else if ( leqi(ln,'tbt.Accuracy') .or. leqi(ln,'Accuracy') ) then
+#else
+       else if ( leqi(ln,'Accuracy') ) then
+#endif
+          call pline_E_parse(pline,1,ln, &
+               val = this%accu, before=3)
+
 #ifdef TBTRANS
        else if ( leqi(ln,'tbt.Eta') .or. leqi(ln,'Eta') ) then
 #else
@@ -2023,6 +2036,7 @@ contains
 #else
     write(*,f9)  '  Electrode self-energy imaginary Eta', this%Eta/eV,' eV'
 #endif
+    write(*,f9)  '  Electrode self-energy accuracy', this%accu/eV,' eV'
     write(*,f6)  '  Electrode inter-layer distance (semi-inf)', this%dINF_layer/Ang,' Ang'
 
 

Src/m_ts_options.F90

@@ -428,6 +428,8 @@ contains
     ! the remaining electrodes have their chemical potential at 0
     ! We should probably warn if +2 electrodes are used and t_dir is the
     ! same for all electrodes... Then the user needs to know what (s)he is doing...
+    ! Accuracy required for self-energy convergence
+    Elecs(:)%accu = fdf_get('TS.Elecs.Accuracy',1e-14_dp*eV,'Ry')
     Elecs(:)%Eta  = fdf_get('TS.Contours.nEq.Eta',0.0001_dp*eV,'Ry')
     Elecs(:)%Eta  = fdf_get('TS.Elecs.Eta',Elecs(1)%Eta,'Ry')
     Elecs(:)%Bulk = fdf_get('TS.Elecs.Bulk',.true.) ! default everything to bulk electrodes

Util/TS/TBtrans/m_tbt_options.F90

@@ -266,6 +266,10 @@ contains
     rtmp = rtmp ** 2
 #endif
     Elecs(:)%Eta = rtmp
+    rtmp = 1.e-14_dp * eV
+    rtmp = fdf_get('TS.Elecs.Accuracy',rtmp,'Ry')
+    rtmp = fdf_get('TBT.Elecs.Accuracy',rtmp,'Ry')
+    Elecs(:)%accu = rtmp
 
     ! whether all calculations should be performed
     ! "out-of-core" i.e. whether the GF files should be created or not

Util/TS/TBtrans/m_tbt_sigma_save.F90

@@ -383,6 +383,9 @@ contains
 
        dic = ('info'.kv.'Imaginary part of self-energy')//('unit'.kv.'Ry')
        call ncdf_def_var(grp,'Eta',NF90_DOUBLE,(/'one'/), atts = dic)
+       
+       dic = ('info'.kv.'Accuracy of the self-energy')//('unit'.kv.'Ry')
+       call ncdf_def_var(grp,'Accuracy',NF90_DOUBLE,(/'one'/), atts = dic)
        call delete(dic)
 
        call ncdf_def_dim(grp,'no_e',Elecs(iEl)%o_inD%n)
@@ -399,6 +402,7 @@ contains
        call delete(dic)
 
        call ncdf_put_var(grp,'Eta',Elecs(iEl)%Eta)
+       call ncdf_put_var(grp,'Accuracy',Elecs(iEl)%accu)
        call ncdf_put_var(grp,'pivot',Elecs(iEl)%o_inD%r)
 
     end do

version.info

-trunk-536
+trunk-536--transiesta-1