PROGRAM mcmtorsion
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      character*80 fn
      parameter (NAT0=10000)
      DIMENSION XX(3),YY(3),A(3),P(3),E0(3),
     1X(NAT0),Y(NAT0),Z(NAT0)
      logical lex
      PI=3.141592653589D0
C     
      WRITE(*,*)'Torsion, angle, distance monitoring'
      WRITE(*,*)'Parameter and list of files defined in LIST.FIL'
      inquire(file='LIST.FIL',exist=lex)
      if(.not.lex)then
       write(6,*)' LIST.FIL not found.'
      else
       open(3,file='LIST.FIL')
       read(3,*)np,L0,L1,L2,L3
       write(6,600)np,L0,L1,L2,L3
600    format(' LIST.FIL found, ',i6,' points, angle ',4i4,'.')
c
       do 1 i=1,np
       read(3,300)fn
300    format(a80)

       open(31,file=fn)

3188   read(31,*,end=3199,err=3199)
       read(31,*,end=3199,err=3199)nat
       if(nat.gt.NAT0)then
        write(6,*)' Too many atoms.'
	close(31)
	stop
       endif
       do 2 ia=1,nat
2      read(31,*)idum,X(ia),Y(ia),Z(ia)

       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)
     1 CA=-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(ANGLE.lt.0.0d0)ANGLE=ANGLE+360.0D0
c      torsion, distance, band angle
       WRITE(6,322)i,AD,BANGLE,ANGLE
       goto 3188
3199   close(31)
1      continue
322    format(i5,3f15.6)
c
       close(31)
      endif
      stop
      end