From 5592b4b7b49efe60a1ba86392630233712777736 Mon Sep 17 00:00:00 2001 From: nw13slx Date: Wed, 22 Apr 2020 22:43:28 -0400 Subject: [PATCH 1/6] fix etype issue --- flare/env.py | 7 ++++--- flare/kernels/mc_sephyps.py | 8 ++++---- flare/kernels/mc_simple.py | 23 ++++++++++++++--------- 3 files changed, 22 insertions(+), 16 deletions(-) diff --git a/flare/env.py b/flare/env.py index f2b1632bb..b72561ac0 100644 --- a/flare/env.py +++ b/flare/env.py @@ -55,7 +55,8 @@ def compute_env(self): # if 3 cutoffs are given, create many-body arrays if len(self.cutoffs) > 2: - self.bond_array_mb, self.neigh_dists_mb, self.num_neighs_mb, self.etype_mb = get_m_body_arrays( + self.bond_array_mb, self.neigh_dists_mb, self.num_neighs_mb, self.etype_mb, \ + self.bond_array_mb_etypes = get_m_body_arrays( self.positions, self.atom, self.cell, self.cutoffs[2], self.species) else: self.bond_array_mb = None @@ -406,7 +407,7 @@ class to allow for njit acceleration with Numba. """ # TODO: this can be probably improved using stored arrays, redundant calls to get_2_body_arrays # Get distances, positions, species and indexes of neighbouring atoms - bond_array_mb, __, _, bond_inds = get_2_body_arrays_ind( + bond_array_mb, __, etypes, bond_inds = get_2_body_arrays_ind( positions, atom, cell, cutoff_mb, species) # For each neighbouring atom, get distances in its neighbourhood @@ -431,7 +432,7 @@ class to allow for njit acceleration with Numba. etypes_mb_array[i, :num_neighs_mb[i]] = neighbouring_etypes[i] - return bond_array_mb, neigh_dists_mb, num_neighs_mb, etypes_mb_array + return bond_array_mb, neigh_dists_mb, num_neighs_mb, etypes_mb_array, etypes if __name__ == '__main__': diff --git a/flare/kernels/mc_sephyps.py b/flare/kernels/mc_sephyps.py index 182f0da0c..ac7c86e46 100644 --- a/flare/kernels/mc_sephyps.py +++ b/flare/kernels/mc_sephyps.py @@ -96,7 +96,7 @@ def two_three_many_body_mc(env1, env2, d1, d2, cutoffs, many_term = many_body_mc_jit(env1.bond_array_mb, env2.bond_array_mb, env1.neigh_dists_mb, env2.neigh_dists_mb, env1.num_neighs_mb, env2.num_neighs_mb, - env1.ctype, env2.ctype, env1.etypes, env2.etypes, + env1.ctype, env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, env1.etype_mb, env2.etype_mb, env1.species, env2.species, d1, d2, sigm, lsm, r_cut_m, cutoff_func) @@ -134,7 +134,7 @@ def two_three_many_body_mc_grad(env1, env2, d1, d2, cutoffs, kern_many, gradm = many_body_mc_grad_jit(env1.bond_array_mb, env2.bond_array_mb, env1.neigh_dists_mb, env2.neigh_dists_mb, env1.num_neighs_mb, env2.num_neighs_mb, env1.ctype, - env2.ctype, env1.etypes, env2.etypes, + env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, env1.etype_mb, env2.etype_mb, env1.species, env2.species, d1, d2, sigm, lsm, r_cut_m, cutoff_func) @@ -174,7 +174,7 @@ def two_three_many_mc_force_en(env1, env2, d1, cutoffs, many_term = many_body_mc_force_en_jit(env1.bond_array_mb, env2.bond_array_mb, env1.neigh_dists_mb, env1.num_neighs_mb, - env1.ctype, env2.ctype, env1.etypes, env2.etypes, + env1.ctype, env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, env1.etype_mb, env1.species, env2.species, d1, sigm, lsm, r_cut_m, cutoff_func) @@ -209,7 +209,7 @@ def two_three_many_mc_en(env1, env2, cutoffs, triplet_mask) many_term = many_body_mc_en_jit(env1.bond_array_2, env2.bond_array_2, env1.ctype, - env2.ctype, env1.etypes, env2.etypes, env1.species, + env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, env1.species, env2.species, sigm, lsm, r_cut_m, cutoff_func) diff --git a/flare/kernels/mc_simple.py b/flare/kernels/mc_simple.py index 47252e8de..e7e3391ec 100644 --- a/flare/kernels/mc_simple.py +++ b/flare/kernels/mc_simple.py @@ -252,7 +252,8 @@ def two_plus_three_plus_many_body_mc(env1: AtomicEnvironment, env2: AtomicEnviro many_term = many_body_mc_jit(env1.bond_array_mb, env2.bond_array_mb, env1.neigh_dists_mb, env2.neigh_dists_mb, env1.num_neighs_mb, env2.num_neighs_mb, - env1.ctype, env2.ctype, env1.etypes, env2.etypes, + env1.ctype, env2.ctype, + env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, env1.etype_mb, env2.etype_mb, env1.species, env2.species, d1, d2, sigm, lsm, r_cut_m, cutoff_func) @@ -305,7 +306,8 @@ def two_plus_three_plus_many_body_mc_grad(env1: AtomicEnvironment, env2: AtomicE kern_many, gradm = many_body_mc_grad_jit(env1.bond_array_mb, env2.bond_array_mb, env1.neigh_dists_mb, env2.neigh_dists_mb, env1.num_neighs_mb, env2.num_neighs_mb, env1.ctype, - env2.ctype, env1.etypes, env2.etypes, + env2.ctype, env1.bond_array_mb_etypes, + env2.bond_array_mb_etypes, env1.etype_mb, env2.etype_mb, env1.species, env2.species, d1, d2, sigm, lsm, r_cut_m, cutoff_func) @@ -359,7 +361,8 @@ def two_plus_three_plus_many_body_mc_force_en(env1: AtomicEnvironment, env2: Ato many_term = many_body_mc_force_en_jit(env1.bond_array_mb, env2.bond_array_mb, env1.neigh_dists_mb, env1.num_neighs_mb, - env1.ctype, env2.ctype, env1.etypes, env2.etypes, + env1.ctype, env2.ctype, env1.bond_array_mb_etypes, + env2.bond_array_mb_etypes, env1.etype_mb, env1.species, env2.species, d1, sigm, lsm, r_cut_m, cutoff_func) @@ -408,7 +411,8 @@ def two_plus_three_plus_many_body_mc_en(env1: AtomicEnvironment, env2: AtomicEnv sig3, ls3, r_cut_3, cutoff_func) many_term = many_body_mc_en_jit(env1.bond_array_2, env2.bond_array_2, env1.ctype, - env2.ctype, env1.etypes, env2.etypes, env1.species, + env2.ctype, env1.bond_array_mb_etypes, + env2.bond_array_mb_etypes, env1.species, env2.species, sigm, lsm, r_cut_m, cutoff_func) @@ -707,7 +711,7 @@ def many_body_mc(env1: AtomicEnvironment, env2: AtomicEnvironment, # Get atomic species of central atom, their neighbours, and their neighbours' neighbours c1, c2 = env1.ctype, env2.ctype - etypes1, etypes2 = env1.etypes, env2.etypes + etypes1, etypes2 = env1.bond_array_mb_etypes, env2.bond_array_mb_etypes etypes_neigh_1, etypes_neigh_2 = env1.etype_mb, env2.etype_mb return many_body_mc_jit(bond_array_1, bond_array_2, neigh_dists_1, neigh_dists_2, num_neigh_1, @@ -735,7 +739,7 @@ def many_body_mc_grad(env1: AtomicEnvironment, env2: AtomicEnvironment, num_neigh_2 = env2.num_neighs_mb c1, c2 = env1.ctype, env2.ctype - etypes1, etypes2 = env1.etypes, env2.etypes + etypes1, etypes2 = env1.bond_array_mb_etypes, env2.bond_array_mb_etypes etypes_neigh_1, etypes_neigh_2 = env1.etype_mb, env2.etype_mb kernel, kernel_grad = many_body_mc_grad_jit(bond_array_1, bond_array_2, neigh_dists_1, @@ -773,7 +777,7 @@ def many_body_mc_force_en(env1, env2, d1, hyps, cutoffs, num_neigh_1 = env1.num_neighs_mb c1, c2 = env1.ctype, env2.ctype - etypes1, etypes2 = env1.etypes, env2.etypes + etypes1, etypes2 = env1.bond_array_mb_etypes, env2.bond_array_mb_etypes etypes_neigh_1 = env1.etype_mb # divide by three to account for triple counting @@ -802,8 +806,9 @@ def many_body_mc_en(env1: AtomicEnvironment, env2: AtomicEnvironment, ls = hyps[1] r_cut = cutoffs[2] - return many_body_mc_en_jit(env1.bond_array_2, env2.bond_array_2, env1.ctype, - env2.ctype, env1.etypes, env2.etypes, env1.species, env2.species, + return many_body_mc_en_jit(env1.bond_array_mb, env2.bond_array_mb, env1.ctype, + env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, + env1.species, env2.species, sig, ls, r_cut, cutoff_func) From 85365a79245b13e95ea57fe557ff08c7497b2c52 Mon Sep 17 00:00:00 2001 From: nw13slx Date: Wed, 22 Apr 2020 23:02:28 -0400 Subject: [PATCH 2/6] update unit test --- tests/fake_gp.py | 84 ++++++- tests/test_mc_mb_kernel.py | 438 ++++++------------------------------- 2 files changed, 144 insertions(+), 378 deletions(-) diff --git a/tests/fake_gp.py b/tests/fake_gp.py index daba2bd85..2816b296a 100644 --- a/tests/fake_gp.py +++ b/tests/fake_gp.py @@ -9,6 +9,7 @@ from flare.struc import Structure from flare.otf_parser import OtfAnalysis + def get_random_structure(cell, unique_species, noa): """Create a random test structure """ np.random.seed(0) @@ -35,11 +36,11 @@ def generate_hm(nbond, ntriplet, cutoffs=[1, 1], constraint=False, multihyps=Tru if (nbond > 0): nbond = 1 hyps_label += ['Length', 'Signal Var.'] - if (ntriplet >0): + if (ntriplet > 0): ntriplet = 1 hyps_label += ['Length', 'Signal Var.'] hyps_label += ['Noise Var.'] - return random((nbond+ntriplet)*2+1), {'hyps_label':hyps_label}, deepcopy(cutoffs) + return random((nbond+ntriplet)*2+1), {'hyps_label': hyps_label}, deepcopy(cutoffs) specs_mask = np.zeros(118, dtype=int) specs_mask[1] = 0 @@ -55,7 +56,7 @@ def generate_hm(nbond, ntriplet, cutoffs=[1, 1], constraint=False, multihyps=Tru cut += [cutoffs[0]] cut += [cutoffs[1]] - if (nbond>1): + if (nbond > 1): sig1 = random(nbond) ls1 = random(nbond) bond_mask = np.ones(nspecs**2, dtype=int) @@ -67,7 +68,7 @@ def generate_hm(nbond, ntriplet, cutoffs=[1, 1], constraint=False, multihyps=Tru bond_mask = np.zeros(nspecs**2, dtype=int) bond_name = ["sig2"]+["ls2"] - if (ntriplet>1): + if (ntriplet > 1): sig2 = random(ntriplet) ls2 = random(ntriplet) triplet_mask = np.ones(nspecs**3, dtype=int) @@ -81,10 +82,10 @@ def generate_hm(nbond, ntriplet, cutoffs=[1, 1], constraint=False, multihyps=Tru sigman = [0.05] - if (nbond>0 and ntriplet>0): + if (nbond > 0 and ntriplet > 0): hyps = np.hstack([sig1, ls1, sig2, ls2, sigman]) hyps_label = np.hstack([bond_name, triplet_name, ['noise']]) - elif (nbond>0): + elif (nbond > 0): hyps = np.hstack([sig1, ls1, sigman]) hyps_label = np.hstack([bond_name, ['noise']]) else: @@ -108,7 +109,7 @@ def generate_hm(nbond, ntriplet, cutoffs=[1, 1], constraint=False, multihyps=Tru count = 0 newhyps = [] newlabel = [] - if (nbond>1): + if (nbond > 1): # fix type 0, and only compute type 1 of bonds hm += [1] newhyps += [hyps[1]] @@ -117,13 +118,13 @@ def generate_hm(nbond, ntriplet, cutoffs=[1, 1], constraint=False, multihyps=Tru newhyps += [hyps[3]] newlabel += [hyps_label[3]] count += 4 - elif (nbond>0): + elif (nbond > 0): # fix sigma, and only vary ls hm += [1] newhyps += [hyps[1]] newlabel += [hyps_label[1]] count += 2 - if (ntriplet>1): + if (ntriplet > 1): # fix type 0, and only compute type 1 of triplets hm += [1+count] newhyps += [hyps[1+count]] @@ -131,7 +132,7 @@ def generate_hm(nbond, ntriplet, cutoffs=[1, 1], constraint=False, multihyps=Tru hm += [3+count] newhyps += [hyps[3+count]] newlabel += [hyps_label[3+count]] - elif (ntriplet>0): + elif (ntriplet > 0): # fix sigma, and only vary ls hm += [1+count] newhyps += [hyps[1+count]] @@ -185,6 +186,7 @@ def get_gp(bodies, kernel_type='mc', multihyps=True) -> GaussianProcess: return gaussian + def get_params(): parameters = {'unique_species': [2, 1], 'cutoff': 0.8, @@ -193,6 +195,7 @@ def get_params(): 'db_pts': 30} return parameters + def get_tstp() -> AtomicEnvironment: """Create test point for kernel to compare against""" # params @@ -209,6 +212,67 @@ def get_tstp() -> AtomicEnvironment: return test_pt +def generate_mb_envs(cutoffs, cell, delt, d1, mask=None): + positions0 = np.array([[0., 0., 0.], + [0., 0.3, 0.], + [0.3, 0., 0.], + [0.0, 0., 0.3], + [1., 1., 0.]]) + positions0[1:] += 0.1*np.random.random([4, 3]) + triplet = [1, 1, 2, 1] + np.random.shuffle(triplet) + species_1 = np.hstack([triplet, randint(1, 2)]) + return generate_mb_envs_pos(positions0, species_1, cutoffs, cell, delt, d1, mask) + + +def generate_mb_twin_envs(cutoffs, cell, delt, d1, mask=None): + + positions0 = np.array([[0., 0., 0.], + [0., 0.3, 0.], + [0.3, 0., 0.], + [0.0, 0., 0.3], + [1., 1., 0.]]) + positions0[1:] += 0.1*np.random.random([4, 3]) + + species_1 = np.array([1, 2, 1, 1, 1], dtype=np.int) + env1 = generate_mb_envs_pos( + positions0, species_1, cutoffs, cell, delt, d1, mask) + + positions1 = positions0[[0, 2, 3, 4]] + species_2 = species_1[[0, 2, 3, 4]] + env2 = generate_mb_envs_pos( + positions1, species_2, cutoffs, cell, delt, d1, mask) + + return env1, env2 + + +def generate_mb_envs_pos(positions0, species_1, cutoffs, cell, delt, d1, mask=None): + + positions = [positions0] + + noa = len(positions0) + + positions_2 = deepcopy(positions0) + positions_2[0][d1-1] = delt + positions += [positions_2] + + positions_3 = deepcopy(positions[0]) + positions_3[0][d1-1] = -delt + positions += [positions_3] + + test_struc = [] + for i in range(3): + test_struc += [struc.Structure(cell, species_1, positions[i])] + + env_0 = [] + env_p = [] + env_m = [] + for i in range(noa): + env_0 += [env.AtomicEnvironment(test_struc[0], i, cutoffs)] + env_p += [env.AtomicEnvironment(test_struc[1], i, cutoffs)] + env_m += [env.AtomicEnvironment(test_struc[2], i, cutoffs)] + + return [env_0, env_p, env_m] def generate_envs(cutoffs, delta): """ create environment with perturbation on diff --git a/tests/test_mc_mb_kernel.py b/tests/test_mc_mb_kernel.py index 5820d1bf7..07508780a 100644 --- a/tests/test_mc_mb_kernel.py +++ b/tests/test_mc_mb_kernel.py @@ -2,15 +2,15 @@ from copy import deepcopy import pytest import numpy as np +from numpy import isclose from numpy.random import random, randint from flare import env, struc, gp -import flare.kernels.mc_simple as mck -from flare.kernels.utils import str_to_kernel_set as stks +from flare.kernels.utils import str_to_kernel_set -from .fake_gp import generate_envs, another_env +from .fake_gp import generate_mb_envs -list_to_test = ['2mc', '3mc', '2+3mc', 'mbmc', '2+3+mbmc'] +list_to_test = ['2', '3', '2+3', 'mb', '2+3+mb'] def generate_hm(kernel_name): @@ -28,16 +28,19 @@ def test_force_en(kernel_name): energy kernel.""" cutoffs = np.ones(3)*1.2 - delta = 1e-8 - env1_1, env1_2, env1_3, \ - env1_2_1, env1_3_1, env1_2_2, env1_3_2, env2_1 \ - = another_env(cutoffs, delta) + delta = 1e-4 + tol = 1e-4 + cell = 1e7 * np.eye(3) # set hyperparameters d1 = 1 + + env1 = generate_mb_envs(cutoffs, cell, delta, d1) + env2 = generate_mb_envs(cutoffs, cell, delta, d1) + hyps = generate_hm(kernel_name) - _, __, en_kernel, force_en_kernel = stks(kernel_name) + _, __, en_kernel, force_en_kernel = str_to_kernel_set(kernel_name) nterm = 0 for term in ['2', '3', 'mb']: @@ -46,50 +49,39 @@ def test_force_en(kernel_name): kern_finite_diff = 0 if ('mb' in kernel_name): - _, __, enm_kernel, ___ = stks('mbmc') + _, __, enm_kernel, ___ = str_to_kernel_set('mb') mhyps = hyps[(nterm-1)*2:] - calc1 = enm_kernel(env1_2, env2_1, mhyps, cutoffs) - calc2 = enm_kernel(env1_3, env2_1, mhyps, cutoffs) - kern_finite_diff_00 = (calc1 - calc2) / (2 * delta) - - calc1 = enm_kernel(env1_2_1, env2_1, mhyps, cutoffs) - calc2 = enm_kernel(env1_3_1, env2_1, mhyps, cutoffs) - kern_finite_diff_10 = (calc1 - calc2) / (2 * delta) - - calc1 = enm_kernel(env1_2_2, env2_1, mhyps, cutoffs) - calc2 = enm_kernel(env1_3_2, env2_1, mhyps, cutoffs) - kern_finite_diff_20 = (calc1 - calc2) / (2 * delta) - - mb_diff = (kern_finite_diff_00 + - kern_finite_diff_10 + kern_finite_diff_20) - + calc = 0 + for i in range(len(env1[0])): + calc += enm_kernel(env1[2][i], env2[0][0], mhyps, cutoffs) + calc -= enm_kernel(env1[1][i], env2[0][0], mhyps, cutoffs) + mb_diff = calc / (2 * delta) kern_finite_diff += mb_diff if ('2' in kernel_name): nbond = 1 - _, __, en2_kernel, ___ = stks('2mc') - calc1 = en2_kernel(env1_2, env2_1, hyps[0:nbond * 2], cutoffs) - calc2 = en2_kernel(env1_1, env2_1, hyps[0:nbond * 2], cutoffs) - diff2b = (calc1 - calc2) / 2.0 / delta + _, __, en2_kernel, ___ = str_to_kernel_set('2') + calc1 = en2_kernel(env1[2][0], env2[0][0], hyps[0:nbond * 2], cutoffs) + calc2 = en2_kernel(env1[1][0], env2[0][0], hyps[0:nbond * 2], cutoffs) + diff2b = (calc1 - calc2) / 2.0 / 2.0 / delta kern_finite_diff += diff2b else: nbond = 0 if ('3' in kernel_name): - _, __, en3_kernel, ___ = stks('3mc') - calc1 = en3_kernel(env1_2, env2_1, hyps[nbond * 2:], cutoffs) - calc2 = en3_kernel(env1_1, env2_1, hyps[nbond * 2:], cutoffs) - diff3b = (calc1 - calc2) / 3.0 / delta + _, __, en3_kernel, ___ = str_to_kernel_set('3') + calc1 = en3_kernel(env1[2][0], env2[0][0], hyps[nbond * 2:], cutoffs) + calc2 = en3_kernel(env1[1][0], env2[0][0], hyps[nbond * 2:], cutoffs) + diff3b = (calc1 - calc2) / 2.0 / 3.0 / delta kern_finite_diff += diff3b - kern_analytical = force_en_kernel(env1_1, env2_1, d1, hyps, cutoffs) + kern_analytical = force_en_kernel(env1[0][0], env2[0][0], d1, hyps, cutoffs) print("\nforce_en", kernel_name, kern_finite_diff, kern_analytical) - tol = 1e-4 - assert (np.isclose(-kern_finite_diff, kern_analytical, atol=tol)) + assert (isclose(kern_finite_diff, kern_analytical, rtol=tol)) @pytest.mark.parametrize('kernel_name', list_to_test) @@ -97,365 +89,75 @@ def test_force(kernel_name): """Check that the analytical force kernel matches finite difference of energy kernel.""" - if ('mb' in kernel_name): - return 0 - # create env 1 + d1 = 1 + d2 = 2 + tol = 1e-4 + cell = 1e7 * np.eye(3) delta = 1e-5 - cutoffs = np.ones(3)*0.8 - env1_1, env1_2, env1_3, env2_1, env2_2, env2_3 \ - = generate_envs(cutoffs, delta) - - # set hyperparameters hyps = generate_hm(kernel_name) - d1 = 1 - d2 = 2 + cutoffs = np.ones(3)*1.2 + kernel, kg, en_kernel, fek = str_to_kernel_set(kernel_name, False) + args = (hyps, cutoffs) - kernel, _, en_kernel, ___ = stks(kernel_name) + np.random.seed(10) + env1 = generate_mb_envs(cutoffs, cell, delta, d1) + env2 = generate_mb_envs(cutoffs, cell, delta, d2) # check force kernel - calc1 = en_kernel(env1_2, env2_2, hyps, cutoffs) - calc2 = en_kernel(env1_3, env2_3, hyps, cutoffs) - calc3 = en_kernel(env1_2, env2_3, hyps, cutoffs) - calc4 = en_kernel(env1_3, env2_2, hyps, cutoffs) - - kern_finite_diff = (calc1 + calc2 - calc3 - calc4) / (4 * delta ** 2) - kern_analytical = kernel(env1_1, env2_1, - d1, d2, hyps, cutoffs) + if ('mb' == kernel_name): + cal = 0 + for i in range(3): + for j in range(len(env1[0])): + cal += en_kernel(env1[1][i], env2[1][j], *args) + cal += en_kernel(env1[2][i], env2[2][j], *args) + cal -= en_kernel(env1[1][i], env2[2][j], *args) + cal -= en_kernel(env1[2][i], env2[1][j], *args) + kern_finite_diff = cal / (4 * delta ** 2) + elif ('mb' not in kernel_name): + calc1 = en_kernel(env1[1][0], env2[1][0], *args) + calc2 = en_kernel(env1[2][0], env2[2][0], *args) + calc3 = en_kernel(env1[1][0], env2[2][0], *args) + calc4 = en_kernel(env1[2][0], env2[1][0], *args) + kern_finite_diff = (calc1 + calc2 - calc3 - calc4) / (4*delta**2) + else: + return + + kern_analytical = kernel(env1[0][0], env2[0][0], + d1, d2, *args) tol = 1e-4 - assert (np.isclose(kern_finite_diff, kern_analytical, atol=tol)) + assert(isclose(kern_finite_diff, kern_analytical, rtol=tol)) @pytest.mark.parametrize('kernel_name', list_to_test) def test_hyps_grad(kernel_name): - if ('mb' in kernel_name): - return 0 - delta = 1e-8 - cutoffs = np.array([1, 1, 1]) - env1_1, env1_2, env1_3, env2_1, env2_2, env2_3 = generate_envs( - cutoffs, delta) + cutoffs = np.ones(3)*1.2 hyps = generate_hm(kernel_name) d1 = randint(1, 3) d2 = randint(1, 3) print("hyps", hyps) + cell = 1e7 * np.eye(3) - kernel, kernel_grad, _, _ = stks(kernel_name, False) + env1 = generate_mb_envs(cutoffs, cell, 0, d1)[0][0] + env2 = generate_mb_envs(cutoffs, cell, 0, d2)[0][0] - grad_test = kernel_grad(env1_1, env2_1, + kernel, kernel_grad, _, _ = str_to_kernel_set(kernel_name, False) + + grad_test = kernel_grad(env1, env2, d1, d2, hyps, cutoffs) - print("analytical gradients", grad_test) tol = 1e-4 - original = kernel(env1_1, env2_1, d1, d2, + original = kernel(env1, env2, d1, d2, hyps, cutoffs) for i in range(len(hyps)-1): newhyps = np.copy(hyps) newhyps[i] += delta - hgrad = (kernel(env1_1, env2_1, d1, d2, newhyps, + hgrad = (kernel(env1, env2, d1, d2, newhyps, cutoffs) - original)/delta print("numerical gradients", hgrad) - assert(np.isclose(grad_test[1][i], hgrad, atol=tol)) - -# ----------------------------------------------------------------------------- -# test many body kernels -# ----------------------------------------------------------------------------- - - -def test_many_body_force(): - """Check that the analytical force kernel matches finite difference of - energy kernel.""" - - # create env 1 - delt = 1e-5 - cell = 1e7 * np.eye(3) - - - d1 = randint(1, 3) - d2 = randint(1, 3) - - cutoffs = np.array([1.2, 1.2, 1.2]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([0., 1., 0.]) + 0.1 * - np.array([random(), random(), random()]), - np.array([1., 0., 0.]) + 0.1 * - np.array([random(), random(), random()]), - np.array([1., 1., 0.]) + 0.1 * np.array([random(), random(), random()])] - - positions_2 = deepcopy(positions_1) - positions_2[0][d1-1] = delt - - positions_3 = deepcopy(positions_1) - positions_3[0][d1-1] = -delt - - species_1 = [randint(1, 2) for i in range(4)] - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - test_structure_3 = struc.Structure(cell, species_1, positions_3) - - env1_1_0 = env.AtomicEnvironment(test_structure_1, 0, cutoffs) - - env1_2_0 = env.AtomicEnvironment(test_structure_2, 0, cutoffs) - env1_2_1 = env.AtomicEnvironment(test_structure_2, 1, cutoffs) - env1_2_2 = env.AtomicEnvironment(test_structure_2, 2, cutoffs) - - env1_3_0 = env.AtomicEnvironment(test_structure_3, 0, cutoffs) - env1_3_1 = env.AtomicEnvironment(test_structure_3, 1, cutoffs) - env1_3_2 = env.AtomicEnvironment(test_structure_3, 2, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([0., 1., 0.]) + 0.1 * - np.array([random(), random(), random()]), - np.array([1., 0., 0.]) + 0.1 * - np.array([random(), random(), random()]), - np.array([1., 1., 0.]) + 0.1 * np.array([random(), random(), random()])] - - positions_2 = deepcopy(positions_1) - positions_2[0][d2-1] = delt - positions_3 = deepcopy(positions_1) - positions_3[0][d2-1] = -delt - - species_2 = [randint(1, 2) for i in range(4)] - test_structure_1 = struc.Structure(cell, species_2, positions_1) - test_structure_2 = struc.Structure(cell, species_2, positions_2) - test_structure_3 = struc.Structure(cell, species_2, positions_3) - - env2_1_0 = env.AtomicEnvironment(test_structure_1, 0, cutoffs) - - env2_2_0 = env.AtomicEnvironment(test_structure_2, 0, cutoffs) - env2_2_1 = env.AtomicEnvironment(test_structure_2, 1, cutoffs) - env2_2_2 = env.AtomicEnvironment(test_structure_2, 2, cutoffs) - - env2_3_0 = env.AtomicEnvironment(test_structure_3, 0, cutoffs) - env2_3_1 = env.AtomicEnvironment(test_structure_3, 1, cutoffs) - env2_3_2 = env.AtomicEnvironment(test_structure_3, 2, cutoffs) - - sig = random()*2 - ls = random()*2 - - hyps = np.array([sig, ls]) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_0, env2_2_0, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_0, env2_3_0, hyps, cutoffs) - calc3 = mck.many_body_mc_en(env1_2_0, env2_3_0, hyps, cutoffs) - calc4 = mck.many_body_mc_en(env1_3_0, env2_2_0, hyps, cutoffs) - - kern_finite_diff_00 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_0, env2_2_1, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_0, env2_3_1, hyps, cutoffs) - calc3 = mck.many_body_mc_en(env1_2_0, env2_3_1, hyps, cutoffs) - calc4 = mck.many_body_mc_en(env1_3_0, env2_2_1, hyps, cutoffs) - - kern_finite_diff_01 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_0, env2_2_2, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_0, env2_3_2, hyps, cutoffs) - calc3 = mck.many_body_mc_en(env1_2_0, env2_3_2, hyps, cutoffs) - calc4 = mck.many_body_mc_en(env1_3_0, env2_2_2, hyps, cutoffs) - - kern_finite_diff_02 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_1, env2_2_0, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_1, env2_3_0, hyps, cutoffs) - calc3 = mck.many_body_mc_en(env1_2_1, env2_3_0, hyps, cutoffs) - calc4 = mck.many_body_mc_en(env1_3_1, env2_2_0, hyps, cutoffs) - - kern_finite_diff_10 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_1, env2_2_1, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_1, env2_3_1, hyps, cutoffs) - calc3 = mck.many_body_mc_en(env1_2_1, env2_3_1, hyps, cutoffs) - calc4 = mck.many_body_mc_en(env1_3_1, env2_2_1, hyps, cutoffs) - - kern_finite_diff_11 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_1, env2_2_2, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_1, env2_3_2, hyps, cutoffs) - calc3 = mck.many_body_mc_en(env1_2_1, env2_3_2, hyps, cutoffs) - calc4 = mck.many_body_mc_en(env1_3_1, env2_2_2, hyps, cutoffs) - - kern_finite_diff_12 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_2, env2_2_0, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_2, env2_3_0, hyps, cutoffs) - calc3 = mck.many_body_mc_en(env1_2_2, env2_3_0, hyps, cutoffs) - calc4 = mck.many_body_mc_en(env1_3_2, env2_2_0, hyps, cutoffs) - - kern_finite_diff_20 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_2, env2_2_1, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_2, env2_3_1, hyps, cutoffs) - calc3 = mck.many_body_mc_en(env1_2_2, env2_3_1, hyps, cutoffs) - calc4 = mck.many_body_mc_en(env1_3_2, env2_2_1, hyps, cutoffs) - - kern_finite_diff_21 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_2, env2_2_2, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_2, env2_3_2, hyps, cutoffs) - calc3 = mck.many_body_mc_en(env1_2_2, env2_3_2, hyps, cutoffs) - calc4 = mck.many_body_mc_en(env1_3_2, env2_2_2, hyps, cutoffs) - - kern_finite_diff_22 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - kern_finite_diff = (kern_finite_diff_00 + kern_finite_diff_01 + kern_finite_diff_02 + - kern_finite_diff_10 + kern_finite_diff_11 + kern_finite_diff_12 + - kern_finite_diff_20 + kern_finite_diff_21 + kern_finite_diff_22) - - kern_analytical = mck.many_body_mc(env1_1_0, env2_1_0, - d1, d2, hyps, cutoffs) - - tol = 1e-4 - - assert (np.isclose(kern_finite_diff, kern_analytical, atol=tol)) - - -def test_many_body_grad(): - # create env 1 - cell = 1e7 * np.eye(3) - cutoffs = np.array([2, 2, 2]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_1 = [randint(1, 2) for i in range(3)] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - env1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_2 = [randint(1, 2) for i in range(3)] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - env2 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - - sig = random()*2 - ls = random()*2 - - d1 = randint(1, 3) - d2 = randint(1, 3) - - hyps = np.array([sig, ls]) - - grad_test = mck.many_body_mc_grad(env1, env2, d1, d2, hyps, cutoffs) - - delta = 1e-8 - new_sig = sig + delta - new_ls = ls + delta - - sig_derv_brute = (mck.many_body_mc(env1, env2, d1, d2, - np.array([new_sig, ls]), - cutoffs) - - mck.many_body_mc(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - l_derv_brute = (mck.many_body_mc(env1, env2, d1, d2, - np.array([sig, new_ls]), - cutoffs) - - mck.many_body_mc(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - tol = 1e-4 - assert (np.isclose(grad_test[1][0], sig_derv_brute, atol=tol)) - assert (np.isclose(grad_test[1][1], l_derv_brute, atol=tol)) - - -def test_many_body_force_en(): - """Check that the analytical force-energy kernel matches finite difference of - energy kernel.""" - - # create env 1 - delt = 1e-5 - cell = 1e7 * np.eye(3) - cutoffs = np.array([1.2, 1.2, 1.2]) - d1 = randint(1, 3) - - positions_1 = [np.array([0., 0., 0.]), - np.array([0., 1., 0.]) + 0.1 * - np.array([random(), random(), random()]), - np.array([1., 0., 0.]) + 0.1 * - np.array([random(), random(), random()]), - np.array([1., 1., 0.]) + 0.1 * np.array([random(), random(), random()])] - - positions_2 = deepcopy(positions_1) - positions_2[0][d1-1] = delt - - positions_3 = deepcopy(positions_1) - positions_3[0][d1-1] = -delt - - species_1 = [randint(1, 2) for i in range(4)] - - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - test_structure_3 = struc.Structure(cell, species_1, positions_3) - - env1_1_0 = env.AtomicEnvironment(test_structure_1, 0, cutoffs) - - env1_2_0 = env.AtomicEnvironment(test_structure_2, 0, cutoffs) - env1_2_1 = env.AtomicEnvironment(test_structure_2, 1, cutoffs) - env1_2_2 = env.AtomicEnvironment(test_structure_2, 2, cutoffs) - - env1_3_0 = env.AtomicEnvironment(test_structure_3, 0, cutoffs) - env1_3_1 = env.AtomicEnvironment(test_structure_3, 1, cutoffs) - env1_3_2 = env.AtomicEnvironment(test_structure_3, 2, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([0., 1., 0.]) + 0.1 * - np.array([random(), random(), random()]), - np.array([1., 0., 0.]) + 0.1 * - np.array([random(), random(), random()]), - np.array([1., 1., 0.]) + 0.1 * np.array([random(), random(), random()])] - - species_2 = [randint(1, 2) for i in range(4)] - - test_structure_1 = struc.Structure(cell, species_2, positions_1) - - env2_1_0 = env.AtomicEnvironment(test_structure_1, 0, cutoffs) - - sig = random()*2 - ls = random()*2 - - hyps = np.array([sig, ls]) - - # check force kernel - calc1 = mck.many_body_mc_en(env1_2_0, env2_1_0, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_0, env2_1_0, hyps, cutoffs) - kern_finite_diff_00 = (calc1 - calc2) / (2 * delt) - - calc1 = mck.many_body_mc_en(env1_2_1, env2_1_0, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_1, env2_1_0, hyps, cutoffs) - kern_finite_diff_10 = (calc1 - calc2) / (2 * delt) - - calc1 = mck.many_body_mc_en(env1_2_2, env2_1_0, hyps, cutoffs) - calc2 = mck.many_body_mc_en(env1_3_2, env2_1_0, hyps, cutoffs) - kern_finite_diff_20 = (calc1 - calc2) / (2 * delt) - - kern_finite_diff = -(kern_finite_diff_00 + - kern_finite_diff_10 + kern_finite_diff_20) - - kern_analytical = mck.many_body_mc_force_en(env1_1_0, env2_1_0, - d1, hyps, cutoffs) - - tol = 1e-4 - - assert (np.isclose(kern_finite_diff, kern_analytical, atol=tol)) + print("analytical gradients", grad_test[1][i]) + assert(isclose(grad_test[1][i], hgrad, rtol=tol)) From c51ddf5343c8feb5bf328533a3414e6e053cadff Mon Sep 17 00:00:00 2001 From: nw13slx Date: Thu, 23 Apr 2020 15:19:41 -0400 Subject: [PATCH 3/6] merge sc and mc unit test --- flare/kernels/sc.py | 8 +- tests/test_kernels.py | 956 ------------------------------------- tests/test_mc_mb_kernel.py | 38 +- 3 files changed, 26 insertions(+), 976 deletions(-) delete mode 100644 tests/test_kernels.py diff --git a/flare/kernels/sc.py b/flare/kernels/sc.py index 78cb474d1..66c3502c1 100644 --- a/flare/kernels/sc.py +++ b/flare/kernels/sc.py @@ -1848,10 +1848,10 @@ def mb_grad_helper_ls(q1, q2, qi, qj, sig, ls): '2+3_grad': two_plus_three_body_grad, '2+3_en': two_plus_three_en, '2+3_force_en': two_plus_three_force_en, - 'many_body': many_body, - 'many_body_en': many_body_en, - 'many_body_grad': many_body_grad, - 'many_body_force_en': many_body_en, + 'many': many_body, + 'many_en': many_body_en, + 'many_grad': many_body_grad, + 'many_force_en': many_body_force_en, 'two_plus_three_plus_many_body': two_plus_three_plus_many_body, 'two_plus_three_plus_many_body_grad': two_plus_three_plus_many_body_grad, 'two_plus_three_plus_many_body_en': two_plus_three_plus_many_body_en, diff --git a/tests/test_kernels.py b/tests/test_kernels.py deleted file mode 100644 index 9d9fb3c7e..000000000 --- a/tests/test_kernels.py +++ /dev/null @@ -1,956 +0,0 @@ -import pytest -import numpy as np -import sys -from random import random, randint -from copy import deepcopy -from flare import env, gp, struc -import flare.kernels.sc as en - -from flare.kernels.utils import from_mask_to_args, from_grad_to_mask - -# ----------------------------------------------------------------------------- -# test two plus three body kernels -# ----------------------------------------------------------------------------- - -# TODO: fix this test to properly account for factors of 2 and 3 -def test_two_plus_three_body_force_en(): - """Check that the analytical force/en kernel matches finite difference of - energy kernel.""" - - # create env 1 - delt = 1e-8 - cell = np.eye(3) - cutoffs = np.array([1, 1]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - positions_2 = deepcopy(positions_1) - positions_2[0][0] = delt - - species_1 = [1, 2, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - - env1_1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - env1_2 = env.AtomicEnvironment(test_structure_2, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_2 = [1, 1, 2] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - env2 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - - # set hyperparameters - sig1 = random() - ls1 = random() - sig2 = random() - ls2 = random() - d1 = 1 - - hyps = np.array([sig1, ls1, sig2, ls2]) - - # check force kernel - calc1 = en.two_body_en(env1_2, env2, hyps[0:2], cutoffs) - calc2 = en.two_body_en(env1_1, env2, hyps[0:2], cutoffs) - calc3 = en.three_body_en(env1_2, env2, hyps[2:4], cutoffs) - calc4 = en.three_body_en(env1_1, env2, hyps[2:4], cutoffs) - - kern_finite_diff = (calc1 - calc2) / (2 * delt) + \ - (calc3 - calc4) / (3 * delt) - kern_analytical = \ - en.two_plus_three_force_en(env1_1, env2, d1, hyps, cutoffs) - - tol = 1e-4 - assert(np.isclose(-kern_finite_diff, kern_analytical, atol=tol)) - - -def test_two_plus_three_body_force(): - """Check that the analytical force kernel matches finite difference of - energy kernel.""" - - # create env 1 - delt = 1e-5 - cell = np.eye(3) - cutoffs = np.array([1, 0.9]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - positions_2 = deepcopy(positions_1) - positions_2[0][0] = delt - - positions_3 = deepcopy(positions_1) - positions_3[0][0] = -delt - - species_1 = [1, 2, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - test_structure_3 = struc.Structure(cell, species_1, positions_3) - - env1_1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - env1_2 = env.AtomicEnvironment(test_structure_2, atom_1, cutoffs) - env1_3 = env.AtomicEnvironment(test_structure_3, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - positions_2 = deepcopy(positions_1) - positions_2[0][1] = delt - positions_3 = deepcopy(positions_1) - positions_3[0][1] = -delt - - species_2 = [1, 1, 2] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - test_structure_2 = struc.Structure(cell, species_2, positions_2) - test_structure_3 = struc.Structure(cell, species_2, positions_3) - - env2_1 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - env2_2 = env.AtomicEnvironment(test_structure_2, atom_2, cutoffs) - env2_3 = env.AtomicEnvironment(test_structure_3, atom_2, cutoffs) - - # set hyperparameters - sig1 = random() - ls1 = random() - sig2 = random() - ls2 = random() - d1 = 1 - d2 = 2 - - hyps = np.array([sig1, ls1, sig2, ls2]) - - # check force kernel - calc1 = en.two_plus_three_en(env1_2, env2_2, hyps, cutoffs) - calc2 = en.two_plus_three_en(env1_3, env2_3, hyps, cutoffs) - calc3 = en.two_plus_three_en(env1_2, env2_3, hyps, cutoffs) - calc4 = en.two_plus_three_en(env1_3, env2_2, hyps, cutoffs) - - kern_finite_diff = (calc1 + calc2 - calc3 - calc4) / (4*delt**2) - kern_analytical = en.two_plus_three_body(env1_1, env2_1, - d1, d2, hyps, cutoffs) - - tol = 1e-4 - assert(np.isclose(kern_finite_diff, kern_analytical, atol=tol)) - - -def test_two_plus_three_body_grad(): - # create env 1 - cell = np.eye(3) - cutoffs = np.array([1, 1]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_1 = [1, 2, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - env1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_2 = [1, 1, 2] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - env2 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - - # set hyperparameters - sig1 = random() - ls1 = random() - sig2 = random() - ls2 = random() - - d1 = randint(1, 3) - d2 = randint(1, 3) - - delta = 1e-8 - - hyps = np.array([sig1, ls1, sig2, ls2]) - hyps1 = np.array([sig1+delta, ls1, sig2, ls2]) - hyps2 = np.array([sig1, ls1+delta, sig2, ls2]) - hyps3 = np.array([sig1, ls1, sig2+delta, ls2]) - hyps4 = np.array([sig1, ls1, sig2, ls2+delta]) - - grad_test = en.two_plus_three_body_grad(env1, env2, d1, d2, hyps, cutoffs) - - sig1_derv_brute = (en.two_plus_three_body(env1, env2, d1, d2, - hyps1, cutoffs) - - en.two_plus_three_body(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - l1_derv_brute = \ - (en.two_plus_three_body(env1, env2, d1, d2, hyps2, cutoffs) - - en.two_plus_three_body(env1, env2, d1, d2, hyps, cutoffs)) / delta - - sig2_derv_brute = \ - (en.two_plus_three_body(env1, env2, d1, d2, - hyps3, cutoffs) - - en.two_plus_three_body(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - l2_derv_brute = \ - (en.two_plus_three_body(env1, env2, d1, d2, hyps4, cutoffs) - - en.two_plus_three_body(env1, env2, d1, d2, hyps, cutoffs)) / delta - - tol = 1e-4 - assert(np.isclose(grad_test[1][0], sig1_derv_brute, atol=tol)) - assert(np.isclose(grad_test[1][1], l1_derv_brute, atol=tol)) - assert(np.isclose(grad_test[1][2], sig2_derv_brute, atol=tol)) - assert(np.isclose(grad_test[1][3], l2_derv_brute, atol=tol)) - - -# ----------------------------------------------------------------------------- -# test two body kernels -# ----------------------------------------------------------------------------- - -def test_two_body_force_en(): - """Check that the analytical force/en kernel matches finite difference of - energy kernel.""" - - # create env 1 - delt = 1e-8 - cell = np.eye(3) - cutoffs = np.array([1]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - positions_2 = deepcopy(positions_1) - positions_2[0][0] = delt - - species_1 = [1, 2, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - - env1_1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - env1_2 = env.AtomicEnvironment(test_structure_2, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_2 = [1, 1, 2] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - env2 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - - sig = random() - ls = random() - d1 = 1 - - hyps = np.array([sig, ls]) - - # check force kernel - calc1 = en.two_body_en(env1_2, env2, hyps, cutoffs) - calc2 = en.two_body_en(env1_1, env2, hyps, cutoffs) - - kern_finite_diff = (calc1 - calc2) / delt - kern_analytical = en.two_body_force_en(env1_1, env2, d1, hyps, cutoffs) - - tol = 1e-4 - assert(np.isclose(-kern_finite_diff/2, kern_analytical, atol=tol)) - - -def test_two_body_force(): - """Check that the analytical force kernel matches finite difference of - energy kernel.""" - - # create env 1 - delt = 1e-5 - cell = np.eye(3) - cutoffs = np.array([1, 1]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - positions_2 = deepcopy(positions_1) - positions_2[0][0] = delt - - positions_3 = deepcopy(positions_1) - positions_3[0][0] = -delt - - species_1 = [1, 2, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - test_structure_3 = struc.Structure(cell, species_1, positions_3) - - env1_1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - env1_2 = env.AtomicEnvironment(test_structure_2, atom_1, cutoffs) - env1_3 = env.AtomicEnvironment(test_structure_3, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - positions_2 = deepcopy(positions_1) - positions_2[0][1] = delt - positions_3 = deepcopy(positions_1) - positions_3[0][1] = -delt - - species_2 = [1, 1, 2] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - test_structure_2 = struc.Structure(cell, species_2, positions_2) - test_structure_3 = struc.Structure(cell, species_2, positions_3) - - env2_1 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - env2_2 = env.AtomicEnvironment(test_structure_2, atom_2, cutoffs) - env2_3 = env.AtomicEnvironment(test_structure_3, atom_2, cutoffs) - - sig = 1 - ls = 0.1 - d1 = 1 - d2 = 2 - - hyps = np.array([sig, ls]) - - # check force kernel - calc1 = en.two_body_en(env1_2, env2_2, hyps, cutoffs) - calc2 = en.two_body_en(env1_3, env2_3, hyps, cutoffs) - calc3 = en.two_body_en(env1_2, env2_3, hyps, cutoffs) - calc4 = en.two_body_en(env1_3, env2_2, hyps, cutoffs) - - kern_finite_diff = (calc1 + calc2 - calc3 - calc4) / (4*delt**2) - kern_analytical = en.two_body(env1_1, env2_1, - d1, d2, hyps, cutoffs) - - tol = 1e-4 - assert(np.isclose(kern_finite_diff, kern_analytical, atol=tol)) - - -def test_two_body_grad(): - # create env 1 - cell = np.eye(3) - cutoffs = np.array([1, 1]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_1 = [1, 2, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - env1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_2 = [1, 1, 2] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - env2 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - - sig = random() - ls = random() - d1 = randint(1, 3) - d2 = randint(1, 3) - - hyps = np.array([sig, ls]) - - grad_test = en.two_body_grad(env1, env2, d1, d2, hyps, cutoffs) - - delta = 1e-8 - new_sig = sig + delta - new_ls = ls + delta - - sig_derv_brute = (en.two_body(env1, env2, d1, d2, - np.array([new_sig, ls]), - cutoffs) - - en.two_body(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - l_derv_brute = (en.two_body(env1, env2, d1, d2, - np.array([sig, new_ls]), - cutoffs) - - en.two_body(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - tol = 1e-4 - assert(np.isclose(grad_test[1][0], sig_derv_brute, atol=tol)) - assert(np.isclose(grad_test[1][1], l_derv_brute, atol=tol)) - - -# ----------------------------------------------------------------------------- -# test three body kernels -# ----------------------------------------------------------------------------- - - -def test_three_body_force_en(): - """Check that the analytical force/en kernel matches finite difference of - energy kernel.""" - - # create env 1 - delt = 1e-8 - cell = np.eye(3) - cutoffs = np.array([1, 1]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - positions_2 = deepcopy(positions_1) - positions_2[0][0] = delt - - species_1 = [1, 2, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - - env1_1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - env1_2 = env.AtomicEnvironment(test_structure_2, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_2 = [1, 1, 2] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - env2 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - - sig = random() - ls = random() - d1 = 1 - - hyps = np.array([sig, ls]) - - # check force kernel - calc1 = en.three_body_en(env1_2, env2, hyps, cutoffs) - calc2 = en.three_body_en(env1_1, env2, hyps, cutoffs) - - kern_finite_diff = (calc1 - calc2) / delt - kern_analytical = en.three_body_force_en(env1_1, env2, d1, hyps, cutoffs) - - tol = 1e-4 - assert(np.isclose(-kern_finite_diff/3, kern_analytical, atol=tol)) - - -def test_three_body_force(): - """Check that the analytical force kernel matches finite difference of - energy kernel.""" - - # create env 1 - delt = 1e-5 - cell = np.eye(3) - cutoffs = np.array([1, 1]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - positions_2 = deepcopy(positions_1) - positions_2[0][0] = delt - - positions_3 = deepcopy(positions_1) - positions_3[0][0] = -delt - - species_1 = [1, 2, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - test_structure_3 = struc.Structure(cell, species_1, positions_3) - - env1_1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - env1_2 = env.AtomicEnvironment(test_structure_2, atom_1, cutoffs) - env1_3 = env.AtomicEnvironment(test_structure_3, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - positions_2 = deepcopy(positions_1) - positions_2[0][1] = delt - positions_3 = deepcopy(positions_1) - positions_3[0][1] = -delt - - species_2 = [1, 1, 2] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - test_structure_2 = struc.Structure(cell, species_2, positions_2) - test_structure_3 = struc.Structure(cell, species_2, positions_3) - - env2_1 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - env2_2 = env.AtomicEnvironment(test_structure_2, atom_2, cutoffs) - env2_3 = env.AtomicEnvironment(test_structure_3, atom_2, cutoffs) - - sig = 1 - ls = 0.1 - d1 = 1 - d2 = 2 - - hyps = np.array([sig, ls]) - - # check force kernel - calc1 = en.three_body_en(env1_2, env2_2, hyps, cutoffs) - calc2 = en.three_body_en(env1_3, env2_3, hyps, cutoffs) - calc3 = en.three_body_en(env1_2, env2_3, hyps, cutoffs) - calc4 = en.three_body_en(env1_3, env2_2, hyps, cutoffs) - - kern_finite_diff = (calc1 + calc2 - calc3 - calc4) / (4*delt**2) - kern_analytical = en.three_body(env1_1, env2_1, - d1, d2, hyps, cutoffs) - - tol = 1e-4 - assert(np.isclose(kern_finite_diff, kern_analytical, atol=tol)) - - -def test_three_body_grad(): - # create env 1 - cell = np.eye(3) - cutoffs = np.array([1, 1]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_1 = [1, 2, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - env1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_2 = [1, 1, 2] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - env2 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - - sig = random() - ls = random() - d1 = randint(1, 3) - d2 = randint(1, 3) - - hyps = np.array([sig, ls]) - - grad_test = en.three_body_grad(env1, env2, d1, d2, hyps, cutoffs) - - delta = 1e-8 - new_sig = sig + delta - new_ls = ls + delta - - sig_derv_brute = (en.three_body(env1, env2, d1, d2, - np.array([new_sig, ls]), - cutoffs) - - en.three_body(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - l_derv_brute = (en.three_body(env1, env2, d1, d2, - np.array([sig, new_ls]), - cutoffs) - - en.three_body(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - tol = 1e-4 - assert(np.isclose(grad_test[1][0], sig_derv_brute, atol=tol)) - assert(np.isclose(grad_test[1][1], l_derv_brute, atol=tol)) - - -def test_masked_hyperparameter_function(): - """ - Test simple input permutations for the from_mask_to_args function - :return: - """ - - cutoffs = [3, 4] - - # Standard sig2, ls2, sig3, ls3, noise hyp array - with pytest.raises(NameError): - from_mask_to_args([], {}, cutoffs) - # ----------------------- - # Test simple input cases - # ----------------------- - hyps_mask = {'nbond': 1, 'nspec':1, 'spec_mask':np.zeros(118)} - hyps = [1,2,5] - c, ns, sm, n2b, bm, n3b, tm, sig2, ls2, sig3, ls3 \ - = from_mask_to_args(hyps, hyps_mask, cutoffs) - assert (np.equal(sig2, [1]).all()) - assert (np.equal(ls2, [2]).all()) - - hyps = [3,4,5] - hyps_mask = {'ntriplet': 1, 'nspec':1, 'spec_mask':np.zeros(118)} - c, ns, sm, n2b, bm, n3b, tm, sig2, ls2, sig3, ls3 \ - = from_mask_to_args(hyps, hyps_mask, cutoffs) - assert (np.equal(sig3, [3]).all()) - assert (np.equal(ls3, [4]).all()) - - hyps = [1, 2, 3, 4, 5] - hyps_mask = {'nbond': 1, 'ntriplet':1, - 'nspec':1, 'spec_mask':np.zeros(118)} - c, ns, sm, n2b, bm, n3b, tm, sig2, ls2, sig3, ls3 \ - = from_mask_to_args(hyps, hyps_mask, cutoffs) - assert (np.equal(sig2, [1]).all()) - assert (np.equal(ls2, [2]).all()) - assert (np.equal(sig3, [3]).all()) - assert (np.equal(ls3, [4]).all()) - - hyps = [1, 2, 3, 4, 5] - hyps_mask['map']=[0, 1, 2, 3, 4] - hyps_mask['original'] = [1, 2, 3, 4, 5, 6] - c, ns, sm, n2b, bm, n3b, tm, sig2, ls2, sig3, ls3 \ - = from_mask_to_args(hyps, hyps_mask, cutoffs) - assert (np.equal(sig2, [1]).all()) - assert (np.equal(ls2, [2]).all()) - assert (np.equal(sig3, [3]).all()) - assert (np.equal(ls3, [4]).all()) - - # ----------------------- - # Test simple 2+3 body input case - # ----------------------- - - # Hyps : sig21, sig22, ls21, ls22, sig31, sig32, ls31, ls32, noise - hyps = [1.1,1.2, 2.1, 2.2, 3.1, 3.2, 4.1, 4.2, 5] - hyps_mask = {'nbond': 2, 'ntriplet': 2, - 'nspec':1, 'spec_mask':np.zeros(118)} - c, ns, sm, n2b, bm, n3b, tm, sig2, ls2, sig3, ls3 \ - = from_mask_to_args(hyps, hyps_mask, cutoffs) - assert (np.equal(sig2, [1.1, 1.2]).all()) - assert (np.equal(ls2, [2.1, 2.2]).all()) - assert (np.equal(sig3, [3.1, 3.2]).all()) - assert (np.equal(ls3, [4.1, 4.2]).all()) - - -def test_grad_mask_function(): - """ - Test simple permutations for the from_grad_to_mask function - :return: - """ - - grad = 'A' - - assert from_grad_to_mask(grad,hyps_mask={}) == 'A' - - # Test map for a standard slate of hyperparameters - grad = [1, 2, 3, 4, 5] - hyps_mask = {'map': [0, 3]} - assert (from_grad_to_mask(grad, hyps_mask) == [1, 4]).all() - - - # Test map when noise variance is included - grad = [1, 2, 3, 4, 5] - hyps_mask = {'map': [0, 3, 4]} - assert (from_grad_to_mask(grad, hyps_mask) == [1, 4, 5]).all() - - # Test map if the last parameter is not sigma_noise - grad = [1, 2, 3, 4, 5] - hyps_mask = {'map': [0, 3, 5]} - assert (from_grad_to_mask(grad, hyps_mask) == [1, 4]).all() - - -# ----------------------------------------------------------------------------- -# test many body kernels -# ----------------------------------------------------------------------------- - - -def test_many_body_force(): - """Check that the analytical force kernel matches finite difference of - energy kernel.""" - - # create env 1 - delt = 1e-5 - cell = 10.0 * np.eye(3) - cutoffs = np.array([1.2, 1.2, 1.2]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([0., 1., 0.]) + 0.1 * np.array([random(), random(), random()]), - np.array([1., 0., 0.]) + 0.1 * np.array([random(), random(), random()]), - np.array([1., 1., 0.]) + 0.1 * np.array([random(), random(), random()])] - - positions_2 = deepcopy(positions_1) - positions_2[0][0] = delt - - positions_3 = deepcopy(positions_1) - positions_3[0][0] = -delt - - species_1 = [1, 1, 1, 1] - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - test_structure_3 = struc.Structure(cell, species_1, positions_3) - - env1_1_0 = env.AtomicEnvironment(test_structure_1, 0, cutoffs) - - env1_2_0 = env.AtomicEnvironment(test_structure_2, 0, cutoffs) - env1_2_1 = env.AtomicEnvironment(test_structure_2, 1, cutoffs) - env1_2_2 = env.AtomicEnvironment(test_structure_2, 2, cutoffs) - - env1_3_0 = env.AtomicEnvironment(test_structure_3, 0, cutoffs) - env1_3_1 = env.AtomicEnvironment(test_structure_3, 1, cutoffs) - env1_3_2 = env.AtomicEnvironment(test_structure_3, 2, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([0., 1., 0.]) + 0.1 * np.array([random(), random(), random()]), - np.array([1., 0., 0.]) + 0.1 * np.array([random(), random(), random()]), - np.array([1., 1., 0.]) + 0.1 * np.array([random(), random(), random()])] - - positions_2 = deepcopy(positions_1) - positions_2[0][1] = delt - positions_3 = deepcopy(positions_1) - positions_3[0][1] = -delt - - species_2 = [1, 1, 1, 1] - test_structure_1 = struc.Structure(cell, species_2, positions_1) - test_structure_2 = struc.Structure(cell, species_2, positions_2) - test_structure_3 = struc.Structure(cell, species_2, positions_3) - - env2_1_0 = env.AtomicEnvironment(test_structure_1, 0, cutoffs) - - env2_2_0 = env.AtomicEnvironment(test_structure_2, 0, cutoffs) - env2_2_1 = env.AtomicEnvironment(test_structure_2, 1, cutoffs) - env2_2_2 = env.AtomicEnvironment(test_structure_2, 2, cutoffs) - - env2_3_0 = env.AtomicEnvironment(test_structure_3, 0, cutoffs) - env2_3_1 = env.AtomicEnvironment(test_structure_3, 1, cutoffs) - env2_3_2 = env.AtomicEnvironment(test_structure_3, 2, cutoffs) - - sig = random() - ls = random() - d1 = 1 - d2 = 2 - - hyps = np.array([sig, ls]) - - # check force kernel - calc1 = en.many_body_en(env1_2_0, env2_2_0, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_0, env2_3_0, hyps, cutoffs) - calc3 = en.many_body_en(env1_2_0, env2_3_0, hyps, cutoffs) - calc4 = en.many_body_en(env1_3_0, env2_2_0, hyps, cutoffs) - - kern_finite_diff_00 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = en.many_body_en(env1_2_0, env2_2_1, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_0, env2_3_1, hyps, cutoffs) - calc3 = en.many_body_en(env1_2_0, env2_3_1, hyps, cutoffs) - calc4 = en.many_body_en(env1_3_0, env2_2_1, hyps, cutoffs) - - kern_finite_diff_01 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = en.many_body_en(env1_2_0, env2_2_2, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_0, env2_3_2, hyps, cutoffs) - calc3 = en.many_body_en(env1_2_0, env2_3_2, hyps, cutoffs) - calc4 = en.many_body_en(env1_3_0, env2_2_2, hyps, cutoffs) - - kern_finite_diff_02 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = en.many_body_en(env1_2_1, env2_2_0, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_1, env2_3_0, hyps, cutoffs) - calc3 = en.many_body_en(env1_2_1, env2_3_0, hyps, cutoffs) - calc4 = en.many_body_en(env1_3_1, env2_2_0, hyps, cutoffs) - - kern_finite_diff_10 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = en.many_body_en(env1_2_1, env2_2_1, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_1, env2_3_1, hyps, cutoffs) - calc3 = en.many_body_en(env1_2_1, env2_3_1, hyps, cutoffs) - calc4 = en.many_body_en(env1_3_1, env2_2_1, hyps, cutoffs) - - kern_finite_diff_11 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = en.many_body_en(env1_2_1, env2_2_2, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_1, env2_3_2, hyps, cutoffs) - calc3 = en.many_body_en(env1_2_1, env2_3_2, hyps, cutoffs) - calc4 = en.many_body_en(env1_3_1, env2_2_2, hyps, cutoffs) - - kern_finite_diff_12 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = en.many_body_en(env1_2_2, env2_2_0, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_2, env2_3_0, hyps, cutoffs) - calc3 = en.many_body_en(env1_2_2, env2_3_0, hyps, cutoffs) - calc4 = en.many_body_en(env1_3_2, env2_2_0, hyps, cutoffs) - - kern_finite_diff_20 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = en.many_body_en(env1_2_2, env2_2_1, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_2, env2_3_1, hyps, cutoffs) - calc3 = en.many_body_en(env1_2_2, env2_3_1, hyps, cutoffs) - calc4 = en.many_body_en(env1_3_2, env2_2_1, hyps, cutoffs) - - kern_finite_diff_21 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - # check force kernel - calc1 = en.many_body_en(env1_2_2, env2_2_2, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_2, env2_3_2, hyps, cutoffs) - calc3 = en.many_body_en(env1_2_2, env2_3_2, hyps, cutoffs) - calc4 = en.many_body_en(env1_3_2, env2_2_2, hyps, cutoffs) - - kern_finite_diff_22 = (calc1 + calc2 - calc3 - calc4) / (4 * delt ** 2) - - kern_finite_diff = (kern_finite_diff_00 + kern_finite_diff_01 + kern_finite_diff_02 + - kern_finite_diff_10 + kern_finite_diff_11 + kern_finite_diff_12 + - kern_finite_diff_20 + kern_finite_diff_21 + kern_finite_diff_22) - - kern_analytical = en.many_body(env1_1_0, env2_1_0, - d1, d2, hyps, cutoffs) - - tol = 1e-4 - - assert (np.isclose(kern_finite_diff, kern_analytical, atol=tol)) - - -def test_many_body_force_en(): - """Check that the analytical force-energy kernel matches finite difference of - energy kernel.""" - - # TODO: why env1_1_1, env1_1_2 and env1_1_3 (and other variables) are never used? - # create env 1 - delt = 1e-5 - cell = 10.0 * np.eye(3) - cutoffs = np.array([1.2, 1.2, 1.2]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([0., 1., 0.]) + 0.1 * np.array([random(), random(), random()]), - np.array([1., 0., 0.]) + 0.1 * np.array([random(), random(), random()]), - np.array([1., 1., 0.]) + 0.1 * np.array([random(), random(), random()])] - - positions_2 = deepcopy(positions_1) - positions_2[0][0] = delt - - positions_3 = deepcopy(positions_1) - positions_3[0][0] = -delt - - species_1 = [1, 1, 1, 1] - - test_structure_1 = struc.Structure(cell, species_1, positions_1) - test_structure_2 = struc.Structure(cell, species_1, positions_2) - test_structure_3 = struc.Structure(cell, species_1, positions_3) - - env1_1_0 = env.AtomicEnvironment(test_structure_1, 0, cutoffs) - - env1_2_0 = env.AtomicEnvironment(test_structure_2, 0, cutoffs) - env1_2_1 = env.AtomicEnvironment(test_structure_2, 1, cutoffs) - env1_2_2 = env.AtomicEnvironment(test_structure_2, 2, cutoffs) - - env1_3_0 = env.AtomicEnvironment(test_structure_3, 0, cutoffs) - env1_3_1 = env.AtomicEnvironment(test_structure_3, 1, cutoffs) - env1_3_2 = env.AtomicEnvironment(test_structure_3, 2, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([0., 1., 0.]) + 0.1 * np.array([random(), random(), random()]), - np.array([1., 0., 0.]) + 0.1 * np.array([random(), random(), random()]), - np.array([1., 1., 0.]) + 0.1 * np.array([random(), random(), random()])] - - species_2 = [1, 1, 1, 1] - - test_structure_1 = struc.Structure(cell, species_2, positions_1) - - env2_1_0 = env.AtomicEnvironment(test_structure_1, 0, cutoffs) - - sig = random() - ls = random() - d1 = 1 - - hyps = np.array([sig, ls]) - - # check force kernel - calc1 = en.many_body_en(env1_2_0, env2_1_0, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_0, env2_1_0, hyps, cutoffs) - kern_finite_diff_00 = (calc1 - calc2) / (2 * delt) - - calc1 = en.many_body_en(env1_2_1, env2_1_0, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_1, env2_1_0, hyps, cutoffs) - kern_finite_diff_10 = (calc1 - calc2) / (2 * delt) - - calc1 = en.many_body_en(env1_2_2, env2_1_0, hyps, cutoffs) - calc2 = en.many_body_en(env1_3_2, env2_1_0, hyps, cutoffs) - kern_finite_diff_20 = (calc1 - calc2) / (2 * delt) - - kern_finite_diff = -(kern_finite_diff_00 + kern_finite_diff_10 + kern_finite_diff_20) - - kern_analytical = en.many_body_force_en(env1_1_0, env2_1_0, - d1, hyps, cutoffs) - - tol = 1e-4 - - assert (np.isclose(kern_finite_diff, kern_analytical, atol=tol)) - - -def test_many_body_grad(): - # create env 1 - cell = np.eye(3) - cutoffs = np.array([1, 1, 1]) - - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_1 = [1, 1, 1] - atom_1 = 0 - test_structure_1 = struc.Structure(cell, species_1, positions_1) - env1 = env.AtomicEnvironment(test_structure_1, atom_1, cutoffs) - - # create env 2 - positions_1 = [np.array([0., 0., 0.]), - np.array([random(), random(), random()]), - np.array([random(), random(), random()])] - - species_2 = [1, 1, 1] - atom_2 = 0 - test_structure_1 = struc.Structure(cell, species_2, positions_1) - env2 = env.AtomicEnvironment(test_structure_1, atom_2, cutoffs) - - sig = random() - ls = random() - - d1 = randint(1, 3) - d2 = randint(1, 3) - - hyps = np.array([sig, ls]) - - grad_test = en.many_body_grad(env1, env2, d1, d2, hyps, cutoffs) - - delta = 1e-8 - new_sig = sig + delta - new_ls = ls + delta - - sig_derv_brute = (en.many_body(env1, env2, d1, d2, - np.array([new_sig, ls]), - cutoffs) - - en.many_body(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - l_derv_brute = (en.many_body(env1, env2, d1, d2, - np.array([sig, new_ls]), - cutoffs) - - en.many_body(env1, env2, d1, d2, - hyps, cutoffs)) / delta - - tol = 1e-4 - assert (np.isclose(grad_test[1][0], sig_derv_brute, atol=tol)) - assert (np.isclose(grad_test[1][1], l_derv_brute, atol=tol)) diff --git a/tests/test_mc_mb_kernel.py b/tests/test_mc_mb_kernel.py index 07508780a..bea3afbdc 100644 --- a/tests/test_mc_mb_kernel.py +++ b/tests/test_mc_mb_kernel.py @@ -11,7 +11,7 @@ from .fake_gp import generate_mb_envs list_to_test = ['2', '3', '2+3', 'mb', '2+3+mb'] - +list_type = ['sc', 'mc'] def generate_hm(kernel_name): hyps = [] @@ -23,24 +23,27 @@ def generate_hm(kernel_name): @pytest.mark.parametrize('kernel_name', list_to_test) -def test_force_en(kernel_name): +@pytest.mark.parametrize('kernel_type', list_type) +def test_force_en(kernel_name, kernel_type): """Check that the analytical force/en kernel matches finite difference of energy kernel.""" cutoffs = np.ones(3)*1.2 - delta = 1e-4 + delta = 1e-5 tol = 1e-4 cell = 1e7 * np.eye(3) # set hyperparameters d1 = 1 + np.random.seed(10) env1 = generate_mb_envs(cutoffs, cell, delta, d1) env2 = generate_mb_envs(cutoffs, cell, delta, d1) hyps = generate_hm(kernel_name) - _, __, en_kernel, force_en_kernel = str_to_kernel_set(kernel_name) + _, __, en_kernel, force_en_kernel = str_to_kernel_set(kernel_name+kernel_type) + print(force_en_kernel.__name__) nterm = 0 for term in ['2', '3', 'mb']: @@ -49,7 +52,7 @@ def test_force_en(kernel_name): kern_finite_diff = 0 if ('mb' in kernel_name): - _, __, enm_kernel, ___ = str_to_kernel_set('mb') + _, __, enm_kernel, ___ = str_to_kernel_set('mb'+kernel_type) mhyps = hyps[(nterm-1)*2:] calc = 0 for i in range(len(env1[0])): @@ -60,7 +63,7 @@ def test_force_en(kernel_name): if ('2' in kernel_name): nbond = 1 - _, __, en2_kernel, ___ = str_to_kernel_set('2') + _, __, en2_kernel, ___ = str_to_kernel_set('2'+kernel_type) calc1 = en2_kernel(env1[2][0], env2[0][0], hyps[0:nbond * 2], cutoffs) calc2 = en2_kernel(env1[1][0], env2[0][0], hyps[0:nbond * 2], cutoffs) diff2b = (calc1 - calc2) / 2.0 / 2.0 / delta @@ -70,7 +73,7 @@ def test_force_en(kernel_name): nbond = 0 if ('3' in kernel_name): - _, __, en3_kernel, ___ = str_to_kernel_set('3') + _, __, en3_kernel, ___ = str_to_kernel_set('3'+kernel_type) calc1 = en3_kernel(env1[2][0], env2[0][0], hyps[nbond * 2:], cutoffs) calc2 = en3_kernel(env1[1][0], env2[0][0], hyps[nbond * 2:], cutoffs) diff3b = (calc1 - calc2) / 2.0 / 3.0 / delta @@ -85,22 +88,24 @@ def test_force_en(kernel_name): @pytest.mark.parametrize('kernel_name', list_to_test) -def test_force(kernel_name): +@pytest.mark.parametrize('kernel_type', list_type) +def test_force(kernel_name, kernel_type): """Check that the analytical force kernel matches finite difference of energy kernel.""" d1 = 1 d2 = 2 - tol = 1e-4 + tol = 1e-3 cell = 1e7 * np.eye(3) - delta = 1e-5 + delta = 1e-4 + cutoffs = np.ones(3)*1.2 + + np.random.seed(10) hyps = generate_hm(kernel_name) - cutoffs = np.ones(3)*1.2 - kernel, kg, en_kernel, fek = str_to_kernel_set(kernel_name, False) + kernel, kg, en_kernel, fek = str_to_kernel_set(kernel_name+kernel_type, False) args = (hyps, cutoffs) - np.random.seed(10) env1 = generate_mb_envs(cutoffs, cell, delta, d1) env2 = generate_mb_envs(cutoffs, cell, delta, d2) @@ -125,12 +130,12 @@ def test_force(kernel_name): kern_analytical = kernel(env1[0][0], env2[0][0], d1, d2, *args) - tol = 1e-4 assert(isclose(kern_finite_diff, kern_analytical, rtol=tol)) @pytest.mark.parametrize('kernel_name', list_to_test) -def test_hyps_grad(kernel_name): +@pytest.mark.parametrize('kernel_type', list_type) +def test_hyps_grad(kernel_name, kernel_type): delta = 1e-8 cutoffs = np.ones(3)*1.2 @@ -141,10 +146,11 @@ def test_hyps_grad(kernel_name): print("hyps", hyps) cell = 1e7 * np.eye(3) + np.random.seed(10) env1 = generate_mb_envs(cutoffs, cell, 0, d1)[0][0] env2 = generate_mb_envs(cutoffs, cell, 0, d2)[0][0] - kernel, kernel_grad, _, _ = str_to_kernel_set(kernel_name, False) + kernel, kernel_grad, _, _ = str_to_kernel_set(kernel_name+kernel_type, False) grad_test = kernel_grad(env1, env2, d1, d2, hyps, cutoffs) From f6ce594edcf99f2d234544927e6b5bafcb24fe65 Mon Sep 17 00:00:00 2001 From: nw13slx Date: Thu, 23 Apr 2020 15:22:07 -0400 Subject: [PATCH 4/6] rename test --- tests/{test_mc_mb_kernel.py => test_kernel.py} | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) rename tests/{test_mc_mb_kernel.py => test_kernel.py} (99%) diff --git a/tests/test_mc_mb_kernel.py b/tests/test_kernel.py similarity index 99% rename from tests/test_mc_mb_kernel.py rename to tests/test_kernel.py index bea3afbdc..1f16fc51f 100644 --- a/tests/test_mc_mb_kernel.py +++ b/tests/test_kernel.py @@ -137,16 +137,15 @@ def test_force(kernel_name, kernel_type): @pytest.mark.parametrize('kernel_type', list_type) def test_hyps_grad(kernel_name, kernel_type): - delta = 1e-8 - cutoffs = np.ones(3)*1.2 - - hyps = generate_hm(kernel_name) d1 = randint(1, 3) d2 = randint(1, 3) - print("hyps", hyps) + tol = 1e-4 cell = 1e7 * np.eye(3) + delta = 1e-8 + cutoffs = np.ones(3)*1.2 np.random.seed(10) + hyps = generate_hm(kernel_name) env1 = generate_mb_envs(cutoffs, cell, 0, d1)[0][0] env2 = generate_mb_envs(cutoffs, cell, 0, d2)[0][0] @@ -155,7 +154,6 @@ def test_hyps_grad(kernel_name, kernel_type): grad_test = kernel_grad(env1, env2, d1, d2, hyps, cutoffs) - tol = 1e-4 original = kernel(env1, env2, d1, d2, hyps, cutoffs) for i in range(len(hyps)-1): From cdee021f7936d65289ba469b8755b2faaceee5c5 Mon Sep 17 00:00:00 2001 From: nw13slx Date: Sat, 25 Apr 2020 12:27:05 -0400 Subject: [PATCH 5/6] fix another wrong argument --- flare/kernels/mc_simple.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/flare/kernels/mc_simple.py b/flare/kernels/mc_simple.py index e7e3391ec..c3cb70ba5 100644 --- a/flare/kernels/mc_simple.py +++ b/flare/kernels/mc_simple.py @@ -410,7 +410,7 @@ def two_plus_three_plus_many_body_mc_en(env1: AtomicEnvironment, env2: AtomicEnv env1.triplet_counts, env2.triplet_counts, sig3, ls3, r_cut_3, cutoff_func) - many_term = many_body_mc_en_jit(env1.bond_array_2, env2.bond_array_2, env1.ctype, + many_term = many_body_mc_en_jit(env1.bond_array_mb, env2.bond_array_mb, env1.ctype, env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, env1.species, env2.species, From 7bdeb6d99e8e08c765118f919da8d7b829e0fd79 Mon Sep 17 00:00:00 2001 From: Shelby Xie Date: Thu, 30 Apr 2020 12:27:16 -0400 Subject: [PATCH 6/6] change format --- flare/kernels/mc_sephyps.py | 37 ++++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/flare/kernels/mc_sephyps.py b/flare/kernels/mc_sephyps.py index ac7c86e46..2ea3cb1bc 100644 --- a/flare/kernels/mc_sephyps.py +++ b/flare/kernels/mc_sephyps.py @@ -94,10 +94,13 @@ def two_three_many_body_mc(env1, env2, d1, d2, cutoffs, d1, d2, sig3, ls3, r_cut_3, cutoff_func, nspec, spec_mask, triplet_mask) - many_term = many_body_mc_jit(env1.bond_array_mb, env2.bond_array_mb, env1.neigh_dists_mb, - env2.neigh_dists_mb, env1.num_neighs_mb, env2.num_neighs_mb, - env1.ctype, env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, - env1.etype_mb, env2.etype_mb, env1.species, env2.species, + many_term = many_body_mc_jit(env1.bond_array_mb, env2.bond_array_mb, + env1.neigh_dists_mb, env2.neigh_dists_mb, + env1.num_neighs_mb, env2.num_neighs_mb, + env1.ctype, env2.ctype, + env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, + env1.etype_mb, env2.etype_mb, + env1.species, env2.species, d1, d2, sigm, lsm, r_cut_m, cutoff_func) return two_term + three_term + many_term @@ -133,11 +136,12 @@ def two_three_many_body_mc_grad(env1, env2, d1, d2, cutoffs, kern_many, gradm = many_body_mc_grad_jit(env1.bond_array_mb, env2.bond_array_mb, env1.neigh_dists_mb, env2.neigh_dists_mb, - env1.num_neighs_mb, env2.num_neighs_mb, env1.ctype, - env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, + env1.num_neighs_mb, env2.num_neighs_mb, + env1.ctype, env2.ctype, + env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, env1.etype_mb, env2.etype_mb, - env1.species, env2.species, d1, d2, sigm, - lsm, r_cut_m, cutoff_func) + env1.species, env2.species, + d1, d2, sigm, lsm, r_cut_m, cutoff_func) g = np.hstack([grad2, grad3, gradm]) @@ -173,11 +177,13 @@ def two_three_many_mc_force_en(env1, env2, d1, cutoffs, triplet_mask) / 3 many_term = many_body_mc_force_en_jit(env1.bond_array_mb, env2.bond_array_mb, - env1.neigh_dists_mb, env1.num_neighs_mb, - env1.ctype, env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, + env1.neigh_dists_mb, + env1.num_neighs_mb, + env1.ctype, env2.ctype, + env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, env1.etype_mb, - env1.species, env2.species, d1, sigm, lsm, r_cut_m, - cutoff_func) + env1.species, env2.species, + d1, sigm, lsm, r_cut_m, cutoff_func) return two_term + three_term + many_term @@ -208,9 +214,10 @@ def two_three_many_mc_en(env1, env2, cutoffs, nspec, spec_mask, triplet_mask) - many_term = many_body_mc_en_jit(env1.bond_array_2, env2.bond_array_2, env1.ctype, - env2.ctype, env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, env1.species, - env2.species, + many_term = many_body_mc_en_jit(env1.bond_array_2, env2.bond_array_2, + env1.ctype, env2.ctype, + env1.bond_array_mb_etypes, env2.bond_array_mb_etypes, + env1.species, env2.species, sigm, lsm, r_cut_m, cutoff_func)