anchor.diff:      double precision, allocatable, dimension(:,:) :: selfl21
       double precision, allocatable, dimension(:,:) :: selfl22
   
  -    ! Pre-computed sine terms
  +    ! Pre-computed terms
       double precision, allocatable, dimension(:,:,:) :: ksi1
       double precision, allocatable, dimension(:,:,:) :: ksi2
  +
  +    ! Pre-computed terms
  +    double precision, allocatable, dimension(:,:,:,:) :: ksi31
  +    double precision, allocatable, dimension(:,:,:,:) :: ksi32
   
       double precision, allocatable, dimension(:,:,:,:,:) :: sinp1
       double precision, allocatable, dimension(:,:,:,:,:) :: sinp2
  
positive.code:     !$OMP PARALLEL DO PRIVATE(ni)
      do i = 1, nm1
          ni = n1(i)
          do m_1 = 1, ni
              do i_1 = 1, ni
                  if (x1(i,i_1,1,m_1) < cut_distance) then
                      ksi1(i, i_1, m_1) = 1.0d0 - sin(x1(i,i_1,1,m_1) * inv_cut)
                  endif
              enddo
          enddo
  
negative.code:     r_width2 = r_width**2
      c_width2 = c_width**2
  
      aadist = pd(int(x1(2,1)), int(x2(2,1)))
  
      ! ang_norm1 = 1.0d0 / (sum(ksi1(:N1)) * sum(ksi2(:N2)))
  
      ! write (*,*) "ANG_NORM1", ang_norm1, sum(ksi1(:N1)), sum(ksi2(:N2))
  
      do m_1 = 2, N1
  
          if (X1(1, m_1) > cut_distance) exit
  
          do m_2 = 2, N2
  
              ! if ((m_1.eq.1).and.(m_2.eq.1)) cycle
  
              if (X2(1, m_2) > cut_distance) exit
  
  