Merge lp:~albertog/siesta/4.0-images into lp:siesta/4.0

[Alberto Garcia]

(See individual commit messages)

Two independent changes:

1. The fix for the "severe folding" leading to wrong norm of the DM initialized by atomic data.
2. A saner hsparse.

I am not sure it is worth detecting the "severe folding" in hsparse (by checking for sameness of io and jo) for informative purposes.

[Alberto Garcia (albertog)]

For Harris-functional calculations, severe folding is just plain crazy, as there is no time to smooth out the initial differences. We should remove the Harris option from the force_2 test.

In general, we should remove tests with non-physically small cells with gamma point (or at least flag them as showcases for the problems).

[Nick Papior (nickpapior)]

Could we force the auxiliary cell for Harris functional?

1) I think very few people are using it,
2) The overhead is negligeble in this case.

[Nick Papior (nickpapior)]

review: Approve

[lp:~albertog/siesta/4.0-images]

576. By Alberto Garcia on 2018-05-07

Check for diagonal folding. Stop Harris runs if found

577. By Alberto Garcia on 2018-05-07

Add diagonalizer for gamma point with auxiliary supercell

The new routine 'diagk0' is a minor modification of 'diagg' that
can work with an auxiliary supercell while still using real
arithmetic. It is meant to be used when 'force-aux-cell T' is
specified in the input file, for gamma point calculations. (This
option can be set to carry out COOP/COHP analyses, or to avoid
severe folding in test cases.)

The rest of the program logic is kept as it was. It might be worth
rethinking the use of the 'gamma' variable.

578. By Alberto Garcia on 2018-05-07

Use 'diagk0' when appropriate in pdos and writewave

This might save time when computing the COOP-related wavefunctions,
and when carrying out a pDOS calculation with the 'gamma' sampling.

Unmerged revisions

578. By Alberto Garcia on 2018-05-07

Use 'diagk0' when appropriate in pdos and writewave

This might save time when computing the COOP-related wavefunctions,
and when carrying out a pDOS calculation with the 'gamma' sampling.

577. By Alberto Garcia on 2018-05-07

Add diagonalizer for gamma point with auxiliary supercell

The new routine 'diagk0' is a minor modification of 'diagg' that
can work with an auxiliary supercell while still using real
arithmetic. It is meant to be used when 'force-aux-cell T' is
specified in the input file, for gamma point calculations. (This
option can be set to carry out COOP/COHP analyses, or to avoid
severe folding in test cases.)

The rest of the program logic is kept as it was. It might be worth
rethinking the use of the 'gamma' variable.

=== modified file 'Src/Makefile'
--- Src/Makefile	2018-04-11 11:00:02 +0000
+++ Src/Makefile	2018-05-08 09:38:49 +0000
@@ -84,7 +84,7 @@
 	chempot.o coceri.o coxmol.o cross.o compute_norm.o\
 	denmat.o denmatlomem.o detover.o dfscf.o m_diagon.o diagon.o digcel.o \
         fft.o dhscf.o constr.o diagk_file.o \
-	diagg.o diagk.o diagkp.o diag2g.o diag2k.o diagpol.o \
+	diagg.o diagk0.o diagk.o diagkp.o diag2g.o diag2k.o diagpol.o \
 	diagsprl.o dipole.o dismin.o dnaefs.o doping_uniform.o dot.o  \
 	m_efield.o egandd.o ener3.o ener3lomem.o errorf.o extrapolon.o \
 	fixed.o interpolation.o gradient.o gradientlomem.o grdsam.o \
@@ -474,6 +474,8 @@
 diagg.o: writewave.o
 diagk.o: compute_norm.o fermid.o parallel.o parallelsubs.o precision.o sys.o
 diagk.o: writewave.o
+diagk0.o: alloc.o fermid.o parallel.o parallelsubs.o precision.o sys.o
+diagk0.o: writewave.o
 diagk_file.o: fermid.o iowfs_netcdf.o parallel.o parallelsubs.o precision.o
 diagk_file.o: sys.o writewave.o
 diagkp.o: alloc.o fermid.o parallel.o parallelsubs.o precision.o sys.o
@@ -528,8 +530,8 @@
 grdsam.o: alloc.o dhscf.o files.o m_mpi_utils.o m_partial_charges.o parallel.o
 grdsam.o: precision.o siesta_geom.o siesta_options.o sys.o units.o
 hsparse.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o mneighb.o
