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#ifndef CDF 9 subroutine dummy_diagk_file() 10 end subroutine dummy_diagk_file 11#else 12 subroutine diagk_file( nspin, nuo, no, maxspn, 13 $ maxnh, maxnd, 14 . maxo, numh, listhptr, listh, numd, listdptr, 15 . listd, H, S, getD, getPSI, fixspin, qtot, qs, 16 . temp, e1, e2, xij, indxuo, nk, kpoint, wk, 17 . eo, qo, Dnew, Enew, ef, efs, Entropy, 18 . Haux, Saux, psi, Dk, Ek, nuotot, occtol, 19 . iscf, n_eigenvectors ) 20C ********************************************************************* 21C Subroutine to calculate the eigenvalues and eigenvectors, density 22C and energy-density matrices, and occupation weights of each 23C eigenvector, for given Hamiltonian and Overlap matrices (including 24C spin polarization). K-sampling version. 25C Writen by J.Soler, August 1998. 26C Modified by A. Garcia, June 2008, to trade file space for speed (!) 27C **************************** INPUT ********************************** 28C integer nspin : Number of spin components (1 or 2) 29C integer nuo : Number of basis orbitals in unit cell 30C local to this processor 31C integer no : Number of basis orbitals in supercell 32C integer maxspn : Second dimension of eo and qo 33C integer maxo : First dimension of eo and qo 34C integer maxnh : Maximum number of orbitals interacting 35C integer maxnd : First dimension of listd and DM 36C integer numh(nuo) : Number of nonzero elements of each row 37C of hamiltonian matrix 38C integer listhptr(nuo) : Pointer to each row (-1) of the 39C hamiltonian matrix 40C integer listh(maxnh) : Nonzero hamiltonian-matrix element 41C column indexes for each matrix row 42C integer numd(nuo) : Number of nonzero elements of each row 43C of density matrix 44C integer listdptr(nuo) : Pointer to each row (-1) of the 45C density matrix 46C integer listd(maxnd) : Nonzero density-matrix element column 47C indexes for each matrix row 48C real*8 H(maxnh,nspin) : Hamiltonian in sparse form 49C real*8 S(maxnh) : Overlap in sparse form 50C logical getD : Find occupations and density matrices? 51C logical getPSI : Find and print wave functions? 52C real*8 qtot : Number of electrons in unit cell 53C real*8 temp : Electronic temperature 54C real*8 e1, e2 : Energy range for density-matrix states 55C (to find local density of states) 56C Not used if e1 > e2 57C real*8 xij(3,maxnh) : Vectors between orbital centers (sparse) 58C (not used if only gamma point) 59C integer indxuo(no) : Index of equivalent orbital in unit cell 60C Unit cell orbitals must be the first in 61C orbital lists, i.e. indxuo.le.nuo, with 62C nuo the number of orbitals in unit cell 63C integer nk : Number of k points 64C real*8 kpoint(3,nk) : k point vectors 65C real*8 wk(nk) : k point weights (must sum one) 66C integer n_eigenvectors : Number of eigenvectors to be saved 67C integer nuotot : total number of orbitals per unit cell 68C over all processors 69C real*8 occtol : Occupancy threshold for DM build 70C *************************** OUTPUT ********************************** 71C real*8 eo(maxo,maxspn,nk) : Eigenvalues 72C ******************** OUTPUT (only if getD=.true.) ******************* 73C real*8 qo(maxo,maxspn,nk) : Occupations of eigenstates 74C real*8 Dnew(maxnd,nspin) : Output Density Matrix 75C real*8 Enew(maxnd,nspin) : Output Energy-Density Matrix 76C real*8 ef : Fermi energy 77C real*8 Entropy : Electronic entropy 78C *************************** AUXILIARY ******************************* 79C real*8 Haux(2,nuotot,nuo) : Auxiliary space for the hamiltonian matrix 80C real*8 Saux(2,nuotot,nuo) : Auxiliary space for the overlap matrix 81C real*8 psi(2,nuotot,nuo) : Auxiliary space for the eigenvectors 82C real*8 Dk(2,nuotot,nuo) : Aux. space that may be the same as Haux 83C real*8 Ek(2,nuotot,nuo) : Aux. space that may be the same as Saux 84C *************************** UNITS *********************************** 85C xij and kpoint must be in reciprocal coordinates of each other. 86C temp and H must be in the same energy units. 87C eo, Enew and ef returned in the units of H. 88C *************************** PARALLEL ******************************** 89C The auxiliary arrays are now no longer symmetric and so the order 90C of referencing has been changed in several places to reflect this. 91C ********************************************************************* 92C 93C Modules 94C 95 use precision 96 use sys 97 use parallel, only : Node, Nodes, BlockSize 98 use parallelsubs, only : LocalToGlobalOrb 99 use writewave, only : writew 100 use m_fermid, only : fermid, fermispin, stepf 101#ifdef MPI 102 use mpi_siesta 103#endif 104 105 use iowfs_netcdf, only : setup_wfs_netcdf_file, write_wfs_netcdf 106 use iowfs_netcdf, only : get_wfs_netcdf, get_wfs_block_netcdf 107 use iowfs_netcdf, only : open_wfs_netcdf_file, 108 $ close_wfs_netcdf_file 109 110 implicit none 111 112#ifdef MPI 113 integer 114 . MPIerror 115#endif 116 117 integer maxnd, maxnh, maxspn, maxo, nk, no, 118 . nspin, nuo, nuotot, iscf, n_eigenvectors 119 integer indxuo(no), listh(maxnh), numh(nuo), 120 . listd(maxnd), numd(nuo), listhptr(nuo), 121 . listdptr(nuo) 122 real(dp) Dnew(maxnd,nspin), 123 . e1, e2, ef, efs(nspin), Enew(maxnd,nspin), 124 . Entropy, eo(maxo,maxspn,nk), H(maxnh,nspin), 125 . kpoint(3,nk), qo(maxo,maxspn,nk), qtot, 126 . S(maxnh), temp, wk(nk), occtol, 127 . xij(3,maxnh), qs(nspin) 128 real(dp) Dk(2,nuotot,nuo), Ek(2,nuotot,nuo), 129 . Haux(2,nuotot,nuo), Saux(2,nuotot,nuo), 130 . psi(2,nuotot,nuo) 131 logical getD, getPSI, fixspin 132 external cdiag 133 134C Internal variables ............................................. 135 integer 136 . BNode, BTest, ie, ierror, iie, ik, ind, io, iio, 137 . ispin, iuo, j, jo, juo, nd, neigneeded 138 139 real(dp) 140 . ckxij, ee, kxij, pipj1, pipj2, qe, skxij, t 141 real(dp) :: qpipj1, qpipj2, epipj1, epipj2 142 143 integer :: ifirst, ilast, nvecs_in_block, nvecs_read, i 144 integer :: max_block_size 145 real(dp), dimension(:,:,:), allocatable, target :: psi_block 146 real(dp), dimension(:,:), pointer :: psi_aux 147C .................... 148 149C Find eigenvalues and eigenvectors 150#ifdef DEBUG 151 call write_debug( ' PRE diagk_file' ) 152#endif 153 154 ! Setup netCDF file 155 ! 156 if (Node == 0) then 157 call setup_wfs_netcdf_file(nuotot, nk, nspin, n_eigenvectors) 158 endif 159 do ik = 1,nk 160 do ispin = 1,nspin 161 call timer( 'c-eigval', 1 ) 162 163 call build_dense_HS(ik,ispin) 164 ! Compute all eigenvectors for now 165 call cdiag(Haux,Saux,nuotot,nuo,nuotot,eo(1,ispin,ik),psi, 166 . nuotot,iscf,ierror, BlockSize) 167 168 if (ierror.gt.0) then 169 call die('Terminating due to failed diagonalisation') 170 elseif (ierror.lt.0) then 171 ! Repeat diagonalisation with increased memory to handle clustering 172 173 call build_dense_HS(ik,ispin) 174 175 call cdiag(Haux,Saux,nuotot,nuo,nuotot,eo(1,ispin,ik),psi, 176 . nuotot,iscf,ierror, BlockSize) 177 endif 178 call timer( 'c-eigval', 2 ) 179! 180! SAVE eigenvectors to file (only up to those specified) 181! 182 call timer( 'c-save', 1 ) 183 call write_wfs_netcdf(nuotot,nuo,ik,ispin, psi, 184 $ n_eigenvectors,eo(1:n_eigenvectors,ispin,ik)) 185 call timer( 'c-save', 2 ) 186 187 enddo 188 enddo 189 190 call close_wfs_netcdf_file() 191 192C Find new Fermi energy and occupation weights ........................ 193 if (fixspin) then 194 call fermispin( nspin, nspin, nk, wk, maxo, nuotot, eo, 195 . temp, qs, qo, efs, Entropy ) 196 else 197 call fermid(nspin, maxspn, nk, wk, maxo, nuotot, eo, 198 . temp, qtot, qo, ef, Entropy ) 199 endif 200 201C Find weights for local density of states ............................ 202 if (e1 .lt. e2) then 203* e1 = e1 - ef 204* e2 = e2 - ef 205 t = max( temp, 1.e-6_dp ) 206 do ik = 1,nk 207 do ispin = 1,nspin 208 do io = 1,nuotot 209 qo(io,ispin,ik) = wk(ik) * 210 . ( stepf((eo(io,ispin,ik)-e2)/t) - 211 . stepf((eo(io,ispin,ik)-e1)/t)) * 2.0_dp/nspin 212 enddo 213 enddo 214 enddo 215 endif 216 217C New density and energy-density matrices of unit-cell orbitals ....... 218 nd = listdptr(nuo) + numd(nuo) 219 Dnew(1:nd,1:nspin) = 0.0_dp 220 Enew(1:nd,1:nspin) = 0.0_dp 221 222 call open_wfs_netcdf_file() 223 224 max_block_size = n_eigenvectors / Nodes 225 allocate(psi_block(2,nuotot,max_block_size)) 226 227C Loop over k points 228 do ik = 1,nk 229 do ispin = 1,nspin 230 231 neigneeded = 0 232 ie = nuotot 233 do while (ie.gt.0.and.neigneeded.eq.0) 234 qe = qo(ie,ispin,ik) 235 if (abs(qe).gt.occtol) neigneeded = ie 236 ie = ie - 1 237 enddo 238 if (neigneeded > n_eigenvectors) then 239 if (Node == 0) then 240 write(6,"(a)") "ERROR: More eigenvectors are needed!" 241 write(6,"(a,2i6)") "Saved, needed: ", n_eigenvectors, 242 $ neigneeded 243 endif 244 call die() 245 endif 246 247 call timer( 'c-buildD', 1 ) 248 249C Global operation to form new density matrix 250 251! Get the eigenvectors in blocks 252 nvecs_read = 0 253 do while (nvecs_read < neigneeded) 254 255 ifirst = nvecs_read + 1 256 ilast = ifirst + max_block_size - 1 257 ilast = min(ilast,neigneeded) 258 nvecs_in_block = ilast - ifirst + 1 259 260 ! Get eigenvectors (will broadcast to all nodes) 261 call get_wfs_block_netcdf(ifirst,nvecs_in_block, 262 $ ik,ispin,nuotot,psi_block) 263 nvecs_read = nvecs_read + nvecs_in_block 264 265 ! Loop over slots in density-matrix in this node 266 267 do iuo = 1,nuo 268 call LocalToGlobalOrb(iuo,Node,Nodes,iio) 269 do j = 1,numd(iuo) 270 ind = listdptr(iuo) + j 271 jo = listd(ind) 272 juo = indxuo(jo) 273 274 kxij = kpoint(1,ik) * xij(1,ind) + 275 . kpoint(2,ik) * xij(2,ind) + 276 . kpoint(3,ik) * xij(3,ind) 277 ckxij = cos(kxij) 278 skxij = sin(kxij) 279 280 qpipj1 = 0.0_dp 281 qpipj2 = 0.0_dp 282 epipj1 = 0.0_dp 283 epipj2 = 0.0_dp 284 do i = 1, nvecs_in_block 285 ie = ifirst + (i-1) 286 psi_aux => psi_block(:,:,i) 287 qe = qo(ie,ispin,ik) 288 ee = qo(ie,ispin,ik) * eo(ie,ispin,ik) 289 pipj1 = psi_aux(1,iio) * psi_aux(1,juo) + 290 . psi_aux(2,iio) * psi_aux(2,juo) 291 pipj2 = psi_aux(1,iio) * psi_aux(2,juo) - 292 . psi_aux(2,iio) * psi_aux(1,juo) 293 qpipj1 = qpipj1 + qe * pipj1 294 qpipj2 = qpipj2 + qe * pipj2 295 epipj1 = epipj1 + ee * pipj1 296 epipj2 = epipj2 + ee * pipj2 297 enddo 298 299 Dnew(ind,ispin)=Dnew(ind,ispin)+ qpipj1*ckxij - 300 . qpipj2*skxij 301 Enew(ind,ispin)=Enew(ind,ispin)+ epipj1*ckxij - 302 . epipj2*skxij 303 enddo 304 enddo 305 306 enddo ! while 307 308 309 call timer( 'c-buildD', 2 ) 310 311 enddo ! ispin 312 enddo ! ik 313 314 deallocate(psi_block) 315 call close_wfs_netcdf_file() 316 317#ifdef DEBUG 318 call write_debug( ' POS diagk_file' ) 319#endif 320 321 CONTAINS 322 323 subroutine build_dense_HS(ik,ispin) 324 integer, intent(in) :: ik, ispin 325 326 ! All other variables taken from parent scope 327 328 call timer( 'c-buildHS', 1 ) 329 Saux = 0.0_dp 330 Haux = 0.0_dp 331 do iuo = 1,nuo 332 do j = 1,numh(iuo) 333 ind = listhptr(iuo) + j 334 jo = listh(ind) 335 juo = indxuo(jo) 336 kxij = kpoint(1,ik) * xij(1,ind) + 337 . kpoint(2,ik) * xij(2,ind) + 338 . kpoint(3,ik) * xij(3,ind) 339 ckxij = cos(kxij) 340 skxij = sin(kxij) 341C Note : sign of complex part changed to match change in order of iuo/juo 342 Saux(1,juo,iuo) = Saux(1,juo,iuo) + S(ind)*ckxij 343 Saux(2,juo,iuo) = Saux(2,juo,iuo) - S(ind)*skxij 344 Haux(1,juo,iuo) = Haux(1,juo,iuo) + H(ind,ispin)*ckxij 345 Haux(2,juo,iuo) = Haux(2,juo,iuo) - H(ind,ispin)*skxij 346 enddo 347 enddo 348 call timer( 'c-buildHS', 2 ) 349 end subroutine build_dense_HS 350 351 end subroutine diagk_file 352#endif 353