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