-hsparse.o: parallel.o parallelsubs.o precision.o siesta_options.o sorting.o
-hsparse.o: sparse_matrices.o sys.o
+hsparse.o: parallel.o parallelsubs.o precision.o sorting.o sparse_matrices.o
+hsparse.o: sys.o
 idiag.o: parallel.o sys.o
 initatom.o: atmparams.o atom.o atom_options.o basis_io.o basis_specs.o
 initatom.o: basis_types.o electrostatic.o old_atmfuncs.o precision.o
@@ -703,6 +705,8 @@
 ofc.o: alloc.o files.o precision.o units.o
 old_atmfuncs.o: alloc.o atmparams.o precision.o sys.o
 on_subs.o: alloc.o onmod.o onmod.o
+one.diagk0.o: compute_norm.o fermid.o parallel.o parallelsubs.o precision.o
+one.diagk0.o: sys.o writewave.o
 onmod.o: precision.o
 optical.o: alloc.o atomlist.o densematrix.o fermid.o files.o parallel.o
 optical.o: parallelsubs.o precision.o sys.o units.o
 
=== added file 'Src/diagk0.F'
--- Src/diagk0.F	1970-01-01 00:00:00 +0000
+++ Src/diagk0.F	2018-05-08 09:38:49 +0000
@@ -0,0 +1,295 @@
+! 
+! Copyright (C) 1996-2016	The SIESTA group
+!  This file is distributed under the terms of the
+!  GNU General Public License: see COPYING in the top directory
+!  or http://www.gnu.org/copyleft/gpl.txt.
+! See Docs/Contributors.txt for a list of contributors.
+!
+      subroutine diagk0( nspin, nuo, maxuo, maxnh, maxnd, 
+     &                  maxo, numh, listhptr, listh, numd, 
+     &                  listdptr, listd, H, S, 
+     &                  getD, getPSI, fixspin, qtot, qs, temp, e1, e2,
+     &                  no_s, indxuo,
+     &                  eo, qo, Dnew, Enew, ef, efs, Entropy,
+     &                  Haux, Saux, psi, nuotot, occtol,
+     &                  iscf, neigwanted )
+C *********************************************************************
+C Subroutine to calculate the eigenvalues and eigenvectors, density
+C and energy-density matrices, and occupation weights of each 
+C eigenvector, for given Hamiltonian and Overlap matrices (including
+C spin polarization). Gamma-point version.
+C Writen by J.Soler, August 1998.
+C **************************** INPUT **********************************
+C integer nspin               : Number of spin components (1 or 2)
+C integer nuo                 : Number of basis orbitals local to node
+C integer maxuo               : Last dimension of xij
+C                               Must be at least max(indxuo)
+C integer maxnh               : Maximum number of orbitals interacting  
+C integer maxnd               : Maximum number of nonzero elements of 
+C                               each row of density matrix
+C integer maxo                : First dimension of eo and qo
+C integer numh(nuo)           : Number of nonzero elements of each row 
+C                               of hamiltonian matrix
+C integer listhptr(nuo)       : Pointer to each row (-1) of the
+C                               hamiltonian matrix
+C integer listh(maxnh)        : Nonzero hamiltonian-matrix element  
+C                               column indexes for each matrix row
+C integer numd(nuo)           : Number of nonzero elements of each row 
+C                               ofdensity matrix
+C integer listdptr(nuo)       : Pointer to each row (-1) of the
+C                               density matrix
+C integer listd(maxnd)        : Nonzero density-matrix element column 
+C                               indexes for each matrix row
+C real*8  H(maxnh,nspin)      : Hamiltonian in sparse form
+C real*8  S(maxnh)            : Overlap in sparse form
+C logical getD                : Find occupations and density matrices?
+C logical getPSI              : Find and print wavefunctions?
+C logical fixspin             : Fix the spin of the system?
+C real*8  qtot                : Number of electrons in unit cell
+C real*8  qs(nspin)           : Number of electrons in unit cell for each
+C                               spin component (if fixed spin option is used)
+C real*8  temp                : Electronic temperature 
+C real*8  e1, e2              : Energy range for density-matrix states
+C                               (to find local density of states)
+C                               Not used if e1 > e2
+C integer no_s                : Number of orbitals in auxiliary supercell
+C integer indxuo(no_s)        : Map from aux supercell to unit cell
+C integer nuotot              : total number of orbitals per unit cell
+C                               over all processors
+C integer iscf                : SCF cycle number
+C real*8  occtol              : Occupancy threshold for DM build
+C integer neigwanted          : Number of eigenvalues wanted
+C *************************** OUTPUT **********************************
+C real*8 eo(maxo,nspn)        : Eigenvalues
+C ******************** OUTPUT (only if getD=.true.) *******************
+C real*8 qo(maxo,nspn)        : Occupations of eigenstates
+C real*8 Dnew(maxnd,nspin)    : Output Density Matrix
+C real*8 Enew(maxnd,nspin)    : Output Energy-Density Matrix
+C real*8 ef                   : Fermi energy
+C real*8 efs(nspin)           : Fermi energy for each spin
+C                               (for fixed spin calculations)
+C real*8 Entropy              : Electronic entropy
+C *************************** AUXILIARY *******************************
+C real*8 Haux(nuotot,nuo)     : Auxiliary space for the hamiltonian matrix
+C real*8 Saux(nuotot,nuo)     : Auxiliary space for the overlap matrix
+C real*8 psi(nuotot,maxuo,nspin) : Auxiliary space for the eigenvectors
+C real*8 aux(nuotot)          : Extra auxiliary space
+C *************************** UNITS ***********************************
+C eo, Enew and ef returned in the units of H.
+C *************************** PARALLEL ********************************
+C The auxiliary arrays are now no longer symmetry and so the order
+C of referencing has been changed in several places to reflect this.
+C *********************************************************************
+C
+C  Modules
+C
+      use precision
+      use sys
+      use parallel,      only : Node, Nodes, BlockSize
+      use parallelsubs,  only : LocalToGlobalOrb
+      use writewave,     only : writew
+      use m_fermid,      only : fermid, fermispin, stepf
+      use alloc
+
+#ifdef MPI
+      use mpi_siesta
+#endif
+
+      implicit none
+
+#ifdef MPI
+      integer 
+     &  MPIerror
+#endif
+
+      integer
+     &  iscf, maxnd, maxnh, maxuo, maxo, nuo, nspin, nuotot,
+     &  neigwanted
+
+      integer 
+     &  listh(maxnh), numh(nuo), listhptr(nuo),
+     &  listd(maxnd), numd(nuo), listdptr(nuo)
+
+      real(dp)
+     &  Dnew(maxnd,nspin), e1, e2, ef, Enew(maxnd,nspin), 
+     &  Entropy, eo(maxo,nspin), H(maxnh,nspin), qo(maxo,nspin), 
+     &  qtot, qs(nspin), S(maxnh), temp, efs(nspin), occtol
+     
+      real(dp)
+     &  Haux(nuotot,nuo), Saux(nuotot,nuo), psi(nuotot,maxuo,nspin) 
+
+      real(dp), dimension(1), parameter :: wk = (/ 1.0_dp /)
+      integer, parameter                :: nk = 1
+
+      integer no_s, indxuo(no_s)
+
+      logical
+     &  getD, getPSI, fixspin
+
+      external rdiag
+
+C  Internal variables .............................................
+      integer           ie, io, iio, ispin, ix, j, jo, BNode, iie, ind,
+     &                  ierror, nd
+      integer           jo_s 
+      real(dp)          ee, eei, qe, qei, rt, t, k(3)
+
+      integer           :: maxnuo, mm
+      integer, pointer  :: nuo_LOC(:)
+      real(dp), pointer :: psi_tmp(:), paux(:)
+
+C Solve eigenvalue problem
+      do ispin = 1,nspin
+        call timer( 'r-eigvec', 1 )
+        call timer( 'r-buildHS', 1 )
+        do io = 1,nuo
+          do jo = 1,nuotot
+            Saux(jo,io) = 0.0d0
+            Haux(jo,io) = 0.0d0
+          enddo
+        enddo
+        do io = 1,nuo
+          do j = 1,numh(io)
+            ind = listhptr(io) + j
+            jo_s = listh(ind)
+            jo = indxuo(jo_s)  ! In unit cell
+            Saux(jo,io) = Saux(jo,io) + S(ind)
+            Haux(jo,io) = Haux(jo,io) + H(ind,ispin)
+          enddo
+        enddo
+        call timer( 'r-buildHS', 2 )
+        call rdiag( Haux, Saux, nuotot, nuo, nuotot, eo(1,ispin), 
+     &              psi(1,1,ispin), neigwanted, iscf, ierror )
+        call timer( 'r-eigvec', 2 )
+
+
+C Check error flag and take appropriate action
+        if (ierror.gt.0) then
+          call die('Terminating due to failed diagonalisation')
+        elseif (ierror.lt.0) then
+C Repeat diagonalisation with increased memory to handle clustering
+          do io = 1,nuo
+            do jo = 1,nuotot
+              Saux(jo,io) = 0.0d0
+              Haux(jo,io) = 0.0d0
+            enddo
+          enddo
+          do io = 1,nuo
+            do j = 1,numh(io)
+              ind = listhptr(io) + j
+              jo_s = listh(ind)
+              jo = indxuo(jo_s) ! In unit cell
+              Saux(jo,io) = Saux(jo,io) + S(ind)
+              Haux(jo,io) = Haux(jo,io) + H(ind,ispin)
+            enddo
+          enddo
+          call rdiag( Haux, Saux, nuotot, nuo, nuotot, eo(1,ispin), 
+     &                psi(1,1,ispin), nuotot, iscf, ierror )
+        endif
+
+        if (getPSI) then
+          do ix = 1,3
+            k(ix) = 0.0d0
+          enddo
+          call writew(nuotot,nuo,1,k,ispin,
+     &                eo(1,ispin),psi(1,1,ispin),.true.)
+        endif
+
+      enddo
+C Check if we are done ................................................
+      if (.not.getD) return
+
+C Find new Fermi energy and occupation weights ........................
+      if (fixspin) then
+        call fermispin( nspin, nspin, nk, wk, maxo, neigwanted, 
+     &                  eo, temp, qs, qo, efs, Entropy )
+      else
+        call fermid( nspin, nspin, nk, wk, maxo, neigwanted, eo, 
+     &               temp, qtot, qo, ef, Entropy )
+      endif
+
+C Find weights for local density of states ............................
+      if (e1 .lt. e2) then
+*       e1 = e1 - ef
+*       e2 = e2 - ef
+        t = max( temp, 1.d-6 )
+        rt = 1.0d0/t
+        do ispin = 1,nspin
+          do io = 1,nuotot
+            qo(io,ispin) = ( stepf((eo(io,ispin)-e2)*rt) -
+     &                       stepf((eo(io,ispin)-e1)*rt)) * 
+     &                       2.0d0/dble(nspin)
+          enddo
+        enddo
+      endif
+      
+      if (nuo.gt.0) then
+C New density and energy-density matrices of unit-cell orbitals .......
+        nd = listdptr(nuo) + numd(nuo)
+        Dnew(1:nd,1:nspin) = 0.0d0
+        Enew(1:nd,1:nspin) = 0.0d0
+      endif
+
+      call timer( 'r-buildD', 1 )
+C Global operation to form new density matrix
+      nullify( nuo_LOC )
+      call re_alloc( nuo_LOC, 0, Nodes-1, 'nuo_LOC', 'diagk0' )
+#ifdef MPI
+      call MPI_Allgather( nuo, 1, MPI_INTEGER, nuo_LOC, 1,
+     &                    MPI_INTEGER, MPI_Comm_World, MPIerror )
+#else
+      nuo_LOC(0) = nuo
+#endif
+      maxnuo = 0
+      do BNode=0, Nodes-1
+        maxnuo = max(nuo_LOC(BNode),maxnuo)
+      enddo
+      nullify( PSI_TMP )
+      call re_alloc( PSI_TMP, 1, nuotot*maxnuo*nspin, 'PSI_TMP',
+     &               'diagk0' )
+
+      do Bnode=0, Nodes-1
+        if (BNode.eq.Node) then
+          mm = 0
+          do ispin= 1, nspin
+            do io= 1, nuo
+              PSI_TMP(mm+1:mm+nuotot) = PSI(1:nuotot,io, ispin)
+              mm = mm + nuotot
+            enddo
+          enddo
+        endif
+#ifdef MPI
+        call MPI_Bcast( PSI_TMP, nuotot*nuo_LOC(Bnode)*nspin,
+     &                  MPI_double_precision, BNode, MPI_Comm_World,
+     &                  MPIerror )
+#endif
+        do ispin = 1,nspin
+          do io = 1,nuo
+            call LocalToGlobalOrb(io,Node,Nodes,iio)
+            mm   = (ispin-1)*nuo_LOC(Bnode)*nuotot
+            do iie = 1,nuo_LOC(Bnode)
+              call  LocalToGlobalOrb(iie,BNode,Nodes,ie)
+              qe = qo(ie,ispin)
+              if (abs(qe).gt.occtol) then
+                paux => PSI_TMP(mm+1:mm+nuotot)
+                ee  = qe*eo(ie,ispin)
+                qei = qe*paux(iio)
+                eei = ee*paux(iio)
+                do j = 1,numd(io)
+                  ind = listdptr(io) + j
+                  jo_s  = listd(ind)
+                  jo = indxuo(jo_s) ! In unit cell
+                  Dnew(ind,ispin) = Dnew(ind,ispin) + qei*paux(jo)
+                  Enew(ind,ispin) = Enew(ind,ispin) + eei*paux(jo)
+                enddo
+              endif
+              mm = mm + nuotot
+            enddo
+          enddo
+        enddo
+      enddo
+
+      call timer( 'r-buildD', 2 )
+      call de_alloc( nuo_LOC, 'nuo_LOC', 'diagk0' )
+      call de_alloc( PSI_TMP, 'PSI_TMP', 'diagk0' )
+      end
 
