      real function dihed (angnum, atoma, atomb, atomc, atomd)          DIHE001
c     routine to calculate the dihedral angle between the 4 atoms       DIHE002
c     given                                                             DIHE003
c                                                                       DIHE004
      include 'brkcln.par'                                              DIHE005
c                                                                       DIHE006
      real    atoma(ncoord), atomb(ncoord), atomc(ncoord), atomd(ncoord)DIHE007
c                                                                       DIHE008
      integer angnum                                                    DIHE009
c                                                                       DIHE010
      real    dihatm(ncoord,dihdat), codist(ncoord,dihdat-1),           DIHE011
     &        atmdst(dihdat-1), dotprd(dihdat-1,dihdat-1), ang,         DIHE012
     &        cosang, detant, atdist                                    DIHE013
c                                                                       DIHE014
      integer i, j, k                                                   DIHE015
c                                                                       DIHE016
      logical xeqy                                                      DIHE017
c                                                                       DIHE018
      do 100 i = 1, ncoord                                              DIHE019
      dihatm(i,1) = atoma(i)                                            DIHE020
      dihatm(i,2) = atomb(i)                                            DIHE021
      dihatm(i,3) = atomc(i)                                            DIHE022
      dihatm(i,4) = atomd(i)                                            DIHE023
CHECK v.3.2-->                                                          DIHE024
      do 99 j=5,dihdat                                                  DIHE025
        dihatm (i,j) = 0.0                                              DIHE026
 99   continue                                                          DIHE027
CHECK v.3.2<--                                                          DIHE028
  100 continue                                                          DIHE029
      do 200 i = 1, dihdat-1                                            DIHE030
      atmdst(i) = atdist(dihatm(1,i),dihatm(1,i+1))                     DIHE031
      do 150 j = 1, ncoord                                              DIHE032
      codist(j,i) = dihatm(j,i+1) - dihatm(j,i)                         DIHE033
  150 continue                                                          DIHE034
  200 continue                                                          DIHE035
      do 300 i = 1, dihdat-1                                            DIHE036
      do 300 j = 1, dihdat-1                                            DIHE037
      dotprd(i,j) = 0.0                                                 DIHE038
  300 continue                                                          DIHE039
      do 400 i = 1, dihdat-2                                            DIHE040
      do 400 j = i+1, dihdat-1                                          DIHE041
      do 400 k = 1, ncoord                                              DIHE042
      dotprd(i,j) = dotprd(i,j) + (codist(k,i) * codist(k,j))           DIHE043
  400 continue                                                          DIHE044
      do 500 i = 1, dihdat-2                                            DIHE045
      do 500 j = i+1, dihdat-1                                          DIHE046
      if (xeqy(atmdst(i),0.0) .or. xeqy(atmdst(j),0.0)) then            DIHE047
         dotprd(i,j) = 1.0                                              DIHE048
      else                                                              DIHE049
         dotprd(i,j) = ( dotprd(i,j) / atmdst(i) ) / atmdst(j)          DIHE050
      end if                                                            DIHE051
  500 continue                                                          DIHE052
      if (angnum .gt. ndihed .and. angnum .le. nmnang) then             DIHE053
         cosang = dotprd(1,3)                                           DIHE054
      else                                                              DIHE055
         if (xeqy(dotprd(1,2),1.0) .or. xeqy(dotprd(2,3),1.0)) then     DIHE056
            cosang = 1.0                                                DIHE057
         else                                                           DIHE058
            cosang = ((dotprd(1,2) * dotprd(2,3)) - dotprd(1,3)) /      DIHE059
     &               (sqrt(1.0 - (dotprd(1,2) ** 2)) *                  DIHE060
     &               sqrt(1.0 - (dotprd(2,3) ** 2)))                    DIHE061
         end if                                                         DIHE062
      end if                                                            DIHE063
      if (abs(cosang) .gt. 1.0) cosang = anint(cosang / abs(cosang))    DIHE064
      ang = acos(cosang) * radian                                       DIHE065
      if (angnum .le. ndihed .or. angnum .gt. nmnang) then              DIHE066
         detant = ((((((codist(1,1) * codist(2,2)) * codist(3,3)) -     DIHE067
     &                ((codist(1,1) * codist(2,3)) * codist(3,2))) +    DIHE068
     &                ((codist(1,2) * codist(2,3)) * codist(3,1))) -    DIHE069
     &                ((codist(1,2) * codist(2,1)) * codist(3,3))) +    DIHE070
     &                ((codist(1,3) * codist(2,1)) * codist(3,2))) -    DIHE071
     &                ((codist(1,3) * codist(2,2)) * codist(3,1))       DIHE072
         if (detant .lt. 0.0) ang = - ang                               DIHE073
      end if                                                            DIHE074
      dihed = ang                                                       DIHE075
      end                                                               DIHE076
