anchor.diff:  
   end function m_dist
   
  -
   function distl2(X1, X2, Z1, Z2, N1, N2, width, &
       & cut_distance, r_width, c_width) result(D12)
   
  
positive.code:     double precision :: c_width2, r_width2
  
      double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  
      double precision :: dd
      double precision :: pair
  
      double precision :: rdist
      double precision :: cdist
  
  
      inv_cut = pi / (2.0d0 * cut_distance)
      inv_width = 1.0d0 / (4.0d0 * width**2)
  
negative.code:     do j_1 = 1, N1
          do j_2 = 1, N2
  
              aadist = 0.0d0
  
              do m_1 = 1, N1
  
                  ! if (X1(1, m_1, j_1) > cut_distance) exit
                  if (X1(j_1, 1,m_1) > cut_distance) exit
  
                  ! sin1 = 1.0d0 - sin(x1(1, m_1, j_1) * inv_cut)
                  sin1 = 1.0d0 - sin(x1(j_1,1,m_1) * inv_cut)
  
                  do m_2 = 1, N2
  
                      ! if (X2(1, m_2, j_2) > cut_distance) exit
                      if (X2(j_2,1, m_2) > cut_distance) exit
  
                      ! if (abs(X2(1, m_2, j_2) - X1(1, m_1, j_1)) < maxgausdist) then
                      if (abs(X2(j_2,1,m_2) - X1(j_1,1,m_1)) < maxgausdist) then
  
                          ! d = exp(-((x1(1, m_1 ,j_1)-x2(1, m_2, j_2))**2) * inv_width ) *  & 
                          !     & sin1 * (1.0d0 - sin(x2(1, m_2, j_2) * inv_cut))
                          d = exp(-((x1(j_1,1,m_1)-x2(j_2,1,m_2))**2) * inv_width ) *  & 
                              & sin1 * (1.0d0 - sin(x2(j_2,1,m_2) * inv_cut))
  
                          ! r_dist = abs(x1(2, m_1, j_1) - x2(2,m_2, j_2))
                          ! c_dist = abs(x1(3, m_1, j_1) - x2(3,m_2, j_2))
                          r_dist = abs(x1(j_1,2,m_1) - x2(j_2,2,m_2))
                          c_dist = abs(x1(j_1,3,m_1) - x2(j_2,3,m_2))
      
                          ! d = d * (r_width2/(r_width2 + r_dist**2) * c_width2/(c_width2 + c_dist**2))
                          d = d * (r_width2/(r_width2 + r_dist**2) * c_width2/(c_width2 + c_dist**2))
  
  
                          ! aadist = aadist + d * (1.0d0 + x1(4, m_1, j_1)*x2(4, m_2, j_2) + & 
                          !     & x1(5, m_1, j_1)*x2(5, m_2, j_2))
                          aadist = aadist + d * (1.0d0 + x1(j_1,4,m_1)*x2(j_2,4,m_2) + & 
                              & x1(j_1,5,m_1)*x2(j_2,5,m_2))
  
  