=== modified file 'Src/diagon.F'
--- Src/diagon.F	2017-02-21 07:36:32 +0000
+++ Src/diagon.F	2018-05-08 09:38:49 +0000
@@ -126,7 +126,7 @@
      .  fixspin, gamma, getD, getPSI
 
       external
-     .  diagg, diagk, diag2g, diag2k, readsp, diagsprl
+     .  diagg, diagk0, diagk, diag2g, diag2k, readsp, diagsprl
 #ifdef CDF
       external diagk_file
 #endif
@@ -255,6 +255,7 @@
      .              eo, qo, Dnew, Enew, ef, efs, Entropy,
      .              Haux, Saux, psi, Haux, Saux, aux, 
      .              nuotot, occtol, iscf )
+          call output_kpoint_occupation()
         else
 #endif
 #ifndef CDF
@@ -279,8 +280,22 @@
      .              eo, qo, Dnew, Enew, ef, efs, Entropy,
      .              Haux, Saux, psi, Haux, Saux, aux, 
      .              nuotot, occtol, iscf, neigwanted )
+              call output_kpoint_occupation()
 #endif
            else
+            if ( (nk==1) .and.
+     $             sum(kpoint(:,1)) == 0.0_dp) then
+!
+!             Use special version  (This is still gamma)
+!
+              call diagk0( nspin, nuo, maxuo, maxnh, maxnd, maxo,
+     .              numh, listhptr, listh, numd, listdptr, listd, 
+     .              H, S, getD, getPSI, fixspin, qtot, qs, temp, e1, e2,
+     .              no, indxuo,
+     .              eo, qo, Dnew, Enew, ef, efs, Entropy,
+     .              Haux, Saux, psi, nuotot, occtol, iscf,
+     .              neigwanted )
+            else     
               call diagk( nspin, nuo, no, maxspn, maxnh, maxnd, 
      .              maxo, numh, listhptr, listh, numd, listdptr,
      .              listd, H, S, getD, getPSI, fixspin, qtot, qs, temp, 
@@ -288,8 +303,9 @@
      .              eo, qo, Dnew, Enew, ef, efs, Entropy,
      .              Haux, Saux, psi, Haux, Saux, aux, 
      .              nuotot, occtol, iscf )
