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 reordpsi(psisave, psi, nuo, nuotot, nocc, ninc)
9C *********************************************************************
10C Reorders the coefficient of the bands, so the bands whose bands index
11C range between 1 and numincbands are the bands included for
12C wannierization.
13C Written by J. Junquera, October 2013, based on subroutine savepsi,
14C by J.D. Gale, March 2000
15C **************************** INPUT **********************************
16C real*8  psi(2,nuotot,nuo)      : Wavefunctions in current k point
17C integer nuo                    : Number of (local) orbitals in the unit cell
18C integer nuotot                 : Number of orbitals in the unit cell
19C integer nocc                   : Number of bands considered for
20C                                  Wannierization (before excluding bands)
21C integer ninc                   : Number of maximum bands whose
22C                                  coefficients will be stored in a node
23C *************************** OUTPUT **********************************
24C complex(dp)  psisave(nuotot,ninc) : Wavefunctions saved in the matrix of
25C                                  coefficients
26C                                  In psi (and last two in psisave):
27C                                  First  index: real or complex value
28C                                  Second index: atomic orbital
29C                                  Third  index: band index
30C *************************** UNITS ***********************************
31C Lengths in atomic units (Bohr).
32C k vectors in reciprocal atomic units.
33C Energies in Rydbergs.
34C *********************************************************************
35      use precision, only: dp
36      use m_siesta2wannier90, only : isexcluded
37#ifdef MPI
38      use alloc,              only : re_alloc, de_alloc
39      use parallel,           only : Node, Nodes, BlockSize
40      use parallelsubs,       only : GetNodeOrbs, LocalToGlobalOrb
41      use parallelsubs,       only : GlobalToLocalOrb
42      use mpi_siesta
43      use m_orderbands,       only : index_included_band
44      use m_orderbands,       only : node_included_band
45      use m_orderbands,       only : band_index_in_node
46      use m_orderbands,       only : which_band_in_node
47#endif
48
49      implicit none
50
51      integer, intent(in) ::  nuo, nuotot, nocc, ninc
52      real(dp), intent(in) ::   psi(2,nuotot,nuo)
53      complex(dp), intent(out) :: psisave(nuotot,ninc)
54
55
56
57C**** Internal variables ***********************************************
58
59      integer :: jband, jbandloc, iuo, juo, iband
60
61#ifdef MPI
62      integer
63     .  MPIerror, iuog, juog, n, noccloc
64      real(dp), dimension(:,:,:), pointer ::  psitmp
65#endif
66
67#ifdef MPI
68! AG
69! Allocate as in detover, using the number of orbitals on the first node,
70! as some of the nodes might have zero orbitals!
71!
72! Actually, this works because the first node handles the largest number
73! of orbitals when using the block-cyclic distribution.
74
75      call GetNodeOrbs(nocc,0,Nodes,noccloc)
76
77      nullify( psitmp )
78      call re_alloc( psitmp, 1, 2, 1, nuotot, 1, noccloc,
79     &                 name='psitmp', routine='reordpsi' )
80
81      jband = 0
82      do n = 0, Nodes-1
83
84!       Broadcast copy of psi on node n to all other nodes
85!       Note how not all the bands are broadcasted.
86!       Only the "occupied" bands in the corresponding node.
87        call GetNodeOrbs(nocc,n,Nodes,noccloc)
88
89        if (Node .eq. n) then
90           psitmp(1:2,1:nuotot,1:noccloc) = psi(1:2,1:nuotot,1:noccloc)
91        endif
92        call MPI_Bcast(psitmp(1,1,1),2*nuotot*noccloc,
93     .    MPI_double_precision,n,MPI_Comm_World,MPIerror)
94
95!       Save local part of psisave
96        do iuo = 1, noccloc
97          call LocalToGlobalOrb(iuo,n,Nodes,iuog)
98
99!         Select here if the local occupied band is included or not.
100!         If it is included, then a new index (jband) that ranges
101!         between 1 and numincbands(ispin) is assigned
102!
103!         This is the same logic as in m_orderbands, so the arrays
104!         can be re-used
105
106          if( .not.  isexcluded(iuog) ) then
107            jband    = jband + 1
108!           The coefficients of the included band will be stored
109!           in Node (node_included_band(jband)),
110!           and will be the band_index_in_node(jband)-th band
111!           stored on that Node
112            jbandloc = band_index_in_node(jband)
113            if( node_included_band(jband) .eq. Node ) then
114              do juo = 1, nuotot
115                psisave(juo,jbandloc) =
116     $             cmplx(psitmp(1,juo,iuo),psitmp(2,juo,iuo),kind=dp)
117              enddo ! atomic orbitals
118            endif   ! If the local Node is the same where the included
119                    !    band should be stored
120          endif     ! Is the band included for wannierization?
121        enddo       ! local bands on a given node
122      enddo         ! nodes
123      call de_alloc( psitmp,  name='psitmp' , routine='reordpsi')
124#else
125!     Straight serial copy
126      jband = 0
127      do iband = 1, nocc
128        if (.not. isexcluded(iband)) then
129          jband = jband + 1
130          psisave(:,jband) =
131     $        cmplx(psi(1,:,iband),psi(2,:,iband),kind=dp)
132        endif
133      enddo
134#endif
135
136
137      end subroutine reordpsi
138
139