anchor.diff:        return
         end
   
  -      subroutine tm_score_all(ndim,nst,n_str,Lch,x,y,z,amat,rmsd_delta,
  +      subroutine tm_score_all(ndim,n_str,Lch,x,y,z,amat,rmsd_delta,
        &                        rmsd_a)
         implicit none
  +
  +      interface
  +      subroutine TMscore(L1,x1,y1,z1,n1,L2,x2,y2,z2,n2,TM,Rcomm,Lcomm)
  +          integer L1, L2
  +          real, allocatable, intent(inout) :: x1(:),y1(:),z1(:)
  +          integer, allocatable, intent(in) :: n1(:)
  +          real, allocatable, intent(inout) :: x2(:),y2(:),z2(:)
  +          integer, allocatable, intent(in) :: n2(:)
  +          real TM, Rcomm
  +          integer Lcomm
  +      end subroutine TMscore
  +      end interface
  +
   !     Intent In
  -      integer ndim ! Longest dimensions of a coordinate array
  -      integer nst ! Maximum number of structures
  +      integer ndim
         integer n_str ! Number of structures in use
         integer Lch ! Chain length of the structures (all same)
  -      real x(ndim,nst),y(ndim,nst),z(ndim,nst)
  +c      real x(ndim,nst),y(ndim,nst),z(ndim,nst)
  +      real, allocatable, intent(in) :: x(:,:),y(:,:), z(:,:)
  +      real, allocatable, intent(inout) :: amat(:,:)
   !     Intent Out
  -      real amat(nst,nst)
  +c      real amat(nst,nst)
         real rmsd_delta
         double precision rmsd_a
   
  
positive.code:       if (score_type .eq. 1) then
            call read_score_matrix(nst,n_str,amat,rmsd_delta, rmsd_a)
        elseif (score_type .eq. 2) then
            call tm_score_all(ndim,nst,n_str,Lch,x,y,z,amat,rmsd_delta,
       &                      rmsd_a)
        else
            call rmsd_score_all(ndim,nst,w,n_str,Lch,x,y,z,amat,
       &                        rmsd_delta,rmsd_a)
        endif
  
  
negative.code:       n_rmsd=0
  !$OMP PARALLEL DEFAULT(NONE)
  !$OMP& SHARED(n_str,Lch,x,y,z,amat,n_rmsd,rmsd_a,rmsd2_a)
  !$OMP& PRIVATE(i,j,k,x1,y1,z1,x2,y2,z2,nn1,nn2,TM,Rcomm,Lcomm,rmsd)
  !$OMP DO
  ! which is the correct place?
  !!$OMP& PRIVATE(i,j,k,x1,y1,z1,x2,y2,z2,nn1,nn2,TM,Rcomm,Lcomm,rmsd)
  !$OMP& REDUCTION(+:n_rmsd,rmsd_a,rmsd2_a)
  !$OMP& COLLAPSE(2)
        do i=1,n_str
           do j=1,n_str
  c OMP the COLLAPSE(2) directive unrolls both loops so we need to loop
  c over all indices. continue doesn't seem to work so we just protect the
  c j lt i with the if statement
              if (j .ge. i) then
                  if (j .eq. i) then
                      rmsd = 0.0
  