+              call output_kpoint_occupation()
+            endif
            endif
-           call output_kpoint_occupation()
 #ifdef MPI
         endif
 #endif
 
=== modified file 'Src/hsparse.F'
--- Src/hsparse.F	2016-01-25 16:00:16 +0000
+++ Src/hsparse.F	2018-05-08 09:38:49 +0000
@@ -48,7 +48,8 @@
       subroutine hsparse( negl, cell, nsc, na, isa, xa,
      &                    lasto, lastkb, iphorb, iphkb,     
      &                    nlhmax, gamma,
-     $                    set_xijo, folding)
+     $                    set_xijo, folding,
+     $                    diagonal_folding, debug_folding)
 
       use precision
       use parallel,        only : Node, Nodes
@@ -61,7 +62,7 @@
       use neighbour,       only : mneighb, reset_neighbour_arrays
       use sys,             only : die
       use alloc,           only : re_alloc, de_alloc
-      use atomlist,        only : no_l
+      use atomlist,        only : no_l, indxuo
       use sparse_matrices, only : listhptr, numh, listh, xijo
       implicit none
 
@@ -76,6 +77,8 @@
       logical, intent(in)           :: gamma
       logical, intent(in)     :: set_xijo
       logical, intent(out)    :: folding
+      logical, intent(out)    :: diagonal_folding
+      logical, intent(in), optional    :: debug_folding
 
       external               timer
 
