C
C $Id: dihedr.F,v 1.3 1998/07/16 16:39:36 jjv5 Exp arjan $
C
!     Modified by Antoine MARION :: 2013-11-08
!          --  add PBC
C------------------------------------------------------------------------
      SUBROUTINE DIHEDRpbc(pbc,XYZ,I,IA,IB,IC,DIH)
C
C     DETERMINES THE I-IA-IB-IC DIHEDRAL ANGLE IN RADIANS.  THE
C     ANGLE DIH IS POSITIVE IF IC IS LOCATED CLOCKWISE FROM I WHEN
C     VIEWING FROM IA THROUGH IB.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XYZ(3,*)
      double precision ix,iy,iz
      double precision iax,iay,iaz
      double precision ibx,iby,ibz
      double precision icx,icy,icz
      logical pbc

      ix = xyz(1,i)
      iy = xyz(2,i)
      iz = xyz(3,i)


      if (pbc) then
       call pbcxyz(i,ia,iax,iay,iaz)
       call pbcxyz(i,ib,ibx,iby,ibz)
       call pbcxyz(i,ic,icx,icy,icz)
      else
       iax = xyz(1,ia)
       iay = xyz(2,ia)
       iaz = xyz(3,ia)

       ibx = xyz(1,ib)
       iby = xyz(2,ib)
       ibz = xyz(3,ib)

       icx = xyz(1,ic)
       icy = xyz(2,ic)
       icz = xyz(3,ic)
      endif

C
C     SHIFT IA-I AND IB-IC BOND VECTORS TO A COMMON ORIGIN.
C

      AIIX = ix - iax
      AIIY = iy - iay
      AIIZ = iz - iaz
      BCX  = icx - ibx
      BCY  = icy - iby
      BCZ  = icz - ibz

!      AIIX = XYZ(1,I) - XYZ(1,IA)
!      AIIY = XYZ(2,I) - XYZ(2,IA)
!      AIIZ = XYZ(3,I) - XYZ(3,IA)
!      BCX = XYZ(1,IC) - XYZ(1,IB)
!      BCY = XYZ(2,IC) - XYZ(2,IB)
!      BCZ = XYZ(3,IC) - XYZ(3,IB)
C
C     FORM THE IA-IB BOND AXIS VECTOR.
C

      ABX = ibx - iax
      ABY = iby - iay
      ABZ = ibz - iaz

!      ABX = XYZ(1,IB) - XYZ(1,IA)
!      ABY = XYZ(2,IB) - XYZ(2,IA)
!      ABZ = XYZ(3,IB) - XYZ(3,IA)
C
C     REMOVE FROM (AIIX,AIIY,AIIZ) AND (BCX,BCY,BCZ) ANY PROJECTION ALONG
C     THE (ABX,ABY,ABZ) AXIS.
C
      DOT1 = AIIX*ABX + AIIY*ABY + AIIZ*ABZ
      ABSQR = ABX**2 + ABY**2 + ABZ**2
      PROJ1 = DOT1/ABSQR
      AIIX = AIIX - PROJ1*ABX
      AIIY = AIIY - PROJ1*ABY
      AIIZ = AIIZ - PROJ1*ABZ
      DOT2 = BCX*ABX + BCY*ABY + BCZ*ABZ
      PROJ2 = DOT2/ABSQR
      BCX = BCX - PROJ2*ABX
      BCY = BCY - PROJ2*ABY
      BCZ = BCZ - PROJ2*ABZ
C
C     COMPUTE THE CROSS-PRODUCT (AIIX,AIIY,AIIZ) X (BCX,BCY,BCZ).  STORE
C     IT IN THE VECTOR (AIBCX,AIBCY,AIBCZ).
C
      AIBCX = AIIY*BCZ - AIIZ*BCY
      AIBCY = AIIZ*BCX - AIIX*BCZ
      AIBCZ = AIIX*BCY - AIIY*BCX
C
C     IF (AIBCX,AIBCY,AIBCZ) POINTS IN THE SAME DIRECTION AS
C     (ABX,ABY,ABZ) THEN IC IS LOCATED CLOCKWISE FROM I WHEN
C     VIEWED FROM IA TOWARD IB.  THUS, IN MOVING ALONG THE PATH
C     I-IA-IB-IC, A CLOCKWISE OR POSITIVE ANGLE IS OBSERVED.
C     TO DETERMINE WHETHER THESE VECTORS POINT IN THE SAME OR
C     OPPOSITE DIRECTIONS, COMPUTE THEIR DOT PRODUCT.
C
      DOT3 = AIBCX*ABX + AIBCY*ABY + AIBCZ*ABZ
      DIREC = SIGN(1.0D0,DOT3)
C
C     COMPUTE THE DIHEDRAL ANGLE DIH.
C
      DOT4 = AIIX*BCX + AIIY*BCY + AIIZ*BCZ
      AILENG = SQRT(AIIX**2 + AIIY**2 + AIIZ**2)
      BCLENG = SQRT(BCX**2 + BCY**2 + BCZ**2)
      IF(ABS(AILENG*BCLENG).LT.1.0D-5)THEN
        COSINE = 1.0D0
      ELSE
        COSINE = DOT4/(AILENG*BCLENG)
      ENDIF
      COSINE = MAX(-1.0D0,COSINE)
      COSINE = MIN(1.0D0,COSINE)
      DIH = ACOS(COSINE)*DIREC
      RETURN
      END