program rau
      implicit none
      real*8 AVO,qw,mw,px,py,pz,dmi,V,rom,dm,dd,amasf
      integer*4 nat,ifa,natu,nats,ia,nmol,iq,nd
      logical lper,lm,ladd
      real*8,allocatable::r(:)
      integer*4,allocatable::iza(:)
      AVO=6.022045d23
c
      write(6,6000)
6000  format(//,' Solvent distribution above a solute surface',/,/,
     1' RAU.OPT    ... option file',/,
     1' FILE.X     ... geometries' ,/,
     1' SOLU.LST   ... list of solute atoms ',/)

      call readopt(ifa,px,py,pz,dmi,dm,nd,lper,dd,natu,nats,lm,ladd)

      open(2,file='FILE.X',status='old')
      read(2,*)
      read(2,*)nat
      allocate(r(3*nat),iza(nat))
      qw=0.0d0
      mw=0.0d0
      do 1 ia=1,nat
      read(2,*)iq
      qw=qw+dble(iq)
1     mw=mw+amasf(iq)
      close(2)
      V=px*py*pz*1.0d-24
      rom=mw/V/AVO
      nmol=(nat-natu)/nats
      write(6,600)nat,nmol,rom
600   format(i6,' atoms,',i6,' molecules,',f10.2,' g/mol')

      call wlist(nat,iza,nd,natu,ifa,lper,r,px,py,pz,dd,dmi,ladd)
      if(lm)call wmol(nat,iza,nd,natu,lper,r,px,py,pz,dd,dmi,nmol,
     1nats,ladd)

      end

      subroutine wmol(nat,iza,nd,natu,lper,r,px,py,pz,dd,dmi,
     1nmol,nats,ladd)
      implicit none
      integer*4 nd,natu,nmc,ii,nat,ia,iza(*),ni,nmol,nats,im,
     1ix,iy,iz,ic,i,ir,sm(3),jj,ja
      real*8 r(*),xa,ya,za,x,y,z,px,py,pz,dd,dmi,di,da,nn,pi,sdd,
     1amasf,rm(3),wm,mm,xi,yi,qi,qj,xii,yii,zii,
     1xci,yci,zci,xcj,ycj,zcj,sx(3),sy(3),sz(3),sp,
     1xj,yj,si,tm,ri,rj,rij2
      logical lper,ladd
      real*8,allocatable::ra(:),rb(:),ci(:)
      integer*4,allocatable::k(:)
      allocate(ra(nd),rb(nd),ci(natu),k(natu))
      call vz(ra,nd)
      call vz(rb,nd)
      call vz(ci,natu)
      k=0
      pi=4.0d0*atan(1.0d0)
      nn=0.0d0
      sdd=1.0d8
      nmc=0
c     if per. box, scan 3x3x3 superbox:
      if(lper)then
       ir=27
      else
       ir=1
      endif
      open(2,file='FILE.X',status='old')
1     continue

      read(2,*,end=500,err=500)
      read(2,*,end=500,err=500)nat
      do 2 ia=1,nat
2     read(2,*)iza(ia),(r(3*(ia-1)+ix),ix=1,3)
      nmc=nmc+1
c
      ix=-1
      iy=-1
      iz=-1
      do 61 ic=1,ir
      do 6 im=1,nmol
      rm=0.0d0
      wm=0.0d0
      do 62 ia=natu+1+(im-1)*nats,natu+im*nats
      mm=amasf(iza(ia))
      wm=wm+mm
      do 62 jj=1,3
62    rm(jj)=rm(jj)+mm*r(3*(ia-1)+jj)
      do 63 jj=1,3
63    rm(jj)=rm(jj)/wm

      x=rm(1)+px*dble(ix)
      y=rm(2)+py*dble(iy)
      z=rm(3)+pz*dble(iz)

c     closest atom and distance from solute:
      di=1.0d8
      ii=0
      do 612 i=1,natu
      xa=r(3*(i-1)+1)-x
      ya=r(3*(i-1)+2)-y
      za=r(3*(i-1)+3)-z
      da=xa**2+ya**2+za**2
      if(da.lt.di)then
       di=da
       ii=i
      endif
612   continue
      xii=r(3*(ii-1)+1)
      yii=r(3*(ii-1)+2)
      zii=r(3*(ii-1)+3)
      sz(1)=x-xii
      sz(2)=y-yii
      sz(3)=z-zii
      call vnorm(sz)

