anchor.diff:      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)
  
positive.code: 
  end function m_dist
  
  
  function distl2(X1, X2, Z1, Z2, N1, N2, width, &
      & cut_distance, r_width, c_width) result(D12)
  
  
negative.code:     double precision, dimension(nmol1) :: D11
      double precision, dimension(nmol2) :: D22
  
      integer :: i, j
  
      do i = 1, nmol1
          ! D11(i) = distl2(X1(:,:,:,i), X1(:,:,:,i), Z1(i,:,:), Z1(i,:,:), N1(i), N1(i), &
          D11(i) = distl2(X1(i,:,:,:), X1(i,:,:,:), Z1(i,:,:), Z1(i,:,:), N1(i), N1(i), &
              & width, cut_distance, r_width, c_width)
      enddo 
      do i = 1, nmol2
          ! D22(i) = distl2(X2(:,:,:,i), X2(:,:,:,i), Z2(i,:,:), Z2(i,:,:), N2(i), N2(i), &
          D22(i) = distl2(X2(i,:,:,:), X2(i,:,:,:), Z2(i,:,:), Z2(i,:,:), N2(i), N2(i), &
              & width, cut_distance, r_width, c_width)
      enddo
  
      do j = 1, nmol2
          do i = 1, nmol1
              ! D12(i,j) = distl2(X1(:,:,:,i), X2(:,:,:,j), Z1(i,:,:), Z2(j,:,:), N1(i), N2(j), &
              D12(i,j) = distl2(X1(i,:,:,:), X2(j,:,:,:), Z1(i,:,:), Z2(j,:,:), N1(i), N2(j), &
                  & width, cut_distance, r_width, c_width)
  
              D12(i,j) = D11(i) + D22(j) - 2.0d0 * D12(i,j)
          enddo
      enddo
  
  end subroutine molecular_arad_l2_distance_all
  
  