Merge lp:~nickpapior/siesta/4.1-superc-gamma into lp:~albertog/siesta/4.1-superc-gamma

[Nick Papior]

Reinsert the Gamma variable in sub-routines for clarity.

I have reworked the patch to also include the Gamma variables for the corresponding k-point lists.
This should make it easier to follow the combinations.

I have also changed the keywords in transiesta, for clarity.

[Alberto Garcia (albertog)]

It looks good. Indeed more clear.

The only remaining problem is that, I think, the proper logic for choosing between dminim and zminim in the OMM method is based on the supercell (no_u == no_s) and not on whether the gamma point only is used. In other words. dminim will not work if we use a supercell with the gamma point.
There are places where it is obvious where to put a call to "MODP", but in other cases I am not sure.

So pending a possible optimization of this, I am going to merge retaining the "supercell" logic for this in compute_dm and post_scf_work.

[Nick Papior (nickpapior)]

Sorry for changing the m_dminim call. I can now see that you are correct, if a supercell is used, then OMM should use zminim.

=== modified file 'Src/Makefile'
--- Src/Makefile	2018-06-20 10:34:33 +0000
+++ Src/Makefile	2018-06-25 11:35:23 +0000
@@ -967,7 +967,7 @@
 m_ts_electrode.o: precision.o sys.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Src/bands.F'
--- Src/bands.F	2018-06-12 22:59:08 +0000
+++ Src/bands.F	2018-06-25 11:35:23 +0000
@@ -395,7 +395,6 @@
       integer :: ik, il, io, ispin, iu, iu_wfs, iuo, naux, nhs, j
 
       logical :: SaveParallelOverK
-      logical :: gamma
       
       real(dp)
      .  Dnew, qs(2), e1, e2, efs(2), emax, emin, Enew, eV, qk, qtot,
@@ -419,9 +418,6 @@
 C Start time counter
       call timer( 'bands', 1 )
 
-      ! This routine always uses diag*k routines
-      gamma = .false.
-      
 C Check parameter maxk
       if (nk .gt. maxk) then
         if (Node.eq.0) then
@@ -451,7 +447,7 @@
 
         rewind (iu_wfs)
 
-        write(iu_wfs) nk, gamma
+        write(iu_wfs) nk, .false. ! nk, Gamma, same file-format in WFS as for Gamma-point
         write(iu_wfs) nspin
         write(iu_wfs) no_u
         write(iu_wfs) (iaorb(j),labelfis(isa(iaorb(j))),
 
=== modified file 'Src/compute_dm.F'
--- Src/compute_dm.F	2018-06-20 10:40:25 +0000
+++ Src/compute_dm.F	2018-06-25 11:35:23 +0000
@@ -66,8 +66,6 @@
       real(dp)            :: buffer1
       integer             :: mpierr
 
-      logical             :: not_using_auxcell
-
 !-------------------------------------------------------------------- BEGIN
 
       if (SIESTA_worker) call timer( 'compute_dm', 1 )
@@ -92,8 +90,6 @@
 !$OMP end parallel
       end if
 
-      not_using_auxcell = (no_u == no_s)
-      
       ! Compute shift in Tr(H*DM) for fermi-level bracketting
       ! Use the current H, the previous iteration H, and the
       ! previous iteration DM
@@ -178,14 +174,13 @@
      &              no_l, maxnh, maxnh, no_u,
      &              numh, listhptr, listh, numh, listhptr, listh, 
      &              H, S, qtot, fixspin, qtots, temp, e1, e2,
-     $              xijo, indxuo, nkpnt, kpoint, kweight,
+     $              xijo, indxuo, gamma_SCF, nkpnt, kpoint, kweight,
      &              eo, qo, Dscf, Escf, ef, efs, Entropy, no_u,
      &              occtol, iscf, neigwanted)
         Ecorrec = 0.0_dp
         PreviousCallDiagon=.true.
       elseif (isolve .eq. SOLVE_ORDERN) then
