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