program rdf3d
      implicit none
      integer*4 np0,nat0,i,ia,ic,icvar,ie2,ifa,iq,ir,is2,ix,iy,iz,
     1ja,jh,jm,jo,nh,nh1,nh2,nmc,nwm,nwmo,j,ica,ici,ii
      parameter (np0=10000000,nat0=10000)
      character*80 fn,sn2
      character*3 sy
      character*5 so
      real*8 rh(np0),ro(np0),r(3*nat0),rm(np0),r0(3),A(3,3),ro0(np0),
     1dr(3),ra(3*nat0),d2,doh,px(3),py(3),pz(3),vol,x,xj,xw,xw1,
     1xw2,y,yj,yw,yw1,yw2,z,zw,zw1,zw2,zj,cuti,ax,ay,az,
     1ppx,ppy,ppz,diff,da,db,cuta,diff2
      integer*4 nd(3),idumm,nc,alist(nat0),nat,nats,iza(nat0),
     1izq(nat0),iaa
      logical lex,lex2,lex3,lper,lsmooth,lamat,lex4,lanal,lskip
c
      write(6,6000)
6000  format(/,'     3D probability distribution',/,/,
     1         '       CAGE.PAR ... cage parameters',/,
     1         '       RDF3D.OPT... list of options',/)
c
      inquire(file='CAGE.PAR',exist=lex)
      if(.not.lex)call report('CAGE.PAR does not exist')
      open(9,file='CAGE.PAR')
      read(9,*)idumm,(r0(i),i=1,3)
      read(9,*)nd(1),dr(1)
      read(9,*)nd(2),dr(2),dr(2)
      read(9,*)nd(3),dr(3),dr(3),dr(3)
      close(9)
      nc=nd(1)*nd(2)*nd(3)
      write(6,6001)nd(1),nd(2),nd(3),nc
6001  format('CAGE.PAR:',i4,' x',i4,' x',i4,' = ',i7,' grid',/)
      if(nc.gt.np0)call report('Too many points')

c     defaults:
c     look at water(otherwise atom type):
c     ifa=0  look for water
c     ifa>0  look for this atom number
c     ifa=-1 look for mass density
c     ifa=-2 look for charge density
c     ifa=-3 look for atomic density
      ifa=0
c     no cut off distances max min:
      cuta=0.0d0
      cuti=0.0d0
c     no periodic conditions:
      lper=.false.
c     no transformation matrices:
      lamat=.false.
c     no post-analysis:
      lanal=.false.
c     use FILE.X:
      fn='FILE.X'
c     sooth the distribution:
      lsmooth=.true.
c     skip solute atoms fro the density:
      lskip=.true.
c     atoms listed in the cube:
      nats=0

      inquire(file='RDF3D.OPT',exist=lex)
      if(lex)then
       write(6,*)'RDF3D.OPT:'
       open(9,file='RDF3D.OPT')
98     read(9,900,err=99,end=99)so
       write(6,900)so
900    format(a5)
       if(so.eq.'ZATOM')read(9,*)ifa
       if(so.eq.'AMATR')read(9,*)lamat
       if(so.eq.'ANALC')read(9,*)lanal
       if(so.eq.'SMOOT')read(9,*)lsmooth
       if(so.eq.'PERIO')read(9,*)lper,ax,ay,az
       if(so.eq.'FILEN')read(9,9001)fn
       if(so.eq.'ATOMS')read(9,*)nats,(alist(i),i=1,nats)
       if(so.eq.'CUTMI')read(9,*)cuti
       if(so.eq.'CUTMA')read(9,*)cuta
       if(so.eq.'LSKIP')read(9,*)lskip
9001   format(a80)
       goto 98
99     close(9)
      endif        
      write(6,*)
      write(6,4003)' Transformation matrices:',lamat
      write(6,4003)'      Analyze deviations:',lanal
      write(6,4003)'               Smoothing:',lsmooth
      write(6,4003)'            Periodic Box:',lper
      write(6,4003)'             Skip solute:',lskip
