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 module m_kinefsm 9 public :: kinefsm 10 CONTAINS 11 subroutine kinefsm(nua, na, no, scell, xa, indxua, rmaxo, 12 . maxnh, maxnd, lasto, iphorb, isa, 13 . numd, listdptr, listd, numh, listhptr, listh, 14 . nspin, Dscf, Ekin, fa, stress, H, 15 . matrix_elements_only ) 16C ********************************************************************* 17C Kinetic contribution to energy, forces, stress and matrix elements. 18C Energies in Ry. Lengths in Bohr. 19C Writen by J.Soler and P.Ordejon, August-October'96, June'98 20C **************************** INPUT ********************************** 21C integer nua : Number of atoms in unit cell 22C integer na : Number of atoms in supercell 23C integer no : Number of orbitals in supercell 24C real*8 scell(3,3) : Supercell vectors SCELL(IXYZ,IVECT) 25C real*8 xa(3,na) : Atomic positions in cartesian coordinates 26c integer indxua(na) : Index of equivalent atom in unit cell 27C real*8 rmaxo : Maximum cutoff for atomic orbitals 28C integer maxnh : First dimension of H and listh 29C integer maxnd : First dimension of Dscf and listd 30C integer lasto(0:na) : Last orbital index of each atom 31C integer iphorb(no) : Orbital index of each orbital in its atom 32C integer isa(na) : Species index of each atom 33C integer numd(nuo) : Number of nonzero elements of each row 34C of the density matrix 35C integer listdptr(nuo) : Pointer to the start of rows (-1) of 36C the density matrix 37C integer listd(maxnh) : Column indexes of the nonzero elements 38C of each row of the density matrix 39C integer numh(nuo) : Number of nonzero elements of each row 40C of the hamiltonian matrix 41C integer listhptr(nuo) : Pointer to the start of rows (-1) of 42C the hamiltonian matrix 43C integer listh(maxnh) : Column indexes of the nonzero elements 44C of each row of the hamiltonian matrix 45C integer nspin : Number of spin components of Dscf and H 46C If computing only matrix elements, it 47C can be set to 1. 48C logical matrix_elements_only: 49C integer Dscf(maxnd,nspin): Density matrix. Not touched if computing 50C only matrix elements. 51C **************************** OUTPUT ********************************* 52C real*8 Ekin : Kinetic energy in unit cell 53C ********************** INPUT and OUTPUT ***************************** 54C real*8 fa(3,nua) : Atomic forces (Kinetic part added to input) 55C real*8 stress(3,3) : Stress tensor (Kinetic part added to input) 56C real*8 H(maxnh,nspin) : Hamiltonian (Kinetic part added to input) 57C ********************************************************************* 58C 59C Modules 60C 61 use precision, only : dp 62 use parallel, only : Node, Nodes 63 use parallelsubs, only : GlobalToLocalOrb 64 use atmfuncs, only : rcut, orb_gindex 65 use neighbour, only : jna=>jan, r2ij, xij, mneighb, 66 & reset_neighbour_arrays 67 use m_new_matel, only : new_matel 68 use alloc, only : re_alloc, de_alloc 69 70 implicit none 71 72 integer :: maxnd, maxnh, na, no, nspin, nua 73 74 integer 75 . indxua(na), iphorb(no), isa(na), lasto(0:na), 76 . listd(maxnd), listh(maxnh), numd(*), numh(*), listdptr(*), 77 . listhptr(*) 78 79 real(dp) 80 . scell(3,3), Dscf(maxnd,nspin), Ekin, 81 . fa(3,nua), H(maxnh,nspin), rmaxo, 82 . stress(3,3), xa(3,na) 83 logical, intent(in) :: matrix_elements_only 84 85C Internal variables .................................................. 86 87 integer 88 . ia, ind, io, iio, ioa, is, ispin, ix, ig, jg, 89 . j, ja, jn, jo, joa, js, jua, jx, nnia 90 91 real(dp) 92 . fij(3), grTij(3) , rij, Tij, volcel, volume 93 94 real(dp), dimension(:), pointer :: Di, Ti 95 96 external :: volcel, timer 97C ...................... 98 99C Start timer 100 call timer( 'kinefsm', 1 ) 101 102C Allocate local memory 103 if (.not. matrix_elements_only) then 104 nullify( Di ) 105 call re_alloc( Di, 1, no, 'Di', 'kinefsm' ) 106 Di(1:no) = 0.0_dp 107 endif 108 nullify( Ti ) 109 call re_alloc( Ti, 1, no, 'Ti', 'kinefsm' ) 110 Ti(1:no) = 0.0_dp 111 112 volume = nua * volcel(scell) / na 113 114 call mneighb( scell, 2.0d0*rmaxo, na, xa, 0, 0, nnia ) 115 116 Ekin = 0.0_dp 117 118 do ia = 1,nua 119 is = isa(ia) 120 call mneighb( scell, 2.0d0*rmaxo, na, xa, ia, 0, nnia ) 121 do io = lasto(ia-1)+1,lasto(ia) 122 123C Is this orbital on this Node? 124 call GlobalToLocalOrb(io,Node,Nodes,iio) 125 if (iio.gt.0) then ! Local orbital 126 ioa = iphorb(io) 127 ig = orb_gindex(is,ioa) 128 129 if (.not. matrix_elements_only) then 130 do j = 1,numd(iio) 131 ind = listdptr(iio) + j 132 jo = listd(ind) 133 do ispin = 1,nspin 134 Di(jo) = Di(jo) + Dscf(ind,ispin) 135 enddo 136 enddo 137 endif 138 do jn = 1,nnia 139 ja = jna(jn) 140 jua = indxua(ja) 141 rij = sqrt( r2ij(jn) ) 142 do jo = lasto(ja-1)+1,lasto(ja) 143 joa = iphorb(jo) 144 js = isa(ja) 145 jg = orb_gindex(js,joa) 146 if (rcut(is,ioa)+rcut(js,joa) .gt. rij) then 147 call new_MATEL( 'T', ig, jg, xij(1:3,jn), 148 . Tij, grTij ) 149 Ti(jo) = Ti(jo) + Tij 150 151 if (.not. matrix_elements_only) then 152 Ekin = Ekin + Di(jo) * Tij 153 do ix = 1,3 154 fij(ix) = Di(jo) * grTij(ix) 155 fa(ix,ia) = fa(ix,ia) + fij(ix) 156 fa(ix,jua) = fa(ix,jua) - fij(ix) 157 do jx = 1,3 158 stress(jx,ix) = stress(jx,ix) + 159 . xij(jx,jn) * fij(ix) / volume 160 enddo 161 enddo 162 endif 163 endif 164 165 enddo 166 enddo 167 if (.not. matrix_elements_only) then 168 do j = 1,numd(iio) 169 jo = listd(listdptr(iio)+j) 170 Di(jo) = 0.0_dp 171 enddo 172 endif 173 do j = 1,numh(iio) 174 ind = listhptr(iio)+j 175 jo = listh(ind) 176 do ispin = 1,nspin 177 H(ind,ispin) = H(ind,ispin) + Ti(jo) 178 enddo 179 Ti(jo) = 0.0_dp 180 enddo 181 endif 182 enddo 183 enddo 184 185C Deallocate local memory 186! call new_MATEL( 'T', 0, 0, 0, 0, xij, Tij, grTij ) 187 call reset_neighbour_arrays( ) 188 call de_alloc( Ti, 'Ti', 'kinefsm' ) 189 if (.not. matrix_elements_only) then 190 call de_alloc( Di, 'Di', 'kinefsm' ) 191 endif 192 193C Finish timer 194 call timer( 'kinefsm', 2 ) 195 196 end subroutine kinefsm 197 end module m_kinefsm 198