anchor.diff:         allocate(c2(n,3))
          allocate(gr(nb))
          allocate(gr0(nb))
  -       write(*,*) "number of atoms", n
  -c       call xyz_read_coords_byname(n,c,atom1_tag,fn_in)
  -c      -----------------------------
  -c      start looping over the frames
  -c      -----------------------------
  +       allocate(atag(n))
  +       allocate(oc1(n))
  +       allocate(oc2(n))
  +c      loop over frames
          nf = 0
          gr = 0
  -30     continue
          do
  -       n1=n
  -       n2=n
  -       call xyz_read_coords_byname2(n1,c1,atom1_tag,
  -     +                              n2,c2,atom2_tag,
  -     +                              fn_in)
  -       if(n1.eq.0.or.n2.eq.0) then
  -          write(*,*) "exiting because there are no coordinates",nf
  -          exit
  -       end if
  -       do i=1,n1
  -         write(34,*) (c1(i,k),k=1,3)
  -       end do
  -       do i=1,n2
  -         write(35,*) (c2(i,k),k=1,3)
  -       end do
  -
  -       gr0 = 0
  -       call rdf_compute(n, n1,c1,atom1_tag,
  -     +                  n2,c2,atom2_tag,
  -     +                  nb,rmax,lat,
  -     +                  gr0,
  -     +                  fn_in)
  -       gr = gr + gr0
  -       nf = nf+1
  -       if(nfrm.gt.0.and.nf.gt.nfrm) then
  -          write(*,*) "exiting because exceeded number of frames",nf,nfrm
  -          exit
  -       end if
  -       end do
  -c       goto 30
  +         call xyz_read(n,c,atag,fn_in)
  +         if(n.eq.0) exit
  +         write(*,*) "number of atoms", n
  +         oc1=.false.
  +         if(atom1_id.gt.0) then
  +           oc1=.false.
  +           oc1(atom1_id) = .true.
  +         else if(atom1_tag.ne." ") then
  +           oc1 = atag.eq.atom1_tag
  +c           call mask_byname(n,c,atag,oc1,atom1_tag)
  +         end if
  +         oc2=.false.
  +         if(atom2_id.gt.0) then
  +           oc2(atom2_id) = .true.
  +         else if(atom2_tag.ne." ") then
  +           oc2 = (atag.eq.atom2_tag).and.(.not.oc1)
  +c           call mask_byname(n,c,atag,oc2,atom2_tag)
  +         else
  +           oc2=oc1
  +         end if
  +         gr0 = 0
  +         call rdf_compute2(n,c,oc1,oc2,
  +     +                    nb,rmax,lat,
  +     +                    gr0,
  +     +                    fn_in)
  +         gr = gr + gr0
  +         nf = nf+1
  +         if(nfrm.gt.0.and.nf.gt.nfrm) then
  +            write(*,*) "exiting because exceeded number of frames",nf,nfrm
  +            exit
  +         end if
  +       end do
   c      ---------------------------
   c      end looping over the frames
   c      ---------------------------
  
positive.code:           call get_command_argument(i,buffer,l,istatus)
            if(istatus.ne.0) goto 18
            if(is_number(buffer)) then
              message = "Only atom tags are supported now"
              goto 911
              read(buffer,*) atom2_id
            else
              atom2_tag = buffer
  
negative.code:        grs = 0.0d0
  c      get reciprocal lattice vectors 
         call smd_lat_invrt(lat,rlat)
  c       write(*,*) "Inverse Lattice vectors"
  c       write(*,*) (rlat(i,1),i=1,3)
  c       write(*,*) (rlat(i,2),i=1,3)
  c       write(*,*) (rlat(i,3),i=1,3)
  c      get lattice volume
         call smd_latt_vol(lat,vol)
  c
  