diff --git a/src/analyz/GNUmakefile b/src/analyz/GNUmakefile
index 20ad45c33b..e02b6deb19 100644
--- a/src/analyz/GNUmakefile
+++ b/src/analyz/GNUmakefile
@@ -1,5 +1,5 @@
 #
-# $Id: GNUmakefile,v 1.20 2003-12-03 05:04:04 d3j191 Exp $
+# $Id: GNUmakefile,v 1.21 2004-03-06 03:41:15 d3j191 Exp $
 #
  OBJ_OPTIMIZE = ana_angle.o\
                 ana_atrad.o\
@@ -19,6 +19,7 @@
                 ana_pltgrd.o\
                 ana_rama.o\
                 ana_rdhdr.o\
+                ana_rdf.o\
                 ana_rdfram.o\
                 ana_rdref.o\
                 ana_rdtop.o\
diff --git a/src/analyz/ana_common.fh b/src/analyz/ana_common.fh
index 92458be89b..cfe7c9ee1a 100644
--- a/src/analyz/ana_common.fh
+++ b/src/analyz/ana_common.fh
@@ -1,5 +1,5 @@
 c
-c $Id: ana_common.fh,v 1.60 2003-12-03 05:04:04 d3j191 Exp $
+c $Id: ana_common.fh,v 1.61 2004-03-06 03:41:15 d3j191 Exp $
 c
       integer mxbnds,mxangs,mxtors,mximps
       parameter(mxbnds=100)
@@ -24,11 +24,11 @@ c
       integer irtdb,me,np,lfnana,lfnref,lfntrj,lfncop,lfnsup,
      + lfnchg,lfnplt,lfnxyz,lfnrms,lfnprj,lfnval,lfnmin,lfnmax,lfncov,
      + lfnvec,lfntop,lfnhis,lfnepz,lfngrp,lfnram,lfnord,lfnloc,lfnsel,
-     + lfnhba,lfnecc,lfndie
+     + lfnhba,lfnecc,lfndie,lfnrdf
       common/ana_rtdb/irtdb,me,np,lfnana,lfnref,lfntrj,lfncop,lfnsup,
      + lfnchg,lfnplt,lfnxyz,lfnrms,lfnprj,lfnval,lfnmin,lfnmax,lfncov,
      + lfnvec,lfntop,lfnhis,lfnepz,lfngrp,lfnram,lfnord,lfnloc,lfnsel,
-     + lfnhba,lfnecc,lfndie
+     + lfnhba,lfnecc,lfndie,lfnrdf
 c
       real*8 time,timr,timoff,box(3),spac(3),rcut,xmax(3),xmin(3)
       real*8 temp,pres,xsmin,xsmax,scale,cpk,stick,rangle
@@ -70,7 +70,8 @@ c
      + i_xp,l_xp,i_ord,l_ord,i_hist,l_hist,i_qdat,l_qdat,i_wsel,l_wsel,
      + i_ndxw,l_ndxw,i_val,l_val,i_iram,l_iram,i_imol,l_imol,
      + i_wrk,l_wrk,i_ssel,l_ssel,i_swt,l_swt,i_hbnd,l_hbnd,
-     + i_sbnd,l_sbnd,i_osel,l_osel,i_owt,l_owt,i_qwdat,l_qwdat
+     + i_sbnd,l_sbnd,i_osel,l_osel,i_owt,l_owt,i_qwdat,l_qwdat,
+     + i_rdf,l_rdf
       common/ana_ptr/i_xref,l_xref,i_xrms,l_xrms,i_wt,l_wt,i_isel,
      + l_isel,i_snam,l_snam,i_xdat,l_xdat,i_idat,l_idat,i_bnd,l_bnd,
      + i_x,l_x,i_q,l_q,i_t,l_t,i_grid,l_grid,i_tag,l_tag,ga_vec,
