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