-        if ((nkpnt > 1) .or. (sum(abs(kpoint(:,1))) /= 0.0))
-     $        call die("Cannot do O(N) with k-points.")
+        if ( .not. gamma_SCF ) call die("Cannot do O(N) with k-points.")
         if ( spin%NCol .or. spin%SO )
      .      call die("Cannot do O(N) with non-coll spins or Spin-orbit")
         call ordern(usesavelwf, ioptlwf, na_u, no_u, no_l, lasto,
@@ -204,7 +199,7 @@
         ! Test based on use of auxiliary supercell
         ! It might still be possible to avoid the complex version
         
-        if (not_using_auxcell) then 
+        if ( gamma_SCF ) then 
           call dminim(.false., PreviousCallDiagon, iscf, istp, no_l,
      &                 spin%H, no_u, maxnh, numh, listhptr, listh, Dscf,
      &                 eta, qtots, H, S, H_kin)
@@ -222,7 +217,7 @@
      &              no_l, maxnh, maxnh, no_u,
      &              numh, listhptr, listh, numh, listhptr, listh,
      &              H, S, qtot, fixspin, qtots, temp, e1, e2,
-     $              xijo, indxuo, nkpnt, kpoint, kweight,
+     $              xijo, indxuo, gamma_SCF, nkpnt, kpoint, kweight,
      &              eo, qo, Dscf, Escf, ef, efs, Entropy, no_u,
      &              occtol, iscf, neigwanted)
 
@@ -231,7 +226,7 @@
       else if (TSrun) then
 
          call transiesta(iscf,spin%H, block_dist, sparse_pattern,
-     &        Not_Using_Auxcell, ucell, nsc, isc_off, no_u, na_u,
+     &        no_u == no_s, ucell, nsc, isc_off, no_u, na_u,
      &        lasto, xa, maxnh,
      &        H, S, Dscf, Escf, Ef, Qtot, .false., DE_NEGF )
 
 
=== modified file 'Src/diag2g.F'
--- Src/diag2g.F	2018-06-20 10:34:33 +0000
+++ Src/diag2g.F	2018-06-25 11:35:23 +0000
@@ -20,7 +20,7 @@
       use parallel,        only : Node, Nodes, BlockSize
       use parallelsubs,    only : LocalToGlobalOrb,GlobalToLocalOrb
       use m_fermid,        only : fermid, stepf
-      use intrinsic_missing, only: MODP1
+      use intrinsic_missing, only: MODP
 
 #ifdef MPI
       use mpi_siesta
@@ -163,7 +163,7 @@
         do j = 1,numh(io)
           ind = listhptr(io) + j
           jo = listh(ind)
-          jo = MODP1(jo,nuotot) ! To allow auxiliary supercells
+          jo = MODP(jo,nuotot) ! To allow auxiliary supercells
           Saux(1,jo,1,io) = Saux(1,jo,1,io) + dcmplx( S(ind), 0.0_dp)
           Saux(2,jo,2,io) = Saux(2,jo,2,io) + dcmplx( S(ind), 0.0_dp)
           Haux(1,jo,1,io) = Haux(1,jo,1,io) + dcmplx(H(ind,1),   0.0_dp)
@@ -189,7 +189,7 @@
           do j = 1,numh(io)
            ind = listhptr(io) + j
            jo = listh(ind)
-           jo = MODP1(jo,nuotot) ! To allow auxiliary supercells
+           jo = MODP(jo,nuotot) ! To allow auxiliary supercells
            Saux(1,jo,1,io) = Saux(1,jo,1,io) +dcmplx( S(ind), 0.0_dp)
            Saux(2,jo,2,io) = Saux(2,jo,2,io) +dcmplx( S(ind), 0.0_dp)
            Haux(1,jo,1,io) = Haux(1,jo,1,io) +dcmplx(H(ind,1),   0.0_dp)
@@ -269,7 +269,7 @@
             do j = 1,numd(io)
               ind = listdptr(io) + j
               jo = listd(ind)
-              jo = MODP1(jo,nuotot) ! To allow auxiliary supercells
+              jo = MODP(jo,nuotot) ! To allow auxiliary supercells
       
 !                 | ------- 1,1 -------     ------- 2,1 ------- |
 !                 | c_{j,up}^* c_{i,up}     c_{j,dn}^* c_{i,up) |
 
=== modified file 'Src/diag3g.F'
--- Src/diag3g.F	2018-06-20 10:34:33 +0000
+++ Src/diag3g.F	2018-06-25 11:35:23 +0000
@@ -17,7 +17,7 @@
       use parallel,        only : Node, Nodes, BlockSize
       use parallelsubs,    only : LocalToGlobalOrb,GlobalToLocalOrb
       use m_fermid,        only : fermid, stepf
-      use intrinsic_missing, only: MODP1
+      use intrinsic_missing, only: MODP
 #ifdef MPI
       use mpi_siesta
 #endif