@@ -99,8 +102,6 @@
      &  rci, rcj, rck, rij, rik, rjk,
      &  rmax, rmaxkb, rmaxo
 
-      logical, save :: warn1 = .false.
-
       real(dp), dimension(:), pointer :: rckb
       logical, dimension(:),  pointer :: conect
       integer, dimension(:),  pointer :: index
@@ -108,14 +109,18 @@
       integer, dimension(:),  pointer :: listhtmp
 
       logical :: connected
-
+      logical :: debug
+      integer :: nprints   ! for debugging
 C -------------------------------------
-#ifdef DEBUG
-      call write_debug( '    PRE hsparse' )
-#endif
+
 C     Start time counter
       call timer( 'hsparse', 1 )
 
+      debug = .false.
+      if (present(debug_folding)) then
+         debug = debug_folding
+      endif
+      
 C     Check size of internal arrays
       ncells = nsc(1) * nsc(2) * nsc(3)
       nua = na / ncells
@@ -192,6 +197,7 @@
       enddo
 
       folding = .false.
+      diagonal_folding = .false.
 
 C------------------------------------C
 C     Find number of non-zeros in H  C
@@ -324,6 +330,7 @@
 C------------------------------------C
 C     Find full H sparsity pattern   C
 C------------------------------------C
+      nprints = 0
 C     Loop on atoms in unit cell
       do ia = 1,nua
 C       Find neighbour atoms within maximum range
