From 2e7a67146399bae041cb64705552bffc6e1c960e Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Fri, 27 Sep 2019 16:15:57 +0200 Subject: [PATCH 1/6] thole damping, screening parameters, and damped T tensors --- cppe/core/math.cc | 75 ++++++++++++++++++++++++++++++++++++++++++---- cppe/core/math.hh | 11 +++++-- tests/test_math.py | 18 +++++++++++ 3 files changed, 96 insertions(+), 8 deletions(-) diff --git a/cppe/core/math.cc b/cppe/core/math.cc index 0daa1c9b..1cbe2be7 100644 --- a/cppe/core/math.cc +++ b/cppe/core/math.cc @@ -19,7 +19,8 @@ Eigen::Vector3d smat_vec(const Eigen::VectorXd& mat, const Eigen::Vector3d& vec, // TODO: add option for damping Eigen::VectorXd Tk_tensor(int k, const Eigen::Vector3d& Rij, - std::vector& Tk_coeffs) { + std::vector& Tk_coeffs, double damping_factor, + double alpha_i, double alpha_j) { int x, y, z; Eigen::VectorXd Tk(multipole_components(k)); int idx = 0; @@ -27,7 +28,7 @@ Eigen::VectorXd Tk_tensor(int k, const Eigen::Vector3d& Rij, for (y = k; y > -1; y--) { for (z = k; z > -1; z--) { if (x + y + z != k) continue; - Tk[idx] = T(Rij, x, y, z, Tk_coeffs); + Tk[idx] = T(Rij, x, y, z, Tk_coeffs, damping_factor, alpha_i, alpha_j); idx++; } } @@ -97,16 +98,57 @@ int xyz2idx(int x, int y, int z) { } double T(const Eigen::Vector3d& Rij, int x, int y, int z, - std::vector& Cijn) { + std::vector& Cijn, double damping_factor, double alpha_i, + double alpha_j) { double t = 0.0; double R = Rij.norm(); double Cx, Cy, Cz; + int k = x + y + z; + + // Thole Damping + std::vector scr_facs; + if (damping_factor != 0.0) { + if (alpha_i <= 0.0 || alpha_j <= 0.0) { + throw std::runtime_error( + "Thole damping only valid for non-zero" + " isotropic polarizabilities."); + } + if (k > 3) { + throw std::runtime_error("Thole damping only implemented up to third order."); + } + // Molecular Simulation, 32:6, 471-484, DOI: 10.1080/08927020600631270 + // v = factor * u, with u = R / (alpha_i * alpha_j)**(1/6) + double v = damping_factor * R / std::pow(alpha_i * alpha_j, 1.0 / 6.0); + scr_facs = thole_screening_factors(v, k); + } + for (size_t l = 0; l <= x; l++) { Cx = Cijn[0](x, l) * pow((Rij(0) / R), l); for (size_t m = 0; m <= y; m++) { Cy = Cx * Cijn[l + x](y, m) * pow((Rij(1) / R), m); for (size_t n = 0; n <= z; n++) { - Cz = Cy * Cijn[l + x + m + y](z, n) * pow((Rij(2) / R), n); + Cz = Cy * Cijn[l + x + m + y](z, n) * pow((Rij(2) / R), n); + int kk = l + m + n; + // Thole damping + if (scr_facs.size()) { + if (kk == 0) { + if (k == 0) { + Cz *= scr_facs[0]; + } else if (k == 2) { + Cz *= scr_facs[1]; + } + } else if (kk == 1) { + if (k == 1) { + Cz *= scr_facs[1]; + } else if (k == 3) { + Cz *= scr_facs[3]; + } + } else if (kk == 2) { + if (k == 2) Cz *= scr_facs[2]; + } else if (kk == 3) { + if (k == 3) Cz *= scr_facs[3]; + } + } t += Cz; } } @@ -115,6 +157,29 @@ double T(const Eigen::Vector3d& Rij, int x, int y, int z, return t; } +std::vector thole_screening_factors(double v, int k) { + std::vector ret; + // screening factors for potential, field, field gradient, field Hessian + double f_v, f_E, f_T, f_D = 1.0; + if (k >= 0) { + f_v = 1.0 - (0.5 * v + 1.0) * std::exp(-v); + ret.push_back(f_v); + } + if (k >= 1) { + f_E = f_v - (0.5 * std::pow(v, 2) + 0.5 * v) * std::exp(-v); + ret.push_back(f_E); + } + if (k >= 2) { + f_T = f_E - 1.0 / 6.0 * std::pow(v, 3) * std::exp(-v); + ret.push_back(f_T); + } + if (k >= 3) { + f_D = f_T - 1.0 / 30.0 * std::pow(v, 4) * std::exp(-v); + ret.push_back(f_D); + } + return ret; +} + std::vector Tk_coefficients(int max_order) { int maxi = 2 * max_order + 3; std::vector Cijn; @@ -219,4 +284,4 @@ std::vector prefactors_nuclei(unsigned k) { int multipole_components(int k) { return (k + 1) * (k + 2) / 2; } -} // namespace libcppe \ No newline at end of file +} // namespace libcppe diff --git a/cppe/core/math.hh b/cppe/core/math.hh index 748d4869..defd4706 100644 --- a/cppe/core/math.hh +++ b/cppe/core/math.hh @@ -13,10 +13,15 @@ Eigen::VectorXd multipole_derivative(int k, int l, const Eigen::Vector3d& Rji, std::vector& Tk_coeffs); double T(const Eigen::Vector3d& Rij, int x, int y, int z, - std::vector& Cijn); + std::vector& Cijn, double damping_factor = 0.0, + double alpha_i = 0.0, double alpha_j = 0.0); Eigen::VectorXd Tk_tensor(int k, const Eigen::Vector3d& Rij, - std::vector& Tk_coeffs); + std::vector& Tk_coeffs, + double damping_factor = 0.0, double alpha_i = 0.0, + double alpha_j = 0.0); + +std::vector thole_screening_factors(double v, int k); std::vector Tk_coefficients(int max_order); @@ -37,4 +42,4 @@ std::vector prefactors_nuclei(unsigned k); int multipole_components(int k); -} // namespace libcppe \ No newline at end of file +} // namespace libcppe diff --git a/tests/test_math.py b/tests/test_math.py index f78a51fc..88b5a56a 100644 --- a/tests/test_math.py +++ b/tests/test_math.py @@ -71,3 +71,21 @@ def test_T_tensors(self): sym_indices = symmetry.get_symm_indices(k) np.testing.assert_almost_equal(actual, ref.take(sym_indices), decimal=10) + + def test_T_tensors_damped(self): + # tests the T tensors against auto-generated Python code + ref_T = tensors.T + + for k in range(4): + R = np.random.random(3) + ref = ref_T[k](*R) + + coeffs = Tk_coefficients(k) + # damped tensor: damping_factor, alpha_i, alpha_j + actual = Tk_tensor(k, R, coeffs, 1.0, 1.0, 1.0) + + # gets the indices of non-redundant components + # e.g., takes only xy from (xy, yx) and so on ... + sym_indices = symmetry.get_symm_indices(k) +# np.testing.assert_almost_equal(actual, +# ref.take(sym_indices), decimal=10) From 51fae746cbdd845663f1483efb895f9ba1b45c77 Mon Sep 17 00:00:00 2001 From: Peter Reinholdt Date: Fri, 27 Sep 2019 16:28:33 +0200 Subject: [PATCH 2/6] Add autogenerated code --- tests/tensors/tensors.py | 2296 +++++++++++++++++++++++++++++++++++++- 1 file changed, 2295 insertions(+), 1 deletion(-) diff --git a/tests/tensors/tensors.py b/tests/tensors/tensors.py index fd10f554..48952f16 100644 --- a/tests/tensors/tensors.py +++ b/tests/tensors/tensors.py @@ -1,5 +1,5 @@ from numba import jit -from numpy import sqrt +from numpy import sqrt, exp import numpy as np @jit(nopython=True, cache=True) def T0(x,y,z): @@ -1477,3 +1477,2297 @@ def T5(x,y,z): arr[(2, 2, 2, 2, 2)] = Tzzzzz(x,y,z) return arr T = [T0, T1, T2, T3, T4, T5, ] + +def T_damp_thole_0(x, y, z, a): + shape = (3,)*0 + result = np.zeros(shape, dtype=np.float64) + x0 = sqrt(x ** 2 + y ** 2 + z ** 2) + x1 = a * x0 + result[...,] = -((x1 + 2) * exp(-x1) / 2 - 1) / x0 + return result + + +def T_damp_thole_1(x, y, z, a): + shape = (3,)*1 + result = np.zeros(shape, dtype=np.float64) + x0 = x ** 2 + y ** 2 + z ** 2 + x1 = a * sqrt(x0) + x2 = -x1 + x3 = exp(x2) + x4 = -a * x3 * (x2 - 1) / (2 * x0) + (x3 * (x1 + 2) - 2) / (2 * x0 ** (3 / 2)) + result[0] = x * x4 + result[1] = x4 * y + result[2] = x4 * z + return result + + +def T_damp_thole_2(x, y, z, a): + shape = (3,)*2 + result = np.zeros(shape, dtype=np.float64) + x0 = x ** 2 + x1 = y ** 2 + x2 = z ** 2 + x3 = x0 + x1 + x2 + x4 = sqrt(x3) + x5 = a * x4 + x6 = exp(-x5) + x7 = a * x6 + x8 = x7 * (x5 + 1) / x3 ** 2 + x9 = 1 / x3 + x10 = 3 * x9 + x11 = x3 ** (-3 / 2) + x12 = x5 + 2 + x13 = x12 * x6 - 2 + x14 = x11 * x13 / 2 + x15 = x0 * x11 + x16 = a * x9 + x17 = 2 * x16 + x18 = x12 * x16 + x19 = 1 / x4 + x20 = -x12 * x19 + x19 + x21 = x19 * x7 / 2 + x22 = x11 * x12 + x23 = 3 * x13 / (2 * x3 ** (5 / 2)) + x21 * (-x11 - x17 + x18 + x22) + x8 + x24 = x * x23 + x25 = -x24 * y + x26 = -x24 * z + x27 = -x23 * y * z + result[0, 0] = ( + -x0 * x8 + - x14 * (x0 * x10 - 1) + - x21 * (-x0 * x17 + x0 * x18 + x12 * x15 - x15 + x20) + ) + result[0, 1] = x25 + result[0, 2] = x26 + result[1, 0] = x25 + result[1, 1] = ( + -x1 * x8 + - x14 * (x1 * x10 - 1) + - x21 * (-x1 * x11 - x1 * x17 + x1 * x18 + x1 * x22 + x20) + ) + result[1, 2] = x27 + result[2, 0] = x26 + result[2, 1] = x27 + result[2, 2] = ( + -x14 * (x10 * x2 - 1) + - x2 * x8 + - x21 * (-x11 * x2 - x17 * x2 + x18 * x2 + x2 * x22 + x20) + ) + return result + + +def T_damp_thole_3(x, y, z, a): + shape = (3,)*3 + result = np.zeros(shape, dtype=np.float64) + x0 = x ** 2 + x1 = y ** 2 + x2 = z ** 2 + x3 = x0 + x1 + x2 + x4 = 1 / x3 + x5 = x0 * x4 + x6 = 3 * x5 - 1 + x7 = sqrt(x3) + x8 = a * x7 + x9 = exp(-x8) + x10 = x8 + 1 + x11 = a / x3 ** 2 + x12 = x10 * x11 + x13 = 3 * x12 * x9 + x14 = x8 + 2 + x15 = x14 * x9 - 2 + x16 = x3 ** (-5 / 2) + x17 = 3 * x16 + x18 = x15 * x17 + x19 = x3 ** (-3 / 2) + x20 = 3 * x19 + x21 = x0 * x19 + x22 = x14 * x21 + x23 = a * x4 + x24 = 2 * x23 + x25 = x14 * x23 + x26 = 1 / x7 + x27 = -x14 * x26 + x26 + x28 = -x0 * x24 + x0 * x25 - x21 + x22 + x27 + x29 = a * x9 + x30 = x28 * x29 + x31 = x0 * x17 + x32 = 6 * x11 + x33 = a ** 2 + x34 = x20 * x33 + x35 = x14 * x31 + x36 = x11 * x14 + x37 = 3 * x36 + x38 = -x14 * x20 + x20 + 6 * x23 - 3 * x25 + x39 = x26 * x29 + x40 = x / 2 + x41 = x10 * x29 / x3 ** 3 + x42 = 3 * x41 + x43 = x15 / x3 ** (7 / 2) + x44 = 3 * x43 + x45 = x9 / 2 + x46 = x12 * x45 + x47 = 3 * x15 * x16 / 2 + x48 = x14 * x19 + x49 = -x19 + x48 + x50 = x29 * (-x24 + x25 + x49) + x51 = x19 / 2 + x52 = x33 * x4 * x45 + x53 = 5 * x11 + x54 = 2 * x36 + x55 = -x23 + x49 + x56 = x39 / 2 + x57 = ( + x0 * x42 + + x0 * x44 + + x21 * x50 + + x28 * x52 + + x30 * x51 + + x46 * x6 + + x47 * x6 + - x56 * (x0 * x53 - x0 * x54 + x21 * x33 + x31 - x35 + x55) + ) + x58 = x57 * y + x59 = x57 * z + x60 = 15 * x43 / 2 + x61 = x19 + x24 - x25 - x48 + x62 = x1 * x17 + x63 = x14 * x62 + x64 = x1 * x48 + x65 = -x1 * x32 - x1 * x34 + x1 * x37 + x33 * x64 - x62 + x63 + x66 = 9 * x41 / 2 + x67 = -x46 - x47 + x68 = x1 * x19 + x69 = -x1 * x24 + x1 * x25 + x27 + x64 - x68 + x70 = x29 * x51 + x71 = x50 * x68 + x69 * x70 + x72 = x * (x1 * x60 + x1 * x66 + x56 * (x61 + x65) + x67 + x71) + x73 = x14 * x17 + x74 = x33 * x48 + x75 = ( + x40 + * y + * z + * (x20 * x50 + x39 * (-x17 - x32 - x34 + x37 + x73 + x74) + 9 * x41 + 15 * x43) + ) + x76 = x19 * x2 + x77 = -x2 * x24 + x2 * x25 + x2 * x48 + x27 - x76 + x78 = -x17 * x2 - x2 * x32 - x2 * x34 + x2 * x37 + x2 * x73 + x2 * x74 + x79 = x2 * x60 + x2 * x66 + x50 * x76 + x56 * (x61 + x78) + x67 + x70 * x77 + x80 = x * x79 + x81 = x1 * x4 + x82 = 3 * x81 - 1 + x83 = x20 * x29 + x84 = z * ( + x1 * x42 + + x1 * x44 + + x46 * x82 + + x47 * x82 + + x52 * x69 + - x56 * (x1 * x53 - x1 * x54 + x33 * x68 + x55 + x62 - x63) + + x71 + ) + x85 = x79 * y + x86 = x2 * x4 + result[0, 0, 0] = x40 * ( + x13 * x6 + + x18 * (5 * x5 - 3) + + x20 * x30 + + x39 * (-x0 * x32 - x0 * x34 + x0 * x37 + x22 * x33 - x31 + x35 + x38) + ) + result[0, 0, 1] = x58 + result[0, 0, 2] = x59 + result[0, 1, 0] = x58 + result[0, 1, 1] = x72 + result[0, 1, 2] = x75 + result[0, 2, 0] = x59 + result[0, 2, 1] = x75 + result[0, 2, 2] = x80 + result[1, 0, 0] = x58 + result[1, 0, 1] = x72 + result[1, 0, 2] = x75 + result[1, 1, 0] = x72 + result[1, 1, 1] = ( + y * (x13 * x82 + x18 * (5 * x81 - 3) + x39 * (x38 + x65) + x69 * x83) / 2 + ) + result[1, 1, 2] = x84 + result[1, 2, 0] = x75 + result[1, 2, 1] = x84 + result[1, 2, 2] = x85 + result[2, 0, 0] = x59 + result[2, 0, 1] = x75 + result[2, 0, 2] = x80 + result[2, 1, 0] = x75 + result[2, 1, 1] = x84 + result[2, 1, 2] = x85 + result[2, 2, 0] = x80 + result[2, 2, 1] = x85 + result[2, 2, 2] = ( + z + * (x13 * (3 * x86 - 1) + x18 * (5 * x86 - 3) + x39 * (x38 + x78) + x77 * x83) + / 2 + ) + return result + + +def T_damp_thole_4(x, y, z, a): + shape = (3,)*4 + result = np.zeros(shape, dtype=np.float64) + x0 = x ** 2 + x1 = y ** 2 + x2 = z ** 2 + x3 = x0 + x1 + x2 + x4 = 1 / x3 + x5 = x0 * x4 + x6 = 5 * x5 - 3 + x7 = a / x3 ** 3 + x8 = x0 * x7 + x9 = sqrt(x3) + x10 = a * x9 + x11 = exp(-x10) + x12 = x10 + 1 + x13 = x11 * x12 + x14 = 6 * x13 + x15 = x ** 4 + x16 = x3 ** (-2) + x17 = 35 * x16 + x18 = x10 + 2 + x19 = x11 * x18 - 2 + x20 = x3 ** (-5 / 2) + x21 = 3 * x20 / 2 + x22 = x19 * x21 + x23 = 3 * x5 - 1 + x24 = x3 ** (-3 / 2) + x25 = 3 * x24 + x26 = x0 * x24 + x27 = x18 * x26 + x28 = a * x4 + x29 = 2 * x28 + x30 = x18 * x28 + x31 = 1 / x9 + x32 = -x18 * x31 + x31 + x33 = -x0 * x29 + x0 * x30 - x26 + x27 + x32 + x34 = x11 * x33 + x35 = a * x34 + x36 = a * x16 + x37 = 6 * x36 + x38 = a ** 2 + x39 = x25 * x38 + x40 = x27 * x38 + x41 = 3 * x36 + x42 = x18 * x41 + x43 = 3 * x20 + x44 = x0 * x43 + x45 = x18 * x44 + x46 = -x44 + x45 + x47 = x18 * x25 + x48 = 6 * x28 + x49 = 3 * x30 + x50 = x25 - x47 + x48 - x49 + x51 = -x0 * x37 - x0 * x39 + x0 * x42 + x40 + x46 + x50 + x52 = a * x11 + x53 = x26 * x52 + x54 = x3 ** (-7 / 2) + x55 = 15 * x54 + x56 = x15 * x55 + x57 = 18 * x20 + x58 = x0 * x57 + x59 = 30 * x7 + x60 = x38 * x57 + x61 = a ** 3 + x62 = x16 * x61 + x63 = x15 * x62 + x64 = x26 * x38 + x65 = x0 * x36 + x66 = 18 * x36 + x67 = x18 * x66 + x68 = x15 * x18 + x69 = x20 * x38 + x70 = 6 * x69 + x71 = 15 * x7 + x72 = -x25 + x47 - x48 + x49 + x73 = x11 / 2 + x74 = a * x31 * x73 + x75 = 9 * x12 + x76 = x3 ** (-4) + x77 = x0 * x52 + x78 = x76 * x77 + x79 = x19 / x3 ** (9 / 2) + x80 = x0 * x79 + x81 = x11 * x7 + x82 = x12 * x81 + x83 = 9 * x82 / 2 + x84 = 3 * x82 / 2 + x85 = x19 * x54 + x86 = 15 * x85 / 2 + x87 = x18 * x24 + x88 = -x24 + x87 + x89 = -x29 + x30 + x88 + x90 = a * x24 + x91 = x11 * x90 + x92 = 3 * x91 / 2 + x93 = x89 * x92 + x94 = 9 * x20 / 2 + x95 = 3 * x16 * x38 / 2 + x96 = x34 * x95 + x97 = 5 * x65 + x98 = 2 * x18 * x65 + x99 = -x28 + x88 + x100 = x11 * (x44 - x45 + x64 + x97 - x98 + x99) + x101 = x100 * x90 + x102 = x91 / 2 + x103 = x38 * x4 + x104 = x103 * x73 + x105 = 12 * x0 + x106 = x18 * x7 + x107 = x0 * x55 + x108 = x107 * x18 + x109 = x107 - x108 + x110 = -x39 + x111 = 9 * x20 + x112 = x111 * x18 + x113 = x110 - x111 + x112 + x18 * x37 - 15 * x36 + x114 = x * ( + -3 * x101 / 2 + + x102 * x51 + + x104 * x51 + + x23 * x83 + + x23 * x93 + + x35 * x94 + + x6 * x84 + + x6 * x86 + - x74 + * (x0 * x62 - x105 * x106 + x105 * x69 + x109 + x113 - x38 * x45 + 27 * x8) + + x75 * x78 + + 15 * x80 + + x96 + ) + x115 = -x114 * y + x116 = -x114 * z + x117 = x1 * x41 + x118 = 5 * x69 + x119 = x0 * x118 + x120 = 23 * x8 + x121 = 8 * x18 + x122 = x121 * x8 + x123 = x24 - x87 + x124 = x1 * x43 + x125 = x124 * x18 + x126 = -x124 + x125 + x127 = x123 + x126 + x28 + x128 = x46 - x64 - x97 + x98 + x129 = 30 * x80 + x130 = x23 * x86 + x131 = x1 * x24 + x132 = a * x100 + x133 = x1 * x4 + x134 = x100 * x38 + x135 = x123 + x29 - x30 + x136 = x1 * x87 + x137 = x136 * x38 + x138 = -x1 * x37 - x1 * x39 + x117 * x18 + x126 + x137 + x139 = x135 + x138 + x140 = x21 * x35 + x141 = -x1 * x29 + x1 * x30 - x131 + x136 + x32 + x142 = x102 * x141 + x143 = x61 * x73 + x144 = x143 * x33 + x145 = 21 * x12 + x146 = x145 * x78 + x147 = x20 * x89 + x148 = 6 * x147 + x149 = x148 * x77 + x150 = x1 * x7 + x151 = 3 * x13 + x152 = x150 * x151 + x153 = 3 * x85 + x154 = x0 * x153 + x102 * x33 + x104 * x33 + x151 * x8 + x22 * x23 + x155 = ( + -x1 * x129 + - x1 * x130 + - x1 * x140 + - x1 * x146 + - x1 * x149 + - x1 * x96 + + x131 * x132 + - x131 * x144 + + x133 * x134 + - x139 * x53 + - x142 * x23 + - x152 * x23 + + x154 + + x74 + * ( + x1 * x107 + - x1 * x108 + + x1 * x119 + + x1 * x120 + - x1 * x122 + - x117 + + x127 + + x128 + ) + ) + x156 = 3 * x82 + x157 = x156 * x23 + x158 = x102 * x89 + x159 = x38 * x87 + x160 = x18 * x43 + x161 = x160 - x43 + x162 = x110 + x159 + x161 - x37 + x42 + x163 = y * z + x164 = -x163 * ( + -x100 * x103 + - x101 + + x129 + + x130 + + x140 + + x144 * x24 + + x146 + + x149 + + x157 + + x158 * x23 + + x162 * x53 + - x74 * (x109 + x119 + x120 - x122 + x161 - x41) + + x96 + ) + x165 = x160 * x2 - x2 * x43 + x166 = x165 - x2 * x41 + x167 = x2 * x24 + x168 = x2 * x4 + x169 = x159 * x2 + x170 = x165 + x169 - x2 * x37 - x2 * x39 + x2 * x42 + x171 = x135 + x170 + x172 = -x167 - x2 * x29 + x2 * x30 + x2 * x87 + x32 + x173 = x102 * x172 + x174 = ( + -x129 * x2 + - x130 * x2 + + x132 * x167 + + x134 * x168 + - x140 * x2 + - x144 * x167 + - x146 * x2 + - x149 * x2 + + x154 + - x157 * x2 + - x171 * x53 + - x173 * x23 + - x2 * x96 + + x74 + * ( + x107 * x2 + - x108 * x2 + + x119 * x2 + + x120 * x2 + - x122 * x2 + + x123 + + x128 + + x166 + + x28 + ) + ) + x175 = x18 * x36 + x176 = x111 - x112 - 9 * x175 + 9 * x24 * x38 - x38 * x47 + x66 + x177 = x1 * x55 + x178 = x1 * x57 + x179 = x1 * x62 + x180 = x177 * x18 + x181 = x1 * x69 + x182 = 15 * x1 + x183 = ( + -x1 * x59 + + x106 * x182 + - x177 + - x178 * x38 + + x179 * x18 + - 4 * x179 + + 6 * x18 * x181 + + x180 + ) + x184 = x139 * x91 + x185 = x1 * x52 + x186 = x89 * x94 + x187 = 105 * x79 / 2 + x188 = x185 * x76 + x189 = 30 * x12 + x190 = x1 * x187 + x188 * x189 + x191 = x138 + x50 + x192 = x141 * x52 + x193 = x102 * x191 + x192 * x94 + x194 = -x75 * x81 - 45 * x85 / 2 - x93 + x195 = x * y + x196 = -x195 * ( + 3 * x184 / 2 + x185 * x186 + x190 + x193 + x194 + x74 * (x176 + x183) + ) + x197 = -x159 - x160 + x37 + x39 - x42 + x43 + x198 = x192 * x21 + x199 = x162 * x52 + x200 = 15 * x147 / 2 + x201 = -x156 - x158 - x86 + x202 = x * z + x203 = -x202 * ( + x131 * x199 + x184 + x185 * x200 + x190 + x198 + x201 + x74 * (x183 + x197) + ) + x204 = x171 * x91 + x205 = x2 * x55 + x206 = x2 * x57 + x207 = x2 * x62 + x208 = x18 * x2 + x209 = ( + x18 * x205 + + x18 * x207 + - x2 * x59 + - x205 + - x206 * x38 + - 4 * x207 + + x208 * x70 + + x208 * x71 + ) + x210 = x172 * x52 + x211 = x2 * x52 + x212 = x187 * x2 + x189 * x211 * x76 + x213 = -x195 * ( + x167 * x199 + + x200 * x211 + + x201 + + x204 + + x21 * x210 + + x212 + + x74 * (x197 + x209) + ) + x214 = x170 + x50 + x215 = ( + x102 * x214 + + x186 * x211 + + x194 + + 3 * x204 / 2 + + x210 * x94 + + x212 + + x74 * (x176 + x209) + ) + x216 = -x202 * x215 + x217 = 5 * x133 - 3 + x218 = y ** 4 + x219 = 3 * x133 - 1 + x220 = x131 * x52 + x221 = x218 * x55 + x222 = x218 * x62 + x223 = x131 * x38 + x224 = x1 * x36 + x225 = x18 * x218 + x226 = 5 * x224 + x227 = 2 * x1 * x175 + x228 = x124 - x125 + x223 + x226 - x227 + x99 + x229 = x11 * x141 * x95 + x230 = -x163 * ( + x104 * x191 + + x182 * x79 + + x188 * x75 + + x193 + + x217 * x84 + + x217 * x86 + + x219 * x83 + + x219 * x93 + - x228 * x92 + + x229 + - x74 + * ( + x113 + - x125 * x38 + - 12 * x150 * x18 + + 27 * x150 + + x177 + + x179 + - x180 + + 12 * x181 + ) + ) + x231 = x1 * x2 + x232 = x231 * x52 + x233 = x2 * x219 + x234 = x167 * x52 + x235 = x150 * x2 + x236 = ( + x1 * x153 + + x104 * x141 + + x11 * x168 * x228 * x38 + - x141 * x143 * x167 + + x142 + - x145 * x232 * x76 + - x148 * x232 + + x152 + - x156 * x233 + - x171 * x220 + - x173 * x219 + - x198 * x2 + - x2 * x229 + + x219 * x22 + + x228 * x234 + - 30 * x231 * x79 + - x233 * x86 + + x74 + * ( + x118 * x231 + - x121 * x235 + + x127 + + x166 + + x177 * x2 + - x180 * x2 + - x223 + - x226 + + x227 + + 23 * x235 + ) + ) + x237 = -x163 * x215 + x238 = z ** 4 + x239 = x238 * x55 + x240 = x238 * x62 + x241 = x18 * x238 + result[0, 0, 0, 0] = ( + -x14 * x6 * x8 + - x22 * (x15 * x17 - 30 * x5 + 3) + - x23 * x25 * x35 + - 2 * x51 * x53 + - x74 + * ( + -x0 * x67 + - x15 * x59 + - x15 * x60 + + x18 * x56 + - x18 * x58 + + x18 * x63 + - 6 * x40 + - x56 + + x58 + - 4 * x63 + + 18 * x64 + + 36 * x65 + + x68 * x70 + + x68 * x71 + + x72 + ) + ) + result[0, 0, 0, 1] = x115 + result[0, 0, 0, 2] = x116 + result[0, 0, 1, 0] = x115 + result[0, 0, 1, 1] = x155 + result[0, 0, 1, 2] = x164 + result[0, 0, 2, 0] = x116 + result[0, 0, 2, 1] = x164 + result[0, 0, 2, 2] = x174 + result[0, 1, 0, 0] = x115 + result[0, 1, 0, 1] = x155 + result[0, 1, 0, 2] = x164 + result[0, 1, 1, 0] = x155 + result[0, 1, 1, 1] = x196 + result[0, 1, 1, 2] = x203 + result[0, 1, 2, 0] = x164 + result[0, 1, 2, 1] = x203 + result[0, 1, 2, 2] = x213 + result[0, 2, 0, 0] = x116 + result[0, 2, 0, 1] = x164 + result[0, 2, 0, 2] = x174 + result[0, 2, 1, 0] = x164 + result[0, 2, 1, 1] = x203 + result[0, 2, 1, 2] = x213 + result[0, 2, 2, 0] = x174 + result[0, 2, 2, 1] = x213 + result[0, 2, 2, 2] = x216 + result[1, 0, 0, 0] = x115 + result[1, 0, 0, 1] = x155 + result[1, 0, 0, 2] = x164 + result[1, 0, 1, 0] = x155 + result[1, 0, 1, 1] = x196 + result[1, 0, 1, 2] = x203 + result[1, 0, 2, 0] = x164 + result[1, 0, 2, 1] = x203 + result[1, 0, 2, 2] = x213 + result[1, 1, 0, 0] = x155 + result[1, 1, 0, 1] = x196 + result[1, 1, 0, 2] = x203 + result[1, 1, 1, 0] = x196 + result[1, 1, 1, 1] = ( + -x14 * x150 * x217 + - 2 * x191 * x220 + - x192 * x219 * x25 + - x22 * (-30 * x133 + x17 * x218 + 3) + - x74 + * ( + -x1 * x67 + - 6 * x137 + - x178 * x18 + + x178 + + x18 * x221 + + x18 * x222 + - x218 * x59 + - x218 * x60 + - x221 + - 4 * x222 + + 18 * x223 + + 36 * x224 + + x225 * x70 + + x225 * x71 + + x72 + ) + ) + result[1, 1, 1, 2] = x230 + result[1, 1, 2, 0] = x203 + result[1, 1, 2, 1] = x230 + result[1, 1, 2, 2] = x236 + result[1, 2, 0, 0] = x164 + result[1, 2, 0, 1] = x203 + result[1, 2, 0, 2] = x213 + result[1, 2, 1, 0] = x203 + result[1, 2, 1, 1] = x230 + result[1, 2, 1, 2] = x236 + result[1, 2, 2, 0] = x213 + result[1, 2, 2, 1] = x236 + result[1, 2, 2, 2] = x237 + result[2, 0, 0, 0] = x116 + result[2, 0, 0, 1] = x164 + result[2, 0, 0, 2] = x174 + result[2, 0, 1, 0] = x164 + result[2, 0, 1, 1] = x203 + result[2, 0, 1, 2] = x213 + result[2, 0, 2, 0] = x174 + result[2, 0, 2, 1] = x213 + result[2, 0, 2, 2] = x216 + result[2, 1, 0, 0] = x164 + result[2, 1, 0, 1] = x203 + result[2, 1, 0, 2] = x213 + result[2, 1, 1, 0] = x203 + result[2, 1, 1, 1] = x230 + result[2, 1, 1, 2] = x236 + result[2, 1, 2, 0] = x213 + result[2, 1, 2, 1] = x236 + result[2, 1, 2, 2] = x237 + result[2, 2, 0, 0] = x174 + result[2, 2, 0, 1] = x213 + result[2, 2, 0, 2] = x216 + result[2, 2, 1, 0] = x213 + result[2, 2, 1, 1] = x236 + result[2, 2, 1, 2] = x237 + result[2, 2, 2, 0] = x216 + result[2, 2, 2, 1] = x237 + result[2, 2, 2, 2] = ( + -6 * x2 * x82 * (5 * x168 - 3) + - x210 * x25 * (3 * x168 - 1) + - 2 * x214 * x234 + - x22 * (-30 * x168 + x17 * x238 + 3) + - x74 + * ( + 18 * x167 * x38 + - 6 * x169 + - x18 * x206 + + x18 * x239 + + x18 * x240 + + 36 * x2 * x36 + - x2 * x67 + + x206 + - x238 * x59 + - x238 * x60 + - x239 + - 4 * x240 + + x241 * x70 + + x241 * x71 + + x72 + ) + ) + return result + + +def T_damp_thole_5(x, y, z, a): + shape = (3,)*5 + result = np.zeros(shape, dtype=np.float64) + x0 = x ** 2 + x1 = y ** 2 + x2 = z ** 2 + x3 = x0 + x1 + x2 + x4 = 1 / x3 + x5 = x0 * x4 + x6 = x ** 4 + x7 = x3 ** (-2) + x8 = x6 * x7 + x9 = -30 * x5 + 35 * x8 + 3 + x10 = x3 ** (-3) + x11 = a * x10 + x12 = sqrt(x3) + x13 = a * x12 + x14 = exp(-x13) + x15 = x13 + 1 + x16 = x14 * x15 + x17 = x11 * x16 + x18 = 15 * x17 / 2 + x19 = x3 ** (-7 / 2) + x20 = x13 + 2 + x21 = x14 * x20 - 2 + x22 = x19 * x21 + x23 = 15 * x22 / 2 + x24 = x3 ** (-3 / 2) + x25 = x0 * x24 + x26 = x20 * x25 + x27 = a * x4 + x28 = 2 * x27 + x29 = x20 * x27 + x30 = 1 / x12 + x31 = -x20 * x30 + x30 + x32 = -x0 * x28 + x0 * x29 - x25 + x26 + x31 + x33 = x14 * x32 + x34 = x3 ** (-5 / 2) + x35 = a * x34 + x36 = 5 * x5 - 3 + x37 = 15 * x36 + x38 = 3 * x5 - 1 + x39 = a * x7 + x40 = 6 * x39 + x41 = a ** 2 + x42 = 3 * x24 + x43 = x41 * x42 + x44 = x26 * x41 + x45 = 3 * x39 + x46 = x20 * x45 + x47 = 3 * x34 + x48 = x0 * x47 + x49 = x20 * x48 + x50 = -x48 + x49 + x51 = x20 * x42 + x52 = 6 * x27 + x53 = 3 * x29 + x54 = x42 - x51 + x52 - x53 + x55 = -x0 * x40 - x0 * x43 + x0 * x46 + x44 + x50 + x54 + x56 = a * x24 + x57 = x14 * x56 + x58 = 5 * x57 + x59 = 15 * x19 + x60 = x59 * x6 + x61 = 18 * x34 + x62 = x0 * x61 + x63 = 30 * x11 + x64 = x41 * x6 + x65 = x20 * x62 + x66 = a ** 3 + x67 = x66 * x8 + x68 = x25 * x41 + x69 = x0 * x39 + x70 = 18 * x39 + x71 = x20 * x70 + x72 = x34 * x41 + x73 = 6 * x20 + x74 = x72 * x73 + x75 = x11 * x20 + x76 = 15 * x75 + x77 = -x42 + x51 - x52 + x53 + x78 = ( + -x0 * x71 + + x20 * x60 + + x20 * x67 + - 6 * x44 + - x6 * x63 + + x6 * x74 + + x6 * x76 + - x60 + - x61 * x64 + + x62 + - x65 + - 4 * x67 + + 18 * x68 + + 36 * x69 + + x77 + ) + x79 = x57 * x78 + x80 = x3 ** (-9 / 2) + x81 = 105 * x80 + x82 = x6 * x81 + x83 = x0 * x19 + x84 = 150 * x83 + x85 = a / x3 ** 4 + x86 = x6 * x85 + x87 = x19 * x64 + x88 = x10 * x66 + x89 = x6 * x88 + x90 = a ** 4 + x91 = x34 * x90 + x92 = x6 * x91 + x93 = x66 * x7 + x94 = x0 * x93 + x95 = x20 * x82 + x96 = x0 * x72 + x97 = x0 * x11 + x98 = x20 * x97 + x99 = 60 * x20 + x100 = 10 * x20 + x101 = 45 * x19 + x102 = x101 * x20 + x103 = x20 * x86 + x104 = 45 * x34 + x105 = x24 * x41 + x106 = x20 * x24 + x107 = x106 * x41 + x108 = x20 * x39 + x109 = x104 * x20 - x104 - 45 * x105 + 15 * x107 + 45 * x108 - 90 * x39 + x110 = x14 / 2 + x111 = a * x30 + x112 = x110 * x111 + x113 = a * x14 + x114 = 60 * x113 + x115 = x15 / x3 ** 5 + x116 = x114 * x115 + x117 = x0 * x85 + x118 = x117 * x16 + x119 = x21 * x80 + x120 = x0 * x119 + x121 = 30 * x120 + x122 = 3 * x14 / 2 + x123 = x11 * x122 * x15 + x124 = x106 - x24 + x125 = x124 - x28 + x29 + x126 = x14 * x35 + x127 = x125 * x126 + x128 = 6 * x127 + x129 = a * x33 + x130 = 18 * x83 + x131 = 9 * x34 + x132 = a * x131 + x133 = x33 * x38 + x134 = x41 * x7 + x135 = 3 * x134 + x136 = 5 * x69 + x137 = 2 * x108 + x138 = x0 * x137 + x139 = x124 - x27 + x140 = x136 - x138 + x139 + x48 - x49 + x68 + x141 = x14 * x140 + x142 = x0 * x55 + x143 = 6 * x126 + x144 = x134 * x14 + x145 = 2 * x144 + x146 = 12 * x96 + x147 = 27 * x97 + x148 = 12 * x98 + x149 = x41 * x49 + x150 = 15 * x83 + x151 = x150 * x20 + x152 = x150 - x151 + x153 = -x43 + x154 = 15 * x39 + x155 = x20 * x40 + x156 = x131 * x20 + x157 = -x131 + x156 + x158 = x153 - x154 + x155 + x157 + x159 = x14 * (x146 + x147 - x148 - x149 + x152 + x158 + x94) + x160 = a * x159 + x161 = x4 * x41 + x162 = x110 * x161 + x163 = 90 * x83 + x164 = 30 * x20 + x165 = 4 * x20 + x166 = x131 - x156 + x167 = x154 - x155 + x43 + x168 = x166 + x167 + x169 = ( + -a * x141 * x38 * x42 + + x0 * x128 * x36 + - x112 + * ( + -90 * x103 + + x163 * x20 + - x163 + - x164 * x87 + - x165 * x89 + + x168 + + x41 * x65 + + x82 + + 195 * x86 + + 105 * x87 + + 22 * x89 + + x92 + - 6 * x94 + - x95 + - 72 * x96 + - 162 * x97 + + 72 * x98 + ) + + x116 * x6 + + 30 * x118 * x36 + + x121 * (7 * x5 - 3) + + x123 * x9 + + x129 * x130 + + x132 * x133 + + x133 * x135 + + x142 * x143 + + x142 * x145 + - 2 * x160 * x25 + + x162 * x78 + + x23 * x9 + + x79 / 2 + ) + x170 = x169 * y + x171 = x169 * z + x172 = x21 / x3 ** (11 / 2) + x173 = x0 * x172 + x174 = 210 * x173 + x175 = x1 * x101 + x176 = 69 * x11 + x177 = 15 * x72 + x178 = x175 * x20 + x179 = 7 * x88 + x180 = x0 * x179 + x181 = 24 * x75 + x182 = x41 * x83 + x183 = 72 * x182 + x184 = 177 * x117 + x185 = x117 * x20 + x186 = 72 * x185 + x187 = x1 * x151 + x188 = x0 * x81 + x189 = x1 * x188 + x190 = -x189 * x20 + x189 + x191 = -x150 + x151 + x192 = -x146 - x147 + x148 + x149 + x168 + x191 - x94 + x193 = 105 * x119 / 2 + x194 = x193 * x36 + x195 = x1 * x24 + x196 = x195 * x66 + x197 = x110 * x55 + x198 = x1 * x4 + x199 = x159 * x41 + x200 = x132 * x141 + x201 = x135 * x141 + x202 = x122 * x55 + x203 = x202 * x35 + x204 = x1 * x106 + x205 = -x1 * x28 + x1 * x29 - x195 + x204 + x31 + x206 = x122 * x35 + x207 = x205 * x206 + x208 = -x106 + x24 + x209 = x208 + x28 - x29 + x210 = x1 * x40 + x211 = x1 * x43 + x212 = x204 * x41 + x213 = x1 * x45 + x214 = x20 * x213 + x215 = x1 * x47 + x216 = x20 * x215 + x217 = -x215 + x216 + x218 = -x210 - x211 + x212 + x214 + x217 + x219 = x209 + x218 + x220 = 3 * x57 / 2 + x221 = x220 * x38 + x222 = x134 * x202 + x223 = 3 * x33 / 2 + x224 = x34 * x66 + x225 = x223 * x224 + x226 = x10 * x41 + x227 = x226 * x33 + x228 = 21 * x227 / 2 + x229 = 45 * x19 / 2 + x230 = x129 * x229 + x231 = x132 * x14 + x232 = x125 * x231 + x233 = x232 * x38 + x234 = x16 * x85 + x235 = x234 * x37 + x236 = x113 * x125 + x237 = x130 * x236 + x238 = x0 * x1 + x239 = x113 * x115 + x240 = 120 * x239 + x241 = -x134 * x223 + x242 = 45 * x234 / 2 + x243 = x242 * x38 + x244 = x1 * x243 + x245 = x241 + x244 + x246 = 5 * x96 + x247 = 23 * x97 + x248 = 8 * x98 + x249 = x208 + x217 + x27 + x250 = -x136 + x138 + x50 - x68 + x251 = x1 * x150 + x1 * x246 + x1 * x247 - x1 * x248 - x187 - x213 + x249 + x250 + x252 = 9 * x33 / 2 + x253 = 9 * x17 / 2 + x254 = -x252 * x35 - x253 * x38 + x255 = -x220 * x251 + x254 + x256 = x57 / 2 + x257 = -9 * x118 - 15 * x120 - x162 * x55 - x23 * x36 - x256 * x55 + x258 = x * ( + x1 * x174 + + x1 * x194 + - x1 * x200 + - x1 * x201 + + x1 * x203 + + x1 * x222 + + x1 * x225 + + x1 * x228 + + x1 * x230 + + x1 * x233 + + x1 * x235 + + x1 * x237 + - x112 + * ( + -x1 * x176 + - x1 * x177 + + x1 * x180 + + x1 * x181 + + x1 * x183 + + x1 * x184 + - x1 * x186 + - x175 + + x178 + - x187 * x41 + + x190 + + x192 + ) + - x160 * x195 + + x196 * x197 + - x198 * x199 + + x207 * x36 + + x219 * x221 + + x238 * x240 + + x245 + + x255 + + x257 + ) + x259 = x0 * x240 + x260 = x20 * x47 + x261 = x260 - x47 + x262 = x107 + x153 + x261 - x40 + x46 + x263 = x220 * x262 + x264 = x152 + x246 + x247 - x248 + x261 - x45 + x265 = x24 * x66 + x266 = x110 * x265 + x267 = x188 * x20 + x268 = x151 * x41 + x269 = x * y * z + x270 = x269 * ( + -x112 + * ( + -x101 + + x102 + - x176 + - x177 + + x180 + + x181 + + x183 + + x184 + - x186 + + x188 + - x267 + - x268 + ) + + x125 * x206 * x36 + - x159 * x161 + - x159 * x56 + + x174 + + x194 + - x200 + - x201 + + x203 + - x220 * x264 + + x222 + + x225 + + x228 + + x230 + + x233 + + x235 + + x237 + + x243 + + x259 + + x263 * x38 + + x266 * x55 + ) + x271 = x188 * x2 - x2 * x267 + x272 = x101 * x2 + x273 = x102 * x2 + x274 = -x176 * x2 - x177 * x2 + x181 * x2 - x272 + x273 + x275 = x2 * x24 + x276 = x275 * x66 + x277 = x2 * x4 + x278 = x106 * x2 - x2 * x28 + x2 * x29 - x275 + x31 + x279 = x206 * x278 + x280 = x2 * x40 + x281 = x2 * x43 + x282 = x107 * x2 + x283 = x2 * x46 + x284 = x2 * x47 + x285 = x2 * x260 + x286 = -x284 + x285 + x287 = -x280 - x281 + x282 + x283 + x286 + x288 = x209 + x287 + x289 = x2 * x243 + x290 = x241 + x289 + x291 = -x2 * x45 + x286 + x292 = ( + x150 * x2 + - x151 * x2 + + x2 * x246 + + x2 * x247 + - x2 * x248 + + x208 + + x250 + + x27 + + x291 + ) + x293 = -x220 * x292 + x254 + x294 = x * ( + -x112 + * ( + x180 * x2 + + x183 * x2 + + x184 * x2 + - x186 * x2 + + x192 + - x2 * x268 + + x271 + + x274 + ) + - x160 * x275 + + x174 * x2 + + x194 * x2 + + x197 * x276 + - x199 * x277 + - x2 * x200 + - x2 * x201 + + x2 * x203 + + x2 * x222 + + x2 * x225 + + x2 * x228 + + x2 * x230 + + x2 * x233 + + x2 * x235 + + x2 * x237 + + x2 * x259 + + x221 * x288 + + x257 + + x279 * x36 + + x290 + + x293 + ) + x295 = 5 * x19 + x296 = x1 * x295 + x297 = x1 * x11 + x298 = 35 * x80 + x299 = x238 * x298 + x300 = x1 * x182 + x301 = x1 * x117 + x302 = x1 * x185 + x303 = -x260 + x47 + x304 = x191 - x246 - x247 + x248 + x45 + x305 = x303 + x304 + x306 = x111 * x122 + x307 = x122 * x161 + x308 = x1 * x59 + x309 = x20 * x308 + x310 = -x308 + x309 + x311 = x166 + x310 + x312 = 9 * x105 + x313 = 9 * x108 + x314 = x41 * x51 + x315 = x312 - x313 - x314 + x70 + x316 = x1 * x63 + x317 = x1 * x61 + x318 = x317 * x41 + x319 = x1 * x93 + x320 = 4 * x319 + x321 = x20 * x319 + x322 = x1 * x72 + x323 = x322 * x73 + x324 = x1 * x76 + x325 = -x316 - x318 - x320 + x321 + x323 + x324 + x326 = x311 + x315 + x325 + x327 = x113 * x25 + x328 = x218 + x54 + x329 = x256 * x328 + x330 = x113 * x83 + x331 = x205 * x330 + x332 = x219 * x231 + x333 = 9 * x205 / 2 + x334 = x126 * x333 + x335 = x0 * x236 + x336 = 45 * x22 / 2 + x337 = 3 * x141 / 2 + x338 = ( + -x0 * x232 + - 54 * x118 + - 90 * x120 + - x134 * x252 + + x140 * x220 + + x161 * x337 + - x223 * x265 + - x336 * x38 + ) + x339 = 315 * x173 + x340 = x1 * x193 + x341 = x110 * x7 * x90 + x342 = x32 * x341 + x343 = x33 * x66 + x344 = 9 * x141 / 2 + x345 = x1 * x344 + x346 = 15 * x1 / 2 + x347 = x129 * x19 + x348 = 195 * x239 + x349 = ( + x1 * x339 + + x1 * x342 + - x134 * x345 + - x196 * x337 + + x215 * x343 + + x227 * x346 + + x238 * x348 + + x340 * x38 + - x345 * x35 + + x346 * x347 + ) + x350 = y * ( + x0 * x332 + + x175 * x335 + + x244 + - x251 * x307 + + x255 + - x306 + * ( + x20 * x296 + - x20 * x299 + - x296 + - 5 * x297 + + x299 + + 11 * x300 + + 51 * x301 + - 16 * x302 + + x305 + ) + + x326 * x327 + + x329 * x38 + + 9 * x331 + + x334 * x38 + + x338 + + x349 + ) + x351 = x303 + x310 + x352 = -x107 + x40 + x43 - x46 + x353 = x325 + x351 + x352 + x354 = x219 * x256 + x355 = x113 * x264 + x356 = x14 * x41 + x357 = x264 * x356 + x358 = x113 * x48 + x359 = x215 * x236 + x360 = x143 * x262 + x361 = 51 * x236 * x83 + x362 = ( + -18 * x118 + - x121 + - x123 * x38 + + x140 * x162 + + x140 * x256 + - x223 * x35 + - x23 * x38 + - x236 * x48 + - x266 * x32 + ) + x363 = z * ( + x1 * x361 + - x112 * (x190 - 15 * x297 + 33 * x300 + 153 * x301 - 48 * x302 + x304 + x351) + - x162 * x251 + - x195 * x355 + - x198 * x357 + + x207 * x38 + + x219 * x358 + + x238 * x360 + + x245 + - x251 * x256 + + x327 * x353 + + 3 * x331 + + x349 + + x354 * x38 + + x359 * x38 + + x362 + ) + x364 = x11 * x2 + x365 = x182 * x2 + x366 = x117 * x2 + x367 = x185 * x2 + x368 = x2 * x59 + x369 = -x368 + x370 = x20 * x368 + x371 = x303 + x369 + x370 + x372 = x2 * x63 + x373 = x2 * x61 + x374 = x373 * x41 + x375 = x2 * x93 + x376 = 4 * x375 + x377 = x20 * x375 + x378 = x2 * x72 + x379 = x378 * x73 + x380 = x2 * x76 + x381 = -x372 - x374 - x376 + x377 + x379 + x380 + x382 = x352 + x371 + x381 + x383 = x256 * x288 + x384 = x278 * x330 + x385 = x236 * x284 + x386 = x143 * x2 + x387 = x193 * x2 + x388 = x2 * x344 + x389 = 15 * x2 / 2 + x390 = ( + x0 * x2 * x348 + - x134 * x388 + + x2 * x339 + + x2 * x342 + + x227 * x389 + - x276 * x337 + + x284 * x343 + + x347 * x389 + - x35 * x388 + + x38 * x387 + ) + x391 = y * ( + x0 * x262 * x386 + - x112 * (x271 + x304 - 15 * x364 + 33 * x365 + 153 * x366 - 48 * x367 + x371) + - x162 * x292 + + x2 * x361 + - x256 * x292 + - x275 * x355 + - x277 * x357 + + x279 * x38 + + x288 * x358 + + x290 + + x327 * x382 + + x362 + + x38 * x383 + + x38 * x385 + + 3 * x384 + + x390 + ) + x392 = x2 * x298 + x393 = x0 * x392 + x394 = x2 * x295 + x395 = x20 * x394 - 5 * x364 - x394 + x396 = x166 + x315 + x369 + x370 + x381 + x397 = x287 + x54 + x398 = x256 * x397 + x399 = x231 * x288 + x400 = 9 * x126 / 2 + x401 = x278 * x400 + x402 = z * ( + x0 * x399 + + x272 * x335 + + x289 + - x292 * x307 + + x293 + - x306 * (-x20 * x393 + x305 + 11 * x365 + 51 * x366 - 16 * x367 + x393 + x395) + + x327 * x396 + + x338 + + x38 * x398 + + x38 * x401 + + 9 * x384 + + x390 + ) + x403 = 315 * x119 + x404 = y ** 4 + x405 = 945 * x172 / 2 + x406 = 90 * x19 + x407 = x1 * x406 + x408 = x20 * x407 + x409 = 108 * x72 + x410 = 90 * x11 + x411 = x1 * x20 + x412 = 36 * x20 + x413 = x157 - x312 + x313 + x314 - x70 + x414 = x404 * x81 + x415 = 210 * x85 + x416 = x404 * x41 + x417 = 135 * x19 + x418 = 40 * x88 + x419 = x404 * x91 + x420 = x20 * x414 + x421 = x404 * x88 + x422 = x404 * x85 + x423 = 105 * x20 + x424 = ( + x100 * x421 + + x102 * x416 + + x20 * x419 + - x404 * x415 + - x404 * x418 + - x414 + - x416 * x417 + - 5 * x419 + + x420 + + x422 * x423 + ) + x425 = x205 * x231 + x426 = x113 * x42 + x427 = 135 * x234 + x428 = x113 * x195 + x429 = 2 * x428 + x430 = x19 * x236 + x431 = 30 * x430 + x432 = x113 * x205 + x433 = 525 * x239 / 2 + x434 = x253 + x336 + x435 = x404 * x59 + x436 = x404 * x93 + x437 = x195 * x41 + x438 = 36 * x39 + x439 = ( + x1 * x438 + - x1 * x71 + - x20 * x317 + + x20 * x435 + + x20 * x436 + - 6 * x212 + + x317 + - x404 * x63 + + x404 * x74 + + x404 * x76 + - x416 * x61 + - x435 + - 4 * x436 + + 18 * x437 + + x77 + ) + x440 = x1 * x328 + x441 = x143 * x440 + x256 * x439 + x442 = x * ( + x1 * x332 + - x1 * x403 + - x1 * x427 + + x112 + * ( + x1 * x409 + + 180 * x297 + + 24 * x319 + - 6 * x321 + - x322 * x412 + + x407 + - x408 + - x410 * x411 + + x413 + + x424 + ) + + x175 * x432 + - x219 * x426 + - x236 * x317 + + x326 * x429 + + x404 * x405 + + x404 * x431 + + x404 * x433 + - x425 + + x434 + + x441 + ) + x443 = x1 * x405 + x444 = x1 * x81 + x445 = x1 * x415 + x446 = x41 * x417 + x447 = x1 * x446 + x448 = x1 * x418 + x449 = x1 * x91 + x450 = 5 * x449 + x451 = x20 * x444 + x452 = x20 * x449 + x453 = x100 * x88 + x454 = x1 * x453 + x455 = x178 * x41 + x456 = x411 * x85 + x457 = 105 * x456 + x458 = x41 * x61 + x459 = x20 * x93 + x460 = x101 - x102 - x20 * x458 + x410 - 3 * x459 + 54 * x72 - 45 * x75 + 12 * x93 + x461 = x206 * x328 + x462 = x113 * x229 + x463 = x205 * x462 + x464 = x1 * x400 + x465 = 105 * x430 / 2 + x466 = -315 * x119 / 2 - 27 * x127 / 2 - 135 * x234 / 2 - x263 + x467 = x269 * ( + x1 * x433 + + x1 * x465 + + x112 + * (-x444 - x445 - x447 - x448 - x450 + x451 + x452 + x454 + x455 + x457 + x460) + + x220 * x353 + + x262 * x464 + + x326 * x57 + + x332 + + x443 + + x461 + + x463 + + x466 + ) + x468 = x2 * x242 + x469 = x1 * x2 + x470 = x19 * x469 + x471 = x113 * x278 + x472 = x19 * x471 + x473 = x19 * x432 + x474 = x389 * x473 + x475 = x2 * x400 + x476 = x113 * x275 + x477 = x2 * x308 + x478 = x2 * x309 + x479 = x215 - x216 + x480 = x2 * x444 + x481 = x2 * x451 + x482 = x308 - x309 + x483 = x * ( + -x1 * x242 + + x112 + * ( + -x2 * x445 + - x2 * x447 + - x2 * x448 + - x2 * x450 + + x2 * x452 + + x2 * x454 + + x2 * x455 + + x2 * x457 + + x262 + + x316 + + x318 + + x320 + - x321 + - x323 + - x324 + + x368 + - x370 + + x372 + + x374 + + x376 + - x377 + - x379 + - x380 + - x480 + + x481 + + x482 + ) + + x114 * x125 * x470 + + x123 + + x2 * x443 + - x207 + + x219 * x475 + + x23 + + x256 + * ( + x125 + - x2 * x316 + - x2 * x318 + - x2 * x320 + + x2 * x321 + + x2 * x323 + + x2 * x324 + + x210 + + x211 + - x212 + - x214 + + x280 + + x281 + - x282 + - x283 + + x284 + - x285 + - x477 + + x478 + + x479 + ) + - x279 + + x288 * x464 + - x340 + + x346 * x472 + + x353 * x476 + - x354 + - x359 + + x360 * x469 + + x382 * x428 + - x383 + - x385 + - x387 + + x433 * x469 + - x468 + + x474 + ) + x484 = x2 * x81 + x485 = x2 * x91 + x486 = x423 * x85 + x487 = x269 * ( + x112 + * ( + -x2 * x415 + - x2 * x418 + - x2 * x446 + + x2 * x453 + + x2 * x486 + + x20 * x484 + + x20 * x485 + + x273 * x41 + + x460 + - x484 + - 5 * x485 + ) + + x2 * x405 + + x2 * x433 + + x2 * x465 + + x206 * x397 + + x220 * x382 + + x262 * x475 + + x278 * x462 + + x396 * x57 + + x399 + + x466 + ) + x488 = z ** 4 + x489 = x488 * x59 + x490 = ( + x2 * x438 + - x2 * x71 + - x20 * x373 + + x20 * x489 + + 18 * x275 * x41 + - 6 * x282 + + x373 + - x458 * x488 + + x459 * x488 + - x488 * x63 + + x488 * x74 + + x488 * x76 + - 4 * x488 * x93 + - x489 + + x77 + ) + x491 = x2 * x406 + x492 = x488 * x81 + x493 = x488 * x91 + x494 = ( + x102 * x41 * x488 + + x20 * x492 + + x20 * x493 + - x415 * x488 + - x418 * x488 + - x446 * x488 + + x453 * x488 + + x486 * x488 + - x492 + - 5 * x493 + ) + x495 = ( + x112 + * ( + -x2 * x20 * x410 + + x2 * x409 + - x20 * x491 + + 180 * x364 + + 24 * x375 + - 6 * x377 + - x378 * x412 + + x413 + + x491 + + x494 + ) + + x2 * x399 + - x2 * x403 + - x2 * x427 + - x231 * x278 + - x236 * x373 + + x256 * x490 + + x272 * x471 + - x288 * x426 + + x386 * x397 + + 2 * x396 * x476 + + x405 * x488 + + x431 * x488 + + x433 * x488 + + x434 + ) + x496 = x * x495 + x497 = x404 * x7 + x498 = -30 * x198 + 35 * x497 + 3 + x499 = 5 * x198 - 3 + x500 = 15 * x126 + x501 = 3 * x198 - 1 + x502 = 5 * x57 / 2 + x503 = 150 * x19 + x504 = x1 * x503 + x505 = x1 * x75 + x506 = x1 * x119 + x507 = x19 * x416 + x508 = 5 * x1 * x39 + x509 = x1 * x137 + x510 = x139 + x437 + x479 + x508 - x509 + x511 = 12 * x322 + x512 = 27 * x297 + x513 = 12 * x505 + x514 = x216 * x41 + x515 = x158 + x319 + x482 + x511 + x512 - x513 - x514 + x516 = x135 * x14 + x517 = x1 * x234 + x518 = z * ( + x1 * x128 * x499 + + 18 * x1 * x473 + - x112 + * ( + -x164 * x507 + - x165 * x421 + + x168 + + x20 * x318 + - 90 * x20 * x422 + - 162 * x297 + - 6 * x319 + - 72 * x322 + - x407 + + x408 + + x414 + + x419 + - x420 + + 22 * x421 + + 195 * x422 + + 72 * x505 + + 105 * x507 + ) + + x116 * x404 + + x123 * x498 + + x145 * x440 + + x162 * x439 + + x205 * x501 * x516 + + x23 * x498 + + x425 * x501 + - x426 * x501 * x510 + - x429 * x515 + + x441 + + 30 * x499 * x517 + + 30 * x506 * (7 * x198 - 3) + ) + x519 = x172 * x469 + x520 = x122 * x205 + x521 = x41 * x470 + x522 = x469 * x85 + x523 = x2 * x456 + x524 = x2 * x510 + x525 = x2 * x205 + x526 = x14 * x226 + x527 = 5 * x322 + x528 = 23 * x297 + x529 = 8 * x505 + x530 = ( + x2 * x527 + + x2 * x528 + - x2 * x529 + + x249 + + x291 + - x437 + + x477 + - x478 + - x508 + + x509 + ) + x531 = -x220 * x530 - x253 * x501 - x334 + x468 * x501 + x532 = y * ( + x110 * x276 * x328 + - x112 + * ( + x167 + + x179 * x469 + + x274 + + x311 + - x319 + - x41 * x478 + + x480 + - x481 + - x511 + - x512 + + x513 + + x514 + + 72 * x521 + + 177 * x522 + - 72 * x523 + ) + + x122 * x134 * x2 * x328 + - x134 * x520 + - x162 * x328 + + x2 * x224 * x520 + + x2 * x232 * x501 + + 15 * x2 * x234 * x499 + + x2 * x461 + + x2 * x463 + + x220 * x288 * x501 + - x23 * x499 + - x231 * x524 + + 18 * x236 * x470 + + x240 * x469 + - x277 * x356 * x515 + + x279 * x499 + - x329 + + x387 * x499 + - x476 * x515 + - 15 * x506 + - x516 * x524 + - 9 * x517 + + 210 * x519 + + 21 * x525 * x526 / 2 + + x531 + ) + x533 = x1 * x392 + x534 = z * ( + -x1 * x232 + + x1 * x399 + + 9 * x1 * x472 + - x122 * x276 * x510 + + x14 * x205 * x284 * x66 + - x144 * x333 + - 9 * x144 * x524 / 2 + + x175 * x2 * x236 + + x205 * x389 * x526 + + x220 * x510 + - x265 * x520 + - x306 + * ( + -x20 * x533 + + x351 + + x395 + + x45 + + 11 * x521 + + 51 * x522 + - 16 * x523 + - x527 + - x528 + + x529 + + x533 + ) + + x307 * x510 + - x307 * x530 + - x336 * x501 + + x341 * x525 + + x348 * x469 + + x387 * x501 + + x396 * x428 + + x398 * x501 + + x401 * x501 + + x474 + - x475 * x510 + - 90 * x506 + - 54 * x517 + + 315 * x519 + + x531 + ) + x535 = x495 * y + x536 = x488 * x7 + x537 = x2 * x503 + result[0, 0, 0, 0, 0] = x * ( + x112 + * ( + x100 * x89 + - x100 * x94 + + x102 * x64 + + 105 * x103 + + x109 + - x20 * x84 + + x20 * x92 + - x82 + + x84 + - 210 * x86 + - 135 * x87 + - 40 * x89 + - 5 * x92 + + 40 * x94 + + x95 + - x96 * x99 + + 180 * x96 + + 300 * x97 + - 150 * x98 + ) + + x18 * x9 + + x23 * (-70 * x5 + 63 * x8 + 15) + + x33 * x35 * x37 + + x38 * x55 * x58 + + 5 * x79 / 2 + ) + result[0, 0, 0, 0, 1] = x170 + result[0, 0, 0, 0, 2] = x171 + result[0, 0, 0, 1, 0] = x170 + result[0, 0, 0, 1, 1] = x258 + result[0, 0, 0, 1, 2] = x270 + result[0, 0, 0, 2, 0] = x171 + result[0, 0, 0, 2, 1] = x270 + result[0, 0, 0, 2, 2] = x294 + result[0, 0, 1, 0, 0] = x170 + result[0, 0, 1, 0, 1] = x258 + result[0, 0, 1, 0, 2] = x270 + result[0, 0, 1, 1, 0] = x258 + result[0, 0, 1, 1, 1] = x350 + result[0, 0, 1, 1, 2] = x363 + result[0, 0, 1, 2, 0] = x270 + result[0, 0, 1, 2, 1] = x363 + result[0, 0, 1, 2, 2] = x391 + result[0, 0, 2, 0, 0] = x171 + result[0, 0, 2, 0, 1] = x270 + result[0, 0, 2, 0, 2] = x294 + result[0, 0, 2, 1, 0] = x270 + result[0, 0, 2, 1, 1] = x363 + result[0, 0, 2, 1, 2] = x391 + result[0, 0, 2, 2, 0] = x294 + result[0, 0, 2, 2, 1] = x391 + result[0, 0, 2, 2, 2] = x402 + result[0, 1, 0, 0, 0] = x170 + result[0, 1, 0, 0, 1] = x258 + result[0, 1, 0, 0, 2] = x270 + result[0, 1, 0, 1, 0] = x258 + result[0, 1, 0, 1, 1] = x350 + result[0, 1, 0, 1, 2] = x363 + result[0, 1, 0, 2, 0] = x270 + result[0, 1, 0, 2, 1] = x363 + result[0, 1, 0, 2, 2] = x391 + result[0, 1, 1, 0, 0] = x258 + result[0, 1, 1, 0, 1] = x350 + result[0, 1, 1, 0, 2] = x363 + result[0, 1, 1, 1, 0] = x350 + result[0, 1, 1, 1, 1] = x442 + result[0, 1, 1, 1, 2] = x467 + result[0, 1, 1, 2, 0] = x363 + result[0, 1, 1, 2, 1] = x467 + result[0, 1, 1, 2, 2] = x483 + result[0, 1, 2, 0, 0] = x270 + result[0, 1, 2, 0, 1] = x363 + result[0, 1, 2, 0, 2] = x391 + result[0, 1, 2, 1, 0] = x363 + result[0, 1, 2, 1, 1] = x467 + result[0, 1, 2, 1, 2] = x483 + result[0, 1, 2, 2, 0] = x391 + result[0, 1, 2, 2, 1] = x483 + result[0, 1, 2, 2, 2] = x487 + result[0, 2, 0, 0, 0] = x171 + result[0, 2, 0, 0, 1] = x270 + result[0, 2, 0, 0, 2] = x294 + result[0, 2, 0, 1, 0] = x270 + result[0, 2, 0, 1, 1] = x363 + result[0, 2, 0, 1, 2] = x391 + result[0, 2, 0, 2, 0] = x294 + result[0, 2, 0, 2, 1] = x391 + result[0, 2, 0, 2, 2] = x402 + result[0, 2, 1, 0, 0] = x270 + result[0, 2, 1, 0, 1] = x363 + result[0, 2, 1, 0, 2] = x391 + result[0, 2, 1, 1, 0] = x363 + result[0, 2, 1, 1, 1] = x467 + result[0, 2, 1, 1, 2] = x483 + result[0, 2, 1, 2, 0] = x391 + result[0, 2, 1, 2, 1] = x483 + result[0, 2, 1, 2, 2] = x487 + result[0, 2, 2, 0, 0] = x294 + result[0, 2, 2, 0, 1] = x391 + result[0, 2, 2, 0, 2] = x402 + result[0, 2, 2, 1, 0] = x391 + result[0, 2, 2, 1, 1] = x483 + result[0, 2, 2, 1, 2] = x487 + result[0, 2, 2, 2, 0] = x402 + result[0, 2, 2, 2, 1] = x487 + result[0, 2, 2, 2, 2] = x496 + result[1, 0, 0, 0, 0] = x170 + result[1, 0, 0, 0, 1] = x258 + result[1, 0, 0, 0, 2] = x270 + result[1, 0, 0, 1, 0] = x258 + result[1, 0, 0, 1, 1] = x350 + result[1, 0, 0, 1, 2] = x363 + result[1, 0, 0, 2, 0] = x270 + result[1, 0, 0, 2, 1] = x363 + result[1, 0, 0, 2, 2] = x391 + result[1, 0, 1, 0, 0] = x258 + result[1, 0, 1, 0, 1] = x350 + result[1, 0, 1, 0, 2] = x363 + result[1, 0, 1, 1, 0] = x350 + result[1, 0, 1, 1, 1] = x442 + result[1, 0, 1, 1, 2] = x467 + result[1, 0, 1, 2, 0] = x363 + result[1, 0, 1, 2, 1] = x467 + result[1, 0, 1, 2, 2] = x483 + result[1, 0, 2, 0, 0] = x270 + result[1, 0, 2, 0, 1] = x363 + result[1, 0, 2, 0, 2] = x391 + result[1, 0, 2, 1, 0] = x363 + result[1, 0, 2, 1, 1] = x467 + result[1, 0, 2, 1, 2] = x483 + result[1, 0, 2, 2, 0] = x391 + result[1, 0, 2, 2, 1] = x483 + result[1, 0, 2, 2, 2] = x487 + result[1, 1, 0, 0, 0] = x258 + result[1, 1, 0, 0, 1] = x350 + result[1, 1, 0, 0, 2] = x363 + result[1, 1, 0, 1, 0] = x350 + result[1, 1, 0, 1, 1] = x442 + result[1, 1, 0, 1, 2] = x467 + result[1, 1, 0, 2, 0] = x363 + result[1, 1, 0, 2, 1] = x467 + result[1, 1, 0, 2, 2] = x483 + result[1, 1, 1, 0, 0] = x350 + result[1, 1, 1, 0, 1] = x442 + result[1, 1, 1, 0, 2] = x467 + result[1, 1, 1, 1, 0] = x442 + result[1, 1, 1, 1, 1] = y * ( + x112 + * ( + x109 + - x20 * x504 + + 300 * x297 + + 40 * x319 + - 10 * x321 + - x322 * x99 + + 180 * x322 + + x424 + + x504 + - 150 * x505 + ) + + x18 * x498 + + x205 * x499 * x500 + + x23 * (-70 * x198 + 63 * x497 + 15) + + x328 * x501 * x58 + + x439 * x502 + ) + result[1, 1, 1, 1, 2] = x518 + result[1, 1, 1, 2, 0] = x467 + result[1, 1, 1, 2, 1] = x518 + result[1, 1, 1, 2, 2] = x532 + result[1, 1, 2, 0, 0] = x363 + result[1, 1, 2, 0, 1] = x467 + result[1, 1, 2, 0, 2] = x483 + result[1, 1, 2, 1, 0] = x467 + result[1, 1, 2, 1, 1] = x518 + result[1, 1, 2, 1, 2] = x532 + result[1, 1, 2, 2, 0] = x483 + result[1, 1, 2, 2, 1] = x532 + result[1, 1, 2, 2, 2] = x534 + result[1, 2, 0, 0, 0] = x270 + result[1, 2, 0, 0, 1] = x363 + result[1, 2, 0, 0, 2] = x391 + result[1, 2, 0, 1, 0] = x363 + result[1, 2, 0, 1, 1] = x467 + result[1, 2, 0, 1, 2] = x483 + result[1, 2, 0, 2, 0] = x391 + result[1, 2, 0, 2, 1] = x483 + result[1, 2, 0, 2, 2] = x487 + result[1, 2, 1, 0, 0] = x363 + result[1, 2, 1, 0, 1] = x467 + result[1, 2, 1, 0, 2] = x483 + result[1, 2, 1, 1, 0] = x467 + result[1, 2, 1, 1, 1] = x518 + result[1, 2, 1, 1, 2] = x532 + result[1, 2, 1, 2, 0] = x483 + result[1, 2, 1, 2, 1] = x532 + result[1, 2, 1, 2, 2] = x534 + result[1, 2, 2, 0, 0] = x391 + result[1, 2, 2, 0, 1] = x483 + result[1, 2, 2, 0, 2] = x487 + result[1, 2, 2, 1, 0] = x483 + result[1, 2, 2, 1, 1] = x532 + result[1, 2, 2, 1, 2] = x534 + result[1, 2, 2, 2, 0] = x487 + result[1, 2, 2, 2, 1] = x534 + result[1, 2, 2, 2, 2] = x535 + result[2, 0, 0, 0, 0] = x171 + result[2, 0, 0, 0, 1] = x270 + result[2, 0, 0, 0, 2] = x294 + result[2, 0, 0, 1, 0] = x270 + result[2, 0, 0, 1, 1] = x363 + result[2, 0, 0, 1, 2] = x391 + result[2, 0, 0, 2, 0] = x294 + result[2, 0, 0, 2, 1] = x391 + result[2, 0, 0, 2, 2] = x402 + result[2, 0, 1, 0, 0] = x270 + result[2, 0, 1, 0, 1] = x363 + result[2, 0, 1, 0, 2] = x391 + result[2, 0, 1, 1, 0] = x363 + result[2, 0, 1, 1, 1] = x467 + result[2, 0, 1, 1, 2] = x483 + result[2, 0, 1, 2, 0] = x391 + result[2, 0, 1, 2, 1] = x483 + result[2, 0, 1, 2, 2] = x487 + result[2, 0, 2, 0, 0] = x294 + result[2, 0, 2, 0, 1] = x391 + result[2, 0, 2, 0, 2] = x402 + result[2, 0, 2, 1, 0] = x391 + result[2, 0, 2, 1, 1] = x483 + result[2, 0, 2, 1, 2] = x487 + result[2, 0, 2, 2, 0] = x402 + result[2, 0, 2, 2, 1] = x487 + result[2, 0, 2, 2, 2] = x496 + result[2, 1, 0, 0, 0] = x270 + result[2, 1, 0, 0, 1] = x363 + result[2, 1, 0, 0, 2] = x391 + result[2, 1, 0, 1, 0] = x363 + result[2, 1, 0, 1, 1] = x467 + result[2, 1, 0, 1, 2] = x483 + result[2, 1, 0, 2, 0] = x391 + result[2, 1, 0, 2, 1] = x483 + result[2, 1, 0, 2, 2] = x487 + result[2, 1, 1, 0, 0] = x363 + result[2, 1, 1, 0, 1] = x467 + result[2, 1, 1, 0, 2] = x483 + result[2, 1, 1, 1, 0] = x467 + result[2, 1, 1, 1, 1] = x518 + result[2, 1, 1, 1, 2] = x532 + result[2, 1, 1, 2, 0] = x483 + result[2, 1, 1, 2, 1] = x532 + result[2, 1, 1, 2, 2] = x534 + result[2, 1, 2, 0, 0] = x391 + result[2, 1, 2, 0, 1] = x483 + result[2, 1, 2, 0, 2] = x487 + result[2, 1, 2, 1, 0] = x483 + result[2, 1, 2, 1, 1] = x532 + result[2, 1, 2, 1, 2] = x534 + result[2, 1, 2, 2, 0] = x487 + result[2, 1, 2, 2, 1] = x534 + result[2, 1, 2, 2, 2] = x535 + result[2, 2, 0, 0, 0] = x294 + result[2, 2, 0, 0, 1] = x391 + result[2, 2, 0, 0, 2] = x402 + result[2, 2, 0, 1, 0] = x391 + result[2, 2, 0, 1, 1] = x483 + result[2, 2, 0, 1, 2] = x487 + result[2, 2, 0, 2, 0] = x402 + result[2, 2, 0, 2, 1] = x487 + result[2, 2, 0, 2, 2] = x496 + result[2, 2, 1, 0, 0] = x391 + result[2, 2, 1, 0, 1] = x483 + result[2, 2, 1, 0, 2] = x487 + result[2, 2, 1, 1, 0] = x483 + result[2, 2, 1, 1, 1] = x532 + result[2, 2, 1, 1, 2] = x534 + result[2, 2, 1, 2, 0] = x487 + result[2, 2, 1, 2, 1] = x534 + result[2, 2, 1, 2, 2] = x535 + result[2, 2, 2, 0, 0] = x402 + result[2, 2, 2, 0, 1] = x487 + result[2, 2, 2, 0, 2] = x496 + result[2, 2, 2, 1, 0] = x487 + result[2, 2, 2, 1, 1] = x534 + result[2, 2, 2, 1, 2] = x535 + result[2, 2, 2, 2, 0] = x496 + result[2, 2, 2, 2, 1] = x535 + result[2, 2, 2, 2, 2] = z * ( + x112 + * ( + x109 + - 150 * x2 * x75 + - x20 * x537 + + 300 * x364 + + 40 * x375 + - 10 * x377 + - x378 * x99 + + 180 * x378 + + x494 + + x537 + ) + + x18 * (-30 * x277 + 35 * x536 + 3) + + x23 * (-70 * x277 + 63 * x536 + 15) + + x278 * x500 * (5 * x277 - 3) + + x397 * x58 * (3 * x277 - 1) + + x490 * x502 + ) + return result + + +T_damp_thole = [ + T_damp_thole_0, + T_damp_thole_1, + T_damp_thole_2, + T_damp_thole_3, + T_damp_thole_4, + T_damp_thole_5, +] From 202e3b52d13d0ccaa5ab2f7deab3573aba540916 Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Fri, 27 Sep 2019 16:47:25 +0200 Subject: [PATCH 3/6] add default parameters to py interface, add tests for damped Tensor (broken) --- cppe/python_iface/export_math.cc | 10 +++++++--- tests/test_math.py | 16 +++++++++++----- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/cppe/python_iface/export_math.cc b/cppe/python_iface/export_math.cc index 00a243d2..d0d5f906 100644 --- a/cppe/python_iface/export_math.cc +++ b/cppe/python_iface/export_math.cc @@ -11,8 +11,12 @@ using namespace libcppe; void export_math(py::module& m) { m.def("smat_vec", &smat_vec); m.def("multipole_derivative", &multipole_derivative); - m.def("T", &T); - m.def("Tk_tensor", &Tk_tensor); + m.def("T", &T, py::arg("Rij"), py::arg("x"), py::arg("y"), py::arg("z"), + py::arg("Cijn"), py::arg("damping_factor") = 0.0, py::arg("alpha_i") = 0.0, + py::arg("alpha_j") = 0.0); + m.def("Tk_tensor", &Tk_tensor, py::arg("k"), py::arg("Rij"), py::arg("Tk_coeffs"), + py::arg("damping_factor") = 0.0, py::arg("alpha_i") = 0.0, + py::arg("alpha_j") = 0.0); m.def("Tk_coefficients", &Tk_coefficients); m.def("xyz2idx", &xyz2idx); m.def("factorial", &factorial); @@ -23,4 +27,4 @@ void export_math(py::module& m) { "Prefactors for the multipole components (alphabetical order)"); m.def("prefactors_nuclei", &prefactors_nuclei); m.def("multipole_components", &multipole_components); -} \ No newline at end of file +} diff --git a/tests/test_math.py b/tests/test_math.py index 88b5a56a..95c9e090 100644 --- a/tests/test_math.py +++ b/tests/test_math.py @@ -74,18 +74,24 @@ def test_T_tensors(self): def test_T_tensors_damped(self): # tests the T tensors against auto-generated Python code - ref_T = tensors.T + ref_T = tensors.T_damp_thole for k in range(4): R = np.random.random(3) - ref = ref_T[k](*R) + damp = 2.0 + a_i = 4.0 + a_j = 10.0 + a = 1/(a_i*a_j)**(1/6)*damp + ref = ref_T[k](*R, a) coeffs = Tk_coefficients(k) # damped tensor: damping_factor, alpha_i, alpha_j - actual = Tk_tensor(k, R, coeffs, 1.0, 1.0, 1.0) + actual = Tk_tensor(k, R, coeffs, damp, a_i, a_j) # gets the indices of non-redundant components # e.g., takes only xy from (xy, yx) and so on ... sym_indices = symmetry.get_symm_indices(k) -# np.testing.assert_almost_equal(actual, -# ref.take(sym_indices), decimal=10) + np.testing.assert_almost_equal( + actual, ref.take(sym_indices), decimal=10, + err_msg="Damped T tensors do not match. Order = {}".format(k) + ) From d5081a3db614c22db705afef95053d2b6c5e6bcd Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Fri, 27 Sep 2019 21:24:36 +0200 Subject: [PATCH 4/6] fixed damped T3 bug --- cppe/core/math.cc | 2 +- tests/test_math.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/cppe/core/math.cc b/cppe/core/math.cc index 1cbe2be7..b02cfef3 100644 --- a/cppe/core/math.cc +++ b/cppe/core/math.cc @@ -141,7 +141,7 @@ double T(const Eigen::Vector3d& Rij, int x, int y, int z, if (k == 1) { Cz *= scr_facs[1]; } else if (k == 3) { - Cz *= scr_facs[3]; + Cz *= scr_facs[2]; } } else if (kk == 2) { if (k == 2) Cz *= scr_facs[2]; diff --git a/tests/test_math.py b/tests/test_math.py index 95c9e090..7fe4ebae 100644 --- a/tests/test_math.py +++ b/tests/test_math.py @@ -70,7 +70,8 @@ def test_T_tensors(self): # e.g., takes only xy from (xy, yx) and so on ... sym_indices = symmetry.get_symm_indices(k) np.testing.assert_almost_equal(actual, - ref.take(sym_indices), decimal=10) + ref.take(sym_indices), decimal=10, + err_msg="T tensors do not match. Order = {}".format(k)) def test_T_tensors_damped(self): # tests the T tensors against auto-generated Python code From a8380c59d06556bd464625169259304eb0436bf0 Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Mon, 30 Sep 2019 18:01:38 +0200 Subject: [PATCH 5/6] implement damping for induced moments and multipole fields --- cppe/core/cppe_state.cc | 2 +- cppe/core/electric_fields.cc | 43 +++++++++++++++++++++-------- cppe/core/electric_fields.hh | 12 +++++--- cppe/core/math.cc | 6 ++-- cppe/core/math.hh | 4 ++- cppe/core/pe_options.hh | 7 +++++ cppe/core/potential.hh | 2 ++ cppe/python_iface/export_fields.cc | 4 +-- cppe/python_iface/export_options.cc | 8 +++++- 9 files changed, 66 insertions(+), 22 deletions(-) diff --git a/cppe/core/cppe_state.cc b/cppe/core/cppe_state.cc index 1dd011f0..64a0c452 100644 --- a/cppe/core/cppe_state.cc +++ b/cppe/core/cppe_state.cc @@ -50,7 +50,7 @@ void CppeState::calculate_static_energies_and_fields() { NuclearFields nfields(m_mol, m_potentials); m_nuc_fields = nfields.compute(); // Multipole fields - MultipoleFields mul_fields(m_potentials); + MultipoleFields mul_fields(m_potentials, m_options); m_multipole_fields = mul_fields.compute(); } diff --git a/cppe/core/electric_fields.cc b/cppe/core/electric_fields.cc index 3584d556..8f1c00f2 100644 --- a/cppe/core/electric_fields.cc +++ b/cppe/core/electric_fields.cc @@ -6,10 +6,7 @@ namespace libcppe { -Eigen::VectorXd NuclearFields::compute(bool damp_core) { - if (damp_core) { - throw std::runtime_error("damping not implemented"); - } +Eigen::VectorXd NuclearFields::compute() { std::vector Tk_coeffs = Tk_coefficients(5); Eigen::VectorXd nuc_fields = Eigen::VectorXd::Zero(3 * m_n_polsites); #pragma omp parallel for firstprivate(Tk_coeffs) @@ -29,10 +26,7 @@ Eigen::VectorXd NuclearFields::compute(bool damp_core) { return nuc_fields; } -Eigen::VectorXd MultipoleFields::compute(bool damp) { - if (damp) { - throw std::runtime_error("damping not implemented"); - } +Eigen::VectorXd MultipoleFields::compute() { std::vector Tk_coeffs = Tk_coefficients(5); Eigen::VectorXd mult_fields = Eigen::VectorXd::Zero(3 * m_n_polsites); // Field at site of potential1 caused by all other sites (also non-polarizable @@ -41,6 +35,8 @@ Eigen::VectorXd MultipoleFields::compute(bool damp) { for (size_t i = 0; i < m_n_polsites; i++) { size_t site_counter = 3 * i; Potential& potential1 = m_polsites[i]; + double alpha_i_isotropic = + potential1.get_polarizabilities()[0].get_isotropic_value(); // for damping for (size_t j = 0; j < m_potentials.size(); j++) { Potential& potential2 = m_potentials[j]; // all other multipoles create el. field at site i @@ -48,14 +44,31 @@ Eigen::VectorXd MultipoleFields::compute(bool damp) { if (potential1.excludes_site(potential2.index)) continue; Eigen::Vector3d diff = potential1.get_site_position() - potential2.get_site_position(); + + // Thole damping: potential2 needs to be polarizable + bool damp_enabled = m_options.damp_multipole; + double alpha_j_isotropic = 0.0; + if (!potential2.is_polarizable()) { + damp_enabled = false; + } else { + alpha_j_isotropic = + potential2.get_polarizabilities()[0].get_isotropic_value(); // for damping + } + // std::cout << "-- created by site " << potential2.index << std::endl; for (auto& mul : potential2.get_multipoles()) { // if (std::all_of(mul.get_values().begin(), mul.get_values().end(), // [](double v) { return std::abs(v) == 0.0; })) { // continue; // } - Eigen::VectorXd Fi = - multipole_derivative(mul.m_k, 1, diff, mul.get_values_vec(), Tk_coeffs); + Eigen::VectorXd Fi; + if (damp_enabled) { + Fi = multipole_derivative(mul.m_k, 1, diff, mul.get_values_vec(), Tk_coeffs, + m_options.damping_factor_multipole, alpha_i_isotropic, + alpha_j_isotropic); + } else { + Fi = multipole_derivative(mul.m_k, 1, diff, mul.get_values_vec(), Tk_coeffs); + } mult_fields(site_counter) += Fi(0); mult_fields(site_counter + 1) += Fi(1); mult_fields(site_counter + 2) += Fi(2); @@ -119,7 +132,15 @@ void InducedMoments::compute(const Eigen::VectorXd& total_fields, continue; } Eigen::Vector3d diff = pot2.get_site_position() - pot1.get_site_position(); - Eigen::VectorXd T2 = Tk_tensor(2, diff, Tk_coeffs); + Eigen::VectorXd T2; + if (m_options.damp_induced) { + Polarizability& alpha_i = pot1.get_polarizabilities()[0]; + Polarizability& alpha_j = pot2.get_polarizabilities()[0]; + T2 = Tk_tensor(2, diff, Tk_coeffs, m_options.damping_factor_induced, + alpha_i.get_isotropic_value(), alpha_j.get_isotropic_value()); + } else { + T2 = Tk_tensor(2, diff, Tk_coeffs); + } Ftmp += smat_vec(T2, induced_moments.segment(m, 3), 1.0); } // keep value to calculate residual diff --git a/cppe/core/electric_fields.hh b/cppe/core/electric_fields.hh index 5da970db..fbd6ce95 100644 --- a/cppe/core/electric_fields.hh +++ b/cppe/core/electric_fields.hh @@ -22,7 +22,7 @@ class ElectricFields { m_n_polsites = m_polsites.size(); }; ~ElectricFields(){}; - virtual Eigen::VectorXd compute(bool damp = false) = 0; + virtual Eigen::VectorXd compute() = 0; }; class NuclearFields : public ElectricFields { @@ -32,13 +32,17 @@ class NuclearFields : public ElectricFields { public: NuclearFields(Molecule mol, std::vector potentials) : ElectricFields(potentials), m_mol(mol){}; - Eigen::VectorXd compute(bool damp = false); + Eigen::VectorXd compute(); }; class MultipoleFields : public ElectricFields { public: - MultipoleFields(std::vector potentials) : ElectricFields(potentials){}; - Eigen::VectorXd compute(bool damp = false); + MultipoleFields(std::vector potentials, PeOptions options) + : ElectricFields(potentials), m_options(options){}; + Eigen::VectorXd compute(); + + private: + PeOptions m_options; }; class InducedMoments { diff --git a/cppe/core/math.cc b/cppe/core/math.cc index b02cfef3..9505234e 100644 --- a/cppe/core/math.cc +++ b/cppe/core/math.cc @@ -40,7 +40,9 @@ Eigen::VectorXd Tk_tensor(int k, const Eigen::Vector3d& Rij, // only 1st derivative supported Eigen::VectorXd multipole_derivative(int k, int l, const Eigen::Vector3d& Rji, Eigen::VectorXd Mkj, - std::vector& Tk_coeffs) { + std::vector& Tk_coeffs, + double damping_factor, double alpha_i, + double alpha_j) { if (l > 1) throw std::runtime_error("Only 1st derivatives supported for multipoles"); Eigen::VectorXd Fi = Eigen::VectorXd::Zero(3); @@ -55,7 +57,7 @@ Eigen::VectorXd multipole_derivative(int k, int l, const Eigen::Vector3d& Rji, double symfac; // std::cout << "mul k = " << k << std::endl; // std::cout << "l = " << l << std::endl; - Eigen::VectorXd Tk = Tk_tensor(k + l, Rji, Tk_coeffs); + Eigen::VectorXd Tk = Tk_tensor(k + l, Rji, Tk_coeffs, damping_factor, alpha_i, alpha_j); for (x = k + l; x > -1; x--) { for (y = k + l; y > -1; y--) { for (z = k + l; z > -1; z--) { diff --git a/cppe/core/math.hh b/cppe/core/math.hh index defd4706..cdf172a8 100644 --- a/cppe/core/math.hh +++ b/cppe/core/math.hh @@ -10,7 +10,9 @@ Eigen::Vector3d smat_vec(const Eigen::VectorXd& mat, const Eigen::Vector3d& vec, Eigen::VectorXd multipole_derivative(int k, int l, const Eigen::Vector3d& Rji, Eigen::VectorXd Mkj, - std::vector& Tk_coeffs); + std::vector& Tk_coeffs, + double damping_factor = 0.0, double alpha_i = 0.0, + double alpha_j = 0.0); double T(const Eigen::Vector3d& Rij, int x, int y, int z, std::vector& Cijn, double damping_factor = 0.0, diff --git a/cppe/core/pe_options.hh b/cppe/core/pe_options.hh index 7cbc4e11..d1a4aa8f 100644 --- a/cppe/core/pe_options.hh +++ b/cppe/core/pe_options.hh @@ -39,6 +39,13 @@ struct PeOptions { double diis_start_norm = 1.0; ///< maximal residual norm for which DIIS is being enabled + bool damp_induced = false; ///< Enable Thole damping for induced moments + double damping_factor_induced = 2.1304; ///< Default damping factor for induced moments + bool damp_multipole = + false; ///< Enable Thole damping for electric fields created by multipole moments + double damping_factor_multipole = 2.1304; ///< Default damping factor for electric + ///< fields created by multipole moments + bool pe_border = false; ///< Activate border options for sites in proximity ///< to the QM/MM border BorderOptions border_options{}; ///< Options for QM/MM border diff --git a/cppe/core/potential.hh b/cppe/core/potential.hh index 3681f0bd..dad0a026 100644 --- a/cppe/core/potential.hh +++ b/cppe/core/potential.hh @@ -59,6 +59,8 @@ class Polarizability { m_values = new_values; } + double get_isotropic_value() { return (m_values[0] + m_values[3] + m_values[5]) / 3.0; } + Eigen::VectorXd get_values_vec() { return Eigen::Map(m_values.data(), m_values.size()); } diff --git a/cppe/python_iface/export_fields.cc b/cppe/python_iface/export_fields.cc index 11bad887..3c9d7961 100644 --- a/cppe/python_iface/export_fields.cc +++ b/cppe/python_iface/export_fields.cc @@ -16,7 +16,7 @@ void export_fields(py::module& m) { py::class_ mul_fields( m, "MultipoleFields", "Electric fields created by multipoles"); - mul_fields.def(py::init>()) + mul_fields.def(py::init, libcppe::PeOptions>()) .def("compute", &libcppe::MultipoleFields::compute); py::class_ ind_moments(m, "InducedMoments"); @@ -28,4 +28,4 @@ void export_fields(py::module& m) { py::arg("total_fields"), py::arg("make_guess")); m.def("get_polarizable_sites", &libcppe::get_polarizable_sites); -} \ No newline at end of file +} diff --git a/cppe/python_iface/export_options.cc b/cppe/python_iface/export_options.cc index c92ddb3b..1f34a977 100644 --- a/cppe/python_iface/export_options.cc +++ b/cppe/python_iface/export_options.cc @@ -29,7 +29,13 @@ void export_options(py::module& m) { .def_readwrite("do_diis", &libcppe::PeOptions::do_diis) .def_readwrite("diis_start_norm", &libcppe::PeOptions::diis_start_norm) .def_readwrite("maxiter", &libcppe::PeOptions::maxiter) + .def_readwrite("damp_induced", &libcppe::PeOptions::damp_induced) + .def_readwrite("damping_factor_induced", + &libcppe::PeOptions::damping_factor_induced) + .def_readwrite("damp_multipole", &libcppe::PeOptions::damp_multipole) + .def_readwrite("damping_factor_multipole", + &libcppe::PeOptions::damping_factor_multipole) .def_readwrite("pe_border", &libcppe::PeOptions::pe_border) .def_readwrite("border_options", &libcppe::PeOptions::border_options); -} \ No newline at end of file +} From 994976e223c98468547ffb58d25f7bef228c39dc Mon Sep 17 00:00:00 2001 From: Maximilian Scheurer Date: Mon, 30 Sep 2019 18:31:50 +0200 Subject: [PATCH 6/6] remove TODO, minor rewording in pe_options --- cppe/core/math.cc | 1 - cppe/core/pe_options.hh | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/cppe/core/math.cc b/cppe/core/math.cc index 9505234e..7a6ba09d 100644 --- a/cppe/core/math.cc +++ b/cppe/core/math.cc @@ -17,7 +17,6 @@ Eigen::Vector3d smat_vec(const Eigen::VectorXd& mat, const Eigen::Vector3d& vec, return result; } -// TODO: add option for damping Eigen::VectorXd Tk_tensor(int k, const Eigen::Vector3d& Rij, std::vector& Tk_coeffs, double damping_factor, double alpha_i, double alpha_j) { diff --git a/cppe/core/pe_options.hh b/cppe/core/pe_options.hh index d1a4aa8f..6fe191d8 100644 --- a/cppe/core/pe_options.hh +++ b/cppe/core/pe_options.hh @@ -40,10 +40,10 @@ struct PeOptions { 1.0; ///< maximal residual norm for which DIIS is being enabled bool damp_induced = false; ///< Enable Thole damping for induced moments - double damping_factor_induced = 2.1304; ///< Default damping factor for induced moments + double damping_factor_induced = 2.1304; ///< damping factor for induced moments bool damp_multipole = false; ///< Enable Thole damping for electric fields created by multipole moments - double damping_factor_multipole = 2.1304; ///< Default damping factor for electric + double damping_factor_multipole = 2.1304; ///< damping factor for electric ///< fields created by multipole moments bool pe_border = false; ///< Activate border options for sites in proximity