PROGRAM xtorsion
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      parameter (ndp0=200,naa0=9)
      character*10 s10
      real*8 XX(3),YY(3),A(3),P(3),E0(3),B0(3),AP(3),
     1b(3),c(3),v1(3),v2(3),rab(3),va(3),vb(3),u(3),va0(3),vb0(3)
      character*10,allocatable:: parname(:)
      real*8,allocatable::X(:,:),xmin(:),xmax(:),distr(:,:),an(:),
     1xt(:,:)
      integer*4,allocatable::ndp(:),itypep(:),L(:,:)
      logical lex,ldistr,lper
      PI=3.141592653589D0
C     
      WRITE(6,*)'FILE.X analysis'
      inquire(file='LIST.PAR',exist=lex)
      if(.not.lex)call report(' LIST.PAR not found.')
      inquire(file='FILE.X',exist=lex)
      if(.not.lex)call report(' FILE.X not found.')
      call readopt(noff,nmol,lper,ldistr)

      open(45,file='LIST.PAR')
c
c     number of coordinates (if negative,calculate distributions):
      read(45,*)nag
      ldistr=ldistr.or.(nag.lt.0)
      nag=abs(nag)
      allocate(parname(nag),itypep(nag),L(nag,naa0),
     1xmin(nag),xmax(nag),ndp(nag),distr(nag,ndp0),an(nag))
      do 211 iag=1,nag
      read(45,4501)parname(iag)
4501  format(a10)
      read(45,*)naa
      naa=abs(naa)
      if(naa.eq. 5)naa=4
      if(naa.eq. 6)naa=4
      if(naa.eq. 7)naa=6
      if(naa.eq. 8)naa=5
      if(naa.eq. 9)naa=5
      if(naa.eq.10)naa=4
      if(naa.gt.naa0)call report('Too many atoms in a list')
      backspace 45
      if(ldistr)then
       read(45,*)itypep(iag),(L(iag,iaa),iaa=1,naa),
     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),(L(iag,iaa),iaa=1,naa)
      endif
      do 25 iaa=naa+1,naa0
25    L(iag,iaa)=L(iag,naa)
211   continue
      close(45)

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

      open(32,file='LIST.PRN')
      is=0
      open(31,file='FILE.X')
991   read(31,*,end=999,err=999)
      read(31,*,end=999,err=999)nat
      if(is.gt.0)deallocate(X,xt)
      allocate(xt(nat,3),X(nat,3))
      do 21 ia=1,nat
21    read(31,*)xt(ia,1),(xt(ia,i),i=1,3)
      is=is+1

c     number of active atoms:
      naa=nat-noff
c     number of atoms in molecule
      nam=naa
      if(nmol.ne.0)nam=nmol
c     number of molecules:
      nm=naa/nam
      if(mod(naa,nam).ne.0)call report('nm x nam <> naa')

      do 111 im=1,nm

      if(nmol.eq.0)then
       write(32,320)is
320    format(i6,$)
      else
       write(32,323)is,im
323    format(2i6,$)
      endif

      do 112 ia=1,nam
      do 112 i=1,3
112   X(ia,i)=xt(noff+ia+nam*(im-1),i)

      do 11 i=1,nag
      itype=abs(itypep(i))
      do 23 ix=1,3
      A(ix) =X(L(i,3),ix)-X(L(i,4),ix)
      P(ix) =X(L(i,1),ix)-X(L(i,2),ix)
23    E0(ix)=X(L(i,2),ix)-X(L(i,3),ix)
      call vp(B0,E0,A)
      SSB0B0=sp(B0,B0)
      SSE0E0=sp(E0,E0)
      DO 8 KK=1,3
      B0(KK)=B0(KK)/DSQRT(SSB0B0)
8     E0(KK)=E0(KK)/DSQRT(SSE0E0)
      SSPP=sp(P,P)
      SSAA=sp(A,A)
      SSAP=sp(A,P)
      AD=DSQRT(SSPP)
      SSAE0=sp(A,E0)
      SSPE0=sp(P,E0)
      DO 9 KK=1,3
      XX(KK)=A(KK)-E0(KK)*SSAE0
