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!
8subroutine diag2kspiral( nuo, no, maxnh, maxnd, maxo, &
9    numh, listhptr, listh, numd, listdptr, &
10    listd, H, S, getD, qtot, temp, e1, e2, &
11    xij, indxuo, nk, kpoint, wk, &
12    eo, qo, Dnew, Enew, ef, Entropy, qw, &
13    Haux, Saux, psi, Dk, Ek, aux, nuotot, &
14    occtol, iscf, neigwanted)
15! *********************************************************************
16! Calculates the eigenvalues and eigenvectors, density
17! and energy-density matrices, and occupation weights of each
18! eigenvector, for given Hamiltonian and Overlap matrices.
19! This version is for non-colinear spin with k-sampling and spiral
20! arrangement of spins.
21! Written by V. M. Garcia-Suarez. June 2002
22! Modified to reduce eigenvector computation by J.D. Gale Nov. 2004
23! Lots of bugs fixed and changed to complex by Nick Papior, June 2020
24! **************************** INPUT **********************************
25! integer nuo                 : Number of basis orbitals in unit cell
26! integer no                  : Number of basis orbitals in supercell
27! integer maxnh               : Maximum number of orbitals interacting
28! integer maxnd               : First dimension of listd / DM
29! integer maxo                : First dimension of eo and qo
30! integer numh(nuo)           : Number of nonzero elements of each row
31!                               of hamiltonian matrix
32! integer listhptr(nuo)       : Pointer to each row (-1) of the
33!                               hamiltonian matrix
34! integer listh(maxnh)        : Nonzero hamiltonian-matrix element
35!                               column indexes for each matrix row
36! integer numd(nuo)           : Number of nonzero elements of each row
37!                               of density matrix
38! integer listdptr(nuo)       : Pointer to each row (-1) of the
39!                               density matrix
40! integer listd(maxnd)        : Nonzero density-matrix element column
41!                               indexes for each matrix row
42! real*8  H(maxnh,4)          : Hamiltonian in sparse form
43! real*8  S(maxnh)            : Overlap in sparse form
44! logical getD                : Find occupations and density matrices?
45! real*8  qtot                : Number of electrons in unit cell
46! real*8  temp                : Electronic temperature
47! real*8  e1, e2              : Energy range for density-matrix states
48!                               (to find local density of states)
49!                               Not used if e1 > e2
50! real*8  xij(3,maxnh)        : Vectors between orbital centers (sparse)
51!                               (not used if only gamma point)
52! integer indxuo(no)          : Index of equivalent orbital in unit cell
53!                               Unit cell orbitals must be the first in
54!                               orbital lists, i.e. indxuo.le.nuo, with
55!                               nuo the number of orbitals in unit cell
56! real*8 qw(3)                : Wave vector for spiral configuration
57! integer nk                  : Number of k points
58! real*8  kpoint(3,nk)        : k point vectors
59! real*8  wk(nk)              : k point weights (must sum one)
60! integer nuotot              : total number of orbitals per unit cell
61!                               over all processors
62! real*8  occtol              : Occupancy threshold for DM build
63! integer iscf                : SCF cycle number
64! integer neigwanted          : Number of eigenvalues wanted
65! *************************** OUTPUT **********************************
66! real*8 eo(maxo*2,nk)        : Eigenvalues
67! real*8 qo(maxo*2,nk)        : Occupations of eigenstates
68! real*8 Dnew(maxnd,4)        : Output Density Matrix
69! real*8 Enew(maxnd,4)        : Output Energy-Density Matrix
70! real*8 ef                   : Fermi energy
71! real*8 Entropy              : Electronic entropy
72! *************************** AUXILIARY *******************************
73! complex*16 Haux(2,nuotot,2,nuo) : Aux. space for the hamiltonian matrix
74! complex*16 Saux(2,nuotot,2,nuo) : Aux. space for the overlap matrix
75! complex*16 psi(2,nuotot,2*nuo)  : Aux. space for the eigenvectors
76! complex*16 aux(2,nuotot)        : Extra auxiliary space
77! complex*16 Dk(2,nuotot,2,nuo)   : Aux. space that may be the same as Haux
78! complex*16 Ek(2,nuotot,2,nuo)   : Aux. space that may be the same as Saux
79! *************************** UNITS ***********************************
80! xij and kpoint must be in reciprocal coordinates of each other.
81! temp and H must be in the same energy units.
82! eo, Enew and ef returned in the units of H.
83! *************************** PARALLEL ********************************
84! The auxiliary arrays are now no longer symmetry and so the order
85! of referencing has been changed in several places to reflect this.
86! *********************************************************************
87!
88!  Modules
89!
90  use precision
91  use sys
92  use parallel,     only : Node, Nodes, BlockSize
93  use parallelsubs, only : LocalToGlobalOrb
94  use m_fermid,     only : fermid, stepf
95#ifdef MPI
96  use mpi_siesta
97#endif
98
99  implicit none
100
101#ifdef MPI
102  integer :: MPIerror
103#endif
104
105  integer :: maxnd, maxnh, maxo, nk, no, nuo, nuotot, iscf, neigwanted
106
107  integer :: indxuo(no), listh(maxnh), numh(nuo), listd(maxnd), numd(nuo), &
108      listhptr(nuo), listdptr(nuo)
109
110  real(dp) :: Dnew(maxnd,4), &
111      e1, e2, ef, Enew(maxnd,4), Entropy, eo(maxo*2,nk), &
112      H(maxnh,4), kpoint(3,nk), qo(maxo*2,nk), qtot, &
113      S(maxnh), temp, wk(nk), xij(3,maxnh), qw(3), occtol
114
115  complex(dp), dimension(2,nuotot,2,nuo) :: Haux, Saux, Dk, Ek
116  complex(dp), dimension(2,nuotot,2*nuo) :: psi
117  complex(dp), dimension(2,nuotot) :: aux
118
119  logical :: getD
120
121  external :: cdiag
122
123!  Internal variables .............................................
124  integer :: BNode, BTest, ie, ierror, iie, ik, ind, io, iio, &
125      iuo, j, jo, juo, neigneeded
126  integer :: nuotot2, nuo2, BlockSize2, maxo2, neigwanted2
127  complex(dp) :: kphs, qphs
128  complex(dp) :: cicj, D11, D22, D12, D21
129  real(dp) ::  kxij, qxij, ee, qe, t, qwh(3)
130
131#ifdef DEBUG
132  call write_debug( '    PRE diag2kspiral' )
133#endif
134
135  ! Easier
136  nuotot2 = nuotot * 2
137  neigwanted2 = neigwanted * 2
138  nuo2 = nuo * 2
139  maxo2 = maxo * 2
140  BlockSize2 = BlockSize * 2
141  qwh(:) = qw(:) * 0.5_dp
142
143! Find eigenvalues at every k point ...............................
144  do ik = 1,nk
145
146    call setup_k()
147
148    ! Find eigenvalues
149    ! Possible memory optimization: equivalence Haux and psi
150    call cdiag(Haux,Saux,nuotot2,nuotot2,nuo2,eo(1,ik),psi, &
151        -neigwanted2,iscf,ierror,BlockSize2)
152    if ( ierror /= 0 ) then
153      call die('Terminating due to failed diagonalisation')
154    end if
155  end do
156
157  ! Check if we are done ................................................
158  if (.not.getD) then
159#ifdef DEBUG
160    call write_debug( '    POS diag2kspiral' )
161#endif
162    return
163  end if
164
165  ! Find new Fermi energy and occupation weights ........................
166  call fermid(2, 1, nk, wk, maxo2, neigwanted2, eo, &
167      temp, qtot, qo, ef, Entropy )
168
169  ! Find weights for local density of states ............................
170  if ( e1 < e2 ) then
171    t = max( temp, 1.d-6 )
172    do ik = 1,nk
173      do io = 1, neigwanted2
174        qo(io,ik) = wk(ik) * &
175            ( stepf( (eo(io,ik)-e2)/t ) - stepf( (eo(io,ik)-e1)/t ) )
176      end do
177    end do
178  end if
179
180! New density and energy-density matrices of unit-cell orbitals .......
181  Dnew(:,:) = 0._dp
182  Enew(:,:) = 0._dp
183
184  do ik = 1, nk
185
186    ! Find maximum eigenvector that is required for this k point
187    neigneeded = 0
188    ie = neigwanted2
189    do while ( ie > 0 .and. neigneeded == 0 )
190      qe = qo(ie,ik)
191      if ( abs(qe) > occtol ) neigneeded = ie
192      ie = ie - 1
193    end do
194
195    call setup_k()
196
197    call cdiag(Haux,Saux,nuotot2,nuotot2,nuo2,eo(1,ik),psi, &
198        neigneeded,iscf,ierror,BlockSize2)
199
200    ! Check error flag and take appropriate action
201    if ( ierror > 0 ) then
202      call die('Terminating due to failed diagonalisation')
203    else if ( ierror < 0 ) then
204
205      call setup_k()
206      call cdiag(Haux,Saux,nuotot2,nuotot2,nuo2,eo(1,ik),psi, &
207          neigneeded,iscf,ierror,BlockSize2)
208
209    end if
210
211! Store the products of eigenvectors in matrices Dk and Ek
212! WARNING: Dk and Ek may be EQUIVALENCE'd to Haux and Saux
213    Dk = cmplx(0.0_dp, 0.0_dp, dp)
214    Ek = cmplx(0.0_dp, 0.0_dp, dp)
215
216    BNode = 0
217    iie = 0
218    do ie = 1, neigneeded
219      if ( Node == BNode ) then
220        iie = iie + 1
221        do j = 1, nuotot
222          aux(1,j) = psi(1,j,iie) ! c_i,up
223          aux(2,j) = psi(2,j,iie) ! c_i,dn
224        end do
225      end if
226#ifdef MPI
227      call MPI_Bcast(aux(1,1),nuotot2,MPI_double_complex,BNode, &
228          MPI_Comm_World,MPIerror)
229#endif
230      qe = qo(ie,ik)
231      ee = qe * eo(ie,ik)
232      do iuo = 1,nuo
233        call LocalToGlobalOrb(iuo,Node,Nodes,iio)
234        do juo = 1, nuotot
235!------- 1,1 -----------------------------------------------------------
236          cicj = conjg(aux(1,juo)) * aux(1,iio)
237          Dk(1,juo,1,iuo) = Dk(1,juo,1,iuo) + qe * cicj
238          Ek(1,juo,1,iuo) = Ek(1,juo,1,iuo) + ee * cicj
239!------- 2,2 -----------------------------------------------------------
240          cicj = conjg(aux(2,juo)) * aux(2,iio)
241          Dk(2,juo,2,iuo) = Dk(2,juo,2,iuo) + qe * cicj
242          Ek(2,juo,2,iuo) = Ek(2,juo,2,iuo) + ee * cicj
243!------- 1,2 -----------------------------------------------------------
244          cicj = conjg(aux(2,juo)) * aux(1,iio)
245          Dk(1,juo,2,iuo) = Dk(1,juo,2,iuo) + qe * cicj
246          Ek(1,juo,2,iuo) = Ek(1,juo,2,iuo) + ee * cicj
247!------- 2,1 -----------------------------------------------------------
248          cicj = conjg(aux(1,juo)) * aux(2,iio)
249          Dk(2,juo,1,iuo) = Dk(2,juo,1,iuo) + qe * cicj
250          Ek(2,juo,1,iuo) = Ek(2,juo,1,iuo) + ee * cicj
251        end do
252      end do
253      BTest = ie/BlockSize2
254      if ( BTest*BlockSize2 == ie ) then
255        BNode = BNode + 1
256        if ( BNode >= Nodes ) BNode = 0
257      end if
258    end do
259
260! Add contribution to density matrices of unit-cell orbitals
261    do iuo = 1,nuo
262      do j = 1,numd(iuo)
263        ind = listdptr(iuo) + j
264        jo = listd(ind)
265        juo = indxuo(jo)
266        kxij = kpoint(1,ik) * xij(1,ind) + &
267            kpoint(2,ik) * xij(2,ind) + &
268            kpoint(3,ik) * xij(3,ind)
269        kphs = exp(cmplx(0.0_dp, - kxij, dp))
270        qxij = qwh(1) * xij(1,ind) + &
271            qwh(2) * xij(2,ind) + &
272            qwh(3) * xij(3,ind)
273        qphs = exp(cmplx(0.0_dp, - qxij, dp))
274
275        D11 = Dk(1,juo,1,iuo) * kphs * qphs
276        D21 = Dk(2,juo,1,iuo) * kphs * conjg(qphs)
277        D12 = Dk(1,juo,2,iuo) * kphs * qphs
278        D22 = Dk(2,juo,2,iuo) * kphs * conjg(qphs)
279
280        ! Make DM Spin-box hermitian
281        D12 = 0.5_dp * (D12 + conjg(D21))
282
283        Dnew(ind,1) = Dnew(ind,1) + real(D11, dp)
284        Dnew(ind,2) = Dnew(ind,2) + real(D22, dp)
285        Dnew(ind,3) = Dnew(ind,3) + real(D12, dp)
286        Dnew(ind,4) = Dnew(ind,4) - aimag(D12)
287
288        D11 = Ek(1,juo,1,iuo) * kphs * qphs
289        D21 = Ek(2,juo,1,iuo) * kphs * conjg(qphs)
290        D12 = Ek(1,juo,2,iuo) * kphs * qphs
291        D22 = Ek(2,juo,2,iuo) * kphs * conjg(qphs)
292
293        ! Make EDM Spin-box hermitian
294        D12 = 0.5_dp * (D12 + conjg(D21))
295
296        Enew(ind,1) = Enew(ind,1) + real(D11, dp)
297        Enew(ind,2) = Enew(ind,2) + real(D22, dp)
298        Enew(ind,3) = Enew(ind,3) + real(D12, dp)
299        Enew(ind,4) = Enew(ind,4) - aimag(D12)
300
301      end do
302    end do
303
304  end do
305
306#ifdef DEBUG
307  call write_debug( '    POS diag2ksipral' )
308#endif
309
310contains
311
312  subroutine setup_k()
313
314    ! Initialize Hamiltonian and overlap matrices in full format
315    ! Index i is for real/imag parts
316    ! Indices is and js are for spin components
317    ! Indices iuo and juo are for orbital components:
318    ! Haux(i,js,juo,is,iuo) = <js,juo|H|is,iuo>
319    Saux = cmplx(0.0_dp, 0.0_dp, dp)
320    Haux = cmplx(0.0_dp, 0.0_dp, dp)
321
322    ! Transfer S,H matrices from sparse format in supercell to
323    ! full format in unit cell
324    ! Convention: ispin=1 => H11, ispin=2 => H22,
325    !             ispin=3 => Real(H12), ispin=4 => Imag(H12)
326    do iuo = 1,nuo
327      do j = 1,numh(iuo)
328        ind = listhptr(iuo) + j
329        jo = listh(ind)
330        juo = indxuo(jo)
331        kxij = kpoint(1,ik) * xij(1,ind) + &
332            kpoint(2,ik) * xij(2,ind) + &
333            kpoint(3,ik) * xij(3,ind)
334        kphs = exp(cmplx(0.0_dp, - kxij, dp))
335        qxij = qwh(1) * xij(1,ind) + &
336            qwh(2) * xij(2,ind) + &
337            qwh(3) * xij(3,ind)
338        qphs = exp(cmplx(0.0_dp, - qxij, dp))
339        D11 = kphs * qphs
340        D22 = kphs * conjg(qphs)
341
342        Saux(1,juo,1,iuo) = Saux(1,juo,1,iuo) + S(ind) * D11
343        Saux(2,juo,2,iuo) = Saux(2,juo,2,iuo) + S(ind) * D22
344
345        Haux(1,juo,1,iuo) = Haux(1,juo,1,iuo) + H(ind,1) * D11
346        Haux(2,juo,1,iuo) = Haux(2,juo,1,iuo) + cmplx(H(ind,3), H(ind,4),dp) * D22
347        Haux(1,juo,2,iuo) = Haux(1,juo,2,iuo) + cmplx(H(ind,3), -H(ind,4),dp) * D11
348        Haux(2,juo,2,iuo) = Haux(2,juo,2,iuo) + H(ind,2) * D22
349
350      end do
351    end do
352
353  end subroutine setup_k
354
355end subroutine diag2kspiral
356