@@ -155,7 +155,7 @@
         do j = 1,numh(io)
           ind = listhptr(io) + j
           jo = listh(ind)
-          jo = MODP1(jo,nuotot) ! To allow auxiliary supercells
+          jo = MODP(jo,nuotot) ! To allow auxiliary supercells
           Saux(1,jo,1,io) = Saux(1,jo,1,io) + dcmplx( S(ind), 0.0_dp)
           Saux(2,jo,2,io) = Saux(2,jo,2,io) + dcmplx( S(ind), 0.0_dp)
           Haux(1,jo,1,io) = Haux(1,jo,1,io) + dcmplx(H(ind,1), H(ind,5))
@@ -180,7 +180,7 @@
         do j = 1,numh(io)
           ind = listhptr(io) + j
           jo = listh(ind)
-          jo = MODP1(jo,nuotot) ! To allow auxiliary supercells
+          jo = MODP(jo,nuotot) ! To allow auxiliary supercells
           Saux(1,jo,1,io) = Saux(1,jo,1,io) + dcmplx( S(ind), 0.0_dp)
           Saux(2,jo,2,io) = Saux(2,jo,2,io) + dcmplx( S(ind), 0.0_dp)
           Haux(1,jo,1,io) = Haux(1,jo,1,io) + dcmplx(H(ind,1), H(ind,5))
@@ -264,7 +264,7 @@
             do j = 1,numd(io)
               ind = listdptr(io) + j
               jo = listd(ind)
-              jo = MODP1(jo,nuotot) ! To allow auxiliary supercells
+              jo = MODP(jo,nuotot) ! To allow auxiliary supercells
 !------- 1,1 -----------------------------------------------------------
               D11  = caux(1,iio) * dconjg(caux(1,jo))
 !------- 2,2 -----------------------------------------------------------
 
=== modified file 'Src/diagg.F'
--- Src/diagg.F	2018-06-20 10:34:33 +0000
+++ Src/diagg.F	2018-06-25 11:35:23 +0000
@@ -90,7 +90,7 @@
       use m_fermid,      only : fermid, fermispin, stepf
       use m_spin,        only : spinor_dim, e_spin_dim
       use alloc
-      use intrinsic_missing, only: MODP1
+      use intrinsic_missing, only: MODP
 
 #ifdef MPI
       use mpi_siesta
@@ -158,7 +158,7 @@
           do j = 1,numh(io)
             ind = listhptr(io) + j
             jo = listh(ind)
-            jo = MODP1(jo,nuotot)   ! To allow auxiliary supercells
+            jo = MODP(jo,nuotot)   ! To allow auxiliary supercells
             Saux(jo,io) = Saux(jo,io) + S(ind)
             Haux(jo,io) = Haux(jo,io) + H(ind,ispin)
           enddo
@@ -184,7 +184,7 @@
             do j = 1,numh(io)
               ind = listhptr(io) + j
               jo = listh(ind)
-              jo = MODP1(jo,nuotot) ! To allow auxiliary supercells
+              jo = MODP(jo,nuotot) ! To allow auxiliary supercells
               Saux(jo,io) = Saux(jo,io) + S(ind)
               Haux(jo,io) = Haux(jo,io) + H(ind,ispin)
             enddo
@@ -301,7 +301,7 @@
                 do j = 1,numd(io)
                   ind = listdptr(io) + j
                   jo  = listd(ind)
-                  jo = MODP1(jo,nuotot) ! To allow auxiliary supercells
+                  jo = MODP(jo,nuotot) ! To allow auxiliary supercells
                   Dnew(ind,ispin) = Dnew(ind,ispin) + qei*paux(jo)
                   Enew(ind,ispin) = Enew(ind,ispin) + eei*paux(jo)
                 enddo
 
=== modified file 'Src/diagon.F'
--- Src/diagon.F	2018-06-12 22:59:08 +0000
+++ Src/diagon.F	2018-06-25 11:35:23 +0000
@@ -15,7 +15,7 @@
       subroutine diagon(no, nspin, maxuo, maxnh, maxnd, 
      .                  maxo, numh, listhptr, listh, numd, 
      .                  listdptr, listd, H, S, qtot, fixspin, 
-     .                  qs, temp, e1, e2, xij, indxuo, nk, 
+     .                  qs, temp, e1, e2, xij, indxuo, gamma, nk, 
      .                  kpoint, wk, eo, qo, Dnew, Enew, ef, efs, 
      .                  Entropy, nuotot, occtol, iscf, neigwanted )
 C *********************************************************************
