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