program ana3x
      implicit none
      integer*4 np,ig,i,n1,n2,nm,im,j,ii,m,ichir(5),ix,ir,is,iin,ic,isu
      integer*4 nat
      real*8 xmin,xmax,px(3),dd,cx(3),dx,v23(3),v34(3),v15(3),u(3),
     1sp,xmi,xma,cau(3),v1(3),v2(3),surf,dm,di,cg(3)
      real*8,allocatable::d(:),rm(:),dr(:),ds(:),xau(:)
      character*80 fn
      logical lchir,lsurf
      open(9,file='ANA.OPT',status='old')
      read(9,*)np,xmin,xmax
      read(9,*)px
      read(9,*)n1,n2,nm
      read(9,*)lchir
      read(9,*)ichir
      read(9,*)lsurf
      read(9,*)surf
      close(9)
      m=(n2-n1+1)/nm
      dx=(xmax-xmin)/np
      allocate(d(np),dr(np),ds(np),rm(3*nm))
      d=0.0d0
      dr=0.0d0
      ds=0.0d0
      write(6,*)m,' molecules'

      ig=0
      open(9,file='FILE.X')
1     read(9,900,end=99,err=99)fn
900   format(a80)
      read(9,*,end=99,err=99)nat
      ig=ig+1
      if(ig.eq.1)allocate(xau(3*nat))
c     geometry ig:      
      ir=0
      is=0
      isu=0
c     iin.. number of molecules inside the wire:
      iin=0

c     geometry
      do 3 i=1,nat
3     read(9,*)xau(1+3*(i-1)),(xau(ix+3*(i-1)),ix=1,3)

c     center of this geometry (only for the first one):
      if(ig.eq.1)then
       do 31 ix=1,3
       cg(ix)=0.0d0
       do 32 i=1,nat
32     cg(ix)=cg(ix)+xau(ix+3*(i-1))
31     cg(ix)=cg(ix)/dble(nat)
      endif

c     mass centers of m molecules
      do 4 im=1,m
c     cx: mass center of molecule
      cx=0.0d0
c     xmi,xma: x-ranges for molecular atoms:
      xmi=0.0d0
      xma=0.0d0
      do 5 j=1,nm
      do 71 ix=1,3
71    rm(ix+3*(j-1))=xau(ix+3*(j-n1+nm*(im-1)))
      if(j.eq.1)then
       xmi=rm(1)
       xma=rm(1)
      endif
      if(rm(1+3*(j-1)).lt.xmi)xmi=rm(1+3*(j-1))
      if(rm(1+3*(j-1)).gt.xma)xma=rm(1+3*(j-1))
      do 5 ix=1,3
5     cx(ix)=cx(ix)+rm(ix+3*(j-1))
      do 51 ix=1,3
      cx(ix)=cx(ix)/dble(nm)
      if(cx(ix)-cg(ix).gt. px(ix)/2.0d0)cx(ix)=cx(ix)-px(ix)
51    if(cx(ix)-cg(ix).lt.-px(ix)/2.0d0)cx(ix)=cx(ix)+px(ix)

c     check that the molecule is not inside:      
c     cau: center of close gold atomsL
      cau=0.0d0
c     close Au atoms along x:
      ic=0
      do 31 i=n2+1,nat
      if(xau(1+3*(i-1)).gt.xmi-2.0d0.and.
     1   xau(1+3*(i-1)).lt.xma+2.0d0)then
       ic=ic+1
       do 311 ix=1,3
311    cau(ix)=cau(ix)+xau(ix+3*(i-1))
      endif
31    continue
      do 313 ix=1,3
313   if(ic.gt.0)cau(ix)=cau(ix)/dble(ic)
      ic=0
      do 314 ix=1,3
314   v2(ix)=cau(ix)-cx(ix)
      do 32 i=1,n1-1
      if(xau(1+3*(i-1)).gt.xmi-2.0d0.and.
     1   xau(1+3*(i-1)).lt.xma+2.0d0)then
       do 312 ix=1,3
312    v1(ix)=xau(ix+3*(i-1))-cx(ix)
       if(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3).lt.0.0d0)ic=ic+1
      endif
32    continue

      if(lsurf)then
c      check closest distance between molecule and gold:
       dm=1.0d10
       do 33 i=n2+1,nat
       di=(xau(1+3*(i-1))-cx(1))**2+(xau(2+3*(i-1))-cx(2))**2
     1   +(xau(3+3*(i-1))-cx(3))**2
33     if(di.lt.dm)dm=di
       if(dsqrt(dm).gt.surf)goto 4
      endif
      

      if(ic.gt.0)then
c      write(6,602)im
c02    format(' Molecule',i6,' inside')
       iin=iin+1
c      write(6,603)cx
c03    format(' cx  :',f10.2,' cy  :',f10.2,' cz  :',f10.2)
      else
c      x-position:
       ii=int((cx(1)-xmin)/dx)+1
       if(ii.ge.1.and.ii.le.np)d(ii)=d(ii)+1.0d0
       if(lchir)then
        do 8 ix=1,3
        v23(ix)=rm(ix+3*(ichir(3)-1))-rm(ix+3*(ichir(2)-1))
        v34(ix)=rm(ix+3*(ichir(4)-1))-rm(ix+3*(ichir(3)-1))
        v15(ix)=rm(ix+3*(ichir(5)-1))-rm(ix+3*(ichir(1)-1))
8       continue
        call vp(u,v23,v34)
        if(sp(u,v15).gt.0.0d0)then
c        chirality "R":
         if(ii.ge.1.and.ii.le.np)dr(ii)=dr(ii)+1.0d0
         ir=ir+1
        else
c        chirality "S":
         if(ii.ge.1.and.ii.le.np)ds(ii)=ds(ii)+1.0d0
         is=is+1
        endif
       endif
      endif
4     continue
      write(6,600)fn(1:40),ir,is,iin
600   format(a40,1x,i6,' R,',i6,' S',i6,' inside')

      goto 1
99    close(9)
      write(6,*)ig,' geometries'

      dd=0.0d0
      do 6 i=1,np
6     dd=dd+d(i)*dx
      do 7 i=1,np
      dr(i)=dr(i)/dd
      ds(i)=ds(i)/dd
7     d(i)=d(i)/dd

      call wd('di.prn',xmin,dx,np,d)
      if(lchir)then
       call wd('dr.prn',xmin,dx,np,dr)
       call wd('ds.prn',xmin,dx,np,ds)
      endif

      end

      subroutine wd(f,xmin,dx,np,d)
      implicit none
      character*(*) f
      real*8 xmin,dx,d(*),x
      integer*4 np,i
      open(9,file=f)
      x=xmin-dx/2.0d0
      do 2 i=1,np
      x=x+dx
2     write(9,902)x,d(i)
902   format(2e12.4)
      close(9)
      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 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