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 diagpol( ispin, nspin, nuo, no, nuotot,
9     .                    maxnh, numh, listhptr, listh, H, S,
10     .                    xij, indxuo, kpoint, eo, psi, ng,
11     .                    Haux, Saux )
12C *********************************************************************
13C Subroutine to calculate the eigenvalues and eigenvectors,
14C for given Hamiltonian and Overlap matrices (including
15C spin polarization), and for a given k_point
16C Written by DSP, March 1999. From the routine diagk of J.Soler.
17C Modified for parallel execution by J.D.Gale, March 2000.
18C **************************** INPUT **********************************
19C integer ispin               : Spin component which will be calculated
20C integer nspin               : Number of spin components (1 or 2)
21C integer nuo                 : Number of basis orbitals in unit cell
22C integer no                  : Number of basis orbitals in supercell
23C integer nuotot              : First dimension of eo, qo, last of xij
24C integer maxnh               : Maximum number of orbitals interacting
25C integer numh(nuo)           : Number of nonzero elements of each row
26C                               of hamiltonian matrix
27C integer listhptr(nuo)       : Pointer to the start of each row
28C                               of hamiltonian matrix
29C integer listh(maxnh)        : Nonzero hamiltonian-matrix element
30C                               column indexes for each matrix row
31C real*8  H(maxnh,nspin)      : Hamiltonian in sparse form
32C real*8  S(maxnh)            : Overlap in sparse form
33C real*8  xij(3,*)            : Vectors between orbital centers (sparse)
34C                               (not used if only gamma point)
35C integer indxuo(no)          : Index of equivalent orbital in unit cell
36C                               Unit cell orbitals must be the first in
37C                               orbital lists, i.e. indxuo.le.nuo, with
38C                               nuo the number of orbitals in unit cell
39C real*8  kpoint(3)           : k point vectors
40C real*8  psi                 : Workspace array
41C integer ng                  : first dimension of Haux, and Saux
42C real*8  Haux(ng,nuotot,nuo) : Workspace for dense H
43C real*8  Saux(ng,nuotot,nuo) : Workspace for dense S
44C *************************** OUTPUT **********************************
45C real*8 eo(nuotot)           : Eigenvalues
46C *************************** UNITS ***********************************
47C xij and kpoint must be in reciprocal coordinates of each other.
48C eo H.
49C *********************************************************************
50
51      use precision
52      use parallel,      only : BlockSize
53      use sys
54
55      implicit          none
56
57      integer
58     .  maxnh, nuotot, no, nspin, nuo, indxuo(no), listh(maxnh),
59     .  listhptr(nuo), numh(nuo), ng
60
61      real(dp)
62     .  eo(nuotot), H(maxnh,nspin), kpoint(3), S(maxnh),
63     .  xij(3,*), psi(ng,nuotot,nuo), Haux(ng,nuotot,nuo),
64     .  Saux(ng,nuotot,nuo)
65      external          rdiag, cdiag
66
67C  Internal variables .............................................
68      integer
69     .  ierror, ind, ispin, iuo, j, jo, juo
70      real(dp)
71     .  ckxij, kxij, skxij
72
73C Solve eigenvalue problem .........................................
74      Saux = 0.0d0
75      Haux = 0.0d0
76      do iuo = 1,nuo
77        do j = 1,numh(iuo)
78          ind = listhptr(iuo) + j
79          jo = listh(ind)
80          juo = indxuo(jo)
81          if(ng.eq.2) then
82           kxij = kpoint(1) * xij(1,ind) +
83     .            kpoint(2) * xij(2,ind) +
84     .            kpoint(3) * xij(3,ind)
85           ckxij = cos(kxij)
86           skxij = sin(kxij)
87           Saux(1,juo,iuo) = Saux(1,juo,iuo) + S(ind)*ckxij
88           Saux(2,juo,iuo) = Saux(2,juo,iuo) - S(ind)*skxij
89           Haux(1,juo,iuo) = Haux(1,juo,iuo) + H(ind,ispin)*ckxij
90           Haux(2,juo,iuo) = Haux(2,juo,iuo) - H(ind,ispin)*skxij
91          else
92           Saux(1,juo,iuo) = Saux(1,juo,iuo) + S(ind)
93           Haux(1,juo,iuo) = Haux(1,juo,iuo) + H(ind,ispin)
94          endif
95        enddo
96      enddo
97      if(ng.eq.2) then
98       call cdiag( Haux, Saux, nuotot, nuo, nuotot, eo, psi,
99     .            nuotot, 1, ierror, BlockSize)
100      else
101       call rdiag( Haux, Saux, nuotot, nuo, nuotot, eo, psi,
102     .            nuotot, 1, ierror, BlockSize)
103      endif
104C Check error flag and take appropriate action
105      if (ierror.gt.0) then
106        call die('Terminating due to failed diagonalisation')
107      elseif (ierror.lt.0) then
108C Repeat diagonalisation with increased memory to handle clustering
109        Saux = 0.0d0
110        Haux = 0.0d0
111        do iuo = 1,nuo
112          do j = 1,numh(iuo)
113            ind = listhptr(iuo) + j
114            jo = listh(ind)
115            juo = indxuo(jo)
116          if(ng.eq.2) then
117           kxij = kpoint(1) * xij(1,ind) +
118     .            kpoint(2) * xij(2,ind) +
119     .            kpoint(3) * xij(3,ind)
120           ckxij = cos(kxij)
121           skxij = sin(kxij)
122           Saux(1,juo,iuo) = Saux(1,juo,iuo) + S(ind)*ckxij
123           Saux(2,juo,iuo) = Saux(2,juo,iuo) - S(ind)*skxij
124           Haux(1,juo,iuo) = Haux(1,juo,iuo) + H(ind,ispin)*ckxij
125           Haux(2,juo,iuo) = Haux(2,juo,iuo) - H(ind,ispin)*skxij
126          else
127           Saux(1,juo,iuo) = Saux(1,juo,iuo) + S(ind)
128           Haux(1,juo,iuo) = Haux(1,juo,iuo) + H(ind,ispin)
129          endif
130        enddo
131      enddo
132      if(ng.eq.2) then
133       call cdiag( Haux, Saux, nuotot, nuo, nuotot, eo, psi,
134     .            nuotot, 1, ierror, BlockSize)
135      else
136       call rdiag( Haux, Saux, nuotot, nuo, nuotot, eo, psi,
137     .            nuotot, 1, ierror, BlockSize)
138      endif
139
140      endif
141
142      end
143