@@ -62,6 +62,7 @@
 C                               Unit cell orbitals must be the first in
 C                               orbital lists, i.e. indxuo.le.nuo, with
 C                               nuo the number of orbitals in unit cell
+C logical Gamma               : Whether only the Gamma point is sampled.
 C integer nk                  : Number of k points
 C real*8  kpoint(3,nk)        : k point vectors
 C real*8  wk(nk)              : k point weights (must sum one)
@@ -115,6 +116,7 @@
       implicit none
 
       real(dp), intent(in) :: H(:,:)
+      logical, intent(in) :: gamma
 
       integer
      .  iscf, maxnd, maxnh, maxuo, maxo, nk, no, nuotot,
@@ -145,7 +147,6 @@
 #endif
       integer            :: io, iuo, naux, nhs, npsi, nuo
       real(dp),  pointer :: aux(:)
-      logical            :: gamma
 C  ....................
 
 C Get Node number and calculate local orbital range
@@ -158,8 +159,6 @@
 C Start time counter ................................................
       call timer( 'diagon', 1 )
 
-      gamma = ((nk == 1) .and. (sum(abs(kpoint(:,1))) == 0.0))
-
 C     Check internal dimensions ..........................................
       
       if ( spin%none .or. spin%Col ) then
 
=== modified file 'Src/intrinsic_missing.F90'
--- Src/intrinsic_missing.F90	2018-06-20 10:34:33 +0000
+++ Src/intrinsic_missing.F90	2018-06-25 11:35:23 +0000
@@ -96,7 +96,6 @@
 
 ! Elemental functions (can be called on arrays)
   public :: MODP
-  public :: MODP1
 
 ! Missing matrix stuff
   public :: EYE
@@ -232,29 +231,9 @@
   elemental function MODP(a,p)
     integer, intent(in) :: a,p
     integer :: MODP
-    if ( a > p ) then
-       MODP = MOD(a,p)
-       if ( MODP == 0 ) MODP = p
-    else
-       MODP = a
-    end if
+    MODP = MOD(a-1,p) + 1
   end function MODP
 
-! Another implementation of the above function
-! NOTE that these functions should only be used on NATURAL numbers
-! The target domain is [1,p] instead of [0,p-1]
-! In particular, modp1(p,p) = p
-! This makes it useful in FORTRAN do-loops which are not zero-based
-! but 1-based.
-
-  elemental function MODP1(a,p)
-    integer, intent(in) :: a,p
-    integer :: MODP1
-
-    MODP1 = MOD(a-1,p) + 1
-    
-  end function MODP1
-
 ! Function to return the unique COUNT of an integer array.
 ! Thus will return how many DIFFERENT entries there exists.
   pure function UNIQC(array)
 
=== modified file 'Src/local_DOS.F'
--- Src/local_DOS.F	2018-06-12 22:59:08 +0000
+++ Src/local_DOS.F	2018-06-25 11:35:23 +0000
@@ -78,7 +78,7 @@
           call diagon(no_s, spinor_dim, no_l, maxnh, maxnh, no_u,
      .                numh, listhptr, listh, numh, listhptr, listh, 
      .                H, S, qtot, fixspin, qtots, temp, e1, e2,
-     .                xijo, indxuo, nkpnt, kpoint, kweight,
+     .                xijo, indxuo, gamma_SCF, nkpnt, kpoint, kweight,
      .                eo, qo, Dscf, Escf, ef, efs, dummy_Entrop, no_u,
      .                occtol, dummy_iscf, neigwanted)
 
 
=== modified file 'Src/m_transiesta.F90'
--- Src/m_transiesta.F90	2018-04-04 13:22:33 +0000
+++ Src/m_transiesta.F90	2018-06-25 11:35:23 +0000
@@ -46,7 +46,7 @@
 
   subroutine transiesta(TSiscf,nspin, &
        sp_dist, sparse_pattern, &
-       Gamma, ucell, nsc, isc_off, no_u, na_u, lasto, xa, n_nzs, &
+       no_aux_cell, ucell, nsc, isc_off, no_u, na_u, lasto, xa, n_nzs, &
        H, S, DM, EDM, Ef, &
        Qtot, Fermi_correct, DE_NEGF)
 
