program dynalook
      implicit none
      integer*4 nat0
      parameter (nat0=100000)
      real*8 r(nat0,3),dist,tstep,p(3),cm(3),cs(3),rd(nat0),dj
      integer*4 isu(nat0),natu,natsolv,nat,as(nat0),nin,nout,nouttot,
     1nintot,list(nat0),list1(nat0),listold(nat0),nlist1,nlist,nlistold,
     1ix,is,k,i,ii,j,nav
      character*50 st
      character*80 s80
      logical lper,lopt,lwr

      lper=.false.
      lwr=.false.
      dist=3.6
      tstep=0.1d0
      nlistold=0
      do 16 i=1,nlistold
16    listold(i)=0
      write(6,600)
600   format(' Looks at changes in a hydration shell, in FILE.X',/,/,
     1' DLOOK.OPT  ....   list of options')
c
c     Options in DLOOK.OPT:
c
c     TSTEP
c     t time step in ps separating the snapshots at FILE.X
c     PERIODIC
c     lper  true for periodic
c     a b c the dimensions of the periodic box
c     SOLUTE
c     nats 
c     is1 ...isnat     number of solute atoms and their list
c     SOLVENT
c     natsolv  number of atoms in ONE solvent molecule
c              regular numbering supposed except for the solute !!
c     DIST
c     dist   dimension of the first solvation shell in A
c     LWR
c     lwr    detailed output
c
      inquire(file='DLOOK.OPT',exist=lopt)
      if(lopt)then
       write(6,*)'DLOOK.OPT found'
       open(44,file='DLOOK.OPT')
101    read(44,440,end=102,err=102)st
440    format(a50)
       if(st(1:5).eq.'TSTEP')read(44,*)tstep
       if(st(1:6).eq.'SOLUTE')then
        read(44,*)natu
        if(natu.gt.nat0)call report('Too many atoms')
        read(44,*)(isu(i),i=1,natu)
       endif
       if(st(1:3).eq.'PER')then
        read(44,*)lper
        read(44,*)(p(ix),ix=1,3)
       endif
       if(st(1:7).eq.'SOLVENT')read(44,*)natsolv
       if(st(1:3).eq.'LWR')read(44,*)lwr
       if(st(1:4).eq.'DIST')read(44,*)dist
       goto 101
102    close(44)
      else
       call report('DLOOK.OPT not found')
      endif
c
c     total change in the solvation shell and in the vicinity of solute atoms:
      open(10,file='FLIP.LST')
c     distances of initial solvent molecules in the shell from the solvent:
      open(11,file='DIST.PRN')
      is=0
      nav=0
      nintot=0
      nouttot=0
      open(3,file='FILE.X',status='old')
1     read(3,300,end=999,err=999)s80
300   format(a80)
      read(3,*,end=999,err=999)nat
      if(nat.gt.nat0)call report('too many atoms')
      is=is+1
      do 2 i=1,nat
2     read(3,*)as(i),(r(i,k),k=1,3)

c     center of solute molecule:
      do 7 ix=1,3
      cs(ix)=0.0d0
      do 8 j=1,natu
8     cs(ix)=cs(ix)+r(isu(j),ix)
7     cs(ix)=cs(ix)/dble(natu)

c     determine solvent molecules in the shell (nsm):
      nlist=0
      k=0
      i=0
c     loop aver all atoms:
3     i=i+1
c     skip solute atoms:
      do 4 j=1,natu
4     if(isu(j).eq.i)goto 3
c     center of solvent molecule number k:
      k=k+1
      do 5 ix=1,3
      cm(ix)=0.0d0
      do 6 ii=i,i+natsolv-1
6     cm(ix)=cm(ix)+r(ii,ix)
5     cm(ix)=cm(ix)/dble(natsolv)
c     correct to periodicity:
      if(lper)then
       do 9 ix=1,3
       if(cm(ix).gt.cs(ix)+p(ix)/2.0d0)cm(ix)=cm(ix)-p(ix)
9      if(cm(ix).lt.cs(ix)-p(ix)/2.0d0)cm(ix)=cm(ix)+p(ix)
      endif
c     distance of molecule k from the solute
      rd(k)=10000.0d0
      do 17 j=1,natu
      dj=dsqrt((cm(1)-r(isu(j),1))**2+(cm(2)-r(isu(j),2))**2
     1+(cm(3)-r(isu(j),3))**2)
17    if(dj.lt.rd(k))rd(k)=dj
c     if in shell, add to list:
      if(rd(k).lt.dist)then
       nlist=nlist+1
       list(nlist)=k
c      write(6,*)k,i,i+natsolv-1,rd(k)
c      write(6,*)(cm(ix),ix=1,3)
c      write(6,*)(r(i,ix),ix=1,3)
c      read(5,*)ii
      endif
      i=i+natsolv-1
      if(i.lt.nat)goto 3

      nav=nav+nlist
      if(lwr)write(6,6001)is,nlist,(list(i),i=1,nlist)
6001  format(' step',i7,',',i3,' shell molecules:',200i4)

c     remember initial list:
      if(is.eq.1)then
       nlist1=nlist
       do 10 i=1,nlist
10     list1(i)=list(i)
      endif
c     take distances of initial molecules and write to DIST.PRN:  
      write(11,1110)is,(rd(list1(i)),i=1,nlist1)
1110  format(i9,40f10.4)

c     compare with old list:
c     nin ... number of new molecules in the shell
c     nout ... number of molecules that left the shell
      nin=0
      do 11 i=1,nlist
      do 12 j=1,nlistold
12    if(listold(j).eq.list(i))goto 11
      nin=nin+1
11    continue
      nout=0
      do 13 i=1,nlistold
      do 14 j=1,nlist
14    if(listold(i).eq.list(j))goto 13
      nout=nout+1
13    continue
      write(10,1000)is,nin,nout
1000  format(6i6)
      nintot=nintot+nin
      nouttot=nouttot+nout

c     save current list
      nlistold=nlist
      do 15 i=1,nlist
15    listold(i)=list(i)

      goto 1

999   close(3)
      close(10)
      close(11)
      write(6,6003)is,nintot,nouttot,dble(nintot)/dble(is)/tstep,
     1dble(nav)/dble(is),
     11.0d0/(dble(nintot)/dble(is)/tstep/dble(nav)*dble(is)),
     1tstep*dble(is),dist
6003  format(i8,' geometries,'i6,' in',i6,' out ',f10.4,'ps-1',
     1/,' average size of the shell :',f10.1,
     1/,' average resident time     :',f10.4,' ps',
     1/,' total simulation time     :',f10.4,' ps',
     1/,' shell radius              :',f10.4,' A')

      stop
      end

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