program rms3d
      implicit none
      integer*4 nat0,nres0
      parameter (nat0=5000,nres0=100)
      integer*4 ngeo,nat,i,is,ix,ntype,it1(nat0),it2(nat0),ia,n1,n2,it,
     1j,ja,jx
      real*8 r(3,nat0),rmean(3,nat0),rms,rs(nat0),cc(nres0,nres0),vi,
     1vj,r2(nres0)
c
      write(6,6004)
6004  format(/,'     RMS MD output deviations',/,/,
     1         '       all.x ... geometries',/,
     1         '       table.txt.. atom groups',/)
c
      ngeo=0
      open(9,file='all.x')
1     read(9,*,end=99,err=99)
      read(9,*,end=99,err=99)nat
      if(nat.gt.nat0)call report('too many atoms')
      ngeo=ngeo+1
      if(ngeo.eq.1)then
       write(6,*)nat,' atoms'
       do 2 i=1,nat
       do 2 ix=1,3
2      rmean(ix,i)=0.0d0
      endif
      do 3 i=1,nat
3     read(9,*)ntype,(r(ix,i),ix=1,3)
      do 4 i=1,nat
      do 4 ix=1,3
4     rmean(ix,i)=rmean(ix,i)+r(ix,i)
      goto 1
99    close(9)
      write(6,*)ngeo,' geometries'
      do 5 i=1,nat
      do 5 ix=1,3
5     rmean(ix,i)=rmean(ix,i)/dble(ngeo)
      write(6,*)'mean geometry calculated'

      open(9,file='table.txt')
      it=0
7     read(9,*,err=88,end=88)n1,n2
      it=it+1
      if(it.gt.nres0)call report('too many residues')
      it1(it)=n1
      it2(it)=n2
      goto 7
88    close(9)
      write(6,*)it,' residues'
      do 8 i=1,it
      rs(i)=0.0d0
      r2(i)=0.0d0
      do 8 j=1,it
8     cc(i,j)=0.0d0

      open(9,file='all.x')
      open(91,file='geo.txt')
      do 6 is=1,ngeo
      read(9,*)
      read(9,*)
      do 61 i=1,nat
61    read(9,*)ntype,(r(ix,i),ix=1,3)

      do 9 i=1,it
      vi=0.0d0
      do 92 ia=it1(i),it2(i)
      do 92 ix=1,3
92    vi=vi+(r(ix,ia)-rmean(ix,ia))**2
      rs(i)=rs(i)+vi
      r2(i)=r2(i)+vi**2
      do 9  j=1,it
      vj=0.0d0
      do 91 ja=it1(j),it2(j)
      do 91 jx=1,3
91    vj=vj+(r(jx,ja)-rmean(jx,ja))**2
9     cc(i,j)=cc(i,j)+vi*vj

      rms=0.0d0
      do 62 i=1,nat
      do 62 ix=1,3
62    rms=rms+(r(ix,i)-rmean(ix,i))**2
      rms=dsqrt(rms/dble(nat))/3.0d0
6     write(91,6000)is,rms
6000  format(i8,f15.3,' A')
      close(9)
      close(91)
      write(6,*)'geometry rms in geo.txt'

      open(9,file='cc.txt')
      do 93 i=1,it
      do 93 j=1,it
      cc(i,j)=cc(i,j)/dsqrt(r2(i)*r2(j))
93    write(9,9000)i,j,cc(i,j)
9000  format(2i6,f12.6)
      close(9)
      write(6,*)'correlation in cc.txt'

      open(9,file='res.txt')
      do 10 i=1,it
      rs(i)=dsqrt(rs(i)/dble(ngeo)/dble(it2(i)-it1(i)+1))/3.0d0
10    write(9,6000)i,rs(i)
      close(9)
      write(6,*)'residue rms in res.txt'

      stop
      end

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