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