4003  format(a25,l2)
      if(lper)write(6,6003)ax,ay,az
6003  format(      '              Dimensions:',3f15.5)
      call ifaifa(ifa,sn2,is2,ie2)

      do 101 i=1,nc
      rh(i)=0.0d0
      rm(i)=0.0d0
101   ro(i)=0.0d0
      jo=0
      jm=0
      jh=0

      inquire(file='STRUCT.LST',exist=lex2)
      if(lex2)then
       open(2,file='STRUCT.LST',status='old')
       write(6,*)'STRUCT.LST found'
      else
       inquire(file=fn,exist=lex3)
       if(lex3)then
        write(6,*)fn
        open(2,file=fn,status='old')
       else
        call report('No input geometries defined')
       endif
      endif

      if(lamat.and.lper)then
       inquire(file='AMAT.LST',exist=lex4)
       if(lex4)then
        open(29,file='AMAT.LST',status='old')
        write(6,*)'AMAT.LST found'
        read(29,*)
       else
        call report('AMAT.LST not found')
       endif
      endif

      if(lanal)then
       if(ifa.eq.0)then
        call rrf('mrdf.cub',nats,izq,ra,dr,r0,nd,ro0)
       else
        call rrf(sn2(is2:ie2)//'rdf.cub',nats,izq,ra,dr,r0,nd,ro0)
       endif
       write(6,*)'Old density loaded'
       open(31,file='PROD.PRN')
      endif

      nmc=0
      nwmo=0
      icvar=0
1     continue
      if(lex2)then
       read(2,200,end=500,err=500)fn
200    format(a80)
       open(4,file=fn)
       read(4,*)nat
      else
       read(2,200,end=500,err=500)fn
       read(2,*,end=500,err=500)nat
      endif

      if(lamat.and.lper)then
       read(29,*)
       read(29,291)A
291    format(3(3F10.4,/))
       do 9 i=1,3
       px(i)=A(i,1)*ax
       py(i)=A(i,2)*ay
9      pz(i)=A(i,3)*az
      else
       do 8 i=1,3
       px(i)=0.0d0
       py(i)=0.0d0
8      pz(i)=0.0d0
       px(1)=ax
       py(2)=ay
       pz(3)=az
      endif

      if(lanal)then
       do 1011 i=1,nc
1011   rm(i)=0.0d0
      endif
       
      if(nat.gt.nat0)call report('Too many atoms') 
      nmc=nmc+1
      if(lex2)then
       do 21 ia=1,nat
       read(4,400)sy,(r(3*(ia-1)+ix),ix=1,3)
400    format(7x,a3,1x,3f12.6)
       iq=0
       if(sy.eq.' H ')iq=1
       if(sy.eq.' N ')iq=7
       if(sy.eq.' O ')iq=8
       if(sy.eq.' C ')iq=6
       if(sy.eq.' F ')iq=9
       if(sy.eq.' Br')iq=35
       if(sy.eq.' Cl')iq=17
       if(iq.eq.0)then
        if(sy(2:2).eq.'H')iq=1
        if(sy(2:2).eq.'O')iq=8
        if(sy(2:2).eq.'C')iq=6
        if(sy(2:2).eq.'N')iq=7
       endif
21     if(iq.eq.0)write(6,*)sy//' unknown atom!'
       close(4)
      else
       do 22 ia=1,nat
22     read(2,*)iza(ia),(r(3*(ia-1)+ix),ix=1,3)
      endif
c

      if(nmc.eq.1)then
       do 3 ia=1,nats
       izq(ia)=iza(alist(ia))
       do 3 ix=1,3
3      ra(ix+3*(ia-1))=r(ix+3*(alist(ia)-1))
       write(6,*)nats,' atoms selected from the first geometry'
      endif
c
c     maximum square of O-H distance:
      doh=1.85d0
c     number of water molecules:
      nwm=0
c
c     if per. box, scan 3x3x3 superbox
      if(lper)then
       ir=27
       ix=-1
       iy=-1
       iz=-1
      else
       ir=1
       ix=0
       iy=0
       iz=0
      endif
c
      do 61 ic=1,ir
      ppx=px(1)*dble(ix)+py(1)*dble(iy)+pz(1)*dble(iz)
      ppy=px(2)*dble(ix)+py(2)*dble(iy)+pz(2)*dble(iz)
      ppz=px(3)*dble(ix)+py(3)*dble(iy)+pz(3)*dble(iz)

      do 6 ia=1,nat
      iq=iza(ia)
      x=r(3*(ia-1)+1)+ppx
      y=r(3*(ia-1)+2)+ppy
      z=r(3*(ia-1)+3)+ppz
c
c     skip solute
      do 63 iaa=1,nats
63    if(alist(iaa).eq.ia.and.lskip)goto 6

      if(ifa.ne.0.and.(ifa.eq.iq.or.ifa.lt.0))then
c      wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
c      specific atom type or global density
       nwm=nwm+1
       if(lsmooth)then
        call addxsmooth(nd,x,y,z,r0,dr,rm,jm,ifa,iq)
       else
        call addx(nd,x,y,z,r0,dr,rm,jm,ifa,iq)
       endif
      else
       if(ifa.eq.0)then
c       wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
c       water molecules
c       
c       decide, whether this is water O:
        if(iq.eq.8)then
c        ia is oxygen
         nh=0
         do 62 ja=1,nat
         if(iza(ja).eq.1)then
          xj=r(3*(ja-1)+1)+ppx
          yj=r(3*(ja-1)+2)+ppy
          zj=r(3*(ja-1)+3)+ppz
          d2=(x-xj)**2+(y-yj)**2+(z-zj)**2
          if(d2.lt.doh)nh=nh+1
          if(d2.lt.doh.and.nh.eq.1)nh1=ja
          if(d2.lt.doh.and.nh.eq.2)nh2=ja
         endif
62       continue
        
         if(nh.eq.2)then
c         ia is bond to two hydrogens nh1 and nh2
          nwm=nwm+1
c       
c         RDF - O
          xw=x
          yw=y
          zw=z
          if(lsmooth)then
           call addxsmooth(nd,xw,yw,zw,r0,dr,ro,jo,ifa,0)
          else
           call addx(nd,xw,yw,zw,r0,dr,ro,jo,ifa,0)
          endif
c       
c         RDF - H
          xw1=r(3*(nh1-1)+1)+ppx
          yw1=r(3*(nh1-1)+2)+ppy
          zw1=r(3*(nh1-1)+3)+ppz
          if(lsmooth)then
           call addxsmooth(nd,xw1,yw1,zw1,r0,dr,rh,jh,ifa,0)
          else
           call addx(nd,xw1,yw1,zw1,r0,dr,rh,jh,ifa,0)
          endif
          xw2=r(3*(nh2-1)+1)+ppx
          yw2=r(3*(nh2-1)+2)+ppy
          zw2=r(3*(nh2-1)+3)+ppz
          if(lsmooth)then
           call addxsmooth(nd,xw2,yw2,zw2,r0,dr,rh,jh,ifa,0)
          else
           call addx(nd,xw2,yw2,zw2,r0,dr,rh,jh,ifa,0)
          endif
c       
c         RDF - mass center
          xw=(xw1+xw2+16.0d0*x)/18.0d0
          yw=(yw1+yw2+16.0d0*y)/18.0d0
          zw=(zw1+zw2+16.0d0*z)/18.0d0
          if(lsmooth)then
           call addxsmooth(nd,xw,yw,zw,r0,dr,rm,jm,ifa,0)
          else
           call addx(nd,xw,yw,zw,r0,dr,rm,jm,ifa,0)
          endif
         endif
c        if(nh.eq.2)
         if(nh.ne.2.and.ia.gt.13)then
          write(6,*)fn
          write(6,*)nmc,ia,nh,ic
         endif
        endif
c       if(iq.eq.8)
c       wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
       endif
      endif
c     if(water RDF
6     continue

      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
      if(nwm.ne.nwmo.and.nwmo.ne.0)icvar=icvar+1
      nwmo=nwm
c       
c     correlation to recorded density:
      if(lanal)then
       call normd(dr(1)*dr(2)*dr(3),nc,rh,rm,ro,1,jh,jm,jo)
       diff=0.0d0
       da=0.0d0
       db=0.0d0
       do 1012 i=1,nc
       da=da+rm(i)**2
       db=db+ro0(i)**2
1012   diff=diff+rm(i)*ro0(i)
       if(da*db.ne.0.0d0)diff=diff/sqrt(da*db)
       if(da*db.ne.0.0d0)diff2=diff-sqrt(db/da)
       write(31,3131)nmc,diff,diff2
3131   format(i7,2G20.8)
      endif

      goto 1
500   close(2)
      if(lamat.and.lper)close(29)
c
      if(icvar.gt.0)write(6,*)'Warning: Variable number of H2O',
     1' molecules ',icvar,' times !'
      write(6,*)nmc,' geometries'
      if(ifa.eq.0)then
       write(6,*)nwm,' waters'
       write(6,*)jo,' os'
       write(6,*)jm,' ms'
       write(6,*)jh,' hs'
      else
       write(6,*)jm,' atoms',ifa
       write(6,*)nwm,' in one snapshot'
      endif

      if(lanal)then
       close(31)
       write(6,*)' PROD.PRN written.'
       stop
      endif
c
c     normalize to grid cell volume and one geometry
      call normd(dr(1)*dr(2)*dr(3),nc,rh,rm,ro,nmc,jh,jm,jo)

      if(cuti.ne.0.0d0.or.cuta.ne.0.0d0)then
       j=0
       write(6,*)'Removing cut-off parts'
       write(6,800)cuti,cuta
800    format(' CUTMIN CUTMAX:',2F7.2)
       do 5 ix=1,nd(1)
       x=r0(1)+dble(ix)*dr(1)
       do 5 iy=1,nd(2)
       y=r0(2)+dble(iy)*dr(2)
       do 5 iz=1,nd(3)
       z=r0(3)+dble(iz)*dr(3)
       ii=iz+(iy-1)*nd(3)+(ix-1)*nd(3)*nd(2)
       ica=0
       ici=0
       do 7 ia=1,nats
       i=1+3*(alist(ia)-1)
       d2=(x-r(i))**2+(y-r(i+1))**2+(z-r(i+2))**2
       if(d2.lt.cuta**2)ica=ica+1
7      if(d2.gt.cuti**2)ici=ici+1
       if(ica.eq.0.or.ica.eq.0)then
        rh(ii)=0.0d0
        ro(ii)=0.0d0
        rm(ii)=0.0d0
        j=j+1
       endif
5      continue
       write(6,*)j,' points removed'
      endif

      if(ifa.eq.0)then
       do 4 i=1,nc
4      rh(i)=-rh(i)
       call wrf('hrdf.cub',nats,izq,ra,dr,r0,nd,rh)
       call wrf('ordf.cub',nats,izq,ra,dr,r0,nd,ro)
       call wrf('mrdf.cub',nats,izq,ra,dr,r0,nd,rm)
      else
       call wrf(sn2(is2:ie2)//'rdf.cub',nats,izq,ra,dr,r0,nd,rm)
      endif

      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      write(6,*)'****** Program stopped. *********'
      stop
      end
      
      subroutine wrf(sn,nats,izq,ra,dr,r0,nd,rh)
      implicit none
      integer*4 nats,izq(*),nd(3),i,ix,iy,iz
      real*8 ra(*),rh(*),r0(3),dr(3),bohr
      character*(*)sn
c      
      bohr=0.529177d0
c
      open(3,file=sn)
      write(3,3000)nats,(r0(ix)/bohr,ix=1,3),nd(1),dr(1)/bohr,0.0,0.0,
     1nd(2),0.0,dr(2)/bohr,0.0,nd(3),0.0,0.0,dr(3)/bohr
3000  format(' Nuclear density',/,' SCF Total Density',/,
     13(i5,3f12.6,/),i5,3f12.6)
      do 2 i=1,nats
2     write(3,3001)izq(i),dble(izq(i)),(ra(ix+3*(i-1))/bohr,ix=1,3)
3001  format(i5,4f12.6)
      do 1 ix=1,nd(1)
      do 1 iy=1,nd(2)
1     write(3,4000)(rh(iz+(iy-1)*nd(3)+(ix-1)*nd(3)*nd(2)),iz=1,nd(3))
4000  format(6E13.5)
      close(3)
      write(6,*)sn
      return
      end

      subroutine rrf(sn,nats,izq,ra,dr,r0,nd,rh)
      implicit none
      integer*4 nats,izq(*),nd(3),i,ix,iy,iz
      real*8 ra(*),rh(*),r0(3),dr(3),bohr,adum
      character*(*)sn
c      
      bohr=0.529177d0
c
      open(3,file=sn)
      read(3,*)
      read(3,*)
      read(3,*)nats,(r0(ix),ix=1,3)
      read(3,*)nd(1),dr(1)
      read(3,*)nd(2),dr(2),dr(2)
      read(3,*)nd(3),dr(3),dr(3),dr(3)
      do 1 ix=1,3
      dr(ix)=dr(ix)*bohr
1     r0(ix)=r0(ix)*bohr
      do 2 i=1,nats
      read(3,*)izq(i),adum,(ra(ix+3*(i-1)),ix=1,3)
      do 2 ix=1,3
2     ra(ix+3*(i-1))=ra(ix+3*(i-1))*bohr
      do 3 ix=1,nd(1)
      do 3 iy=1,nd(2)
3     read(3,4000)(rh(iz+(iy-1)*nd(3)+(ix-1)*nd(3)*nd(2)),iz=1,nd(3))
4000  format(6E13.5)
      close(3)
      write(6,*)sn//' read in'
      return
      end

      subroutine addxsmooth(nd,xw,yw,zw,r0,dr,r,jo,ifa,iq)
      implicit none
      real*8 r(*),r0(3),dr(3),xw,yw,zw,vx,vy,vz,wx,wy,wz,
     1w(27),a14,dr1,dr2,dr3,an,p,weight
      integer*4 j,nd(3),ix,i,k,iy,iz,ni,jo,ifa,iq
      real*4 amas(89)
      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/     
c
      dr1=dr(1)
      dr2=dr(2)
      dr3=dr(3)
      ix=int((xw-r0(1))/dr1)+1
      if(ix.le.1.or.ix.gt.nd(1))return
      wx=xw-dr1*(0.50d0+dble(ix-1))-r0(1)
      iy=int((yw-r0(2))/dr2)+1
      if(iy.le.1.or.iy.gt.nd(2))return
      wy=yw-dr2*(0.50d0+dble(iy-1))-r0(2)
      iz=int((zw-r0(3))/dr3)+1
      if(iz.le.1.or.iz.gt.nd(3))return
      wz=zw-dr3*(0.50d0+dble(iz-1))-r0(3)
c
      weight=1.0d0
c     mass:
      if(ifa.eq.-1)weight=dble(amas(iq))
c     charge
      if(ifa.eq.-2)weight=dble(iq)

c     zero cube contributions
      do 2 i=1,27
2     w(i)=0.d0
      a14=0.0d0

c     loop over the all 27 cubes
      vx=-dr1-dr1
      do 1 i=-1,1
      vx=vx+dr1
      vy=-dr2-dr2
      do 1 j=-1,1
      vy=vy+dr2
      vz=-dr3-dr3
      do 1 k=-1,1
      vz=vz+dr3
      if(k.ne.0.or.j.ne.0.or.i.ne.0)then
       p=(wx*vx+wy*vy+wz*vz)/(vx*vx+vy*vy+vz*vz)
       if(dabs(p).gt.1.d0)then
        write(6,*)xw,yw,zw,ix,iy,iz,dr1,dr2,dr3,r0,vx,vy,vz
        write(6,*)p
        write(6,*)wx,wy,wz
        call report('error')
       endif
       if(p.gt.0.d0)then
        w(k+2+(j+1)*3+(i+1)*9)=p
        w(14)=w(14)+(1.0d0-p)
        a14=a14+1.0d0
       endif
      endif
1     continue
c
c     adjust central cube:
      w(14)=w(14)/a14
c
c     normalize total contribution to the weight:
      an=0.0d0
      do 3 i=1,27
3     an=an+w(i)
      if(an.gt.0.0d0)then
       do 4 i=1,27
4      w(i)=w(i)/an*weight
      endif
c
c     increment RDF, or difference for analysis:
      do 5 i=-1,1
      do 5 j=-1,1
      do 5 k=-1,1
      if(ix+i.gt.0.and.ix+i.le.nd(1).and.
     1   iy+j.gt.0.and.iy+j.le.nd(2).and.
     1   iz+k.gt.0.and.iz+k.le.nd(3))then
       ni=k+iz+(j+iy-1)*nd(3)+(i+ix-1)*nd(3)*nd(2)
       r(ni)=r(ni)+w(k+2+(j+1)*3+(i+1)*9)
      endif
5     continue
      jo=jo+1
      return
      end

      subroutine addx(nd,xw,yw,zw,r0,dr,r,j,ifa,iq)
      implicit none
      real*8 r(*),r0(3),dr(3),xw,yw,zw,weight
      integer*4 j,nd(3),ix,iy,iz,ni,ifa,iq
      real*4 amas(89)
      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/     
      
      weight=1.0d0
c     mass:
      if(ifa.eq.-1)weight=dble(amas(iq))
c     charge
      if(ifa.eq.-2)weight=dble(iq)

      ix=int((xw-r0(1))/dr(1))+1
      if(ix.le.0.or.ix.gt.nd(1))return
      iy=int((yw-r0(2))/dr(2))+1
      if(iy.le.0.or.iy.gt.nd(2))return
      iz=int((zw-r0(3))/dr(3))+1
      if(iz.le.0.or.iz.gt.nd(3))return
      ni=iz+(iy-1)*nd(3)+(ix-1)*nd(3)*nd(2)
      r(ni)=r(ni)+weight
      j=j+1
      return
      end

      subroutine ifaifa(ifa,sn2,is2,ie2)
      implicit none
      integer*4 ifa,ie2,is2
      character*80 sn2
      write(sn2,309)iabs(ifa)
309   format(i80)
      if(ifa.gt.0)then
       write(6,*)'Distribution of atom type ',ifa
      else
       if(ifa.eq.0)then
        write(6,*)' Water RDF'
       else
        if(ifa.eq.-1)then
         write(6,*)' Mass density'
         sn2='romass'
        else
         if(ifa.eq.-2)then
          write(6,*)' Charge density'
          sn2='rocharge'
         else
          if(ifa.eq.-3)then
           write(6,*)' Atomic density'
           sn2='roatom'
          else
           call report('Unknown ZATOM option.')
          endif
         endif
        endif
       endif
      endif
      do 111 is2=1,80
111   if(sn2(is2:is2).ne.' ')goto 121
121   do 131 ie2=80,1,-1
131   if(sn2(ie2:ie2).ne.' ')goto 141
141   return
      end

      subroutine normd(v,nc,rh,rm,ro,n,jh,jm,jo)
      implicit none
      real*8 v,rh(*),rm(*),ro(*),f
      integer*4 n,nc,jh,jm,jo,i
      f=1.0d0/v/dble(n)
      if(jh.gt.0)then
       do 1 i=1,nc
1      rh(i)=rh(i)*f
      endif
      if(jm.gt.0)then
       do 2 i=1,nc
2      rm(i)=rm(i)*f
      endif
      if(jo.gt.0)then
       do 3 i=1,nc
3      ro(i)=ro(i)*f
      endif
      return
      end