anchor.diff:         write(*,*) message
          stop
          end program
  +
  +      subroutine trj_read_coords(fn,nv,nu,cv,cu,ok)
  +      implicit none
  +      integer, intent(in)    :: fn  
  +      integer, intent(in)   :: nv
  +      integer, intent(in)   :: nu
  +      logical, intent(out)  :: ok
  +      double precision, intent(out)   :: cv(3,nv)
  +      double precision, intent(out)   :: cu(3,nu)
  +c      
  +      character*80 card
  +      character*80 message
  +      character*30 token
  +      character*1 sep
  +      character*1 a1
  +      integer i0
  +      integer i,k,lb,lt
  +      logical ostatus
  +c
  +      ok = .false.
  +      sep = " "
  +c      
  +c     look for the start of the fame record
  +c     -------------------------------------
  +  100 continue
  +      read(fn,1000,end=11,err=911) card
  +c      write(*,*) card
  + 1000 format(a)
  +      if(card(1:5).ne.'frame') goto 100
  +c
  +c     skip five lines here
  +c     --------------------
  +      do i=1,5
  +        read(fn,1000,err=911,end=911) card
  +      end do
  +      do i=1,nv
  +        read(fn,1000,err=911,end=911) card
  +        i0 = 1
  +        do k=1,3
  +          call get_next_token(i0,card,sep,token,ostatus)
  +          read(token,*,ERR=911,END=911) cv(k,i)
  +        end do
  +      end do
  +      do i=1,nu
  +        read(fn,1000,err=911,end=911) card
  +        i0 = 1
  +        do k=1,3
  +          call get_next_token(i0,card,sep,token,ostatus)
  +          read(token,*,ERR=911,END=911) cu(k,i)
  +        end do
  +      end do
  +      ok = .true.
  +11    continue
  +      return
  +911    continue       
  +c      if you reach this you are in trouble
  +       write(*,*) "found error"
  +       write(*,*) message
  +       stop
  +
  +      end subroutine
  +
  +      subroutine trj_read_header(fn,nvm,nva,nu,nprec)
  +      implicit none
  +      integer, intent(in)    :: fn  
  +      integer, intent(out)   :: nvm
  +      integer, intent(out)   :: nva
  +      integer, intent(out)   :: nu
  +      integer, intent(out)   :: nprec
  +c      
  +      character*80 card
  +      character*80 message
  +      character*30 token
  +      character*1 sep
  +      character*1 a1
  +      integer i0
  +      integer i,k,lb,lt
  +      logical ostatus
  +c
  +      rewind(fn)
  +c      
  +c     look for the "header"
  +c     ---------------------
  +  100 continue
  +      read(fn,1000,end=911) card
  +c      write(*,*) card
  + 1000 format(a)
  +      if(card(1:6).ne.'header') goto 100
  +c     parse the line      
  +      read(fn,1000,err=911,end=911) card
  +      i0 = 1
  +      sep = " "
  +      call get_next_token(i0,card,sep,token,ostatus)
  +      read(token,*) nva
  +      call get_next_token(i0,card,sep,token,ostatus)
  +      read(token,*) nu
  +      call get_next_token(i0,card,sep,token,ostatus)
  +      call get_next_token(i0,card,sep,token,ostatus)
  +      call get_next_token(i0,card,sep,token,ostatus)
  +      read(token,*) nvm
  +      call get_next_token(i0,card,sep,token,ostatus)
  +      read(token,*) nprec
  +      return
  +911    continue       
  +c      if you reach this you are in trouble
  +       write(*,*) "found error"
  +       write(*,*) message
  +       stop
  +
  +      end subroutine
  +
  +      subroutine get_next_token(i0,buffer,sep,token,ostatus)
  +      implicit none
  +      integer ,         intent(inout) :: i0
  +      character(LEN=*), intent(in)    :: buffer
  +      character(LEN=*), intent(out)   :: token
  +      character*1,      intent(in)    :: sep
  +      logical,          intent(out)   :: ostatus
  +c      
  +      character*1 a1
  +      integer i,k,lb,lt
  +c      
  +      token = " "
  +      ostatus = .false.
  +      lb=len(buffer)
  +      lt=len(token)
  +      k=0
  +      do i=i0,lb
  +       a1=buffer(i:i)
  +       if(a1.eq." ") then
  +          if(k.gt.0.and.sep.eq." ") goto 30
  +       else if(a1.eq.sep.and.sep.ne." ") then
  +          goto 30
  +       else
  +         k = k+1
  +         if(k.gt.lt) goto 30
  +         token(k:k) = a1
  +       end if
  +      end do
  +30    continue
  +      i0=i+1
  +      ostatus = k.gt.0
  +      end subroutine
  +
  +       subroutine lattice_read_file(fname,lat)
  +       implicit none
  +       character*(*) fname
  +       double precision lat(3,3)
  +c       
  +       integer ifn
  +       logical ofile
  +       character*180 bigbuf
  +c       
  +       character*180 message
  +       character*30 pname
  +       integer i,k
  +       logical is_integer
  +       external is_integer
  +c
  +       pname = "lattice_read_file"
  +       write(*,*) "reading lattice from file"
  +       lat = 0.0d0
  +       ifn = 13
  +       inquire(file=fname,exist=ofile)
  +       if(ofile) then
  +               open(ifn,file=fname,
  +     $            form='formatted',status='old',err=911)
  +       else
  +           message = "no file found: "//fname
  +           goto 911
  +       end if 
  +c 
  +       i = 0
  +       do 
  +         message = "reading lattice file"
  +         read(ifn,20,end=10,err=911) bigbuf
  +c        increment line number
  +         if(bigbuf.eq." ") cycle
  +         i = i+1
  +         read(bigbuf,*,err=10) (lat(k,i),k=1,3)
  +         if(i.eq.3) exit
  +       end do
  +10     continue
  +       if(i.eq.0.or.i.eq.2) then
  +         message = "something wrong with format"
  +         goto 911
  +       end if
  +       if(i.eq.1) then
  +          lat(2,2) = lat(2,1)
  +          lat(3,3) = lat(3,1)
  +          lat(2,1) = 0.00d0
  +          lat(3,1) = 0.00d0
  +       end if
  +       return
  +c
  +911    continue       
  +c      if you reach this you are in trouble
  +       write(*,*) "Emergency STOP in"
  +       write(*,*) pname(1:len_trim(pname))
  +       write(*,*) message(1:len_trim(message))
  +       stop
  +c
  +c      format statements
  +c      -----------------
  +30     FORMAT(180A1)
  +20     FORMAT(A180)
  +       end
   
          subroutine xyz_read_natoms(n,fn_in)
          implicit none
          integer n
          integer fn_in
   c       
  +       character*30 pname
          character*180 bigbuf
  -       character*180 message
  -       integer l
          logical is_integer
          external is_integer
  +c
  +       pname = "xyz_read_natoms"
  +       n = 0
   c      ------------------------------------------
   c      get number of atoms (skipping empty lines)
   c      ------------------------------------------
  -11     continue
  -       message = "looking for number of atoms"
  -       read(fn_in,20,end=138) bigbuf
  -c      increment line number
  -       l=l+1
  -       if(.not.is_integer(bigbuf)) goto 11
  +       bigbuf = " "
  +       do 
  +         read(fn_in,20,end=10) bigbuf
  +         if(is_integer(bigbuf)) exit
  +       end do
          read(bigbuf,*,err=136) n
  +10     continue
          return
   c      -------------
   c      error section
   c      -------------
   136    continue
  -       write(*,*) "error number of atoms field"
  -       write(*,*) "line number: ", l
  -       write(*,*) "buffer: ",bigbuf
  -       stop
  -       return
  -138    continue
  -       write(*,*) message
  +       write(*,*) "error in ",pname(1:len_trim(pname))
  +       write(*,*) "current buffer: ",bigbuf(1:len_trim(bigbuf))
          stop
   c
   c      format statements
   c      -----------------
  -30     FORMAT(180A1)
   20     FORMAT(A180)
          end
   
  
positive.code:        integer nf
         integer nfrm
  c      allocatable arrays
         double precision, dimension(:,:), allocatable :: c
         double precision, dimension(:,:), allocatable :: c1
         double precision, dimension(:,:), allocatable :: c2
         double precision, dimension(:), allocatable :: gr0
         double precision, dimension(:), allocatable :: gr
  c
  c      --------------------------------------------      
  c      beging parsing command line arguments if any
  
negative.code:           call get_command_argument(i,buffer,l,istatus)
            if(istatus.ne.0) goto 18
            if(is_number(buffer)) then
  c            message = "Only atom tags are supported now"
  c            goto 911
              read(buffer,*) atom2_id
            else
              atom2_tag = buffer
  