diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 93c5b499..5b1f1f46 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -3,7 +3,7 @@ repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v3.4.0 + rev: v5.0.0 hooks: - id: trailing-whitespace exclude: ^tests/resources/ @@ -23,7 +23,7 @@ repos: args: ['--maxkb=3000'] - repo: https://github.com/myint/autoflake - rev: v2.2.0 + rev: v2.3.1 hooks: - id: autoflake name: Removes unused variables @@ -48,7 +48,7 @@ repos: ] - repo: https://github.com/psf/black - rev: 22.3.0 + rev: 24.10.0 hooks: - id: black name: Fixes formatting @@ -56,7 +56,7 @@ repos: args: ["--line-length=99"] - repo: https://github.com/pycqa/flake8 - rev: 7.0.0 + rev: 7.1.1 hooks: - id: flake8 name: Checks pep8 style @@ -71,7 +71,7 @@ repos: "--ignore=E501,E203,W503,E741", ] - - repo: https://github.com/pseewald/fprettify + - repo: https://github.com/fortran-lang/fprettify rev: v0.3.7 hooks: - id: fprettify diff --git a/Makefile b/Makefile index dbbe9eac..27b0ce3b 100644 --- a/Makefile +++ b/Makefile @@ -31,7 +31,7 @@ test: types: ${python} -m monkeytype run $$(which ${pytest}) ./tests - ${python} -m monkeytype list-modules | grep ${pkg} | parallel -j${j} "${python} -m monkeytype apply {}" + ${python} -m monkeytype list-modules | grep ${pkg} | parallel -j${j} "${python} -m monkeytype apply {} > /dev/null && echo {}" cov: ${python} -m pytest -vrs --cov=${pkg} --cov-report html tests diff --git a/README.rst b/README.rst index 04930da5..5dd5bd81 100644 --- a/README.rst +++ b/README.rst @@ -23,11 +23,18 @@ Breaking changes from ``qml``: * FCHL representations callable interface to be consistent with other representations (e.i. atoms, coordinates) - ============== How to install ============== +You need a fortran compiler and math library. Default is `gfortran` and `openblas`. + + +.. code-block:: bash + + sudo apt install libopenblas-dev gcc + + A proper pip-package is on the way, for now .. code-block:: bash diff --git a/_compile.py b/_compile.py index 01256fb2..b2b4f92e 100644 --- a/_compile.py +++ b/_compile.py @@ -5,67 +5,85 @@ import sys from pathlib import Path +import numpy as np + f90_modules = { "representations/frepresentations": ["frepresentations.f90"], "representations/facsf": ["facsf.f90"], "representations/fslatm": ["fslatm.f90"], "representations/arad/farad_kernels": ["farad_kernels.f90"], "representations/fchl/ffchl_module": [ - "ffchl_module.f90", - "ffchl_scalar_kernels.f90", "ffchl_kernel_types.f90", + "ffchl_module.f90", + "ffchl_module_ef.f90", "ffchl_kernels.f90", - "ffchl_electric_field_kernels.f90", + "ffchl_scalar_kernels.f90", + "ffchl_kernels_ef.f90", "ffchl_force_kernels.f90", ], "solvers/fsolvers": ["fsolvers.f90"], "kernels/fdistance": ["fdistance.f90"], - "kernels/fkernels": ["fkernels.f90", "fkpca.f90"], + "kernels/fkernels": [ + "fkernels.f90", + "fkpca.f90", + "fkwasserstein.f90", + ], "kernels/fgradient_kernels": ["fgradient_kernels.f90"], + "utils/fsettings": ["fsettings.f90"], } def find_mkl(): - - return + raise NotImplementedError() -def find_flags(fcc: str): +def find_env() -> dict[str, str]: """Find compiler flags""" + # TODO Find ifort flags, choose from FCC # TODO Find math lib # TODO Find os - # -lgomp", "-lpthread", "-lm", "-ldl - # ["-L${MKLROOT}/lib/intel64", "-lmkl_rt"] - - # COMPILER_FLAGS = ["-O3", "-fopenmp", "-m64", "-march=native", "-fPIC", - # "-Wno-maybe-uninitialized", "-Wno-unused-function", "-Wno-cpp"] - # LINKER_FLAGS = ["-lgomp"] + COMPILER_FLAGS = [ + "-O3", + "-fopenmp", + "-m64", + "-march=native", + "-fPIC", + "-Wno-maybe-uninitialized", + "-Wno-unused-function", + "-Wno-cpp", + ] extra_flags = ["-lgomp", "-lpthread", "-lm", "-ldl"] + math_flags = ["-L/usr/lib/", "-lblas", "-llapack"] - flags = ["-L/usr/lib/", "-lblas", "-llapack"] + extra_flags - - return flags - - -def find_fcc(): - """Find the fortran compiler. Either gnu or intel""" - - # fcc = "ifort" + fflags = [] + COMPILER_FLAGS + ldflags = [] + extra_flags + math_flags fcc = "gfortran" - return fcc + env = {"FFLAGS": " ".join(fflags), "LDFLAGS": " ".join(ldflags), "FCC": fcc} + + return env def main(): """Compile f90 in src/qmllib""" - fcc = find_fcc() - flags = find_flags(fcc) + print(f"Using numpy {np.__version__}") + + # Find and set Fortran compiler, compiler flags and linker flags + env = find_env() + for key, value in env.items(): + print(f"export {key}='{value}'") + os.environ[key] = value + + f2py = [sys.executable, "-m", "numpy.f2py"] - os.environ["FCC"] = fcc + meson_flags = [ + "--backend", + "meson", + ] for module_name, module_sources in f90_modules.items(): @@ -74,9 +92,7 @@ def main(): stem = path.stem cwd = Path("src/qmllib") / parent - cmd = ( - [sys.executable, "-m", "numpy.f2py", "-c"] + flags + module_sources + ["-m", str(stem)] - ) + cmd = f2py + ["-c"] + module_sources + ["-m", str(stem)] + meson_flags print(cwd, " ".join(cmd)) proc = subprocess.run(cmd, cwd=cwd, capture_output=True, text=True) @@ -85,9 +101,9 @@ def main(): exitcode = proc.returncode if exitcode > 0: - print(stdout) - print() print(stderr) + print() + print(stdout) exit(exitcode) diff --git a/environment_dev.yaml b/environment_dev.yaml index bd3e993a..26f71544 100644 --- a/environment_dev.yaml +++ b/environment_dev.yaml @@ -4,19 +4,22 @@ channels: - defaults dependencies: - python==3.12 - - jupytext - - monkeytype - - numpy<2.0.0 - - pandas - pip - - pre-commit - - pyarrow - - pytest - - scikit-learn - - scipy - # build - - build - - meson - - ninja - # publish - - twine + - pip: + - jupytext + - monkeytype + - numpy<2.0.0 + - pandas + - pandas + - pip + - pre-commit + - pyarrow + - pytest + - scikit-learn + - scipy + # build + - build + - meson + - ninja + # publish + - twine diff --git a/src/qmllib/kernels/fkernels.f90 b/src/qmllib/kernels/fkernels.f90 index d3a9d4ab..44e8dc39 100644 --- a/src/qmllib/kernels/fkernels.f90 +++ b/src/qmllib/kernels/fkernels.f90 @@ -1,83 +1,4 @@ -module searchtools - - implicit none - -contains - - function searchsorted(all_values, sorted) result(cdf_idx) - - implicit none - - double precision, dimension(:), intent(in) :: all_values - double precision, dimension(:), intent(in) :: sorted - - integer, allocatable, dimension(:) :: cdf_idx - - double precision :: val - - integer :: i, j, n, m - - n = size(all_values) - 1 - m = size(sorted) - - allocate (cdf_idx(n)) - - cdf_idx(:) = 0 - - do i = 1, n - - val = all_values(i) - - do j = 1, m - - !write (*,*) i, j, sorted(j), val - - ! if ((sorted(j) <= val) .and. (val < sorted(j+1))) then - if (sorted(j) > val) then - - cdf_idx(i) = j - 1 - !write(*,*) "found" - exit - - ! endif - else !iif (val > maxval(sorted)) then - cdf_idx(i) = m - end if - - end do - - end do - - end function searchsorted - - recursive subroutine quicksort(a, first, last) - implicit none - double precision :: a(*), x, t - integer first, last - integer i, j - - x = a((first + last)/2) - i = first - j = last - do - do while (a(i) < x) - i = i + 1 - end do - do while (x < a(j)) - j = j - 1 - end do - if (i >= j) exit - t = a(i); a(i) = a(j); a(j) = t - i = i + 1 - j = j - 1 - end do - if (first < i - 1) call quicksort(a, first, i - 1) - if (j + 1 < last) call quicksort(a, j + 1, last) - end subroutine quicksort - -end module searchtools - subroutine fget_local_kernels_gaussian(q1, q2, n1, n2, sigmas, & & nm1, nm2, nsigmas, kernels) @@ -796,95 +717,3 @@ subroutine fsargan_kernel(a, na, b, nb, k, sigma, gammas, ng) deallocate (prefactor) end subroutine fsargan_kernel - -subroutine fwasserstein_kernel(a, na, b, nb, k, sigma, p, q) - - use searchtools - implicit none - - double precision, dimension(:, :), intent(in) :: a - double precision, dimension(:, :), intent(in) :: b - - double precision, allocatable, dimension(:, :) :: asorted - double precision, allocatable, dimension(:, :) :: bsorted - - double precision, allocatable, dimension(:) :: rep - - integer, intent(in) :: na, nb - - double precision, dimension(:, :), intent(inout) :: k - double precision, intent(in) :: sigma - - integer, intent(in) :: p - integer, intent(in) :: q - - double precision :: inv_sigma - - integer :: i, j, l - integer :: rep_size - - double precision, allocatable, dimension(:) :: deltas - double precision, allocatable, dimension(:) :: all_values - - double precision, allocatable, dimension(:) :: a_cdf - double precision, allocatable, dimension(:) :: b_cdf - integer, allocatable, dimension(:) :: a_cdf_idx - integer, allocatable, dimension(:) :: b_cdf_idx - - rep_size = size(a, dim=1) - allocate (asorted(rep_size, na)) - allocate (bsorted(rep_size, nb)) - allocate (rep(rep_size)) - - allocate (all_values(rep_size*2)) - allocate (deltas(rep_size*2 - 1)) - - allocate (a_cdf(rep_size*2 - 1)) - allocate (b_cdf(rep_size*2 - 1)) - - allocate (a_cdf_idx(rep_size*2 - 1)) - allocate (b_cdf_idx(rep_size*2 - 1)) - - asorted(:, :) = a(:, :) - bsorted(:, :) = b(:, :) - - do i = 1, na - rep(:) = asorted(:, i) - call quicksort(rep, 1, rep_size) - asorted(:, i) = rep(:) - end do - - do i = 1, nb - rep(:) = bsorted(:, i) - call quicksort(rep, 1, rep_size) - bsorted(:, i) = rep(:) - end do - - !$OMP PARALLEL DO PRIVATE(all_values,a_cdf_idx,b_cdf_idx,a_cdf,b_cdf,deltas) - do j = 1, nb - do i = 1, na - - all_values(:rep_size) = asorted(:, i) - all_values(rep_size + 1:) = bsorted(:, j) - - call quicksort(all_values, 1, 2*rep_size) - - do l = 1, 2*rep_size - 1 - deltas(l) = all_values(l + 1) - all_values(l) - end do - - a_cdf_idx = searchsorted(all_values, asorted(:, i)) - b_cdf_idx = searchsorted(all_values, bsorted(:, j)) - - a_cdf(:) = a_cdf_idx(:) - b_cdf(:) = b_cdf_idx(:) - a_cdf(:) = a_cdf(:)/rep_size - b_cdf(:) = b_cdf(:)/rep_size - - ! k(i,j) = exp(-sum(abs(a_cdf-b_cdf)*deltas)/sigma) - k(i, j) = exp(-(sum((abs(a_cdf - b_cdf)**p)*deltas)**(1.0d0/p))**q/sigma) - - end do - end do - !$OMP END PARALLEL DO -end subroutine fwasserstein_kernel diff --git a/src/qmllib/kernels/fkwasserstein.f90 b/src/qmllib/kernels/fkwasserstein.f90 new file mode 100644 index 00000000..73820f54 --- /dev/null +++ b/src/qmllib/kernels/fkwasserstein.f90 @@ -0,0 +1,177 @@ + +module searchtools + + implicit none + +contains + + subroutine searchsorted(all_values, sorted, cdf_idx) + !function searchsorted(all_values, sorted) result(cdf_idx) + + implicit none + + double precision, dimension(:), intent(in) :: all_values + double precision, dimension(:), intent(in) :: sorted + + integer, dimension(:), intent(out) :: cdf_idx + !integer, allocatable, dimension(:) :: cdf_idx +! integer, allocatable, dimension(:) :: searchsorted + + double precision :: val + + integer :: i, j, n, m + + n = size(all_values) - 1 + m = size(sorted) + + !allocate (cdf_idx(n)) + + cdf_idx(:) = 0 + + do i = 1, n + + val = all_values(i) + + do j = 1, m + + !write (*,*) i, j, sorted(j), val + + ! if ((sorted(j) <= val) .and. (val < sorted(j+1))) then + if (sorted(j) > val) then + + cdf_idx(i) = j - 1 + !write(*,*) "found" + exit + + ! endif + else !iif (val > maxval(sorted)) then + cdf_idx(i) = m + end if + + end do + + end do + + end subroutine searchsorted + !end function searchsorted + + recursive subroutine quicksort(a, first, last) + implicit none + double precision :: a(*), x, t + integer first, last + integer i, j + + x = a((first + last)/2) + i = first + j = last + do + do while (a(i) < x) + i = i + 1 + end do + do while (x < a(j)) + j = j - 1 + end do + if (i >= j) exit + t = a(i); a(i) = a(j); a(j) = t + i = i + 1 + j = j - 1 + end do + if (first < i - 1) call quicksort(a, first, i - 1) + if (j + 1 < last) call quicksort(a, j + 1, last) + end subroutine quicksort + +end module searchtools + +subroutine fwasserstein_kernel(a, na, b, nb, k, sigma, p, q) + + use searchtools + implicit none + + double precision, dimension(:, :), intent(in) :: a + double precision, dimension(:, :), intent(in) :: b + + double precision, allocatable, dimension(:, :) :: asorted + double precision, allocatable, dimension(:, :) :: bsorted + + double precision, allocatable, dimension(:) :: rep + + integer, intent(in) :: na, nb + + double precision, dimension(:, :), intent(inout) :: k + double precision, intent(in) :: sigma + + integer, intent(in) :: p + integer, intent(in) :: q + + double precision :: inv_sigma + + integer :: i, j, l + integer :: rep_size + + double precision, allocatable, dimension(:) :: deltas + double precision, allocatable, dimension(:) :: all_values + + double precision, allocatable, dimension(:) :: a_cdf + double precision, allocatable, dimension(:) :: b_cdf + integer, allocatable, dimension(:) :: a_cdf_idx + integer, allocatable, dimension(:) :: b_cdf_idx + + rep_size = size(a, dim=1) + allocate (asorted(rep_size, na)) + allocate (bsorted(rep_size, nb)) + allocate (rep(rep_size)) + + allocate (all_values(rep_size*2)) + allocate (deltas(rep_size*2 - 1)) + + allocate (a_cdf(rep_size*2 - 1)) + allocate (b_cdf(rep_size*2 - 1)) + + allocate (a_cdf_idx(rep_size*2 - 1)) + allocate (b_cdf_idx(rep_size*2 - 1)) + + asorted(:, :) = a(:, :) + bsorted(:, :) = b(:, :) + + do i = 1, na + rep(:) = asorted(:, i) + call quicksort(rep, 1, rep_size) + asorted(:, i) = rep(:) + end do + + do i = 1, nb + rep(:) = bsorted(:, i) + call quicksort(rep, 1, rep_size) + bsorted(:, i) = rep(:) + end do + + !$OMP PARALLEL DO PRIVATE(all_values,a_cdf_idx,b_cdf_idx,a_cdf,b_cdf,deltas) + do j = 1, nb + do i = 1, na + + all_values(:rep_size) = asorted(:, i) + all_values(rep_size + 1:) = bsorted(:, j) + + call quicksort(all_values, 1, 2*rep_size) + + do l = 1, 2*rep_size - 1 + deltas(l) = all_values(l + 1) - all_values(l) + end do + + !a_cdf_idx = searchsorted(all_values, asorted(:, i)) + !b_cdf_idx = searchsorted(all_values, bsorted(:, j)) + call searchsorted(all_values, asorted(:, i), a_cdf_idx) + call searchsorted(all_values, bsorted(:, j), b_cdf_idx) + + a_cdf(:) = a_cdf_idx(:) + b_cdf(:) = b_cdf_idx(:) + a_cdf(:) = a_cdf(:)/rep_size + b_cdf(:) = b_cdf(:)/rep_size + + ! k(i,j) = exp(-sum(abs(a_cdf-b_cdf)*deltas)/sigma) + k(i, j) = exp(-(sum((abs(a_cdf - b_cdf)**p)*deltas)**(1.0d0/p))**q/sigma) + + end do + end do + !$OMP END PARALLEL DO +end subroutine fwasserstein_kernel diff --git a/src/qmllib/kernels/gradient_kernels.py b/src/qmllib/kernels/gradient_kernels.py index 59f7f056..fa105d48 100644 --- a/src/qmllib/kernels/gradient_kernels.py +++ b/src/qmllib/kernels/gradient_kernels.py @@ -295,7 +295,6 @@ def get_atomic_local_kernel( 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_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -361,7 +360,6 @@ def get_atomic_local_gradient_kernel( 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} = \\frac{\\part}{\\part x}\\sum_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -443,7 +441,6 @@ def get_atomic_local_gradient_kernel( 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}`: :math:`K_{ij} = \\frac{\\part}{\\part x}\\sum_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -517,7 +514,6 @@ def get_gdml_kernel( Q2: List[List[int]], SIGMA: float, ) -> ndarray: - """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: :math:`K_{Ij} = \\frac{\\part^2}{\\part x_i\\part x_j}\\sum_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -601,7 +597,6 @@ def get_gdml_kernel( 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}`: :math:`K_{Ij} = \\frac{\\part^2}{\\part x_i\\part x_j}\\sum_{J\\in j}\\exp \\big( -\\frac{\\|X_I - X_J\\|_2^2}{2\\sigma^2} \\big)` @@ -662,7 +657,6 @@ def get_gp_kernel( Q2: List[Union[ndarray, List[int]]], SIGMA: float, ) -> ndarray: - """Calculates the Gaussian kernel matrix K with the local decomposition where :math:`K_{ij}`: This kernel corresponds to a Gaussian process regression (GPR) approach. @@ -743,7 +737,6 @@ def get_gp_kernel( 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. The kernel has four blocks, consisting of the 0'th, 1st and 2nd derivatives. diff --git a/src/qmllib/representations/fchl/ffchl_force_kernels.f90 b/src/qmllib/representations/fchl/ffchl_force_kernels.f90 index d939c0a4..e7d9062a 100644 --- a/src/qmllib/representations/fchl/ffchl_force_kernels.f90 +++ b/src/qmllib/representations/fchl/ffchl_force_kernels.f90 @@ -118,6 +118,9 @@ subroutine fget_gaussian_process_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, integer :: maxneigh1 integer :: maxneigh2 + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + ! Angular normalization constant ang_norm2 = get_angular_norm2(t_width) @@ -132,8 +135,14 @@ subroutine fget_gaussian_process_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, pmax2 = get_pmax_displaced(x2, n2) ! Get two-body weight function - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), 3, size(x2, dim=3), maxval(n2), maxval(n2), maxval(nneigh2))) + + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + + ! ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + ! ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, maxval(n1), pmax1, order, maxneigh1)) @@ -152,14 +161,21 @@ subroutine fget_gaussian_process_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, & cut_distance, cosp2, sinp2, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + allocate (self_scalar2(nm2, 3, size(x2, dim=3), maxval(n2), maxval(n2))) + + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + call get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & - & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + ! & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,s12) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,s12,ktmp) do a = 1, nm1 na = n1(a) do j1 = 1, na @@ -176,9 +192,14 @@ subroutine fget_gaussian_process_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernels(:, a, b) = kernels(:, a, b) & - & + kernel(self_scalar1(a, j1), self_scalar1(b, j2), s12, & - & kernel_idx, parameters) + ktmp = 0.0d0 + call kernel(self_scalar1(a, j1), self_scalar1(b, j2), s12, & + & kernel_idx, parameters, ktmp) + + kernels(:, a, b) = kernels(:, a, b) + ktmp + ! kernels(:, a, b) = kernels(:, a, b) & + ! & + kernel(self_scalar1(a, j1), self_scalar1(b, j2), s12, & + ! & kernel_idx, parameters) end do end do @@ -186,7 +207,7 @@ subroutine fget_gaussian_process_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, end do !$OMP END PARALLEL do - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,xyz_pm2,s12),& + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,xyz_pm2,s12,ktmp),& !$OMP& PRIVATE(idx1,idx2) do a = 1, nm1 na = n1(a) @@ -210,19 +231,25 @@ subroutine fget_gaussian_process_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + ktmp = 0.0d0 + call kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & + & kernel_idx, parameters, ktmp) + if (pm2 == 2) then - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) + ktmp + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & + ! & kernel_idx, parameters) kernels(:, idx2, idx1) = kernels(:, idx1, idx2) else - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & - kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) - ktmp + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & - kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & + ! & kernel_idx, parameters) kernels(:, idx2, idx1) = kernels(:, idx1, idx2) @@ -240,7 +267,7 @@ subroutine fget_gaussian_process_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, kernels(:, :nm1, nm1 + 1:) = kernels(:, :nm1, nm1 + 1:)/(2*dx) kernels(:, nm1 + 1:, :nm1) = kernels(:, nm1 + 1:, :nm1)/(2*dx) - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,xyz_pm1,xyz_pm2,s12),& + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,xyz_pm1,xyz_pm2,s12,ktmp),& !$OMP& PRIVATE(idx1,idx2) do a = 1, nm1 na = n1(a) @@ -268,27 +295,39 @@ subroutine fget_gaussian_process_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + ktmp = 0.0d0 + call kernel(self_scalar2(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + & kernel_idx, parameters, ktmp) + if (pm1 == pm2) then - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & + kernel(self_scalar2(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) + ktmp + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & + kernel(self_scalar2(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) if (a /= b) then - kernels(:, idx2, idx1) = kernels(:, idx2, idx1) & - & + kernel(self_scalar2(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx2, idx1) = kernels(:, idx2, idx1) + ktmp + + ! kernels(:, idx2, idx1) = kernels(:, idx2, idx1) & + ! & + kernel(self_scalar2(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) end if else - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & - kernel(self_scalar2(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) - ktmp + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & - kernel(self_scalar2(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) if (a /= b) then - kernels(:, idx2, idx1) = kernels(:, idx2, idx1) & - & - kernel(self_scalar2(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx2, idx1) = kernels(:, idx2, idx1) - ktmp + + ! kernels(:, idx2, idx1) = kernels(:, idx2, idx1) & + ! & - kernel(self_scalar2(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) end if end if @@ -307,6 +346,16 @@ subroutine fget_gaussian_process_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, kernels(:, nm1 + 1:, nm1 + 1:) = kernels(:, nm1 + 1:, nm1 + 1:)/(4*dx**2) + deallocate (ktmp) + deallocate (ksi1) + deallocate (ksi2) + deallocate (cosp1) + deallocate (sinp1) + deallocate (cosp2) + deallocate (sinp2) + deallocate (self_scalar1) + deallocate (self_scalar2) + end subroutine fget_gaussian_process_kernels_fchl subroutine fget_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, & @@ -426,6 +475,10 @@ subroutine fget_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nn integer :: maxneigh1 integer :: maxneigh2 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + kernels = 0.0d0 ! Angular normalization constant @@ -440,8 +493,14 @@ subroutine fget_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nn pmax2 = get_pmax_displaced(x2, n2) ! Get two-body weight function - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), 3, size(x2, dim=3), maxval(n2), maxval(n2), maxval(nneigh2))) + + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + + ! ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + ! ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, maxval(n1), pmax1, order, maxneigh1)) @@ -460,12 +519,18 @@ subroutine fget_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nn & cut_distance, cosp2, sinp2, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + allocate (self_scalar2(nm2, 3, size(x2, dim=3), maxval(n2), maxval(n2))) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + call get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & - & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + ! & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,xyz_pm2,s12),& !$OMP& PRIVATE(idx1,idx2) @@ -490,15 +555,24 @@ subroutine fget_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nn & cosp1(a, j1, :, :, :), cosp2(b, xyz_pm2, i2, j2, :, :, :), & & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + + ktmp = 0.0d0 + call kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & + & kernel_idx, parameters, ktmp) + if (pm2 == 2) then - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) + ktmp + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & + ! & kernel_idx, parameters) else - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & - kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) - ktmp + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & - kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & + ! & kernel_idx, parameters) end if end do @@ -512,6 +586,16 @@ subroutine fget_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nn kernels = kernels/(2*dx) + deallocate (ktmp) + deallocate (ksi1) + deallocate (ksi2) + deallocate (cosp1) + deallocate (sinp1) + deallocate (cosp2) + deallocate (sinp2) + deallocate (self_scalar1) + deallocate (self_scalar2) + end subroutine fget_local_gradient_kernels_fchl subroutine fget_local_hessian_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, & @@ -631,6 +715,10 @@ subroutine fget_local_hessian_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nne integer :: maxneigh1 integer :: maxneigh2 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + ! Angular normalization constant ang_norm2 = get_angular_norm2(t_width) @@ -645,8 +733,13 @@ subroutine fget_local_hessian_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nne pmax2 = get_pmax_displaced(x2, n2) ! Get two-body weight function - ksi1 = get_ksi_displaced(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), 3, size(x1, dim=3), maxval(n1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), 3, size(x2, dim=3), maxval(n2), maxval(n2), maxval(nneigh2))) + call get_ksi_displaced(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + + ! ksi1 = get_ksi_displaced(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + ! ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, 3*2, maxval(n1), maxval(n1), pmax1, order, maxval(nneigh1))) @@ -665,12 +758,18 @@ subroutine fget_local_hessian_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nne & cut_distance, cosp2, sinp2, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar_displaced(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, & - & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, 3, size(x1, dim=3), maxval(n1), maxval(n1))) + allocate (self_scalar2(nm2, 3, size(x2, dim=3), maxval(n2), maxval(n2))) + call get_selfscalar_displaced(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, & + & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + call get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & - & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar1 = get_selfscalar_displaced(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, & + ! & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + ! & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,xyz_pm1,xyz_pm2,s12),& !$OMP& PRIVATE(idx1,idx2) @@ -700,16 +799,24 @@ subroutine fget_local_hessian_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nne & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + ktmp = 0.0d0 + call kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + & kernel_idx, parameters, ktmp) + if (pm1 == pm2) then - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & + kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) + ktmp + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & + kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) else - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & - kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) - ktmp + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & - kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) end if @@ -727,6 +834,16 @@ subroutine fget_local_hessian_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nne kernels = kernels/(4*dx**2) + deallocate (ktmp) + deallocate (ksi1) + deallocate (ksi2) + deallocate (cosp1) + deallocate (sinp1) + deallocate (cosp2) + deallocate (sinp2) + deallocate (self_scalar1) + deallocate (self_scalar2) + end subroutine fget_local_hessian_kernels_fchl subroutine fget_local_symmetric_hessian_kernels_fchl(x1, verbose, n1, nneigh1, & @@ -834,6 +951,10 @@ subroutine fget_local_symmetric_hessian_kernels_fchl(x1, verbose, n1, nneigh1, & ! Max number of neighbors integer :: maxneigh1 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + ! Angular normalization constant ang_norm2 = get_angular_norm2(t_width) @@ -846,7 +967,9 @@ subroutine fget_local_symmetric_hessian_kernels_fchl(x1, verbose, n1, nneigh1, & pmax1 = get_pmax_displaced(x1, n1) ! Get two-body weight function - ksi1 = get_ksi_displaced(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), 3, size(x1, dim=3), maxval(n1), maxval(n1), maxval(nneigh1))) + call get_ksi_displaced(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + ! ksi1 = get_ksi_displaced(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, 3*2, maxval(n1), maxval(n1), pmax1, order, maxval(nneigh1))) @@ -857,8 +980,11 @@ subroutine fget_local_symmetric_hessian_kernels_fchl(x1, verbose, n1, nneigh1, & & cosp1, sinp1, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar_displaced(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width,& - & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, 3, size(x1, dim=3), maxval(n1), maxval(n1))) + call get_selfscalar_displaced(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width,& + & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + ! self_scalar1 = get_selfscalar_displaced(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width,& + ! & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,xyz_pm1,xyz_pm2,s12),& !$OMP& PRIVATE(idx1,idx2) @@ -888,26 +1014,38 @@ subroutine fget_local_symmetric_hessian_kernels_fchl(x1, verbose, n1, nneigh1, & & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + ktmp = 0.0d0 + call kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar1(b, xyz2, pm2, i2, j2), s12,& + & kernel_idx, parameters, ktmp) + if (pm1 == pm2) then - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & + kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar1(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) + ktmp + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & + kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar1(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) + if (a /= b) then - kernels(:, idx2, idx1) = kernels(:, idx2, idx1) & - & + kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar1(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx2, idx1) = kernels(:, idx2, idx1) + ktmp + + ! kernels(:, idx2, idx1) = kernels(:, idx2, idx1) & + ! & + kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar1(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) end if else - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & - kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar1(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) - ktmp + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & - kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar1(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) if (a /= b) then - kernels(:, idx2, idx1) = kernels(:, idx2, idx1) & - & - kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar1(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx2, idx1) = kernels(:, idx2, idx1) - ktmp + + ! kernels(:, idx2, idx1) = kernels(:, idx2, idx1) & + ! & - kernel(self_scalar1(a, xyz1, pm1, i1, j1), self_scalar1(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) end if end if @@ -926,6 +1064,12 @@ subroutine fget_local_symmetric_hessian_kernels_fchl(x1, verbose, n1, nneigh1, & kernels = kernels/(4*dx**2) + deallocate (ktmp) + deallocate (ksi1) + deallocate (cosp1) + deallocate (sinp1) + deallocate (self_scalar1) + end subroutine fget_local_symmetric_hessian_kernels_fchl subroutine fget_force_alphas_fchl(x1, x2, verbose, forces, energies, n1, n2, & @@ -1068,6 +1212,10 @@ subroutine fget_force_alphas_fchl(x1, x2, verbose, forces, energies, n1, n2, & ! Kernel between molecules and atom double precision, allocatable, dimension(:, :, :) :: kernel_ma + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + inv_2dx = 1.0d0/(2.0d0*dx) ! Angular normalization constant @@ -1082,8 +1230,13 @@ subroutine fget_force_alphas_fchl(x1, x2, verbose, forces, energies, n1, n2, & pmax2 = get_pmax_displaced(x2, n2) ! Get two-body weight function - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), 3, size(x2, dim=3), maxval(n2), maxval(n2), maxval(nneigh2))) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + + ! ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + ! ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, maxval(n1), pmax1, order, maxneigh1)) @@ -1102,12 +1255,18 @@ subroutine fget_force_alphas_fchl(x1, x2, verbose, forces, energies, n1, n2, & & cut_distance, cosp2, sinp2, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + allocate (self_scalar2(nm2, 3, size(x2, dim=3), maxval(n2), maxval(n2))) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + call get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & - & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + ! & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) allocate (kernel_delta(na1, na1, nsigmas)) allocate (y(na1, nsigmas)) @@ -1147,15 +1306,25 @@ subroutine fget_force_alphas_fchl(x1, x2, verbose, forces, energies, n1, n2, & & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + ktmp = 0.0d0 + call kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & + kernel_idx, parameters, ktmp) + if (pm2 == 2) then - kernel_delta(idx1, idx2, :) = kernel_delta(idx1, idx2, :) & - & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & - kernel_idx, parameters)*inv_2dx + kernel_delta(idx1, idx2, :) = kernel_delta(idx1, idx2, :) + ktmp*inv_2dx + + ! kernel_delta(idx1, idx2, :) = kernel_delta(idx1, idx2, :) & + ! & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & + ! kernel_idx, parameters)*inv_2dx + else - kernel_delta(idx1, idx2, :) = kernel_delta(idx1, idx2, :) & - & - kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & - kernel_idx, parameters)*inv_2dx + + kernel_delta(idx1, idx2, :) = kernel_delta(idx1, idx2, :) - ktmp*inv_2dx + + ! kernel_delta(idx1, idx2, :) = kernel_delta(idx1, idx2, :) & + ! & - kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12, & + ! kernel_idx, parameters)*inv_2dx end if @@ -1212,9 +1381,15 @@ subroutine fget_force_alphas_fchl(x1, x2, verbose, forces, energies, n1, n2, & & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernel_MA(b, idx1, :) = kernel_MA(b, idx1, :) & - & + kernel(self_scalar1(a, i), self_scalar1(b, j), s12, & - kernel_idx, parameters) + ktmp = 0.0d0 + call kernel(self_scalar1(a, i), self_scalar1(b, j), s12, & + kernel_idx, parameters, ktmp) + + kernel_MA(b, idx1, :) = kernel_MA(b, idx1, :) + ktmp + + ! kernel_MA(b, idx1, :) = kernel_MA(b, idx1, :) & + ! & + kernel(self_scalar1(a, i), self_scalar1(b, j), s12, & + ! kernel_idx, parameters) end do end do @@ -1278,6 +1453,16 @@ subroutine fget_force_alphas_fchl(x1, x2, verbose, forces, energies, n1, n2, & deallocate (y) deallocate (kernel_scratch) + deallocate (ktmp) + deallocate (ksi1) + deallocate (ksi2) + deallocate (cosp1) + deallocate (sinp1) + deallocate (cosp2) + deallocate (sinp2) + deallocate (self_scalar1) + deallocate (self_scalar2) + end subroutine fget_force_alphas_fchl subroutine fget_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, & @@ -1398,6 +1583,10 @@ subroutine fget_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nnei integer :: maxneigh1 integer :: maxneigh2 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + kernels = 0.0d0 ! Angular normalization constant @@ -1412,8 +1601,13 @@ subroutine fget_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nnei pmax2 = get_pmax_displaced(x2, n2) ! Get two-body weight function - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), 3, size(x2, dim=3), maxval(n2), maxval(n2), maxval(nneigh2))) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + + ! ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + ! ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, maxval(n1), pmax1, order, maxneigh1)) @@ -1432,12 +1626,18 @@ subroutine fget_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nnei & cut_distance, cosp2, sinp2, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + allocate (self_scalar2(nm2, 3, size(x2, dim=3), maxval(n2), maxval(n2))) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + call get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & - & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + ! & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,xyz_pm2,s12),& !$OMP& PRIVATE(idx1,idx2,idx1_start,idx1_end,idx2_start,idx2_end) @@ -1473,15 +1673,23 @@ subroutine fget_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nnei & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + ktmp = 0.0d0 + call kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + & kernel_idx, parameters, ktmp) + if (pm2 == 2) then - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) + ktmp + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) else - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & - kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters) + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) - ktmp + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & - kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters) end if @@ -1496,6 +1704,16 @@ subroutine fget_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nnei kernels = kernels/(2*dx) + deallocate (ktmp) + deallocate (ksi1) + deallocate (ksi2) + deallocate (cosp1) + deallocate (sinp1) + deallocate (cosp2) + deallocate (sinp2) + deallocate (self_scalar1) + deallocate (self_scalar2) + end subroutine fget_atomic_local_gradient_kernels_fchl subroutine fget_atomic_local_gradient_5point_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, & @@ -1619,6 +1837,10 @@ subroutine fget_atomic_local_gradient_5point_kernels_fchl(x1, x2, verbose, n1, n ! For numerical differentiation double precision, parameter, dimension(5) :: fact = (/1.0d0, -8.0d0, 0.0d0, 8.0d0, -1.0d0/) + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + ! fact(1) = 1.0d0 ! fact(2) = -8.0d0 ! fact(3) = 0.0d0 @@ -1639,8 +1861,12 @@ subroutine fget_atomic_local_gradient_5point_kernels_fchl(x1, x2, verbose, n1, n pmax2 = get_pmax_displaced(x2, n2) ! Get two-body weight function - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), 3, size(x2, dim=3), maxval(n2), maxval(n2), maxval(nneigh2))) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + ! ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + ! ksi2 = get_ksi_displaced(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, maxval(n1), pmax1, order, maxneigh1)) @@ -1659,12 +1885,18 @@ subroutine fget_atomic_local_gradient_5point_kernels_fchl(x1, x2, verbose, n1, n & cut_distance, cosp2, sinp2, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + allocate (self_scalar2(nm2, 3, size(x2, dim=3), maxval(n2), maxval(n2))) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + call get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & - & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + ! self_scalar2 = get_selfscalar_displaced(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, & + ! & d_width, cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) !$OMP PARALLEL DO schedule(dynamic) PRIVATE(na,nb,xyz_pm2,s12),& !$OMP& PRIVATE(idx1,idx2,idx1_start,idx1_end,idx2_start,idx2_end) @@ -1704,9 +1936,15 @@ subroutine fget_atomic_local_gradient_5point_kernels_fchl(x1, x2, verbose, n1, n & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & - & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& - & kernel_idx, parameters)*fact(pm2) + ktmp = 0.0d0 + call kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + & kernel_idx, parameters, ktmp) + + kernels(:, idx1, idx2) = kernels(:, idx1, idx2) + ktmp*fact(pm2) + + ! kernels(:, idx1, idx2) = kernels(:, idx1, idx2) & + ! & + kernel(self_scalar1(a, j1), self_scalar2(b, xyz2, pm2, i2, j2), s12,& + ! & kernel_idx, parameters)*fact(pm2) end do end do @@ -1722,4 +1960,14 @@ subroutine fget_atomic_local_gradient_5point_kernels_fchl(x1, x2, verbose, n1, n kernels = kernels/(12*dx) + deallocate (ktmp) + deallocate (ksi1) + deallocate (ksi2) + deallocate (cosp1) + deallocate (sinp1) + deallocate (cosp2) + deallocate (sinp2) + deallocate (self_scalar1) + deallocate (self_scalar2) + end subroutine fget_atomic_local_gradient_5point_kernels_fchl diff --git a/src/qmllib/representations/fchl/ffchl_kernels.f90 b/src/qmllib/representations/fchl/ffchl_kernels.f90 index b0984c15..b13d88d4 100644 --- a/src/qmllib/representations/fchl/ffchl_kernels.f90 +++ b/src/qmllib/representations/fchl/ffchl_kernels.f90 @@ -1,6 +1,3 @@ - -! Inspiration from: -! http://crsouza.com/2010/03/17/kernel-functions-for-machine-learning-applications/#kernel_functions module ffchl_kernels implicit none @@ -251,7 +248,8 @@ subroutine polynomial2_kernel(s12, parameters, k) end subroutine polynomial2_kernel - function kernel(s11, s22, s12, kernel_idx, parameters) result(k) + !function kernel(s11, s22, s12, kernel_idx, parameters) result(k) + subroutine kernel(s11, s22, s12, kernel_idx, parameters, k) use ffchl_kernel_types @@ -264,10 +262,11 @@ function kernel(s11, s22, s12, kernel_idx, parameters) result(k) double precision, intent(in), dimension(:, :) :: parameters integer :: n - double precision, allocatable, dimension(:) :: k + !double precision, allocatable, dimension(:) :: k + double precision, dimension(:) :: k n = size(parameters, dim=1) - allocate (k(n)) + !allocate (k(n)) if (kernel_idx == GAUSSIAN) then call gaussian_kernel(s11, s22, s12, parameters, k) @@ -307,6 +306,6 @@ function kernel(s11, s22, s12, kernel_idx, parameters) result(k) stop end if - end function kernel + end subroutine kernel end module ffchl_kernels diff --git a/src/qmllib/representations/fchl/ffchl_electric_field_kernels.f90 b/src/qmllib/representations/fchl/ffchl_kernels_ef.f90 similarity index 83% rename from src/qmllib/representations/fchl/ffchl_electric_field_kernels.f90 rename to src/qmllib/representations/fchl/ffchl_kernels_ef.f90 index 320b9f35..2fbf58a9 100644 --- a/src/qmllib/representations/fchl/ffchl_electric_field_kernels.f90 +++ b/src/qmllib/representations/fchl/ffchl_kernels_ef.f90 @@ -3,8 +3,9 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 & distance_scale, angular_scale, alchemy, two_body_power, three_body_power, ef_scale,& & df, kernel_idx, parameters, kernels) - use ffchl_module, only: scalar, get_angular_norm2, get_pmax, get_ksi_ef_field, init_cosp_sinp_ef_field, & - & get_selfscalar, get_ksi_ef, init_cosp_sinp_ef + use ffchl_module, only: scalar, get_angular_norm2, get_pmax, get_selfscalar + + use ffchl_module_ef, only: get_ksi_ef_field, init_cosp_sinp_ef_field, get_ksi_ef, init_cosp_sinp_ef use ffchl_kernels, only: kernel @@ -119,6 +120,13 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 integer :: maxneigh1 integer :: maxneigh2 + ! Work for kernel + integer :: n + double precision, allocatable, dimension(:) :: ktmp + + n = size(parameters, dim=1) + allocate (ktmp(n)) + kernels(:, :, :) = 0.0d0 ! Get max number of neighbors @@ -133,14 +141,25 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 pmax2 = get_pmax(x2, n2) ! Get two-body weight function - ksi1 = get_ksi_ef_field(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, & - & f1, ef_scale, verbose) - ksi1_ef = get_ksi_ef(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, ef_scale, df, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + call get_ksi_ef_field(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, & + & f1, ef_scale, verbose, ksi1) + !ksi1 = get_ksi_ef_field(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, & + ! & f1, ef_scale, verbose) + + allocate (ksi1_ef(size(x1, dim=1), 3, 2, maxval(n1), maxval(nneigh1))) + call get_ksi_ef(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, ef_scale, df, verbose, ksi1_ef) + !ksi1_ef = get_ksi_ef(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, ef_scale, df, verbose) ! Get two-body weight function - ksi2 = get_ksi_ef_field(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, & - & f2, ef_scale, verbose) - ksi2_ef = get_ksi_ef(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, ef_scale, df, verbose) + allocate (ksi2(size(x2, dim=1), maxval(n2), maxval(nneigh2))) + allocate (ksi2_ef(size(x2, dim=1), 3, 2, maxval(n2), maxval(nneigh2))) + call get_ksi_ef_field(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, & + & f2, ef_scale, verbose, ksi2) + call get_ksi_ef(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, ef_scale, df, verbose, ksi2_ef) + !ksi2 = get_ksi_ef_field(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, & + ! & f2, ef_scale, verbose) + !ksi2_ef = get_ksi_ef(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, ef_scale, df, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, maxval(n1), pmax1, order, maxneigh1)) @@ -174,13 +193,20 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 call init_cosp_sinp_ef(x2, n2, nneigh2, three_body_power, order, cut_start, cut_distance, & & cosp2_ef, sinp2_ef, ef_scale, df, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + allocate (self_scalar2(nm2, maxval(n2))) + ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + !self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + call get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) + !self_scalar2 = get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) ! Self-scalar derivatives allocate (self_scalar1_ef(nm1, 3, 2, maxval(n1))) @@ -238,9 +264,13 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernels(:, a, b) = kernels(:, a, b) & - & + kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & - & kernel_idx, parameters) + ktmp = 0.0d0 + call kernel(self_scalar1(a, i), self_scalar2(b, j), s12, kernel_idx, parameters, ktmp) + kernels(:, a, b) = kernels(:, a, b) + ktmp + + !kernels(:, a, b) = kernels(:, a, b) & + ! & + kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & + ! & kernel_idx, parameters) end do end do @@ -248,7 +278,7 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 end do !$OMP END PARALLEL DO - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,idx_a,idx_b) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,idx_a,idx_b,ktmp) do a = 1, nm1 ni = n1(a) do i = 1, ni @@ -271,17 +301,24 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + ktmp = 0.0d0 + call kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, kernel_idx, parameters, ktmp) + if (pm == 1) then - kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & - & + kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, & - & kernel_idx, parameters)/(2*df) + !kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & + ! & + kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, & + ! & kernel_idx, parameters)/(2*df) + + kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) + ktmp/(2*df) else - kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & - & - kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, & - & kernel_idx, parameters)/(2*df) + !kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & + ! & - kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, & + ! & kernel_idx, parameters)/(2*df) + + kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) - ktmp/(2*df) end if @@ -293,7 +330,7 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 end do !$OMP END PARALLEL DO - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,idx_a,idx_b) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,idx_a,idx_b,ktmp) do a = 1, nm2 ni = n2(a) @@ -317,19 +354,27 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + ktmp = 0.0d0 + call kernel(self_scalar2(a, i), self_scalar1_ef(b, xyz, pm, j), s12, & + & kernel_idx, parameters, ktmp) + if (pm == 1) then ! kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & - kernels(:, idx_b, idx_a) = kernels(:, idx_b, idx_a) & - & + kernel(self_scalar2(a, i), self_scalar1_ef(b, xyz, pm, j), s12, & - & kernel_idx, parameters)/(2*df) + kernels(:, idx_b, idx_a) = kernels(:, idx_b, idx_a) + ktmp/(2*df) + + !kernels(:, idx_b, idx_a) = kernels(:, idx_b, idx_a) & + ! & + kernel(self_scalar2(a, i), self_scalar1_ef(b, xyz, pm, j), s12, & + ! & kernel_idx, parameters)/(2*df) else ! kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & - kernels(:, idx_b, idx_a) = kernels(:, idx_b, idx_a) & - & - kernel(self_scalar2(a, i), self_scalar1_ef(b, xyz, pm, j), s12, & - & kernel_idx, parameters)/(2*df) + kernels(:, idx_b, idx_a) = kernels(:, idx_b, idx_a) - ktmp/(2*df) + + !kernels(:, idx_b, idx_a) = kernels(:, idx_b, idx_a) & + ! & - kernel(self_scalar2(a, i), self_scalar1_ef(b, xyz, pm, j), s12, & + ! & kernel_idx, parameters)/(2*df) end if @@ -363,19 +408,30 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 s12 = scalar(x1(a, i, :, :), x2(b, j, :, :), & & nneigh1(a, i), nneigh2(b, j), & & ksi1_ef(a, xyz1, pm1, i, :), ksi2_ef(b, xyz2, pm2, j, :), & - & sinp1_ef(a, xyz1, pm1, i, :, :, :), sinp2_ef(b, xyz2, pm2, j, :, :, :), & - & cosp1_ef(a, xyz1, pm1, i, :, :, :), cosp2_ef(b, xyz2, pm2, j, :, :, :), & + & sinp1_ef(a, xyz1, pm1, i, :, :, :), sinp2_ef(b, xyz2, pm2, j, :, :, :), & + & cosp1_ef(a, xyz1, pm1, i, :, :, :), cosp2_ef(b, xyz2, pm2, j, :, :, :), & & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + + ktmp = 0.0d0 + call kernel(self_scalar1_ef(a, xyz1, pm1, i), self_scalar2_ef(b, xyz2, pm2, j), s12, & + & kernel_idx, parameters, ktmp) + if (pm1 == pm2) then - kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & - & + kernel(self_scalar1_ef(a, xyz1, pm1, i), self_scalar2_ef(b, xyz2, pm2, j), s12, & - & kernel_idx, parameters)/(4*df**2) + kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) + ktmp/(4*df**2) + + ! kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & + !& + kernel(self_scalar1_ef(a, xyz1, pm1, i), self_scalar2_ef(b, xyz2, pm2, j), s12, & + ! & kernel_idx, parameters)/(4*df**2) + else - kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & - & - kernel(self_scalar1_ef(a, xyz1, pm1, i), self_scalar2_ef(b, xyz2, pm2, j), s12, & - & kernel_idx, parameters)/(4*df**2) + + kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) - ktmp/(4*df**2) + + ! kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & + !& - kernel(self_scalar1_ef(a, xyz1, pm1, i), self_scalar2_ef(b, xyz2, pm2, j), s12, & + ! & kernel_idx, parameters)/(4*df**2) end if @@ -389,6 +445,7 @@ subroutine fget_ef_gaussian_process_kernels_fchl(x1, x2, verbose, f1, f2, n1, n2 end do !$OMP END PARALLEL DO + deallocate (ktmp) deallocate (self_scalar1) deallocate (self_scalar1_ef) deallocate (ksi1) @@ -405,9 +462,11 @@ subroutine fget_ef_atomic_local_kernels_fchl(x1, x2, verbose, f2, n1, n2, nneigh & distance_scale, angular_scale, alchemy, two_body_power, three_body_power, ef_scale,& & kernel_idx, parameters, kernels) - use ffchl_module, only: scalar, get_angular_norm2, get_pmax, get_ksi_ef_field, init_cosp_sinp_ef_field, & + use ffchl_module, only: scalar, get_angular_norm2, get_pmax, & & get_selfscalar, get_ksi, init_cosp_sinp + use ffchl_module_ef, only: get_ksi_ef_field, init_cosp_sinp_ef_field + use ffchl_kernels, only: kernel implicit none @@ -509,6 +568,13 @@ subroutine fget_ef_atomic_local_kernels_fchl(x1, x2, verbose, f2, n1, n2, nneigh integer :: maxneigh1 integer :: maxneigh2 + ! Work for kernel + integer :: n + double precision, allocatable, dimension(:) :: ktmp + + n = size(parameters, dim=1) + allocate (ktmp(n)) + kernels(:, :, :) = 0.0d0 ! Get max number of neighbors @@ -522,12 +588,18 @@ subroutine fget_ef_atomic_local_kernels_fchl(x1, x2, verbose, f2, n1, n2, nneigh pmax1 = get_pmax(x1, n1) pmax2 = get_pmax(x2, n2) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), maxval(n2), maxval(nneigh2))) + ! Get two-body weight function - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + ! ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) ! Get two-body weight function - ksi2 = get_ksi_ef_field(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, & - & f2, ef_scale, verbose) + call get_ksi_ef_field(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, & + & f2, ef_scale, verbose, ksi2) + !ksi2 = get_ksi_ef_field(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, & + ! & f2, ef_scale, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, maxval(n1), pmax1, order, maxneigh1)) @@ -546,16 +618,23 @@ subroutine fget_ef_atomic_local_kernels_fchl(x1, x2, verbose, f2, n1, n2, nneigh & cosp2, sinp2, f2, ef_scale, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + allocate (self_scalar2(nm2, maxval(n2))) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) + call get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) + + !self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + !self_scalar2 = get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) kernels(:, :, :) = 0.0d0 - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(ni,nj,idx1,s12) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(ni,nj,idx1,s12,ktmp) do a = 1, nm1 ni = n1(a) do i = 1, ni @@ -573,9 +652,15 @@ subroutine fget_ef_atomic_local_kernels_fchl(x1, x2, verbose, f2, n1, n2, nneigh & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernels(:, idx1, b) = kernels(:, idx1, b) & - & + kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & - & kernel_idx, parameters) + ktmp = 0.0d0 + call kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & + & kernel_idx, parameters, ktmp) + + kernels(:, idx1, b) = kernels(:, idx1, b) + ktmp + + !kernels(:, idx1, b) = kernels(:, idx1, b) & + ! & + kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & + ! & kernel_idx, parameters) end do end do @@ -584,6 +669,7 @@ subroutine fget_ef_atomic_local_kernels_fchl(x1, x2, verbose, f2, n1, n2, nneigh end do !$OMP END PARALLEL DO + deallocate (ktmp) deallocate (self_scalar1) deallocate (self_scalar2) deallocate (ksi1) @@ -592,6 +678,7 @@ subroutine fget_ef_atomic_local_kernels_fchl(x1, x2, verbose, f2, n1, n2, nneigh deallocate (cosp2) deallocate (sinp1) deallocate (sinp2) + end subroutine fget_ef_atomic_local_kernels_fchl subroutine fget_ef_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, nm1, nm2, na1, nsigmas, & @@ -599,8 +686,9 @@ subroutine fget_ef_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, n & distance_scale, angular_scale, alchemy, two_body_power, three_body_power, ef_scale,& & df, kernel_idx, parameters, kernels) - use ffchl_module, only: scalar, get_angular_norm2, get_pmax, get_ksi, init_cosp_sinp, & - & get_selfscalar, get_ksi_ef, init_cosp_sinp_ef + use ffchl_module, only: scalar, get_angular_norm2, get_pmax, get_ksi, init_cosp_sinp, get_selfscalar + + use ffchl_module_ef, only: get_ksi_ef, init_cosp_sinp_ef use ffchl_kernels, only: kernel @@ -704,6 +792,13 @@ subroutine fget_ef_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, n integer :: maxneigh1 integer :: maxneigh2 + ! Work for kernel + integer :: n + double precision, allocatable, dimension(:) :: ktmp + + n = size(parameters, dim=1) + allocate (ktmp(n)) + kernels(:, :, :) = 0.0d0 ! Get max number of neighbors @@ -718,8 +813,12 @@ subroutine fget_ef_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, n pmax2 = get_pmax(x2, n2) ! Get two-body weight function - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2_ef = get_ksi_ef(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, ef_scale, df, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2_ef(size(x2, dim=1), 3, 2, maxval(n2), maxval(nneigh2))) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi_ef(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, ef_scale, df, verbose, ksi2_ef) + !ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + !ksi2_ef = get_ksi_ef(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, ef_scale, df, verbose) ! Allocate three-body Fourier terms allocate (cosp1(nm1, maxval(n1), pmax1, order, maxneigh1)) @@ -738,8 +837,11 @@ subroutine fget_ef_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, n & cosp2_ef, sinp2_ef, ef_scale, df, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + !self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) allocate (self_scalar2_ef(nm2, 3, 2, maxval(n2))) do a = 1, nm2 @@ -760,7 +862,7 @@ subroutine fget_ef_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, n end do end do - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,idx_a,idx_b) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,idx_a,idx_b,ktmp) do a = 1, nm1 ni = n1(a) @@ -784,17 +886,25 @@ subroutine fget_ef_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, n & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) + ktmp = 0.0d0 + call kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, & + & kernel_idx, parameters, ktmp) + if (pm == 1) then - kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & - & + kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, & - & kernel_idx, parameters) + kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) + ktmp + + !kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & + ! & + kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, & + ! & kernel_idx, parameters) else - kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & - & - kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, & - & kernel_idx, parameters) + kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) - ktmp + + !kernels(:, idx_a, idx_b) = kernels(:, idx_a, idx_b) & + ! & - kernel(self_scalar1(a, i), self_scalar2_ef(b, xyz, pm, j), s12, & + ! & kernel_idx, parameters) end if @@ -808,6 +918,7 @@ subroutine fget_ef_atomic_local_gradient_kernels_fchl(x1, x2, verbose, n1, n2, n kernels = kernels/(2*df) + deallocate (ktmp) deallocate (self_scalar1) deallocate (self_scalar2_ef) deallocate (ksi1) diff --git a/src/qmllib/representations/fchl/ffchl_module.f90 b/src/qmllib/representations/fchl/ffchl_module.f90 index f894ab98..0fb63a62 100644 --- a/src/qmllib/representations/fchl/ffchl_module.f90 +++ b/src/qmllib/representations/fchl/ffchl_module.f90 @@ -2,29 +2,6 @@ module ffchl_module implicit none - ! For polarization - ! double precision, parameter, dimension(19, 3) :: field_displacement = reshape( & - ! & (/ -1, -1, 0, & ! 1 x- y- - ! & -1, 0, -1, & ! 2 x- z- - ! & -1, 0, 0, & ! 3 x- - ! & -1, 0, 1, & ! 4 x- z+ - ! & -1, 1, 0, & ! 5 x- y+ - ! & 0, -1, -1, & ! 6 y- z- - ! & 0, -1, 0, & ! 7 y- - ! & 0, -1, 1, & ! 8 y- z+ - ! & 0, 0, -1, & ! 9 z- - ! & 0, 0, 0, & ! 10 unperturbed = 10 - ! & 0, 0, 1, & ! 11 z+ - ! & 0, 1, -1, & ! 12 y+ z- - ! & 0, 1, 0, & ! 13 y+ - ! & 0, 1, 1, & ! 14 y+ z+ - ! & 1, -1, 0, & ! 15 x+ y- - ! & 1, 0, -1, & ! 16 x+ z- - ! & 1, 0, 0, & ! 17 x+ - ! & 1, 0, 1, & ! 18 x+ z+ - ! & 1, 1, 0 & ! 19 x+ y+ - ! &/), (/ 19, 3 /)) - contains function get_pmax(x, na) result(pmax) @@ -700,7 +677,7 @@ pure function scalar_alchemy(X1, X2, N1, N2, ksi1, ksi2, sin1, sin2, cos1, cos2, deallocate (mask2) end function scalar_alchemy - function get_ksi(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose) result(ksi) + subroutine get_ksi(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose, ksi)! result(ksi) implicit none @@ -724,7 +701,8 @@ function get_ksi(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose logical, intent(in) :: verbose ! Pre-computed two-body weights - double precision, allocatable, dimension(:, :, :) :: ksi + !double precision, allocatable, dimension(:, :, :) :: ksi + double precision, dimension(:, :, :) :: ksi ! Internal counters integer :: maxneigh, maxatoms, nm, a, ni, i @@ -733,7 +711,7 @@ function get_ksi(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose maxatoms = maxval(na) nm = size(x, dim=1) - allocate (ksi(nm, maxatoms, maxneigh)) + !allocate (ksi(nm, maxatoms, maxneigh)) ksi = 0.0d0 @@ -747,9 +725,11 @@ function get_ksi(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose end do !$OMP END PARALLEL do - end function get_ksi + !end function get_ksi + end subroutine get_ksi - function get_ksi_displaced(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose) result(ksi) + !function get_ksi_displaced(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose) result(ksi) + subroutine get_ksi_displaced(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose, ksi) implicit none @@ -773,7 +753,8 @@ function get_ksi_displaced(x, na, nneigh, two_body_power, cut_start, cut_distanc logical, intent(in) :: verbose ! Pre-computed two-body weights - double precision, allocatable, dimension(:, :, :, :, :, :) :: ksi + !double precision, allocatable, dimension(:, :, :, :, :, :) :: ksi + double precision, dimension(:, :, :, :, :, :), intent(out) :: ksi ! Internal counters integer :: maxneigh, maxatoms, nm, a, ni, i, j, pm, xyz, ndisp @@ -783,7 +764,7 @@ function get_ksi_displaced(x, na, nneigh, two_body_power, cut_start, cut_distanc nm = size(x, dim=1) ndisp = size(x, dim=3) - allocate (ksi(nm, 3, ndisp, maxatoms, maxatoms, maxneigh)) + !allocate (ksi(nm, 3, ndisp, maxatoms, maxatoms, maxneigh)) ksi = 0.0d0 @@ -804,9 +785,11 @@ function get_ksi_displaced(x, na, nneigh, two_body_power, cut_start, cut_distanc end do !$OMP END PARALLEL do - end function get_ksi_displaced + !end function get_ksi_displaced + end subroutine get_ksi_displaced - function get_ksi_atomic(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose) result(ksi) + !function get_ksi_atomic(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose) result(ksi) + subroutine get_ksi_atomic(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose, ksi) implicit none @@ -830,7 +813,8 @@ function get_ksi_atomic(x, na, nneigh, two_body_power, cut_start, cut_distance, logical, intent(in) :: verbose ! Pre-computed two-body weights - double precision, allocatable, dimension(:, :) :: ksi + ! double precision, allocatable, dimension(:, :) :: ksi + double precision, dimension(:, :), intent(out) :: ksi ! Internal counters integer :: maxneigh, nm, i @@ -838,7 +822,7 @@ function get_ksi_atomic(x, na, nneigh, two_body_power, cut_start, cut_distance, maxneigh = maxval(nneigh) nm = size(x, dim=1) - allocate (ksi(na, maxneigh)) + ! allocate (ksi(na, maxneigh)) ksi = 0.0d0 @@ -849,7 +833,8 @@ function get_ksi_atomic(x, na, nneigh, two_body_power, cut_start, cut_distance, end do !$OMP END PARALLEL do - end function get_ksi_atomic + !end function get_ksi_atomic + end subroutine get_ksi_atomic subroutine init_cosp_sinp(x, na, nneigh, three_body_power, order, cut_start, cut_distance, cosp, sinp, verbose) @@ -1043,8 +1028,10 @@ subroutine init_cosp_sinp_atomic(x, na, nneigh, three_body_power, order, cut_sta end subroutine init_cosp_sinp_atomic - function get_selfscalar(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) result(self_scalar) + !function get_selfscalar(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) result(self_scalar) + subroutine get_selfscalar(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar) implicit none @@ -1094,12 +1081,13 @@ function get_selfscalar(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d_width, & ! Whether to be verbose with output logical, intent(in) :: verbose - double precision, allocatable, dimension(:, :) :: self_scalar + !double precision, allocatable, dimension(:, :) :: self_scalar + double precision, dimension(:, :) :: self_scalar ! Internal counters integer :: a, ni, i - allocate (self_scalar(nm, maxval(na))) + !allocate (self_scalar(nm, maxval(na))) !$OMP PARALLEL DO PRIVATE(ni) do a = 1, nm @@ -1115,10 +1103,13 @@ function get_selfscalar(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d_width, & end do !$OMP END PARALLEL DO - end function get_selfscalar + !end function get_selfscalar + end subroutine get_selfscalar - function get_selfscalar_displaced(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) result(self_scalar) + !function get_selfscalar_displaced(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) result(self_scalar) + subroutine get_selfscalar_displaced(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar) implicit none @@ -1168,13 +1159,14 @@ function get_selfscalar_displaced(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d ! Whether to be verbose with output logical, intent(in) :: verbose - double precision, allocatable, dimension(:, :, :, :, :) :: self_scalar + !double precision, allocatable, dimension(:, :, :, :, :) :: self_scalar + double precision, dimension(:, :, :, :, :) :: self_scalar ! Internal counters integer :: a, ni, i, j, pm, xyz, xyz_pm, ndisp ndisp = size(x, dim=3) - allocate (self_scalar(nm, 3, ndisp, maxval(na), maxval(na))) + ! allocate (self_scalar(nm, 3, ndisp, maxval(na), maxval(na))) self_scalar = 0.0d0 !$OMP PARALLEL DO PRIVATE(ni, xyz_pm) schedule(dynamic) @@ -1202,751 +1194,7 @@ function get_selfscalar_displaced(x, nm, na, nneigh, ksi, sinp, cosp, t_width, d end do !$OMP END PARALLEL do - end function get_selfscalar_displaced - - function get_ksi_ef(x, na, nneigh, two_body_power, cut_start, cut_distance, ef_scale, df, verbose) result(ksi) - - implicit none - - ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) - double precision, dimension(:, :, :, :), intent(in) :: x - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: na - - ! Number of neighbors for each atom in each compound - integer, dimension(:, :), intent(in) :: nneigh - - ! Decaying powerlaws for two-body term - double precision, intent(in) :: two_body_power - - ! Fraction of cut_distance at which cut-off starts - double precision, intent(in) :: cut_start - double precision, intent(in) :: cut_distance - - ! Electric field displacement - double precision, intent(in) :: ef_scale - double precision, intent(in) :: df - - ! Display output - logical, intent(in) :: verbose - - ! Pre-computed two-body weights - double precision, allocatable, dimension(:, :, :, :, :) :: ksi - - ! Internal counters - integer :: maxneigh, maxatoms, nm, a, ni, i, xyz, pm - - ! Electric field - double precision, dimension(3) :: field - - maxneigh = maxval(nneigh) - maxatoms = maxval(na) - nm = size(x, dim=1) - - allocate (ksi(nm, 3, 2, maxatoms, maxneigh)) - - ksi = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(ni, field) - do a = 1, nm - ni = na(a) - do xyz = 1, 3 - do pm = 1, 2 - - field = 0.0d0 - field(xyz) = (pm - 1.5d0)*2.0d0*df - - ! write(*,*) xyz, pm, field - - do i = 1, ni - - ksi(a, xyz, pm, i, :) = get_twobody_weights_ef(x(a, i, :, :), field, nneigh(a, i), & - & two_body_power, cut_start, cut_distance, maxneigh, ef_scale) - - end do - end do - end do - end do - !$OMP END PARALLEL do - - end function get_ksi_ef - - function get_ksi_ef_field(x, na, nneigh, two_body_power, cut_start, cut_distance, fields, ef_scale, verbose) result(ksi) - - implicit none - - ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) - double precision, dimension(:, :, :, :), intent(in) :: x - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: na - - ! Number of neighbors for each atom in each compound - integer, dimension(:, :), intent(in) :: nneigh - - ! Decaying powerlaws for two-body term - double precision, intent(in) :: two_body_power - - ! Fraction of cut_distance at which cut-off starts - double precision, intent(in) :: cut_start - double precision, intent(in) :: cut_distance - - ! Electric fields for each representation - double precision, dimension(:, :), intent(in) :: fields - - ! Display output - logical, intent(in) :: verbose - - ! Pre-computed two-body weights - double precision, allocatable, dimension(:, :, :) :: ksi - - ! Internal counters - integer :: maxneigh, maxatoms, nm, a, ni, i - - double precision, intent(in) :: ef_scale - - ! Electric field displacement - ! double precision, intent(in) :: df - - maxneigh = maxval(nneigh) - maxatoms = maxval(na) - nm = size(x, dim=1) - - allocate (ksi(nm, maxatoms, maxneigh)) - - ksi = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(ni) - do a = 1, nm - ni = na(a) - - do i = 1, ni - - ksi(a, i, :) = get_twobody_weights_ef(x(a, i, :, :), fields(a, :), nneigh(a, i), & - & two_body_power, cut_start, cut_distance, maxneigh, ef_scale) - end do - - end do - !$OMP END PARALLEL do - - end function get_ksi_ef_field - - subroutine init_cosp_sinp_ef(x, na, nneigh, three_body_power, order, cut_start, cut_distance, & - & cosp, sinp, ef_scale, df, verbose) - - implicit none - - ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) - double precision, dimension(:, :, :, :), intent(in) :: x - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: na - - ! Number of neighbors for each atom in each compound - integer, dimension(:, :), intent(in) :: nneigh - - ! Decaying powerlaws for two-body term - double precision, intent(in) :: three_body_power - - integer, intent(in) :: order - - ! Fraction of cut_distance at which cut-off starts - double precision, intent(in) :: cut_start - double precision, intent(in) :: cut_distance - double precision, intent(in) :: ef_scale - - ! Electric field displacement - double precision, intent(in) :: df - - ! Cosine and sine terms for each atomtype - double precision, dimension(:, :, :, :, :, :, :), intent(out) :: cosp - double precision, dimension(:, :, :, :, :, :, :), intent(out) :: sinp - - ! Display output - logical, intent(in) :: verbose - - ! Internal counters - integer :: maxneigh, maxatoms, pmax, nm, a, ni, i - - double precision, allocatable, dimension(:, :, :, :) :: fourier - - ! Internal counters - integer :: xyz, pm - - ! Electric field - double precision, dimension(3) :: field - - maxneigh = maxval(nneigh) - maxatoms = maxval(na) - nm = size(x, dim=1) - - pmax = get_pmax(x, na) - - cosp = 0.0d0 - sinp = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(ni, fourier, field) schedule(dynamic) - do a = 1, nm - ni = na(a) - do xyz = 1, 3 - do pm = 1, 2 - - field = 0.0d0 - field(xyz) = (pm - 1.5d0)*2.0d0*df - - do i = 1, ni - - fourier = get_threebody_fourier_ef(x(a, i, :, :), field, & - & nneigh(a, i), order, three_body_power, cut_start, cut_distance, pmax, order, maxneigh, ef_scale) - - cosp(a, xyz, pm, i, :, :, :) = fourier(1, :, :, :) - sinp(a, xyz, pm, i, :, :, :) = fourier(2, :, :, :) - - end do - end do - end do - end do - !$OMP END PARALLEL DO - - end subroutine init_cosp_sinp_ef - - subroutine init_cosp_sinp_ef_field(x, na, nneigh, three_body_power, & - & order, cut_start, cut_distance, cosp, sinp, fields, ef_scale, verbose) - - implicit none - - ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) - double precision, dimension(:, :, :, :), intent(in) :: x - - ! List of numbers of atoms in each molecule - integer, dimension(:), intent(in) :: na - - ! Number of neighbors for each atom in each compound - integer, dimension(:, :), intent(in) :: nneigh - - ! Decaying powerlaws for two-body term - double precision, intent(in) :: three_body_power - - integer, intent(in) :: order - - ! Fraction of cut_distance at which cut-off starts - double precision, intent(in) :: cut_start - double precision, intent(in) :: cut_distance - double precision, intent(in) :: ef_scale - - ! Electric field displacement - ! double precision, intent(in) :: df - - ! Display output - logical, intent(in) :: verbose - - ! Cosine and sine terms for each atomtype - double precision, dimension(:, :, :, :, :), intent(out) :: cosp - double precision, dimension(:, :, :, :, :), intent(out) :: sinp - - ! Electric fields for each representation - double precision, dimension(:, :), intent(in) :: fields - - ! Internal counters - integer :: maxneigh, maxatoms, pmax, nm, a, ni, i - - double precision, allocatable, dimension(:, :, :, :) :: fourier - - maxneigh = maxval(nneigh) - maxatoms = maxval(na) - nm = size(x, dim=1) - - pmax = get_pmax(x, na) - - cosp = 0.0d0 - sinp = 0.0d0 - - !$OMP PARALLEL DO PRIVATE(ni, fourier) schedule(dynamic) - do a = 1, nm - - ni = na(a) - - do i = 1, ni - - fourier = get_threebody_fourier_ef(x(a, i, :, :), fields(a, :), & - & nneigh(a, i), order, three_body_power, cut_start, cut_distance, pmax, order, maxneigh, ef_scale) - - cosp(a, i, :, :, :) = fourier(1, :, :, :) - sinp(a, i, :, :, :) = fourier(2, :, :, :) - end do - - end do - !$OMP END PARALLEL DO - - end subroutine init_cosp_sinp_ef_field - -! Calculate the Fourier terms for the FCHL three-body expansion - function get_threebody_fourier_ef(x, field, neighbors, order, power, cut_start, cut_distance, & - & dim1, dim2, dim3, ef_scale) result(fourier) - - implicit none - - ! Input representation, dimension=(5,n). - double precision, dimension(:, :), intent(in) :: x - - double precision, dimension(3), intent(in) :: field - - ! Number of neighboring atoms to iterate over. - integer, intent(in) :: neighbors - - ! Fourier-expansion order. - integer, intent(in) :: order - - ! Power law - double precision, intent(in) :: power - - ! Lower limit of damping function - double precision, intent(in) :: cut_start - - ! Upper limit of damping function - double precision, intent(in) :: cut_distance - - double precision, intent(in) :: ef_scale - - ! Dimensions or the output array. - integer, intent(in) :: dim1, dim2, dim3 - - ! dim(1,:,:,:) are cos terms, dim(2,:,:,:) are sine terms. - double precision, dimension(2, dim1, dim2, dim3) :: fourier - - ! Pi at double precision. - double precision, parameter :: pi = 4.0d0*atan(1.0d0) - - ! Internal counters. - integer :: j, k, m - - ! Indexes for the periodic-table distance matrix. - integer :: pj, pk - - ! Angle between atoms for the three-body term. - double precision :: theta - - ! Three-body weight - double precision :: ksi3 - - ! Temporary variables for cos and sine Fourier terms. - double precision :: cos_m, sin_m - - fourier = 0.0d0 - - do j = 2, neighbors - do k = j + 1, neighbors - - ksi3 = calc_ksi3_ef(X(:, :), field, j, k, neighbors, power, cut_start, cut_distance, ef_scale) - theta = calc_angle(x(3:5, j), x(3:5, 1), x(3:5, k)) - - pj = int(x(2, k)) - pk = int(x(2, j)) - - do m = 1, order - - cos_m = (cos(m*theta) - cos((theta + pi)*m))*ksi3 - sin_m = (sin(m*theta) - sin((theta + pi)*m))*ksi3 - - fourier(1, pj, m, j) = fourier(1, pj, m, j) + cos_m - fourier(2, pj, m, j) = fourier(2, pj, m, j) + sin_m - - fourier(1, pk, m, k) = fourier(1, pk, m, k) + cos_m - fourier(2, pk, m, k) = fourier(2, pk, m, k) + sin_m - - end do - - end do - end do - - return - - end function get_threebody_fourier_ef - - function calc_ksi3_ef(X, field, j, k, neighbors, power, cut_start, cut_distance, ef_scale) result(ksi3) - - implicit none - - double precision, dimension(:, :), intent(in) :: X - - double precision, dimension(3), intent(in) :: field - - integer, intent(in) :: j - integer, intent(in) :: k - - ! Number of neighboring atoms to iterate over. - integer, intent(in) :: neighbors - - double precision, intent(in) :: power - double precision, intent(in) :: cut_start - double precision, intent(in) :: cut_distance - double precision, intent(in) :: ef_scale - - double precision :: cos_i, cos_j, cos_k - double precision :: di, dj, dk - - double precision :: ksi3 - double precision :: cut - - ! Center of nuclear charge - double precision, dimension(3) :: coz - double precision, dimension(3) :: dipole - - double precision :: total_charge - - integer :: i - ! coz - ! dipole - - cos_i = calc_cos_angle(x(3:5, k), x(3:5, 1), x(3:5, j)) - cos_j = calc_cos_angle(x(3:5, j), x(3:5, k), x(3:5, 1)) - cos_k = calc_cos_angle(x(3:5, 1), x(3:5, j), x(3:5, k)) - - dk = x(1, j) - dj = x(1, k) - di = norm2(x(3:5, j) - x(3:5, k)) - - cut = cut_function(dk, cut_start, cut_distance)* & - & cut_function(dj, cut_start, cut_distance)* & - & cut_function(di, cut_start, cut_distance) - - total_charge = 0.0d0 - - coz = 0.0d0 - - do i = 1, neighbors - coz = coz + x(3:5, i)*x(2, i) - total_charge = total_charge + x(2, i) - end do - - coz = coz/total_charge - - dipole = (x(3:5, 1) - coz)*x(6, 1) + (x(3:5, j) - coz)*x(6, j) & - & + (x(3:5, k) - coz)*x(6, k) - - ! ksi3 = cut * (1.0d0 + 3.0d0 * cos_i*cos_j*cos_k) / (di * dj * dk)**power - - ksi3 = cut*((1.0d0 + 3.0d0*cos_i*cos_j*cos_k)/(di*dj*dk)**power & - & + ef_scale*dot_product(dipole, field)) - - ! ksi3 = cut * dot_product(dipole, field) - - end function calc_ksi3_ef - - function get_twobody_weights_ef(x, field, neighbors, power, cut_start, cut_distance, dim1, ef_scale) result(ksi) - - implicit none - - double precision, dimension(:, :), intent(in) :: x - double precision, dimension(3), intent(in) :: field - integer, intent(in) :: neighbors - double precision, intent(in) :: power - double precision, intent(in) :: cut_start - double precision, intent(in) :: cut_distance - double precision, intent(in) :: ef_scale - integer, intent(in) :: dim1 - - double precision, dimension(dim1) :: ksi - - ! Electric field displacement - ! double precision, intent(in) :: df - - integer :: i - - double precision, dimension(3) :: dipole - - ! Center of nuclear charge - double precision, dimension(3) :: coz - double precision :: total_charge - - ksi = 0.0d0 - - coz = 0.0d0 - total_charge = 0.0d0 - do i = 1, neighbors - coz = coz + x(3:5, i)*x(2, i) - total_charge = total_charge + x(2, i) - end do - - coz = coz/total_charge - - do i = 2, neighbors - - dipole = (x(3:5, 1) - coz)*x(6, 1) + (x(3:5, i) - coz)*x(6, i) - - ksi(i) = cut_function(x(1, i), cut_start, cut_distance) & - & *(1.0d0/x(1, i)**power + ef_scale*dot_product(dipole, field)) - - ! ksi(i) = cut_function(x(1, i), cut_start, cut_distance) & - ! & * (dot_product(dipole, field)) - - end do - - end function get_twobody_weights_ef - -! Code for polarization -! function get_twobody_weights_pol(x, field, neighbors, power, cut_start, cut_distance, dim1, ef_scale, qidx, df) result(ksi) -! -! implicit none -! -! double precision, dimension(:,:), intent(in) :: x -! double precision, dimension(3), intent(in) :: field -! integer, intent(in) :: neighbors -! double precision, intent(in) :: power -! double precision, intent(in) :: cut_start -! double precision, intent(in) :: cut_distance -! double precision, intent(in) :: ef_scale -! integer, intent(in) :: dim1 -! -! integer, intent(in) :: qidx -! -! ! Electric field displacement -! double precision, intent(in) :: df -! -! double precision, dimension(dim1) :: ksi -! -! -! integer :: i -! -! double precision, dimension(3) :: dipole -! -! ! Center of nuclear charge -! double precision, dimension(3) :: coz -! double precision :: total_charge -! -! ksi = 0.0d0 -! -! coz = 0.0d0 -! total_charge = 0.0d0 -! do i = 1, neighbors -! coz = coz + x(3:5, i) * x(2, i) -! total_charge = total_charge + x(2, i) -! enddo -! -! coz = coz / total_charge -! -! write(*,*) "SHAPE", size(x, dim=1), size(x, dim=2) -! do i = 2, neighbors -! -! write(*,*) i, qidx, x(qidx + 5, 1), x(qidx + 5, i), sum(x(qidx + 5, :)) -! -! dipole = (x(3:5, 1) - coz) * x(qidx + 5, 1) + (x(3:5, i) - coz) * x(qidx + 5, i) -! -! ! ksi(i) = cut_function(x(1, i), cut_start, cut_distance) & -! ! & * (1.0d0/ x(1, i)**power + ef_scale * dot_product(dipole, field)) -! -! ksi(i) = cut_function(x(1, i), cut_start, cut_distance) & -! & * (dot_product(dipole, field)) -! -! enddo -! -! end function get_twobody_weights_pol -! -! function get_ksi_pol(x, na, nneigh, two_body_power, cut_start, cut_distance, ef_scale, df) result(ksi) -! -! use omp_lib, only: omp_get_wtime -! -! implicit none -! -! ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) -! double precision, dimension(:,:,:,:), intent(in) :: x -! -! ! List of numbers of atoms in each molecule -! integer, dimension(:), intent(in) :: na -! -! ! Number of neighbors for each atom in each compound -! integer, dimension(:,:), intent(in) :: nneigh -! -! ! Decaying powerlaws for two-body term -! double precision, intent(in) :: two_body_power -! -! ! Fraction of cut_distance at which cut-off starts -! double precision, intent(in) :: cut_start -! double precision, intent(in) :: cut_distance -! -! ! Electric field displacement -! double precision, intent(in) :: ef_scale -! double precision, intent(in) :: df -! -! ! Pre-computed two-body weights -! double precision, allocatable, dimension(:,:,:,:) :: ksi -! -! ! Internal counters -! integer :: maxneigh, maxatoms, nm, a, ni, i, xyz -! -! ! Electric field -! double precision, dimension(3) :: field -! -! double precision :: t_start, t_end -! -! write (*,"(A)", advance="no") "TWO-BODY TERMS POLARIZABLE" -! t_start = omp_get_wtime() -! -! maxneigh = maxval(nneigh) -! maxatoms = maxval(na) -! nm = size(x, dim=1) -! -! allocate(ksi(nm, 19, maxatoms, maxneigh)) -! -! ksi = 0.0d0 -! -! ! !$OMP PARALLEL DO PRIVATE(ni, field) -! do a = 1, nm -! ni = na(a) -! do xyz = 1, 19 -! -! field(1) = field_displacement(xyz,1) * df -! field(2) = field_displacement(xyz,2) * df -! field(3) = field_displacement(xyz,3) * df -! -! do i = 1, ni -! -! ksi(a, xyz, i, :) = get_twobody_weights_pol(x(a,i,:,:), field, nneigh(a, i), & -! & two_body_power, cut_start, cut_distance, maxneigh, ef_scale, xyz, df) -! -! enddo -! enddo -! enddo -! ! !$OMP END PARALLEL do -! -! t_end = omp_get_wtime() -! write (*,"(A,F12.4,A)") " Time = ", t_end - t_start, " s" -! -! end function get_ksi_pol -! -! subroutine init_cosp_sinp_pol(x, na, nneigh, three_body_power, order, cut_start, cut_distance, cosp, sinp, ef_scale, df) -! -! use omp_lib, only: omp_get_wtime -! -! implicit none -! -! ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) -! double precision, dimension(:,:,:,:), intent(in) :: x -! -! ! List of numbers of atoms in each molecule -! integer, dimension(:), intent(in) :: na -! -! ! Number of neighbors for each atom in each compound -! integer, dimension(:,:), intent(in) :: nneigh -! -! ! Decaying powerlaws for two-body term -! double precision, intent(in) :: three_body_power -! -! integer, intent(in) :: order -! -! ! Fraction of cut_distance at which cut-off starts -! double precision, intent(in) :: cut_start -! double precision, intent(in) :: cut_distance -! double precision, intent(in) :: ef_scale -! -! ! Electric field displacement -! double precision, intent(in) :: df -! -! ! Cosine and sine terms for each atomtype -! double precision, dimension(:,:,:,:,:,:), intent(out) :: cosp -! double precision, dimension(:,:,:,:,:,:), intent(out) :: sinp -! -! -! ! Internal counters -! integer :: maxneigh, maxatoms, pmax, nm, a, ni, i -! -! double precision, allocatable, dimension(:,:,:,:) :: fourier -! -! double precision :: t_start, t_end -! -! ! Internal counters -! integer :: xyz -! -! ! Electric field -! double precision, dimension(3) :: field -! -! write (*,"(A)", advance="no") "THREE-BODY TERMS POLARIZABLE" -! t_start = omp_get_wtime() -! -! maxneigh = maxval(nneigh) -! maxatoms = maxval(na) -! nm = size(x, dim=1) -! -! pmax = get_pmax(x, na) -! -! cosp = 0.0d0 -! sinp = 0.0d0 -! -! ! !$OMP PARALLEL DO PRIVATE(ni, fourier, field) schedule(dynamic) -! do a = 1, nm -! ni = na(a) -! do xyz = 1, 19 -! -! field(1) = field_displacement(xyz,1) * df -! field(2) = field_displacement(xyz,2) * df -! field(3) = field_displacement(xyz,3) * df -! ! write(*,*) xyz, field_displacement(xyz,:), field -! -! do i = 1, ni -! -! fourier = get_threebody_fourier_ef(x(a,i,:,:), field, & -! & nneigh(a, i), order, three_body_power, cut_start, cut_distance, pmax, order, maxneigh, ef_scale) -! -! cosp(a,xyz,i,:,:,:) = fourier(1,:,:,:) -! sinp(a,xyz,i,:,:,:) = fourier(2,:,:,:) -! -! enddo -! enddo -! enddo -! ! !$OMP END PARALLEL DO -! -! t_end = omp_get_wtime() -! write (*,"(A,F12.4,A)") " Time = ", t_end - t_start, " s" -! -! end subroutine init_cosp_sinp_pol - - function get_twobody_weights_ef_field(x, field, neighbors, power, cut_start, cut_distance, dim1, ef_scale) result(ksi) - - implicit none - - double precision, dimension(:, :), intent(in) :: x - double precision, dimension(3), intent(in) :: field - integer, intent(in) :: neighbors - double precision, intent(in) :: power - double precision, intent(in) :: cut_start - double precision, intent(in) :: cut_distance - double precision, intent(in) :: ef_scale - integer, intent(in) :: dim1 - - double precision, dimension(dim1) :: ksi - - ! Electric field displacement - ! double precision, intent(in) :: df - - integer :: i - - double precision, dimension(3) :: dipole - - ! Center of nuclear charge - double precision, dimension(3) :: coz - double precision :: total_charge - - ksi = 0.0d0 - - coz = 0.0d0 - total_charge = 0.0d0 - do i = 1, neighbors - coz = coz + x(3:5, i)*x(2, i) - total_charge = total_charge + x(2, i) - end do - - coz = coz/total_charge - - do i = 2, neighbors - - dipole = (x(3:5, 1) - coz)*x(6, 1) + (x(3:5, i) - coz)*x(6, i) - - ksi(i) = cut_function(x(1, i), cut_start, cut_distance) & - & *(1.0d0/x(1, i)**power + ef_scale*dot_product(dipole, field)) - - end do - - end function get_twobody_weights_ef_field + !end function get_selfscalar_displaced + end subroutine get_selfscalar_displaced end module ffchl_module diff --git a/src/qmllib/representations/fchl/ffchl_module_ef.f90 b/src/qmllib/representations/fchl/ffchl_module_ef.f90 new file mode 100644 index 00000000..775a86b4 --- /dev/null +++ b/src/qmllib/representations/fchl/ffchl_module_ef.f90 @@ -0,0 +1,554 @@ +module ffchl_module_ef + + implicit none + +contains + + !function get_ksi_ef(x, na, nneigh, two_body_power, cut_start, cut_distance, ef_scale, df, verbose) result(ksi) + subroutine get_ksi_ef(x, na, nneigh, two_body_power, cut_start, cut_distance, ef_scale, df, verbose, ksi) + + implicit none + + ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) + double precision, dimension(:, :, :, :), intent(in) :: x + + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: na + + ! Number of neighbors for each atom in each compound + integer, dimension(:, :), intent(in) :: nneigh + + ! Decaying powerlaws for two-body term + double precision, intent(in) :: two_body_power + + ! Fraction of cut_distance at which cut-off starts + double precision, intent(in) :: cut_start + double precision, intent(in) :: cut_distance + + ! Electric field displacement + double precision, intent(in) :: ef_scale + double precision, intent(in) :: df + + ! Display output + logical, intent(in) :: verbose + + ! Pre-computed two-body weights + ! double precision, allocatable, dimension(:, :, :, :, :) :: ksi + double precision, dimension(:, :, :, :, :) :: ksi + + ! Internal counters + integer :: maxneigh, maxatoms, nm, a, ni, i, xyz, pm + + ! Electric field + double precision, dimension(3) :: field + + maxneigh = maxval(nneigh) + maxatoms = maxval(na) + nm = size(x, dim=1) + + !allocate (ksi(nm, 3, 2, maxatoms, maxneigh)) + + ksi = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(ni, field) + do a = 1, nm + ni = na(a) + do xyz = 1, 3 + do pm = 1, 2 + + field = 0.0d0 + field(xyz) = (pm - 1.5d0)*2.0d0*df + + ! write(*,*) xyz, pm, field + + do i = 1, ni + + ksi(a, xyz, pm, i, :) = get_twobody_weights_ef(x(a, i, :, :), field, nneigh(a, i), & + & two_body_power, cut_start, cut_distance, maxneigh, ef_scale) + + end do + end do + end do + end do + !$OMP END PARALLEL do + + ! end function get_ksi_ef + end subroutine get_ksi_ef + + !function get_ksi_ef_field(x, na, nneigh, two_body_power, cut_start, cut_distance, fields, ef_scale, verbose) result(ksi) + subroutine get_ksi_ef_field(x, na, nneigh, two_body_power, cut_start, cut_distance, fields, ef_scale, verbose, ksi) + + implicit none + + ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) + double precision, dimension(:, :, :, :), intent(in) :: x + + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: na + + ! Number of neighbors for each atom in each compound + integer, dimension(:, :), intent(in) :: nneigh + + ! Decaying powerlaws for two-body term + double precision, intent(in) :: two_body_power + + ! Fraction of cut_distance at which cut-off starts + double precision, intent(in) :: cut_start + double precision, intent(in) :: cut_distance + + ! Electric fields for each representation + double precision, dimension(:, :), intent(in) :: fields + + ! Display output + logical, intent(in) :: verbose + + ! Pre-computed two-body weights + !double precision, allocatable, dimension(:, :, :) :: ksi + double precision, dimension(:, :, :), intent(out) :: ksi + + ! Internal counters + integer :: maxneigh, maxatoms, nm, a, ni, i + + double precision, intent(in) :: ef_scale + + ! Electric field displacement + ! double precision, intent(in) :: df + + maxneigh = maxval(nneigh) + maxatoms = maxval(na) + nm = size(x, dim=1) + + !allocate (ksi(nm, maxatoms, maxneigh)) + + ksi = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(ni) + do a = 1, nm + ni = na(a) + + do i = 1, ni + + ksi(a, i, :) = get_twobody_weights_ef(x(a, i, :, :), fields(a, :), nneigh(a, i), & + & two_body_power, cut_start, cut_distance, maxneigh, ef_scale) + end do + + end do + !$OMP END PARALLEL do + + !end function get_ksi_ef_field + end subroutine get_ksi_ef_field + + subroutine init_cosp_sinp_ef(x, na, nneigh, three_body_power, order, cut_start, cut_distance, & + & cosp, sinp, ef_scale, df, verbose) + + use ffchl_module, only: get_pmax + + implicit none + + ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) + double precision, dimension(:, :, :, :), intent(in) :: x + + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: na + + ! Number of neighbors for each atom in each compound + integer, dimension(:, :), intent(in) :: nneigh + + ! Decaying powerlaws for two-body term + double precision, intent(in) :: three_body_power + + integer, intent(in) :: order + + ! Fraction of cut_distance at which cut-off starts + double precision, intent(in) :: cut_start + double precision, intent(in) :: cut_distance + double precision, intent(in) :: ef_scale + + ! Electric field displacement + double precision, intent(in) :: df + + ! Cosine and sine terms for each atomtype + double precision, dimension(:, :, :, :, :, :, :), intent(out) :: cosp + double precision, dimension(:, :, :, :, :, :, :), intent(out) :: sinp + + ! Display output + logical, intent(in) :: verbose + + ! Internal counters + integer :: maxneigh, maxatoms, pmax, nm, a, ni, i + + double precision, allocatable, dimension(:, :, :, :) :: fourier + + ! Internal counters + integer :: xyz, pm + + ! Electric field + double precision, dimension(3) :: field + + maxneigh = maxval(nneigh) + maxatoms = maxval(na) + nm = size(x, dim=1) + + pmax = get_pmax(x, na) + + cosp = 0.0d0 + sinp = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(ni, fourier, field) schedule(dynamic) + do a = 1, nm + ni = na(a) + do xyz = 1, 3 + do pm = 1, 2 + + field = 0.0d0 + field(xyz) = (pm - 1.5d0)*2.0d0*df + + do i = 1, ni + + fourier = get_threebody_fourier_ef(x(a, i, :, :), field, & + & nneigh(a, i), order, three_body_power, cut_start, cut_distance, pmax, order, maxneigh, ef_scale) + + cosp(a, xyz, pm, i, :, :, :) = fourier(1, :, :, :) + sinp(a, xyz, pm, i, :, :, :) = fourier(2, :, :, :) + + end do + end do + end do + end do + !$OMP END PARALLEL DO + + end subroutine init_cosp_sinp_ef + + subroutine init_cosp_sinp_ef_field(x, na, nneigh, three_body_power, & + & order, cut_start, cut_distance, cosp, sinp, fields, ef_scale, verbose) + + use ffchl_module, only: get_pmax + + implicit none + + ! FCHL descriptors for the training set, format (i,maxatoms,5,maxneighbors) + double precision, dimension(:, :, :, :), intent(in) :: x + + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: na + + ! Number of neighbors for each atom in each compound + integer, dimension(:, :), intent(in) :: nneigh + + ! Decaying powerlaws for two-body term + double precision, intent(in) :: three_body_power + + integer, intent(in) :: order + + ! Fraction of cut_distance at which cut-off starts + double precision, intent(in) :: cut_start + double precision, intent(in) :: cut_distance + double precision, intent(in) :: ef_scale + + ! Electric field displacement + ! double precision, intent(in) :: df + + ! Display output + logical, intent(in) :: verbose + + ! Cosine and sine terms for each atomtype + double precision, dimension(:, :, :, :, :), intent(out) :: cosp + double precision, dimension(:, :, :, :, :), intent(out) :: sinp + + ! Electric fields for each representation + double precision, dimension(:, :), intent(in) :: fields + + ! Internal counters + integer :: maxneigh, maxatoms, pmax, nm, a, ni, i + + double precision, allocatable, dimension(:, :, :, :) :: fourier + + maxneigh = maxval(nneigh) + maxatoms = maxval(na) + nm = size(x, dim=1) + + pmax = get_pmax(x, na) + + cosp = 0.0d0 + sinp = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(ni, fourier) schedule(dynamic) + do a = 1, nm + + ni = na(a) + + do i = 1, ni + + fourier = get_threebody_fourier_ef(x(a, i, :, :), fields(a, :), & + & nneigh(a, i), order, three_body_power, cut_start, cut_distance, pmax, order, maxneigh, ef_scale) + + cosp(a, i, :, :, :) = fourier(1, :, :, :) + sinp(a, i, :, :, :) = fourier(2, :, :, :) + end do + + end do + !$OMP END PARALLEL DO + + end subroutine init_cosp_sinp_ef_field + + ! Calculate the Fourier terms for the FCHL three-body expansion + function get_threebody_fourier_ef(x, field, neighbors, order, power, cut_start, cut_distance, & + & dim1, dim2, dim3, ef_scale) result(fourier) + + use ffchl_module, only: calc_angle + + implicit none + + ! Input representation, dimension=(5,n). + double precision, dimension(:, :), intent(in) :: x + + double precision, dimension(3), intent(in) :: field + + ! Number of neighboring atoms to iterate over. + integer, intent(in) :: neighbors + + ! Fourier-expansion order. + integer, intent(in) :: order + + ! Power law + double precision, intent(in) :: power + + ! Lower limit of damping function + double precision, intent(in) :: cut_start + + ! Upper limit of damping function + double precision, intent(in) :: cut_distance + + double precision, intent(in) :: ef_scale + + ! Dimensions or the output array. + integer, intent(in) :: dim1, dim2, dim3 + + ! dim(1,:,:,:) are cos terms, dim(2,:,:,:) are sine terms. + double precision, dimension(2, dim1, dim2, dim3) :: fourier + + ! Pi at double precision. + double precision, parameter :: pi = 4.0d0*atan(1.0d0) + + ! Internal counters. + integer :: j, k, m + + ! Indexes for the periodic-table distance matrix. + integer :: pj, pk + + ! Angle between atoms for the three-body term. + double precision :: theta + + ! Three-body weight + double precision :: ksi3 + + ! Temporary variables for cos and sine Fourier terms. + double precision :: cos_m, sin_m + + fourier = 0.0d0 + + do j = 2, neighbors + do k = j + 1, neighbors + + ksi3 = calc_ksi3_ef(X(:, :), field, j, k, neighbors, power, cut_start, cut_distance, ef_scale) + theta = calc_angle(x(3:5, j), x(3:5, 1), x(3:5, k)) + + pj = int(x(2, k)) + pk = int(x(2, j)) + + do m = 1, order + + cos_m = (cos(m*theta) - cos((theta + pi)*m))*ksi3 + sin_m = (sin(m*theta) - sin((theta + pi)*m))*ksi3 + + fourier(1, pj, m, j) = fourier(1, pj, m, j) + cos_m + fourier(2, pj, m, j) = fourier(2, pj, m, j) + sin_m + + fourier(1, pk, m, k) = fourier(1, pk, m, k) + cos_m + fourier(2, pk, m, k) = fourier(2, pk, m, k) + sin_m + + end do + + end do + end do + + return + + end function get_threebody_fourier_ef + + function calc_ksi3_ef(X, field, j, k, neighbors, power, cut_start, cut_distance, ef_scale) result(ksi3) + + use ffchl_module, only: cut_function, calc_cos_angle + + implicit none + + double precision, dimension(:, :), intent(in) :: X + + double precision, dimension(3), intent(in) :: field + + integer, intent(in) :: j + integer, intent(in) :: k + + ! Number of neighboring atoms to iterate over. + integer, intent(in) :: neighbors + + double precision, intent(in) :: power + double precision, intent(in) :: cut_start + double precision, intent(in) :: cut_distance + double precision, intent(in) :: ef_scale + + double precision :: cos_i, cos_j, cos_k + double precision :: di, dj, dk + + double precision :: ksi3 + double precision :: cut + + ! Center of nuclear charge + double precision, dimension(3) :: coz + double precision, dimension(3) :: dipole + + double precision :: total_charge + + integer :: i + ! coz + ! dipole + + cos_i = calc_cos_angle(x(3:5, k), x(3:5, 1), x(3:5, j)) + cos_j = calc_cos_angle(x(3:5, j), x(3:5, k), x(3:5, 1)) + cos_k = calc_cos_angle(x(3:5, 1), x(3:5, j), x(3:5, k)) + + dk = x(1, j) + dj = x(1, k) + di = norm2(x(3:5, j) - x(3:5, k)) + + cut = cut_function(dk, cut_start, cut_distance)* & + & cut_function(dj, cut_start, cut_distance)* & + & cut_function(di, cut_start, cut_distance) + + total_charge = 0.0d0 + + coz = 0.0d0 + + do i = 1, neighbors + coz = coz + x(3:5, i)*x(2, i) + total_charge = total_charge + x(2, i) + end do + + coz = coz/total_charge + + dipole = (x(3:5, 1) - coz)*x(6, 1) + (x(3:5, j) - coz)*x(6, j) & + & + (x(3:5, k) - coz)*x(6, k) + + ! ksi3 = cut * (1.0d0 + 3.0d0 * cos_i*cos_j*cos_k) / (di * dj * dk)**power + + ksi3 = cut*((1.0d0 + 3.0d0*cos_i*cos_j*cos_k)/(di*dj*dk)**power & + & + ef_scale*dot_product(dipole, field)) + + ! ksi3 = cut * dot_product(dipole, field) + + end function calc_ksi3_ef + + function get_twobody_weights_ef(x, field, neighbors, power, cut_start, cut_distance, dim1, ef_scale) result(ksi) + + use ffchl_module, only: cut_function + + implicit none + + double precision, dimension(:, :), intent(in) :: x + double precision, dimension(3), intent(in) :: field + integer, intent(in) :: neighbors + double precision, intent(in) :: power + double precision, intent(in) :: cut_start + double precision, intent(in) :: cut_distance + double precision, intent(in) :: ef_scale + integer, intent(in) :: dim1 + + double precision, dimension(dim1) :: ksi + + ! Electric field displacement + ! double precision, intent(in) :: df + + integer :: i + + double precision, dimension(3) :: dipole + + ! Center of nuclear charge + double precision, dimension(3) :: coz + double precision :: total_charge + + ksi = 0.0d0 + + coz = 0.0d0 + total_charge = 0.0d0 + do i = 1, neighbors + coz = coz + x(3:5, i)*x(2, i) + total_charge = total_charge + x(2, i) + end do + + coz = coz/total_charge + + do i = 2, neighbors + + dipole = (x(3:5, 1) - coz)*x(6, 1) + (x(3:5, i) - coz)*x(6, i) + + ksi(i) = cut_function(x(1, i), cut_start, cut_distance) & + & *(1.0d0/x(1, i)**power + ef_scale*dot_product(dipole, field)) + + ! ksi(i) = cut_function(x(1, i), cut_start, cut_distance) & + ! & * (dot_product(dipole, field)) + + end do + + end function get_twobody_weights_ef + + function get_twobody_weights_ef_field(x, field, neighbors, power, cut_start, cut_distance, dim1, ef_scale) result(ksi) + + use ffchl_module, only: cut_function + + implicit none + + double precision, dimension(:, :), intent(in) :: x + double precision, dimension(3), intent(in) :: field + integer, intent(in) :: neighbors + double precision, intent(in) :: power + double precision, intent(in) :: cut_start + double precision, intent(in) :: cut_distance + double precision, intent(in) :: ef_scale + integer, intent(in) :: dim1 + + double precision, dimension(dim1) :: ksi + + ! Electric field displacement + ! double precision, intent(in) :: df + + integer :: i + + double precision, dimension(3) :: dipole + + ! Center of nuclear charge + double precision, dimension(3) :: coz + double precision :: total_charge + + ksi = 0.0d0 + + coz = 0.0d0 + total_charge = 0.0d0 + do i = 1, neighbors + coz = coz + x(3:5, i)*x(2, i) + total_charge = total_charge + x(2, i) + end do + + coz = coz/total_charge + + do i = 2, neighbors + + dipole = (x(3:5, 1) - coz)*x(6, 1) + (x(3:5, i) - coz)*x(6, i) + + ksi(i) = cut_function(x(1, i), cut_start, cut_distance) & + & *(1.0d0/x(1, i)**power + ef_scale*dot_product(dipole, field)) + + end do + + end function get_twobody_weights_ef_field + +end module ffchl_module_ef diff --git a/src/qmllib/representations/fchl/ffchl_scalar_kernels.f90 b/src/qmllib/representations/fchl/ffchl_scalar_kernels.f90 index a354fd2c..d2251f7f 100644 --- a/src/qmllib/representations/fchl/ffchl_scalar_kernels.f90 +++ b/src/qmllib/representations/fchl/ffchl_scalar_kernels.f90 @@ -98,6 +98,10 @@ subroutine fget_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, nm1, nm2 integer :: maxneigh1 integer :: maxneigh2 + ! Work kernel space + integer :: n + double precision, allocatable, dimension(:) :: ktmp + kernels(:, :, :) = 0.0d0 ! Get max number of neighbors @@ -112,8 +116,20 @@ subroutine fget_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, nm1, nm2 pmax2 = get_pmax(x2, n2) ! Get two-body weight function - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + !ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + !ksi2 = get_ksi(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + + ! subroutine get_ksi(x, na, nneigh, two_body_power, cut_start, cut_distance, verbose, ksi)! result(ksi) + ! maxneigh = maxval(nneigh) + ! maxatoms = maxval(na) + ! nm = size(x, dim=1) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), maxval(n2), maxval(nneigh2))) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + + n = size(parameters, dim=1) + allocate (ktmp(n)) ! Allocate three-body Fourier terms allocate (cosp1(nm1, maxval(n1), pmax1, order, maxneigh1)) @@ -132,14 +148,22 @@ subroutine fget_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, nm1, nm2 & cosp2, sinp2, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + !self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose,& + &self_scalar1) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) - - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj) + !self_scalar2 = get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar2(nm2, maxval(n2))) + call get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose,& + &self_scalar2) + + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,ktmp) do b = 1, nm2 nj = n2(b) do a = 1, nm1 @@ -155,9 +179,14 @@ subroutine fget_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, nm1, nm2 & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernels(:, a, b) = kernels(:, a, b) & - & + kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & - & kernel_idx, parameters) + !kernels(:, a, b) = kernels(:, a, b) & + ! & + kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & + ! & kernel_idx, parameters) + + ktmp = 0.0d0 + call kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & + & kernel_idx, parameters, ktmp) + kernels(:, a, b) = kernels(:, a, b) + ktmp end do end do @@ -166,6 +195,7 @@ subroutine fget_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, nm1, nm2 end do !$OMP END PARALLEL DO + deallocate (ktmp) deallocate (self_scalar1) deallocate (self_scalar2) deallocate (ksi1) @@ -250,6 +280,9 @@ subroutine fget_symmetric_kernels_fchl(x1, verbose, n1, nneigh1, nm1, nsigmas, & integer :: maxneigh1 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + kernels(:, :, :) = 0.0d0 ang_norm2 = get_angular_norm2(t_width) @@ -257,7 +290,9 @@ subroutine fget_symmetric_kernels_fchl(x1, verbose, n1, nneigh1, nm1, nsigmas, & maxneigh1 = maxval(nneigh1) pmax1 = get_pmax(x1, n1) - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + !ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) allocate (cosp1(nm1, maxval(n1), pmax1, order, maxval(nneigh1))) allocate (sinp1(nm1, maxval(n1), pmax1, order, maxval(nneigh1))) @@ -265,10 +300,15 @@ subroutine fget_symmetric_kernels_fchl(x1, verbose, n1, nneigh1, nm1, nsigmas, & call init_cosp_sinp(x1, n1, nneigh1, three_body_power, order, cut_start, cut_distance, & & cosp1, sinp1, verbose) - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + !self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + + allocate (ktmp(size(parameters, dim=1))) - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,ktmp) do b = 1, nm1 nj = n1(b) do a = b, nm1 @@ -284,9 +324,14 @@ subroutine fget_symmetric_kernels_fchl(x1, verbose, n1, nneigh1, nm1, nsigmas, & & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernels(:, a, b) = kernels(:, a, b) & - & + kernel(self_scalar1(a, i), self_scalar1(b, j), s12, & - & kernel_idx, parameters) + ktmp(:) = 0.0d0 + call kernel(self_scalar1(a, i), self_scalar1(b, j), s12, & + & kernel_idx, parameters, ktmp) + + kernels(:, a, b) = kernels(:, a, b) + ktmp + !kernels(:, a, b) = kernels(:, a, b) & + ! & + kernel(self_scalar1(a, i), self_scalar1(b, j), s12, & + ! & kernel_idx, parameters) kernels(:, b, a) = kernels(:, a, b) @@ -297,6 +342,7 @@ subroutine fget_symmetric_kernels_fchl(x1, verbose, n1, nneigh1, nm1, nsigmas, & end do !$OMP END PARALLEL DO + deallocate (ktmp) deallocate (self_scalar1) deallocate (ksi1) deallocate (cosp1) @@ -377,13 +423,19 @@ subroutine fget_global_symmetric_kernels_fchl(x1, verbose, n1, nneigh1, nm1, nsi integer :: maxneigh1 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + maxneigh1 = maxval(nneigh1) ang_norm2 = get_angular_norm2(t_width) pmax1 = get_pmax(x1, n1) - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + !ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) allocate (cosp1(nm1, maxval(n1), pmax1, order, maxval(nneigh1))) allocate (sinp1(nm1, maxval(n1), pmax1, order, maxval(nneigh1))) @@ -414,7 +466,7 @@ subroutine fget_global_symmetric_kernels_fchl(x1, verbose, n1, nneigh1, nm1, nsi kernels(:, :, :) = 0.0d0 - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,mol_dist) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,mol_dist,ktmp) do b = 1, nm1 nj = n1(b) do a = b, nm1 @@ -437,8 +489,12 @@ subroutine fget_global_symmetric_kernels_fchl(x1, verbose, n1, nneigh1, nm1, nsi end do end do - kernels(:, a, b) = kernel(self_scalar1(a), self_scalar1(b), mol_dist, & - & kernel_idx, parameters) + ktmp = 0.0d0 + call kernel(self_scalar1(a), self_scalar1(b), mol_dist, & + & kernel_idx, parameters, ktmp) + kernels(:, a, b) = ktmp + !kernels(:, a, b) = kernel(self_scalar1(a), self_scalar1(b), mol_dist, & + ! & kernel_idx, parameters) kernels(:, b, a) = kernels(:, a, b) @@ -446,6 +502,7 @@ subroutine fget_global_symmetric_kernels_fchl(x1, verbose, n1, nneigh1, nm1, nsi end do !$OMP END PARALLEL DO + deallocate (ktmp) deallocate (self_scalar1) deallocate (ksi1) deallocate (cosp1) @@ -539,6 +596,10 @@ subroutine fget_global_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, & integer :: maxneigh1 integer :: maxneigh2 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + maxneigh1 = maxval(nneigh1) maxneigh2 = maxval(nneigh2) @@ -547,8 +608,12 @@ subroutine fget_global_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, & pmax1 = get_pmax(x1, n1) pmax2 = get_pmax(x2, n2) - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), maxval(n2), maxval(nneigh2))) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + !ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + !ksi2 = get_ksi(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) allocate (cosp1(nm1, maxval(n1), pmax1, order, maxval(nneigh1))) allocate (sinp1(nm1, maxval(n1), pmax1, order, maxval(nneigh1))) @@ -604,7 +669,7 @@ subroutine fget_global_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, & kernels(:, :, :) = 0.0d0 - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,mol_dist) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12,ni,nj,mol_dist,ktmp) do b = 1, nm2 nj = n2(b) do a = 1, nm1 @@ -626,13 +691,19 @@ subroutine fget_global_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nneigh2, & end do end do - kernels(:, a, b) = kernel(self_scalar1(a), self_scalar2(b), mol_dist, & - & kernel_idx, parameters) + + ktmp = 0.0d0 + call kernel(self_scalar1(a), self_scalar2(b), mol_dist, & + & kernel_idx, parameters, ktmp) + kernels(:, a, b) = ktmp + !kernels(:, a, b) = kernel(self_scalar1(a), self_scalar2(b), mol_dist, & + ! & kernel_idx, parameters) end do end do !$OMP END PARALLEL DO + deallocate (ktmp) deallocate (self_scalar1) deallocate (self_scalar2) deallocate (ksi1) @@ -722,6 +793,10 @@ subroutine fget_atomic_kernels_fchl(x1, x2, verbose, nneigh1, nneigh2, & integer :: maxneigh1 integer :: maxneigh2 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + maxneigh1 = maxval(nneigh1) maxneigh2 = maxval(nneigh2) @@ -730,8 +805,12 @@ subroutine fget_atomic_kernels_fchl(x1, x2, verbose, nneigh1, nneigh2, & pmax1 = get_pmax_atomic(x1, nneigh1) pmax2 = get_pmax_atomic(x2, nneigh2) - ksi1 = get_ksi_atomic(x1, na1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi_atomic(x2, na2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(na1, maxval(nneigh1))) + allocate (ksi2(na2, maxval(nneigh2))) + call get_ksi_atomic(x1, na1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi_atomic(x2, na2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + !ksi1 = get_ksi_atomic(x1, na1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + !ksi2 = get_ksi_atomic(x2, na2, nneigh2, two_body_power, cut_start, cut_distance, verbose) allocate (cosp1(na1, pmax1, order, maxneigh1)) allocate (sinp1(na1, pmax1, order, maxneigh1)) @@ -775,7 +854,7 @@ subroutine fget_atomic_kernels_fchl(x1, x2, verbose, nneigh1, nneigh2, & kernels(:, :, :) = 0.0d0 - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12, ktmp) do i = 1, na1 do j = 1, na2 @@ -786,13 +865,18 @@ subroutine fget_atomic_kernels_fchl(x1, x2, verbose, nneigh1, nneigh2, & & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernels(:, i, j) = kernel(self_scalar1(i), self_scalar2(j), s12, & - & kernel_idx, parameters) + ktmp = 0.0d0 + call kernel(self_scalar1(i), self_scalar2(j), s12, & + & kernel_idx, parameters, ktmp) + kernels(:, i, j) = ktmp + !kernels(:, i, j) = kernel(self_scalar1(i), self_scalar2(j), s12, & + ! & kernel_idx, parameters) end do end do !$OMP END PARALLEL DO + deallocate (ktmp) deallocate (self_scalar1) deallocate (self_scalar2) deallocate (ksi1) @@ -871,13 +955,19 @@ subroutine fget_atomic_symmetric_kernels_fchl(x1, verbose, nneigh1, na1, nsigmas integer :: maxneigh1 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + maxneigh1 = maxval(nneigh1) ang_norm2 = get_angular_norm2(t_width) pmax1 = get_pmax_atomic(x1, nneigh1) - ksi1 = get_ksi_atomic(x1, na1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(na1, maxval(nneigh1))) + call get_ksi_atomic(x1, na1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + !ksi1 = get_ksi_atomic(x1, na1, nneigh1, two_body_power, cut_start, cut_distance, verbose) allocate (cosp1(na1, pmax1, order, maxneigh1)) allocate (sinp1(na1, pmax1, order, maxneigh1)) @@ -902,7 +992,7 @@ subroutine fget_atomic_symmetric_kernels_fchl(x1, verbose, nneigh1, na1, nsigmas kernels(:, :, :) = 0.0d0 - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(s12, ktmp) do i = 1, na1 do j = i, na1 @@ -913,14 +1003,21 @@ subroutine fget_atomic_symmetric_kernels_fchl(x1, verbose, nneigh1, na1, nsigmas & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernels(:, i, j) = kernel(self_scalar1(i), self_scalar1(j), s12, & - & kernel_idx, parameters) + !kernels(:, i, j) = kernel(self_scalar1(i), self_scalar1(j), s12, & + ! & kernel_idx, parameters) + ktmp = 0.0d0 + call kernel(self_scalar1(i), self_scalar1(j), s12, & + & kernel_idx, parameters, ktmp) + + kernels(:, i, j) = ktmp kernels(:, j, i) = kernels(:, i, j) + end do end do !$OMP END PARALLEL DO + deallocate (ktmp) deallocate (self_scalar1) deallocate (ksi1) deallocate (cosp1) @@ -1023,6 +1120,10 @@ subroutine fget_atomic_local_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nnei integer :: maxneigh1 integer :: maxneigh2 + ! Work kernel + double precision, allocatable, dimension(:) :: ktmp + allocate (ktmp(size(parameters, dim=1))) + maxneigh1 = maxval(nneigh1) maxneigh2 = maxval(nneigh2) @@ -1031,8 +1132,12 @@ subroutine fget_atomic_local_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nnei pmax1 = get_pmax(x1, n1) pmax2 = get_pmax(x2, n2) - ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) - ksi2 = get_ksi(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) + allocate (ksi1(size(x1, dim=1), maxval(n1), maxval(nneigh1))) + allocate (ksi2(size(x2, dim=1), maxval(n2), maxval(nneigh2))) + call get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose, ksi1) + call get_ksi(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose, ksi2) + !ksi1 = get_ksi(x1, n1, nneigh1, two_body_power, cut_start, cut_distance, verbose) + !ksi2 = get_ksi(x2, n2, nneigh2, two_body_power, cut_start, cut_distance, verbose) allocate (cosp1(nm1, maxval(n1), pmax1, order, maxval(nneigh1))) allocate (sinp1(nm1, maxval(n1), pmax1, order, maxval(nneigh1))) @@ -1047,16 +1152,22 @@ subroutine fget_atomic_local_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nnei & cosp2, sinp2, verbose) ! Pre-calculate self-scalar terms - self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + allocate (self_scalar1(nm1, maxval(n1))) + allocate (self_scalar2(nm2, maxval(n2))) + call get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar1) + call get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & + & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose, self_scalar2) + !self_scalar1 = get_selfscalar(x1, nm1, n1, nneigh1, ksi1, sinp1, cosp1, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) ! Pre-calculate self-scalar terms - self_scalar2 = get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & - & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) + !self_scalar2 = get_selfscalar(x2, nm2, n2, nneigh2, ksi2, sinp2, cosp2, t_width, d_width, & + ! & cut_distance, order, pd, ang_norm2, distance_scale, angular_scale, alchemy, verbose) kernels(:, :, :) = 0.0d0 - !$OMP PARALLEL DO schedule(dynamic) PRIVATE(ni,nj,idx1,s12) + !$OMP PARALLEL DO schedule(dynamic) PRIVATE(ni,nj,idx1,s12,ktmp) do a = 1, nm1 ni = n1(a) do i = 1, ni @@ -1074,9 +1185,14 @@ subroutine fget_atomic_local_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nnei & t_width, d_width, cut_distance, order, & & pd, ang_norm2, distance_scale, angular_scale, alchemy) - kernels(:, idx1, b) = kernels(:, idx1, b) & - & + kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & - & kernel_idx, parameters) + ktmp = 0.0d0 + call kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & + & kernel_idx, parameters, ktmp) + kernels(:, idx1, b) = kernels(:, idx1, b) + ktmp + + !kernels(:, idx1, b) = kernels(:, idx1, b) & + ! & + kernel(self_scalar1(a, i), self_scalar2(b, j), s12, & + ! & kernel_idx, parameters) end do end do @@ -1085,6 +1201,7 @@ subroutine fget_atomic_local_kernels_fchl(x1, x2, verbose, n1, n2, nneigh1, nnei end do !$OMP END PARALLEL DO + deallocate (ktmp) deallocate (self_scalar1) deallocate (self_scalar2) deallocate (ksi1) diff --git a/src/qmllib/representations/slatm.py b/src/qmllib/representations/slatm.py index 7de31c32..98422660 100644 --- a/src/qmllib/representations/slatm.py +++ b/src/qmllib/representations/slatm.py @@ -169,7 +169,6 @@ def get_sbot( dgrid: float = 0.0262, pbc: str = "000", ) -> ndarray: - """ sigma -- standard deviation of gaussian distribution centered on a specific angle defaults to 0.05 (rad), approximately 3 degree diff --git a/src/qmllib/utils/__init__.py b/src/qmllib/utils/__init__.py index e69de29b..b5964d44 100644 --- a/src/qmllib/utils/__init__.py +++ b/src/qmllib/utils/__init__.py @@ -0,0 +1 @@ +from .fsettings import check_openmp, get_threads diff --git a/src/qmllib/utils/fib.f90 b/src/qmllib/utils/fib.f90 new file mode 100644 index 00000000..ced53fd4 --- /dev/null +++ b/src/qmllib/utils/fib.f90 @@ -0,0 +1,902 @@ +! fib.f90 +subroutine fib(a, n) + USE OMP_LIB + use iso_c_binding + integer(c_int), intent(in) :: n + integer(c_int), intent(out) :: a(n) + + !$OMP PARALLEL + PRINT *, "Hello from process: ", OMP_GET_THREAD_NUM() +!$OMP END PARALLEL + + do i = 1, n + if (i .eq. 1) then + a(i) = 0.0d0 + elseif (i .eq. 2) then + a(i) = 1.0d0 + else + a(i) = a(i - 1) + a(i - 2) + end if + end do +end + +function searchsorted(all_values, sorted) result(cdf_idx) + + implicit none + + double precision, dimension(:), intent(in) :: all_values + double precision, dimension(:), intent(in) :: sorted + + integer, allocatable, dimension(:) :: cdf_idx + + double precision :: val + + integer :: i, j, n, m + + n = size(all_values) - 1 + m = size(sorted) + + allocate (cdf_idx(n)) + + cdf_idx(:) = 0 + + do i = 1, n + + val = all_values(i) + + do j = 1, m + + !write (*,*) i, j, sorted(j), val + + ! if ((sorted(j) <= val) .and. (val < sorted(j+1))) then + if (sorted(j) > val) then + + cdf_idx(i) = j - 1 + !write(*,*) "found" + exit + + ! endif + else !iif (val > maxval(sorted)) then + cdf_idx(i) = m + end if + + end do + + end do + +end function searchsorted + +recursive subroutine quicksort(a, first, last) + implicit none + double precision :: a(*), x, t + integer first, last + integer i, j + + x = a((first + last)/2) + i = first + j = last + do + do while (a(i) < x) + i = i + 1 + end do + do while (x < a(j)) + j = j - 1 + end do + if (i >= j) exit + t = a(i); a(i) = a(j); a(j) = t + i = i + 1 + j = j - 1 + end do + if (first < i - 1) call quicksort(a, first, i - 1) + if (j + 1 < last) call quicksort(a, j + 1, last) +end subroutine quicksort + +subroutine fget_local_kernels_gaussian(q1, q2, n1, n2, sigmas, & + & nm1, nm2, nsigmas, kernels) + + implicit none + + double precision, dimension(:, :), intent(in) :: q1 + double precision, dimension(:, :), intent(in) :: q2 + + ! 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 + + ! Resulting alpha vector + double precision, dimension(nsigmas, nm1, nm2), intent(out) :: kernels + + ! Internal counters + integer :: a, b, i, j, k, ni, nj + + ! Temporary variables necessary for parallelization + double precision, allocatable, dimension(:, :) :: atomic_distance + + integer, allocatable, dimension(:) :: i_starts + integer, allocatable, dimension(:) :: j_starts + + allocate (i_starts(nm1)) + allocate (j_starts(nm2)) + + !$OMP PARALLEL DO + do i = 1, nm1 + i_starts(i) = sum(n1(:i)) - n1(i) + end do + !$OMP END PARALLEL DO + + !$OMP PARALLEL DO + do j = 1, nm2 + j_starts(j) = sum(n2(:j)) - n2(j) + end do + !$OMP END PARALLEL DO + + inv_sigma2(:) = -0.5d0/(sigmas(:))**2 + kernels(:, :, :) = 0.0d0 + + allocate (atomic_distance(maxval(n1), maxval(n2))) + atomic_distance(:, :) = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj) SCHEDULE(dynamic) COLLAPSE(2) + do a = 1, nm1 + do b = 1, nm2 + nj = n2(b) + ni = n1(a) + + atomic_distance(:, :) = 0.0d0 + do i = 1, ni + do j = 1, nj + + atomic_distance(i, j) = sum((q1(:, i + i_starts(a)) - q2(:, j + j_starts(b)))**2) + + end do + end do + + do k = 1, nsigmas + kernels(k, a, b) = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) + end do + + end do + end do + !$OMP END PARALLEL DO + + deallocate (atomic_distance) + deallocate (i_starts) + deallocate (j_starts) + +end subroutine fget_local_kernels_gaussian + +subroutine fget_local_kernels_laplacian(q1, q2, n1, n2, sigmas, & + & nm1, nm2, nsigmas, kernels) + + implicit none + + double precision, dimension(:, :), intent(in) :: q1 + double precision, dimension(:, :), intent(in) :: q2 + + ! 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 + + ! Resulting alpha vector + double precision, dimension(nsigmas, nm1, nm2), intent(out) :: kernels + + ! Internal counters + integer :: a, b, i, j, k, ni, nj + + ! Temporary variables necessary for parallelization + double precision, allocatable, dimension(:, :) :: atomic_distance + + integer, allocatable, dimension(:) :: i_starts + integer, allocatable, dimension(:) :: j_starts + + allocate (i_starts(nm1)) + allocate (j_starts(nm2)) + + !$OMP PARALLEL DO + do i = 1, nm1 + i_starts(i) = sum(n1(:i)) - n1(i) + end do + !$OMP END PARALLEL DO + + !$OMP PARALLEL DO + do j = 1, nm2 + j_starts(j) = sum(n2(:j)) - n2(j) + end do + !$OMP END PARALLEL DO + + inv_sigma2(:) = -1.0d0/sigmas(:) + kernels(:, :, :) = 0.0d0 + + allocate (atomic_distance(maxval(n1), maxval(n2))) + atomic_distance(:, :) = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj) SCHEDULE(dynamic) COLLAPSE(2) + do a = 1, nm1 + do b = 1, nm2 + nj = n2(b) + ni = n1(a) + + atomic_distance(:, :) = 0.0d0 + do i = 1, ni + do j = 1, nj + + atomic_distance(i, j) = sum(abs(q1(:, i + i_starts(a)) - q2(:, j + j_starts(b)))) + + end do + end do + + do k = 1, nsigmas + kernels(k, a, b) = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) + end do + + end do + end do + !$OMP END PARALLEL DO + + deallocate (atomic_distance) + deallocate (i_starts) + deallocate (j_starts) + +end subroutine fget_local_kernels_laplacian + +subroutine fget_vector_kernels_laplacian(q1, q2, n1, n2, sigmas, & + & nm1, nm2, nsigmas, kernels) + + implicit none + + ! Descriptors for the training set + double precision, dimension(:, :, :), intent(in) :: q1 + double precision, dimension(:, :, :), intent(in) :: q2 + + ! 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_sigma + + ! Resulting alpha vector + double precision, dimension(nsigmas, nm1, nm2), intent(out) :: kernels + + ! Internal counters + integer :: i, j, k, ni, nj, ia, ja + + ! Temporary variables necessary for parallelization + double precision, allocatable, dimension(:, :) :: atomic_distance + + inv_sigma(:) = -1.0d0/sigmas(:) + + kernels(:, :, :) = 0.0d0 + + allocate (atomic_distance(maxval(n1), maxval(n2))) + atomic_distance(:, :) = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj) SCHEDULE(dynamic) COLLAPSE(2) + do j = 1, nm2 + do i = 1, nm1 + ni = n1(i) + nj = n2(j) + + atomic_distance(:, :) = 0.0d0 + + do ja = 1, nj + do ia = 1, ni + + atomic_distance(ia, ja) = sum(abs(q1(:, ia, i) - q2(:, ja, j))) + + end do + end do + + do k = 1, nsigmas + kernels(k, i, j) = sum(exp(atomic_distance(:ni, :nj)*inv_sigma(k))) + end do + + end do + end do + !$OMP END PARALLEL DO + + deallocate (atomic_distance) + +end subroutine fget_vector_kernels_laplacian + +subroutine fget_vector_kernels_gaussian(q1, q2, n1, n2, sigmas, & + & nm1, nm2, nsigmas, kernels) + + implicit none + + ! Representations (n_samples, n_max_atoms, rep_size) + double precision, dimension(:, :, :), intent(in) :: q1 + double precision, dimension(:, :, :), intent(in) :: q2 + + ! 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 + + ! Resulting alpha vector + double precision, dimension(nsigmas, nm1, nm2), intent(out) :: kernels + + ! Internal counters + integer :: i, j, k, ni, nj, ia, ja + + ! Temporary variables necessary for parallelization + double precision, allocatable, dimension(:, :) :: atomic_distance + + inv_sigma2(:) = -0.5d0/(sigmas(:))**2 + + kernels(:, :, :) = 0.0d0 + + allocate (atomic_distance(maxval(n1), maxval(n2))) + atomic_distance(:, :) = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj,ja,ia) SCHEDULE(dynamic) COLLAPSE(2) + do j = 1, nm2 + do i = 1, nm1 + ni = n1(i) + nj = n2(j) + + atomic_distance(:, :) = 0.0d0 + + do ja = 1, nj + do ia = 1, ni + + atomic_distance(ia, ja) = sum((q1(:, ia, i) - q2(:, ja, j))**2) + + 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) + +end subroutine fget_vector_kernels_gaussian + +subroutine fget_vector_kernels_gaussian_symmetric(q, n, sigmas, & + & nm, nsigmas, kernels) + + implicit none + + ! Representations (rep_size, n_samples, n_max_atoms) + double precision, dimension(:, :, :), intent(in) :: q + + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: n + + ! Sigma in the Gaussian kernel + double precision, dimension(:), intent(in) :: sigmas + + ! Number of molecules + integer, intent(in) :: nm + + ! Number of sigmas + integer, intent(in) :: nsigmas + + ! Resulting kernels + double precision, dimension(nsigmas, nm, nm), intent(out) :: kernels + + ! Temporary variables necessary for parallelization + double precision, allocatable, dimension(:, :) :: atomic_distance + double precision, allocatable, dimension(:) :: inv_sigma2 + + ! Internal counters + integer :: i, j, k, ni, nj, ia, ja + double precision :: val + + allocate (inv_sigma2(nsigmas)) + + inv_sigma2 = -0.5d0/(sigmas)**2 + + kernels = 1.0d0 + + i = size(q, dim=3) + allocate (atomic_distance(i, i)) + atomic_distance(:, :) = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj,ja,ia,val) SCHEDULE(dynamic) COLLAPSE(2) + do j = 1, nm + do i = 1, nm + if (i .lt. j) cycle + ni = n(i) + nj = n(j) + + atomic_distance(:, :) = 0.0d0 + + do ja = 1, nj + do ia = 1, ni + + atomic_distance(ia, ja) = sum((q(:, ia, i) - q(:, ja, j))**2) + + end do + end do + + do k = 1, nsigmas + val = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) + kernels(k, i, j) = val + kernels(k, j, i) = val + end do + + end do + end do + !$OMP END PARALLEL DO + + deallocate (atomic_distance) + deallocate (inv_sigma2) + +end subroutine fget_vector_kernels_gaussian_symmetric + +subroutine fget_vector_kernels_laplacian_symmetric(q, n, sigmas, & + & nm, nsigmas, kernels) + + implicit none + + ! Representations (rep_size, n_samples, n_max_atoms) + double precision, dimension(:, :, :), intent(in) :: q + + ! List of numbers of atoms in each molecule + integer, dimension(:), intent(in) :: n + + ! Sigma in the Laplacian kernel + double precision, dimension(:), intent(in) :: sigmas + + ! Number of molecules + integer, intent(in) :: nm + + ! Number of sigmas + integer, intent(in) :: nsigmas + + ! Resulting kernels + double precision, dimension(nsigmas, nm, nm), intent(out) :: kernels + + ! Temporary variables necessary for parallelization + double precision, allocatable, dimension(:, :) :: atomic_distance + double precision, allocatable, dimension(:) :: inv_sigma2 + + ! Internal counters + integer :: i, j, k, ni, nj, ia, ja + double precision :: val + + allocate (inv_sigma2(nsigmas)) + + inv_sigma2 = -1.0d0/sigmas + + kernels = 1.0d0 + + i = size(q, dim=3) + allocate (atomic_distance(i, i)) + atomic_distance(:, :) = 0.0d0 + + !$OMP PARALLEL DO PRIVATE(atomic_distance,ni,nj,ja,ia,val) SCHEDULE(dynamic) COLLAPSE(2) + do j = 1, nm + do i = 1, nm + if (i .lt. j) cycle + ni = n(i) + nj = n(j) + + atomic_distance(:, :) = 0.0d0 + + do ja = 1, nj + do ia = 1, ni + + atomic_distance(ia, ja) = sum(abs(q(:, ia, i) - q(:, ja, j))) + + end do + end do + + do k = 1, nsigmas + val = sum(exp(atomic_distance(:ni, :nj)*inv_sigma2(k))) + kernels(k, i, j) = val + kernels(k, j, i) = val + end do + + end do + end do + !$OMP END PARALLEL DO + + deallocate (atomic_distance) + deallocate (inv_sigma2) + +end subroutine fget_vector_kernels_laplacian_symmetric + +subroutine fgaussian_kernel(a, na, b, nb, k, sigma) + + 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 + double precision, intent(in) :: sigma + + double precision, allocatable, dimension(:) :: temp + + double precision :: inv_sigma + integer :: i, j + + inv_sigma = -0.5d0/(sigma*sigma) + + allocate (temp(size(a, dim=1))) + + !$OMP PARALLEL DO PRIVATE(temp) COLLAPSE(2) + do i = 1, nb + do j = 1, na + temp(:) = a(:, j) - b(:, i) + k(j, i) = exp(inv_sigma*dot_product(temp, temp)) + end do + end do + !$OMP END PARALLEL DO + + deallocate (temp) + +end subroutine fgaussian_kernel + +subroutine fgaussian_kernel_symmetric(x, n, k, sigma) + + implicit none + + double precision, dimension(:, :), intent(in) :: x + + integer, intent(in) :: n + + double precision, dimension(:, :), intent(inout) :: k + double precision, intent(in) :: sigma + + double precision, allocatable, dimension(:) :: temp + double precision :: val + + double precision :: inv_sigma + integer :: i, j + + inv_sigma = -0.5d0/(sigma*sigma) + + k = 1.0d0 + + allocate (temp(size(x, dim=1))) + + !$OMP PARALLEL DO PRIVATE(temp, val) SCHEDULE(dynamic) + do i = 1, n + do j = i, n + temp = x(:, j) - x(:, i) + val = exp(inv_sigma*dot_product(temp, temp)) + k(j, i) = val + k(i, j) = val + end do + end do + !$OMP END PARALLEL DO + + deallocate (temp) + +end subroutine fgaussian_kernel_symmetric + +subroutine flaplacian_kernel(a, na, b, nb, k, sigma) + + 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 + double precision, intent(in) :: sigma + + double precision :: inv_sigma + + integer :: i, j + + inv_sigma = -1.0d0/sigma + + !$OMP PARALLEL DO COLLAPSE(2) + do i = 1, nb + do j = 1, na + k(j, i) = exp(inv_sigma*sum(abs(a(:, j) - b(:, i)))) + end do + end do + !$OMP END PARALLEL DO + +end subroutine flaplacian_kernel + +subroutine flaplacian_kernel_symmetric(x, n, k, sigma) + + implicit none + + double precision, dimension(:, :), intent(in) :: x + + integer, intent(in) :: n + + double precision, dimension(:, :), intent(inout) :: k + double precision, intent(in) :: sigma + + double precision :: val + + double precision :: inv_sigma + integer :: i, j + + inv_sigma = -1.0d0/sigma + + k = 1.0d0 + + !$OMP PARALLEL DO PRIVATE(val) SCHEDULE(dynamic) + do i = 1, n + do j = i, n + val = exp(inv_sigma*sum(abs(x(:, j) - x(:, i)))) + k(j, i) = val + k(i, j) = val + end do + end do + !$OMP END PARALLEL DO + +end subroutine flaplacian_kernel_symmetric + +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 COLLAPSE(2) + do i = 1, nb + do j = 1, na + k(j, i) = dot_product(a(:, j), b(:, i)) + end do + end do +!$OMP END PARALLEL DO + +end subroutine flinear_kernel + +subroutine fmatern_kernel_l2(a, na, b, nb, k, sigma, order) + + 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 + double precision, intent(in) :: sigma + integer, intent(in) :: order + + double precision, allocatable, dimension(:) :: temp + + double precision :: inv_sigma, inv_sigma2, d, d2 + integer :: i, j + + allocate (temp(size(a, dim=1))) + + if (order == 0) then + inv_sigma = -1.0d0/sigma + + !$OMP PARALLEL DO PRIVATE(temp) COLLAPSE(2) + do i = 1, nb + do j = 1, na + temp(:) = a(:, j) - b(:, i) + k(j, i) = exp(inv_sigma*sqrt(sum(temp*temp))) + end do + end do + !$OMP END PARALLEL DO + else if (order == 1) then + inv_sigma = -sqrt(3.0d0)/sigma + + !$OMP PARALLEL DO PRIVATE(temp, d) COLLAPSE(2) + do i = 1, nb + do j = 1, na + temp(:) = a(:, j) - b(:, i) + d = sqrt(sum(temp*temp)) + k(j, i) = exp(inv_sigma*d)*(1.0d0 - inv_sigma*d) + end do + end do + !$OMP END PARALLEL DO + else + inv_sigma = -sqrt(5.0d0)/sigma + inv_sigma2 = 5.0d0/(3.0d0*sigma*sigma) + + !$OMP PARALLEL DO PRIVATE(temp, d, d2) COLLAPSE(2) + do i = 1, nb + do j = 1, na + temp(:) = a(:, j) - b(:, i) + d2 = sum(temp*temp) + d = sqrt(d2) + k(j, i) = exp(inv_sigma*d)*(1.0d0 - inv_sigma*d + inv_sigma2*d2) + end do + end do + !$OMP END PARALLEL DO + end if + + deallocate (temp) + +end subroutine fmatern_kernel_l2 + +subroutine fsargan_kernel(a, na, b, nb, k, sigma, gammas, ng) + + implicit none + + double precision, dimension(:, :), intent(in) :: a + double precision, dimension(:, :), intent(in) :: b + double precision, dimension(:), intent(in) :: gammas + + integer, intent(in) :: na, nb, ng + + double precision, dimension(:, :), intent(inout) :: k + double precision, intent(in) :: sigma + + double precision, allocatable, dimension(:) :: prefactor + double precision :: inv_sigma + double precision :: d + + integer :: i, j, m + + inv_sigma = -1.0d0/sigma + + ! Allocate temporary + allocate (prefactor(ng)) + + !$OMP PARALLEL DO PRIVATE(d, prefactor) SCHEDULE(dynamic) COLLAPSE(2) + do i = 1, nb + do j = 1, na + d = sum(abs(a(:, j) - b(:, i))) + do m = 1, ng + prefactor(m) = gammas(m)*(-inv_sigma*d)**m + end do + k(j, i) = exp(inv_sigma*d)*(1 + sum(prefactor(:))) + end do + end do + !$OMP END PARALLEL DO + + ! Clean up + deallocate (prefactor) + +end subroutine fsargan_kernel + +subroutine fwasserstein_kernel(a, na, b, nb, k, sigma, p, q) + + implicit none + + double precision, dimension(:, :), intent(in) :: a + double precision, dimension(:, :), intent(in) :: b + + double precision, allocatable, dimension(:, :) :: asorted + double precision, allocatable, dimension(:, :) :: bsorted + + double precision, allocatable, dimension(:) :: rep + + integer, intent(in) :: na, nb + + double precision, dimension(:, :), intent(inout) :: k + double precision, intent(in) :: sigma + + integer, intent(in) :: p + integer, intent(in) :: q + + double precision :: inv_sigma + + integer :: i, j, l + integer :: rep_size + + double precision, allocatable, dimension(:) :: deltas + double precision, allocatable, dimension(:) :: all_values + + double precision, allocatable, dimension(:) :: a_cdf + double precision, allocatable, dimension(:) :: b_cdf + integer, allocatable, dimension(:) :: a_cdf_idx + integer, allocatable, dimension(:) :: b_cdf_idx + + rep_size = size(a, dim=1) + allocate (asorted(rep_size, na)) + allocate (bsorted(rep_size, nb)) + allocate (rep(rep_size)) + + allocate (all_values(rep_size*2)) + allocate (deltas(rep_size*2 - 1)) + + allocate (a_cdf(rep_size*2 - 1)) + allocate (b_cdf(rep_size*2 - 1)) + + allocate (a_cdf_idx(rep_size*2 - 1)) + allocate (b_cdf_idx(rep_size*2 - 1)) + + asorted(:, :) = a(:, :) + bsorted(:, :) = b(:, :) + + do i = 1, na + rep(:) = asorted(:, i) + call quicksort(rep, 1, rep_size) + asorted(:, i) = rep(:) + end do + + do i = 1, nb + rep(:) = bsorted(:, i) + call quicksort(rep, 1, rep_size) + bsorted(:, i) = rep(:) + end do + + !$OMP PARALLEL DO PRIVATE(all_values,a_cdf_idx,b_cdf_idx,a_cdf,b_cdf,deltas) + do j = 1, nb + do i = 1, na + + all_values(:rep_size) = asorted(:, i) + all_values(rep_size + 1:) = bsorted(:, j) + + call quicksort(all_values, 1, 2*rep_size) + + do l = 1, 2*rep_size - 1 + deltas(l) = all_values(l + 1) - all_values(l) + end do + + a_cdf_idx = searchsorted(all_values, asorted(:, i)) + b_cdf_idx = searchsorted(all_values, bsorted(:, j)) + + a_cdf(:) = a_cdf_idx(:) + b_cdf(:) = b_cdf_idx(:) + a_cdf(:) = a_cdf(:)/rep_size + b_cdf(:) = b_cdf(:)/rep_size + + ! k(i,j) = exp(-sum(abs(a_cdf-b_cdf)*deltas)/sigma) + k(i, j) = exp(-(sum((abs(a_cdf - b_cdf)**p)*deltas)**(1.0d0/p))**q/sigma) + + end do + end do + !$OMP END PARALLEL DO +end subroutine fwasserstein_kernel diff --git a/src/qmllib/utils/fsettings.f90 b/src/qmllib/utils/fsettings.f90 new file mode 100644 index 00000000..2a43a145 --- /dev/null +++ b/src/qmllib/utils/fsettings.f90 @@ -0,0 +1,21 @@ + + subroutine check_openmp(compiled_with_openmp) + + implicit none + logical, intent(out):: compiled_with_openmp + + compiled_with_openmp = .false. + +!$ compiled_with_openmp = .true. + + end subroutine check_openmp + + function get_threads() result(nt) +!$ use omp_lib + implicit none + integer :: nt + + nt = 0 +!$ nt = omp_get_max_threads() + + end function get_threads diff --git a/tests/test_fchl_scalar.py b/tests/test_fchl_scalar.py index 9685c2e1..04b36d0c 100644 --- a/tests/test_fchl_scalar.py +++ b/tests/test_fchl_scalar.py @@ -2,9 +2,8 @@ from conftest import ASSETS, get_energies, shuffle_arrays from scipy.special import binom, factorial, jn -from qmllib.representations.fchl import generate_representation -from qmllib.representations.fchl import generate_representation as generate_fchl_representation from qmllib.representations.fchl import ( + generate_representation, get_atomic_kernels, get_atomic_symmetric_kernels, get_global_kernels, @@ -41,7 +40,7 @@ def _get_training_data(n_points, representation_options={}): # Associate a property (heat of formation) with the object all_properties.append(data[xyz_file]) - representation = generate_fchl_representation(coord, atoms, **_representation_options) + representation = generate_representation(coord, atoms, **_representation_options) assert ( representation.shape[0] == _representation_options["max_size"] @@ -1252,6 +1251,8 @@ def test_fchl_l2(): K_test = np.zeros((n_points, n_points)) + print(K) + # UNUSED sigma = 2.0 # UNUSED v = 3 # UNUSED n = 2 @@ -1268,6 +1269,9 @@ def test_fchl_l2(): K_test[i, j] += np.exp(Sij[ii, jj] * inv_sigma) + print(K_test) + print(np.max(K - K_test)) + assert np.allclose(K, K_test), "Error in FCHL l2 kernels"