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