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 pdos3g( nuo, no, maxuo, maxnh,
9     .                  maxo, numh, listhptr, listh, H, S,
10     .                  E1, E2, nhist, sigma, indxuo, eo,
11     .                  haux, saux, psi, dtot, dpr, nuotot )
12
13C **********************************************************************
14C Find the density of states projected onto the atomic orbitals
15C     D_mu(E) = Sum(n,k,nu) C(mu,n,k) C(nu,n,k) S(mu,nu,k) Delta(E-E(n,k))
16C where n run over all the bands between two given energies
17C Written by J. Junquera and E. Artacho. Nov' 99
18C Gamma point version adapted from PDOSK by Julian Gale. Feb' 03
19C Spin-orbit coupling version by J. Ferrer, October 2007
20C ****  INPUT  *********************************************************
21C integer nuo               : Number of atomic orbitals in the unit cell
22C integer no                : Number of atomic orbitals in the supercell
23C                             (maximum number of differents spin polarizations)
24C integer maxuo             : Maximum number of atomic orbitals in the unit cell
25C integer maxnh             : Maximum number of orbitals interacting
26C                             with any orbital
27C integer maxo              : First dimension of eo
28C integer numh(nuo)         : Number of nonzero elements of each row
29C                             of hamiltonian matrix
30C integer listhptr(nuo)     : Pointer to each row (-1) of the
31C                             hamiltonian matrix
32C integer listh(maxnh)      : Nonzero hamiltonian-matrix element
33C                             column indexes for each matrix row
34C real*8  H(maxnh,4)    : Hamiltonian in sparse format
35C real*8  S(maxnh)          : Overlap in sparse format
36C real*8  E1, E2            : Energy range for density-matrix states
37C                             (to find local density of states)
38C                             Not used if e1 > e2
39C integer nhist             : Number of the subdivisions of the histogram
40C real*8  sigma             : Width of the gaussian to expand the eigenvectors
41C integer indxuo(no)        : Index of equivalent orbital in unit cell
42C real*8  eo(maxo,2)   : Eigenvalues
43C integer nuotot            : Total number of orbitals per unit cell
44C ****  AUXILIARY  *****************************************************
45C real*8  haux(nuo,nuo)     : Auxiliary space for the hamiltonian matrix
46C real*8  saux(nuo,nuo)     : Auxiliary space for the overlap matrix
47C real*8  psi(nuo,nuo)      : Auxiliary space for the eigenvectors
48C ****  OUTPUT  ********************************************************
49C real*8  dtot(nhist,4)   : Total density of states
50C real*8  dpr(nhist,nuo,4): Proyected density of states
51C **********************************************************************
52
53      use precision
54      use parallel,     only : Node, Nodes, BlockSize
55      use parallelsubs, only : GetNodeOrbs, LocalToGlobalOrb
56      use units,        only : pi
57      use alloc,        only : re_alloc, de_alloc
58#ifdef MPI
59      use mpi_siesta
60#endif
61      use sys,          only : die
62
63      implicit none
64
65      integer
66     .  nuo, no, maxuo, maxnh,
67     .  maxo, nhist, nuotot
68
69      integer
70     .  numh(nuo), listhptr(nuo), listh(maxnh), indxuo(no)
71
72      real(dp)
73     .  H(maxnh,8), S(maxnh), E1, E2, sigma, eo(maxo*2),
74     .  dtot(nhist,4), dpr(nhist,nuotot,4)
75      complex(dp), target :: psi(2,nuotot,2*nuo)
76      complex(dp) Haux(2,nuotot,2,nuo), Saux(2,nuotot,2,nuo)
77      complex(dp), pointer :: caux(:,:)
78      external  cdiag
79
80C Internal variables ---------------------------------------------------
81      integer
82     .  ispin, iuo, juo, io, j, jo, ihist, iband, ind, ierror
83
84      real(dp)
85     .  delta, ener, diff, rpipj, ipipj, gauss, norm
86      real(dp), pointer :: Spr(:,:) => null()
87
88#ifdef MPI
89      integer :: BNode, Bnuo, ibandg, maxnuo, MPIerror
90      real(dp), pointer :: Sloc(:,:)
91      real(dp)          :: tmp(nhist,nuotot,4)
92#endif
93
94C Initialize some variables
95      delta = (E2 - E1)/nhist
96
97C Initialize auxiliary variables
98
99      call re_alloc(Spr, 1, nuotot, 1, nuo, name='Spr',
100     &     routine='pdos3g')
101
102      do io = 1,nuo
103        Saux(:,:,:,io) = 0._dp
104        Haux(:,:,:,io) = 0._dp
105        do j = 1,numh(io)
106          ind = listhptr(io) + j
107          jo = listh(ind)
108          jo = indxuo(jo)
109          Saux(1,jo,1,io) = Saux(1,jo,1,io) + dcmplx( S(ind), 0.0_dp)
110          Saux(2,jo,2,io) = Saux(2,jo,2,io) + dcmplx( S(ind), 0.0_dp)
111          Haux(1,jo,1,io) = Haux(1,jo,1,io) + dcmplx(H(ind,1), H(ind,5))
112          Haux(2,jo,2,io) = Haux(2,jo,2,io) + dcmplx(H(ind,2), H(ind,6))
113          Haux(1,jo,2,io) = Haux(1,jo,2,io) + dcmplx(H(ind,3),-H(ind,4))
114          Haux(2,jo,1,io) = Haux(2,jo,1,io) + dcmplx(H(ind,7), H(ind,8))
115        enddo
116      enddo
117C Diagonalize at the Gamma point
118      call cdiag( Haux, Saux, 2*nuotot, 2*nuo, 2*nuotot,
119     .            eo, psi, 2*nuotot, 1, ierror, 2*BlockSize )
120      if (ierror.gt.0) then
121        call die('Terminating due to failed diagonalisation')
122      elseif (ierror .lt. 0) then
123C Repeat diagonalisation with increased memory to handle clustering
124      do io = 1,nuo
125        Saux(:,:,:,io) = 0._dp
126        Haux(:,:,:,io) = 0._dp
127        do j = 1,numh(io)
128          ind = listhptr(io) + j
129          jo = listh(ind)
130          jo = indxuo(jo)
131          Saux(1,jo,1,io) = Saux(1,jo,1,io) + dcmplx( S(ind), 0.0_dp)
132          Saux(2,jo,2,io) = Saux(2,jo,2,io) + dcmplx( S(ind), 0.0_dp)
133          Haux(1,jo,1,io) = Haux(1,jo,1,io) + dcmplx(H(ind,1), H(ind,5))
134          Haux(2,jo,2,io) = Haux(2,jo,2,io) + dcmplx(H(ind,2), H(ind,6))
135          Haux(1,jo,2,io) = Haux(1,jo,2,io) + dcmplx(H(ind,3),-H(ind,4))
136          Haux(2,jo,1,io) = Haux(2,jo,1,io) + dcmplx(H(ind,7), H(ind,8))
137        enddo
138      enddo
139        call cdiag(Haux,Saux,2*nuotot,2*nuo,2*nuotot,eo,psi,
140     .                         2*nuotot,1,ierror,2*BlockSize)
141      endif
142
143C     Rebuild the full overlap matrix
144      Spr = 0._dp
145      do io = 1, nuo
146        do j = 1, numh(io)
147          ind = listhptr(io) + j
148          jo = listh(ind)
149          jo = indxuo(jo)
150          Spr(jo,io) = Spr(jo,io) + S(ind)
151        enddo
152      enddo
153
154#ifdef MPI
155C Find maximum number of orbitals per node
156        call MPI_AllReduce(nuo,maxnuo,1,MPI_integer,MPI_max,
157     .                   MPI_Comm_World,MPIerror)
158
159C Allocate workspace array for broadcast overlap matrix
160        nullify( Sloc )
161        call re_alloc( Sloc, 1, nuotot, 1, maxnuo,
162     .                 name='Sloc', routine='pdos3g' )
163
164C Loop over nodes broadcasting overlap matrix
165        do BNode = 0,Nodes-1
166
167C Find out how many orbitals there are on the broadcast node
168          call GetNodeOrbs(nuotot,BNode,Nodes,Bnuo)
169
170C Transfer data
171          if (Node.eq.BNode) then
172            Sloc(1:nuotot,1:Bnuo) = Spr(1:nuotot,1:Bnuo)
173          endif
174          call MPI_Bcast(Sloc(1,1),nuotot*Bnuo,
175     .      MPI_double_precision,BNode,MPI_Comm_World,MPIerror)
176
177          do ihist = 1, nhist
178            ener = E1 + (ihist - 1) * delta
179            do iband = 1, nuo*2
180              call LocalToGlobalOrb((iband+1)/2,Node,Nodes,ibandg)
181              ibandg = ibandg * 2 - mod(iband, 2)
182              diff = (ener - eo(ibandg))**2 / (sigma ** 2)
183              if (diff .gt. 15.0d0) cycle
184              gauss = exp(-diff)
185              caux => psi(:,:,iband) ! c_{up,j}, c_{down,j}
186              do jo = 1, Bnuo
187                call LocalToGlobalOrb(jo,BNode,Nodes,juo)
188                do io = 1, nuotot
189                  rpipj = dreal(caux(1,io)*dconjg(caux(1,juo)))*gauss
190                  dpr(ihist,juo,1)=dpr(ihist,juo,1)+rpipj*Sloc(io,jo)
191
192                  rpipj = dreal(caux(2,io)*dconjg(caux(2,juo)))*gauss
193                  dpr(ihist,juo,2)=dpr(ihist,juo,2)+rpipj*Sloc(io,jo)
194
195                  rpipj = dreal(caux(1,io)*dconjg(caux(2,juo)))*gauss
196                  dpr(ihist,juo,3)=dpr(ihist,juo,3)+rpipj*Sloc(io,jo)
197
198                  ipipj = dimag(caux(1,io)*dconjg(caux(2,juo)))*gauss
199                  dpr(ihist,juo,4)=dpr(ihist,juo,4)-ipipj*Sloc(io,jo)
200                enddo
201              enddo
202            enddo
203          enddo
204
205        enddo  !BNode
206
207C Free workspace array for overlap
208        call de_alloc(Sloc, 'Sloc', 'pdos3g')
209
210#else
211C Loop over all the energy range
212
213        do ihist = 1, nhist
214          ener = E1 + (ihist - 1) * delta
215          do iband = 1, nuo*2
216            diff = (ener - eo(iband))**2 / (sigma ** 2)
217            if (diff .gt. 15.0d0) cycle
218            gauss = exp(-diff)
219            caux => psi(:,:,iband) ! c_{up,j}, c_{down,j}
220            do jo = 1, nuotot
221              do io = 1, nuotot
222                rpipj = dreal(caux(1,io)*dconjg(caux(1,jo)))*gauss
223                dpr(ihist,jo,1)=dpr(ihist,jo,1)+rpipj*Spr(io,jo)
224
225                rpipj = dreal(caux(2,io)*dconjg(caux(2,jo)))*gauss
226                dpr(ihist,jo,2)=dpr(ihist,jo,2)+rpipj*Spr(io,jo)
227
228                rpipj = dreal(caux(1,io)*dconjg(caux(2,jo)))*gauss
229                dpr(ihist,jo,3)=dpr(ihist,jo,3)+rpipj*Spr(io,jo)
230
231                ipipj = dimag(caux(1,io)*dconjg(caux(2,jo)))*gauss
232                dpr(ihist,jo,4)=dpr(ihist,jo,4)-ipipj*Spr(io,jo)
233              enddo
234            enddo
235          enddo
236        enddo
237#endif
238
239      call de_alloc(Spr, 'Spr', 'pdos3g')
240
241#ifdef MPI
242
243C Global reduction of dpr matrix
244      tmp = 0.0d0
245      call MPI_AllReduce(dpr(1,1,1),tmp(1,1,1),nhist*nuotot*4,
246     .  MPI_double_precision,MPI_sum,MPI_Comm_World,MPIerror)
247      dpr = tmp
248
249#endif
250
251      norm = sigma * sqrt(pi)
252      dpr = dpr / norm
253
254      do ihist = 1, nhist
255        dtot(ihist,1) = sum(dpr(ihist,:,1))
256        dtot(ihist,2) = sum(dpr(ihist,:,2))
257        dtot(ihist,3) = sum(dpr(ihist,:,3))
258        dtot(ihist,4) = sum(dpr(ihist,:,4))
259      enddo
260
261      return
262      end
263