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