9     YY(KK)=P(KK)-E0(KK)*SSPE0
      SSXXXX=sp(XX,XX)
      SSYYYY=sp(YY,YY)
      SSXXYY=sp(XX,YY)
      call vp(AP,XX,YY)
      SSE0AP=sp(AP,E0)
      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(lper)then
       xper=xmax(i)-xmin(i)
c      write(6,*)i,ANGLE,xper,xmin(i),xmax(i)
1001   if(ANGLE.lt.xmin(i))then
        ANGLE=ANGLE+xper
        goto 1001
       endif
1002   if(ANGLE.gt.xmax(i))then
        ANGLE=ANGLE-xper
        goto 1002
       endif
      else
       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
      endif
      par=-1.0d10
c     2 - bond length:
      if(itype.eq.2)par=AD
c     3 - bond angle:
      if(itype.eq.3)par=BANGLE
c     4 - dihedral (torsional) angle:
      if(itype.eq.4)par=ANGLE
c     5 - out of plane deviation:
      if(itype.eq.5)par=sp(B0,P)
c     6 - deviation of two vectors:
      if(itype.eq.6)then
       par=0.0d0
       if(dabs(SSPP*SSAA).gt.0.0d0)par=dacos(SSAP/dsqrt(SSPP*SSAA))*
     1 180.0d0/PI
      endif
c     7 - deviation of two planes:
      if(itype.eq.7)then
       do 24 ix=1,3
       b(ix)=X(L(i,5),ix)-X(L(i,4),ix)
24     c(ix)=X(L(i,5),ix)-X(L(i,6),ix)
       call vp(v1,P,E0)
       call vp(v2,b,c)
       par=0.0d0
       if(dabs(sp(v1,v1)*sp(v2,v2)).gt.0.0d0)par=dacos(sp(v1,v2)
     1 /dsqrt(sp(v1,v1)*sp(v2,v2)))*180.0d0/PI
      endif
c     8,9 - puckering,amplitude of 5-member ring:
      if(itype.eq.8.or.itype.eq.9)then
c                 a1
c            t4  /  \ t0
c              a5    a2
c           t3 |     | t1
c              a4 -  a3
c                 t2
       t0=tau(nat,L(i,5),L(i,1),L(i,2),L(i,3),X)
       t1=tau(nat,L(i,1),L(i,2),L(i,3),L(i,4),X)
       t2=tau(nat,L(i,2),L(i,3),L(i,4),L(i,5),X)
       t3=tau(nat,L(i,3),L(i,4),L(i,5),L(i,1),X)
       t4=tau(nat,L(i,4),L(i,5),L(i,1),L(i,2),X)
       P0=atan2(t4+t1-t3-t0,2.0d0*t2*1.5388d0)
       tm=t2/cos(P0)
       P0=180.0d0/pi*P0
       if(P0.lt.0.0d0)P0=P0+360.0d0
       if(itype.eq.8)par=P0
       if(itype.eq.9)par=tm
      endif
c     distance beteen two lines given by atoms l1-l2 and l3-l4:
      if(itype.eq.10)then
c                l3--------------------l4
c                        |.
c                        |
c                        | d
c                        |.  
c          l1---------------------l2  
c     
       do 26 ix=1,3
       rab(ix)=X(L(i,3),ix)-X(L(i,1),ix)
       va( ix)=X(L(i,2),ix)-X(L(i,1),ix)
26     vb( ix)=X(L(i,4),ix)-X(L(i,3),ix)
       ao=dsqrt(sp(va,va))
       bo=dsqrt(sp(vb,vb))
       do 27 ix=1,3
       va0( ix)=va(ix)/ao
27     vb0( ix)=vb(ix)/bo
       sab=sp(va0,vb0)
       if(dabs(1.0d0-sab**2).lt.1.0d-4)then
        par=dsqrt(sp(rab,rab))
       else
        apar=0.0d0
        bpar=0.0d0
        do 28 ix=1,3
        apar=apar+rab(ix)*(va0(ix)-sab*vb0(ix))
28      bpar=bpar-rab(ix)*(vb0(ix)-sab*va0(ix))
        apar=apar/ao/(1.0d0-sab**2)
        bpar=bpar/bo/(1.0d0-sab**2)
        do 29 ix=1,3
29      u(ix)=rab(ix)+bpar*vb(ix)-apar*va(ix)
        par=dsqrt(sp(u,u))
       endif
      endif
      if(par.eq.-1.0d10)then
       write(6,*)itype
       call report('Unknown coordinate type')
      endif

      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
      if(itype.eq.2)then
       write(32,322)par
