Linear kernel

qml/arad.py

@@ -24,6 +24,9 @@
 
 import numpy as np
 
+from .farad_kernels import fget_global_kernels_arad
+from .farad_kernels import fget_global_symmetric_kernels_arad
+
 from .farad_kernels import fget_local_kernels_arad
 from .farad_kernels import fget_local_symmetric_kernels_arad
 
@@ -131,6 +134,85 @@ def generate_arad_representation(coordinates, nuclear_charges, size=23, cut_dist
 
     return M
 
+def get_global_kernels_arad(X1, X2, sigmas, 
+        width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5):
+    """ Calculates the global Gaussian kernel matrix K for atomic ARAD
+        descriptors for a list of different sigmas. Each kernel element
+        is the sum of all kernel elements between pairs of atoms in two molecules.
+
+        K is calculated using an OpenMP parallel Fortran routine.
+
+        :param X1: ARAD descriptors for molecules in set 1.
+        :type X1: numpy array
+        :param X2: Array of ARAD descriptors for molecules in set 2.
+        :type X2: numpy array
+        :param sigmas: List of sigmas for which to calculate the Kernel matrices.
+        :type sigmas: list
+
+        :return: The kernel matrices for each sigma - shape (number_sigmas, number_molecules1, number_molecules2)
+        :rtype: numpy array
+    """
+
+    amax = X1.shape[1]
+
+    assert X1.shape[3] == amax, "ERROR: Check ARAD decriptor sizes! code = 1"
+    assert X2.shape[1] == amax, "ERROR: Check ARAD decriptor sizes! code = 2"
+    assert X2.shape[3] == amax, "ERROR: Check ARAD decriptor sizes! code = 3"
+
+    nm1 = X1.shape[0]
+    nm2 = X2.shape[0]
+
+    N1 = np.empty(nm1, dtype = np.int32)
+    Z1_arad = np.zeros((nm1, amax, 2))
+    for i in range(nm1):
+        N1[i] = len(np.where(X1[i,:,2,0] > 0)[0])
+        Z1_arad[i] = X1[i,:,1:3,0]
+
+    N2 = np.empty(nm2, dtype = np.int32)
+    Z2_arad = np.zeros((nm2, amax, 2))
+    for i in range(nm2):
+        N2[i] = len(np.where(X2[i,:,2,0] > 0)[0])
+        Z2_arad[i] = X2[i,:,1:3,0]
+
+    sigmas = np.array(sigmas)
+    nsigmas = sigmas.size
+
+    return fget_global_kernels_arad(X1, X2, Z1_arad, Z2_arad, N1, N2, sigmas, 
+                nm1, nm2, nsigmas, width, cut_distance, r_width, c_width)
+
+
+def get_global_symmetric_kernels_arad(X1, sigmas, 
+    width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5):
+    """ Calculates the global Gaussian kernel matrix K for atomic ARAD
+        descriptors for a list of different sigmas. Each kernel element
+        is the sum of all kernel elements between pairs of atoms in two molecules.
+
+        K is calculated using an OpenMP parallel Fortran routine.
+
+        :param X1: ARAD descriptors for molecules in set 1.
+        :type X1: numpy array
+        :param sigmas: List of sigmas for which to calculate the Kernel matrices.
+        :type sigmas: list
+
+        :return: The kernel matrices for each sigma - shape (number_sigmas, number_molecules1, number_molecules1)
+        :rtype: numpy array
+    """
+
+    nm1 = X1.shape[0]
+    amax = X1.shape[1]
+
+    N1 = np.empty(nm1, dtype = np.int32)
+    Z1_arad = np.zeros((nm1, amax, 2))
+    for i in range(nm1):
+        N1[i] = len(np.where(X1[i,:,2,0] > 0)[0])
+        Z1_arad[i] = X1[i,:,1:3,0]
+
+    sigmas = np.array(sigmas)
+    nsigmas = sigmas.size
+
+    return fget_global_symmetric_kernels_arad(X1, Z1_arad, N1, sigmas, 
+                nm1, nsigmas, width, cut_distance, r_width, c_width)
+
 
 def get_local_kernels_arad(X1, X2, sigmas, 
         width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5):

qml/farad_kernels.f90

@@ -82,7 +82,6 @@ function atomic_distl2(X1, X2, N1, N2, sin1, sin2, width, cut_distance, r_width,
 
                 aadist = aadist + d * (1.0d0 + x1(4,m_1)*x2(4,m_2) + x1(5,m_1)*x2(5,m_2))
 
-                ! write (*,*) m_1, m_2, x1(4,m_1), x2(4,m_2)
             end if
         end do
     end do
@@ -157,9 +156,6 @@ subroutine fget_local_kernels_arad(q1, q2, z1, z2, n1, n2, sigmas, nm1, nm2, nsi
     ! 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
 
@@ -172,8 +168,6 @@ subroutine fget_local_kernels_arad(q1, q2, z1, z2, n1, n2, sigmas, nm1, nm2, nsi
     sin1 = 0.0d0
     sin2 = 0.0d0
 
-    ! write (*,*) "INV_CUT", inv_cut
-
     !$OMP PARALLEL DO PRIVATE(ni)
     do i = 1, nm1
         ni = n1(i)
@@ -247,10 +241,7 @@ subroutine fget_local_kernels_arad(q1, q2, z1, z2, n1, n2, sigmas, nm1, nm2, nsi
                         & * (r_width2/(r_width2 + (z1(i,i_1,1) - z2(j,j_1,1))**2) * &
                         & c_width2/(c_width2 + (z1(i,i_1,2) - z2(j,j_1,2))**2))
 
-                    if (abs(l2dist) < eps) l2dist = 0.0d0
-
                     atomic_distance(i_1,j_1) = l2dist
-                    ! write (*,*) i_1, j_1, l2dist
 
                 enddo
             enddo
@@ -331,9 +322,6 @@ subroutine fget_local_symmetric_kernels_arad(q1, z1, n1, sigmas, nm1, nsigmas, &
     ! 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
 
@@ -388,8 +376,6 @@ subroutine fget_local_symmetric_kernels_arad(q1, z1, n1, sigmas, nm1, nsigmas, &
                         & * (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
-
                     atomic_distance(i_1,j_1) = l2dist
 
                 enddo
@@ -474,9 +460,6 @@ subroutine fget_atomic_kernels_arad(q1, q2, z1, z2, n1, n2, sigmas, na1, na2, ns
     ! 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
 
@@ -534,8 +517,6 @@ subroutine fget_atomic_kernels_arad(q1, q2, z1, z2, n1, n2, sigmas, na1, na2, ns
                 & * (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
-
             do k = 1, nsigmas
                 kernels(k, i, j) =  exp(l2dist * inv_sigma2(k))
             enddo
@@ -610,9 +591,6 @@ subroutine fget_atomic_symmetric_kernels_arad(q1, z1, n1, sigmas, na1, nsigmas,
     ! 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
 
@@ -652,8 +630,6 @@ subroutine fget_atomic_symmetric_kernels_arad(q1, z1, n1, sigmas, na1, nsigmas,
                 & * (r_width2/(r_width2 + (z1(i,1) - z1(j,1))**2) * &
                 & c_width2/(c_width2 + (z1(i,2) - z1(j,2))**2))
 
-            if (abs(l2dist) < eps) l2dist = 0.0d0
-
             do k = 1, nsigmas
                 kernels(k, i, j) =  exp(l2dist * inv_sigma2(k))
                 kernels(k, j, i) =  exp(l2dist * inv_sigma2(k))
@@ -667,3 +643,309 @@ subroutine fget_atomic_symmetric_kernels_arad(q1, z1, n1, sigmas, na1, nsigmas,
     deallocate(sin1)
 
 end subroutine fget_atomic_symmetric_kernels_arad
+
+
+subroutine fget_global_symmetric_kernels_arad(q1, z1, n1, sigmas, nm1, nsigmas, &
+        & width, cut_distance, r_width, c_width, kernels)
+
+    use arad, only: atomic_distl2
+
+    implicit none
+
+    ! ARAD descriptors for the training set, format (i,j_1,5,m_1)
+    double precision, dimension(:,:,:,:), intent(in) :: q1
+
+    ! ARAD atom-types for each atom in each molecule, format (i, j_1, 2)
+    double precision, dimension(:,:,:), intent(in) :: z1
+
+    ! List of numbers of atoms in each molecule
+    integer, dimension(:), intent(in) :: n1
+
+    ! Sigma in the Gaussian kernel
+    double precision, dimension(:), intent(in) :: sigmas
+
+    ! Number of molecules
+    integer, intent(in) :: nm1
+
+    ! Number of sigmas
+    integer, intent(in) :: nsigmas
+
+    ! -1.0 / sigma^2 for use in the kernel
+    double precision, dimension(nsigmas) :: inv_sigma2
+
+    ! 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(nsigmas,nm1,nm1), intent(out) :: kernels
+
+    ! Internal counters
+    integer :: i, j, k, ni, nj
+    integer :: m_1, i_1, j_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
+
+    ! Pre-computed sine terms
+    double precision, allocatable, dimension(:,:,:) :: sin1
+
+    ! Value of PI at full FORTRAN precision.
+    double precision, parameter :: pi = 4.0d0 * atan(1.0d0)
+
+    r_width2 = r_width**2
+    c_width2 = c_width**2
+
+    inv_cut = pi / (2.0d0 * cut_distance)
+    inv_sigma2(:) = -1.0d0 / (sigmas(:))**2
+
+    allocate(sin1(nm1, maxval(n1), maxval(n1)))
+
+    !$OMP PARALLEL DO PRIVATE(ni)
+    do i = 1, nm1
+        ni = n1(i)
+        do m_1 = 1, ni
+            do i_1 = 1, ni
+                sin1(i, i_1, m_1) = 1.0d0 - sin(q1(i,i_1,1,m_1) * inv_cut)
+            enddo
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    allocate(selfl21(nm1, maxval(n1)))
+
+    !$OMP PARALLEL DO PRIVATE(ni)
+    do i = 1, nm1
+        ni = n1(i)
+        do i_1 = 1, ni
+            selfl21(i,i_1) = atomic_distl2(q1(i,i_1,:,:), q1(i,i_1,:,:), n1(i), n1(i), &
+                & sin1(i,i_1,:), sin1(i,i_1,:), width, cut_distance, r_width, c_width)
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    allocate(atomic_distance(maxval(n1), maxval(n1)))
+
+    kernels(:,:,:) = 0.0d0
+    atomic_distance(:,:) = 0.0d0
+
+    !$OMP PARALLEL DO PRIVATE(l2dist,atomic_distance,ni,nj) schedule(dynamic)
+    do j = 1, nm1
+        nj = n1(j)
+        do i = 1, j
+            ni = n1(i)
+
+            atomic_distance(:,:) = 0.0d0
+
+            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))
+
+                    atomic_distance(i_1,j_1) = l2dist
+
+                enddo
+            enddo
+
+            do k = 1, nsigmas
+                kernels(k, i, j) =  exp(sum(atomic_distance(:ni,:nj)) * inv_sigma2(k))
+                kernels(k, j, i) =  kernels(k, i, j)
+            enddo
+
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    deallocate(atomic_distance)
+    deallocate(selfl21)
+    deallocate(sin1)
+
+end subroutine fget_global_symmetric_kernels_arad
+
+
+subroutine fget_global_kernels_arad(q1, q2, z1, z2, n1, n2, sigmas, nm1, nm2, nsigmas, &
+        & width, cut_distance, r_width, c_width, kernels)
+
+    use arad, only: atomic_distl2
+
+    implicit none
+
+    ! ARAD descriptors for the training set, format (i,j_1,5,m_1)
+    double precision, dimension(:,:,:,:), intent(in) :: q1
+    double precision, dimension(:,:,:,:), intent(in) :: q2
+
+    ! ARAD atom-types for each atom in each molecule, format (i, j_1, 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
+
+    ! Sigma in the Gaussian kernel
+    double precision, dimension(:), intent(in) :: sigmas
+
+    ! Number of molecules
+    integer, intent(in) :: nm1
+    integer, intent(in) :: nm2
+
+    ! Number of sigmas
+    integer, intent(in) :: nsigmas
+
+    ! -1.0 / sigma^2 for use in the kernel
+    double precision, dimension(nsigmas) :: inv_sigma2
+
+    ! 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(nsigmas,nm1,nm2), intent(out) :: kernels
+
+    ! Internal counters
+    integer :: i, j, k, ni, nj
+    integer :: m_1, i_1, j_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)
+
+    r_width2 = r_width**2
+    c_width2 = c_width**2
+
+    inv_cut = pi / (2.0d0 * cut_distance)
+    inv_sigma2(:) = -1.0d0 / (sigmas(:))**2
+
+    allocate(sin1(nm1, maxval(n1), maxval(n1)))
+    allocate(sin2(nm2, maxval(n2), maxval(n2)))
+
+    sin1 = 0.0d0
+    sin2 = 0.0d0
+
+    !$OMP PARALLEL DO PRIVATE(ni)
+    do i = 1, nm1
+        ni = n1(i)
+        do m_1 = 1, ni
+            do i_1 = 1, ni
+                if (q1(i,i_1,1,m_1) < cut_distance) then
+                    sin1(i, i_1, m_1) = 1.0d0 - sin(q1(i,i_1,1,m_1) * inv_cut)
+                endif
+            enddo
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    !$OMP PARALLEL DO PRIVATE(ni)
+    do i = 1, nm2
+        ni = n2(i)
+        do m_1 = 1, ni
+            do i_1 = 1, ni
+                if (q2(i,i_1,1,m_1) < cut_distance) then
+                    sin2(i, i_1, m_1) = 1.0d0 - sin(q2(i,i_1,1,m_1) * inv_cut)
+                endif
+            enddo
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    allocate(selfl21(nm1, maxval(n1)))
+    allocate(selfl22(nm2, maxval(n2)))
+
+    !$OMP PARALLEL DO PRIVATE(ni)
+    do i = 1, nm1
+        ni = n1(i)
+        do i_1 = 1, ni
+            selfl21(i,i_1) = atomic_distl2(q1(i,i_1,:,:), q1(i,i_1,:,:), n1(i), n1(i), &
+                & sin1(i,i_1,:), sin1(i,i_1,:), width, cut_distance, r_width, c_width)
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    !$OMP PARALLEL DO PRIVATE(ni)
+    do i = 1, nm2
+        ni = n2(i)
+        do i_1 = 1, ni
+            selfl22(i,i_1) = atomic_distl2(q2(i,i_1,:,:), q2(i,i_1,:,:), n2(i), n2(i), &
+                & sin2(i,i_1,:), sin2(i,i_1,:), width, cut_distance, r_width, c_width)
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+
+    allocate(atomic_distance(maxval(n1), maxval(n2)))
+
+    kernels(:,:,:) = 0.0d0
+    atomic_distance(:,:) = 0.0d0
+
+    !$OMP PARALLEL DO PRIVATE(l2dist,atomic_distance,ni,nj) schedule(dynamic)
+    do j = 1, nm2
+        nj = n2(j)
+        do i = 1, nm1
+            ni = n1(i)
+
+            atomic_distance(:,:) = 0.0d0
+
+            do i_1 = 1, ni
+                do j_1 = 1, nj
+
+                    l2dist = atomic_distl2(q1(i,i_1,:,:), q2(j,j_1,:,:), n1(i), n2(j), &
+                        & sin1(i,i_1,:), sin2(j,j_1,:), width, cut_distance, r_width, c_width)
+
+                    l2dist = selfl21(i,i_1) + selfl22(j,j_1) - 2.0d0 * l2dist &
+                        & * (r_width2/(r_width2 + (z1(i,i_1,1) - z2(j,j_1,1))**2) * &
+                        & c_width2/(c_width2 + (z1(i,i_1,2) - z2(j,j_1,2))**2))
+
+                    atomic_distance(i_1,j_1) = l2dist
+
+                enddo
+            enddo
+
+            do k = 1, nsigmas
+                kernels(k, i, j) =  exp(sum(atomic_distance(:ni,:nj)) * inv_sigma2(k))
+            enddo
+
+        enddo
+    enddo
+    !$OMP END PARALLEL DO
+
+    deallocate(atomic_distance)
+    deallocate(selfl21)
+    deallocate(selfl22)
+    deallocate(sin1)
+    deallocate(sin2)
+
+end subroutine fget_global_kernels_arad

qml/fkernels.f90

@@ -1,6 +1,6 @@
 ! MIT License
 !
-! Copyright (c) 2016 Anders Steen Christensen
+! Copyright (c) 2016 Anders Steen Christensen, Lars A. Bratholm, Felix A. Faber
 !
 ! Permission is hereby granted, free of charge, to any person obtaining a copy
 ! of this software and associated documentation files (the "Software"), to deal
@@ -225,6 +225,30 @@ subroutine flaplacian_kernel(a, na, b, nb, k, sigma)
 end subroutine flaplacian_kernel
 
 
+subroutine flinear_kernel(a, na, b, nb, k)
+
+    implicit none
+
+    double precision, dimension(:,:), intent(in) :: a
+    double precision, dimension(:,:), intent(in) :: b
+
+    integer, intent(in) :: na, nb
+
+    double precision, dimension(:,:), intent(inout) :: k
+
+    integer :: i, j
+
+!$OMP PARALLEL DO
+    do i = 1, nb
+        do j = 1, na
+            k(j,i) = dot_product(a(:,j), b(:,i))
+        enddo
+    enddo
+!$OMP END PARALLEL DO
+
+end subroutine flinear_kernel
+
+
 subroutine fmatern_kernel_l2(a, na, b, nb, k, sigma, order)
 
     implicit none

qml/frepresentations.f90

@@ -230,7 +230,6 @@ subroutine fgenerate_local_coulomb_matrix(atomic_charges, coordinates, natoms, n
 
     double precision, allocatable, dimension(:, :, :) :: pair_distance_matrix
     double precision, allocatable, dimension(:, :) :: distance_matrix
-    double precision, allocatable, dimension(:, :) :: distance_matrix_tmp
 
     integer i, j, m, n, k
 

qml/kernels.py

@@ -1,6 +1,6 @@
 # MIT License
 #
-# Copyright (c) 2016 Anders Steen Christensen, Felix Faber
+# Copyright (c) 2016 Anders Steen Christensen, Felix A. Faber, Lars A. Bratholm
 #
 # Permission is hereby granted, free of charge, to any person obtaining a copy
 # of this software and associated documentation files (the "Software"), to deal
@@ -26,6 +26,7 @@
 
 from .fkernels import fgaussian_kernel
 from .fkernels import flaplacian_kernel
+from .fkernels import flinear_kernel
 from .fkernels import fsargan_kernel
 from .fkernels import fmatern_kernel_l2
 
@@ -37,9 +38,9 @@ def laplacian_kernel(A, B, sigma):
         Where :math:`A_{i}` and :math:`B_{j}` are representation vectors.
         K is calculated using an OpenMP parallel Fortran routine.
 
-        :param A: 2D array of descriptors - shape (N, representation size).
+        :param A: 2D array of representations - shape (N, representation size).
         :type A: numpy array
-        :param B: 2D array of descriptors - shape (M, representation size).
+        :param B: 2D array of representations - shape (M, representation size).
         :type B: numpy array
         :param sigma: The value of sigma in the kernel matrix.
         :type sigma: float
@@ -66,9 +67,9 @@ def gaussian_kernel(A, B, sigma):
         Where :math:`A_{i}` and :math:`B_{j}` are representation vectors.
         K is calculated using an OpenMP parallel Fortran routine.
 
-        :param A: 2D array of descriptors - shape (N, representation size).
+        :param A: 2D array of representations - shape (N, representation size).
         :type A: numpy array
-        :param B: 2D array of descriptors - shape (M, representation size).
+        :param B: 2D array of representations - shape (M, representation size).
         :type B: numpy array
         :param sigma: The value of sigma in the kernel matrix.
         :type sigma: float
@@ -87,6 +88,34 @@ def gaussian_kernel(A, B, sigma):
 
     return K
 
+def linear_kernel(A, B):
+    """ Calculates the linear kernel matrix K, where :math:`K_{ij}`:
+
+            :math:`K_{ij} = A_i \cdot B_j`
+
+        VWhere :math:`A_{i}` and :math:`B_{j}` are  representation vectors. 
+
+        K is calculated using an OpenMP parallel Fortran routine.
+
+        :param A: 2D array of representations - shape (N, representation size).
+        :type A: numpy array
+        :param B: 2D array of representations - shape (M, representation size).
+        :type B: numpy array
+
+        :return: The Gaussian kernel matrix - shape (N, M)
+        :rtype: numpy array
+    """
+
+    na = A.shape[0]
+    nb = B.shape[0]
+
+    K = np.empty((na, nb), order='F')
+
+    # Note: Transposed for Fortran
+    flinear_kernel(A.T, na, B.T, nb, K)
+
+    return K
+
 def sargan_kernel(A, B, sigma, gammas):
     """ Calculates the Sargan kernel matrix K, where :math:`K_{ij}`:
 
@@ -95,9 +124,9 @@ def sargan_kernel(A, B, sigma, gammas):
         Where :math:`A_{i}` and :math:`B_{j}` are representation vectors.
         K is calculated using an OpenMP parallel Fortran routine.
 
-        :param A: 2D array of descriptors - shape (N, representation size).
+        :param A: 2D array of representations - shape (N, representation size).
         :type A: numpy array
-        :param B: 2D array of descriptors - shape (M, representation size).
+        :param B: 2D array of representations - shape (M, representation size).
         :type B: numpy array
         :param sigma: The value of sigma in the kernel matrix.
         :type sigma: float
@@ -137,9 +166,9 @@ def matern_kernel(A, B, sigma, order = 0, metric = "l1"):
 
         K is calculated using an OpenMP parallel Fortran routine.
 
-        :param A: 2D array of descriptors - shape (N, representation size).
+        :param A: 2D array of representations - shape (N, representation size).
         :type A: numpy array
-        :param B: 2D array of descriptors - shape (M, representation size).
+        :param B: 2D array of representations - shape (M, representation size).
         :type B: numpy array
         :param sigma: The value of sigma in the kernel matrix.
         :type sigma: float

tests/test_arad.py

@@ -13,6 +13,9 @@
 from qml.arad import get_local_kernels_arad
 from qml.arad import get_local_symmetric_kernels_arad
 
+from qml.arad import get_global_kernels_arad
+from qml.arad import get_global_symmetric_kernels_arad
+
 from qml.arad import get_atomic_kernels_arad
 from qml.arad import get_atomic_symmetric_kernels_arad
 
@@ -70,7 +73,11 @@ def test_arad():
     assert np.allclose(K_local_symm, K_local_asymm), "Symmetry error in local kernels"
     assert np.invert(np.all(np.isnan(K_local_asymm))), "ERROR: ARAD local symmetric kernel contains NaN"
 
-    K_local_asymm = get_local_kernels_arad(X1[-4:], X1[:6], sigmas)
+    K_global_asymm = get_global_kernels_arad(X1, X1, sigmas)
+    K_global_symm = get_global_symmetric_kernels_arad(X1, sigmas)
+
+    assert np.allclose(K_global_symm, K_global_asymm), "Symmetry error in global kernels"
+    assert np.invert(np.all(np.isnan(K_global_asymm))), "ERROR: ARAD global symmetric kernel contains NaN"
 
     molid = 5
     X1 = generate_arad_representation(mols[molid].coordinates,

tests/test_kernels.py

@@ -25,11 +25,16 @@
 import sys
 import numpy as np
 import qml
-from qml.kernels import laplacian_kernel, gaussian_kernel, \
-    matern_kernel, sargan_kernel
+from qml.kernels import laplacian_kernel
+from qml.kernels import gaussian_kernel
+from qml.kernels import linear_kernel
+from qml.kernels import matern_kernel
+from qml.kernels import sargan_kernel
 
 def test_laplacian_kernel():
 
+    np.random.seed(666)
+
     n_train = 25
     n_test = 20
 
@@ -57,6 +62,8 @@ def test_laplacian_kernel():
 
 def test_gaussian_kernel():
 
+    np.random.seed(666)
+
     n_train = 25
     n_test = 20
 
@@ -82,7 +89,40 @@ def test_gaussian_kernel():
     # Check for symmetry:
     assert np.allclose(Ksymm, Ksymm.T), "Error in Gaussian kernel"
 
+
+def test_linear_kernel():
+
+    np.random.seed(666)
+
+    n_train = 25
+    n_test = 20
+
+    # List of dummy representations
+    X = np.random.rand(n_train, 1000)
+    Xs = np.random.rand(n_test, 1000)
+
+    sigma = 100.0
+
+    Ktest = np.zeros((n_train, n_test))
+
+    for i in range(n_train):
+        for j in range(n_test):
+            Ktest[i,j] = np.dot(X[i], Xs[j])
+
+    K = linear_kernel(X, Xs)
+
+    # Compare two implementations:
+    assert np.allclose(K, Ktest), "Error in linear kernel"
+
+    Ksymm = linear_kernel(X, X)
+
+    # Check for symmetry:
+    assert np.allclose(Ksymm, Ksymm.T), "Error in linear kernel"
+
 def test_matern_kernel():
+
+    np.random.seed(666)
+
     for metric in ("l1", "l2"):
         for order in (0, 1, 2):
             print(metric,order)
@@ -129,6 +169,9 @@ def matern(metric, order):
     assert np.allclose(Ksymm, Ksymm.T), "Error in Matern kernel"
 
 def test_sargan_kernel():
+
+    np.random.seed(666)
+
     for ngamma in (0, 1, 2):
         sargan(ngamma)
 
@@ -173,5 +216,6 @@ def sargan(ngamma):
 if __name__ == "__main__":
     test_laplacian_kernel()
     test_gaussian_kernel()
+    test_linear_kernel()
     test_matern_kernel()
     test_sargan_kernel()