@@ -82,7 +82,7 @@
     integer, intent(in)  :: nspin
     type(OrbitalDistribution), intent(inout) :: sp_dist
     type(Sparsity), intent(inout) :: sparse_pattern
-    logical, intent(in)  :: Gamma
+    logical, intent(in)  :: no_aux_cell
     real(dp), intent(in) :: ucell(3,3)
     integer, intent(in)  :: nsc(3), no_u, na_u
     integer, intent(in) :: isc_off(3,product(nsc))
@@ -134,7 +134,7 @@
        ! local sparsity pattern...
        converged = IsVolt .or. TS_RHOCORR_METHOD == TS_RHOCORR_FERMI
        call ts_sparse_init(slabel,converged, N_Elec, Elecs, &
-            ucell, nsc, na_u, xa, lasto, sp_dist, sparse_pattern, Gamma, &
+            ucell, nsc, na_u, xa, lasto, sp_dist, sparse_pattern, no_aux_cell, &
             isc_off)
 
        if ( ts_method == TS_BTD ) then
 
=== modified file 'Src/pdos.F'
--- Src/pdos.F	2018-06-12 22:59:08 +0000
+++ Src/pdos.F	2018-06-25 11:35:23 +0000
@@ -8,7 +8,7 @@
       subroutine pdos( NO, nspin, maxspn, NO_L, MAXNH, 
      .                 MAXO, NUMH, LISTHPTR, LISTH, H, S, 
      .                 E1, E2, SIGMA, NHIST,
-     .                 XIJ, INDXUO, NK, KPOINT, WK, EO, NO_U )
+     .                 XIJ, INDXUO, GAMMA, NK, KPOINT, WK, EO, NO_U )
 C **********************************************************************
 C Subroutine to calculate the projected density of states on the
 C atomic orbitals for a given eigenvalue spectra
@@ -43,6 +43,7 @@
 C                               Unit cell orbitals must be the first in
 C                               orbital lists, i.e. indxuo.le.nuo, with
 C                               nuo the nuber of orbitals in the unit cell
+C logical Gamma               : whether only the Gamma point is sampled
 C INTEGER NK                  : Number of k points
 C REAL*8  KPOINT(3,NK)        : k point vectors
 C REAL*8  WK(NK)              : k point weights (must sum one)
@@ -79,6 +80,7 @@
      .  NO, NSPIN, MAXSPN, NO_L, MAXNH, NK, NHIST, 
      .  MAXO, NO_U
 
+      logical, intent(in) :: Gamma
       integer 
      .  NUMH(*), LISTH(MAXNH), LISTHPTR(*), INDXUO(NO)
 
@@ -98,7 +100,6 @@
       integer iat, spec, ii, iorb
 
       logical :: orig_ParallelOverK, orig_Serial
-      logical :: gamma
 
       real(dp), dimension(:), pointer :: tmp
 
@@ -121,8 +122,6 @@
 
       call timer( 'pdos', 1)
 
-      gamma = ((nk == 1) .and. (sum(abs(kpoint(:,1))) == 0.0))
-
       orig_Serial = Serial
       orig_ParallelOverK = ParallelOverK
 
 
=== modified file 'Src/post_scf_work.F'
--- Src/post_scf_work.F	2018-06-20 10:34:33 +0000
+++ Src/post_scf_work.F	2018-06-25 11:35:23 +0000
@@ -50,12 +50,9 @@
       character(len=20) :: msg
       type(Pair_Geometry_dSpData2D)  :: pair
       type(Geometry)                :: geom
-      logical                        :: not_using_auxcell
 
       call timer( 'PostSCF', 1 )
 
-      not_using_auxcell = (no_s == no_u)
-!     
 !     If converged, make one last iteration to find forces and stress
 
 !     If we use the minimization routine, the energy-density
@@ -79,12 +76,12 @@
      &                no_l, maxnh, maxnh, no_u,
      &                numh, listhptr, listh, numh, listhptr, listh, 
      &                H, S, qtot, fixspin, qtots, temp, 1.0_dp, -1.0_dp,
-     &                xijo, indxuo, nkpnt, kpoint, kweight,
+     &                xijo, indxuo, gamma_SCF, nkpnt, kpoint, kweight,
      &                eo, qo, Dscf, Escf, ef, efs, Entropy, no_u,
      &                occtol, iscf, neigwanted)
           Ecorrec = 0.0_dp
         else