@@ -396,62 +403,53 @@
                     endif
 
                    if (connected) then
-                     if (conect(jo)) then
-                        folding = .true.
-
-                       ! This test is now deferred to be able
-                       ! to catch multiple images while avoiding
-                       ! false positives (i.e., we test first
-                       ! whether there is indeed a connection).
-
-                      if (.true.) then
-                         ! If already connected and using supercell, 
-                         ! the latter might not be big enough...
-                         ! We warn the user and keep the first instance
+                      if (conect(jo)) then
+                         folding = .true.
+                         if (io == indxuo(jo)) then
+                            diagonal_folding = .true.
+                         endif
+
+                         ! If already connected and using supercell, the
+                         ! latter might not be big enough...  We defer
+                         ! the decision on whether to kill the program
+                         ! to the caller.  Here keep the first instance
                          ! of xij (same behavior as the old xijorb, as
                          ! earlier jnats are closer)
-                   ! Warn also if Gamma-point calculation, just
-                         ! in case
-                         if (.not.warn1) then
-                            if (Node.eq.0) then
-                              if (gamma) then
-                                 call check_cohp(io,jo)
-                              else
-                                write(6,'(/,a,2i6,a,/,a)')
-     .                           'WARNING: orbital pair ',io,jo,
+                         if (debug) then
+                            ! This might not be pretty in
+                            ! parallel. Better to build a list of a few
+                            ! cases, and pass it to the caller.
+                            if (nprints <= 20) then
+                               print "(a,2i6,a)",
+     .                           'WARNING: orbital pair ',io,indxuo(jo),
      .                           ' is multiply connected'
-                              endif
+                               nprints = nprints + 1
                             endif
-                            warn1 = .true.
-                         endif
-                      endif
-
-                   else
-
-                      conect(jo) = .true.
-                      numh(iio) = numh(iio) + 1
-                      ind = listhptr(iio)+numh(iio)
-                      listh(ind) = jo
-                      if (set_xijo) then
-                         xijo(1:3,ind) = xij(1:3,jnat)
+                         endif
+
+                      else
+
+                         conect(jo) = .true.
+                         numh(iio) = numh(iio) + 1
+                         ind = listhptr(iio)+numh(iio)
+                         listh(ind) = jo
+                         if (set_xijo) then
+                            xijo(1:3,ind) = xij(1:3,jnat)
+                         endif
                       endif
                    endif
-                endif
+                enddo
              enddo
-          enddo
 
 C           Restore conect array for next orbital io
-            do j = 1,numh(iio)
-               jo = listh(listhptr(iio)+j)
-               conect(jo) = .false.
-            enddo
-          endif     ! iio > 0
-        enddo       ! io
-      enddo         ! ia
+             do j = 1,numh(iio)
+                jo = listh(listhptr(iio)+j)
+                conect(jo) = .false.
+             enddo
+          endif                 ! iio > 0
+       enddo                    ! io
+      enddo                     ! ia
 
-!!        print "(a5,i3,a40,3i8)",
-!!     $         "Node: ", Node, "in hsparse nuo, nuotot, nlhmax: ",
-!!     $         nuo, nuotot, nlhmax
 
 C     Initialize listsc
       call LISTSC_INIT( nsc, nuotot )
@@ -466,28 +464,7 @@
       call de_alloc( rckb, 'rckb', 'hsparse' )
 
       call timer( 'hsparse', 2 )
-#ifdef DEBUG
-      call write_debug( '    POS hsparse' )
-#endif
+
       end subroutine hsparse
-      
-      subroutine check_cohp(io,jo)
-      use siesta_options, only: write_coop
-
-      integer, intent(in) :: io, jo
- 
-      if (write_coop) then
-         write(6,'(/,a,2i6,a,/,a)')
-     .     'NOTE: orbital pair ',io,jo,
-     .     ' (at least) is multiply connected.',
-     .     'NOTE: Your COOP/COHP analysis might ' //
-     $     'be affected by folding.'
-         write(0,'(/,a,2i6,a,/,a)')
-     .     'NOTE: orbital pair ',io,jo,
-     .     ' (at least) is multiply connected.',
-     .     'NOTE: Your COOP/COHP analysis might ' //
-     $     'be affected by folding.'
-      endif
-      end subroutine check_cohp
 
       end module m_hsparse
 
