1!
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8      subroutine diagg( h_spin_dim, nuo, maxuo, maxnh, maxnd,
9     &                  maxo, numh, listhptr, listh, numd,
10     &                  listdptr, listd, H, S,
11     &                  getD, getPSI, fixspin, qtot, qs, temp, e1, e2,
12     &                  eo, qo, Dnew, Enew, ef, efs, Entropy,
13     &                  Haux, Saux, psi, nuotot, occtol,
14     &                  iscf, neigwanted )
15C *********************************************************************
16C Subroutine to calculate the eigenvalues and eigenvectors, density
17C and energy-density matrices, and occupation weights of each
18C eigenvector, for given Hamiltonian and Overlap matrices (including
19C spin polarization). Gamma-point version.
20C Writen by J.Soler, August 1998.
21C **************************** INPUT **********************************
22C integer h_spin_dim    : Number of spin components of H and D
23C integer spinor_dim          : # of spin components of eo and qo, qs, efs
24C integer e_spin_dim          : Number of spin components of E_dm
25C integer nuo                 : Number of basis orbitals local to node
26C integer maxuo               : Last dimension of xij
27C                               Must be at least max(indxuo)
28C integer maxnh               : Maximum number of orbitals interacting
29C integer maxnd               : Maximum number of nonzero elements of
30C                               each row of density matrix
31C integer maxo                : First dimension of eo and qo
32C integer numh(nuo)           : Number of nonzero elements of each row
33C                               of hamiltonian matrix
34C integer listhptr(nuo)       : Pointer to each row (-1) of the
35C                               hamiltonian matrix
36C integer listh(maxnh)        : Nonzero hamiltonian-matrix element
37C                               column indexes for each matrix row
38C integer numd(nuo)           : Number of nonzero elements of each row
39C                               ofdensity matrix
40C integer listdptr(nuo)       : Pointer to each row (-1) of the
41C                               density matrix
42C integer listd(maxnd)        : Nonzero density-matrix element column
43C                               indexes for each matrix row
44C real*8  H(maxnh,h_spin_dim) : Hamiltonian in sparse form
45C real*8  S(maxnh)            : Overlap in sparse form
46C logical getD                : Find occupations and density matrices?
47C logical getPSI              : Find and print wavefunctions?
48C logical fixspin             : Fix the spin of the system?
49C real*8  qtot                : Number of electrons in unit cell
50C real*8  qs(spinor_dim)      : Number of electrons in unit cell for each
51C                               spin component (if fixed spin option is used)
52C real*8  temp                : Electronic temperature
53C real*8  e1, e2              : Energy range for density-matrix states
54C                               (to find local density of states)
55C                               Not used if e1 > e2
56C integer nuotot              : total number of orbitals per unit cell
57C                               over all processors
58C integer iscf                : SCF cycle number
59C real*8  occtol              : Occupancy threshold for DM build
60C integer neigwanted          : Number of eigenvalues wanted
61C *************************** OUTPUT **********************************
62C real*8 eo(maxo,spinor_dim)  : Eigenvalues
63C ******************** OUTPUT (only if getD=.true.) *******************
64C real*8 qo(maxo,spinor_dim)   : Occupations of eigenstates
65C real*8 Dnew(maxnd,h_spin_dim): Output Density Matrix
66C real*8 Enew(maxnd,e_spin_dim): Output Energy-Density Matrix
67C real*8 ef                    : Fermi energy
68C real*8 efs(spinor_dim)       : Fermi energy for each spin
69C                                 (for fixed spin calculations)
70C real*8 Entropy               : Electronic entropy
71C *************************** AUXILIARY *******************************
72C real*8 Haux(nuotot,nuo)     : Auxiliary space for the hamiltonian matrix
73C real*8 Saux(nuotot,nuo)     : Auxiliary space for the overlap matrix
74C real*8 psi(nuotot,maxuo,spinor_dim) : Auxiliary space for the eigenvectors
75C real*8 aux(nuotot)          : Extra auxiliary space
76C *************************** UNITS ***********************************
77C eo, Enew and ef returned in the units of H.
78C *************************** PARALLEL ********************************
79C The auxiliary arrays are now no longer symmetry and so the order
80C of referencing has been changed in several places to reflect this.
81C *********************************************************************
82C
83C  Modules
84C
85      use precision
86      use sys
87      use parallel,      only : Node, Nodes, BlockSize
88      use parallelsubs,  only : LocalToGlobalOrb
89      use writewave,     only : writew
90      use m_fermid,      only : fermid, fermispin, stepf
91      use m_spin,        only : spinor_dim, e_spin_dim
92      use alloc
93      use intrinsic_missing, only: MODP
94      use iso_c_binding, only : c_f_pointer, c_loc
95#ifdef MPI
96      use mpi_siesta
97#endif
98
99      implicit none
100
101#ifdef MPI
102      integer
103     &  MPIerror
104#endif
105
106      integer
107     &  iscf, maxnd, maxnh, maxuo, maxo, nuo, nuotot,
108     &  neigwanted, h_spin_dim
109
110      integer
111     &  listh(maxnh), numh(nuo), listhptr(nuo),
112     &  listd(maxnd), numd(nuo), listdptr(nuo)
113
114      real(dp)
115     &  Dnew(maxnd,h_spin_dim), e1, e2, ef, Enew(maxnd,e_spin_dim),
116     &  Entropy, eo(maxo,spinor_dim), H(maxnh,h_spin_dim),
117     &  qo(maxo,spinor_dim), qtot, qs(spinor_dim), S(maxnh), temp,
118     &  efs(spinor_dim), occtol
119
120      real(dp) Haux(nuotot,nuo), Saux(nuotot,nuo)
121      real(dp), target :: psi(nuotot,maxuo,spinor_dim)
122      real(dp), pointer :: psi_real_1d(:)
123      real(dp), dimension(1), parameter :: wk = (/ 1.0_dp /)
124      integer, parameter                :: nk = 1
125
126      logical
127     .  getD, getPSI, fixspin
128
129      external rdiag
130
131C  Internal variables .............................................
132      integer           ie, io, iio, ispin, ix, j, jo, BNode, iie, ind,
133     &                  ierror, nd
134      real(dp)          ee, eei, qe, qei, rt, t, k(3)
135
136      integer           :: maxnuo, mm
137      integer, pointer  :: nuo_LOC(:)
138      real(dp), pointer :: psi_tmp(:), paux(:)
139
140#ifdef DEBUG
141      call write_debug( '    PRE diagg' )
142#endif
143      if (h_spin_dim /= spinor_dim) then
144         call die("Spin size mismatch in diagg")
145      endif
146
147C Solve eigenvalue problem
148
149      do ispin = 1,spinor_dim
150
151        call timer( 'r-eigvec', 1 )
152        call timer( 'r-buildHS', 1 )
153!$OMP parallel do default(shared), private(io,j,ind,jo)
154        do io = 1,nuo
155          Saux(:,io) = 0._dp
156          Haux(:,io) = 0._dp
157          do j = 1,numh(io)
158            ind = listhptr(io) + j
159            jo = listh(ind)
160            jo = MODP(jo,nuotot)   ! To allow auxiliary supercells
161            Saux(jo,io) = Saux(jo,io) + S(ind)
162            Haux(jo,io) = Haux(jo,io) + H(ind,ispin)
163          enddo
164        enddo
165!$OMP end parallel do
166        call timer( 'r-buildHS', 2 )
167        call rdiag( Haux, Saux, nuotot, nuo, nuotot,
168     $       eo(1,ispin),
169     .       psi(1,1,ispin), neigwanted, iscf, ierror,
170     .       BlockSize)
171        call timer( 'r-eigvec', 2 )
172
173
174C Check error flag and take appropriate action
175        if (ierror.gt.0) then
176          call die('Terminating due to failed diagonalisation')
177        elseif (ierror.lt.0) then
178C Repeat diagonalisation with increased memory to handle clustering
179!$OMP parallel do default(shared), private(io,j,ind,jo)
180          do io = 1,nuo
181            Saux(:,io) = 0._dp
182            Haux(:,io) = 0._dp
183            do j = 1,numh(io)
184              ind = listhptr(io) + j
185              jo = listh(ind)
186              jo = MODP(jo,nuotot) ! To allow auxiliary supercells
187              Saux(jo,io) = Saux(jo,io) + S(ind)
188              Haux(jo,io) = Haux(jo,io) + H(ind,ispin)
189            enddo
190          enddo
191!$OMP end parallel do
192          call rdiag( Haux, Saux, nuotot, nuo, nuotot, eo(1,ispin),
193     &         psi(1,1,ispin), neigwanted, iscf, ierror,
194     .         BlockSize )
195        endif
196
197        if (getPSI) then
198          ! We do not have info about occupation
199          ! In the current implementation, neigneeded
200          ! is ALWAYS nuotot.
201          ! This is because writewave expects the full spectrum
202          do ix = 1,3
203            k(ix) = 0.0d0
204          enddo
205          call c_f_pointer(c_loc(psi(:,:,ispin)),psi_real_1d,
206     $                     [size(psi,1)*size(psi,2)])
207          call writew(nuotot,nuo,1,k,ispin,
208     &                eo(1,ispin),psi_real_1d,gamma=.true.,
209     $                non_coll=.false.,blocksize=BlockSize)
210        endif
211
212      enddo
213C Check if we are done ................................................
214      if (.not.getD) then
215#ifdef DEBUG
216         call write_debug( '    POS diagg' )
217#endif
218         return
219      end if
220
221C Find new Fermi energy and occupation weights ........................
222      if (fixspin) then
223        call fermispin( spinor_dim, spinor_dim, nk, wk, maxo,
224     &                  neigwanted, eo, temp, qs, qo, efs, Entropy )
225      else
226        call fermid(spinor_dim, spinor_dim, nk, wk, maxo, neigwanted,
227     &               eo, temp, qtot, qo, ef, Entropy )
228      endif
229
230
231C Find weights for local density of states ............................
232      if (e1 .lt. e2) then
233*       e1 = e1 - ef
234*       e2 = e2 - ef
235        t = max( temp, 1.d-6 )
236        rt = 1.0d0/t
237!$OMP parallel do default(shared), private(ispin,io)
238        do ispin = 1,spinor_dim
239          do io = 1,nuotot
240            qo(io,ispin) = ( stepf((eo(io,ispin)-e2)*rt) -
241     &                       stepf((eo(io,ispin)-e1)*rt)) *
242     &                       2.0d0/dble(spinor_dim)
243          enddo
244        enddo
245!$OMP end parallel do
246      endif
247
248      if (nuo.gt.0) then
249C New density and energy-density matrices of unit-cell orbitals .......
250        nd = listdptr(nuo) + numd(nuo)
251        Dnew(1:nd,1:h_spin_dim) = 0.0d0
252        Enew(1:nd,1:e_spin_dim) = 0.0d0
253      endif
254
255      call timer( 'r-buildD', 1 )
256C Global operation to form new density matrix
257      nullify( nuo_LOC )
258      call re_alloc( nuo_LOC, 0, Nodes-1, 'nuo_LOC', 'diagg' )
259#ifdef MPI
260      call MPI_Allgather( nuo, 1, MPI_INTEGER, nuo_LOC, 1,
261     &                    MPI_INTEGER, MPI_Comm_World, MPIerror )
262#else
263      nuo_LOC(0) = nuo
264#endif
265      maxnuo = 0
266      do BNode=0, Nodes-1
267        maxnuo = max(nuo_LOC(BNode),maxnuo)
268      enddo
269      nullify( PSI_TMP )
270      call re_alloc( PSI_TMP, 1, nuotot*maxnuo*spinor_dim, 'PSI_TMP',
271     &               'diagg' )
272
273      do Bnode=0, Nodes-1
274        if (BNode.eq.Node) then
275          mm = 0
276          do ispin= 1, spinor_dim
277            do io= 1, nuo
278              PSI_TMP(mm+1:mm+nuotot) = PSI(1:nuotot,io, ispin)
279              mm = mm + nuotot
280            enddo
281          enddo
282        endif
283#ifdef MPI
284        call MPI_Bcast( PSI_TMP, nuotot*nuo_LOC(Bnode)*spinor_dim,
285     &                  MPI_double_precision, BNode, MPI_Comm_World,
286     &                  MPIerror )
287#endif
288        do ispin = 1,h_spin_dim
289          do io = 1,nuo
290            call LocalToGlobalOrb(io,Node,Nodes,iio)
291            mm   = (ispin-1)*nuo_LOC(Bnode)*nuotot
292            do iie = 1,nuo_LOC(Bnode)
293              call  LocalToGlobalOrb(iie,BNode,Nodes,ie)
294              qe = qo(ie,ispin)
295              if (abs(qe).gt.occtol) then
296                paux => PSI_TMP(mm+1:mm+nuotot)
297                qei = qe*paux(iio)
298                eei = qei*eo(ie,ispin)
299                do j = 1,numd(io)
300                  ind = listdptr(io) + j
301                  jo  = listd(ind)
302                  jo = MODP(jo,nuotot) ! To allow auxiliary supercells
303                  Dnew(ind,ispin) = Dnew(ind,ispin) + qei*paux(jo)
304                  Enew(ind,ispin) = Enew(ind,ispin) + eei*paux(jo)
305                enddo
306              endif
307              mm = mm + nuotot
308            enddo
309          enddo
310        enddo
311      enddo
312
313      call timer( 'r-buildD', 2 )
314      call de_alloc( nuo_LOC, 'nuo_LOC', 'diagg' )
315      call de_alloc( PSI_TMP, 'PSI_TMP', 'diagg' )
316
317#ifdef DEBUG
318      call write_debug( '    POS diagg' )
319#endif
320      end
321