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 pdos( NO, nspin, maxspn, NO_L, MAXNH, 9 . MAXO, NUMH, LISTHPTR, LISTH, H, S, 10 . E1, E2, SIGMA, NHIST, 11 . XIJ, INDXUO, GAMMA, NK, KPOINT, WK, EO, NO_U ) 12C ********************************************************************** 13C Subroutine to calculate the projected density of states on the 14C atomic orbitals for a given eigenvalue spectra 15C Written by J. Junquera and E. Artacho, November 1999. 16C *********** INPUT ************************************************** 17C INTEGER NO : Number of basis orbitals in the supercell 18C integer nspin=h_spin_dim : Number of spin components of H and D 19C integer maxspn=spinor_dim : Number of spin components of eo 20C INTEGER NO_L : Maximum number of atomic orbitals in the unit 21C cell. First dimension of eo, qo, last of xij 22C Must be at least max(indxuo) 23! IN THIS PROCESSOR 24C INTEGER MAXNH : Maximum number of orbitals interacting 25C with any orbital 26C INTEGER MAXO : First dimension of eo 27C INTEGER NUMH(NUO) : Number of nonzero elements of each row 28C of hamiltonian matrix 29C INTEGER LISTH(MAXNH) : Nonzero hamiltonian-matrix element 30C column indexes for each matrix row 31C INTEGER LISTHPTR(NUO) : Pointer to each row (-1) of the 32C density matrix 33C REAL*8 H(MAXNH,nspin) : Hamiltonian in sparse format 34C REAL*8 S(MAXNH) : Overlap in sparse format 35C REAL*8 E1, E2 : Energy range for density-matrix states 36C (to find local density of states) 37C Not used if e1 > e2 38C REAL*8 SIGMA : Width of the gaussian to expand the eigenvalues 39C INTEGER NHIST : Number of subdivisions of the histogram 40C REAL*8 XIJ(3,MAXNH) : Vectors between orbital centers (sparse) 41C (not used if only gamma point) 42C INTEGER INDXUO(NO) : Index of equivalent orbital in unit cell 43C Unit cell orbitals must be the first in 44C orbital lists, i.e. indxuo.le.nuo, with 45C nuo the nuber of orbitals in the unit cell 46C logical Gamma : whether only the Gamma point is sampled 47C INTEGER NK : Number of k points 48C REAL*8 KPOINT(3,NK) : k point vectors 49C REAL*8 WK(NK) : k point weights (must sum one) 50C REAL*8 EO(NO_L,maxspn,NK) : Eigenvalues 51C INTEGER NO_U : Total number of orbitals in unit cell 52C ********************************************************************** 53 54 use precision, only : dp 55 use parallel, only : Node, Nodes, IOnode 56 use fdf 57 use siesta_geom, only : xa, isa 58 use m_spin, only : spin 59 use atomlist, only : iphorb, iaorb 60 use atmfuncs, only : zetafio, mofio, lofio, cnfigfio, labelfis 61 use atmfuncs, only : pol, izofis 62 use flib_wxml, only : xmlf_t, xml_OpenFile, xml_NewElement, 63 . xml_AddArray, xml_AddPCData, 64 $ xml_AddAttribute, 65 . xml_EndElement, xml_Close 66 use xml, only : str, xml_dump_attribute 67 use densematrix, only : allocDenseMatrix, resetDenseMatrix 68 use densematrix, only : Haux, Saux, psi 69 use alloc, only : re_alloc, de_alloc 70 use units, only : eV 71 use files, only : slabel, label_length 72 use m_energies, only: Ef 73#ifdef MPI 74 use parallelsubs, only : GetNodeOrbs 75#endif 76 use m_diag_option, only: ParallelOverK, Serial 77 78 implicit none 79 80 integer 81 . NO, NSPIN, MAXSPN, NO_L, MAXNH, NK, NHIST, 82 . MAXO, NO_U 83 84 logical, intent(in) :: Gamma 85 integer 86 . NUMH(*), LISTH(MAXNH), LISTHPTR(*), INDXUO(NO) 87 88 real(dp) 89 . H(MAXNH,NSPIN), S(MAXNH), E1, E2, SIGMA, 90 . XIJ(3,MAXNH), KPOINT(3,NK), WK(NK), EO(MAXO,MAXSPN,NK) 91 92C Dynamic arrays ------------------------------------------------------- 93 real(dp), dimension(:,:) , pointer :: DTOT 94 real(dp), dimension(:,:,:), pointer :: DPR 95 96C Internal variables --------------------------------------------------- 97 integer 98 . nuo, nhs, npsi, iuo, ihist, ispin, 99 . iunit1, iunit2, i 100 101 integer iat, spec, ii, iorb 102 103 logical :: orig_ParallelOverK, orig_Serial 104 105 real(dp), dimension(:), pointer :: tmp 106 107 character*40 pos_string 108 109 character(len=label_length+5) :: fnamepro 110 character(len=label_length+4) :: fnametot 111 112 real(dp) delta, ener 113 114 external 115 . io_assign, io_close, 116 . pdosk, timer 117 118 type(xmlf_t) :: xf ! For new XML output 119 120#ifdef DEBUG 121 call write_debug( ' PRE pdos' ) 122#endif 123 124 call timer( 'pdos', 1) 125 126 orig_Serial = Serial 127 orig_ParallelOverK = ParallelOverK 128 129C Find the intervals between the subdivisions in the energy scale ------ 130 delta = (E2 - E1) / (NHIST-1) 131 132 ! Ensure we only use paralleloverk for 133 ! spin-(un)polarized calculations 134 if ( nspin > 2 ) ParallelOverK = .false. 135 ! Reset for Gamma-only PDOS calculations 136 if ( Gamma ) ParallelOverK = .false. 137 138 if ( Nodes > 1 .and. .not. ParallelOverK ) then 139 Serial = .false. 140 end if 141 142C Find number of orbitals per unit cell 143#ifdef MPI 144 call GetNodeOrbs(no_u,Node,Nodes,nuo) 145#else 146 nuo = no_u 147#endif 148 149C Check internal dimensions -------------------------------------------- 150 if ( nspin.le.2 .and. gamma) then 151 nhs = no_u * nuo 152 npsi = no_u * no_l * maxspn 153 elseif ( nspin.le.2 .and. .not.gamma) then 154 if (ParallelOverK) then 155 nhs = 2 * no_u * no_u 156 npsi = 2 * no_u * no_u 157 else 158 nhs = 2 * no_u * nuo 159 npsi = 2 * no_u * nuo 160 endif 161 elseif (nspin.ge.4) then 162 nhs = 2 * (2*no_u) * (2*nuo) 163 npsi = 2 * (2*no_u) * (2*nuo) 164 else 165 call die('diagon: ERROR: incorrect value of nspin') 166 endif 167 168 169C Allocate local arrays ------------------------------------------- 170C ----- 171 call allocDenseMatrix(nhs, nhs, npsi) 172 nullify( dtot, dpr ) 173 call re_alloc( dtot, 1, nhist, 1, spin%EDM, 'dtot', 'pdos' ) 174 call re_alloc( dpr, 1, nhist, 1, no_u, 1, spin%EDM, 175 & 'dpr', 'pdos' ) 176 177C Initialize the projected density of states --------------------------- 178 do ispin = 1, spin%EDM 179 do ihist = 1,nhist 180 dtot(ihist,ispin) = 0.d0 181 do iuo = 1,no_u 182 dpr(ihist,iuo,ispin) = 0.d0 183 enddo 184 enddo 185 enddo 186 187C Call appropiate routine ---------------------------------------------- 188 if (nspin.le.2 .and. gamma) then 189 call pdosg( nspin, NUO, NO, maxspn, MAXNH, 190 . MAXO, NUMH, LISTHPTR, LISTH, H, S, 191 . E1, E2, NHIST, SIGMA, INDXUO, EO, 192 . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) 193 elseif ( nspin.le.2 .and. .not.gamma) then 194 if (ParallelOverK) then 195 call pdoskp( nspin, NUO, NO, maxspn, MAXNH, 196 . MAXO, NUMH, LISTHPTR, LISTH, H, S, 197 . E1, E2, NHIST, SIGMA, 198 . XIJ, INDXUO, NK, KPOINT, WK, EO, 199 . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) 200 else 201 call pdosk( nspin, NUO, NO, maxspn, MAXNH, 202 . MAXO, NUMH, LISTHPTR, LISTH, H, S, 203 . E1, E2, NHIST, SIGMA, 204 . XIJ, INDXUO, NK, KPOINT, WK, EO, 205 . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) 206 endif 207 elseif (nspin == 4 .and. gamma) then 208 call pdos2g( NUO, NO, NO_L, MAXNH, 209 . MAXO, NUMH, LISTHPTR, LISTH, H, S, 210 . E1, E2, NHIST, SIGMA, INDXUO, EO, 211 . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) 212 elseif (nspin == 4 .and. .not. gamma) then 213 call pdos2k( NUO, NO, NO_L, MAXNH, 214 . MAXO, NUMH, LISTHPTR, LISTH, H, S, 215 . E1, E2, NHIST, SIGMA, 216 . XIJ, INDXUO, NK, KPOINT, WK, EO, 217 . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) 218 elseif (nspin == 8 .and. gamma) then 219 call pdos3g( NUO, NO, NO_L, MAXNH, 220 . MAXO, NUMH, LISTHPTR, LISTH, H, S, 221 . E1, E2, NHIST, SIGMA, INDXUO, EO, 222 . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) 223 elseif (nspin == 8 .and. .not. gamma) then 224 call pdos3k( NUO, NO, NO_L, MAXNH, 225 . MAXO, NUMH, LISTHPTR, LISTH, H, S, 226 . E1, E2, NHIST, SIGMA, 227 . XIJ, INDXUO, NK, KPOINT, WK, EO, 228 . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) 229 endif 230 231 232 233 234 if (IOnode) then 235C Open file for write on I/O node 236 fnametot = trim(slabel)//'.DOS' 237 call io_assign(iunit1) 238 open(unit=iunit1, file=fnametot, form='formatted', 239 . status='unknown') 240 241C Output histogram 242 select case ( spin%EDM ) 243 case ( 1 ) 244 do ihist = 1,nhist 245 ENER = E1 + (ihist-1) * delta 246 write(iunit1,'(2f20.5)') ener/ev,dtot(ihist,1)*eV 247 enddo 248 case ( 2 ) 249 do ihist = 1,nhist 250 ENER = E1 + (ihist-1) * delta 251 write(iunit1,'(3f20.5)') ener/ev,dtot(ihist,1)*eV, 252 . dtot(ihist,2)*eV 253 enddo 254 case default 255 do ihist = 1,nhist 256 ENER = E1 + (ihist-1) * delta 257 write(iunit1,'(5f20.5)') ener/ev,dtot(ihist,1)*eV, 258 . dtot(ihist,2)*eV,2.0_dp*dtot(ihist,3)*eV,2.0_dp*dtot(ihist,4)*eV 259 enddo 260 end select 261 262C Close file for write 263 call io_close(iunit1) 264 endif 265 266C New writing 267 if (IOnode) then 268 call xml_OpenFile(trim(slabel)//".PDOS.xml",xf,indent=.false.) 269 270 fnamepro = trim(slabel)//'.PDOS' 271 call io_assign(iunit2) 272 open(iunit2,file=fnamepro,form='formatted',status='unknown') 273 call xml_NewElement(xf,"pdos") 274 call xml_NewElement(xf,"nspin") 275 call xml_AddPCData(xf,spin%EDM) 276 call xml_EndElement(xf,"nspin") 277 278 write(iunit2,'(a)') '<pdos>' 279 write(iunit2,'(a,i1,a)') '<nspin>', spin%EDM, '</nspin>' 280 write(iunit2,'(a,i0,a)') '<norbitals>', no_u, '</norbitals>' 281 282 ! Write fermi-level 283 call xml_NewElement(xf,"fermi_energy") 284 call xml_AddAttribute(xf,"units","eV") 285 call xml_AddPCData(xf, Ef / eV) 286 call xml_EndElement(xf,"fermi_energy") 287 write(iunit2,'(a,e17.9,a)') '<fermi_energy units="eV">', 288 & Ef / eV, '</fermi_energy>' 289 290 ! Write sampled energies 291 write(iunit2,'(a)') '<energy_values units="eV">' 292 call xml_NewElement(xf,"energy_values") 293 call xml_AddAttribute(xf,"units","eV") 294 295 nullify( tmp ) 296 call re_alloc( tmp, 1, spin%EDM*nhist, 'tmp', 'pdos' ) 297 do ihist = 1,nhist 298 ENER = E1 + (ihist-1) * delta 299 tmp(ihist) = ener/eV 300 write(iunit2,'(e17.9)') tmp(ihist) 301 enddo 302 write(iunit2,'(a)') '</energy_values>' 303 304 call xml_AddArray(xf,tmp(1:nhist)) 305 call xml_EndElement(xf,"energy_values") 306 307 do i = 1, no_u 308 iat = iaorb(i) 309 iorb = iphorb(i) 310 spec = isa(iat) 311 312 call xml_NewElement(xf,"orbital") 313 write(iunit2,'(a)') '<orbital ' 314 call xml_AddAttribute(xf,"index",i) 315 call xml_dump_attribute(iunit2,"index",str(i)) 316 call xml_AddAttribute(xf,"atom_index",iat) 317 call xml_dump_attribute(iunit2,"atom_index",str(iat)) 318 call xml_AddAttribute(xf,"species",trim(labelfis(spec))) 319 call xml_dump_attribute(iunit2,"species",labelfis(spec)) 320 call xml_AddAttribute(xf,"Z",izofis(spec)) 321 call xml_dump_attribute(iunit2,"Z",str(izofis(spec))) 322 write(pos_string,'(3f11.6)') xa(1:3,iat) 323 call xml_AddAttribute(xf,"position",trim(pos_string)) 324 call xml_dump_attribute(iunit2,"position",pos_string) 325 call xml_AddAttribute(xf,"n",cnfigfio(spec,iorb)) 326 call xml_dump_attribute(iunit2,"n",str(cnfigfio(spec,iorb))) 327 call xml_AddAttribute(xf,"l",lofio(spec,iorb)) 328 call xml_dump_attribute(iunit2,"l",str(lofio(spec,iorb))) 329 call xml_AddAttribute(xf,"m",mofio(spec,iorb)) 330 call xml_dump_attribute(iunit2,"m",str(mofio(spec,iorb))) 331 call xml_AddAttribute(xf,"z",zetafio(spec,iorb)) 332 call xml_dump_attribute(iunit2,"z",str(zetafio(spec,iorb))) 333 call xml_AddAttribute(xf,"P",pol(spec,iorb)) 334 call xml_dump_attribute(iunit2,"P",str(pol(spec,iorb))) 335 write(iunit2,'(a)') '>' 336 337!------------------------------------------------------------ 338 call xml_NewElement(xf,"data") 339 write(iunit2,'(a)') '<data>' 340 select case ( spin%EDM ) 341 case ( 1 ) 342 do ihist = 1, nhist 343 tmp(ihist) = dpr(ihist,i,1) * eV 344 write(iunit2,'(e17.9)') tmp(ihist) 345 end do 346 call xml_AddArray(xf,tmp(1:nhist)) 347 case ( 2 ) 348 do ihist = 1, nhist 349 tmp(ihist*2-1) = dpr(ihist,i,1) * eV 350 tmp(ihist*2) = dpr(ihist,i,2) * eV 351 write(iunit2,'(2(tr1,e17.9))') tmp(ihist*2-1:ihist*2) 352 end do 353 call xml_AddArray(xf,tmp(1:2*nhist)) 354 case ( 4 ) 355 do ihist = 1, nhist 356 tmp(ihist*4-3) = dpr(ihist,i,1) * eV 357 tmp(ihist*4-2) = dpr(ihist,i,2) * eV 358 tmp(ihist*4-1) = 2 * dpr(ihist,i,3) * eV 359 tmp(ihist*4 ) = 2 * dpr(ihist,i,4) * eV 360 write(iunit2,'(4f10.5)') tmp(ihist*4-3:ihist*4) 361 end do 362 call xml_AddArray(xf,tmp(1:4*nhist)) 363 end select 364 call xml_EndElement(xf,"data") 365 write(iunit2,'(a)') '</data>' 366 call xml_EndElement(xf,"orbital") 367 write(iunit2,'(a)') '</orbital>' 368 369 end do 370 371C Close file 372 call xml_EndElement(xf,"pdos") 373 call xml_Close(xf) 374 write(iunit2,'(a)') '</pdos>' 375 call io_close(iunit2) 376 377C Free local workspace array 378 call de_alloc( tmp, 'tmp', 'pdos' ) 379 380 endif 381 382C Free local arrays ---------------------------------------------------- 383 call de_alloc( dpr, 'dpr', 'pdos' ) 384 call de_alloc( dtot, 'dtot', 'pdos' ) 385 call resetDenseMatrix() 386 387 ! Restore the paralleloverk options 388 ParallelOverK = orig_ParallelOverK 389 Serial = orig_Serial 390 391 call timer( 'pdos', 2) 392 393#ifdef DEBUG 394 call write_debug( ' POS pdos' ) 395#endif 396 397 end 398