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! 8subroutine diag2kspiral( nuo, no, maxnh, maxnd, maxo, & 9 numh, listhptr, listh, numd, listdptr, & 10 listd, H, S, getD, qtot, temp, e1, e2, & 11 xij, indxuo, nk, kpoint, wk, & 12 eo, qo, Dnew, Enew, ef, Entropy, qw, & 13 Haux, Saux, psi, Dk, Ek, aux, nuotot, & 14 occtol, iscf, neigwanted) 15! ********************************************************************* 16! Calculates the eigenvalues and eigenvectors, density 17! and energy-density matrices, and occupation weights of each 18! eigenvector, for given Hamiltonian and Overlap matrices. 19! This version is for non-colinear spin with k-sampling and spiral 20! arrangement of spins. 21! Written by V. M. Garcia-Suarez. June 2002 22! Modified to reduce eigenvector computation by J.D. Gale Nov. 2004 23! Lots of bugs fixed and changed to complex by Nick Papior, June 2020 24! **************************** INPUT ********************************** 25! integer nuo : Number of basis orbitals in unit cell 26! integer no : Number of basis orbitals in supercell 27! integer maxnh : Maximum number of orbitals interacting 28! integer maxnd : First dimension of listd / DM 29! integer maxo : First dimension of eo and qo 30! integer numh(nuo) : Number of nonzero elements of each row 31! of hamiltonian matrix 32! integer listhptr(nuo) : Pointer to each row (-1) of the 33! hamiltonian matrix 34! integer listh(maxnh) : Nonzero hamiltonian-matrix element 35! column indexes for each matrix row 36! integer numd(nuo) : Number of nonzero elements of each row 37! of density matrix 38! integer listdptr(nuo) : Pointer to each row (-1) of the 39! density matrix 40! integer listd(maxnd) : Nonzero density-matrix element column 41! indexes for each matrix row 42! real*8 H(maxnh,4) : Hamiltonian in sparse form 43! real*8 S(maxnh) : Overlap in sparse form 44! logical getD : Find occupations and density matrices? 45! real*8 qtot : Number of electrons in unit cell 46! real*8 temp : Electronic temperature 47! real*8 e1, e2 : Energy range for density-matrix states 48! (to find local density of states) 49! Not used if e1 > e2 50! real*8 xij(3,maxnh) : Vectors between orbital centers (sparse) 51! (not used if only gamma point) 52! integer indxuo(no) : Index of equivalent orbital in unit cell 53! Unit cell orbitals must be the first in 54! orbital lists, i.e. indxuo.le.nuo, with 55! nuo the number of orbitals in unit cell 56! real*8 qw(3) : Wave vector for spiral configuration 57! integer nk : Number of k points 58! real*8 kpoint(3,nk) : k point vectors 59! real*8 wk(nk) : k point weights (must sum one) 60! integer nuotot : total number of orbitals per unit cell 61! over all processors 62! real*8 occtol : Occupancy threshold for DM build 63! integer iscf : SCF cycle number 64! integer neigwanted : Number of eigenvalues wanted 65! *************************** OUTPUT ********************************** 66! real*8 eo(maxo*2,nk) : Eigenvalues 67! real*8 qo(maxo*2,nk) : Occupations of eigenstates 68! real*8 Dnew(maxnd,4) : Output Density Matrix 69! real*8 Enew(maxnd,4) : Output Energy-Density Matrix 70! real*8 ef : Fermi energy 71! real*8 Entropy : Electronic entropy 72! *************************** AUXILIARY ******************************* 73! complex*16 Haux(2,nuotot,2,nuo) : Aux. space for the hamiltonian matrix 74! complex*16 Saux(2,nuotot,2,nuo) : Aux. space for the overlap matrix 75! complex*16 psi(2,nuotot,2*nuo) : Aux. space for the eigenvectors 76! complex*16 aux(2,nuotot) : Extra auxiliary space 77! complex*16 Dk(2,nuotot,2,nuo) : Aux. space that may be the same as Haux 78! complex*16 Ek(2,nuotot,2,nuo) : Aux. space that may be the same as Saux 79! *************************** UNITS *********************************** 80! xij and kpoint must be in reciprocal coordinates of each other. 81! temp and H must be in the same energy units. 82! eo, Enew and ef returned in the units of H. 83! *************************** PARALLEL ******************************** 84! The auxiliary arrays are now no longer symmetry and so the order 85! of referencing has been changed in several places to reflect this. 86! ********************************************************************* 87! 88! Modules 89! 90 use precision 91 use sys 92 use parallel, only : Node, Nodes, BlockSize 93 use parallelsubs, only : LocalToGlobalOrb 94 use m_fermid, only : fermid, stepf 95#ifdef MPI 96 use mpi_siesta 97#endif 98 99 implicit none 100 101#ifdef MPI 102 integer :: MPIerror 103#endif 104 105 integer :: maxnd, maxnh, maxo, nk, no, nuo, nuotot, iscf, neigwanted 106 107 integer :: indxuo(no), listh(maxnh), numh(nuo), listd(maxnd), numd(nuo), & 108 listhptr(nuo), listdptr(nuo) 109 110 real(dp) :: Dnew(maxnd,4), & 111 e1, e2, ef, Enew(maxnd,4), Entropy, eo(maxo*2,nk), & 112 H(maxnh,4), kpoint(3,nk), qo(maxo*2,nk), qtot, & 113 S(maxnh), temp, wk(nk), xij(3,maxnh), qw(3), occtol 114 115 complex(dp), dimension(2,nuotot,2,nuo) :: Haux, Saux, Dk, Ek 116 complex(dp), dimension(2,nuotot,2*nuo) :: psi 117 complex(dp), dimension(2,nuotot) :: aux 118 119 logical :: getD 120 121 external :: cdiag 122 123! Internal variables ............................................. 124 integer :: BNode, BTest, ie, ierror, iie, ik, ind, io, iio, & 125 iuo, j, jo, juo, neigneeded 126 integer :: nuotot2, nuo2, BlockSize2, maxo2, neigwanted2 127 complex(dp) :: kphs, qphs 128 complex(dp) :: cicj, D11, D22, D12, D21 129 real(dp) :: kxij, qxij, ee, qe, t, qwh(3) 130 131#ifdef DEBUG 132 call write_debug( ' PRE diag2kspiral' ) 133#endif 134 135 ! Easier 136 nuotot2 = nuotot * 2 137 neigwanted2 = neigwanted * 2 138 nuo2 = nuo * 2 139 maxo2 = maxo * 2 140 BlockSize2 = BlockSize * 2 141 qwh(:) = qw(:) * 0.5_dp 142 143! Find eigenvalues at every k point ............................... 144 do ik = 1,nk 145 146 call setup_k() 147 148 ! Find eigenvalues 149 ! Possible memory optimization: equivalence Haux and psi 150 call cdiag(Haux,Saux,nuotot2,nuotot2,nuo2,eo(1,ik),psi, & 151 -neigwanted2,iscf,ierror,BlockSize2) 152 if ( ierror /= 0 ) then 153 call die('Terminating due to failed diagonalisation') 154 end if 155 end do 156 157 ! Check if we are done ................................................ 158 if (.not.getD) then 159#ifdef DEBUG 160 call write_debug( ' POS diag2kspiral' ) 161#endif 162 return 163 end if 164 165 ! Find new Fermi energy and occupation weights ........................ 166 call fermid(2, 1, nk, wk, maxo2, neigwanted2, eo, & 167 temp, qtot, qo, ef, Entropy ) 168 169 ! Find weights for local density of states ............................ 170 if ( e1 < e2 ) then 171 t = max( temp, 1.d-6 ) 172 do ik = 1,nk 173 do io = 1, neigwanted2 174 qo(io,ik) = wk(ik) * & 175 ( stepf( (eo(io,ik)-e2)/t ) - stepf( (eo(io,ik)-e1)/t ) ) 176 end do 177 end do 178 end if 179 180! New density and energy-density matrices of unit-cell orbitals ....... 181 Dnew(:,:) = 0._dp 182 Enew(:,:) = 0._dp 183 184 do ik = 1, nk 185 186 ! Find maximum eigenvector that is required for this k point 187 neigneeded = 0 188 ie = neigwanted2 189 do while ( ie > 0 .and. neigneeded == 0 ) 190 qe = qo(ie,ik) 191 if ( abs(qe) > occtol ) neigneeded = ie 192 ie = ie - 1 193 end do 194 195 call setup_k() 196 197 call cdiag(Haux,Saux,nuotot2,nuotot2,nuo2,eo(1,ik),psi, & 198 neigneeded,iscf,ierror,BlockSize2) 199 200 ! Check error flag and take appropriate action 201 if ( ierror > 0 ) then 202 call die('Terminating due to failed diagonalisation') 203 else if ( ierror < 0 ) then 204 205 call setup_k() 206 call cdiag(Haux,Saux,nuotot2,nuotot2,nuo2,eo(1,ik),psi, & 207 neigneeded,iscf,ierror,BlockSize2) 208 209 end if 210 211! Store the products of eigenvectors in matrices Dk and Ek 212! WARNING: Dk and Ek may be EQUIVALENCE'd to Haux and Saux 213 Dk = cmplx(0.0_dp, 0.0_dp, dp) 214 Ek = cmplx(0.0_dp, 0.0_dp, dp) 215 216 BNode = 0 217 iie = 0 218 do ie = 1, neigneeded 219 if ( Node == BNode ) then 220 iie = iie + 1 221 do j = 1, nuotot 222 aux(1,j) = psi(1,j,iie) ! c_i,up 223 aux(2,j) = psi(2,j,iie) ! c_i,dn 224 end do 225 end if 226#ifdef MPI 227 call MPI_Bcast(aux(1,1),nuotot2,MPI_double_complex,BNode, & 228 MPI_Comm_World,MPIerror) 229#endif 230 qe = qo(ie,ik) 231 ee = qe * eo(ie,ik) 232 do iuo = 1,nuo 233 call LocalToGlobalOrb(iuo,Node,Nodes,iio) 234 do juo = 1, nuotot 235!------- 1,1 ----------------------------------------------------------- 236 cicj = conjg(aux(1,juo)) * aux(1,iio) 237 Dk(1,juo,1,iuo) = Dk(1,juo,1,iuo) + qe * cicj 238 Ek(1,juo,1,iuo) = Ek(1,juo,1,iuo) + ee * cicj 239!------- 2,2 ----------------------------------------------------------- 240 cicj = conjg(aux(2,juo)) * aux(2,iio) 241 Dk(2,juo,2,iuo) = Dk(2,juo,2,iuo) + qe * cicj 242 Ek(2,juo,2,iuo) = Ek(2,juo,2,iuo) + ee * cicj 243!------- 1,2 ----------------------------------------------------------- 244 cicj = conjg(aux(2,juo)) * aux(1,iio) 245 Dk(1,juo,2,iuo) = Dk(1,juo,2,iuo) + qe * cicj 246 Ek(1,juo,2,iuo) = Ek(1,juo,2,iuo) + ee * cicj 247!------- 2,1 ----------------------------------------------------------- 248 cicj = conjg(aux(1,juo)) * aux(2,iio) 249 Dk(2,juo,1,iuo) = Dk(2,juo,1,iuo) + qe * cicj 250 Ek(2,juo,1,iuo) = Ek(2,juo,1,iuo) + ee * cicj 251 end do 252 end do 253 BTest = ie/BlockSize2 254 if ( BTest*BlockSize2 == ie ) then 255 BNode = BNode + 1 256 if ( BNode >= Nodes ) BNode = 0 257 end if 258 end do 259 260! Add contribution to density matrices of unit-cell orbitals 261 do iuo = 1,nuo 262 do j = 1,numd(iuo) 263 ind = listdptr(iuo) + j 264 jo = listd(ind) 265 juo = indxuo(jo) 266 kxij = kpoint(1,ik) * xij(1,ind) + & 267 kpoint(2,ik) * xij(2,ind) + & 268 kpoint(3,ik) * xij(3,ind) 269 kphs = exp(cmplx(0.0_dp, - kxij, dp)) 270 qxij = qwh(1) * xij(1,ind) + & 271 qwh(2) * xij(2,ind) + & 272 qwh(3) * xij(3,ind) 273 qphs = exp(cmplx(0.0_dp, - qxij, dp)) 274 275 D11 = Dk(1,juo,1,iuo) * kphs * qphs 276 D21 = Dk(2,juo,1,iuo) * kphs * conjg(qphs) 277 D12 = Dk(1,juo,2,iuo) * kphs * qphs 278 D22 = Dk(2,juo,2,iuo) * kphs * conjg(qphs) 279 280 ! Make DM Spin-box hermitian 281 D12 = 0.5_dp * (D12 + conjg(D21)) 282 283 Dnew(ind,1) = Dnew(ind,1) + real(D11, dp) 284 Dnew(ind,2) = Dnew(ind,2) + real(D22, dp) 285 Dnew(ind,3) = Dnew(ind,3) + real(D12, dp) 286 Dnew(ind,4) = Dnew(ind,4) - aimag(D12) 287 288 D11 = Ek(1,juo,1,iuo) * kphs * qphs 289 D21 = Ek(2,juo,1,iuo) * kphs * conjg(qphs) 290 D12 = Ek(1,juo,2,iuo) * kphs * qphs 291 D22 = Ek(2,juo,2,iuo) * kphs * conjg(qphs) 292 293 ! Make EDM Spin-box hermitian 294 D12 = 0.5_dp * (D12 + conjg(D21)) 295 296 Enew(ind,1) = Enew(ind,1) + real(D11, dp) 297 Enew(ind,2) = Enew(ind,2) + real(D22, dp) 298 Enew(ind,3) = Enew(ind,3) + real(D12, dp) 299 Enew(ind,4) = Enew(ind,4) - aimag(D12) 300 301 end do 302 end do 303 304 end do 305 306#ifdef DEBUG 307 call write_debug( ' POS diag2ksipral' ) 308#endif 309 310contains 311 312 subroutine setup_k() 313 314 ! Initialize Hamiltonian and overlap matrices in full format 315 ! Index i is for real/imag parts 316 ! Indices is and js are for spin components 317 ! Indices iuo and juo are for orbital components: 318 ! Haux(i,js,juo,is,iuo) = <js,juo|H|is,iuo> 319 Saux = cmplx(0.0_dp, 0.0_dp, dp) 320 Haux = cmplx(0.0_dp, 0.0_dp, dp) 321 322 ! Transfer S,H matrices from sparse format in supercell to 323 ! full format in unit cell 324 ! Convention: ispin=1 => H11, ispin=2 => H22, 325 ! ispin=3 => Real(H12), ispin=4 => Imag(H12) 326 do iuo = 1,nuo 327 do j = 1,numh(iuo) 328 ind = listhptr(iuo) + j 329 jo = listh(ind) 330 juo = indxuo(jo) 331 kxij = kpoint(1,ik) * xij(1,ind) + & 332 kpoint(2,ik) * xij(2,ind) + & 333 kpoint(3,ik) * xij(3,ind) 334 kphs = exp(cmplx(0.0_dp, - kxij, dp)) 335 qxij = qwh(1) * xij(1,ind) + & 336 qwh(2) * xij(2,ind) + & 337 qwh(3) * xij(3,ind) 338 qphs = exp(cmplx(0.0_dp, - qxij, dp)) 339 D11 = kphs * qphs 340 D22 = kphs * conjg(qphs) 341 342 Saux(1,juo,1,iuo) = Saux(1,juo,1,iuo) + S(ind) * D11 343 Saux(2,juo,2,iuo) = Saux(2,juo,2,iuo) + S(ind) * D22 344 345 Haux(1,juo,1,iuo) = Haux(1,juo,1,iuo) + H(ind,1) * D11 346 Haux(2,juo,1,iuo) = Haux(2,juo,1,iuo) + cmplx(H(ind,3), H(ind,4),dp) * D22 347 Haux(1,juo,2,iuo) = Haux(1,juo,2,iuo) + cmplx(H(ind,3), -H(ind,4),dp) * D11 348 Haux(2,juo,2,iuo) = Haux(2,juo,2,iuo) + H(ind,2) * D22 349 350 end do 351 end do 352 353 end subroutine setup_k 354 355end subroutine diag2kspiral 356