PROGRAM LIQUID_PAR
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      parameter (NAT0=10000,ndp0=900)
      character*80 s80
      DIMENSION X(NAT0),Y(NAT0),Z(NAT0),
     2it(NAT0),distr(ndp0)
      logical lex,lper,LANGLE,LDIST
      PI=3.141592653589D0
C     
      WRITE(6,*)
      WRITE(6,6000)
6000  format('Structural analysis of MD-modelled liquids ',/,/,
     1       'Geometries in FILE.X',/,/,
     2       'In LIQUID.OPT define:',/,
     3       '   ANGLE true - analyze angle',/,
     3       '   DISTA true - analyze distance',/,
     3       '   PER   true - consider periodic box',/,
     3       '   IJK   atomic numbers of interest, always',
     3       ' 3, even if last not used',/,
     3       '   XPER  periodic box dimension',/,
     3       '   YPER  periodic box dimension',/,
     3       '   ZPER  periodic box dimension',/,
     3       '   PMIN  parameter minimum value',/,
     3       '   PMAX  parameter maximum value',/,
     3       '   XMIN  distribution minimum value',/,
     3       '   XMAX  distribution maximum value',/,
     3       '   NDP   number of distribution points',/,
     3       '   DMIN  minimum distance value (also for angle)',/,
     3       '   DMAX  maximum distance value (also for angle)',/)
      inquire(file='LIQUID.OPT',exist=lex)
      if(.not.lex)call report(' LIQUID.OPT not found.')
      inquire(file='FILE.X',exist=lex)
      if(.not.lex)call report(' FILE.X not found.')

      LANGLE=.true.
      LDIST=.false.
      lper=.false.
      XPER=0.0d0
      YPER=0.0d0
      ZPER=0.0d0
      DMIN=0.7d0
      DMAX=2.4d0
      PMIN=80.0d0
      PMAX=160.0d0
      ndp=100
      XMIN=80.0d0
      XMAX=160.0d0
      I1=1
      J1=8
      K1=1
      open(45,file='LIQUID.OPT')
451   read(45,45)s80
45    format(a80)
      if(s80(1:5).eq.'ANGLE')read(45,*)LANGLE
      if(s80(1:4).eq.'DIST' )read(45,*)LDIST
      if(s80(1:3).eq.'PER'  )read(45,*)lper
      if(s80(1:4).eq.'XPER' )read(45,*)XPER
      if(s80(1:4).eq.'YPER' )read(45,*)YPER
      if(s80(1:4).eq.'ZPER' )read(45,*)ZPER
      if(s80(1:4).eq.'DMIN' )read(45,*)DMIN
      if(s80(1:4).eq.'DMAX' )read(45,*)DMAX
      if(s80(1:4).eq.'PMIN' )read(45,*)PMIN
      if(s80(1:4).eq.'PMAX' )read(45,*)PMAX
      if(s80(1:3).eq.'NDP'  )read(45,*)ndp
      if(s80(1:4).eq.'XMIN' )read(45,*)XMIN
      if(s80(1:4).eq.'XMAX' )read(45,*)XMAX
      if(s80(1:3).eq.'IJK'  )read(45,*)I1,J1,K1
      if(s80(1:3).eq.'END'  )goto 998
      goto 451
998   close(45)
      XPER2=XPER/2.0d0
      YPER2=YPER/2.0d0
      ZPER2=ZPER/2.0d0
      d2min=DMIN**2
      d2max=DMAX**2
      dd=(xmax-xmin)/dble(ndp)

      if(ndp.gt.ndp0)call report('too many distribution points')
      an=0.0d0
      do 10 i=1,ndp
10    distr(i)=0.0d0

      is=0
      open(31,file='FILE.X')
991   read(31,*,end=999,err=999)
      read(31,*,end=999,err=999)nat
      if(nat.gt.NAT0)call report('Too many atoms')
      do 21 ia=1,nat
21    read(31,*)it(ia),X(ia),Y(ia),Z(ia)
      is=is+1

      do 211 ia=1,nat
      if(it(ia).eq.I1)then
       xa=X(ia)
       ya=Y(ia)
       za=Z(ia)
       do 212 ib=1,nat
       if(it(ib).eq.J1.and.ib.ne.ia)then
        xb=X(ib)
        yb=Y(ib)
        zb=Z(ib)
        if(lper)then
         if(xb.gt.xa+XPER2)xb=xb-XPER
         if(yb.gt.ya+YPER2)yb=yb-YPER
         if(zb.gt.za+ZPER2)zb=zb-ZPER
         if(xb.lt.xa-XPER2)xb=xb+XPER
         if(yb.lt.ya-YPER2)yb=yb+YPER
         if(zb.lt.za-ZPER2)zb=zb+ZPER
        endif
        d2=(xa-xb)**2+(ya-yb)**2+(za-zb)**2
        if(d2.gt.d2min.and.d2.lt.d2max)then
         if(LDIST)then       
          par=dsqrt(d2)
          if(par.le.PMAX.and.par.ge.PMIN)then
           an=an+1.0d0
           j=int((par-xmin+dd)/dd)
           if(j.lt.1)j=1
           if(j.gt.ndp)j=ndp
           distr(j)=distr(j)+1.0d0
          endif
         endif
         if(LANGLE)then
          do 213 ic=1,nat
          if(it(ic).eq.K1.and.ic.gt.ia)then
           xc=X(ic)
           yc=Y(ic)
           zc=Z(ic)
           if(lper)then
            if(xc.gt.xb+XPER2)xc=xc-XPER
            if(yc.gt.yb+YPER2)yc=yc-YPER
            if(zc.gt.zb+ZPER2)zc=zc-ZPER
            if(xc.lt.xb-XPER2)xc=xc+XPER
            if(yc.lt.yb-YPER2)yc=yc+YPER
            if(zc.lt.zb-ZPER2)zc=zc+ZPER
           endif
           d2=(xc-xb)**2+(yc-yb)**2+(zc-zb)**2
           if(d2.gt.d2min.and.d2.lt.d2max)then
            an=an+1.0d0
            PX=xa-xb
            PY=ya-yb
            PZ=za-zb
            EX=xb-xc
            EY=yb-yc
            EZ=zb-zc
            CB=-(PX*EX+PY*EY+PZ*EZ)/DSQRT(PX*PX+PY*PY+PZ*PZ)
     1      /dsqrt(EX*EX+EY*EY+EZ*EZ)
            IF(CB.LE.-1.0D0)then
             BANGLE=PI
            ELSE
             IF(CB.LT.1.0D0)THEN
              BANGLE=ACOS(CB)
             ELSE
              BANGLE=0.0D0
             ENDIF
            ENDIF
            par=180.0D0*BANGLE/PI
            if(par.le.PMAX.and.par.ge.PMIN)then
             an=an+1.0d0
             j=int((par-xmin+dd)/dd)
             if(j.lt.1)j=1
             if(j.gt.ndp)j=ndp
             distr(j)=distr(j)+1.0d0
            endif
           endif
          endif
213       continue
         endif
        endif
       endif
212    continue
      endif
211   continue
c
c     goto next geometry:
      goto 991
999   close(31)
c
      write(6,6009)is,an
6009  format(i10,' geometries',/,
     1       f10.0,' parameters')
      open(33,file='DISTR.PRN')
      fac=1.0d0/an/dd
      do 17 j=1,ndp
17    write(33,3333)xmin+dd*dble(2*j-1)/2.0d0,distr(j)*fac
3333  format(2f12.6)
      close(33)
      call report('Distribution file DISTR.PRN written')
      end

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