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 module m_writedelk 9 public :: write_delk 10 CONTAINS 11 subroutine write_delk( gamma, no_u, no_s, nspin, indxuo, maxnh, 12 . numh, listhptr, listh, delkmat, S, qtot, temp, xij) 13C ********************************************************************* 14C Saves the hamiltonian and overlap matrices, and other data required 15C to obtain the bands and density of states 16C Writen by J.Soler July 1997. 17C Note because of the new more compact method of storing H and S 18C this routine is NOT backwards compatible 19C *************************** INPUT ********************************** 20C logical gamma : Is only gamma point used? 21C ******************** INPUT or OUTPUT (depending on task) *********** 22C integer no_u : Number of basis orbitals per unit cell 23C integer no_s : Number of basis orbitals per supercell 24C integer nspin : Spin polarization (1 or 2) 25C integer indxuo(no_s) : Index of orbitals in supercell 26C integer maxnh : First dimension of listh, delkmat, S and 27C second of xij 28C integer numh(nuo) : Number of nonzero elements of each row 29C of hamiltonian matrix 30C integer listhptr(nuo) : Pointer to the start of each row (-1) 31C of hamiltonian matrix 32C integer listh(maxnh) : Nonzero hamiltonian-matrix element column 33C indexes for each matrix row 34C real*8 delkmat(maxnh) : Hamiltonian in sparse form 35C real*8 S(maxnh) : Overlap in sparse form 36C real*8 qtot : Total number of electrons 37C real*8 temp : Electronic temperature for Fermi smearing 38C real*8 xij(3,maxnh) : Vectors between orbital centers (sparse) 39C (not read/written if only gamma point) 40 41C 42C Modules 43C 44 use precision, only: dp, sp 45 use parallel, only : Node, Nodes 46 use parallelsubs, only : WhichNodeOrb, LocalToGlobalOrb, 47 . GlobalToLocalOrb, GetNodeOrbs 48 use atm_types, only : nspecies 49 use atomlist, only : iphorb, iaorb 50 use siesta_geom, only : na_u, isa 51 use atmfuncs, only : nofis, labelfis, zvalfis 52 use atmfuncs, only : cnfigfio, lofio, zetafio 53 use fdf 54 use files, only : slabel, label_length 55 use sys, only : die 56#ifdef MPI 57 use mpi_siesta 58#endif 59 60 implicit none 61 62 logical gamma 63 integer maxnh, no_u, no_s, nspin 64 integer indxuo(no_s), listh(maxnh), numh(*), listhptr(*) 65 real(dp) S(maxnh), 66 . qtot, temp, xij(3,maxnh) 67 complex(dp) delkmat(maxnh) 68 69 external io_assign, io_close 70 71C Internal variables and arrays 72 character(len=label_length+4) :: fname 73 integer im, is, iu, ju, k, mnh, ns, ia, io 74 integer ih,hl,nuo,maxnhtot,maxhg 75 integer, dimension(:), allocatable :: numhg 76#ifdef MPI 77 integer MPIerror, Request, Status(MPI_Status_Size), BNode 78 integer, dimension(:), allocatable :: ibuffer 79 real(dp), dimension(:), allocatable :: buffer 80 complex(dp), dimension(:), allocatable :: buffer3 81 real(dp), dimension(:,:), allocatable :: buffer2 82#endif 83 logical baddim, gammaonfile, write_xijk 84 85C Find name of file 86 fname = trim(slabel) // '.dek' 87 88C Find total numbers over all Nodes 89#ifdef MPI 90 call MPI_AllReduce(maxnh,maxnhtot,1,MPI_integer,MPI_sum, 91 . MPI_Comm_World,MPIerror) 92#else 93 maxnhtot = maxnh 94#endif 95 96 if (Node.eq.0) then 97C Open file 98 call io_assign( iu ) 99 open( iu, file=fname, form='formatted', status='unknown' ) 100 101C Write overall data 102 write(iu,*) " no_u, no_s, nspin, maxnhtot = ", 103 . no_u, no_s, nspin, maxnhtot 104 105C Read logical 106 write(iu,*) " gamma = ", gamma 107 108C Allocate local array for global numh 109 allocate(numhg(no_u)) 110 call memory('A','I',no_u,'iohs') 111 112C Write out indxuo 113 if (.not.gamma) then 114 do ih = 1, no_s 115 write(iu,*)" ih, indxuo(ih) = ", ih, indxuo(ih) 116 enddo 117 endif 118 119 endif 120 121C Create globalised numh 122 do ih = 1,no_u 123#ifdef MPI 124 call WhichNodeOrb(ih,Nodes,BNode) 125 if (BNode.eq.0.and.Node.eq.BNode) then 126 call GlobalToLocalOrb(ih,Node,Nodes,hl) 127#else 128 hl = ih 129#endif 130 numhg(ih) = numh(hl) 131#ifdef MPI 132 elseif (Node.eq.BNode) then 133 call GlobalToLocalOrb(ih,Node,Nodes,hl) 134 call MPI_ISend(numh(hl),1,MPI_integer, 135 . 0,1,MPI_Comm_World,Request,MPIerror) 136 call MPI_Wait(Request,Status,MPIerror) 137 elseif (Node.eq.0) then 138 call MPI_IRecv(numhg(ih),1,MPI_integer, 139 . BNode,1,MPI_Comm_World,Request,MPIerror) 140 call MPI_Wait(Request,Status,MPIerror) 141 endif 142 if (BNode.ne.0) then 143 call MPI_Barrier(MPI_Comm_World,MPIerror) 144 endif 145#endif 146 enddo 147 148 if (Node.eq.0) then 149C Write numh 150 maxhg = 0 151 do ih = 1,no_u 152 maxhg = max(maxhg,numhg(ih)) 153 enddo 154 do ih = 1, no_u 155 write(iu,*)" ih, numhg(ih) = ", ih, numhg(ih) 156 enddo 157#ifdef MPI 158 allocate(buffer(maxhg)) 159 call memory('A','D',maxhg,'iohs') 160 allocate(buffer3(maxhg)) 161 call memory('A','D',maxhg,'iohs') 162 allocate(ibuffer(maxhg)) 163 call memory('A','I',maxhg,'iohs') 164#endif 165 endif 166 167C Write listh 168 do ih = 1,no_u 169#ifdef MPI 170 call WhichNodeOrb(ih,Nodes,BNode) 171 if (BNode.eq.0.and.Node.eq.BNode) then 172 call GlobalToLocalOrb(ih,Node,Nodes,hl) 173#else 174 hl = ih 175#endif 176 do im = 1, numh(hl) 177 write(iu,*) ' im, listh(listhptr(hl)+im) = ', 178 . im, listh(listhptr(hl)+im) 179 enddo 180#ifdef MPI 181 elseif (Node.eq.0) then 182 call MPI_IRecv(ibuffer,numhg(ih),MPI_integer,BNode,1, 183 . MPI_Comm_World,Request,MPIerror) 184 call MPI_Wait(Request,Status,MPIerror) 185 elseif (Node.eq.BNode) then 186 call GlobalToLocalOrb(ih,Node,Nodes,hl) 187 call MPI_ISend(listh(listhptr(hl)+1),numh(hl),MPI_integer, 188 . 0,1,MPI_Comm_World,Request,MPIerror) 189 call MPI_Wait(Request,Status,MPIerror) 190 endif 191 if (BNode.ne.0) then 192 call MPI_Barrier(MPI_Comm_World,MPIerror) 193 if (Node.eq.0) then 194 do im = 1, numhg(ih) 195 write(iu,*) " im, ibuffer(im) = ", 196 . im, ibuffer(im) 197 enddo 198 endif 199 endif 200#endif 201 enddo 202 203#ifdef MPI 204 if (Node.eq.0) then 205 call memory('D','I',size(ibuffer),'iohs') 206 deallocate(ibuffer) 207 endif 208#endif 209 210 211C Write delkmat 212 do ih=1,no_u 213#ifdef MPI 214 call WhichNodeOrb(ih,Nodes,BNode) 215 if (BNode.eq.0.and.Node.eq.BNode) then 216 call GlobalToLocalOrb(ih,Node,Nodes,hl) 217#else 218 hl = ih 219#endif 220 do im = 1, numh(hl) 221! write(iu,'(a,2i5,2f12.5)') 222! . ' im, delkmat(listhptr(hl)+im) = ', 223! . im, listhptr(hl)+im, delkmat(listhptr(hl)+im) 224 write(iu,'(2f12.5)') 225 . delkmat(listhptr(hl)+im) 226 enddo 227#ifdef MPI 228 elseif (Node.eq.0) then 229 call MPI_IRecv(buffer3,numhg(ih),MPI_double_complex, 230 . BNode,1,MPI_Comm_World,Request,MPIerror) 231 call MPI_Wait(Request,Status,MPIerror) 232 elseif (Node.eq.BNode) then 233 call GlobalToLocalOrb(ih,Node,Nodes,hl) 234 call MPI_ISend(delkmat(listhptr(hl)+1),numh(hl), 235 . MPI_double_complex,0,1,MPI_Comm_World, 236 . Request,MPIerror) 237 call MPI_Wait(Request,Status,MPIerror) 238 endif 239 if (BNode.ne.0) then 240 call MPI_Barrier(MPI_Comm_World,MPIerror) 241 if (Node.eq.0) then 242 do im = 1, numhg(ih) 243! write(iu,'(a,i5,2f12.5)') ' im, buffer3(im) = ', 244! . im, buffer3(im) 245 write(iu,'(2f12.5)') buffer3(im) 246 enddo 247 endif 248 endif 249#endif 250 enddo 251 252C Write Overlap matrix 253 do ih = 1,no_u 254#ifdef MPI 255 call WhichNodeOrb(ih,Nodes,BNode) 256 if (BNode.eq.0.and.Node.eq.BNode) then 257 call GlobalToLocalOrb(ih,Node,Nodes,hl) 258#else 259 hl = ih 260#endif 261 do im = 1, numh(hl) 262 write(iu,'(a,2i5,2f12.5)') ' im, S(listhptr(hl)+im) = ', 263 . im, listhptr(hl)+im, (real(S(listhptr(hl)+im),kind=sp)) 264 enddo 265#ifdef MPI 266 elseif (Node.eq.0) then 267 call MPI_IRecv(buffer,numhg(ih),MPI_double_precision, 268 . BNode,1,MPI_Comm_World,Request,MPIerror) 269 call MPI_Wait(Request,Status,MPIerror) 270 elseif (Node.eq.BNode) then 271 call GlobalToLocalOrb(ih,Node,Nodes,hl) 272 call MPI_ISend(S(listhptr(hl)+1),numh(hl), 273 . MPI_double_precision,0,1,MPI_Comm_World, 274 . Request,MPIerror) 275 call MPI_Wait(Request,Status,MPIerror) 276 endif 277 if (BNode.ne.0) then 278 call MPI_Barrier(MPI_Comm_World,MPIerror) 279 if (Node.eq.0) then 280 do im = 1, numhg(ih) 281 write(iu,'(a,i5,2f12.5)') 282 . ' im, (real(buffer(im),kind=sp) = ', 283 . im, (real(buffer(im),kind=sp)) 284 enddo 285 endif 286 endif 287#endif 288 enddo 289 290#ifdef MPI 291 if (Node .eq. 0) then 292C Free buffer array 293 call memory('D','D',size(buffer),'iohs') 294 deallocate(buffer) 295 call memory('D','D',size(buffer3),'iohs') 296 deallocate(buffer3) 297 endif 298#endif 299 300 if (Node.eq.0) then 301 write(iu,*) ' qtot,temp = ', qtot,temp 302 endif 303 304 write_xijk = .TRUE. 305 306 if (write_xijk) then 307#ifdef MPI 308C Allocate buffer array 309 if (Node .eq. 0) then 310 allocate(buffer2(3,maxhg)) 311 call memory('A','D',3*maxhg,'iohs') 312 endif 313#endif 314 do ih = 1,no_u 315#ifdef MPI 316 call WhichNodeOrb(ih,Nodes,BNode) 317 if (BNode.eq.0.and.Node.eq.BNode) then 318 call GlobalToLocalOrb(ih,Node,Nodes,hl) 319#else 320 hl = ih 321#endif 322 do im = 1, numh(hl) 323 write(iu,*) im, (real(xij(k,listhptr(hl)+im),kind=sp), 324 $ k=1,3) 325 enddo 326#ifdef MPI 327 elseif (Node.eq.0) then 328 call MPI_IRecv(buffer2(1,1),3*numhg(ih), 329 . MPI_double_precision,BNode,1,MPI_Comm_World, 330 . Request,MPIerror) 331 call MPI_Wait(Request,Status,MPIerror) 332 elseif (Node.eq.BNode) then 333 call GlobalToLocalOrb(ih,Node,Nodes,hl) 334 call MPI_ISend(xij(1,listhptr(hl)+1),3*numh(hl), 335 . MPI_double_precision,0,1,MPI_Comm_World, 336 . Request,MPIerror) 337 call MPI_Wait(Request,Status,MPIerror) 338 endif 339 if (BNode.ne.0) then 340 call MPI_Barrier(MPI_Comm_World,MPIerror) 341 if (Node.eq.0) then 342 do im = 1, numhg(ih) 343 write(iu,*) im, (real(buffer2(k,im),kind=sp), 344 $ k=1,3) 345 enddo 346 endif 347 endif 348#endif 349 enddo 350#ifdef MPI 351 if (Node .eq. 0) then 352C Free buffer array 353 call memory('D','D',size(buffer2),'iohs') 354 deallocate(buffer2) 355 endif 356#endif 357 endif ! write_xijk 358 359 if (Node.eq.0) then 360 361! 362! Write other useful info 363! 364 write(iu,*)' nspecies = ', nspecies 365 write(iu,*)(labelfis(is), zvalfis(is),nofis(is),is=1,nspecies) 366 do is = 1, nspecies 367 do io=1,nofis(is) 368 write(iu,*)cnfigfio(is,io),lofio(is,io),zetafio(is,io) 369 enddo 370 enddo 371 write(iu,*) na_u 372 write(iu,*) (isa(ia),ia=1,na_u) 373 write(iu,*) (iaorb(io), iphorb(io), io=1,no_u) 374 375C Deallocate local array for global numh 376 call memory('D','I',size(numhg),'iohs') 377 deallocate(numhg) 378C Close file 379 call io_close( iu ) 380 endif 381 382 end subroutine write_delk 383 end module m_writedelk 384