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 meshdscf
9C
10C Stores quantities that are connected with Dscf in mesh local
11C form when data is distributed for parallel execution
12C
13      use precision, only: dp
14      implicit none
15
16C ----------------------------------------------------------------------
17C Dscf related variables for parallel distributed form
18C ----------------------------------------------------------------------
19C integer listdl(sum(numdl))           : List of non-zero elements in
20C                                      : a row of DscfL
21C integer listdlptr(nrowsDscfL)        : Pointer to row in listdl
22C integer NeedDscfL(nuotot)            : Pointer as to whether a row of
23C                                      : Dscf is needed in DscfL
24C integer nrowsDscfL                   : Number of rows of DscfL
25C integer numdl(nrowsDscfL)            : Number of non-zero elements in
26C                                      : a row of DscfL
27C real(dp)  DscfL(maxndl,nrowsDscfL)     : Local copy of Dscf elements
28C                                      : needed for the local mesh
29C ----------------------------------------------------------------------
30      integer, public            :: nrowsDscfL
31      integer,  pointer, public  :: listdl(:), listdlptr(:),
32     $                              NeedDscfL(:), numdl(:)
33      real(dp), pointer, public  :: DscfL(:,:)
34      logical           :: first_time = .true.
35
36      public ::  matrixOtoM, matrixMtoO, matrixMtoOC
37      public ::  resetDscfPointers
38      public ::  CreateLocalDscfPointers
39
40      private
41
42      CONTAINS
43
44      subroutine resetDscfPointers( )
45      use alloc, only : de_alloc
46      use m_dscfcomm, only: resetdscfComm
47      implicit none
48      call resetdscfComm( )
49      call de_alloc( listdl,    'listdl',    'meshdscf' )
50      call de_alloc( listdlptr, 'listdlptr', 'meshdscf' )
51      call de_alloc( NeedDscfL, 'NeedDscfL', 'meshdscf' )
52      call de_alloc( numdl,     'numdl',     'meshdscf' )
53      call de_alloc( DscfL,     'DscfL',     'meshdscf' )
54      end subroutine resetDscfPointers
55
56      subroutine CreateLocalDscfPointers( nmpl, nuotot, numd,
57     &                                    listdptr, listd )
58C
59C Calculates the values of the orbitals at the mesh points
60C
61C Update: Computes the communications needed to move data ordered
62C by orbitals to data ordered by mesh (function dscfComm). All-to-all
63C communications has been substituted by point-to-point communications.
64C Written by Rogeli Grima (BSC) Dec.2007
65C
66C ----------------------------------------------------------------------
67C Input :
68C ----------------------------------------------------------------------
69C integer nmpl          : Number of mesh points in unit cell locally
70C integer nuotot        : Total number of basis orbitals in unit cell
71C integer numd(nuo)     : Number of nonzero density-matrix
72C                       : elements for each matrix row
73C integer listdptr(nuo) : Pointer to start of rows of density-matrix
74C integer listd(maxnh)  : Nonzero-density-matrix-element column
75C                       : indexes for each matrix row
76C ----------------------------------------------------------------------
77C Output :
78C ----------------------------------------------------------------------
79C All output quantities are in the module meshdscf
80C ----------------------------------------------------------------------
81
82C
83C Modules
84C
85      use atomlist,     only : indxuo
86      use meshphi,      only : endpht, lstpht
87      use parallel,     only : Node, Nodes
88      use parallelsubs, only : GlobalToLocalOrb, WhichNodeOrb
89      use precision
90      use alloc
91      use m_dscfComm,   only : dscfcomm
92      use m_dscfComm,   only : DCsize, DCself, DCtotal, DCmax
93      use m_dscfComm,   only : DCmaxnd
94      use m_dscfComm,   only : DCsrc, DCdst, DCsic, DCinvp
95
96
97#ifdef MPI
98      use mpi_siesta
99#endif
100      implicit none
101C
102C Passed arguments
103C
104      integer, intent(in) :: nmpl
105      integer, intent(in) :: nuotot
106      integer, intent(in) :: numd(*)
107      integer, intent(in) :: listdptr(*)
108      integer, intent(in) :: listd(*)
109C
110C Local variables
111C
112      integer          :: i, ii, io, iio, ip, imp, iu, numdele,
113     &                    maxndmax, nsize, nn, mm
114      integer, pointer :: ibuffer(:)
115#ifdef MPI
116      integer          :: MPIerror, Status(MPI_Status_Size)
117#endif
118
119#ifdef DEBUG
120      call write_debug( '      PRE CreateLocalDscfPointers' )
121#endif
122
123      if (first_time) then
124        first_time = .false.
125      else
126        call resetDscfPointers( )
127      endif
128
129      nullify( NeedDscfL, listdl, listdlptr, numdl )
130
131C     Create pointer as to whether a given row of DscfL is needed in NeedDscfL
132C     This pointer is never deallocated...
133      call re_alloc( NeedDscfL, 1, nuotot, 'NeedDscfL', 'meshdscf' )
134
135      do io= 1, nuotot
136        NeedDscfL(io) = 0
137      enddo
138
139      do ip = 1,nmpl
140        do imp = 1+endpht(ip-1), endpht(ip)
141          i = lstpht(imp)
142          iu = indxuo(i)
143          NeedDscfL(iu) = 1
144        enddo
145      enddo
146      nrowsDscfL = 0
147      do i = 1,nuotot
148        if (NeedDscfL(i).eq.1) then
149          nrowsDscfL = nrowsDscfL + 1
150          NeedDscfL(i) = nrowsDscfL
151        endif
152      enddo
153
154C     Computes the communications needed to move data ordered
155C     by orbitals to data ordered by mesh.
156      call dscfComm( nuotot, nrowsDscfL, NeedDscfL )
157
158C     Allocate/reallocate memory for numdl and listdlptr
159      call re_alloc( numdl, 1, max(1,nrowsDscfL), 'numdl', 'meshdscf' )
160      call re_alloc( listdlptr, 1, max(1,nrowsDscfL),
161     &               'listdlptr', 'meshdscf' )
162
163C     Distribute information about numd globally
164C     Use communications scheduling precomputed in dscfComm
165      nullify( ibuffer )
166      call re_alloc( ibuffer, 1, DCmax, 'ibuffer', 'meshdscf' )
167
168C     This data is in the current process. No communications needed
169      do ii= 1, DCself
170        io = DCinvp(ii)
171        call GlobalToLocalOrb(io,Node,Nodes,iio)
172        numdele              = numd(iio)
173        numdl(NeedDscfL(io)) = numdele
174      enddo
175
176      nn       = DCself + 1
177      maxndmax = 0
178#ifdef MPI
179      do i= 1, DCsize
180C       For all the needed communications
181        numdele = 0
182        if (DCsrc(i).eq.Node) then
183C         If this is the source node, copy data into a buffer
184C         and send it to the receiver.
185          do ii= 1, DCsic(i)
186            io = DCinvp(nn)
187            call GlobalToLocalOrb(io,Node,Nodes,iio)
188            ibuffer(ii) = numd(iio)
189            numdele = numdele + ibuffer(ii)
190            nn = nn + 1
191          enddo
192          call MPI_Send( ibuffer, DCsic(i), MPI_integer,
193     &                   DCdst(i), 1, MPI_Comm_World, MPIerror )
194        else
195C         If this is the destination node, receive data from source node.
196          call mpi_recv( ibuffer, DCsic(i), MPI_integer,
197     &                   DCsrc(i),  1, MPI_Comm_world, Status,
198     &                   MPIerror )
199          do ii= 1, DCsic(i)
200            io                   = DCinvp(nn)
201            numdl(NeedDscfL(io)) = ibuffer(ii)
202            numdele              = numdele + ibuffer(ii)
203            nn                   = nn + 1
204          enddo
205        endif
206        if (numdele.gt.maxndmax) maxndmax = numdele
207      enddo
208#endif
209      call de_alloc( ibuffer, 'ibuffer', 'meshdscf' )
210      DCmaxnd = maxndmax
211
212C     Create listdlptr using numdl
213      listdlptr(1) = 0
214      do io = 2,nrowsDscfL
215        listdlptr(io) = listdlptr(io-1) + numdl(io-1)
216      enddo
217
218C     Allocate/reallocate listdl
219      if (nrowsDscfL.gt.0) then
220        nsize = listdlptr(nrowsDscfL)+numdl(nrowsDscfL)
221      else
222        nsize = 1
223      endif
224      call re_alloc( listdl, 1, nsize, 'listdl', 'meshdscf' )
225
226C     Distribute information about listd globally
227C     Use communications scheduling precomputed in dscfComm
228      nullify( ibuffer )
229      call re_alloc( ibuffer, 1, DCmaxnd, 'ibuffer', 'meshdscf' )
230
231C     This data is in the current process. No communications needed
232      do ii= 1, DCself
233        io = DCinvp(ii)
234        call GlobalToLocalOrb(io,Node,Nodes,iio)
235        iu = NeedDscfL(io)
236        do i = 1,numd(iio)
237          listdl(listdlptr(iu)+i) = listd(listdptr(iio)+i)
238        enddo
239      enddo
240
241      nn = DCself + 1
242#ifdef MPI
243      do i= 1, DCsize
244C       For all the needed communications
245        if (DCsrc(i).eq.Node) then
246C         If this is the source node, copy data into a buffer
247C         and send it to the receiver.
248          mm = 0
249          do ii= 1, DCsic(i)
250            io = DCinvp(nn)
251            call GlobalToLocalOrb(io,Node,Nodes,iio)
252            ibuffer(mm+1:mm+numd(iio)) = listd(listdptr(iio)+1:
253     &                                         listdptr(iio)+numd(iio))
254            mm = mm + numd(iio)
255            nn = nn + 1
256          enddo
257          call MPI_Send( ibuffer, mm, MPI_integer, DCdst(i), 1,
258     &                   MPI_Comm_World, MPIerror )
259        else
260C         If this is the destination node, receive data from source node.
261          mm = 0
262          do ii= 0, DCsic(i)-1
263            io = DCinvp(nn+ii)
264            mm = mm + numdl(NeedDscfL(io))
265          enddo
266          call mpi_recv( ibuffer, mm, MPI_integer, DCsrc(i), 1,
267     &                   MPI_Comm_world, Status, MPIerror )
268          mm = 0
269          do ii= 1, DCsic(i)
270            io = DCinvp(nn)
271            iu = NeedDscfL(io)
272            listdl(listdlptr(iu)+1:listdlptr(iu)+numdl(iu)) =
273     &              ibuffer(mm+1:mm+numdl(iu))
274            mm = mm + numdl(iu)
275            nn = nn + 1
276          enddo
277        endif
278      enddo
279#endif
280      call de_alloc( ibuffer, 'ibuffer', 'meshdscf' )
281
282#ifdef DEBUG
283      call write_debug( '      POS CreateLocalDscfPointers' )
284#endif
285      end subroutine CreateLocalDscfPointers
286
287      subroutine matrixOtoM( maxnd, numd, listdptr, maxndl, nuo,
288     &                       nspin, Dscf, DscfL )
289C ********************************************************************
290C Transforms a matrix which is distributed by block cyclic
291C distribution of orbitals to a matrix that contains all
292C the orbital rows needed for a mesh point distribution
293C over the nodes.
294C Created by J.D.Gale, February 2000
295C
296C Update: All-to-all communications has been substituted by
297C point-to-point communications. These have been precomputed in
298C dscfComm.
299C Written by Rogeli Grima (BSC) Dec.2007
300C
301C *********************** INPUT **************************************
302C integer maxnd         : First dimension of Dscf
303C integer numd(nuo)     : Number of non-zero elements in row of Dscf
304C integer listdptr(nuo) : Pointer to start of rows in Dscf
305C integer maxndl        : First dimension of DscfL
306C integer nuo           : Local no. of orbitals in unit cell
307C integer nspin         : Number of spin components
308C real*8  Dscf(maxnd,nspin) : Matrix in orbital distributed form
309C *********************** OUTPUT *************************************
310C real*8  DscfL(maxndl,nspin) : Matrix in mesh distributed form
311C ********************************************************************
312
313C  Modules
314      use precision
315      use alloc,        only : re_alloc, de_alloc
316#ifdef MPI
317      use mpi_siesta
318      use parallel,     only : Node, Nodes
319      use parallelsubs, only : WhichNodeOrb, GlobalToLocalOrb
320      use m_dscfComm
321#endif
322      implicit none
323
324C Argument types and dimensions
325      integer           :: maxnd, maxndl, nspin, nuo, numd(nuo),
326     &                     listdptr(nuo)
327      real(dp)          :: Dscf(maxnd,nspin), DscfL(maxndl,nspin)
328      external             memory
329
330C Internal variables and arrays
331      integer           :: io, ispin
332#ifdef MPI
333      integer           :: i, ii, iu, mm, nn, iio, MPIerror,
334     &                     Status(MPI_Status_Size)
335      real(dp), pointer :: buffer(:)
336#else
337      integer           :: il
338#endif
339
340C***********************
341C  Parallel execution  *
342C***********************
343#ifdef MPI
344C Allocate local Dscf storage array
345      nullify(buffer)
346      call re_alloc( buffer, 1, DCmaxnd*nspin, 'buffer', 'meshdscf' )
347
348C     Copy the data that is in the current process
349      do ii= 1, DCself
350        io = DCinvp(ii)
351        call GlobalToLocalOrb(io,Node,Nodes,iio)
352        iu = NeedDscfL(io)
353        do ispin = 1,nspin
354          DscfL(listdlptr(iu)+1:listdlptr(iu)+numdl(iu),ispin) =
355     &    Dscf(listdptr(iio)+1:listdptr(iio)+numd(iio),ispin)
356        enddo
357      enddo
358      nn = DCself
359
360C     Send or receive the data from/to other processes
361      do i= 1, DCsize
362        if (DCsrc(i).eq.Node) then
363          mm = 0
364          do ispin = 1,nspin
365            do ii= 1, DCsic(i)
366              io = DCinvp(nn+ii)
367              call GlobalToLocalOrb(io,Node,Nodes,iio)
368              buffer(mm+1:mm+numd(iio)) = Dscf(listdptr(iio)+1:
369     &                           listdptr(iio)+numd(iio),ispin)
370              mm = mm + numd(iio)
371            enddo
372          enddo
373          call MPI_Send( buffer, mm, MPI_double_precision, DCdst(i), 1,
374     &                   MPI_Comm_World, MPIerror )
375        else
376          mm = 0
377          do ii= 1, DCsic(i)
378            io = DCinvp(nn+ii)
379            iu = NeedDscfL(io)
380            mm = mm + numdl(iu)
381          enddo
382          mm = mm*nspin
383          call mpi_recv( buffer, mm, MPI_double_precision, DCsrc(i), 1,
384     &                   MPI_Comm_world, Status, MPIerror )
385          mm = 0
386          do ispin = 1,nspin
387            do ii= 1, DCsic(i)
388              io = DCinvp(nn+ii)
389              iu = NeedDscfL(io)
390              DscfL(listdlptr(iu)+1:listdlptr(iu)+numdl(iu),ispin) =
391     &              buffer(mm+1:mm+numdl(iu))
392              mm = mm + numdl(iu)
393            enddo
394          enddo
395        endif
396        nn = nn + DCsic(i)
397      enddo
398      call de_alloc( buffer, 'buffer', 'meshdscf' )
399#else
400C*********************
401C  Serial execution  *
402C*********************
403C Loop over rows of Dscf checking to see if they are in DscfL
404      do ispin = 1,nspin
405        do io = 1,nuo
406
407C Get pointer for this row of Dscf and see if it is needed for DscfL
408          il = NeedDscfL(io)
409          if (il.gt.0) then
410            DscfL(listdlptr(il)+1:listdlptr(il)+numdl(il),ispin) =
411     &        Dscf(listdptr(io)+1:listdptr(io)+numdl(il),ispin)
412          endif
413
414        enddo
415      enddo
416#endif
417      end subroutine matrixOtoM
418
419      subroutine matrixMtoO( maxnvl, maxnv, numVs, listVsptr, nuo,
420     &                       nspin, VsL, Vs )
421C ********************************************************************
422C Transforms a matrix which is distributed by mesh points to a matrix
423C that is distributed by a block cyclic distribution over the orbitals
424C and adds it to an existing array of this form.
425C Created by J.D.Gale, February 2000
426C
427C Update: All-to-all communications has been substituted by
428C point-to-point communications. These have been precomputed in
429C dscfComm.
430C Written by Rogeli Grima (BSC) Dec.2007
431C
432C *********************** INPUT **************************************
433C integer maxnvl          : First dimension of VsL and maximum number
434C                           of nonzero elements in VsL
435C integer maxnv           : First dimension of Vs and maximum number
436C                           of nonzero elements in Vs
437C integer numVs(nuo)      : Number of non-zero elements in row of Vs
438C integer listVsptr(nuo)  : Pointer to start of rows in Vs
439C integer nuo             : Local no. of orbitals in unit cell
440C integer nspin           : Number of spin components
441C real*8  VsL(maxnvl,nspin) : Mesh contribution to be added to Vs
442C ******************** INPUT AND OUTPUT *******************************
443C real*8  Vs(maxnv,nspin) : Value of nonzero elements in each row of Vs
444C                           to which the potential matrix elements are
445C                           summed up
446C *********************************************************************
447
448C  Modules
449      use precision
450      use alloc, only: re_alloc, de_alloc
451#ifdef MPI
452      use mpi_siesta
453      use parallel,     only : Node, Nodes
454      use parallelsubs, only : WhichNodeOrb, GlobalToLocalOrb,
455     &                         LocalToGlobalOrb
456      use m_dscfComm
457#endif
458
459      implicit none
460
461C Argument types and dimensions
462      integer
463     &   maxnv, maxnvl, nspin, nuo, numVs(nuo),
464     &   listVsptr(nuo)
465      real(dp)
466     &   Vs(maxnv,nspin), VsL(maxnvl,nspin)
467
468C Internal variables and arrays
469      integer
470     &  i, iu, ispin
471
472#ifdef MPI
473      integer           :: ii, MPIerror, Status(MPI_Status_Size), nn,
474     &                     mm, io, iio
475      real(dp), pointer :: buffer(:)
476#endif
477
478C***********************
479C  Parallel execution  *
480C***********************
481#ifdef MPI
482C     Allocate a buffer to Send/Receive data
483      nullify(buffer)
484      call re_alloc( buffer, 1, DCmaxnd*nspin, 'buffer', 'meshdscf' )
485
486C     Copy the data that is in the current process from VsL to Vs
487      do ii= 1, DCself
488        io = DCinvp(ii)
489        call GlobalToLocalOrb(io,Node,Nodes,iio)
490        iu = NeedDscfL(io)
491        do ispin = 1,nspin
492          Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio),ispin) =
493     &    Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio),ispin) +
494     &    VsL(listdlptr(iu)+1:listdlptr(iu)+numdl(iu),ispin)
495        enddo
496      enddo
497      nn = DCself
498
499C     Send or receive the data from/to other processes
500      do i= 1, DCsize
501C       The communication is from dst to src
502        if (DCsrc(i).eq.Node) then
503          mm = 0
504          do ii= 1, DCsic(i)
505            io = DCinvp(nn+ii)
506            call GlobalToLocalOrb(io,Node,Nodes,iio)
507            mm = mm + numVs(iio)
508          enddo
509          mm = mm*nspin
510
511          call mpi_recv( buffer, mm, MPI_double_precision, DCdst(i), 1,
512     &                   MPI_Comm_world, Status, MPIerror )
513          mm = 0
514          do ispin = 1,nspin
515            do ii= 1, DCsic(i)
516              io = DCinvp(nn+ii)
517              call GlobalToLocalOrb(io,Node,Nodes,iio)
518              Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio),ispin) =
519     &        Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio),ispin) +
520     &        buffer(mm+1:mm+numVs(iio))
521              mm = mm + numVs(iio)
522            enddo
523          enddo
524        else
525          mm = 0
526          do ispin = 1,nspin
527            do ii= 1, DCsic(i)
528              io = DCinvp(nn+ii)
529              iio = NeedDscfL(io)
530              buffer(mm+1:mm+numdl(iio)) =
531     &        VsL(listdlptr(iio)+1:listdlptr(iio)+numdl(iio),ispin)
532              mm = mm + numdl(iio)
533            enddo
534          enddo
535          call MPI_Send( buffer, mm, MPI_double_precision, DCsrc(i), 1,
536     &                   MPI_Comm_World, MPIerror )
537        endif
538        nn = nn + DCsic(i)
539      enddo
540      call de_alloc( buffer, 'buffer', 'meshdscf' )
541#else
542C*********************
543C  Serial execution  *
544C*********************
545C Add those elements that are needed locally to the values already
546C stored in the orbital oriented array
547      do ispin = 1,nspin
548        do i = 1,nuo
549          iu = NeedDscfL(i)
550          if (iu.gt.0) then
551            Vs(listVsptr(i)+1:listVsptr(i)+numVs(i),ispin) =
552     &        Vs(listVsptr(i)+1:listVsptr(i)+numVs(i),ispin) +
553     &        VsL(listdlptr(iu)+1:listdlptr(iu)+numVs(i),ispin)
554          endif
555        enddo
556      enddo
557#endif
558      end subroutine matrixMtoO
559
560      subroutine matrixMtoOC( maxnvl, maxnv, numVs, listVsptr, nuo,
561     &                        VsL, Vs )
562C ********************************************************************
563C Transforms a matrix which is distributed by mesh points to a matrix
564C that is distributed by a block cyclic distribution over the orbitals
565C and adds it to an existing array of this form.
566C Created by J.D.Gale, February 2000
567C
568C Update: All-to-all communications has been substituted by
569C point-to-point communications. These have been precomputed in
570C dscfComm.
571C Written by Rogeli Grima (BSC) Dec.2007
572C Generalized by Javier Junquera in May. 2012 to the case where the
573C potential is a complex variable.
574C
575C *********************** INPUT **************************************
576C integer maxnvl          : First dimension of VsL and maximum number
577C                           of nonzero elements in VsL
578C integer maxnv           : First dimension of Vs and maximum number
579C                           of nonzero elements in Vs
580C integer numVs(nuo)      : Number of non-zero elements in row of Vs
581C integer listVsptr(nuo)  : Pointer to start of rows in Vs
582C integer nuo             : Local no. of orbitals in unit cell
583C real*8  VsL(maxnvl,2)   : Mesh contribution to be added to Vs
584C ******************** INPUT AND OUTPUT *******************************
585C complex*8  Vs(maxnv)    : Value of nonzero elements in each row of Vs
586C                           to which the matrix elements are
587C                           summed up
588C *********************************************************************
589
590C  Modules
591      use precision
592      use alloc, only: re_alloc, de_alloc
593#ifdef MPI
594      use mpi_siesta
595      use parallel,     only : Node, Nodes
596      use parallelsubs, only : WhichNodeOrb, GlobalToLocalOrb,
597     &                         LocalToGlobalOrb
598      use m_dscfComm
599#endif
600
601      implicit none
602
603C Argument types and dimensions
604      integer
605     &   maxnv, maxnvl, nuo, numVs(nuo),
606     &   listVsptr(nuo)
607      real(dp)
608     &   VsL(maxnvl,2)
609      complex(dp)
610     &   Vs(maxnv)
611
612C Internal variables and arrays
613      integer
614     &  i, iu, irelim
615
616#ifdef MPI
617      integer           :: ii, MPIerror, Status(MPI_Status_Size), nn,
618     &                     mm, io, iio
619      real(dp), pointer :: buffer(:)
620#endif
621
622C***********************
623C  Parallel execution  *
624C***********************
625#ifdef MPI
626C     Allocate a buffer to Send/Receive data
627      nullify(buffer)
628      call re_alloc( buffer, 1, DCmaxnd*2, 'buffer', 'meshdscf' )
629
630C     Copy the data that is in the current process from VsL to Vs
631      do ii= 1, DCself
632        io = DCinvp(ii)
633        call GlobalToLocalOrb(io,Node,Nodes,iio)
634        iu = NeedDscfL(io)
635        Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio)) =
636     &    Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio)) +
637     &    cmplx(VsL(listdlptr(iu)+1:listdlptr(iu)+numdl(iu),1),
638     &          VsL(listdlptr(iu)+1:listdlptr(iu)+numdl(iu),2), kind=dp)
639      enddo
640      nn = DCself
641
642C     Send or receive the data from/to other processes
643      do i= 1, DCsize
644C       The communication is from dst to src
645        if (DCsrc(i).eq.Node) then
646          mm = 0
647          do ii= 1, DCsic(i)
648            io = DCinvp(nn+ii)
649            call GlobalToLocalOrb(io,Node,Nodes,iio)
650            mm = mm + numVs(iio)
651          enddo
652          mm = mm*2
653
654          call mpi_recv( buffer, mm, MPI_double_precision, DCdst(i), 1,
655     &                   MPI_Comm_world, Status, MPIerror )
656          mm = 0
657          do irelim = 1, 2
658            do ii= 1, DCsic(i)
659              io = DCinvp(nn+ii)
660              call GlobalToLocalOrb(io,Node,Nodes,iio)
661              if( irelim .eq. 1) then
662                Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio)) =
663     &          Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio)) +
664     &          buffer(mm+1:mm+numVs(iio)) * cmplx(1.0,0.0,kind=dp)
665              elseif( irelim .eq. 2) then
666                Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio)) =
667     &          Vs(listVsptr(iio)+1:listVsptr(iio)+numVs(iio)) +
668     &          buffer(mm+1:mm+numVs(iio)) * cmplx(0.0,1.0,kind=dp)
669              endif
670              mm = mm + numVs(iio)
671            enddo
672          enddo
673        else
674          mm = 0
675          do irelim = 1, 2
676            do ii= 1, DCsic(i)
677              io = DCinvp(nn+ii)
678              iio = NeedDscfL(io)
679              buffer(mm+1:mm+numdl(iio)) =
680     &        VsL(listdlptr(iio)+1:listdlptr(iio)+numdl(iio),irelim)
681              mm = mm + numdl(iio)
682            enddo
683          enddo
684          call MPI_Send( buffer, mm, MPI_double_precision, DCsrc(i), 1,
685     &                   MPI_Comm_World, MPIerror )
686        endif
687        nn = nn + DCsic(i)
688      enddo
689      call de_alloc( buffer, 'buffer', 'meshdscf' )
690#else
691C*********************
692C  Serial execution  *
693C*********************
694C Add those elements that are needed locally to the values already
695C stored in the orbital oriented array
696      do i = 1,nuo
697        iu = NeedDscfL(i)
698        if (iu.gt.0) then
699          Vs(listVsptr(i)+1:listVsptr(i)+numVs(i)) =
700     &      Vs(listVsptr(i)+1:listVsptr(i)+numVs(i)) +
701     &      cmplx(VsL(listdlptr(iu)+1:listdlptr(iu)+numVs(i),1),
702     &            VsL(listdlptr(iu)+1:listdlptr(iu)+numVs(i),2),kind=dp)
703        endif
704      enddo
705#endif
706      end subroutine matrixMtoOC
707
708      end module meshdscf
709
710