PROGRAM tinkertorsion
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      parameter (NAT0=10000,np0=100,ndp0=200)
      character*80 fn
      character*10 parname(np0),s10
      DIMENSION XX(3),YY(3),A(3),P(3),E0(3),X(NAT0),Y(NAT0),Z(NAT0),
     2itypep(np0),L0p(np0),L1p(np0),L2p(np0),L3p(np0),
     1xmin(np0),xmax(np0),ndp(np0),distr(np0,ndp0),an(np0)
      logical lex,ldistr
      PI=3.141592653589D0
C     
      WRITE(6,*)'Tinker torsion angle trajectory.'
      inquire(file='LIST.PAR',exist=lex)
      if(.not.lex)call report(' LIST.PAR not found.')
      inquire(file='STRUCT.LST',exist=lex)
      if(.not.lex)call report(' STRUCT.LST not found.')

      open(45,file='LIST.PAR')
c
c     number of coordinates (if negative,calculate distributions):
      read(45,*)nag
      ldistr=nag.lt.0
      nag=abs(nag)
      if(nag.gt.np0)call report(' Too many parameters')
      do 211 iag=1,nag
      read(45,4501)parname(iag)
4501  format(a10)
      if(ldistr)then
       read(45,*)itypep(iag),L0p(iag),L1p(iag),L2p(iag),L3p(iag),
     1 ndp(iag),xmin(iag),xmax(iag)
       if(ndp(iag).gt.ndp0)call report('too many distribution points')
       an(iag)=0.0d0
       do 10 i=1,ndp(iag)
10     distr(iag,i)=0.0d0
      else
       read(45,*)itypep(iag),L0p(iag),L1p(iag),L2p(iag),L3p(iag)
      endif
211   continue
      close(45)

      open(32,file='CONSIGNE.TXT')
      do 1 i=1,nag
      write(6,*)i,parname(i)
1     write(32,324)parname(i)
324   format(a10,$)
      write(32,*)
      close(32)

      open(32,file='LIST.PRN')
      open(33,file='STRUCT.LST')
      is=0
991   read(33,450,end=99,err=99)fn
450   format(a80)
      is=is+1

      open(31,file=fn)
c
c     Is it an MCM file?
      do 300 i=len(fn),1,-1
300   if(fn(i:i).ne.' ')goto 301
301   if(fn(i-1:i).eq.'.x'.or.fn(i-1:i).eq.'.X')then
       read(31,*)
       read(31,*)nat
       if(nat.gt.NAT0)call report('Too many atoms')
       do 21 ia=1,nat
21     read(31,*)idumm,X(ia),Y(ia),Z(ia)
      else
       read(31,*)nat
       if(nat.gt.NAT0)call report('Too many atoms')
       do 2 ia=1,nat
2      read(31,310)X(ia),Y(ia),Z(ia)
310    format(11x,3f12.6)
      endif
      close(31)
      
      write(32,320)is
320   format(i6,$)

      do 11 i=1,nag
      itype=abs(itypep(i))
      L0=L0p(i)
      L1=L1p(i)
      L2=L2p(i)
      L3=L3p(i)
      A(1)=X(L2)-X(L3)
      A(2)=Y(L2)-Y(L3)
      A(3)=Z(L2)-Z(L3)
      P(1)=X(L0)-X(L1)
      P(2)=Y(L0)-Y(L1)
      P(3)=Z(L0)-Z(L1)
      E0(1)=X(L1)-X(L2)
      E0(2)=Y(L1)-Y(L2)
      E0(3)=Z(L1)-Z(L2)
      SSE0E0=E0(1)*E0(1)+E0(2)*E0(2)+E0(3)*E0(3)
      DO 8 KK=1,3
8     E0(KK)=E0(KK)/DSQRT(SSE0E0)
      SSPP=P(1)*P(1)+P(2)*P(2)+P(3)*P(3)
      AD=DSQRT(SSPP)
      SSAE0=A(1)*E0(1)+A(2)*E0(2)+A(3)*E0(3)
      SSPE0=P(1)*E0(1)+P(2)*E0(2)+P(3)*E0(3)
      DO 9 KK=1,3
      XX(KK)=A(KK)-E0(KK)*SSAE0
9     YY(KK)=P(KK)-E0(KK)*SSPE0
      SSXXXX=XX(1)*XX(1)+XX(2)*XX(2)+XX(3)*XX(3)
      SSYYYY=YY(1)*YY(1)+YY(2)*YY(2)+YY(3)*YY(3)
      SSXXYY=XX(1)*YY(1)+XX(2)*YY(2)+XX(3)*YY(3)
      APX=XX(2)*YY(3)-XX(3)*YY(2)
      APY=XX(3)*YY(1)-XX(1)*YY(3)
      APZ=XX(1)*YY(2)-XX(2)*YY(1)
      SSE0AP=APX*E0(1)+APY*E0(2)+APZ*E0(3)
      SIGN=1.0D0
      IF (SSE0AP.NE.0.0D0)SIGN=-SSE0AP/ABS(SSE0AP)
      CB=0.0D0
      IF (AD.NE.0.0D0)CB=-SSPE0/AD
      CA=0.0D0
      IF (SSXXXX.NE.0.0D0.AND.SSYYYY.NE.0.0D0)
     1CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY)
      BANGLE=0.0D0
      ANGLE=0.0D0
      IF (CB.EQ.-1.0D0)BANGLE=PI
      IF (ABS(CB).LT.1.0D0)BANGLE=ACOS(CB)
      IF (CA.LE.-1.0D0) ANGLE=PI
      IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA)
      ANGLE=180.0D0* ANGLE/PI
      BANGLE=180.0D0*BANGLE/PI
      ANGLE=SIGN*ANGLE
      if(itypep(i).lt.0)then
       if(ANGLE.lt.-180.0d0)ANGLE= 360.0d0+ANGLE
       if(ANGLE.gt.+180.0d0)ANGLE=-360.0d0+ANGLE
      else
       if(ANGLE.lt.0.0d0)ANGLE=360.0d0+ANGLE
      endif
      if(itype.eq.2)par=AD
      if(itype.eq.3)par=BANGLE
      if(itype.eq.4)par=ANGLE
      if(ldistr)then
       an(i)=an(i)+1.0d0
       if(par.le.xmax(i).and.par.ge.xmin(i))then
        dd=(xmax(i)-xmin(i))/dble(ndp(i))
        j=int((par-xmin(i)+dd)/dd)
        if(j.lt.1)j=1
        if(j.gt.ndp(i))j=ndp(i)
        distr(i,j)=distr(i,j)+1.0d0
       endif
      endif
11    write(32,321)par
321   format(f10.3,$)
      write(32,*)

      goto 991
99    close(32)
      close(33)
c
      if(ldistr)then
       do 12 ip=1, nag
       write(s10,4502)ip
4502   format(i10)
       do 13 is=1,10
13     if(s10(is:is).ne.' ')goto 14
14     do 15 ie=10,1,-1
15     if(s10(ie:ie).ne.' ')goto 16
16     write(6,*)is,ie,s10
       open(33,file='DISTR.'//s10(is:ie)//'.PRN')
       dd=(xmax(ip)-xmin(ip))/dble(ndp(ip))
       fac=1.0d0/an(ip)/dd
       do 17 j=1,ndp(ip)
17     write(33,3333)xmin(ip)+dd*dble(2*j-1)/2.0d0,distr(ip,j)*fac
3333   format(2f12.6)
12     close(33)
       write(6,*)'Distribution files written'
      endif

      call report('LIST.PRN written')

      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end