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