@@ -78,7 +79,8 @@ c
      + i_xp,l_xp,i_ord,l_ord,i_hist,l_hist,i_qdat,l_qdat,i_wsel,l_wsel,
      + i_ndxw,l_ndxw,i_val,l_val,i_iram,l_iram,i_imol,l_imol,
      + i_wrk,l_wrk,i_ssel,l_ssel,i_swt,l_swt,i_hbnd,l_hbnd,
-     + i_sbnd,l_sbnd,i_osel,l_osel,i_owt,l_owt,i_qwdat,l_qwdat
+     + i_sbnd,l_sbnd,i_osel,l_osel,i_owt,l_owt,i_qwdat,l_qwdat,
+     + i_rdf,l_rdf
 c
       integer nsa,msa,nwm,mwm,nwa,mwa,ndata,nq,nave,nsel,npov,nwsel,
      + nwrit,indx,msgm,nsgm,msb,mwb,nselo,ifr,ito
@@ -86,14 +88,15 @@ c
      + nwsel,nwrit,indx,msgm,nsgm,msb,mwb,nselo,ifr,ito
 c
       logical lxw,lvw,lfw,lpw,lsx,lsv,lfs,lps,lsel,lana,lrms,lsonly,
-     + ltop,ldcd,lrama,lesd,lloc,lsels,lhbond,lselo,ldiel
+     + ltop,ldcd,lrama,lesd,lloc,lsels,lhbond,lselo,ldiel,lrdf
       common/ana_log/lxw,lvw,lfw,lpw,lsx,lsv,lfs,lps,lsel,lana,lrms,
-     + lsonly,ltop,ldcd,lrama,lesd,lloc,lsels,lhbond,lselo,ldiel
+     + lsonly,ltop,ldcd,lrama,lesd,lloc,lsels,lhbond,lselo,ldiel,lrdf
 c
       character*24 wtag(100,2)
       character*16 wnam(100)
       character*10 datum,tijd
-      common/ana_chr/wtag,wnam,datum,tijd
+      character*3 cnum
+      common/ana_chr/wtag,wnam,datum,tijd,cnum
 c
       real*8 wval(100,2)
       common/ana_val/wval
@@ -103,13 +106,13 @@ c
      + ihbndw(100),ibndw(100,2)
       real*8 rgroups(maxgrp,2),rgroup(maxgrp,2)
       integer nhis,lhis,ihis(mlhis,mhis),idhis(mhis,3)
-      integer nord,idord(mord,4),iord(mord,2)
+      integer nord,idord(mord,4),iord(mord,2),nrdf,numrdf
       real*8 rord(mord,2)
-      real*8 rhis,dhis
-      common/ana_def/rord,rgroups,rgroup,rhis,dhis,
+      real*8 rhis,dhis,rrdf
+      common/ana_def/rord,rgroups,rgroup,rhis,dhis,rrdf,
      + ndef,ldef,idef,ngroups,igroups,ngroup,igroup,
      + ihbndw,ibndw,nhis,lhis,ihis,idhis,
-     + nord,idord,iord
+     + nord,idord,iord,nrdf,numrdf
 c
       integer nrot,irot(100)
       real*8 arot(100)
diff --git a/src/analyz/ana_input.F b/src/analyz/ana_input.F
index ee3fc3859e..13871c503b 100644
--- a/src/analyz/ana_input.F
+++ b/src/analyz/ana_input.F
@@ -1,6 +1,6 @@
       subroutine ana_input(irtdb)
 c
-c $Id: ana_input.F,v 1.67 2003-12-03 05:04:04 d3j191 Exp $
+c $Id: ana_input.F,v 1.68 2004-03-06 03:41:15 d3j191 Exp $
 c
       implicit none
 c
@@ -16,9 +16,9 @@ c
       integer lfncmd,numcmd,len,ivec,nspac,ibond,iangl,itors,iplan,icd
       integer ifrfr,ifrto,ifrsk,ifrst,ilast,isel,jsel,ksel,lsel
       integer iopt,ivctr,iorder
