diff --git a/Makefile b/Makefile index 302b49a..97ef22b 100644 --- a/Makefile +++ b/Makefile @@ -2,6 +2,7 @@ python=./env/bin/python mamba=mamba pkg=qmllib pip=./env/bin/pip +pytest=pytest j=1 .PHONY: build @@ -24,7 +25,7 @@ test: ${python} -m pytest -rs ./tests types: - ${python} -m monkeytype run $(which pytest) ./tests/ + ${python} -m monkeytype run $$(which ${pytest}) ./tests ${python} -m monkeytype list-modules | grep ${pkg} | parallel -j${j} "${python} -m monkeytype apply {}" cov: diff --git a/src/qmllib/kernels/distance.py b/src/qmllib/kernels/distance.py index a0dae08..5507ac0 100644 --- a/src/qmllib/kernels/distance.py +++ b/src/qmllib/kernels/distance.py @@ -1,9 +1,12 @@ +from typing import Union + import numpy as np +from numpy import ndarray from .fdistance import fl2_distance, fmanhattan_distance, fp_distance_double, fp_distance_integer -def manhattan_distance(A, B): +def manhattan_distance(A: ndarray, B: ndarray) -> ndarray: """Calculates the Manhattan distances, D, between two Numpy arrays of representations. @@ -37,7 +40,7 @@ def manhattan_distance(A, B): return D -def l2_distance(A, B): +def l2_distance(A: ndarray, B: ndarray) -> ndarray: """Calculates the L2 distances, D, between two Numpy arrays of representations. @@ -71,7 +74,7 @@ def l2_distance(A, B): return D -def p_distance(A, B, p=2): +def p_distance(A: ndarray, B: ndarray, p: Union[int, float] = 2) -> ndarray: """Calculates the p-norm distances between two Numpy arrays of representations. The value of the keyword argument ``p =`` sets the norm order. diff --git a/src/qmllib/kernels/gradient_kernels.py b/src/qmllib/kernels/gradient_kernels.py index d1a8675..59f7f05 100644 --- a/src/qmllib/kernels/gradient_kernels.py +++ b/src/qmllib/kernels/gradient_kernels.py @@ -1,4 +1,7 @@ +from typing import List, Union + import numpy as np +from numpy import ndarray from qmllib.utils.environment_manipulation import ( mkl_get_num_threads, @@ -22,7 +25,9 @@ ) -def get_global_kernel(X1, X2, Q1, Q2, SIGMA): +def get_global_kernel( + X1: ndarray, X2: ndarray, Q1: List[List[int]], Q2: List[List[int]], SIGMA: float +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: :math:`K_{ij} = \\sum_{I\\in i} \\sum_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -75,7 +80,9 @@ def get_global_kernel(X1, X2, Q1, Q2, SIGMA): return K -def get_local_kernels(X1, X2, Q1, Q2, SIGMAS): +def get_local_kernels( + X1: ndarray, X2: ndarray, Q1: List[List[int]], Q2: List[List[int]], SIGMAS: List[float] +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: :math:`K_{ij} = \\sum_{I\\in i} \\sum_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -131,7 +138,13 @@ def get_local_kernels(X1, X2, Q1, Q2, SIGMAS): return K -def get_local_kernel(X1, X2, Q1, Q2, SIGMA): +def get_local_kernel( + X1: ndarray, + X2: ndarray, + Q1: List[Union[ndarray, List[int]]], + Q2: List[Union[ndarray, List[int]]], + SIGMA: float, +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: :math:`K_{ij} = \\sum_{I\\in i} \\sum_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -184,7 +197,7 @@ def get_local_kernel(X1, X2, Q1, Q2, SIGMA): return K -def get_local_symmetric_kernels(X1, Q1, SIGMAS): +def get_local_symmetric_kernels(X1: ndarray, Q1: List[List[int]], SIGMAS: List[float]) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: :math:`K_{ij} = \\sum_{I\\in i} \\sum_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -229,7 +242,9 @@ def get_local_symmetric_kernels(X1, Q1, SIGMAS): return K -def get_local_symmetric_kernel(X1, Q1, SIGMA): +def get_local_symmetric_kernel( + X1: ndarray, Q1: List[Union[ndarray, List[int]]], SIGMA: float +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: :math:`K_{ij} = \\sum_{I\\in i} \\sum_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -273,7 +288,13 @@ def get_local_symmetric_kernel(X1, Q1, SIGMA): return K -def get_atomic_local_kernel(X1, X2, Q1, Q2, SIGMA): +def get_atomic_local_kernel( + X1: ndarray, + X2: ndarray, + Q1: List[Union[ndarray, List[int]]], + Q2: List[Union[ndarray, List[int]]], + SIGMA: float, +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: @@ -332,7 +353,14 @@ def get_atomic_local_kernel(X1, X2, Q1, Q2, SIGMA): return K -def get_atomic_local_gradient_kernel(X1, X2, dX2, Q1, Q2, SIGMA): +def get_atomic_local_gradient_kernel( + X1: ndarray, + X2: ndarray, + dX2: ndarray, + Q1: List[Union[ndarray, List[int]]], + Q2: List[Union[ndarray, List[int]]], + SIGMA: float, +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: @@ -412,7 +440,9 @@ def get_atomic_local_gradient_kernel(X1, X2, dX2, Q1, Q2, SIGMA): return K -def get_local_gradient_kernel(X1, X2, dX2, Q1, Q2, SIGMA): +def get_local_gradient_kernel( + X1: ndarray, X2: ndarray, dX2: ndarray, Q1: List[List[int]], Q2: List[List[int]], SIGMA: float +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: @@ -478,7 +508,15 @@ def get_local_gradient_kernel(X1, X2, dX2, Q1, Q2, SIGMA): return K -def get_gdml_kernel(X1, X2, dX1, dX2, Q1, Q2, SIGMA): +def get_gdml_kernel( + X1: ndarray, + X2: ndarray, + dX1: ndarray, + dX2: ndarray, + Q1: List[List[int]], + Q2: List[List[int]], + SIGMA: float, +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: @@ -560,7 +598,9 @@ def get_gdml_kernel(X1, X2, dX1, dX2, Q1, Q2, SIGMA): return K -def get_symmetric_gdml_kernel(X1, dX1, Q1, SIGMA): +def get_symmetric_gdml_kernel( + X1: ndarray, dX1: ndarray, Q1: List[List[int]], SIGMA: float +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: @@ -613,7 +653,15 @@ def get_symmetric_gdml_kernel(X1, dX1, Q1, SIGMA): return K -def get_gp_kernel(X1, X2, dX1, dX2, Q1, Q2, SIGMA): +def get_gp_kernel( + X1: ndarray, + X2: ndarray, + dX1: ndarray, + dX2: ndarray, + Q1: List[Union[ndarray, List[int]]], + Q2: List[Union[ndarray, List[int]]], + SIGMA: float, +) -> ndarray: """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: @@ -692,7 +740,9 @@ def get_gp_kernel(X1, X2, dX1, dX2, Q1, Q2, SIGMA): return K -def get_symmetric_gp_kernel(X1, dX1, Q1, SIGMA): +def get_symmetric_gp_kernel( + X1: ndarray, dX1: ndarray, Q1: List[Union[ndarray, List[int]]], SIGMA: float +) -> ndarray: """ This symmetric kernel corresponds to a Gaussian process regression (GPR) approach. diff --git a/src/qmllib/kernels/kernels.py b/src/qmllib/kernels/kernels.py index e567efb..6064bb0 100644 --- a/src/qmllib/kernels/kernels.py +++ b/src/qmllib/kernels/kernels.py @@ -1,4 +1,7 @@ +from typing import List, Union + import numpy as np +from numpy import float64, ndarray from .fkernels import ( fgaussian_kernel, @@ -15,7 +18,7 @@ ) -def wasserstein_kernel(A, B, sigma, p=1, q=1): +def wasserstein_kernel(A: ndarray, B: ndarray, sigma: float, p: int = 1, q: int = 1) -> ndarray: """Calculates the Wasserstein kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{(W_p(A_i, B_i))^q}{\\sigma} \\big)` @@ -45,7 +48,7 @@ def wasserstein_kernel(A, B, sigma, p=1, q=1): return K -def laplacian_kernel(A, B, sigma): +def laplacian_kernel(A: ndarray, B: ndarray, sigma: float) -> ndarray: """Calculates the Laplacian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_1}{\\sigma} \\big)` @@ -75,7 +78,7 @@ def laplacian_kernel(A, B, sigma): return K -def laplacian_kernel_symmetric(A, sigma): +def laplacian_kernel_symmetric(A: ndarray, sigma: float) -> ndarray: """Calculates the symmetric Laplacian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - A_j\\|_1}{\\sigma} \\big)` @@ -102,7 +105,7 @@ def laplacian_kernel_symmetric(A, sigma): return K -def gaussian_kernel(A, B, sigma): +def gaussian_kernel(A: ndarray, B: ndarray, sigma: float) -> ndarray: """Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\\sigma^2} \\big)` @@ -132,7 +135,7 @@ def gaussian_kernel(A, B, sigma): return K -def gaussian_kernel_symmetric(A, sigma): +def gaussian_kernel_symmetric(A: ndarray, sigma: float) -> ndarray: """Calculates the symmetric Gaussian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - A_j\\|_2^2}{2\\sigma^2} \\big)` @@ -159,7 +162,7 @@ def gaussian_kernel_symmetric(A, sigma): return K -def linear_kernel(A, B): +def linear_kernel(A: ndarray, B: ndarray) -> ndarray: """Calculates the linear kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = A_i \\cdot B_j` @@ -188,7 +191,12 @@ def linear_kernel(A, B): return K -def sargan_kernel(A, B, sigma, gammas): +def sargan_kernel( + A: ndarray, + B: ndarray, + sigma: Union[float, float64], + gammas: Union[ndarray, List[Union[int, float]], List[int]], +) -> ndarray: """Calculates the Sargan kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\| A_i - B_j \\|_1)}{\\sigma} \\big) \\big(1 + \\sum_{k} \\frac{\\gamma_{k} \\| A_i - B_j \\|_1^k}{\\sigma^k} \\big)` @@ -225,7 +233,9 @@ def sargan_kernel(A, B, sigma, gammas): return K -def matern_kernel(A, B, sigma, order=0, metric="l1"): +def matern_kernel( + A: ndarray, B: ndarray, sigma: float, order: int = 0, metric: str = "l1" +) -> ndarray: """Calculates the Matern kernel matrix K, where :math:`K_{ij}`: for order = 0: @@ -286,7 +296,9 @@ def matern_kernel(A, B, sigma, order=0, metric="l1"): return K -def get_local_kernels_gaussian(A, B, na, nb, sigmas): +def get_local_kernels_gaussian( + A: ndarray, B: ndarray, na: ndarray, nb: ndarray, sigmas: List[float] +) -> ndarray: """Calculates the Gaussian kernel matrix K, for a local representation where :math:`K_{ij}`: :math:`K_{ij} = \\sum_{a \\in i} \\sum_{b \\in j} \\exp \\big( -\\frac{\\|A_a - B_b\\|_2^2}{2\\sigma^2} \\big)` @@ -328,7 +340,9 @@ def get_local_kernels_gaussian(A, B, na, nb, sigmas): return fget_local_kernels_gaussian(A.T, B.T, na, nb, sigmas, nma, nmb, nsigmas) -def get_local_kernels_laplacian(A, B, na, nb, sigmas): +def get_local_kernels_laplacian( + A: ndarray, B: ndarray, na: ndarray, nb: ndarray, sigmas: List[float] +) -> ndarray: """Calculates the Local Laplacian kernel matrix K, for a local representation where :math:`K_{ij}`: :math:`K_{ij} = \\sum_{a \\in i} \\sum_{b \\in j} \\exp \\big( -\\frac{\\|A_a - B_b\\|_1}{\\sigma} \\big)` @@ -370,7 +384,7 @@ def get_local_kernels_laplacian(A, B, na, nb, sigmas): return fget_local_kernels_laplacian(A.T, B.T, na, nb, sigmas, nma, nmb, nsigmas) -def kpca(K, n=2, centering=True): +def kpca(K: ndarray, n: int = 2, centering: bool = True) -> ndarray: """Calculates `n` first principal components for the kernel :math:`K`. The PCA is calculated using an OpenMP parallel Fortran routine. diff --git a/src/qmllib/representations/arad/arad.py b/src/qmllib/representations/arad/arad.py index 7f76381..221fa5c 100644 --- a/src/qmllib/representations/arad/arad.py +++ b/src/qmllib/representations/arad/arad.py @@ -1,4 +1,7 @@ +from typing import List + import numpy as np +from numpy import ndarray from qmllib.utils.alchemy import PTP @@ -12,7 +15,7 @@ ) -def getAngle(sp, norms): +def getAngle(sp: ndarray, norms: ndarray) -> ndarray: epsilon = 10.0 * np.finfo(float).eps angles = np.zeros(sp.shape) mask1 = np.logical_and(np.abs(sp - norms) > epsilon, np.abs(norms) > epsilon) @@ -20,7 +23,9 @@ def getAngle(sp, norms): return angles -def generate_arad_representation(coordinates, nuclear_charges, size=23, cut_distance=5.0): +def generate_arad_representation( + coordinates: ndarray, nuclear_charges: ndarray, size: int = 23, cut_distance: float = 5.0 +) -> ndarray: """Generates a representation for the ARAD kernel module. :param coordinates: Input coordinates. @@ -116,7 +121,15 @@ 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): +def get_global_kernels_arad( + X1: ndarray, + X2: ndarray, + sigmas: List[float], + width: float = 0.2, + cut_distance: float = 5.0, + r_width: float = 1.0, + c_width: float = 0.5, +) -> ndarray: """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. @@ -177,8 +190,13 @@ def get_global_kernels_arad(X1, X2, sigmas, width=0.2, cut_distance=5.0, r_width def get_global_symmetric_kernels_arad( - X1, sigmas, width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5 -): + X1: ndarray, + sigmas: List[float], + width: float = 0.2, + cut_distance: float = 5.0, + r_width: float = 1.0, + c_width: float = 0.5, +) -> ndarray: """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. @@ -211,7 +229,15 @@ def get_global_symmetric_kernels_arad( ) -def get_local_kernels_arad(X1, X2, sigmas, width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5): +def get_local_kernels_arad( + X1: ndarray, + X2: ndarray, + sigmas: List[float], + width: float = 0.2, + cut_distance: float = 5.0, + r_width: float = 1.0, + c_width: float = 0.5, +) -> ndarray: """Calculates the 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. @@ -272,8 +298,13 @@ def get_local_kernels_arad(X1, X2, sigmas, width=0.2, cut_distance=5.0, r_width= def get_local_symmetric_kernels_arad( - X1, sigmas, width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5 -): + X1: ndarray, + sigmas: List[float], + width: float = 0.2, + cut_distance: float = 5.0, + r_width: float = 1.0, + c_width: float = 0.5, +) -> ndarray: """Calculates the 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. @@ -306,7 +337,15 @@ def get_local_symmetric_kernels_arad( ) -def get_atomic_kernels_arad(X1, X2, sigmas, width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5): +def get_atomic_kernels_arad( + X1: ndarray, + X2: ndarray, + sigmas: List[float], + width: float = 0.2, + cut_distance: float = 5.0, + r_width: float = 1.0, + c_width: float = 0.5, +) -> ndarray: """Calculates the Gaussian kernel matrix K for atomic ARAD descriptors for a list of different sigmas. For atomic properties, e.g. partial charges, chemical shifts, etc. @@ -366,8 +405,13 @@ def get_atomic_kernels_arad(X1, X2, sigmas, width=0.2, cut_distance=5.0, r_width def get_atomic_symmetric_kernels_arad( - X1, sigmas, width=0.2, cut_distance=5.0, r_width=1.0, c_width=0.5 -): + X1: ndarray, + sigmas: List[float], + width: float = 0.2, + cut_distance: float = 5.0, + r_width: float = 1.0, + c_width: float = 0.5, +) -> ndarray: """Calculates the Gaussian kernel matrix K for atomic ARAD descriptors for a list of different sigmas. For atomic properties, e.g. partial charges, chemical shifts, etc. diff --git a/src/qmllib/representations/arad/farad_kernels.f90 b/src/qmllib/representations/arad/farad_kernels.f90 index 0a441dd..9ec5da8 100644 --- a/src/qmllib/representations/arad/farad_kernels.f90 +++ b/src/qmllib/representations/arad/farad_kernels.f90 @@ -22,967 +22,956 @@ module arad - implicit none + implicit none contains -function atomic_distl2(X1, X2, N1, N2, sin1, sin2, width, cut_distance, r_width, c_width) result(aadist) + function atomic_distl2(X1, X2, N1, N2, sin1, sin2, width, cut_distance, r_width, c_width) result(aadist) - implicit none + implicit none - double precision, dimension(:,:), intent(in) :: X1 - double precision, dimension(:,:), intent(in) :: X2 + double precision, dimension(:, :), intent(in) :: X1 + double precision, dimension(:, :), intent(in) :: X2 - integer, intent(in) :: N1 - integer, intent(in) :: N2 + integer, intent(in) :: N1 + integer, intent(in) :: N2 - double precision, dimension(:), intent(in) :: sin1 - double precision, dimension(:), intent(in) :: sin2 + double precision, dimension(:), intent(in) :: sin1 + double precision, dimension(:), intent(in) :: sin2 - double precision, intent(in) :: width - double precision, intent(in) :: cut_distance - double precision, intent(in) :: r_width - double precision, intent(in) :: c_width + double precision, intent(in) :: width + double precision, intent(in) :: cut_distance + double precision, intent(in) :: r_width + double precision, intent(in) :: c_width - double precision :: aadist + double precision :: aadist - double precision :: d + double precision :: d - integer :: m_1, m_2 + integer :: m_1, m_2 - double precision :: maxgausdist2 + double precision :: maxgausdist2 - double precision :: inv_width - double precision :: c_width2, r_width2, r2 + double precision :: inv_width + double precision :: c_width2, r_width2, r2 - inv_width = -1.0d0 / (4.0d0 * width**2) + inv_width = -1.0d0/(4.0d0*width**2) - maxgausdist2 = (8.0d0 * width)**2 - r_width2 = r_width**2 - c_width2 = c_width**2 + maxgausdist2 = (8.0d0*width)**2 + r_width2 = r_width**2 + c_width2 = c_width**2 - aadist = 0.0d0 + aadist = 0.0d0 - do m_1 = 1, N1 + do m_1 = 1, N1 - if (X1(1, m_1) > cut_distance) exit + if (X1(1, m_1) > cut_distance) exit - do m_2 = 1, N2 + do m_2 = 1, N2 if (X2(1, m_2) > cut_distance) exit - r2 = (X2(1,m_2) - X1(1,m_1))**2 + r2 = (X2(1, m_2) - X1(1, m_1))**2 if (r2 < maxgausdist2) then - d = exp(r2 * inv_width ) * sin1(m_1) * sin2(m_2) + d = exp(r2*inv_width)*sin1(m_1)*sin2(m_2) - d = d * (r_width2/(r_width2 + (x1(2,m_1) - x2(2,m_2))**2) * & - & c_width2/(c_width2 + (x1(3,m_1) - x2(3,m_2))**2)) + d = d*(r_width2/(r_width2 + (x1(2, m_1) - x2(2, m_2))**2)* & + & c_width2/(c_width2 + (x1(3, m_1) - x2(3, m_2))**2)) - aadist = aadist + d * (1.0d0 + x1(4,m_1)*x2(4,m_2) + x1(5,m_1)*x2(5,m_2)) + aadist = aadist + d*(1.0d0 + x1(4, m_1)*x2(4, m_2) + x1(5, m_1)*x2(5, m_2)) end if - end do - end do + end do + end do -end function atomic_distl2 + end function atomic_distl2 end module arad - subroutine fget_local_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) = sum(exp(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) + use arad, only: atomic_distl2 -end subroutine fget_local_kernels_arad + 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) + end if + end do + end do + end do + !$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) + end if + end do + end do + end do + !$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) + end do + end do + !$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) + end do + end do + !$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 + + end do + end do + + do k = 1, nsigmas + kernels(k, i, j) = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) + end do + + end do + end do + !$OMP END PARALLEL DO + + deallocate (atomic_distance) + deallocate (selfl21) + deallocate (selfl22) + deallocate (sin1) + deallocate (sin2) + +end subroutine fget_local_kernels_arad subroutine fget_local_symmetric_kernels_arad(q1, z1, n1, sigmas, nm1, nsigmas, & & width, cut_distance, r_width, c_width, kernels) - use arad, only: atomic_distl2 + use arad, only: atomic_distl2 - implicit none + implicit none - ! ARAD descriptors for the training set, format (i,j_1,5,m_1) - double precision, dimension(:,:,:,:), intent(in) :: q1 + ! 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 + ! 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 + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: n1 - ! Sigma in the Gaussian kernel - double precision, dimension(:), intent(in) :: sigmas + ! Sigma in the Gaussian kernel + double precision, dimension(:), intent(in) :: sigmas - ! Number of molecules - integer, intent(in) :: nm1 + ! Number of molecules + integer, intent(in) :: nm1 - ! Number of sigmas - integer, intent(in) :: nsigmas + ! Number of sigmas + integer, intent(in) :: nsigmas - ! -1.0 / sigma^2 for use in the kernel - double precision, dimension(nsigmas) :: inv_sigma2 + ! -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 + ! 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 + ! 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 + ! 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 + ! 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 + ! 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 terms in the full distance matrix + double precision, allocatable, dimension(:, :) :: selfl21 - ! Pre-computed sine terms - double precision, allocatable, dimension(:,:,:) :: sin1 + ! 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) + ! 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 + r_width2 = r_width**2 + c_width2 = c_width**2 - inv_cut = pi / (2.0d0 * cut_distance) - inv_sigma2(:) = -1.0d0 / (sigmas(:))**2 + inv_cut = pi/(2.0d0*cut_distance) + inv_sigma2(:) = -1.0d0/(sigmas(:))**2 - allocate(sin1(nm1, maxval(n1), maxval(n1))) + 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 + !$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) + end do + end do + end do + !$OMP END PARALLEL DO - allocate(selfl21(nm1, maxval(n1))) + 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 + !$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) + end do + end do + !$OMP END PARALLEL DO - allocate(atomic_distance(maxval(n1), maxval(n1))) + allocate (atomic_distance(maxval(n1), maxval(n1))) - kernels(:,:,:) = 0.0d0 - atomic_distance(:,:) = 0.0d0 + 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) + !$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 + atomic_distance(:, :) = 0.0d0 - do i_1 = 1, ni - do j_1 = 1, nj + 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 = 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)) + 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 + atomic_distance(i_1, j_1) = l2dist - enddo - enddo + end do + end do - do k = 1, nsigmas - kernels(k, i, j) = sum(exp(atomic_distance(:ni,:nj) * inv_sigma2(k))) - kernels(k, j, i) = kernels(k, i, j) - enddo + do k = 1, nsigmas + kernels(k, i, j) = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) + kernels(k, j, i) = kernels(k, i, j) + end do - enddo - enddo - !$OMP END PARALLEL DO + end do + end do + !$OMP END PARALLEL DO - deallocate(atomic_distance) - deallocate(selfl21) - deallocate(sin1) + deallocate (atomic_distance) + deallocate (selfl21) + deallocate (sin1) end subroutine fget_local_symmetric_kernels_arad subroutine fget_atomic_kernels_arad(q1, q2, z1, z2, n1, n2, sigmas, na1, na2, nsigmas, & & width, cut_distance, r_width, c_width, kernels) - use arad, only: atomic_distl2 + 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 + 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 + ! ARAD atom-types for each atom, format (i, 2) + double precision, dimension(:, :), intent(in) :: z1 + double precision, dimension(:, :), intent(in) :: z2 - ! Sigma in the Gaussian kernel - double precision, dimension(:), intent(in) :: sigmas + ! Sigma in the Gaussian kernel + double precision, dimension(:), intent(in) :: sigmas - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: n1 - integer, dimension(:), intent(in) :: n2 + ! 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 + ! Number of atom + integer, intent(in) :: na1 + integer, intent(in) :: na2 - ! Number of sigmas - integer, intent(in) :: nsigmas + ! Number of sigmas + integer, intent(in) :: nsigmas - ! -1.0 / sigma^2 for use in the kernel - double precision, dimension(nsigmas) :: inv_sigma2 + ! -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 + ! 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,na1,na2), intent(out) :: kernels + ! Resulting alpha vector + double precision, dimension(nsigmas, na1, na2), intent(out) :: kernels - ! Internal counters - integer :: i, j, k, ni - integer :: m_1 + ! Internal counters + integer :: i, j, k, ni + integer :: m_1 - ! Pre-computed constants - double precision :: r_width2 - double precision :: c_width2 - double precision :: inv_cut + ! Pre-computed constants + double precision :: r_width2 + double precision :: c_width2 + double precision :: inv_cut - ! Temporary variables necessary for parallelization - double precision :: l2dist + ! Temporary variables necessary for parallelization + double precision :: l2dist - ! Pre-computed terms in the full distance matrix - double precision, allocatable, dimension(:) :: selfl21 - double precision, allocatable, dimension(:) :: selfl22 + ! 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 + ! 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) + ! 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 + r_width2 = r_width**2 + c_width2 = c_width**2 - inv_cut = pi / (2.0d0 * cut_distance) - inv_sigma2(:) = -1.0d0 / (sigmas(:))**2 + inv_cut = pi/(2.0d0*cut_distance) + inv_sigma2(:) = -1.0d0/(sigmas(:))**2 - allocate(sin1(na1, maxval(n1))) - allocate(sin2(na2, maxval(n2))) + 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, na1 + ni = n1(i) + do m_1 = 1, ni + sin1(i, m_1) = 1.0d0 - sin(q1(i, 1, m_1)*inv_cut) + end do + end do + !$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 + !$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) + end do + end do + !$OMP END PARALLEL DO - allocate(selfl21(na1)) - allocate(selfl22(na2)) + 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, na1 + selfl21(i) = atomic_distl2(q1(i, :, :), q1(i, :, :), n1(i), n1(i), & + & sin1(i, :), sin1(i, :), width, cut_distance, r_width, c_width) + end do + !$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 - - kernels(:,:,:) = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(l2dist) schedule(dynamic) - 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)) - - do k = 1, nsigmas - kernels(k, i, j) = exp(l2dist * inv_sigma2(k)) - enddo - - enddo - enddo - !$OMP END PARALLEL DO - - deallocate(selfl21) - deallocate(selfl22) - deallocate(sin1) - deallocate(sin2) + !$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) + end do + !$OMP END PARALLEL DO + + kernels(:, :, :) = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(l2dist) schedule(dynamic) + 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)) + + do k = 1, nsigmas + kernels(k, i, j) = exp(l2dist*inv_sigma2(k)) + end do + + end do + end do + !$OMP END PARALLEL DO + + deallocate (selfl21) + deallocate (selfl22) + deallocate (sin1) + deallocate (sin2) end subroutine fget_atomic_kernels_arad - subroutine fget_atomic_symmetric_kernels_arad(q1, z1, n1, sigmas, na1, nsigmas, & & width, cut_distance, r_width, c_width, kernels) - use arad, only: atomic_distl2 + use arad, only: atomic_distl2 - implicit none + implicit none - ! ARAD descriptors for each atom in the training set, format (i,5,m_1) - double precision, dimension(:,:,:), intent(in) :: q1 + ! ARAD descriptors for each atom in the training set, format (i,5,m_1) + double precision, dimension(:, :, :), intent(in) :: q1 - ! ARAD atom-types for each atom, format (i, 2) - double precision, dimension(:,:), intent(in) :: z1 + ! ARAD atom-types for each atom, format (i, 2) + double precision, dimension(:, :), intent(in) :: z1 - ! Sigma in the Gaussian kernel - double precision, dimension(:), intent(in) :: sigmas + ! Sigma in the Gaussian kernel + double precision, dimension(:), intent(in) :: sigmas - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: n1 + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: n1 - ! Number of atom - integer, intent(in) :: na1 + ! Number of atom + integer, intent(in) :: na1 - ! Number of sigmas - integer, intent(in) :: nsigmas + ! Number of sigmas + integer, intent(in) :: nsigmas - ! -1.0 / sigma^2 for use in the kernel - double precision, dimension(nsigmas) :: inv_sigma2 + ! -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 + ! 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,na1,na1), intent(out) :: kernels + ! Resulting alpha vector + double precision, dimension(nsigmas, na1, na1), intent(out) :: kernels - ! Internal counters - integer :: i, j, k, ni - integer :: m_1 + ! Internal counters + integer :: i, j, k, ni + integer :: m_1 - ! Pre-computed constants - double precision :: r_width2 - double precision :: c_width2 - double precision :: inv_cut + ! Pre-computed constants + double precision :: r_width2 + double precision :: c_width2 + double precision :: inv_cut - ! Temporary variables necessary for parallelization - double precision :: l2dist + ! Temporary variables necessary for parallelization + double precision :: l2dist - ! Pre-computed terms in the full distance matrix - double precision, allocatable, dimension(:) :: selfl21 + ! 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) + ! Pre-computed sine terms + double precision, allocatable, dimension(:, :) :: sin1 - r_width2 = r_width**2 - c_width2 = c_width**2 + ! Value of PI at full FORTRAN precision. + double precision, parameter :: pi = 4.0d0*atan(1.0d0) - inv_cut = pi / (2.0d0 * cut_distance) - inv_sigma2(:) = -1.0d0 / (sigmas(:))**2 + r_width2 = r_width**2 + c_width2 = c_width**2 - allocate(sin1(na1, maxval(n1))) + inv_cut = pi/(2.0d0*cut_distance) + inv_sigma2(:) = -1.0d0/(sigmas(:))**2 - !$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 + allocate (sin1(na1, maxval(n1))) - allocate(selfl21(na1)) + !$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) + end do + end do + !$OMP END PARALLEL DO - !$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 + allocate (selfl21(na1)) - kernels(:,:,:) = 0.0d0 + !$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) + end do + !$OMP END PARALLEL DO - !$OMP PARALLEL DO PRIVATE(l2dist) schedule(dynamic) - do j = 1, na1 - do i = j, na1 + kernels(:, :, :) = 0.0d0 - l2dist = atomic_distl2(q1(i,:,:), q1(j,:,:), n1(i), n1(j), & - & sin1(i,:), sin1(j,:), width, cut_distance, r_width, c_width) + !$OMP PARALLEL DO PRIVATE(l2dist) schedule(dynamic) + do j = 1, na1 + do i = j, na1 - l2dist = selfl21(i) + selfl21(j) - 2.0d0 * l2dist & - & * (r_width2/(r_width2 + (z1(i,1) - z1(j,1))**2) * & - & c_width2/(c_width2 + (z1(i,2) - z1(j,2))**2)) + l2dist = atomic_distl2(q1(i, :, :), q1(j, :, :), n1(i), n1(j), & + & sin1(i, :), sin1(j, :), width, cut_distance, r_width, c_width) - do k = 1, nsigmas - kernels(k, i, j) = exp(l2dist * inv_sigma2(k)) - kernels(k, j, i) = exp(l2dist * inv_sigma2(k)) - enddo + l2dist = selfl21(i) + selfl21(j) - 2.0d0*l2dist & + & *(r_width2/(r_width2 + (z1(i, 1) - z1(j, 1))**2)* & + & c_width2/(c_width2 + (z1(i, 2) - z1(j, 2))**2)) - enddo - enddo - !$OMP END PARALLEL DO + do k = 1, nsigmas + kernels(k, i, j) = exp(l2dist*inv_sigma2(k)) + kernels(k, j, i) = exp(l2dist*inv_sigma2(k)) + end do - deallocate(selfl21) - deallocate(sin1) + end do + end do + !$OMP END PARALLEL DO -end subroutine fget_atomic_symmetric_kernels_arad + deallocate (selfl21) + 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 + use arad, only: atomic_distl2 - implicit none + implicit none - ! ARAD descriptors for the training set, format (i,j_1,5,m_1) - double precision, dimension(:,:,:,:), intent(in) :: q1 + ! 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 + ! 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 + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: n1 - ! Sigma in the Gaussian kernel - double precision, dimension(:), intent(in) :: sigmas + ! Sigma in the Gaussian kernel + double precision, dimension(:), intent(in) :: sigmas - ! Number of molecules - integer, intent(in) :: nm1 + ! Number of molecules + integer, intent(in) :: nm1 - ! Number of sigmas - integer, intent(in) :: nsigmas + ! Number of sigmas + integer, intent(in) :: nsigmas - ! -1.0 / sigma^2 for use in the kernel - double precision, dimension(nsigmas) :: inv_sigma2 + ! -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 + ! 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 + ! 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 + ! 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 + ! 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 + ! 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 terms in the full distance matrix + double precision, allocatable, dimension(:) :: selfl21 - ! Pre-computed sine terms - double precision, allocatable, dimension(:,:,:) :: sin1 + ! 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) - double precision :: mol_dist + ! Value of PI at full FORTRAN precision. + double precision, parameter :: pi = 4.0d0*atan(1.0d0) + double precision :: mol_dist - r_width2 = r_width**2 - c_width2 = c_width**2 + r_width2 = r_width**2 + c_width2 = c_width**2 - inv_cut = pi / (2.0d0 * cut_distance) - inv_sigma2(:) = -1.0d0 / (sigmas(:))**2 + inv_cut = pi/(2.0d0*cut_distance) + inv_sigma2(:) = -1.0d0/(sigmas(:))**2 - allocate(sin1(nm1, maxval(n1), maxval(n1))) + 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 + !$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) + end do + end do + end do + !$OMP END PARALLEL DO - allocate(selfl21(nm1)) + allocate (selfl21(nm1)) - selfl21 = 0.0d0 + selfl21 = 0.0d0 - !$OMP PARALLEL DO PRIVATE(ni) REDUCTION(+:selfl21) - do i = 1, nm1 - ni = n1(i) - do i_1 = 1, ni - do j_1 = 1, ni + !$OMP PARALLEL DO PRIVATE(ni) REDUCTION(+:selfl21) + do i = 1, nm1 + ni = n1(i) + do i_1 = 1, ni + do j_1 = 1, ni - selfl21(i) = selfl21(i) + atomic_distl2(q1(i,i_1,:,:), q1(i,j_1,:,:), n1(i), n1(i), & - & sin1(i,i_1,:), sin1(i,j_1,:), width, cut_distance, r_width, c_width) & - & * (r_width2/(r_width2 + (z1(i,i_1,1) - z1(i,j_1,1))**2) * & - & c_width2/(c_width2 + (z1(i,i_1,2) - z1(i,j_1,2))**2)) + selfl21(i) = selfl21(i) + atomic_distl2(q1(i, i_1, :, :), q1(i, j_1, :, :), n1(i), n1(i), & + & sin1(i, i_1, :), sin1(i, j_1, :), width, cut_distance, r_width, c_width) & + & *(r_width2/(r_width2 + (z1(i, i_1, 1) - z1(i, j_1, 1))**2)* & + & c_width2/(c_width2 + (z1(i, i_1, 2) - z1(i, j_1, 2))**2)) - enddo - enddo - enddo - !$OMP END PARALLEL DO + end do + end do + end do + !$OMP END PARALLEL DO - allocate(atomic_distance(maxval(n1), maxval(n1))) + allocate (atomic_distance(maxval(n1), maxval(n1))) - kernels(:,:,:) = 0.0d0 - atomic_distance(:,:) = 0.0d0 + kernels(:, :, :) = 0.0d0 + atomic_distance(:, :) = 0.0d0 - !$OMP PARALLEL DO PRIVATE(l2dist,atomic_distance,ni,nj,mol_dist) schedule(dynamic) - do j = 1, nm1 - nj = n1(j) - do i = 1, j! nm1 + !$OMP PARALLEL DO PRIVATE(l2dist,atomic_distance,ni,nj,mol_dist) schedule(dynamic) + do j = 1, nm1 + nj = n1(j) + do i = 1, j! nm1 - ni = n1(i) + ni = n1(i) - atomic_distance(:,:) = 0.0d0 + atomic_distance(:, :) = 0.0d0 - do i_1 = 1, ni - do j_1 = 1, nj + 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 = 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 = 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)) + L2dist = 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 + atomic_distance(i_1, j_1) = l2dist - enddo - enddo + end do + end do - mol_dist = selfl21(i) + selfl21(j) - 2.0d0 * sum(atomic_distance(:ni,:nj)) + mol_dist = selfl21(i) + selfl21(j) - 2.0d0*sum(atomic_distance(:ni, :nj)) - do k = 1, nsigmas - kernels(k, i, j) = exp(mol_dist * inv_sigma2(k)) - kernels(k, j, i) = kernels(k, i, j) - enddo + do k = 1, nsigmas + kernels(k, i, j) = exp(mol_dist*inv_sigma2(k)) + kernels(k, j, i) = kernels(k, i, j) + end do - enddo - enddo - !$OMP END PARALLEL DO + end do + end do + !$OMP END PARALLEL DO - deallocate(atomic_distance) - deallocate(selfl21) - deallocate(sin1) + 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 + use arad, only: atomic_distl2 + + implicit none - 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 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 + ! 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 - ! -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 + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: n1 + integer, dimension(:), intent(in) :: n2 - ! Internal counters - integer :: i, j, k, ni, nj - integer :: m_1, i_1, j_1 + ! Sigma in the Gaussian kernel + double precision, dimension(:), intent(in) :: sigmas - ! Pre-computed constants - double precision :: r_width2 - double precision :: c_width2 - double precision :: inv_cut + ! Number of molecules + integer, intent(in) :: nm1 + integer, intent(in) :: nm2 - ! Temporary variables necessary for parallelization - double precision :: l2dist - double precision, allocatable, dimension(:,:) :: atomic_distance + ! Number of sigmas + integer, intent(in) :: nsigmas - ! Pre-computed terms in the full distance matrix - double precision, allocatable, dimension(:) :: selfl21 - double precision, allocatable, dimension(:) :: selfl22 + ! -1.0 / sigma^2 for use in the kernel + double precision, dimension(nsigmas) :: inv_sigma2 - ! Pre-computed sine terms - double precision, allocatable, dimension(:,:,:) :: sin1 - double precision, allocatable, dimension(:,:,:) :: sin2 + ! 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 - ! Value of PI at full FORTRAN precision. - double precision, parameter :: pi = 4.0d0 * atan(1.0d0) + ! Resulting alpha vector + double precision, dimension(nsigmas, nm1, nm2), intent(out) :: kernels - double precision :: mol_dist + ! Internal counters + integer :: i, j, k, ni, nj + integer :: m_1, i_1, j_1 - r_width2 = r_width**2 - c_width2 = c_width**2 + ! Pre-computed constants + double precision :: r_width2 + double precision :: c_width2 + double precision :: inv_cut - inv_cut = pi / (2.0d0 * cut_distance) - inv_sigma2(:) = -1.0d0 / (sigmas(:))**2 + ! Temporary variables necessary for parallelization + double precision :: l2dist + double precision, allocatable, dimension(:, :) :: atomic_distance - allocate(sin1(nm1, maxval(n1), maxval(n1))) - allocate(sin2(nm2, maxval(n2), maxval(n2))) + ! Pre-computed terms in the full distance matrix + double precision, allocatable, dimension(:) :: selfl21 + double precision, allocatable, dimension(:) :: selfl22 - sin1 = 0.0d0 - sin2 = 0.0d0 + ! Pre-computed sine terms + double precision, allocatable, dimension(:, :, :) :: sin1 + double precision, allocatable, dimension(:, :, :) :: sin2 - !$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 + ! Value of PI at full FORTRAN precision. + double precision, parameter :: pi = 4.0d0*atan(1.0d0) - !$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 + double precision :: mol_dist - allocate(selfl21(nm1)) - allocate(selfl22(nm2)) + r_width2 = r_width**2 + c_width2 = c_width**2 - selfl21 = 0.0d0 - selfl22 = 0.0d0 + inv_cut = pi/(2.0d0*cut_distance) + inv_sigma2(:) = -1.0d0/(sigmas(:))**2 - !$OMP PARALLEL DO PRIVATE(ni) REDUCTION(+:selfl21) - do i = 1, nm1 - ni = n1(i) - do i_1 = 1, ni - do j_1 = 1, ni + allocate (sin1(nm1, maxval(n1), maxval(n1))) + allocate (sin2(nm2, maxval(n2), maxval(n2))) - selfl21(i) = selfl21(i) + atomic_distl2(q1(i,i_1,:,:), q1(i,j_1,:,:), n1(i), n1(i), & - & sin1(i,i_1,:), sin1(i,j_1,:), width, cut_distance, r_width, c_width) & - & * (r_width2/(r_width2 + (z1(i,i_1,1) - z1(i,j_1,1))**2) * & - & c_width2/(c_width2 + (z1(i,i_1,2) - z1(i,j_1,2))**2)) + sin1 = 0.0d0 + sin2 = 0.0d0 - enddo + !$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) + end if + end do + end do + end do + !$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) + end if + end do + end do + end do + !$OMP END PARALLEL DO - enddo - enddo - !$OMP END PARALLEL DO + allocate (selfl21(nm1)) + allocate (selfl22(nm2)) - !$OMP PARALLEL DO PRIVATE(ni) REDUCTION(+:selfl22) - do i = 1, nm2 - ni = n2(i) - do i_1 = 1, ni - do j_1 = 1, ni + selfl21 = 0.0d0 + selfl22 = 0.0d0 - selfl22(i) = selfl22(i) + atomic_distl2(q2(i,i_1,:,:), q2(i,j_1,:,:), n2(i), n2(i), & - & sin2(i,i_1,:), sin2(i,j_1,:), width, cut_distance, r_width, c_width) & - &* (r_width2/(r_width2 + (z2(i,i_1,1) - z2(i,j_1,1))**2) * & - & c_width2/(c_width2 + (z2(i,i_1,2) - z2(i,j_1,2))**2)) + !$OMP PARALLEL DO PRIVATE(ni) REDUCTION(+:selfl21) + do i = 1, nm1 + ni = n1(i) + do i_1 = 1, ni + do j_1 = 1, ni + selfl21(i) = selfl21(i) + atomic_distl2(q1(i, i_1, :, :), q1(i, j_1, :, :), n1(i), n1(i), & + & sin1(i, i_1, :), sin1(i, j_1, :), width, cut_distance, r_width, c_width) & + & *(r_width2/(r_width2 + (z1(i, i_1, 1) - z1(i, j_1, 1))**2)* & + & c_width2/(c_width2 + (z1(i, i_1, 2) - z1(i, j_1, 2))**2)) - enddo + end do - enddo - enddo - !$OMP END PARALLEL DO + end do + end do + !$OMP END PARALLEL DO + !$OMP PARALLEL DO PRIVATE(ni) REDUCTION(+:selfl22) + do i = 1, nm2 + ni = n2(i) + do i_1 = 1, ni + do j_1 = 1, ni - allocate(atomic_distance(maxval(n1), maxval(n2))) + selfl22(i) = selfl22(i) + atomic_distl2(q2(i, i_1, :, :), q2(i, j_1, :, :), n2(i), n2(i), & + & sin2(i, i_1, :), sin2(i, j_1, :), width, cut_distance, r_width, c_width) & + &*(r_width2/(r_width2 + (z2(i, i_1, 1) - z2(i, j_1, 1))**2)* & + & c_width2/(c_width2 + (z2(i, i_1, 2) - z2(i, j_1, 2))**2)) - kernels(:,:,:) = 0.0d0 - atomic_distance(:,:) = 0.0d0 + end do + end do + end do + !$OMP END PARALLEL DO - !$OMP PARALLEL DO PRIVATE(l2dist,atomic_distance,ni,nj,mol_dist) schedule(dynamic) + allocate (atomic_distance(maxval(n1), maxval(n2))) - do j = 1, nm2 - nj = n2(j) - do i = 1, nm1 - ni = n1(i) + kernels(:, :, :) = 0.0d0 + atomic_distance(:, :) = 0.0d0 - atomic_distance(:,:) = 0.0d0 + !$OMP PARALLEL DO PRIVATE(l2dist,atomic_distance,ni,nj,mol_dist) schedule(dynamic) - do i_1 = 1, ni - do j_1 = 1, nj + do j = 1, nm2 + nj = n2(j) + do i = 1, nm1 + ni = n1(i) - 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) + atomic_distance(:, :) = 0.0d0 + do i_1 = 1, ni + do j_1 = 1, nj - L2dist = 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)) + 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 = 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 + atomic_distance(i_1, j_1) = l2dist - enddo - enddo + end do + end do - mol_dist = selfl21(i) + selfl22(j) - 2.0d0 * sum(atomic_distance(:ni,:nj)) + mol_dist = selfl21(i) + selfl22(j) - 2.0d0*sum(atomic_distance(:ni, :nj)) - do k = 1, nsigmas - kernels(k, i, j) = exp(mol_dist * inv_sigma2(k)) + do k = 1, nsigmas + kernels(k, i, j) = exp(mol_dist*inv_sigma2(k)) - enddo + end do - enddo - enddo - !$OMP END PARALLEL DO + end do + end do + !$OMP END PARALLEL DO - deallocate(atomic_distance) - deallocate(selfl21) - deallocate(selfl22) - deallocate(sin1) - deallocate(sin2) + deallocate (atomic_distance) + deallocate (selfl21) + deallocate (selfl22) + deallocate (sin1) + deallocate (sin2) end subroutine fget_global_kernels_arad diff --git a/src/qmllib/representations/bob/__init__.py b/src/qmllib/representations/bob/__init__.py index 8f3abfc..38d46fc 100644 --- a/src/qmllib/representations/bob/__init__.py +++ b/src/qmllib/representations/bob/__init__.py @@ -2,7 +2,10 @@ Bag of bonds utils functions """ +from typing import List + import numpy as np +from numpy import ndarray from qmllib.constants.periodic_table import ELEMENT_NAME @@ -23,7 +26,7 @@ def get_natypes(nuclear_charges: np.ndarray) -> dict[str, int]: return natypes -def get_asize(list_nuclear_charges, pad) -> dict[str, int]: +def get_asize(list_nuclear_charges: List[ndarray], pad: int) -> dict[str, int]: """ example: diff --git a/src/qmllib/representations/fchl/fchl_kernel_functions.py b/src/qmllib/representations/fchl/fchl_kernel_functions.py index 785cf36..1af3cc7 100644 --- a/src/qmllib/representations/fchl/fchl_kernel_functions.py +++ b/src/qmllib/representations/fchl/fchl_kernel_functions.py @@ -1,12 +1,17 @@ from __future__ import division, print_function +from typing import Dict, List, Optional, Tuple, Union + import numpy as np +from numpy import ndarray from scipy.special import binom, factorial from .ffchl_module import ffchl_kernel_types as kt -def get_gaussian_parameters(tags): +def get_gaussian_parameters( + tags: Optional[Dict[str, List[float]]] +) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = { @@ -25,7 +30,7 @@ def get_gaussian_parameters(tags): return kt.gaussian, parameters, n_kernels -def get_linear_parameters(tags): +def get_linear_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = { @@ -41,7 +46,7 @@ def get_linear_parameters(tags): return kt.linear, parameters, n_kernels -def get_polynomial_parameters(tags): +def get_polynomial_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = {"alpha": [1.0], "c": [0.0], "d": [1.0]} @@ -54,7 +59,7 @@ def get_polynomial_parameters(tags): return kt.polynomial, parameters, n_kernels -def get_sigmoid_parameters(tags): +def get_sigmoid_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = { @@ -74,7 +79,7 @@ def get_sigmoid_parameters(tags): return kt.sigmoid, parameters, n_kernels -def get_multiquadratic_parameters(tags): +def get_multiquadratic_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = { @@ -93,7 +98,9 @@ def get_multiquadratic_parameters(tags): return kt.multiquadratic, parameters, n_kernels -def get_inverse_multiquadratic_parameters(tags): +def get_inverse_multiquadratic_parameters( + tags: Dict[str, List[float]] +) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = { @@ -112,7 +119,7 @@ def get_inverse_multiquadratic_parameters(tags): return kt.inv_multiquadratic, parameters, n_kernels -def get_bessel_parameters(tags): +def get_bessel_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = {"sigma": [1.0], "v": [1.0], "n": [1.0]} @@ -126,7 +133,7 @@ def get_bessel_parameters(tags): return kt.bessel, parameters, n_kernels -def get_l2_parameters(tags): +def get_l2_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = { @@ -146,7 +153,7 @@ def get_l2_parameters(tags): return kt.l2, parameters, n_kernels -def get_matern_parameters(tags): +def get_matern_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = { @@ -174,7 +181,7 @@ def get_matern_parameters(tags): return kt.matern, parameters, n_kernels -def get_cauchy_parameters(tags): +def get_cauchy_parameters(tags: Dict[str, List[float]]) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = { @@ -193,7 +200,7 @@ def get_cauchy_parameters(tags): return kt.cauchy, parameters, n_kernels -def get_polynomial2_parameters(tags): +def get_polynomial2_parameters(tags: Dict[str, List[List[float]]]) -> Tuple[ndarray, ndarray, int]: if tags is None: tags = { @@ -211,7 +218,9 @@ def get_polynomial2_parameters(tags): return kt.polynomial2, parameters, n_kernels -def get_kernel_parameters(name, tags): +def get_kernel_parameters( + name: str, tags: Optional[Union[Dict[str, List[List[float]]], Dict[str, List[float]]]] +) -> Tuple[ndarray, ndarray, int]: parameters = None idx = kt.gaussian diff --git a/src/qmllib/representations/fchl/fchl_representations.py b/src/qmllib/representations/fchl/fchl_representations.py index 6b21e59..f8d9588 100644 --- a/src/qmllib/representations/fchl/fchl_representations.py +++ b/src/qmllib/representations/fchl/fchl_representations.py @@ -1,11 +1,18 @@ import copy +from typing import List, Optional, Union import numpy as np +from numpy import ndarray def generate_representation( - coordinates, nuclear_charges, max_size=23, neighbors=23, cut_distance=5.0, cell=None -): + coordinates: Union[ndarray, List[List[float]]], + nuclear_charges: ndarray, + max_size: int = 23, + neighbors: int = 23, + cut_distance: float = 5.0, + cell: Optional[ndarray] = None, +) -> ndarray: """Generates a representation for the FCHL kernel module. :param coordinates: Input coordinates. @@ -180,13 +187,13 @@ def generate_displaced_representations_5point( def generate_representation_electric_field( - coordinates, - nuclear_charges, - fictitious_charges="gasteiger", - max_size=23, - neighbors=23, - cut_distance=5.0, -): + coordinates: ndarray, + nuclear_charges: ndarray, + fictitious_charges: Union[ndarray, List[float]] = "gasteiger", + max_size: int = 23, + neighbors: int = 23, + cut_distance: float = 5.0, +) -> ndarray: """Generates a representation for the FCHL kernel module, including fictitious partial charges. For use with fist-order electric field-dependent properties, such as dipole moments and IR intensity. diff --git a/src/qmllib/representations/fchl/fchl_scalar_kernels.py b/src/qmllib/representations/fchl/fchl_scalar_kernels.py index 58ecc06..9ae0511 100644 --- a/src/qmllib/representations/fchl/fchl_scalar_kernels.py +++ b/src/qmllib/representations/fchl/fchl_scalar_kernels.py @@ -1,6 +1,9 @@ from __future__ import absolute_import +from typing import Dict, List, Optional, Union + import numpy as np +from numpy import float64, ndarray from qmllib.utils.alchemy import get_alchemy @@ -17,24 +20,24 @@ def get_local_kernels( - A, - B, - verbose=False, - two_body_scaling=np.sqrt(8), - three_body_scaling=1.6, - two_body_width=0.2, - three_body_width=np.pi, - two_body_power=4.0, - three_body_power=2.0, - cut_start=1.0, - cut_distance=5.0, - fourier_order=1, - alchemy="periodic-table", - alchemy_period_width=1.6, - alchemy_group_width=1.6, - kernel="gaussian", - kernel_args=None, -): + A: ndarray, + B: ndarray, + verbose: bool = False, + two_body_scaling: float = np.sqrt(8), + three_body_scaling: float = 1.6, + two_body_width: float = 0.2, + three_body_width: float = np.pi, + two_body_power: float = 4.0, + three_body_power: float = 2.0, + cut_start: float = 1.0, + cut_distance: float = 5.0, + fourier_order: int = 1, + alchemy: str = "periodic-table", + alchemy_period_width: float = 1.6, + alchemy_group_width: float = 1.6, + kernel: str = "gaussian", + kernel_args: Optional[Dict[str, List[float]]] = None, +) -> ndarray: """Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\\sigma^2} \\big)` @@ -147,23 +150,23 @@ def get_local_kernels( def get_local_symmetric_kernels( - A, - verbose=False, - two_body_scaling=np.sqrt(8), - three_body_scaling=1.6, - two_body_width=0.2, - three_body_width=np.pi, - two_body_power=4.0, - three_body_power=2.0, - cut_start=1.0, - cut_distance=5.0, - fourier_order=1, - alchemy="periodic-table", - alchemy_period_width=1.6, - alchemy_group_width=1.6, - kernel="gaussian", - kernel_args=None, -): + A: ndarray, + verbose: bool = False, + two_body_scaling: Union[float, float64] = np.sqrt(8), + three_body_scaling: float = 1.6, + two_body_width: float = 0.2, + three_body_width: float = np.pi, + two_body_power: float = 4.0, + three_body_power: float = 2.0, + cut_start: float = 1.0, + cut_distance: float = 5.0, + fourier_order: int = 1, + alchemy: Union[ndarray, str] = "periodic-table", + alchemy_period_width: float = 1.6, + alchemy_group_width: float = 1.6, + kernel: str = "gaussian", + kernel_args: Optional[Union[Dict[str, List[List[float]]], Dict[str, List[float]]]] = None, +) -> ndarray: """Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - A_j\\|_2^2}{2\\sigma^2} \\big)` @@ -254,23 +257,23 @@ def get_local_symmetric_kernels( def get_global_symmetric_kernels( - A, - verbose=False, - two_body_scaling=np.sqrt(8), - three_body_scaling=1.6, - two_body_width=0.2, - three_body_width=np.pi, - two_body_power=4.0, - three_body_power=2.0, - cut_start=1.0, - cut_distance=5.0, - fourier_order=1, - alchemy="periodic-table", - alchemy_period_width=1.6, - alchemy_group_width=1.6, - kernel="gaussian", - kernel_args=None, -): + A: ndarray, + verbose: bool = False, + two_body_scaling: float64 = np.sqrt(8), + three_body_scaling: float = 1.6, + two_body_width: float = 0.2, + three_body_width: float = np.pi, + two_body_power: float = 4.0, + three_body_power: float = 2.0, + cut_start: float = 1.0, + cut_distance: float = 5.0, + fourier_order: int = 1, + alchemy: str = "periodic-table", + alchemy_period_width: float = 1.6, + alchemy_group_width: float = 1.6, + kernel: str = "gaussian", + kernel_args: Optional[Dict[str, List[float]]] = None, +) -> ndarray: """Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - A_j\\|_2^2}{2\\sigma^2} \\big)` @@ -361,24 +364,24 @@ def get_global_symmetric_kernels( def get_global_kernels( - A, - B, - verbose=False, - two_body_scaling=np.sqrt(8), - three_body_scaling=1.6, - two_body_width=0.2, - three_body_width=np.pi, - two_body_power=4.0, - three_body_power=2.0, - cut_start=1.0, - cut_distance=5.0, - fourier_order=1, - alchemy="periodic-table", - alchemy_period_width=1.6, - alchemy_group_width=1.6, - kernel="gaussian", - kernel_args=None, -): + A: ndarray, + B: ndarray, + verbose: bool = False, + two_body_scaling: float64 = np.sqrt(8), + three_body_scaling: float = 1.6, + two_body_width: float = 0.2, + three_body_width: float = np.pi, + two_body_power: float = 4.0, + three_body_power: float = 2.0, + cut_start: float = 1.0, + cut_distance: float = 5.0, + fourier_order: int = 1, + alchemy: str = "periodic-table", + alchemy_period_width: float = 1.6, + alchemy_group_width: float = 1.6, + kernel: str = "gaussian", + kernel_args: Optional[Dict[str, List[float]]] = None, +) -> ndarray: """Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\\sigma^2} \\big)` @@ -490,24 +493,24 @@ def get_global_kernels( def get_atomic_kernels( - A, - B, - verbose=False, - two_body_scaling=np.sqrt(8), - three_body_scaling=1.6, - two_body_width=0.2, - three_body_width=np.pi, - two_body_power=4.0, - three_body_power=2.0, - cut_start=1.0, - cut_distance=5.0, - fourier_order=1, - alchemy="periodic-table", - alchemy_period_width=1.6, - alchemy_group_width=1.6, - kernel="gaussian", - kernel_args=None, -): + A: ndarray, + B: ndarray, + verbose: bool = False, + two_body_scaling: float64 = np.sqrt(8), + three_body_scaling: float = 1.6, + two_body_width: float = 0.2, + three_body_width: float = np.pi, + two_body_power: float = 4.0, + three_body_power: float = 2.0, + cut_start: float = 1.0, + cut_distance: float = 5.0, + fourier_order: int = 1, + alchemy: str = "periodic-table", + alchemy_period_width: float = 1.6, + alchemy_group_width: float = 1.6, + kernel: str = "gaussian", + kernel_args: Optional[Dict[str, List[float]]] = None, +) -> ndarray: """Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\\sigma^2} \\big)` @@ -606,23 +609,23 @@ def get_atomic_kernels( def get_atomic_symmetric_kernels( - A, - verbose=False, - two_body_scaling=np.sqrt(8), - three_body_scaling=1.6, - two_body_width=0.2, - three_body_width=np.pi, - two_body_power=4.0, - three_body_power=2.0, - cut_start=1.0, - cut_distance=5.0, - fourier_order=1, - alchemy="periodic-table", - alchemy_period_width=1.6, - alchemy_group_width=1.6, - kernel="gaussian", - kernel_args=None, -): + A: ndarray, + verbose: bool = False, + two_body_scaling: float64 = np.sqrt(8), + three_body_scaling: float = 1.6, + two_body_width: float = 0.2, + three_body_width: float = np.pi, + two_body_power: float = 4.0, + three_body_power: float = 2.0, + cut_start: float = 1.0, + cut_distance: float = 5.0, + fourier_order: int = 1, + alchemy: str = "periodic-table", + alchemy_period_width: float = 1.6, + alchemy_group_width: float = 1.6, + kernel: str = "gaussian", + kernel_args: Optional[Dict[str, List[float]]] = None, +) -> ndarray: """Calculates the Gaussian kernel matrix K, where :math:`K_{ij}`: :math:`K_{ij} = \\exp \\big( -\\frac{\\|A_i - B_j\\|_2^2}{2\\sigma^2} \\big)` diff --git a/src/qmllib/representations/representations.py b/src/qmllib/representations/representations.py index fa47e4b..78bf8b3 100644 --- a/src/qmllib/representations/representations.py +++ b/src/qmllib/representations/representations.py @@ -1,6 +1,8 @@ import itertools +from typing import Dict, List, Optional, Tuple, Union import numpy as np +from numpy import int64, ndarray from qmllib.constants.periodic_table import NUCLEAR_CHARGE @@ -50,7 +52,9 @@ def vector_to_matrix(v): return M -def generate_coulomb_matrix(nuclear_charges, coordinates, size=23, sorting="row-norm"): +def generate_coulomb_matrix( + nuclear_charges: ndarray, coordinates: ndarray, size: int = 23, sorting: str = "row-norm" +) -> ndarray: """ Creates a Coulomb Matrix representation of a molecule. Sorting of the elements can either be done by ``sorting="row-norm"`` or ``sorting="unsorted"``. A matrix :math:`M` is constructed with elements @@ -102,16 +106,16 @@ def generate_coulomb_matrix(nuclear_charges, coordinates, size=23, sorting="row- def generate_atomic_coulomb_matrix( - nuclear_charges, - coordinates, - size=23, - sorting="distance", - central_cutoff=1e6, - central_decay=-1, - interaction_cutoff=1e6, - interaction_decay=-1, - indices=None, -): + nuclear_charges: ndarray, + coordinates: ndarray, + size: int = 23, + sorting: str = "distance", + central_cutoff: float = 1e6, + central_decay: Union[float, int] = -1, + interaction_cutoff: float = 1e6, + interaction_decay: Union[float, int] = -1, + indices: Optional[List[int]] = None, +) -> ndarray: """ Creates a Coulomb Matrix representation of the local environment of a central atom. For each central atom :math:`k`, a matrix :math:`M` is constructed with elements @@ -240,7 +244,9 @@ def generate_atomic_coulomb_matrix( raise SystemExit -def generate_eigenvalue_coulomb_matrix(nuclear_charges, coordinates, size=23): +def generate_eigenvalue_coulomb_matrix( + nuclear_charges: ndarray, coordinates: ndarray, size: int = 23 +) -> ndarray: """ Creates an eigenvalue Coulomb Matrix representation of a molecule. A matrix :math:`M` is constructed with elements @@ -271,12 +277,12 @@ def generate_eigenvalue_coulomb_matrix(nuclear_charges, coordinates, size=23): def generate_bob( - nuclear_charges, - coordinates, - atomtypes, - size=23, - asize={"O": 3, "C": 7, "N": 3, "H": 16, "S": 1}, -): + nuclear_charges: ndarray, + coordinates: ndarray, + atomtypes: ndarray, + size: int = 23, + asize: Dict[str, Union[int64, int]] = {"O": 3, "C": 7, "N": 3, "H": 16, "S": 1}, +) -> ndarray: """Creates a Bag of Bonds (BOB) representation of a molecule. The representation expands on the coulomb matrix representation. For each element a bag (vector) is constructed for self interactions @@ -327,7 +333,7 @@ def generate_bob( return fgenerate_bob(nuclear_charges, coordinates, nuclear_charges, ids, nmax, n) -def get_slatm_mbtypes(nuclear_charges, pbc="000"): +def get_slatm_mbtypes(nuclear_charges: List[ndarray], pbc: str = "000") -> List[List[int64]]: """ Get the list of minimal types of many-body terms in a dataset. This resulting list is necessary as input in the ``generate_slatm_representation()`` function. @@ -395,18 +401,18 @@ def get_slatm_mbtypes(nuclear_charges, pbc="000"): def generate_slatm( - coordinates, - nuclear_charges, - mbtypes, - unit_cell=None, - local=False, - sigmas=[0.05, 0.05], - dgrids=[0.03, 0.03], - rcut=4.8, - alchemy=False, - pbc="000", - rpower=6, -): + coordinates: ndarray, + nuclear_charges: ndarray, + mbtypes: List[List[int64]], + unit_cell: None = None, + local: bool = False, + sigmas: List[float] = [0.05, 0.05], + dgrids: List[float] = [0.03, 0.03], + rcut: float = 4.8, + alchemy: bool = False, + pbc: str = "000", + rpower: int = 6, +) -> Union[ndarray, List[ndarray]]: """ Generate Spectrum of London and Axillrod-Teller-Muto potential (SLATM) representation. Both global (``local=False``) and local (``local=True``) SLATM are available. @@ -613,21 +619,21 @@ def generate_slatm( def generate_acsf( - nuclear_charges, - coordinates, - elements=[1, 6, 7, 8, 16], - nRs2=3, - nRs3=3, - nTs=3, - eta2=1, - eta3=1, - zeta=1, - rcut=5, - acut=5, - bin_min=0.8, - gradients=False, - pad=None, -): + nuclear_charges: List[int], + coordinates: ndarray, + elements: List[int] = [1, 6, 7, 8, 16], + nRs2: int = 3, + nRs3: int = 3, + nTs: int = 3, + eta2: int = 1, + eta3: int = 1, + zeta: int = 1, + rcut: int = 5, + acut: int = 5, + bin_min: float = 0.8, + gradients: bool = False, + pad: Optional[int] = None, +) -> Union[Tuple[ndarray, ndarray], ndarray]: """ Generate the variant of atom-centered symmetry functions used in https://doi.org/10.1039/C7SC04934J @@ -730,23 +736,23 @@ def generate_acsf( def generate_fchl_acsf( - nuclear_charges, - coordinates, - elements=[1, 6, 7, 8, 16], - nRs2=24, - nRs3=20, - nFourier=1, - eta2=0.32, - eta3=2.7, - zeta=np.pi, - rcut=8.0, - acut=8.0, - two_body_decay=1.8, - three_body_decay=0.57, - three_body_weight=13.4, - pad=False, - gradients=False, -): + nuclear_charges: ndarray, + coordinates: ndarray, + elements: List[int] = [1, 6, 7, 8, 16], + nRs2: int = 24, + nRs3: int = 20, + nFourier: int = 1, + eta2: float = 0.32, + eta3: float = 2.7, + zeta: float = np.pi, + rcut: float = 8.0, + acut: float = 8.0, + two_body_decay: float = 1.8, + three_body_decay: float = 0.57, + three_body_weight: float = 13.4, + pad: Union[int, bool] = False, + gradients: bool = False, +) -> Union[Tuple[ndarray, ndarray], ndarray]: """ FCHL-ACSF diff --git a/src/qmllib/representations/slatm.py b/src/qmllib/representations/slatm.py index c40ecb6..7de31c3 100644 --- a/src/qmllib/representations/slatm.py +++ b/src/qmllib/representations/slatm.py @@ -1,4 +1,7 @@ +from typing import List, Optional + import numpy as np +from numpy import int64, ndarray from .fslatm import fget_sbop, fget_sbop_local, fget_sbot, fget_sbot_local @@ -97,7 +100,7 @@ def update_m(obj, ia, rcut=9.0, pbc=None): # return obj_u -def get_boa(z1, zs_): +def get_boa(z1: int64, zs_: ndarray) -> ndarray: return z1 * np.array( [ (zs_ == z1).sum(), @@ -107,17 +110,17 @@ def get_boa(z1, zs_): def get_sbop( - mbtype, - obj, - iloc=False, - ia=None, - normalize=True, - sigma=0.05, - rcut=4.8, - dgrid=0.03, - pbc="000", - rpower=6, -): + mbtype: List[int64], + obj: List[ndarray], + iloc: bool = False, + ia: Optional[int] = None, + normalize: bool = True, + sigma: float = 0.05, + rcut: float = 4.8, + dgrid: float = 0.03, + pbc: str = "000", + rpower: int = 6, +) -> ndarray: """ two-body terms @@ -156,8 +159,16 @@ def get_sbop( def get_sbot( - mbtype, obj, iloc=False, ia=None, normalize=True, sigma=0.05, rcut=4.8, dgrid=0.0262, pbc="000" -): + mbtype: List[int64], + obj: List[ndarray], + iloc: bool = False, + ia: Optional[int] = None, + normalize: bool = True, + sigma: float = 0.05, + rcut: float = 4.8, + dgrid: float = 0.0262, + pbc: str = "000", +) -> ndarray: """ sigma -- standard deviation of gaussian distribution centered on a specific angle diff --git a/src/qmllib/solvers/__init__.py b/src/qmllib/solvers/__init__.py index 3f9bf1d..eab39b5 100644 --- a/src/qmllib/solvers/__init__.py +++ b/src/qmllib/solvers/__init__.py @@ -1,3 +1,5 @@ +from typing import Optional + import numpy as np from numpy import ndarray @@ -158,7 +160,7 @@ def bkf_solve(A: ndarray, y: ndarray) -> ndarray: return x -def svd_solve(A, y, rcond=None): +def svd_solve(A: ndarray, y: ndarray, rcond: Optional[float] = None) -> ndarray: """Solves the equation :math:`A x = y` diff --git a/src/qmllib/utils/alchemy.py b/src/qmllib/utils/alchemy.py index 005c02c..fad9610 100644 --- a/src/qmllib/utils/alchemy.py +++ b/src/qmllib/utils/alchemy.py @@ -1,6 +1,8 @@ from copy import copy +from typing import Any, Dict, Tuple, Union import numpy as np +from numpy import float64, ndarray # Periodic table indexes PTP = { @@ -260,16 +262,16 @@ def get_alchemy( - alchemy, - emax=100, - r_width=0.001, - c_width=0.001, - elemental_vectors={}, - n_width=0.001, - m_width=0.001, - l_width=0.001, - s_width=0.001, -): + alchemy: Union[ndarray, str], + emax: int = 100, + r_width: float = 0.001, + c_width: float = 0.001, + elemental_vectors: Dict[Any, Any] = {}, + n_width: float = 0.001, + m_width: float = 0.001, + l_width: float = 0.001, + s_width: float = 0.001, +) -> Tuple[bool, ndarray]: if isinstance(alchemy, np.ndarray): @@ -353,7 +355,7 @@ def gen_QNum_distances(emax=100, n_width=0.001, m_width=0.001, l_width=0.001, s_ return pd -def periodic_distance(a, b, r_width, c_width): +def periodic_distance(a: int, b: int, r_width: float, c_width: float) -> float64: """Calculate stochiometric distance a -- nuclear charge of element a @@ -372,7 +374,7 @@ def periodic_distance(a, b, r_width, c_width): return np.exp(-((ra - rb) ** 2) / (4 * r_width**2) - (ca - cb) ** 2 / (4 * c_width**2)) -def gen_pd(emax=100, r_width=0.001, c_width=0.001): +def gen_pd(emax: int = 100, r_width: float = 0.001, c_width: float = 0.001) -> ndarray: """Generate stochiometric ditance matrix emax -- Largest element diff --git a/src/qmllib/utils/xyz_format.py b/src/qmllib/utils/xyz_format.py index d722098..3c46d92 100644 --- a/src/qmllib/utils/xyz_format.py +++ b/src/qmllib/utils/xyz_format.py @@ -1,11 +1,13 @@ from pathlib import Path +from typing import Tuple import numpy as np +from numpy import ndarray from qmllib.constants.periodic_table import NUCLEAR_CHARGE -def read_xyz(filename: str | Path): +def read_xyz(filename: str | Path) -> Tuple[ndarray, ndarray]: """(Re-)initializes the Compound-object with data from an xyz-file. :param filename: Input xyz-filename or file-like obejct