c     some other atom
      do 613 i=1,natu
      if(i.ne.ii)then
       sx(1)=r(3*(i-1)+1)-xii
       sx(2)=r(3*(i-1)+2)-yii
       sx(3)=r(3*(i-1)+3)-zii
       call vp(sz,sx,sy)
       if(sp(sy,sy).lt.1.0d-30)goto 613
       call vnorm(sy)
       call vp(sy,sz,sx)
       call vnorm(sx)
       goto 614
      endif
613   continue
      call report('plane vectors not found')
c
c     plane chirality index:
614   si=0.0d0
      tm=0.0d0
      do 64 ia=natu+1+(im-1)*nats,natu+im*nats
c     mass center coordinates      
      xci=r(3*(ia-1)+1)-rm(1)
      yci=r(3*(ia-1)+2)-rm(2)
      zci=r(3*(ia-1)+3)-rm(3)
c     plane coordinates:
      xi=xci*sx(1)+yci*sx(2)+zci*sx(3)
      yi=xci*sy(1)+yci*sy(2)+zci*sy(3)
      ri=xi**2+yi**2
      qi=dble(iza(ia))
      tm=tm+qi**2
      do 64 ja=natu+1+(im-1)*nats,natu+im*nats
      xcj=r(3*(ja-1)+1)-rm(1)
      ycj=r(3*(ja-1)+2)-rm(2)
      zcj=r(3*(ja-1)+3)-rm(3)
      xj=xcj*sx(1)+ycj*sx(2)+zcj*sx(3)
      yj=xcj*sy(1)+ycj*sy(2)+zcj*sy(3)
      rj=xj**2+yj**2
      qj=dble(iza(ja))
      rij2=(xi-xj)**2+(yi-yj)**2
64    if(rj.gt.ri.and.ia.ne.ja)si=si+qi*qj*(xi*yj-yi*xj)/rij2
      si=si/tm

 
      ni=int((sqrt(di)-dmi)/dd)+1
      if(ni.le.nd.and.ni.ge.1)ra(ni)=ra(ni)+1.0d0/di
      if(ni.le.nd.and.ni.ge.1)rb(ni)=rb(ni)+si/di
      ci(ii)=ci(ii)+si/di
      k(ii)=k(ii)+1
      nn=nn+1.0d0
      if(sdd.gt.di)then
       sdd=di
       sm(1)=nmc
       sm(2)=im
       sm(3)=ii
      endif
6     continue
c     loop over molecules

      ix=ix+1
      if(ix.gt.1)then
       ix=-1
       iy=iy+1
       if(iy.gt.1)then
        iy=-1
        iz=iz+1
       endif
      endif
61    continue
c     loop over periodic boxes
c       
      goto 1
500   close(2)
c
      write(6,800)nn
800   format(f20.0,' instances of counted molecules')
      write(6,801)sdd,sm(1),sm(2),sm(3)
801   format(' Smallest distance:',f6.2,':',
     1' Geometry:',i5,', Molecule:',i5,', Atom:',i6) 

      do 103 i=1,nd
      ra(i)=ra(i)*px*py*pz/(4.0d0*pi*dd*nn)
103   rb(i)=rb(i)*px*py*pz/(4.0d0*pi*dd*nn)
c
      call wrf('mrdf.prn',nd,dmi,ra,dd,nmc,ladd)
      call wrf('brdf.prn',nd,dmi,rb,dd,nmc,ladd)

      open(8,file='SOLI.X')
      write(8,*)'solute with chirality powers'
      write(8,*)natu
      do 104 ia=1,natu
      if(k(ia).ne.0)ci(ia)=ci(ia)/dble(k(ia))
104   write(8,802)iza(ia),(r(3*(ia-1)+ix),ix=1,3),ci(ia)
802   format(i6,3f12.6,' 0 0 0 0 0 0 0',f8.4)
      close(8)

      return
      end

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

      subroutine wlist(nat,iza,nd,natu,ifa,lper,r,px,py,pz,dd,dmi,
     1ladd)
      implicit none
      integer*4 nd,natu,nmc,ii,nat,ia,iza(*),ni,
     1ix,iy,iz,ic,i,iq,ifa,is2,ie2,ir,sm(3)
      character*80 sn2
      real*8 r(*),xa,ya,za,x,y,z,px,py,pz,dd,dmi,di,da,nn,pi,sdd
      logical lper,ladd
      real*8,allocatable::ra(:)
      allocate(ra(nd))
      call vz(ra,nd)
      pi=4.0d0*atan(1.0d0)
      nn=0.0d0
      sdd=1.0d8
      nmc=0