-      integer itag,iatag,jatag,ipbc
+      integer itag,iatag,jatag,ipbc,nrdf
       real*8 timoff,rsel,rtag,rcut,scale,cpk,stick,rval,rmin,rmax,rang
-      real*8 arota
+      real*8 arota,rrdf
       logical lref,lfil,lsol,lsuper,lselo
       integer iesppb,lent,idpdb,indx,icent,jcent,nwhb
       character*8 option
@@ -35,6 +35,8 @@ c
       numcmd=0
       lref=.false.
       lfil=.false.
+      nrdf=0
+      rrdf=1.0d0
 c
       timoff=0.0d0
       scale=1.0d0
@@ -450,6 +452,19 @@ c
       goto 2
       endif
 c
+c     radial distribution function
+c     ----------------------------
+c
+      if(inp_compare(.false.,'rdf',item)) then
+      if(.not.inp_i(nrdf))
+     + call md_abort('No rdf number of bins specified',0)
+      if(.not.inp_f(rrdf)) call md_abort('No rdf range specified',0)
+      write(lfncmd,1032) nrdf,rrdf
+ 1032 format('rdf    ',i7,f12.6)
+      numcmd=numcmd+1
+      goto 2
+      endif
+c
 c     group analysis
 c     --------------
 c
diff --git a/src/analyz/ana_rdf.F b/src/analyz/ana_rdf.F
new file mode 100644
index 0000000000..d7b7243a85
--- /dev/null
+++ b/src/analyz/ana_rdf.F
@@ -0,0 +1,137 @@
+      subroutine ana_rdf_init()
+c
+c $Id: ana_rdf.F,v 1.1 2004-03-06 03:41:15 d3j191 Exp $
+c
+      implicit none
+c
+#include "ana_common.fh"
+#include "global.fh"
+#include "mafdecls.fh"
+#include "msgids.fh"
+#include "rtdb.fh"
+c
+      if(.not.ma_push_get(mt_int,nrdf*nsel*mwa,'irdf',l_rdf,i_rdf))
+     + call md_abort('Could not allocate irdf',0)
+      print*,'rdf allocated in rdfhdr ',nrdf*nsel*mwa
+c
+      return
+      end
+      subroutine ana_rdfhdr(irdf)
+c
+      implicit none
+c
+#include "ana_common.fh"
+#include "global.fh"
+#include "mafdecls.fh"
+#include "msgids.fh"
+#include "rtdb.fh"
+c
+      integer irdf(nsel,mwa,nrdf)
+c
+      character*255 fname
+      integer i,j,k,lq,l,m
+c
+      fname=filtrj
+      lq=index(filtrj,'?')
+      if(lq.gt.0) then
+      fname=filtrj(1:lq-1)//cnum//filtrj(lq+1:index(filtrj,' ')-1)
+      endif
+      lq=index(fname,'.trj')
+      fname(lq:lq+3)='.rdf'
+c
+      open(unit=lfnrdf,file=fname(1:index(fname,' ')-1),
+     + status='unknown')
+c
+      write(*,3333) fname(1:index(fname,' ')-1)
+ 3333 format(/,' Opening rdf file ',a)
+      rewind(lfnrdf)
+c
+      numrdf=0
+c
+      do 1 i=1,nrdf
+      do 2 j=1,mwa
+      do 3 k=1,nsel
+      irdf(k,j,i)=0
+    3 continue
+    2 continue
+    1 continue
+c
+      print*,'rdf set to zero in rdfhdr'
+c
+      return
+      end
+      subroutine ana_rdf(isel,xs,xw,irdf)
+c
+      implicit none
+c
+#include "ana_common.fh"
+#include "global.fh"
+#include "mafdecls.fh"
+#include "msgids.fh"
+#include "rtdb.fh"
+c
+      integer isel(msa),irdf(nsel,mwa,nrdf)
+      real*8 xs(msa,3),xw(mwm,mwa,3)
+      integer i,j,k,l,m
+      real*8 d
+c
+      i=0
+      do 1 j=1,nsa
+      if(isel(j).gt.0) then
+      i=i+1
+      do 2 k=1,nwm
+      do 3 l=1,nwa
+      d=sqrt((xs(j,1)-xw(k,l,1))**2+(xs(j,2)-xw(k,l,2))**2+
+     + (xs(j,3)-xw(k,l,3))**2)
+      m=int(dble(nrdf*rrdf)/d)
+      if(m.le.nrdf) irdf(i,l,m)=irdf(i,l,m)+1
+    3 continue
+    2 continue
+      endif
+    1 continue
+c
+      numrdf=numrdf+1
+c
+      return
+      end
+      subroutine ana_rdfwrt(isel,irdf)
+c
+      implicit none
+c
+#include "ana_common.fh"
+#include "global.fh"
+#include "mafdecls.fh"
+#include "msgids.fh"
+#include "rtdb.fh"
+c
+      integer isel(msa)
+      integer irdf(nsel,mwa,nrdf)
+c
+      integer i,j,k,l,m
+      real*8 d
+c
+      write(*,'(a)') 'ANA_RDFWRT'
+c
+      do 4 i=1,nsel
+      do 5 m=1,nrdf
+      write(*,'(2i5,10i10)') i,m,(irdf(i,l,m),l=1,nwa)
+    5 continue
+    4 continue
+c
+c      i=0
+c      do 1 j=1,nsa
+c      if(isel(j).gt.0) then
+c      i=i+1
+c      do 2 k=1,nrdf
+c      d=dble(k)*rrdf/dble(nrdf)
+c      write(*,1000) d,(irdf(i,l,k),l=1,nwa)
+c 1000 format(f12.6,10i5)
+c    2 continue
+c      endif
+c    1 continue
+c
+      close(unit=lfnrdf,status='keep')
+      write(*,'(a)') ' Closing rdf file '
+c
+      return
+      end
diff --git a/src/analyz/ana_rdfram.F b/src/analyz/ana_rdfram.F
index 1e0d177383..2a297464e0 100644
--- a/src/analyz/ana_rdfram.F
+++ b/src/analyz/ana_rdfram.F
@@ -1,6 +1,6 @@
       logical function ana_rdfram(x,ix,w)
 c
