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 detover(psiprev, psi, S, Sr,
9     .               numh, listhptr, listh, indxuo,
10     .               no, nuo, xij, maxnh, nuotot, nocc,
11     .               kpoint, dk, detr, deti )
12C *********************************************************************
13C Finds the determinant of the overlap matrix
14C between the periodic Bloch functions corresponding to neighboring
15C k points
16C Written by DSP. March 1999
17C Modified for parallel execution by J.D. Gale, March 2000
18C **************************** INPUT **********************************
19C real*8  psi(2,nuotot,nuo)   : Wavefunctions in current k point
20C real*8  S(maxnh)            : Overlap in sparse form
21C real*8  Sr(maxnh)           : Position operator matrix elements (sparse)
22C integer numh(nuo)           : Number of nonzero elements of each row
23C                               of hamiltonian matrix
24C integer listhptr(nuo)       : Pointer to the start of each row
25C                               of hamiltonian matrix
26C integer listh(maxnh)        : Nonzero hamiltonian-matrix element
27C                               column indexes for each matrix row
28C integer indxuo(no)          : Index of equivalent orbital in unit cell
29C                               Unit cell orbitals must be the first in
30C                               orbital lists, i.e. indxuo.le.nuo, with
31C                               nuo the number of orbitals in unit cell
32C integer no                  : Number of basis orbitals in supercell
33C integer nuo                 : Number of basis orbitals in the unit cell
34C integer maxnh               : Maximum number of orbitals interacting
35C integer nuotot              : Third dimension of xij
36C real*8  xij(3,maxnh)        : Vectors between orbital centers (sparse)
37C                               (not used if only gamma point)
38C integer nocc                : number of occupied states
39C real*8  kpoint(3)           : Current kpoint
40C real*8  dk(3)               : Vector joining the previous and current
41C                               kpoint
42C *************************** INPUT/OUTPUT ****************************
43C real*8  psiprev(2,nuo,nuotot) : Wavefunctions in previous k point
44C real*8  detr                  : Real part of the determinant
45C real*8  deti                  : Imaginary part of the determinant
46C *************************** UNITS ***********************************
47C Lengths in atomic units (Bohr).
48C k vectors in reciprocal atomic units.
49C Energies in Rydbergs.
50C *********************************************************************
51      use precision
52      use parallel,     only : Node, Nodes
53      use parallelsubs, only : GetNodeOrbs, LocalToGlobalOrb
54      use alloc,        only : re_alloc, de_alloc
55      use m_linpack,    only : zgedi, zgefa
56#ifdef MPI
57      use mpi_siesta
58#endif
59      implicit none
60
61      integer
62     .  nuo, maxnh, nuotot, no, nocc,
63     .  listh(maxnh), listhptr(nuo), numh(nuo), indxuo(no),
64     .  info, job
65      parameter(job=10)
66
67      real(dp)
68     .  psiprev(2,nuo,nuotot), dk(3),detr, deti,
69     .  xij(3,maxnh), S(maxnh), Sr(maxnh),
70     .  psi(2,nuotot,nuo), kpoint(3)
71
72      complex(dp) ::   det(2), pipj, determ
73
74C**** Internal variables ***********************************************
75
76      integer :: iuo, juo, j, ie, je, jo, ind
77      integer, dimension(:), pointer ::  Auxint
78#ifdef MPI
79      integer :: MPIerror, jeg, n, noccloc, noccmax
80      real(dp), dimension(:,:,:), pointer :: psitmp
81#endif
82      real(dp) ::  kxij, skxij, ckxij,pipj1, pipj2
83
84      complex(dp), dimension(:,:), pointer :: Aux, Aux2
85
86C Allocate local memory
87      nullify( Aux )
88      call re_alloc( Aux, 1, nocc, 1, nocc, name='Aux',
89     &               routine='detover' )
90      nullify( Aux2 )
91      call re_alloc( Aux2, 1, nuotot, 1, nuo, name='Aux2',
92     &               routine='detover' )
93      nullify( Auxint )
94      call re_alloc( Auxint, 1, nocc, name='Auxint', routine='detover' )
95
96      do iuo = 1,nuo
97        do j = 1,numh(iuo)
98          ind = listhptr(iuo)+j
99          jo = listh(ind)
100          juo = indxuo(jo)
101          kxij = (kpoint(1)-0.5d0*dk(1)) * xij(1,ind) +
102     .           (kpoint(2)-0.5d0*dk(2)) * xij(2,ind) +
103     .           (kpoint(3)-0.5d0*dk(3)) * xij(3,ind)
104          ckxij = dcos(kxij)
105          skxij = dsin(kxij)
106          Aux2(juo,iuo)=Aux2(juo,iuo)+
107     .      cmplx(  S(ind)*ckxij + Sr(ind)*skxij,
108     .      S(ind)*skxij - Sr(ind)*ckxij , dp )
109
110        enddo
111      enddo
112
113#ifdef MPI
114      call GetNodeOrbs(nuotot,0,Nodes,noccmax)
115      nullify( psitmp )
116      call re_alloc( psitmp, 1, 2, 1, nuotot, 1, noccmax, name='psitmp',
117     &               routine='detover' )
118
119C Ultimately this needs modifying so that Aux is distributed
120      do n = 0,Nodes-1
121
122C Broadcast copy of psi on node n to all other nodes
123        call GetNodeOrbs(nocc,n,Nodes,noccloc)
124        call GetNodeOrbs(nuotot,n,Nodes,noccmax)
125        if (Node .eq. n) then
126          do iuo = 1,noccmax
127            do juo = 1,nuotot
128              psitmp(1,juo,iuo) = psi(1,juo,iuo)
129              psitmp(2,juo,iuo) = psi(2,juo,iuo)
130            enddo
131          enddo
132        endif
133        call MPI_Bcast(psitmp(1,1,1),2*nuotot*noccmax,
134     .    MPI_double_precision,n,MPI_Comm_World,MPIerror)
135
136        do ie = 1,nocc
137          do je = 1, noccloc
138            call LocalToGlobalOrb(je,n,Nodes,jeg)
139
140            do iuo=1,nuo
141              do juo=1,nuotot
142
143                pipj1 = psiprev(1,iuo,ie) * psitmp(1,juo,je) +
144     .                  psiprev(2,iuo,ie) * psitmp(2,juo,je)
145                pipj2 = psiprev(1,iuo,ie) * psitmp(2,juo,je) -
146     .                  psiprev(2,iuo,ie) * psitmp(1,juo,je)
147                pipj=cmplx(pipj1,pipj2, dp)
148
149                Aux(jeg,ie) = Aux(jeg,ie) + pipj*Aux2(juo,iuo)
150
151              enddo
152            enddo
153
154          enddo
155        enddo
156
157      enddo
158
159      call de_alloc( psitmp, name='psitmp' , routine='detover')
160#else
161      do ie = 1,nocc
162        do je = 1, nocc
163
164          do iuo=1,nuo
165            do juo=1,nuotot
166
167              pipj1 = psiprev(1,iuo,ie) * psi(1,juo,je) +
168     .                psiprev(2,iuo,ie) * psi(2,juo,je)
169              pipj2 = psiprev(1,iuo,ie) * psi(2,juo,je) -
170     .                psiprev(2,iuo,ie) * psi(1,juo,je)
171              pipj=cmplx(pipj1,pipj2, dp)
172
173              Aux(je,ie) = Aux(je,ie) + pipj*Aux2(juo,iuo)
174
175            enddo
176          enddo
177
178        enddo
179      enddo
180#endif
181
182C Resize Aux2 for re-use
183      call re_alloc( Aux2, 1, nocc, 1, nocc, name='Aux2',
184     &               routine='detover', copy=.false. )
185
186#ifdef MPI
187      call MPI_AllReduce(Aux(1,1),Aux2(1,1),nocc*nocc,
188     .  MPI_double_complex,MPI_sum,MPI_Comm_World,MPIerror)
189      Aux = Aux2
190#endif
191
192      call zgefa(Aux,nocc,nocc,Auxint,info)
193      call zgedi(Aux,nocc,nocc,Auxint,det,Aux2,job)
194
195      determ=det(1)
196      deti=aimag(determ)
197      detr=real(determ)
198
199C Deallocate local memory
200      call de_alloc( Auxint, name='Auxint' )
201      call de_alloc( Aux2, name='Aux2' )
202      call de_alloc( Aux, name='Aux' )
203
204      end
205