c     if per. box, scan 3x3x3 superbox:
      if(lper)then
       ir=27
      else
       ir=1
      endif
      open(2,file='FILE.X',status='old')
1     continue

      read(2,*,end=500,err=500)
      read(2,*,end=500,err=500)nat
      do 2 ia=1,nat
2     read(2,*)iza(ia),(r(3*(ia-1)+ix),ix=1,3)
      nmc=nmc+1
c
c     surface for solute atoms:
c     call dsv(s,natu
c
      ix=-1
      iy=-1
      iz=-1
      do 61 ic=1,ir
      do 6 ia=natu+1,nat
      iq=iza(ia)
      x=r(3*(ia-1)+1)+px*dble(ix)
      y=r(3*(ia-1)+2)+py*dble(iy)
      z=r(3*(ia-1)+3)+pz*dble(iz)

c     specific atom type or all:
      if(ifa.eq.iq.or.ifa.eq.0)then
c      shortest distance from the solute:
       di=1.0d8
       ii=0
       do 612 i=1,natu
       xa=r(3*(i-1)+1)-x
       ya=r(3*(i-1)+2)-y
       za=r(3*(i-1)+3)-z
       da=xa**2+ya**2+za**2
       if(da.lt.di)then
        di=da
        ii=i
       endif
612    continue
       ni=int((sqrt(di)-dmi)/dd)+1
       if(ni.le.nd.and.ni.ge.1)ra(ni)=ra(ni)+1.0d0/di
       nn=nn+1.0d0
       if(sdd.gt.di)then
        sdd=di
        sm(1)=nmc
        sm(2)=ia
        sm(3)=ii
       endif
      endif
6     continue
c     loop over atoms

      ix=ix+1
      if(ix.gt.1)then
       ix=-1
       iy=iy+1
       if(iy.gt.1)then
        iy=-1
        iz=iz+1
       endif
      endif
61    continue
c     loop over periodic boxes
c       
      goto 1
500   close(2)
c
      write(6,*)nmc,' geometries'
      write(6,800)nn
800   format(f20.0,' instances of counted atoms')
      write(6,801)sdd,sm(1),sm(2),sm(3)
801   format(' Smallest distance:',f10.2,':',
     1' Geometry:',i6,', atom:',i6,', solute atom:',i6) 

      do 103 i=1,nd
103   ra(i)=ra(i)*px*py*pz/(4.0d0*pi*dd*nn)
c
      write(sn2,309)ifa
309   format(i80)
      call dse(sn2,is2,ie2)
      call wrf(sn2(is2:ie2)//'rdf.prn',nd,dmi,ra,dd,nmc,ladd)

      return
      end

      subroutine wrf(sn,nd,dmi,r,dd,nmc,ladd)
      implicit none
      integer*4 nd,nmc,nmco,i
      real*8 r(*),dmi,dd,a,ao,s,x
      character*(*)sn
      logical lex,ladd
      write(6,*)sn
      if(ladd)then
       inquire(file=sn//'.N',exist=lex)
       open(3,file=sn//'.N')
       if(lex)then
        read(3,*)nmco
        rewind 3
        write(3,*)nmco+nmc
        close(3)

        a=dble(nmc)
        ao=dble(nmco)
        s=a+ao
        open(3,file=sn)
        do 2 i=1,nd
        read(3,*)x,x
2       r(i)=(r(i)*a+x*ao)/s
        close(3)
        write(6,602)nmco,nmc
602     format(i10,' old geometries,',i10,' new ones')
       else
        nmco=0
        write(3,*)nmc
        close(3)
        write(6,603)
603     format(' starting averaging')
       endif

      endif

      open(3,file=sn)
      do 1 i=1,nd
1     write(3,4000)dmi+dd/2.0d0+dble(i-1)*dd,r(i)
4000  format(f10.5,' ',g12.6)
      close(3)
      return
      end

      function amasf(q)
      implicit none
      real*4 amas(89)
      real*8 amasf
      integer*4 q
      data amas/1.008,4.003,
     2  6.941, 9.012,   10.810,12.011,14.007,15.999,18.998,20.179,
     3 22.990,24.305,   26.981,28.086,30.974,32.060,35.453,39.948,
     4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847,
     4 58.933,58.700,63.546,65.380,
     4                  69.720,72.590,74.922,78.960,79.904,83.800,
     5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070,
     5 102.906,106.400,107.868,112.410,
     5                 114.82,118.69,121.75,127.600,126.905,131.300,
     6 132.905,137.330,138.906,
     6                 140.120,140.908,144.240,145.000,150.400,
     6 151.960,157.250,158.925,162.500,164.930,167.260,168.934,
     6 173.040,174.970,
     6 178.490,180.948,183.850,186.207,190.200,192.220,195.090,
     6 196.967,207.590,204.370,207.200,208.980,210.000,210.001,
     6 222.02,
     7 223.000,226.025,227.028/     
      amasf=dble(amas(q))
      return
      end

      subroutine vz(v,n)
      real*8 v(*)
      integer*4 n,i
      do 1 i=1,n
1     v(i)=0.0d0
      return
      end

      subroutine dse(sn,is,ie)
      character*(*) sn
      integer*4 is,ie
      do 11 is=1,len(sn)
11    if(sn(is:is).ne.' ')goto 12
12    do 13 ie=len(sn),1,-1
13    if(sn(ie:ie).ne.' ')goto 14
14    return
      end


      subroutine readopt(ifa,px,py,pz,dmi,dm,nd,lper,dd,natu,nats,lm,
     1ladd)
      implicit none
      integer*4 ifa,nd,natu,nats
      real*8 px,py,pz,dmi,dm,dd
      logical lper,lm,ladd
      character*80 s80,s81
      nats=0
      natu=0
      ifa=0
      px=0.0d0
      py=0.0d0
      pz=0.0d0
      dmi=0.0d0
      dm=10.0d0
      nd=200
      lper=.false.
      ladd=.false.
      lm=.false.
      open(9,file='RAU.OPT',status='old')
1000  read(9,900,end=99,err=99)s80
900   format(a80)
c     minimal distance:
      if(s80(1:4).eq.'RMIN')read(9,*)dmi
c     maximal distance:
      if(s80(1:4).eq.'RMAX')read(9,*)dm
c     number of poins:
      if(s80(1:2).eq.'NP')read(9,*)nd
c     follow atom, zero for all:
      if(s80(1:6).eq.'FOLLOW')read(9,*)ifa
c     number of solute atoms:
      if(s80(1:4).eq.'NSOL')read(9,*)natu
c     number of atoms in molecule:
      if(s80(1:4).eq.'NATS')read(9,*)nats
c     periodic boundary conditions:
      if(s80(1:4).eq.'LPER')read(9,*)lper
c     periodic box dimensions:
      if(s80(1:4).eq.'PERI')read(9,*)px,py,pz
c     add to previous distributions:
      if(s80(1:4).eq.'LADD')read(9,*)ladd
c     chirality index:
      if(s80(1:4).eq.'CHIR')read(9,*)lm
      if(s80(1:3).eq.'END')goto 99
      backspace 9
      read(9,900)s81
      if(s80.eq.s81)then
       write(6,900)s80
       write(6,*)'unknown option'
       stop
      endif
      goto 1000
99    close(9)
      dd=(dm-dmi)/dble(nd)
      write(6,6002)dd
6002  format(' Tick: ',f6.3,' A')
      if(.not.lper)call report('Lper not defined')
      if(nats.eq.0)call report('NATS not defined')
      if(natu.eq.0)call report('NSOL not defined')
      if(px.eq.0.0d0)call report('px=0')
      return
      end

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

      subroutine vp(a,b,c)
      implicit none
      real*8 a(3),b(3),c(3)
      c(1)=a(2)*b(3)-a(3)*b(2)
      c(2)=a(3)*b(1)-a(1)*b(3)
      c(3)=a(1)*b(2)-a(2)*b(1)
      return
      end



      subroutine vnorm(v)
      real*8 v(3),n,sp
      n=dsqrt(sp(v,v))
      v(1)=v(1)/n
      v(2)=v(2)/n
      v(3)=v(3)/n
      return
      end