-c $Id: ana_rdfram.F,v 1.32 2003-12-03 05:04:04 d3j191 Exp $
+c $Id: ana_rdfram.F,v 1.33 2004-03-06 03:41:15 d3j191 Exp $
 c
       implicit none
 c
@@ -15,7 +15,6 @@ c
       character*108 card
       integer i,j,k,lq
       character*255 fname
-      character*3 cnum
 c
       logical succes
       real*8 timp
@@ -141,6 +140,7 @@ c
 c
       close(unit=lfntrj)
       write(*,'(a)') ' Closing trj file '
+      if(lrdf) call ana_rdfwrt()
 c
       fname=filtrj
 c
@@ -156,6 +156,8 @@ c
      + status='old',err=9998)
       write(*,3333) fname(1:index(fname,' ')-1)
  3333 format(/,' Opening trj file ',a)
+c
+      if(lrdf) call ana_rdfhdr(int_mb(i_rdf))
       goto 100
 c
  9998 continue
diff --git a/src/analyz/ana_rdhdr.F b/src/analyz/ana_rdhdr.F
index 9478095479..2131a40303 100644
--- a/src/analyz/ana_rdhdr.F
+++ b/src/analyz/ana_rdhdr.F
@@ -1,16 +1,16 @@
       subroutine ana_rdhdr(sgmnam)
 c
-c $Id: ana_rdhdr.F,v 1.15 2003-10-19 03:31:01 d3j191 Exp $
+c $Id: ana_rdhdr.F,v 1.16 2004-03-06 03:41:15 d3j191 Exp $
 c
       implicit none
 c
 #include "ana_common.fh"