=== modified file 'Src/new_dm.F'
--- Src/new_dm.F	2018-04-06 12:32:04 +0000
+++ Src/new_dm.F	2018-05-08 09:38:49 +0000
@@ -485,6 +485,7 @@
       use mpi_siesta
 #endif
       use units, only : pi
+      use sparse_matrices, only: S
 
       implicit          none
 
@@ -655,7 +656,7 @@
         endif
 
 
-      else
+      else  ! Initialize with neutral atoms
 
 C See whether specific initial spins are given in a DM.InitSpin block
 C and read them in a loop on atoms where lines are read and parsed
@@ -903,6 +904,24 @@
 
         endif
 
+      ! We have initialized with atomic information. Correct in case we
+      ! are using such a small cell that there are direct interactions
+      ! of orbitals with their own images, and we insist on using the
+      ! Gamma-point only. Otherwise S(diagonal) is always 1.0 and the
+      ! simple atomic-orbital superposition works as intended.
+        
+
+        do io = 1, no_l
+            call LocalToGlobalOrb(io,Node,Nodes,iio)
+            do in = 1,numh(io)
+               ind = listhptr(io)+in
+               jo = listh(ind)
+               if (iio .eq. jo) then  ! diagonal element
+                  Dscf(ind,:) = Dscf(ind,:) / S(ind)
+               endif
+            enddo
+         enddo
+
       endif
 
       end subroutine initdm
 
=== modified file 'Src/pdos.F'
--- Src/pdos.F	2018-05-07 07:36:34 +0000
+++ Src/pdos.F	2018-05-08 09:38:49 +0000
@@ -180,11 +180,19 @@
      .                XIJ, INDXUO, NK, KPOINT, WK, EO,
      .                HAUX, SAUX, PSI, DTOT, DPR, NO_U )
         else
-          call pdosk( NSPIN, NUO, NO, MAXSPN, MAXNH,
+          if ((nk==1) .and. (sum(kpoint(:,1)) == 0.0d0)) then
+           ! pdosg already setup for auxiliary supercell!!
+           call pdosg( NSPIN, NUO, NO, MAXSPN, MAXNH,
+     .              MAXO, NUMH, LISTHPTR, LISTH, H, S,
+     .              E1, E2, NHIST, SIGMA, INDXUO, EO,
+     .              HAUX, SAUX, PSI, DTOT, DPR, NO_U )
+          else
+           call pdosk( NSPIN, NUO, NO, MAXSPN, MAXNH,
      .                MAXO, NUMH, LISTHPTR, LISTH, H, S,
      .                E1, E2, NHIST, SIGMA, 
      .                XIJ, INDXUO, NK, KPOINT, WK, EO,
      .                HAUX, SAUX, PSI, DTOT, DPR, NO_U )
+          endif
         endif
       endif
 
 
=== modified file 'Src/state_init.F'
--- Src/state_init.F	2018-04-11 09:23:27 +0000
+++ Src/state_init.F	2018-05-08 09:38:49 +0000
@@ -48,7 +48,7 @@
       use m_mpi_utils,       only: globalize_or
       use m_mpi_utils,       only: globalize_sum
       use domain_decom,      only: domainDecom, use_dd, use_dd_perm
-
+      use fdf,               only: fdf_get
 #ifdef TRANSIESTA
       use m_ts_options, only : onlyS
       use sys, only : bye
@@ -69,6 +69,7 @@
       integer            :: i, ix, iadispl, ixdispl
       logical            :: auxchanged   ! Auxiliary supercell changed?
       logical            :: folding, folding1
+      logical            :: diag_folding, diag_folding1
       logical            :: foundxv  ! dummy for call to ioxv
       real(dp)           :: dummy_qspin(8)
       external           ::  madelung, timer
@@ -255,7 +256,39 @@
 !           analyses.
       call hsparse( negl, scell, nsc, na_s, isa, xa, lasto,
      &              lastkb, iphorb, iphKB, maxnh, gamma,
-     $                 set_xijo=.true., folding=folding1)
+     $              set_xijo=.true., folding=folding1,
+     $              diagonal_folding=diag_folding1,
+     $              debug_folding=fdf_get('debug-folding',.false.))
+!
+      call globalize_or(diag_folding1,diag_folding)
+      if (diag_folding .and. gamma) then
+         if (IOnode) then
+            write(6,"(a)") "Gamma-point calculation " //
+     $         "with interaction between periodic images!"
+            write(6,"(a)") "Some features might not work optimally"
+         endif
+         if (harrisfun) call die("Harris functional run needs " //
+     $                              "'force-aux-cell T'")
+      endif
+      !
+      call globalize_or(folding1,folding)
+      if (folding) then
+         if (gamma) then
+            if (IOnode) then
+               write(6,"(a)") "Gamma-point calculation " //
+     $                        "with multiply-connected orbital pairs"
+               write(6,"(a)") "Folding of H and S implicitly performed"
+               call check_cohp()
+            endif
+         else
+            write(6,"(a,/,a)") "Non Gamma-point calculation " //
+     $           "with multiply-connected orbital pairs " //
+     $           "in auxiliary supercell.",
+     $           "Possible internal error. " //
+     $           "Use 'debug-folding T' to debug."
+            call die("Inadequate auxiliary supercell")
+         endif
+      endif
 !
       call globalize_sum(maxnh,nnz)
       if (cml_p) then