-          if (not_using_auxcell) then
+          if ( gamma_SCF ) then
             call dminim(.true., .false., iscf, istp, no_l, nspin, no_u,
      &                  maxnh, numh, listhptr, listh, Escf, eta, qtots)
           else
 
=== modified file 'Src/projected_DOS.F'
--- Src/projected_DOS.F	2018-06-12 22:59:08 +0000
+++ Src/projected_DOS.F	2018-06-25 11:35:23 +0000
@@ -121,7 +121,7 @@
          call pdos( no_s, h_spin_dim, spinor_dim, no_l, 
      .        maxnh,
      .        no_u, numh, listhptr, listh, H, S,
-     .        e1, e2, sigma, nhist, xijo, indxuo, 
+     .        e1, e2, sigma, nhist, xijo, indxuo, gamma_PDOS,
      .        nkpnt_pdos, kpoints_pdos, kweight_pdos, eo,
      .        no_u)
       else
@@ -130,7 +130,8 @@
      .        maxnh,
      .        no_u, numh, listhptr, listh, H, S,
      .        e1, e2, sigma, nhist,
-     .        xijo, indxuo, nkpnt, kpoint, kweight, eo,
+     .        xijo, indxuo, gamma_SCF,
+     .        nkpnt, kpoint, kweight, eo,
      .        no_u)
       endif
 
 
=== modified file 'Src/siesta_analysis.F'
--- Src/siesta_analysis.F	2018-06-12 22:59:08 +0000
+++ Src/siesta_analysis.F	2018-06-25 11:35:23 +0000
@@ -15,7 +15,7 @@
 
       subroutine siesta_analysis( relaxd )
       USE band, only: nbk, bk, maxbk, bands
-      USE writewave, only: nwk, wfk, wwave
+      USE writewave, only: nwk, wfk, wwave, gamma_wavefunctions
       USE writewave, only: setup_wfs_list, wfs_filename
       USE m_ksvinit, only: nkpol, kpol, wgthpol
       use m_ksv
@@ -35,7 +35,6 @@
      &                    iphkb, no_u, no_s, iza, iphorb, rmaxo, indxua
       use atomlist, only: qtot
       use fdf
-      use writewave,    only: wwave
       use siesta_cml
       use files,        only : slabel
       use files,        only : filesOut_t   ! derived type for output file names
@@ -274,7 +273,7 @@
         call wwave( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, 
      &              nwk,
      &              numh, listhptr, listh, H, S, Ef, xijo, indxuo,
-     &              nwk, wfk, no_u, occtol )
+     &              gamma_wavefunctions, nwk, wfk, no_u, occtol )
       endif
 
 
@@ -294,7 +293,7 @@
          call wwave( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, 
      .        nkpnt,
      .        numh, listhptr, listh, H, S, Ef, xijo, indxuo,
-     .        nkpnt, kpoint, no_u, occtol)
+     .        gamma_SCF, nkpnt, kpoint, no_u, occtol)
       endif
 
 !     Compute bands
 
=== modified file 'Src/writewave.F'
--- Src/writewave.F	2018-06-12 22:59:08 +0000
+++ Src/writewave.F	2018-06-25 11:35:23 +0000
@@ -386,7 +386,7 @@
       subroutine wwave( no, nspin, maxspn, maxo, maxuo, maxnh, 
      .                  maxk,
      .                  numh, listhptr, listh, H, S, ef, xij, indxuo,
-     .                  nk, kpoint, nuotot, occtol)
+     .                  gamma, nk, kpoint, nuotot, occtol)
 C *********************************************************************
 C Finds wavefunctions at selected k-points.
 C Written by P. Ordejon, June 2003
@@ -415,6 +415,7 @@
 C                               Unit cell orbitals must be the first in
 C                               orbital lists, i.e. indxuo.le.nuo, with
 C                               nuo the number of orbitals in unit cell
+C logical Gamma               : whether only the Gamma-point is sampled
 C integer nk                  : Number of band k points
 C real*8  kpoint(3,maxk)      : k point vectors
 C integer nuotot              : Total number of orbitals in unit cell
@@ -447,6 +448,7 @@
 
       implicit          none
 
+      logical, intent(in) :: Gamma
       integer           maxk, maxnh, maxo, maxuo, nk, no, 
      .                  h_spin_dim, spinor_dim, nspin, maxspn,
      .                  nuotot, indxuo(no), listh(maxnh), numh(*), 