+#include "mafdecls.fh"
 c
       character*16 sgmnam(msa)
       character*80 card
       character*255 fname
       integer i,lq
-      character*3 cnum
 c
       integer ib,jb,nsb,nwb
 c
@@ -34,6 +34,7 @@ c
       write(*,3333) fname(1:index(fname,' ')-1)
  3333 format(/,' Opening trj file ',a)
       rewind(lfntrj)
+      if(lrdf) call ana_rdfhdr(int_mb(i_rdf))
     1 continue
       read(lfntrj,1000,err=9998,end=9997) card
  1000 format(a)
diff --git a/src/analyz/ana_rdtop.F b/src/analyz/ana_rdtop.F
index 2de598ed21..61db2c68d8 100644
--- a/src/analyz/ana_rdtop.F
+++ b/src/analyz/ana_rdtop.F
@@ -5,6 +5,7 @@ c
 #include "mafdecls.fh"
 #include "global.fh"
 #include "msgids.fh"
+#include "util.fh"
 #include "ana_common.fh"
 c
       character*16 sgmnam(msa)
@@ -12,14 +13,14 @@ c
       integer iram(msgm,7),imol(msa),ibnd(msb,2)
 c
       character*1 cdummy
-      integer i,j,nat,naq,i_tmp,l_tmp,i_itmp,l_itmp
+      integer i,j,k,nat,naq,nparm,nseq,i_tmp,l_tmp,i_itmp,l_itmp
       integer naw,nbw,nhw,ndw,now,ntw,nnw
       integer nas,nbs,num
       character*5 sname,aname
 c
       ltop=.false.
 c
-      do 9 i=1,nsgm
+      do 1 i=1,nsgm
       iram(i,1)=0
       iram(i,2)=0
       iram(i,3)=0
@@ -27,37 +28,47 @@ c
       iram(i,5)=0
       iram(i,6)=0
       iram(i,7)=0
-    9 continue
+    1 continue
 c
       if(me.eq.0) then
 c
       open(unit=lfntop,file=filtop(1:index(filtop,' ')-1),
      + form='formatted',status='old',err=999)
 c
+c      print*,'TOPOLOGY FILE ',filtop(1:index(filtop,' ')-1)
       read(lfntop,1000,end=999,err=999) cdummy
       read(lfntop,1000) cdummy
       read(lfntop,1000) cdummy
       read(lfntop,1000) cdummy
-      read(lfntop,1000) cdummy
  1000 format(a1)
+      read(lfntop,1001) nparm
       read(lfntop,1001) nat
       read(lfntop,1001) naq
+      read(lfntop,1001) nseq
+      read(lfntop,1000) cdummy
  1001 format(i5)
-      do 1 i=1,nat
+      do 2 i=1,nat*nparm
       read(lfntop,1000) cdummy
-    1 continue
-      do 2 i=1,nat
-      do 3 j=i,nat
-      read(lfntop,1000) cdummy
-      read(lfntop,1000) cdummy
-    3 continue
     2 continue
+      do 3 i=1,nat
+      do 4 j=i,nat
+      do 5 k=1,nparm
+      read(lfntop,1000) cdummy
+    5 continue
+    4 continue
+    3 continue
       if(.not.ma_push_get(mt_dbl,naq,'tmp',l_tmp,i_tmp))
      + call md_abort('Failed to allocate tmp',0)
-      do 4 i=1,naq
+      do 6 i=1,naq
       read(lfntop,1002) dbl_mb(i_tmp-1+i)
- 1002 format(f12.6)
-    4 continue
+ 1002 format(5x,f12.6)
+      do 7 j=1,nparm-1
+      read(lfntop,1000) cdummy
+    7 continue
+    6 continue
+      do 8 i=1,nseq
+      read(lfntop,1000) cdummy
+    8 continue
       read(lfntop,1003) naw,nbw,nhw,ndw,now,ntw,nnw
  1003 format(5i7,2i10)
       read(lfntop,1003) nas,nbs
