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_intramol_pressure
9CONTAINS
10subroutine remove_intramol_pressure(ucell,stress,na_u,xa,fa,mstress)
11 use precision, only: dp
12 use zmatrix, only: nZmol, nZmolStartAtom
13 use zmatrix, only: lUseZmatrix, nZmolAtoms
14
15 implicit none
16
17 real(dp), external :: volcel
18
19 integer, intent(in)   :: na_u
20 real(dp), intent(in)  :: ucell(3,3), xa(3,na_u), fa(3,na_u)
21 real(dp), intent(in)  :: stress(3,3)
22 real(dp), intent(out) :: mstress(3,3)
23
24 real(dp) :: virial, fmean, volume
25 integer  :: ifirst, ilast, natoms_mol, imol, ix, ia
26
27! Find intramolecular contributions to the pressure
28        volume = volcel(ucell)
29        virial = 0.0_dp
30        if (lUseZmatrix) then
31           ! Molecule by molecule
32           do imol = 1, nZmol
33              ifirst = nZmolStartAtom(imol)
34              natoms_mol = nZmolAtoms(imol)
35              ilast  = ifirst + natoms_mol - 1
36              do ix = 1,3
37               fmean = sum(fa(ix,ifirst:ilast)) / natoms_mol
38               do ia = ifirst, ilast
39                 virial = virial + xa(ix,ia) * (fa(ix,ia) - fmean)
40               enddo
41              enddo
42           enddo
43        else
44           ! No Zmatrix
45           ! Consider the whole system as a molecule, even if it is meaningless
46           do ix = 1,3
47              fmean = sum(fa(ix,1:na_u)) / na_u
48              do ia = 1,na_u
49                 virial = virial + xa(ix,ia) * (fa(ix,ia) - fmean)
50              enddo
51           enddo
52        endif
53!
54        mstress(1:3,1:3) = stress(1:3,1:3)
55        do ix = 1, 3
56           mstress(ix,ix) = mstress(ix,ix) + virial / volume / 3.0_dp
57        enddo
58
59end subroutine remove_intramol_pressure
60end module m_intramol_pressure
61