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!
8module m_overfsm
9
10  use precision,     only : dp
11  use parallel,      only : Node, Nodes
12  use parallelsubs,  only : GlobalToLocalOrb
13  use atmfuncs,      only : rcut, orb_gindex
14  use neighbour,     only : jna=>jan, r2ij, xij, mneighb, reset_neighbour_arrays
15  use alloc,         only : re_alloc, de_alloc
16  use m_new_matel,   only : new_matel
17  use t_spin, only: tSpin
18
19  implicit none
20
21  private
22  public :: overfsm
23
24contains
25
26  subroutine overfsm(na_u, na_s, no_l, no_s, scell, xa, indxua, rmaxo, &
27      lasto, iphorb, isa, maxnd, numd, listdptr, listd, &
28      spin, Escf, fa, stress)
29
30    ! *********************************************************************
31    ! Overlap matrix and orthonormalization contribution to forces and stress.
32    ! Energies in Ry. Lengths in Bohr.
33    ! Writen by J.Soler and P.Ordejon, August-October'96, June'98
34    ! **************************** INPUT **********************************
35    ! integer na_u             : Number of atoms in unit cell
36    ! integer na_s             : Number of atoms in supercell
37    ! integer no_l             : Number of orbitals in unit-cell (local)
38    ! integer no_s             : Number of orbitals in supercell
39    ! real*8  scell(3,3)       : Supercell vectors SCELL(IXYZ,IVECT)
40    ! real*8  xa(3,na_s)       : Atomic positions in cartesian coordinates
41    ! integer indxua(na_s)     : Index of equivalent atom in unit cell
42    ! real*8  rmaxo            : Maximum cutoff for atomic orbitals
43    ! integer maxnd            : First dimension of Escf and listd
44    ! integer lasto(0:na_s)    : Last orbital index of each atom
45    ! integer iphorb(no_s)     : Orbital index of each orbital in its atom
46    ! integer isa(na_s)        : Species index of each atom
47    ! integer numd(no_l)       : Number of nonzero elements of each row
48    !                            of Escf
49    ! integer listdptr(no_l)   : Pointer to start of rows (-1) of Escf
50    ! integer listd(maxnd)     : Column indexes of the nonzero elements
51    !                            of each row of Escf
52    ! type(tSpin) spin         : Spin configuration
53    ! integer Escf(maxnd,spin%EDM): Energy-Density matrix
54    ! ********************** INPUT and OUTPUT *****************************
55    ! real*8  fa(3,na_u)       : Atomic forces (Orthog. part added to input)
56    ! real*8  stress(3,3)      : Stress tensor (Orthog. part added to input)
57    ! *********************************************************************
58
59
60    integer, intent(in) :: na_u, na_s, no_l, no_s
61    integer, intent(in) :: indxua(na_s), iphorb(no_s), isa(na_s), lasto(0:na_s)
62    integer, intent(in) :: maxnd
63    integer, intent(in) :: listd(maxnd), numd(no_l), listdptr(no_l)
64
65    real(dp), intent(inout) :: fa(3,na_u), stress(3,3)
66    type(tSpin), intent(in) :: spin
67    real(dp), intent(in) :: scell(3,3), Escf(maxnd,spin%EDM), rmaxo, xa(3,na_s)
68
69    ! Internal variables ......................................................
70
71    integer :: ia, ind, io, ioa, is, ispin, ix, iio, ig, jg
72    integer :: j, ja, jn, jo, joa, js, jua, jx, nnia
73
74    real(dp) :: fij(3), grSij(3), rij, Sij, volcel, volume
75
76    real(dp), dimension(:), pointer ::  Di => null()
77
78    external :: timer
79
80    ! Start timer
81    call timer( 'overfsm', 1 )
82
83    volume = na_u * volcel(scell) / na_s
84
85    ! Initialize neighb subroutine
86    call mneighb( scell, 2.d0*rmaxo, na_s, xa, 0, 0, nnia)
87
88    ! Allocate local memory
89    call re_alloc( Di, 1, no_s, 'Di', 'overfsm' )
90
91    Di(1:no_s) = 0.0d0
92
93    do ia = 1, na_u
94      is = isa(ia)
95      call mneighb( scell, 2.d0*rmaxo, na_s, xa, ia, 0, nnia )
96      do io = lasto(ia-1)+1,lasto(ia)
97
98        ! Is this orbital on this Node?
99        call GlobalToLocalOrb(io,Node,Nodes,iio)
100        if ( iio > 0 ) then
101
102          ! Valid orbital
103          ioa = iphorb(io)
104          ig = orb_gindex(is,ioa)
105          do j = 1,numd(iio)
106            ind = listdptr(iio)+j
107            jo = listd(ind)
108            do ispin = 1, spin%spinor
109              Di(jo) = Di(jo) + Escf(ind,ispin)
110            end do
111          end do
112
113          do jn = 1,nnia
114            ja = jna(jn)
115            jua = indxua(ja)
116            rij = sqrt( r2ij(jn) )
117            do jo = lasto(ja-1) + 1 , lasto(ja)
118              joa = iphorb(jo)
119              js = isa(ja)
120
121              if ( rcut(is,ioa) + rcut(js,joa) > rij ) then
122                jg = orb_gindex(js,joa)
123                call new_MATEL( 'S', ig, jg, xij(1:3,jn), Sij, grSij )
124                do ix = 1,3
125                  fij(ix) = - Di(jo) * grSij(ix)
126                  fa(ix,ia)  = fa(ix,ia)  + fij(ix)
127                  fa(ix,jua) = fa(ix,jua) - fij(ix)
128                  do jx = 1,3
129                    stress(jx,ix) = stress(jx,ix) + xij(jx,jn) * fij(ix) / volume
130                  end do
131                end do
132              end if
133
134            end do
135          end do
136
137          ! Reset arrays for next orbital interaction
138          do j = 1,numd(iio)
139            jo = listd(listdptr(iio)+j)
140            Di(jo) = 0._dp
141          end do
142
143        end if
144
145      end do
146    end do
147
148    ! Deallocate local memory
149    call reset_neighbour_arrays( )
150    call de_alloc( Di, 'Di', 'overfsm' )
151
152    ! Finish timer
153    call timer( 'overfsm', 2 )
154
155  end subroutine overfsm
156
157end module m_overfsm
158
159
160
161
162