Commit a81e35a0 authored by Artem Pulkin's avatar Artem Pulkin

potentials: add behler4 and fix issues in potential tests

parent 1315f616
Pipeline #43465 passed with stage
in 2 minutes and 30 seconds
This diff is collapsed.
......@@ -8,7 +8,8 @@ from scipy.sparse import csr_matrix
from ._potentials import kernel_general_2, kernel_g_general_2, kernel_general_3, kernel_g_general_3, \
kernel_lj, kernel_g_lj, kernel_sw_phi2, kernel_g_sw_phi2, \
kernel_sw_phi3, kernel_g_sw_phi3, kernel_mlsf_g2, kernel_g_mlsf_g2, kernel_harmonic_repulsion, \
kernel_g_harmonic_repulsion, kernel_sigmoid, kernel_g_sigmoid, kernel_mlsf_g5, kernel_g_mlsf_g5
kernel_g_harmonic_repulsion, kernel_sigmoid, kernel_g_sigmoid, kernel_mlsf_g5, kernel_g_mlsf_g5, \
kernel_mlsf_g4, kernel_g_mlsf_g4
from ._util import calc_sparse_distances
from .util import num_grad, ParameterUnits
......@@ -795,6 +796,16 @@ behler5_descriptor_family = LocalPotentialFamily(
tag="Behler type 5",
pre_compute_r=[sine_cutoff_fn, sine_cutoff_fp, exp_fn],
)
behler4_descriptor_family = LocalPotentialFamily(
coordination_number=3,
parameters=dict(a=(0, np.inf), eta=(0, np.inf), l=(-1, 1), zeta=(0, np.inf)),
parameter_units=ParameterUnits(dict(a="angstrom", eta="1/angstrom**2")),
cutoff="a",
kernel=kernel_mlsf_g4,
kernel_gradient=kernel_g_mlsf_g4,
tag="Behler type 4",
pre_compute_r=[sine_cutoff_fn, sine_cutoff_fp, exp_fn],
)
def behler2_p2(r, a, eta, r_sphere):
......
......@@ -170,17 +170,33 @@ class PotentialTestMixin:
self.potential("kernel_gradient", self.r, self.cartesian_row, self.cartesian_col, out=out)
testing.assert_allclose(out, val_ref)
def test_numgrad(self, atol=1e-10, rtol=1e-3):
testing.assert_allclose(
self.potential("kernel_gradient", self.r, self.cartesian_row, self.cartesian_col),
self.potential("kernel_gradient_d", self.r, self.cartesian_row, self.cartesian_col, out=np.zeros(
potentials.shape_for_kernel("kernel_gradient", len(self.cartesian_row)), dtype=float,
)),
atol=atol, rtol=rtol,
)
def test_numgrad(self, atol=1e-10, rtol=1e-3, debug=False):
test = self.potential("kernel_gradient", self.r, self.cartesian_row, self.cartesian_col)
ref = self.potential("kernel_gradient_d", self.r, self.cartesian_row, self.cartesian_col, out=np.zeros(
potentials.shape_for_kernel("kernel_gradient", len(self.cartesian_row)), dtype=float,
))
if debug:
test_nz = np.nonzero(test)
ref_nz = np.nonzero(ref)
print(test_nz, ref_nz)
print(test[test_nz], ref[ref_nz])
testing.assert_allclose(test, ref, atol=atol, rtol=rtol)
class RefPotentialMixin:
def _test_ref_numgrad(self, atol=1e-10, rtol=1e-3, debug=False): # This test is not stable because of units
self.potential_ref.kernels["kernel_gradient_d"] = potentials.num_grad_potential(self.potential_ref.kernels["kernel"])
test = self.potential_ref("kernel_gradient", self.r, self.cartesian_row, self.cartesian_col)
ref = self.potential_ref("kernel_gradient_d", self.r, self.cartesian_row, self.cartesian_col, out=np.zeros(
potentials.shape_for_kernel("kernel_gradient", len(self.cartesian_row)), dtype=float,
))
if debug:
test_nz = np.nonzero(test)
ref_nz = np.nonzero(ref)
print(test_nz, ref_nz)
print(test[test_nz], ref[ref_nz])
testing.assert_allclose(test, ref, atol=atol, rtol=rtol)
def test_val_ref(self):
testing.assert_allclose(
self.potential("kernel", self.r, self.cartesian_row, self.cartesian_col),
......@@ -324,6 +340,7 @@ class AnglePotentialTestMixin(PotentialTestMixin, RefPotentialMixin):
_i = np.arange(r1.shape[1])
_r1[:, _i, _i] = 0
_r2[:, _i, _i] = 0
_cos[:, _i, _i] = 0
result_1 = _r1[..., np.newaxis] * uvecs[:, :, np.newaxis] + _cos[..., np.newaxis] * dc_dr1 # [n_rows, n_columns, n_columns, 3]
result_2 = _r2[..., np.newaxis] * uvecs[:, np.newaxis, :] + _cos[..., np.newaxis] * dc_dr2 # [n_rows, n_columns, n_columns, 3]
......@@ -546,6 +563,56 @@ class TestBehlerSF5(AnglePotentialTestMixin, SerializationMixin, TestCase):
return a
class TestBehlerSF4(AnglePotentialTestMixin, SerializationMixin, TestCase):
@classmethod
def setUpClassPotential(cls):
cls.potential = potentials.behler4_descriptor_family(a=.19, eta=5, l=0.1, zeta=4)
cls.degenerate = False
@classmethod
def __r3__(cls, r1, r2, r12_cos):
return (r1 ** 2 + r2 ** 2 - 2 * r1 * r2 * r12_cos) ** .5
@classmethod
def __r3_factor__(cls, r1, r2, r12_cos, a, eta):
r3 = cls.__r3__(r1, r2, r12_cos)
return np.exp(- eta * r3 ** 2) * (.5 + np.cos(np.pi * r3 / a) / 2) * (r3 < a)
@classmethod
def __r3_factor_prime__(cls, r1, r2, r12_cos, a, eta):
r3 = cls.__r3__(r1, r2, r12_cos)
return - np.exp(- eta * r3 ** 2) * (2 * eta * r3 * (.5 + np.cos(np.pi * r3 / a) / 2) + np.pi / a * np.sin(np.pi * r3 / a) / 2) * (r3 < a)
@classmethod
def __raw_fun__(cls, r1, r2, r12_cos, a, eta, l, zeta):
return TestBehlerSF5.__raw_fun__(r1, r2, r12_cos, a, eta, l, zeta) * cls.__r3_factor__(r1, r2, r12_cos, a, eta)
@classmethod
def __raw_fun_prime_r1__(cls, r1, r2, r12_cos, a, eta, l, zeta):
r3 = cls.__r3__(r1, r2, r12_cos)
dr3_dr1 = (r1 - r2 * r12_cos) / r3
return TestBehlerSF5.__raw_fun_prime_r1__(r1, r2, r12_cos, a, eta, l, zeta) * cls.__r3_factor__(r1, r2, r12_cos, a, eta) + \
TestBehlerSF5.__raw_fun__(r1, r2, r12_cos, a, eta, l, zeta) * cls.__r3_factor_prime__(r1, r2, r12_cos, a, eta) * dr3_dr1
@classmethod
def __raw_fun_prime_r2__(cls, r1, r2, r12_cos, a, eta, l, zeta):
r3 = cls.__r3__(r1, r2, r12_cos)
dr3_dr2 = (r2 - r1 * r12_cos) / r3
return TestBehlerSF5.__raw_fun_prime_r2__(r1, r2, r12_cos, a, eta, l, zeta) * cls.__r3_factor__(r1, r2, r12_cos, a, eta) + \
TestBehlerSF5.__raw_fun__(r1, r2, r12_cos, a, eta, l, zeta) * cls.__r3_factor_prime__(r1, r2, r12_cos, a, eta) * dr3_dr2
@classmethod
def __raw_fun_prime_cos__(cls, r1, r2, r12_cos, a, eta, l, zeta):
r3 = cls.__r3__(r1, r2, r12_cos)
dr3_dcos = -r1 * r2 / r3
return TestBehlerSF5.__raw_fun_prime_cos__(r1, r2, r12_cos, a, eta, l, zeta) * cls.__r3_factor__(r1, r2, r12_cos, a, eta) + \
TestBehlerSF5.__raw_fun__(r1, r2, r12_cos, a, eta, l, zeta) * cls.__r3_factor_prime__(r1, r2, r12_cos, a, eta) * dr3_dcos
@classmethod
def __cutoff__(cls, a, eta, l, zeta):
return a
# TODO: add gradients here
class NWPotentialTestMixin:
......
......@@ -80,6 +80,50 @@
"grad_cosine": "zeta * l * pow(1 + l * r12_cos, zeta - 1) * _fn_exponent * _fn_cutoff1 * _fn_cutoff2",
"r12_symmetry_allowed": "1"
},
{
"name": "mlsf_g4",
"type": "potential-3",
"parameters": "double eta, double l, double zeta, double[:, ::1] pre_compute_r, int[::1] pre_compute_r_handles",
"preamble": [
"cdef double _fn_cutoff1, _fn_cutoff2, _fn_cutoff3, _fn_exponent, _fn_pw",
"cdef double _prefactor = pow(2, 1 - zeta)",
"cdef int pre_compute_r_cutoff_fn_handle = pre_compute_r_handles[0], pre_compute_r_exp_fn_handle = pre_compute_r_handles[2]",
"cdef double _r3, _r3_factor, g5_fun"
],
"preamble_grad": [
"cdef double _fp_cutoff1, _fp_cutoff2",
"cdef int pre_compute_r_cutoff_fp_handle = pre_compute_r_handles[1]",
"cdef double _r3_factor_p, _r3_grad_r1, _r3_grad_r2, _r3_grad_cosine, g5_grad_r1, g5_grad_r2, g5_grad_cosine"
],
"before1": "_fn_cutoff1 = pre_compute_r[ptr1, pre_compute_r_cutoff_fn_handle]",
"before1_grad": "_fp_cutoff1 = pre_compute_r[ptr1, pre_compute_r_cutoff_fp_handle]",
"before": [
"_r3 = sqrt(r1 * r1 + r2 * r2 - 2 * r1 * r2 * r12_cos)"
],
"final_filter": "_r3 < a",
"before_inner": [
"_fn_cutoff2 = pre_compute_r[ptr2, pre_compute_r_cutoff_fn_handle]",
"_fn_exponent = _prefactor * pre_compute_r[ptr1, pre_compute_r_exp_fn_handle] * pre_compute_r[ptr2, pre_compute_r_exp_fn_handle]",
"_fn_pw = pow(1 + l * r12_cos, zeta)",
"g5_fun = _fn_pw * _fn_exponent * _fn_cutoff1 * _fn_cutoff2",
"_r3_factor = exp(-eta * _r3 * _r3) * (.5 + cos(pi * _r3 / a) / 2)"
],
"before_inner_grad": [
"_fp_cutoff2 = pre_compute_r[ptr2, pre_compute_r_cutoff_fp_handle]",
"g5_grad_r1 = - 2 * eta * r1 * g5_fun + _fn_pw * _fn_exponent * _fn_cutoff2 * _fp_cutoff1",
"g5_grad_r2 = - 2 * eta * r2 * g5_fun + _fn_pw * _fn_exponent * _fn_cutoff1 * _fp_cutoff2",
"g5_grad_cosine = zeta * l * pow(1 + l * r12_cos, zeta - 1) * _fn_exponent * _fn_cutoff1 * _fn_cutoff2",
"_r3_factor_p = - 2 * eta * _r3 * _r3_factor - exp(- eta * _r3 * _r3) * sin(pi * _r3 / a) * pi / a / 2",
"_r3_grad_r1 = (r1 - r2 * r12_cos) / _r3",
"_r3_grad_r2 = (r2 - r1 * r12_cos) / _r3",
"_r3_grad_cosine = - r1 * r2 / _r3"
],
"kernel": "g5_fun * _r3_factor",
"grad_r1": "g5_grad_r1 * _r3_factor + g5_fun * _r3_factor_p * _r3_grad_r1",
"grad_r2": "g5_grad_r2 * _r3_factor + g5_fun * _r3_factor_p * _r3_grad_r2",
"grad_cosine": "g5_grad_cosine * _r3_factor + g5_fun * _r3_factor_p * _r3_grad_cosine",
"r12_symmetry_allowed": "1"
},
{
"name": "sigmoid",
"type": "potential-2",
......
......@@ -73,8 +73,9 @@ class BuildPotentials(Command):
r12_symmetry_allowed="0",
parameters='',
degenerate=False,
final_filter=True,
)
for k in "preamble", "preamble_grad", "before", "before_grad", "before1", "before1_grad", "after", "after_grad":
for k in "preamble", "preamble_grad", "before", "before_grad", "before1", "before1_grad", "before_inner", "before_inner_grad":
defaults[k] = f"# (no '{k}' statements)"
for p in potentials:
_p = defaults.copy()
......
......@@ -43,8 +43,6 @@ def kernel_g_$name(
function_value = $kernel
# --- grad ---
g = $grad
# --- after ---
$after
# --------------
for dim in range(3):
x = (cartesian_row[row, dim] - cartesian_col[col, dim]) / r * g
......
......@@ -39,6 +39,3 @@ def kernel_$name(
$before
# --- kernel ---
out[row] += $kernel
# --- after ---
$after
# --------------
......@@ -62,30 +62,30 @@ def kernel_g_$name(
# --- before ---
$before
$before_grad
# --- kernel ---
function_value = $kernel
# --- grad ---
dfunc_dr1 = $grad_r1
dfunc_dr2 = $grad_r2
dfunc_dct = $grad_cosine
# --- after ---
$after
$after_grad
# --------------
if r12_symmetry_allowed and ptr1 != ptr2:
dfunc_dr1 *= 2
dfunc_dr2 *= 2
dfunc_dct *= 2
if $final_filter:
# --- before ---
$before_inner
$before_inner_grad
# --- kernel ---
function_value = $kernel
# --- grad ---
dfunc_dr1 = $grad_r1
dfunc_dr2 = $grad_r2
dfunc_dct = $grad_cosine
if r12_symmetry_allowed and ptr1 != ptr2:
dfunc_dr1 *= 2
dfunc_dr2 *= 2
dfunc_dct *= 2
# Derivatives
# Derivatives
for dim in range(3):
nx1 = (cartesian_row[row, dim] - cartesian_col[col1, dim]) / r1
nx2 = (cartesian_row[row, dim] - cartesian_col[col2, dim]) / r2
# d(cos θ) / dr = 1/r (n_s - n_r cos θ)
cx1 = (nx2 - nx1 * r12_cos) / r1
cx2 = (nx1 - nx2 * r12_cos) / r2
for dim in range(3):
nx1 = (cartesian_row[row, dim] - cartesian_col[col1, dim]) / r1
nx2 = (cartesian_row[row, dim] - cartesian_col[col2, dim]) / r2
# d(cos θ) / dr = 1/r (n_s - n_r cos θ)
cx1 = (nx2 - nx1 * r12_cos) / r1
cx2 = (nx1 - nx2 * r12_cos) / r2
out[row, row, dim] += nx1 * dfunc_dr1 + cx1 * dfunc_dct + nx2 * dfunc_dr2 + cx2 * dfunc_dct # df_self / dr_self
out[row, col1_, dim] -= nx1 * dfunc_dr1 + cx1 * dfunc_dct # df_self / dr_n1
out[row, col2_, dim] -= nx2 * dfunc_dr2 + cx2 * dfunc_dct # df_self / dr_n2
out[row, row, dim] += nx1 * dfunc_dr1 + cx1 * dfunc_dct + nx2 * dfunc_dr2 + cx2 * dfunc_dct # df_self / dr_self
out[row, col1_, dim] -= nx1 * dfunc_dr1 + cx1 * dfunc_dct # df_self / dr_n1
out[row, col2_, dim] -= nx2 * dfunc_dr2 + cx2 * dfunc_dct # df_self / dr_n2
......@@ -55,8 +55,8 @@ def kernel_$name(
r12_cos /= r1 * r2
# --- before ---
$before
# --- kernel ---
out[row] += (1 + r12_symmetry_allowed * (ptr1 != ptr2)) * ($kernel)
# --- after ---
$after
# --------------
if $final_filter:
# --- before ---
$before_inner
# --- kernel ---
out[row] += (1 + r12_symmetry_allowed * (ptr1 != ptr2)) * ($kernel)
......@@ -2,6 +2,6 @@
import numpy as np
cimport numpy as np
cimport cython
from libc.math cimport exp, cos, sin, pi, pow
from libc.math cimport exp, cos, sin, pi, pow, sqrt
$kernels
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment