anchor.diff:      c_width2 = c_width**2
   
       inv_cut = pi / (2.0d0 * cut_distance)
  -    inv_sigma2(:) = -1.0d0 / (sigmas(:))**2
   
       allocate(sin1(na1, maxval(n1)))
       allocate(sin2(na2, maxval(n2)))
  
positive.code:     deallocate(sin1)
      deallocate(sin2)
  
  end subroutine fget_atomic_kernels_arad
  
  
  subroutine fget_atomic_distance_arad(q1, q2, z1, z2, n1, n2, na1, na2, &
          & width, cut_distance, r_width, c_width, kernel)
  
      use arad, only: atomic_distl2
  
      implicit none
  
      ! ARAD descriptors for each atom in the training set, format (i,5,m_1)
      double precision, dimension(:,:,:), intent(in) :: q1
      double precision, dimension(:,:,:), intent(in) :: q2
  
      ! ARAD atom-types for each atom, format (i, 2)
      double precision, dimension(:,:), intent(in) :: z1
      double precision, dimension(:,:), intent(in) :: z2
  
      ! List of numbers of atoms in each molecule
      integer, dimension(:), intent(in) :: n1
      integer, dimension(:), intent(in) :: n2
  
      ! Number of atom
      integer, intent(in) :: na1
      integer, intent(in) :: na2
  
      ! ARAD parameters
      double precision, intent(in) :: width
      double precision, intent(in) :: cut_distance
      double precision, intent(in) :: r_width
      double precision, intent(in) :: c_width
  
      ! Resulting alpha vector
      double precision, dimension(na1,na2), intent(out) :: kernel
  
      ! Internal counters
      integer :: i, j, k, ni, nj
      integer :: m_1
  
      ! Pre-computed constants
      double precision :: r_width2
      double precision :: c_width2
      double precision :: inv_cut
  
      ! Temporary variables necessary for parallelization
      double precision :: l2dist
      double precision, allocatable, dimension(:,:) :: atomic_distance
  
      ! Pre-computed terms in the full distance matrix
      double precision, allocatable, dimension(:) :: selfl21
      double precision, allocatable, dimension(:) :: selfl22
  
      ! Pre-computed sine terms
      double precision, allocatable, dimension(:,:) :: sin1
      double precision, allocatable, dimension(:,:) :: sin2
  
      ! Value of PI at full FORTRAN precision.
      double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
  
      ! Small number // to test for numerical stability
      double precision, parameter :: eps = 5.0e-12
  
      r_width2 = r_width**2
      c_width2 = c_width**2
  
      inv_cut = pi / (2.0d0 * cut_distance)
  
      allocate(sin1(na1, maxval(n1)))
      allocate(sin2(na2, maxval(n2)))
  
      !$OMP PARALLEL DO PRIVATE(ni)
      do i = 1, na1
          ni = n1(i)
          do m_1 = 1, ni
              sin1(i, m_1) = 1.0d0 - sin(q1(i,1,m_1) * inv_cut)
          enddo
      enddo
      !$OMP END PARALLEL DO
  
      !$OMP PARALLEL DO PRIVATE(ni)
      do i = 1, na2
          ni = n2(i)
          do m_1 = 1, ni
              sin2(i, m_1) = 1.0d0 - sin(q2(i,1,m_1) * inv_cut)
          enddo
      enddo
      !$OMP END PARALLEL DO
  
      allocate(selfl21(na1))
      allocate(selfl22(na2))
  
      !$OMP PARALLEL DO PRIVATE(ni)
      do i = 1, na1
          selfl21(i) = atomic_distl2(q1(i,:,:), q1(i,:,:), n1(i), n1(i), &
              & sin1(i,:), sin1(i,:), width, cut_distance, r_width, c_width)
      enddo
      !$OMP END PARALLEL DO
  
      !$OMP PARALLEL DO PRIVATE(ni)
      do i = 1, na2
          selfl22(i) = atomic_distl2(q2(i,:,:), q2(i,:,:), n2(i), n2(i), &
              & sin2(i,:), sin2(i,:), width, cut_distance, r_width, c_width)
      enddo
      !$OMP END PARALLEL DO
  
      kernel(:,:) = 0.0d0
  
      !$OMP PARALLEL DO PRIVATE(l2dist)
      do j = 1, na2
          do i = 1, na1
  
              l2dist = atomic_distl2(q1(i,:,:), q2(j,:,:), n1(i), n2(j), &
                  & sin1(i,:), sin2(j,:), width, cut_distance, r_width, c_width)
  
              l2dist = selfl21(i) + selfl22(j) - 2.0d0 * l2dist &
                  & * (r_width2/(r_width2 + (z1(i,1) - z2(j,1))**2) * &
                  & c_width2/(c_width2 + (z1(i,2) - z2(j,2))**2))
  
              if (abs(l2dist) < eps) l2dist = 0.0d0
  
              kernel(j, i) =  l2dist
  
          enddo
      enddo
      !$OMP END PARALLEL DO
  
      deallocate(selfl21)
      deallocate(selfl22)
      deallocate(sin1)
      deallocate(sin2)
  
  end subroutine fget_atomic_distance_arad
  
negative.code:             do i_1 = 1, ni
                  do j_1 = 1, nj
  
                      l2dist = atomic_distl2(q1(i,i_1,:,:), q1(j,j_1,:,:), n1(i), n1(j), &
                          & sin1(i,i_1,:), sin1(j,j_1,:), width, cut_distance, r_width, c_width)
  
                      l2dist = selfl21(i,i_1) + selfl21(j,j_1) - 2.0d0 * l2dist &
                          & * (r_width2/(r_width2 + (z1(i,i_1,1) - z1(j,j_1,1))**2) &
                          & * c_width2/(c_width2 + (z1(i,i_1,2) - z1(j,j_1,2))**2))
  
                      if (abs(l2dist) < eps) l2dist = 0.0d0
  
  