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 diag3g( nuo, no, maxnh, maxnd, maxo, 9 . numh, listhptr, listh, numd, listdptr, 10 . listd, H, S, getD, getPSI, qtot, temp, e1, e2, 11 . eo, qo, Dnew, Enew, ef, Entropy, 12 . Haux, Saux, psi, caux, 13 . nuotot, occtol, iscf, neigwanted) 14 15 use precision 16 use sys 17 use parallel, only : Node, Nodes, BlockSize 18 use parallelsubs, only : LocalToGlobalOrb,GlobalToLocalOrb 19 use writewave, only : writew 20 use m_fermid, only : fermid, stepf 21 use intrinsic_missing, only: MODP 22 use iso_c_binding, only : c_loc, c_f_pointer 23#ifdef MPI 24 use mpi_siesta 25#endif 26 27 implicit none 28C ********************************************************************* 29C Subroutine to calculate the eigenvalues and eigenvectors, density 30C and energy-density matrices, and occupation weights of each 31C eigenvector, for given Hamiltonian and Overlap matrices. 32C This version is for spin-orbit at gamma point. 33C Writen by J.Soler, May and August 1998. 34C Reduced memory requirements, Nick, Aug. 2017 35C **************************** INPUT ********************************** 36C integer nuo : Number of basis orbitals on local node 37C integer no : Number of basis orbitals 38C integer maxnh : Maximum number of orbitals interacting 39C integer maxnd : First dimension of Dnew and Enew 40C integer maxo : First dimension of eo and qo 41C integer numh(nuo) : Number of nonzero elements of each row 42C of hamiltonian matrix 43C integer listhptr(nuo) : Pointer to each row (-1) of the 44C hamiltonian matrix 45C integer listh(maxnh) : Nonzero hamiltonian-matrix element 46C column indexes for each matrix row 47C integer numd(nuo) : Number of nonzero elements of each row 48C of density matrix 49C integer listdptr(nuo) : Pointer to each row (-1) of the 50C density matrix 51C integer listd(maxnd) : Nonzero density-matrix element column 52C indexes for each matrix row 53C real*8 H(maxnh,4) : Hamiltonian in sparse form 54C real*8 S(maxnh) : Overlap in sparse form 55C logical getD : Find occupations and density matrices? 56C real*8 qtot : Number of electrons in unit cell 57C real*8 temp : Electronic temperature 58C real*8 e1, e2 : Energy range for density-matrix states 59C (to find local density of states) 60C Not used if e1 > e2 61C integer nuotot : total number of orbitals per unit cell 62C over all processors 63C real*8 occtol : Occupancy threshold for DM build 64C integer iscf : SCF cycle number 65C integer neigwanted : Number of eigenvalues wanted 66C *************************** OUTPUT ********************************** 67C real*8 eo(maxo*2) : Eigenvalues 68C real*8 qo(maxo*2) : Occupations of eigenstates 69C real*8 Dnew(maxnd,4) : Output Density Matrix 70C real*8 Enew(maxnd,4) : Output Energy-Density Matrix 71C real*8 ef : Fermi energy 72C real*8 Entropy : Electronic entropy 73C *************************** AUXILIARY ******************************* 74C complex*16 Haux(2,nuotot,2,nuo): Auxiliary space for the hamiltonian matrix 75C complex*16 Saux(2,nuotot,2,nuo): Auxiliary space for the overlap matrix 76C complex*16 psi(2,nuotot,2*nuo) : Auxiliary space for the eigenvectors 77C complex*16 caux(2,nuotot) : Extra auxiliary space 78C *************************** UNITS *********************************** 79C xij and kpoint must be in reciprocal coordinates of each other. 80C temp and H must be in the same energy units. 81C eo, Enew and ef returned in the units of H. 82C *************************** PARALLEL ******************************** 83C The auxiliary arrays are now no longer symmetry and so the order 84C of referencing has been changed in several places to reflect this. 85! ********************************************************************* 86! Haux(js,juo,is,iuo) = <js,juo|H|is,iuo> 87! Indices is and js are for spin components 88! Indices iuo and juo are for orbital components: 89! ********************************************************************* 90 integer maxo, maxnd, maxnh 91 integer no, nuo, nuotot, iscf, neigwanted 92 93 integer listh(maxnh), numh(nuo), listhptr(nuo) 94 integer listd(maxnd), numd(nuo), listdptr(nuo) 95 96 real(dp) Dnew(maxnd,8), e1, e2, ef, Enew(maxnd,4) 97 real(dp) Entropy, eo(maxo*2), H(maxnh,8), qo(maxo*2) 98 real(dp) qtot, S(maxnh), temp, occtol 99 100 complex(dp), dimension(2,nuotot,2*nuo), target :: psi 101 real(dp), pointer :: psi_real_1d(:) 102 complex(dp), dimension(2,nuotot,2,nuo) :: Haux, Saux 103 complex(dp), dimension(2,nuotot) :: caux 104 logical getD, getPSI 105 106! Internal variables ............................................. 107 108 integer BNode, BTest, ie, ierror, iie, iio 109 integer ind, io, j, jo, nd, iuo, juo 110 real(dp) ee, qe, t, k(3) 111 complex(dp) D11, D22, D12, D21 112#ifdef MPI 113 integer MPIerror 114#endif 115 external cdiag 116!*********************************************************************** 117! B E G I N 118!*********************************************************************** 119 120!*********************************************************************** 121! BUILD HAMILTONIAN 122!*********************************************************************** 123! The different subroutines build H_{j,i}^{js,is} = <js,j|H|is,i> 124! The spin notation is as follows: 125! 126! | H_{j,i}^{u,u} H_{j,i}^{u,d} | 127! H_(j,i} = | | 128! | H_{j,i}^{d,u} H_{j,i}^{d,d} | 129! 130! | H(ind,1) + i H(ind,5) H(ind,3) - i H(ind,4) | 131! = | | 132! | H(ind,7) + i H(ind,8) H(ind,2) + i H(ind,6) | 133! 134! 1. Hermiticity imposes H_{i,j}^{is,js}=H_{j,i}^{js,is}^* 135! 2. Since wave functions are real, if there are no single P or L operators: 136! (a) H_{i,j}^{is,js}=H_{j,i}^{is,js} 137! (b) These imply spin-box hermiticity: H_{i,j}^{is,js}=H_{i,j}^{js,is}^* 138! 139! (c) Hence 140! 141! | H(ind,1) H(ind,3) - i H(ind,4) | 142! H_{j,i} = H_{i,j} = | | 143! | H(ind,3) + i H(ind,4) H(ind,2) | 144! 145! 146! The spin-orbit interaction and the orbital part of the Zeeman interaction 147! break the "spin-box hermiticity" of H. Hence the full format of H is needed 148! in those cases. 149! 150! Spin-box symmetrization of D: 151! D_{i,j}(1,2) = 0.5 ( D_{i,j}(1,2) + D_{i,j}(2,1)^* ) 152! can not be enforced here. 153 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(1,jo,1,io) = Saux(1,jo,1,io) + dcmplx( S(ind), 0.0_dp) 162 Saux(2,jo,2,io) = Saux(2,jo,2,io) + dcmplx( S(ind), 0.0_dp) 163 Haux(1,jo,1,io) = Haux(1,jo,1,io) + dcmplx(H(ind,1), H(ind,5)) 164 Haux(2,jo,2,io) = Haux(2,jo,2,io) + dcmplx(H(ind,2), H(ind,6)) 165 Haux(1,jo,2,io) = Haux(1,jo,2,io) + dcmplx(H(ind,3),-H(ind,4)) 166 Haux(2,jo,1,io) = Haux(2,jo,1,io) + dcmplx(H(ind,7), H(ind,8)) 167 enddo 168 enddo 169 170! Solve the eigenvalue problem 171 call cdiag(Haux,Saux,2*nuotot,2*nuo,2*nuotot,eo,psi, 172 . 2*neigwanted,iscf,ierror,2*BlockSize) 173 174! 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 178! Repeat diagonalisation with increased memory to handle clustering 179 do io = 1,nuo 180 Saux(:,:,:,io) = 0._dp 181 Haux(:,:,:,io) = 0._dp 182 do j = 1,numh(io) 183 ind = listhptr(io) + j 184 jo = listh(ind) 185 jo = MODP(jo,nuotot) ! To allow auxiliary supercells 186 Saux(1,jo,1,io) = Saux(1,jo,1,io) + dcmplx( S(ind), 0.0_dp) 187 Saux(2,jo,2,io) = Saux(2,jo,2,io) + dcmplx( S(ind), 0.0_dp) 188 Haux(1,jo,1,io) = Haux(1,jo,1,io) + dcmplx(H(ind,1), H(ind,5)) 189 Haux(2,jo,2,io) = Haux(2,jo,2,io) + dcmplx(H(ind,2), H(ind,6)) 190 Haux(1,jo,2,io) = Haux(1,jo,2,io) + dcmplx(H(ind,3),-H(ind,4)) 191 Haux(2,jo,1,io) = Haux(2,jo,1,io) + dcmplx(H(ind,7), H(ind,8)) 192 enddo 193 enddo 194 195 call cdiag(Haux,Saux,2*nuotot,2*nuo,2*nuotot,eo,psi, 196 . 2*neigwanted,iscf,ierror,2*BlockSize) 197 endif 198 199 if (getPSI) then 200 k(1:3) = 0.0d0 201 202 ! We do not have info about occupation 203 ! In the current implementation, neigneeded 204 ! is ALWAYS nuotot. 205 ! This is because writewave expects the full spectrum 206 call c_f_pointer( c_loc(psi), psi_real_1d, [size(psi)*2] ) 207 call writew(nuotot,nuo,1,k,1, 208 & eo,psi_real_1d,gamma=.true., 209 $ non_coll=.true.,blocksize=2*BlockSize) 210 endif 211 212! Check if we are done 213 if (.not.getD) return 214 215! Find new Fermi energy and occupation weights 216 217 call fermid(2, 1, 1, (/1._dp/), 2*maxo, 2*neigwanted, 218 . eo, temp, qtot, qo, ef, Entropy) 219 220! write(6,'(a,f12.6)') ' Ef = ', ef 221 222! Find weights for local density of states ............................ 223 if (e1 .lt. e2) then 224 t = max( temp, 1.d-6 ) 225 do io = 1, nuotot*2 226 qo(io) = ( stepf((eo(io)-e2)/t) - stepf((eo(io)-e1)/t)) 227 enddo 228 endif 229 230!*********************************************************************** 231! BUILD NEW DENSITY MATRIX 232!*********************************************************************** 233! 234! | ------- 1,1 ------- ------- 1,2 --------- | 235! | c_{i,up} c_{j,up}^* c_{i,up} c_{j,down)^* | 236! D_{i,j} = | | 237! | ------- 2,1 ------- ------- 2,2 --------- | 238! | c_{i,down} c_{j,up}^* c_{i,dn} c_{j,dn)^* | 239! 240! = | D_{i,j}(1) D_{i,j}(3)-i D_{i,j}(4) | 241! | D_{i,j}(7)+i D_{i,j}(8) D_{i,j}(2) | 242 243! The Energy is computed as E = Tr [ D_{i,j} H_{j,i} ] 244! 245! The Density Matrix is not "spin box hermitian" even if H does not contain P or 246! L operators. 247! 248! If H contains P or L operators, then spin-box "symetrization" of D leads to erroneous 249! results for E 250 251 Dnew(:,:) = 0.0_dp 252 Enew(:,:) = 0.0_dp 253 254 BNode = 0 255 iie = 0 256 257 do ie = 1,2*nuotot 258 qe = qo(ie) 259 if (Node.eq.BNode) then 260 iie = iie + 1 261 endif 262 263 caux(:,:) = dcmplx(0.0_dp,0.0_dp) 264 if (abs(qe).gt.occtol) then 265 if (Node.eq.BNode) then 266 do j = 1,nuotot 267 caux(1,j) = psi(1,j,iie) ! c_{i,up} 268 caux(2,j) = psi(2,j,iie) ! c_{i,dn} 269 enddo 270 endif 271#ifdef MPI 272 call MPI_Bcast(caux(1,1),2*nuotot,MPI_double_complex, 273 . BNode,MPI_Comm_World,MPIerror) 274#endif 275 ee = qo(ie) * eo(ie) 276 do io = 1,nuo 277 call LocalToGlobalOrb(io,Node,Nodes,iio) 278 do j = 1,numd(io) 279 ind = listdptr(io) + j 280 jo = listd(ind) 281 jo = MODP(jo,nuotot) ! To allow auxiliary supercells 282!------- 1,1 ----------------------------------------------------------- 283 D11 = caux(1,iio) * dconjg(caux(1,jo)) 284!------- 2,2 ----------------------------------------------------------- 285 D22 = caux(2,iio) * dconjg(caux(2,jo)) 286!------- 1,2 ----------------------------------------------------------- 287 D12 = caux(1,iio) * dconjg(caux(2,jo)) 288!------- 2,1 ----------------------------------------------------------- 289 D21 = caux(2,iio) * dconjg(caux(1,jo)) 290 291 Dnew(ind,1) = Dnew(ind,1) + dreal(D11) * qe 292 Dnew(ind,2) = Dnew(ind,2) + dreal(D22) * qe 293 Dnew(ind,3) = Dnew(ind,3) + dreal(D12) * qe 294 Dnew(ind,4) = Dnew(ind,4) - dimag(D12) * qe 295 Dnew(ind,5) = Dnew(ind,5) + dimag(D11) * qe 296 Dnew(ind,6) = Dnew(ind,6) + dimag(D22) * qe 297 Dnew(ind,7) = Dnew(ind,7) + dreal(D21) * qe 298 Dnew(ind,8) = Dnew(ind,8) + dimag(D21) * qe 299 300 Enew(ind,1) = Enew(ind,1) + dreal(D11) * ee 301 Enew(ind,2) = Enew(ind,2) + dreal(D22) * ee 302 Enew(ind,3) = Enew(ind,3) + dreal(D12) * ee 303 Enew(ind,4) = Enew(ind,4) - dimag(D12) * ee 304 305 enddo 306 enddo 307 endif 308 BTest = ie/(2*BlockSize) 309 if (BTest*2*BlockSize.eq.ie) then 310 BNode = BNode + 1 311 if (BNode .gt. Nodes-1) BNode = 0 312 endif 313 enddo 314 315!*********************************************************************** 316 return 317 end 318!********************************************************************** 319 320