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