@@ -477,7 +479,6 @@
       data Dnew, Enew, e1, e2, qk, qtot, temp, wk /8*0.d0/
 
       logical :: SaveParallelOverK
-      logical           gamma
 
       h_spin_dim=size(H,dim=2)
 
@@ -492,8 +493,6 @@
 C Start time counter 
       call timer( 'writewave', 1 )
 
-      gamma = ((nk == 1) .and. (sum(abs(kpoint(:,1))) == 0.0))
-      
 C Check parameter maxk 
       if (nk .gt. maxk) then
         if (Node.eq.0) then
 
=== modified file 'Util/COOP/Makefile'
--- Util/COOP/Makefile	2018-06-20 10:34:33 +0000
+++ Util/COOP/Makefile	2018-06-25 11:35:23 +0000
@@ -483,7 +483,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/Denchar/Src/Makefile'
--- Util/Denchar/Src/Makefile	2018-06-20 10:34:33 +0000
+++ Util/Denchar/Src/Makefile	2018-06-25 11:35:23 +0000
@@ -557,7 +557,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/Gen-basis/Makefile'
--- Util/Gen-basis/Makefile	2018-06-20 10:34:33 +0000
+++ Util/Gen-basis/Makefile	2018-06-25 11:35:23 +0000
@@ -547,7 +547,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/Grimme/Makefile'
--- Util/Grimme/Makefile	2018-06-20 10:34:33 +0000
+++ Util/Grimme/Makefile	2018-06-25 11:35:23 +0000
@@ -485,7 +485,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/Helpers/Makefile'
--- Util/Helpers/Makefile	2018-06-20 10:34:33 +0000
+++ Util/Helpers/Makefile	2018-06-25 11:35:23 +0000
@@ -488,7 +488,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/STM/ol-stm/Src/Makefile'
--- Util/STM/ol-stm/Src/Makefile	2018-06-20 10:34:33 +0000
+++ Util/STM/ol-stm/Src/Makefile	2018-06-25 11:35:23 +0000
@@ -570,7 +570,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/SpPivot/Makefile'
--- Util/SpPivot/Makefile	2018-06-20 10:34:33 +0000
+++ Util/SpPivot/Makefile	2018-06-25 11:35:23 +0000
@@ -499,7 +499,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/TS/TBtrans/Makefile'
--- Util/TS/TBtrans/Makefile	2018-06-20 10:34:33 +0000
+++ Util/TS/TBtrans/Makefile	2018-06-25 11:35:23 +0000
@@ -755,7 +755,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/TS/ts2ts/Makefile'
--- Util/TS/ts2ts/Makefile	2018-06-20 10:34:33 +0000
+++ Util/TS/ts2ts/Makefile	2018-06-25 11:35:23 +0000
@@ -492,7 +492,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/TS/tshs2tshs/Makefile'
--- Util/TS/tshs2tshs/Makefile	2018-06-20 10:34:33 +0000
+++ Util/TS/tshs2tshs/Makefile	2018-06-25 11:35:23 +0000
@@ -558,7 +558,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'Util/VCA/Makefile'
--- Util/VCA/Makefile	2018-06-20 10:34:33 +0000
+++ Util/VCA/Makefile	2018-06-25 11:35:23 +0000
@@ -518,7 +518,7 @@
 m_ts_electrode.o: parallel.o precision.o units.o
 m_ts_electype.o: class_OrbitalDistribution.o class_SpData1D.o class_SpData2D.o
 m_ts_electype.o: class_Sparsity.o create_Sparsity_SC.o geom_helper.o
-m_ts_electype.o: intrinsic_missing.o m_geom_box.o m_geom_plane.o
+m_ts_electype.o: intrinsic_missing.o m_char.o m_geom_box.o m_geom_plane.o
 m_ts_electype.o: m_handle_sparse.o m_iodm.o m_os.o m_region.o m_ts_chem_pot.o
 m_ts_electype.o: m_ts_io.o m_ts_io_ctype.o m_ts_iodm.o parallel.o precision.o
 m_ts_electype.o: units.o
 
=== modified file 'version.info'
--- version.info	2018-06-20 10:43:13 +0000
+++ version.info	2018-06-25 11:35:23 +0000
@@ -1,1 +1,1 @@
-siesta-4.1--933--gamma-931
+siesta-4.1--933--gamma-931--nick-1