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