@@ -66,20 +77,22 @@ c
       rewind(44)
       write(44,2003) naw,nas,nbs
  2003 format(3i10)
-      do 7 i=1,naw
+      do 9 i=1,naw
       read(lfntop,1005) sname,aname,num,num,j
       write(44,2004) sname,aname
  2004 format(2a5,i6)
       qw(i)=dbl_mb(i_tmp-1+j)
-    7 continue
-      do 19 i=1,nbw
+    9 continue
+      do 10 i=1,nbw
       read(lfntop,'(2i7)') ibndw(i,1),ibndw(i,2)
+      do 11 j=1,nparm
       read(lfntop,1000) cdummy
-   19 continue
-      nat=2*(nhw+ndw+now)
-      do 5 i=1,nat
+   11 continue
+   10 continue
+      nat=4*(nhw+ndw+now)
+      do 12 i=1,nat
       read(lfntop,1000) cdummy
-    5 continue
+   12 continue
       if(ntw.gt.0) then
       read(lfntop,1004) (j,i=1,ntw)
       read(lfntop,1004) (j,i=1,ntw)
@@ -90,9 +103,12 @@ c
       read(lfntop,1004) (j,i=1,nnw)
       endif
       read(lfntop,1000) cdummy
+      do 13 i=1,nparm
+      read(lfntop,1000) cdummy
+   13 continue
       if(.not.ma_push_get(mt_int,nas,'itmp',l_itmp,i_itmp))
      + call md_abort('Failed to allocate itmp',0)
-      do 6 i=1,nas
+      do 14 i=1,nas
       read(lfntop,1005) sname,aname,imol(i),num,j
  1005 format(a5,5x,a5,6x,2i5,15x,i5)
       write(44,2004) sname,aname,num
@@ -104,12 +120,14 @@ c
       if(aname.eq.' C   ') iram(num,4)=i
       if(aname.eq.' H   ') iram(num,6)=i
       if(aname.eq.' O   ') iram(num,7)=i
-    6 continue
-      do 8 i=1,nbs
+   14 continue
+      do 15 i=1,nbs
       read(lfntop,1006) j,num
       ibnd(i,1)=j
       ibnd(i,2)=num
+      do 16 k=1,nparm
       read(lfntop,1000) cdummy
+   16 continue
  1006 format(2i7)
       write(44,2005) j,num
  2005 format(2i8)
@@ -125,7 +143,7 @@ c
       iram(int_mb(i_itmp-1+num),5)=j
       endif
       endif
-    8 continue
+   15 continue
       if(.not.ma_pop_stack(l_itmp))
      + call md_abort('Failed to deallocate itmp',0)
       close(unit=44)
@@ -146,6 +164,7 @@ c
       endif
       endif
 c
+c      print*,'ANA_RDTOP DONE'
       return
       end
       subroutine ana_siztop()
@@ -160,10 +179,11 @@ c
       character*1 cdummy
       real*8 rdummy
 c
-      integer i,j,nat,naq
+      integer i,j,k,nat,naq,nparm,nseq
       integer naw,nbw,nhw,ndw,now,ntw,nnw
       integer nas,nbs,num
       character*5 sname,aname
+      character*80 string
 c
       ltop=.false.
       nsgm=0
@@ -177,24 +197,31 @@ c
       read(lfntop,1000) cdummy
       read(lfntop,1000) cdummy
       read(lfntop,1000) cdummy
-      read(lfntop,1000) cdummy
  1000 format(a1)
+ 5000 format(a)
+      read(lfntop,1001) nparm
       read(lfntop,1001) nat
       read(lfntop,1001) naq
+      read(lfntop,1001) nseq
  1001 format(i5)