322    format(f10.5,$)
      else
       write(32,321)par
321    format(f10.3,$)
      endif
11    continue
      write(32,*)
111   continue

      goto 991
999   close(32)
      close(31)
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     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 readopt(noff,nmol,lper,ldistr)
      implicit none
      integer*4 noff,nmol
      logical lex,lper,ldistr
      character*4 key
      lper=.false.
      ldistr=.false.
      noff=0
      nmol=0
      inquire(file='XTORSION.OPT',exist=lex)
      if(lex)then
       write(6,600)
 600   format(' XTORSION.OPT found')
       open(9,file='XTORSION.OPT')
1      read(9,40,end=99,err=99)key
40     format(a4)
c      offset of atom numbers:
       if(key.eq.'NOFF')read(9,*)noff
c      periodicity in atom numbers:
c      (e.g. solvent of nmol atoms aroudn solute of noff atoms)
       if(key.eq.'NMOL')read(9,*)nmol
c      if ldistr (nag>0), use xper=xmax-amin to put values between:
       if(key.eq.'LPER')read(9,*)lper
c      make distribution (trigger also ny nag<0):
       if(key.eq.'DIST')read(9,*)ldistr
       goto 1
99     close(9)
      endif
      return
      end

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

      function sp(a,b)
      IMPLICIT none
      real*8 a(*),b(*),sp
      sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
      return
      end

      subroutine vp(A,X,Y)
      implicit none
      real*8 A(*),X(*),Y(*)
      A(1)=X(2)*Y(3)-X(3)*Y(2)
      A(2)=X(3)*Y(1)-X(1)*Y(3)
      A(3)=X(1)*Y(2)-X(2)*Y(1)
      return
      end

      function tau(nat,L1,L2,l3,L4,X)
c     torsion angle between -180 ... 180 deg
      implicit none
      integer*4 L1,L2,L3,L4,ix,KK,nat
      real*8 X(nat,3),A(3),P(3),E0(3),B0(3),sp,SSB0B0,SSE0E0,
     1SSAE0,SSPE0,XX(3),YY(3),SSXXXX,SSYYYY,SSXXYY,SSE0AP,AP(3),
     1SIGN,CA,PI,ANGLE,tau
      PI=3.141592653589D0
      do 23 ix=1,3
      A(ix) =X(L3,ix)-X(L4,ix)
      P(ix) =X(L1,ix)-X(L2,ix)
23    E0(ix)=X(L2,ix)-X(L3,ix)
      call vp(B0,E0,A)
      SSB0B0=sp(B0,B0)
      SSE0E0=sp(E0,E0)
      DO 8 KK=1,3
      B0(KK)=B0(KK)/DSQRT(SSB0B0)
8     E0(KK)=E0(KK)/DSQRT(SSE0E0)
      SSAE0=sp(A,E0)
      SSPE0=sp(P,E0)
      DO 9 KK=1,3
      XX(KK)=A(KK)-E0(KK)*SSAE0
9     YY(KK)=P(KK)-E0(KK)*SSPE0
      SSXXXX=sp(XX,XX)
      SSYYYY=sp(YY,YY)
      SSXXYY=sp(XX,YY)
      call vp(AP,XX,YY)
      SSE0AP=sp(AP,E0)
      SIGN=1.0D0
      IF (SSE0AP.NE.0.0D0)SIGN=-SSE0AP/ABS(SSE0AP)
      CA=0.0D0
      IF (SSXXXX.NE.0.0D0.AND.SSYYYY.NE.0.0D0)
     1CA=-SSXXYY/DSQRT(SSXXXX)/DSQRT(SSYYYY)
      ANGLE=0.0D0
      IF (CA.LE.-1.0D0) ANGLE=PI
      IF (ABS(CA).LT.1.0D0) ANGLE=ACOS(CA)
      ANGLE=180.0D0* ANGLE/PI
      ANGLE=SIGN*ANGLE
      if(ANGLE.lt.-180.0d0)ANGLE= 360.0d0+ANGLE
      if(ANGLE.gt.+180.0d0)ANGLE=-360.0d0+ANGLE
      tau=ANGLE
      return
      end