@@ -269,12 +302,6 @@
          call cmlEndPropertyList(mainXML)
       endif
       !
-      call globalize_or(folding1,folding)
-      if (folding) then
-         if (IOnode) then
-            print *, "Folding of H and S is implicitly performed"
-         endif
-      endif
       !
       ! If using domain decomposition, redistribute orbitals
       ! for this geometry, based on the hsparse info. 
@@ -387,4 +414,22 @@
 
 !--------------------------------------------------------------------------- END
       END subroutine state_init
+
+      subroutine check_cohp()
+      use siesta_options, only: write_coop
+
+      if (write_coop) then
+         write(6,'(/,a,a,/,a)')
+     .     'NOTE: There are multiply-connected orbitals.',
+     .     'NOTE: Your COOP/COHP analysis might ' //
+     $     'be affected by folding.',
+     .     'NOTE: Use "force-aux-cell T" or k-point sampling'
+         write(0,'(/,a,a,/,a)')
+     .     'NOTE: There are multiply-connected orbitals.',
+     .     'NOTE: Your COOP/COHP analysis might ' //
+     $     'be affected by folding.',
+     .     'NOTE: Use "force-aux-cell T" or k-point sampling'
+      endif
+      end subroutine check_cohp
+
       END module m_state_init
 
=== modified file 'Src/writewave.F'
--- Src/writewave.F	2017-10-04 09:48:27 +0000
+++ Src/writewave.F	2018-05-08 09:38:49 +0000
@@ -471,6 +471,7 @@
       data Dnew, Enew, e1, e2, qk, qtot, temp, wk /8*0.d0/
 
       logical  SaveParallelOverK
+      external :: diagg, diagk0, diagk
 
 C Get local number of orbitals
 #ifdef MPI
@@ -600,8 +601,17 @@
          ParallelOverK = .false.
          ResetFirstCall = .true.
 
-            
-         call diagk( nspin, nuo, no, maxspn, maxnh, maxnh, 
+         if ( (nk==1) .and. (sum(kpoint(:,1)) == 0.0d0))  then
+           call diagk0( nspin, nuo, maxuo, maxnh, maxnh, 
+     .               maxo, numh, listhptr, listh, numh, listhptr, 
+     .               listh, H, S, getD, getPSI,
+     .               fixspin, qtot, qs, temp,
+     .               e1, e2,
+     $               no, indxuo,
+     $               ek, qk, Dnew, Enew, ef, efs, Entropy,
+     .               Haux, Saux, psi, nuotot, occtol, 1, nuotot )
+         else
+           call diagk( nspin, nuo, no, maxspn, maxnh, maxnh, 
      .               maxo, numh, listhptr, listh, numh, listhptr, 
      .               listh, H, S, getD, getPSI,
      .               fixspin, qtot, qs, temp,
@@ -609,6 +619,7 @@
      .               ek, qk, Dnew, Enew, ef, efs, Entropy,
      .               Haux, Saux, psi, Haux, Saux, aux, nuotot,
      .               occtol, 1 )
+         endif
          ParallelOverK = SaveParallelOverK
          ResetFirstCall = .false.
 
@@ -663,7 +674,7 @@
       real(SP), dimension(:,:), allocatable :: aux   !! NOTE SP
 
       external io_assign, io_close
-
+      
 C ...................
 
 C Allocate auxiliary arrays
 
=== modified file 'version.info'
--- version.info	2018-05-07 12:07:29 +0000
+++ version.info	2018-05-08 09:38:49 +0000
@@ -1,1 +1,10 @@
+<<<<<<< TREE
 siesta-4.0--575
+=======
+siesta-4.0--573--folding-578
+
+
+
+
+
+>>>>>>> MERGE-SOURCE