-      do 1 i=1,nat
+      read(lfntop,1000) cdummy
+      do 1 i=1,nat*nparm
       read(lfntop,1000) cdummy
     1 continue
       do 2 i=1,nat
       do 3 j=i,nat
+      do 4 k=1,nparm
       read(lfntop,1000) cdummy
-      read(lfntop,1000) cdummy
+    4 continue
     3 continue
     2 continue
-      do 4 i=1,naq
+      do 5 i=1,naq*nparm
       read(lfntop,1002) rdummy
- 1002 format(f12.6)
-    4 continue
+ 1002 format(5x,f12.6)
+    5 continue
+      do 6 i=1,nseq
+      read(lfntop,1000) cdummy
+    6 continue
       read(lfntop,1003) naw,nbw,nhw,ndw,now,ntw,nnw
  1003 format(5i7,2i10)
       mwb=nbw
@@ -204,10 +231,10 @@ c
       read(lfntop,1005) sname,aname,num,j
  2004 format(2a5,i6)
     7 continue
-      nat=2*(nbw+nhw+ndw+now)
-      do 5 i=1,nat
+      nat=4*(nbw+nhw+ndw+now)
+      do 8 i=1,nat
       read(lfntop,1000) cdummy
-    5 continue
+    8 continue
       if(ntw.gt.0) then
       read(lfntop,1004) (j,i=1,ntw)
       read(lfntop,1004) (j,i=1,ntw)
@@ -218,17 +245,20 @@ c
       read(lfntop,1004) (j,i=1,nnw)
       endif
       read(lfntop,1000) cdummy
-      do 6 i=1,nas
+      do 9 i=1,nparm
+      read(lfntop,1000) cdummy
+    9 continue
+      do 10  i=1,nas
       read(lfntop,1005) sname,aname,num,j
  1005 format(a5,5x,a5,11x,i5,15x,i5)
       nsgm=max(nsgm,num)
-    6 continue
-      do 8 i=1,nbs
+   10 continue
+      do 11 i=1,nbs
       read(lfntop,1006) j,num
       read(lfntop,1000) cdummy
  1006 format(2i7)
  2005 format(2i8)
-    8 continue
+   11 continue
       close(unit=lfntop)
 c
       msgm=nsgm
diff --git a/src/analyz/ana_task.F b/src/analyz/ana_task.F
index f02348fd15..0013142e9e 100644
--- a/src/analyz/ana_task.F
+++ b/src/analyz/ana_task.F
@@ -1,6 +1,6 @@
       subroutine ana_task
 c
-c $Id: ana_task.F,v 1.94 2003-12-03 05:04:05 d3j191 Exp $
+c $Id: ana_task.F,v 1.95 2004-03-06 03:41:15 d3j191 Exp $
 c
       implicit none
 c
@@ -41,6 +41,7 @@ c
       lfnloc=70
       lfnhba=71
       lfnecc=72
+      lfnrdf=73
       lfnref=76
       lfntrj=77
       lfncop=78
@@ -105,6 +106,7 @@ c
       lesd=.false.
       lloc=.false.
       lhbond=.false.
+      lrdf=.false.
 c
       nave=0
       ndata=0
@@ -481,6 +483,15 @@ c
       goto 1
       endif
 c
+c     radial distribution function
+c     ----------------------------
+c
+      if(card(1:6).eq.'rdf   ') then
+      read(card,'(7x,i7,f12.6)') nrdf,rrdf
+      lrdf=.true.
+      goto 1
+      endif
+c
 c     order parameters
 c     ----------------
 c
@@ -861,6 +872,11 @@ c
       filhba=string(1:index(string,' ')-1)//'.hba'
       fildie=string(1:index(string,' ')-1)//'.die'
       time=0.0d0
+c
+      if(lrdf) then
+      call ana_rdf_init()
+      print*,'rdf_init'
+      endif
 c
 c     read file header
 c
@@ -1017,6 +1033,11 @@ c
      + dbl_mb(i_wt),dbl_mb(i_xdat),i,int_mb(i_wrk))
  3307 continue
       endif
+c
+      if(lrdf) then
+      call ana_rdf(int_mb(i_isel),dbl_mb(i_xdat),dbl_mb(i_wdat),
+     + int_mb(i_rdf))
+      endif
 c
       if(ngroups.gt.0) then
       if(me.eq.0) write(lfngrp,4004) time
@@ -1056,15 +1077,11 @@ c
      + call md_abort('Failed to deallocate hist',0)
       close(unit=lfnhis)
       endif
-c
-      if(lhbond) then
-      if(.not.ma_pop_stack(l_temp))
-     + call md_abort('Failed to deallocate temp',0)
-      endif
 c
       if(me.eq.0) then
       close(unit=lfntrj)
       write(*,'(a)') ' Closing trj file '
+      if(lrdf) call ana_rdfwrt(int_mb(i_isel),int_mb(i_rdf))
       if(lana) then
       close(unit=lfnana)
       write(*,'(a)') ' Closing scn file '
@@ -1105,6 +1122,17 @@ c
  3304 continue
       endif
       endif
+c
+      if(lrdf) then
+      print*,'deallocating rdf'
+      if(.not.ma_pop_stack(l_rdf))
+     + call md_abort('Failed to deallocate rdf',0)
+      endif
+c
+      if(lhbond) then
+      if(.not.ma_pop_stack(l_temp))
+     + call md_abort('Failed to deallocate temp',0)
+      endif
 c
       nbnds=0
       nangs=0
@@ -1457,6 +1485,7 @@ c
       endif
 c
       if(me.eq.0) close(unit=lfncmd,status='delete')
+c      if(me.eq.0) close(unit=lfncmd,status='keep')
 c
       if(active) call ana_edfinal()
 c
diff --git a/src/analyz/ana_wthdr.F b/src/analyz/ana_wthdr.F
index c76e2ab9cc..256c75b132 100644
--- a/src/analyz/ana_wthdr.F
+++ b/src/analyz/ana_wthdr.F
@@ -1,6 +1,6 @@
       subroutine ana_wthdr(iunit,fmt,sgmnam,tag,isel,logw)
 c
-c $Id: ana_wthdr.F,v 1.18 2003-10-19 03:31:02 d3j191 Exp $
+c $Id: ana_wthdr.F,v 1.19 2004-03-06 03:41:15 d3j191 Exp $
 c
       implicit none
 c
@@ -10,7 +10,7 @@ c
       integer isel(msa)
       integer iunit
       character*3 fmt
-      character*4 cnum
+      character*4 cdnum
       character*255 fname
       character*24 tag(msa,2)
       logical logw
@@ -42,8 +42,8 @@ c
       fname=filcop
       if(mcopf.gt.0) then
       lq=index(filcop,'.')
-      write(cnum,'(a1,i3.3)') '_',icopf
-      fname=filcop(1:lq-1)//cnum//filcop(lq:index(filcop,' ')-1)
+      write(cdnum,'(a1,i3.3)') '_',icopf
+      fname=filcop(1:lq-1)//cdnum//filcop(lq:index(filcop,' ')-1)
       endif
       if(binary) then
       open(unit=lfncop,file=fname(1:index(fname,' ')-1),
@@ -60,8 +60,8 @@ c
       fname=filsup
       if(msupf.gt.0) then
       lq=index(filsup,'.')
-      write(cnum,'(a1,i3.3)') '_',isupf
-      fname=filsup(1:lq-1)//cnum//filsup(lq:index(filsup,' ')-1)
+      write(cdnum,'(a1,i3.3)') '_',isupf
+      fname=filsup(1:lq-1)//cdnum//filsup(lq:index(filsup,' ')-1)
       endif
       if(binary) then
       open(unit=lfnsup,file=fname(1:index(fname,' ')-1),