From 769870cfc9579d4946d8836077b93d54d69f499b Mon Sep 17 00:00:00 2001 From: Markus Hoehnerbach Date: Thu, 29 Jun 2017 15:50:54 +0200 Subject: [PATCH 1/9] Proper spline coefficient calculation for AIREBO --- src/MANYBODY/pair_airebo.cpp | 326 +++++++++++++++++++++++++++++++++++ src/MANYBODY/pair_airebo.h | 6 + 2 files changed, 332 insertions(+) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 0ca80c6b76..b8ef4db48c 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -4066,6 +4066,276 @@ double PairAIREBO::Sptricubic(double x, double y, double z, return f; } +/* ---------------------------------------------------------------------- + spline coefficient matrix python script +------------------------------------------------------------------------- + +import numpy as np +import numpy.linalg as lin + +# Generate all the derivatives that are spline conditions +# Ordered such that df / dx_i / d_xj i < j. +# Gives the derivatives at which the spline's values are prescribed. +def generate_derivs(n): + def generate_derivs_order(n, m): + if m == 0: + return [tuple()] + if m == 1: + return [tuple([i]) for i in range(n)] + rec = generate_derivs_order(n, m - 1) + return [tuple([i]+list(j)) for i in range(n) for j in rec if j[0] > i] + ret = [] + m = 0 + while m <= n: + ret += generate_derivs_order(n, m) + m += 1 + return ret + +# Generate all the points in an n-dimensional unit cube. +# Gives the points at which the spline's values are prescribed. +def generate_points(n): + if n == 1: + return [(0,), (1,)] + rec = generate_points(n - 1) + return [tuple([j]+list(i)) for j in range(2) for i in rec] + +# Generate all the coefficients in the order later expected. +def generate_coeffs(n): + if n == 1: + return [tuple([i]) for i in range(4)] # cubic + rec = generate_coeffs(n-1) + return [tuple([i]+list(j)) for i in range(4) for j in rec] + +# Evaluate the `deriv`'s derivative at `point` symbolically +# with respect to the coefficients `coeffs`. +def eval_at(n, coeffs, deriv, point): + def eval_single(order, value, the_deriv): + if the_deriv: + if order == 0: + return 0 + if order == 1: + return 1 + return order * value + else: + if order == 0: + return 1 + else: + return value + result = {} + for c in coeffs: + result[c] = 1 + for i in range(n): + result[c] *= eval_single(c[i], point[i], i in deriv) + return result + +# Build the matrix transforming prescribed values to coefficients. +def get_matrix(n): + coeffs = generate_coeffs(n) + points = generate_points(n) + derivs = generate_derivs(n) + assert(len(coeffs) == len(points)*len(derivs)) + i = 0 + A = np.zeros((len(coeffs), len(points)*len(derivs))) + for d in derivs: + for p in points: + coeff = eval_at(n, coeffs, d, p) + for j, c in enumerate(coeffs): + A[i, j] = coeff[c] + i += 1 + return lin.inv(A) + +# Output the first k values with padding n from A. +def output_matrix(n, k, A): + print('\n'.join([''.join([("%{}d,".format(n+1)) % i for i in j[:k]]) for j in A])) + +*/ + +/* ---------------------------------------------------------------------- + tricubic spline coefficient calculation +------------------------------------------------------------------------- */ + +void PairAIREBO::Sptricubic_patch_adjust(double * dl, double wid, double lo, + char dir) { + int rowOuterL = 16, rowInnerL = 1, colL; + if (dir == 'R') { + rowOuterL = 4; + colL = 16; + } else if (dir == 'M') { + colL = 4; + } else if (dir == 'L') { + rowInnerL = 4; + colL = 1; + } + double binomial[5] = {1, 1, 2, 6}; + for (int rowOuter = 0; rowOuter < 4; rowOuter++) { + for (int rowInner = 0; rowInner < 4; rowInner++) { + for (int col = 0; col < 4; col++) { + double acc = 0; + for (int k = col; k < 4; k++) { + acc += dl[rowOuterL * rowOuter + rowInnerL * rowInner + colL * k] + * pow(wid, -k) * pow(-lo, k - col) * binomial[k] / binomial[col] + / binomial[k - col]; + } + dl[rowOuterL * rowOuter + rowInnerL * rowInner + colL * col] = acc; + } + } + } +} + +void PairAIREBO::Sptricubic_patch_coeffs( + double xmin, double xmax, double ymin, double ymax, double zmin, double zmax, + double * y, double * y1, double * y2, double * y3, double * dl +) { + const double C_inv[64][32] = { + // output_matrix(2, 8*4, get_matrix(3)) + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, + -3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, -1, 0, 0, 0, 0, 0, 0, + 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + -3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, + 9, -9, -9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, -6, 3, -3, 0, 0, 0, 0, 6, 3, -6, -3, 0, 0, 0, 0, + -6, 6, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, 4, -2, 2, 0, 0, 0, 0, -3, -3, 3, 3, 0, 0, 0, 0, + 2, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, + -6, 6, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, -3, 3, 0, 0, 0, 0, -4, -2, 4, 2, 0, 0, 0, 0, + 4, -4, -4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 2, -2, 0, 0, 0, 0, 2, 2, -2, -2, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, -3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 9, -9, -9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, -6, 6, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, -6, 6, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 4, -4, -4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + -3, 0, 0, 0, 3, 0, 0, 0, -2, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 3, 0, 0, 0, + 9, -9, 0, 0, -9, 9, 0, 0, 6, -6, 0, 0, 3, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0, -6, -3, 0, 0, + -6, 6, 0, 0, 6, -6, 0, 0, -4, 4, 0, 0, -2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, -3, 0, 0, 3, 3, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, -9, 0, 0, -9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 6, 0, 0, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 9, 0, -9, 0, -9, 0, 9, 0, 6, 0, -6, 0, 3, 0, -3, 0, 6, 0, 3, 0, -6, 0, -3, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9, 0, -9, 0, -9, 0, 9, 0, + -27, 27, 27,-27, 27,-27,-27, 27,-18, 18, 18,-18, -9, 9, 9, -9,-18, 18, -9, 9, 18,-18, 9, -9,-18, -9, 18, 9, 18, 9,-18, -9, + 18,-18,-18, 18,-18, 18, 18,-18, 12,-12,-12, 12, 6, -6, -6, 6, 12,-12, 6, -6,-12, 12, -6, 6, 9, 9, -9, -9, -9, -9, 9, 9, + -6, 0, 6, 0, 6, 0, -6, 0, -4, 0, 4, 0, -2, 0, 2, 0, -3, 0, -3, 0, 3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, 6, 0, 6, 0, -6, 0, + 18,-18,-18, 18,-18, 18, 18,-18, 12,-12,-12, 12, 6, -6, -6, 6, 9, -9, 9, -9, -9, 9, -9, 9, 12, 6,-12, -6,-12, -6, 12, 6, + -12, 12, 12,-12, 12,-12,-12, 12, -8, 8, 8, -8, -4, 4, 4, -4, -6, 6, -6, 6, 6, -6, 6, -6, -6, -6, 6, 6, 6, 6, -6, -6, + 2, 0, 0, 0, -2, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, -2, 0, 0, 0, + -6, 6, 0, 0, 6, -6, 0, 0, -3, 3, 0, 0, -3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -4, -2, 0, 0, 4, 2, 0, 0, + 4, -4, 0, 0, -4, 4, 0, 0, 2, -2, 0, 0, 2, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, -2, -2, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 6, 0, 0, 6, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, -4, 0, 0, -4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + -6, 0, 6, 0, 6, 0, -6, 0, -3, 0, 3, 0, -3, 0, 3, 0, -4, 0, -2, 0, 4, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, 6, 0, 6, 0, -6, 0, + 18,-18,-18, 18,-18, 18, 18,-18, 9, -9, -9, 9, 9, -9, -9, 9, 12,-12, 6, -6,-12, 12, -6, 6, 12, 6,-12, -6,-12, -6, 12, 6, + -12, 12, 12,-12, 12,-12,-12, 12, -6, 6, 6, -6, -6, 6, 6, -6, -8, 8, -4, 4, 8, -8, 4, -4, -6, -6, 6, 6, 6, 6, -6, -6, + 4, 0, -4, 0, -4, 0, 4, 0, 2, 0, -2, 0, 2, 0, -2, 0, 2, 0, 2, 0, -2, 0, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, -4, 0, -4, 0, 4, 0, + -12, 12, 12,-12, 12,-12,-12, 12, -6, 6, 6, -6, -6, 6, 6, -6, -6, 6, -6, 6, 6, -6, 6, -6, -8, -4, 8, 4, 8, 4, -8, -4, + 8, -8, -8, 8, -8, 8, 8, -8, 4, -4, -4, 4, 4, -4, -4, 4, 4, -4, 4, -4, -4, 4, -4, 4, 4, 4, -4, -4, -4, -4, 4, 4, + }; + double dx = xmax - xmin; + double dy = ymax - ymin; + double dz = zmax - zmin; + double x[32]; + for (int i = 0; i < 8; i++) { + x[i+0*8] = y[i]; + x[i+1*8] = y1[i] * dx; + x[i+2*8] = y2[i] * dy; + x[i+3*8] = y3[i] * dz; + } + for (int i = 0; i < 64; i++) { + dl[i] = 0; + for (int k = 0; k < 32; k++) { + dl[i] += x[k] * C_inv[i][k]; + } + } + Sptricubic_patch_adjust(dl, dx, xmin, 'R'); + Sptricubic_patch_adjust(dl, dy, ymin, 'M'); + Sptricubic_patch_adjust(dl, dz, zmin, 'L'); +} + +/* ---------------------------------------------------------------------- + bicubic spline coefficient calculation +------------------------------------------------------------------------- */ + +void PairAIREBO::Spbicubic_patch_adjust(double * dl, double wid, double lo, + char dir) { + int rowL = dir == 'R' ? 1 : 4; + int colL = dir == 'L' ? 1 : 4; + double binomial[5] = {1, 1, 2, 6}; + for (int row = 0; row < 4; row++) { + for (int col = 0; col < 4; col++) { + double acc = 0; + for (int k = col; k < 4; k++) { + acc += dl[rowL * row + colL * k] * pow(wid, -k) * pow(-lo, k - col) + * binomial[k] / binomial[col] / binomial[k - col]; + } + dl[rowL * row + colL * col] = acc; + } + } +} + +void PairAIREBO::Spbicubic_patch_coeffs( + double xmin, double xmax, double ymin, double ymax, double * y, + double * y1, double * y2, double * dl +) { + const double C_inv[16][12] = { + // output_matrix(1, 4*3, get_matrix(2)) + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, + -3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, + 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, + 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, + -3, 0, 3, 0,-2, 0,-1, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, + 9,-9,-9, 9, 6,-6, 3,-3, 6, 3,-6,-3, + -6, 6, 6,-6,-4, 4,-2, 2,-3,-3, 3, 3, + 2, 0,-2, 0, 1, 0, 1, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, + -6, 6, 6,-6,-3, 3,-3, 3,-4,-2, 4, 2, + 4,-4,-4, 4, 2,-2, 2,-2, 2, 2,-2,-2, + }; + double dx = xmax - xmin; + double dy = ymax - ymin; + double x[12]; + for (int i = 0; i < 4; i++) { + x[i+0*4] = y[i]; + x[i+1*4] = y1[i] * dx; + x[i+2*4] = y2[i] * dy; + } + for (int i = 0; i < 16; i++) { + dl[i] = 0; + for (int k = 0; k < 12; k++) { + dl[i] += x[k] * C_inv[i][k]; + } + } + Spbicubic_patch_adjust(dl, dx, xmin, 'R'); + Spbicubic_patch_adjust(dl, dy, ymin, 'L'); +} + /* ---------------------------------------------------------------------- initialize spline knot values ------------------------------------------------------------------------- */ @@ -4107,6 +4377,22 @@ void PairAIREBO::spline_init() PCHf[2][1] = -0.300529172; PCHf[3][0] = -0.307584705; + for (int nH = 0; nH < 4; nH++) { + for (int nC = 0; nC < 4; nC++) { + double y[4] = {0}, y1[4] = {0}, y2[4] = {0}; + y[0] = PCCf[nC][nH]; + y[1] = PCCf[nC][nH+1]; + y[2] = PCCf[nC+1][nH]; + y[3] = PCCf[nC+1][nH+1]; + Spbicubic_patch_coeffs(nC, nC+1, nH, nH+1, y, y1, y2, &pCC[nC][nH][0]); + y[0] = PCHf[nC][nH]; + y[1] = PCHf[nC][nH+1]; + y[2] = PCHf[nC+1][nH]; + y[3] = PCHf[nC+1][nH+1]; + Spbicubic_patch_coeffs(nC, nC+1, nH, nH+1, y, y1, y2, &pCH[nC][nH][0]); + } + } + for (i = 0; i < 5; i++) { for (j = 0; j < 5; j++) { for (k = 0; k < 10; k++) { @@ -4271,6 +4557,46 @@ void PairAIREBO::spline_init() Tf[2][2][1] = -0.035140; for (i = 2; i < 10; i++) Tf[2][2][i] = -0.0040480; + + for (int nH = 0; nH < 4; nH++) { + for (int nC = 0; nC < 4; nC++) { + // Note: Spline knot values exist up to "10", but are never used because + // they are clamped down to 9. + for (int nConj = 0; nConj < 9; nConj++) { + double y[8] = {0}, y1[8] = {0}, y2[8] = {0}, y3[8] = {0}; + #define FILL_KNOTS_TRI(dest, src) \ + dest[0] = src[nC+0][nH+0][nConj+0]; \ + dest[1] = src[nC+0][nH+0][nConj+1]; \ + dest[2] = src[nC+0][nH+1][nConj+0]; \ + dest[3] = src[nC+0][nH+1][nConj+1]; \ + dest[4] = src[nC+1][nH+0][nConj+0]; \ + dest[5] = src[nC+1][nH+0][nConj+1]; \ + dest[6] = src[nC+1][nH+1][nConj+0]; \ + dest[7] = src[nC+1][nH+1][nConj+1]; + FILL_KNOTS_TRI(y, piCCf) + FILL_KNOTS_TRI(y1, piCCdfdx) + FILL_KNOTS_TRI(y2, piCCdfdy) + FILL_KNOTS_TRI(y3, piCCdfdz) + Sptricubic_patch_coeffs(nC, nC+1, nH, nH+1, nConj, nConj+1, y, y1, y2, y3, &piCC[nC][nH][nConj][0]); + FILL_KNOTS_TRI(y, piCHf) + FILL_KNOTS_TRI(y1, piCHdfdx) + FILL_KNOTS_TRI(y2, piCHdfdy) + FILL_KNOTS_TRI(y3, piCHdfdz) + Sptricubic_patch_coeffs(nC, nC+1, nH, nH+1, nConj, nConj+1, y, y1, y2, y3, &piCH[nC][nH][nConj][0]); + FILL_KNOTS_TRI(y, piHHf) + FILL_KNOTS_TRI(y1, piHHdfdx) + FILL_KNOTS_TRI(y2, piHHdfdy) + FILL_KNOTS_TRI(y3, piHHdfdz) + Sptricubic_patch_coeffs(nC, nC+1, nH, nH+1, nConj, nConj+1, y, y1, y2, y3, &piHH[nC][nH][nConj][0]); + FILL_KNOTS_TRI(y, Tf) + FILL_KNOTS_TRI(y1, Tdfdx) + FILL_KNOTS_TRI(y2, Tdfdy) + FILL_KNOTS_TRI(y3, Tdfdz) + Sptricubic_patch_coeffs(nC, nC+1, nH, nH+1, nConj, nConj+1, y, y1, y2, y3, &Tijc[nC][nH][nConj][0]); + #undef FILL_KNOTS_TRI + } + } + } } /* ---------------------------------------------------------------------- diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index e927794ec0..aabc3fe17a 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -113,6 +113,12 @@ class PairAIREBO : public Pair { double Sp5th(double, double *, double *); double Spbicubic(double, double, double *, double *); double Sptricubic(double, double, double, double *, double *); + void Sptricubic_patch_adjust(double *, double, double, char); + void Sptricubic_patch_coeffs(double, double, double, double, double, double, + double*, double*, double*, double*, double*); + void Spbicubic_patch_adjust(double *, double, double, char); + void Spbicubic_patch_coeffs(double, double, double, double, double *, + double *, double *, double *); void spline_init(); void allocate(); From 74d63c24fd349c06594b94b6e51415677b1e98b2 Mon Sep 17 00:00:00 2001 From: Markus Hoehnerbach Date: Wed, 5 Jul 2017 14:51:34 +0200 Subject: [PATCH 2/9] Fix AIREBO missing derivative in bondorderLJ This change replaces the bondorderLJ() function with code provided by Github user CF17, which is based on the bondorder() code. It could be fixed with a shorter patch [1], but layering fix upon fix seems to be unwise in this case. While the code at this point departs from following the Fortran code closely, the reason is that the bug is present in the Fortran code as well. Instead, the new code follows closely the bondorder() code that already exists, which should be easier to maintain in the future. This patch makes the two functions consistent with each other, and makes outside contributions easier. Since it uses a different approach to compute its value, some explanation of that reasoning has been added on top. 1: https://github.com/v0i0/lammps/commit/e8c5c662b28914de6b77e0ca49010eddfe1423e9 --- src/MANYBODY/pair_airebo.cpp | 401 +++++++++++++++++------------------ 1 file changed, 200 insertions(+), 201 deletions(-) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index b8ef4db48c..bc3af2d332 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -2081,44 +2081,63 @@ double PairAIREBO::bondorder(int i, int j, double rij[3], /* ---------------------------------------------------------------------- Bij* function -------------------------------------------------------------------------- */ +------------------------------------------------------------------------- -double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, - double VA, double rij0[3], double rij0mag, +This function calculates S(t_b(b_ij*)) as specified in the AIREBO paper. +To do so, it needs to compute b_ij*, i.e. the bondorder given that the +atoms i and j are placed a ficticious distance rijmag_mod apart. +Now there are two approaches to calculate the resulting forces: +1. Carry through the ficticious distance and corresponding vector + rij_mod, correcting afterwards using the derivative of r/|r|. +2. Perform all the calculations using the real distance, and do not + use a correction, only using rijmag_mod where necessary. +This code opts for (2). Mathematically, the approaches are equivalent +if implemented correctly, since in terms where only the normalized +vector is used, both calculations necessarily lead to the same result +since if f(x) = g(x/|x|) then for x = y/|y| f(x) = g(y/|y|/1). +The actual modified distance is only used in the lamda terms. +Note that these do not contribute to the forces between i and j, since +rijmag_mod is a constant and the corresponding derivatives are +accordingly zero. +This function should be kept in sync with bondorder(), i.e. changes +there probably also need to be performed here. + +*/ + +double PairAIREBO::bondorderLJ(int i, int j, double rij_mod[3], double rijmag_mod, + double VA, double rij[3], double rijmag, double **f, int vflag_atom) { - int k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; - int atomi,atomj,itype,jtype,ktype,ltype,ntype; - double rik[3], rjl[3], rkn[3],rknmag,dNki; + int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; + int itype,jtype,ktype,ltype,ntype; + double rik[3],rjl[3],rkn[3],rji[3],rki[3],rlj[3],rknmag,dNki,dwjl,bij; double NijC,NijH,NjiC,NjiH,wik,dwik,dwkn,wjl; double rikmag,rjlmag,cosjik,cosijl,g,tmp2,tmp3; - double Etmp,pij,tmp,wij,dwij,NconjtmpI,NconjtmpJ; - double Nki,Nlj,dS,lamdajik,lamdaijl,dgdc,dgdN,pji,Nijconj,piRC; - double dN2[2],dN3[3],dwjl; - double Tij,crosskij[3],crosskijmag; - double crossijl[3],crossijlmag,omkijl; - double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3]; - double bij,tmp3pij,tmp3pji,Stb,dStb; - double r32[3],r32mag,cos321; + double Etmp,pij,tmp,wij,dwij,NconjtmpI,NconjtmpJ,Nki,Nlj,dS; + double lamdajik,lamdaijl,dgdc,dgdN,pji,Nijconj,piRC; + double dcosjikdri[3],dcosijldri[3],dcosjikdrk[3]; + double dN2[2],dN3[3]; + double dcosjikdrj[3],dcosijldrj[3],dcosijldrl[3]; + double Tij; + double r32[3],r32mag,cos321,r43[3],r13[3]; + double dNlj; double om1234,rln[3]; double rlnmag,dwln,r23[3],r23mag,r21[3],r21mag; double w21,dw21,r34[3],r34mag,cos234,w34,dw34; double cross321[3],cross234[3],prefactor,SpN; - double fcikpc,fcjlpc,fcjkpc,fcilpc; - double dt2dik[3],dt2djl[3],aa,aaa2,at2,cw,cwnum,cwnom; + double fcikpc,fcjlpc,fcjkpc,fcilpc,fcijpc; + double dt2dik[3],dt2djl[3],dt2dij[3],aa,aaa1,aaa2,at2,cw,cwnum,cwnom; double sin321,sin234,rr,rijrik,rijrjl,rjk2,rik2,ril2,rjl2; - double dctik,dctjk,dctjl,dctil,rik2i,rjl2i,sink2i,sinl2i; - double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil; - double dNlj; - double PijS,PjiS; + double dctik,dctjk,dctjl,dctij,dctji,dctil,rik2i,rjl2i,sink2i,sinl2i; + double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil,dt1dij; + double F23[3],F12[3],F34[3],F31[3],F24[3],fi[3],fj[3],fk[3],fl[3]; + double f1[3],f2[3],f3[3],f4[4]; + double dcut321,PijS,PjiS; double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp; int *REBO_neighs,*REBO_neighs_i,*REBO_neighs_j,*REBO_neighs_k,*REBO_neighs_l; - double F12[3],F34[3],F31[3],F24[3]; - double fi[3],fj[3],fk[3],fl[3],f1[3],f2[3],f3[3],f4[4]; - double rji[3],rki[3],rlj[3],r13[3],r43[3]; - double realrij[3], realrijmag; - double rjkmag, rilmag, dctdjk, dctdik, dctdil, dctdjl; - double fjk[3], fil[3], rijmbr; + double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3]; + double tmp3pij,tmp3pji,Stb,dStb; + double **x = atom->x; int *type = atom->type; @@ -2127,12 +2146,11 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, atomj = j; itype = map[type[atomi]]; jtype = map[type[atomj]]; - wij = Sp(rij0mag,rcmin[itype][jtype],rcmax[itype][jtype],dwij); + wij = Sp(rijmag,rcmin[itype][jtype],rcmax[itype][jtype],dwij); NijC = nC[atomi]-(wij*kronecker(jtype,0)); NijH = nH[atomi]-(wij*kronecker(jtype,1)); NjiC = nC[atomj]-(wij*kronecker(itype,0)); NjiH = nH[atomj]-(wij*kronecker(itype,1)); - bij = 0.0; tmp = 0.0; tmp2 = 0.0; @@ -2145,12 +2163,6 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, Stb = 0.0; dStb = 0.0; - realrij[0] = x[atomi][0] - x[atomj][0]; - realrij[1] = x[atomi][1] - x[atomj][1]; - realrij[2] = x[atomi][2] - x[atomj][2]; - realrijmag = sqrt(realrij[0] * realrij[0] + realrij[1] * realrij[1] - + realrij[2] * realrij[2]); - REBO_neighs = REBO_firstneigh[i]; for (k = 0; k < REBO_numneigh[i]; k++) { atomk = REBO_neighs[k]; @@ -2161,9 +2173,9 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, rik[2] = x[atomi][2]-x[atomk][2]; rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2])); lamdajik = 4.0*kronecker(itype,1) * - ((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag)); + ((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag_mod)); wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dS); - Nki = nC[atomk]-(wik*kronecker(itype,0)) + + Nki = nC[atomk]-(wik*kronecker(itype,0))+ nH[atomk]-(wik*kronecker(itype,1)); cosjik = ((rij[0]*rik[0])+(rij[1]*rik[1])+(rij[2]*rik[2])) / (rijmag*rikmag); @@ -2186,6 +2198,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, pij = 1.0/sqrt(1.0+Etmp+PijS); tmppij = -.5*cube(pij); tmp3pij = tmp3; + tmp = 0.0; tmp2 = 0.0; tmp3 = 0.0; @@ -2201,7 +2214,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, rjl[2] = x[atomj][2]-x[atoml][2]; rjlmag = sqrt((rjl[0]*rjl[0])+(rjl[1]*rjl[1])+(rjl[2]*rjl[2])); lamdaijl = 4.0*kronecker(jtype,1) * - ((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag)); + ((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag_mod)); wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dS); Nlj = nC[atoml]-(wjl*kronecker(jtype,0))+nH[atoml] - (wjl*kronecker(jtype,1)); @@ -2231,82 +2244,80 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, Nijconj = 1.0+(NconjtmpI*NconjtmpI)+(NconjtmpJ*NconjtmpJ); piRC = piRCSpline(NijC+NijH,NjiC+NjiH,Nijconj,itype,jtype,dN3piRC); + Tij = 0.0; dN3Tij[0] = 0.0; dN3Tij[1] = 0.0; dN3Tij[2] = 0.0; if (itype == 0 && jtype == 0) Tij=TijSpline((NijC+NijH),(NjiC+NjiH),Nijconj,dN3Tij); - Etmp = 0.0; + if (fabs(Tij) > TOL) { + atom2 = atomi; + atom3 = atomj; + r32[0] = x[atom3][0]-x[atom2][0]; + r32[1] = x[atom3][1]-x[atom2][1]; + r32[2] = x[atom3][2]-x[atom2][2]; + r32mag = sqrt((r32[0]*r32[0])+(r32[1]*r32[1])+(r32[2]*r32[2])); + r23[0] = -r32[0]; + r23[1] = -r32[1]; + r23[2] = -r32[2]; + r23mag = r32mag; REBO_neighs_i = REBO_firstneigh[i]; for (k = 0; k < REBO_numneigh[i]; k++) { atomk = REBO_neighs_i[k]; + atom1 = atomk; ktype = map[type[atomk]]; if (atomk != atomj) { - rik[0] = x[atomi][0]-x[atomk][0]; - rik[1] = x[atomi][1]-x[atomk][1]; - rik[2] = x[atomi][2]-x[atomk][2]; - rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2])); - cos321 = ((rij[0]*rik[0])+(rij[1]*rik[1])+(rij[2]*rik[2])) / - (rijmag*rikmag); + r21[0] = x[atom2][0]-x[atom1][0]; + r21[1] = x[atom2][1]-x[atom1][1]; + r21[2] = x[atom2][2]-x[atom1][2]; + r21mag = sqrt(r21[0]*r21[0] + r21[1]*r21[1] + r21[2]*r21[2]); + cos321 = -1.0*((r21[0]*r32[0])+(r21[1]*r32[1])+(r21[2]*r32[2])) / + (r21mag*r32mag); cos321 = MIN(cos321,1.0); cos321 = MAX(cos321,-1.0); + sin321 = sqrt(1.0 - cos321*cos321); + if ((sin321 > TOL) && (r21mag > TOL)) { // XXX was sin321 != 0.0 + w21 = Sp(r21mag,rcmin[itype][ktype],rcmaxp[itype][ktype],dw21); + tspjik = Sp2(cos321,thmin,thmax,dtsjik); - rjk[0] = rik[0]-rij[0]; - rjk[1] = rik[1]-rij[1]; - rjk[2] = rik[2]-rij[2]; - rjk2 = (rjk[0]*rjk[0])+(rjk[1]*rjk[1])+(rjk[2]*rjk[2]); - rij2 = rijmag*rijmag; - rik2 = rikmag*rikmag; - costmp = 0.5*(rij2+rik2-rjk2)/rijmag/rikmag; - tspjik = Sp2(costmp,thmin,thmax,dtsjik); - - if (sqrt(1.0 - cos321*cos321) > sqrt(TOL)) { - wik = Sp(rikmag,rcmin[itype][ktype],rcmaxp[itype][ktype],dwik); REBO_neighs_j = REBO_firstneigh[j]; for (l = 0; l < REBO_numneigh[j]; l++) { atoml = REBO_neighs_j[l]; + atom4 = atoml; ltype = map[type[atoml]]; if (!(atoml == atomi || atoml == atomk)) { - rjl[0] = x[atomj][0]-x[atoml][0]; - rjl[1] = x[atomj][1]-x[atoml][1]; - rjl[2] = x[atomj][2]-x[atoml][2]; - rjlmag = sqrt(rjl[0]*rjl[0] + rjl[1]*rjl[1] + rjl[2]*rjl[2]); - cos234 = -((rij[0]*rjl[0])+(rij[1]*rjl[1])+(rij[2]*rjl[2])) / - (rijmag*rjlmag); + r34[0] = x[atom3][0]-x[atom4][0]; + r34[1] = x[atom3][1]-x[atom4][1]; + r34[2] = x[atom3][2]-x[atom4][2]; + r34mag = sqrt((r34[0]*r34[0])+(r34[1]*r34[1])+(r34[2]*r34[2])); + cos234 = (r32[0]*r34[0] + r32[1]*r34[1] + r32[2]*r34[2]) / + (r32mag*r34mag); cos234 = MIN(cos234,1.0); cos234 = MAX(cos234,-1.0); + sin234 = sqrt(1.0 - cos234*cos234); - ril[0] = rij[0]+rjl[0]; - ril[1] = rij[1]+rjl[1]; - ril[2] = rij[2]+rjl[2]; - ril2 = (ril[0]*ril[0])+(ril[1]*ril[1])+(ril[2]*ril[2]); - rjl2 = rjlmag*rjlmag; - costmp = 0.5*(rij2+rjl2-ril2)/rijmag/rjlmag; - tspijl = Sp2(costmp,thmin,thmax,dtsijl); + if ((sin234 > TOL) && (r34mag > TOL)) { // XXX was sin234 != 0.0 + w34 = Sp(r34mag,rcmin[jtype][ltype],rcmaxp[jtype][ltype],dw34); + tspijl = Sp2(cos234,thmin,thmax,dtsijl); - if (sqrt(1.0 - cos234*cos234) > sqrt(TOL)) { - wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmaxp[jtype][ltype],dS); - crosskij[0] = (rij[1]*rik[2]-rij[2]*rik[1]); - crosskij[1] = (rij[2]*rik[0]-rij[0]*rik[2]); - crosskij[2] = (rij[0]*rik[1]-rij[1]*rik[0]); - crosskijmag = sqrt(crosskij[0]*crosskij[0] + - crosskij[1]*crosskij[1] + - crosskij[2]*crosskij[2]); - crossijl[0] = (rij[1]*rjl[2]-rij[2]*rjl[1]); - crossijl[1] = (rij[2]*rjl[0]-rij[0]*rjl[2]); - crossijl[2] = (rij[0]*rjl[1]-rij[1]*rjl[0]); - crossijlmag = sqrt(crossijl[0]*crossijl[0] + - crossijl[1]*crossijl[1] + - crossijl[2]*crossijl[2]); - omkijl = -1.0*(((crosskij[0]*crossijl[0]) + - (crosskij[1]*crossijl[1]) + - (crosskij[2]*crossijl[2])) / - (crosskijmag*crossijlmag)); - Etmp += ((1.0-square(omkijl))*wik*wjl) * + cross321[0] = (r32[1]*r21[2])-(r32[2]*r21[1]); + cross321[1] = (r32[2]*r21[0])-(r32[0]*r21[2]); + cross321[2] = (r32[0]*r21[1])-(r32[1]*r21[0]); + cross234[0] = (r23[1]*r34[2])-(r23[2]*r34[1]); + cross234[1] = (r23[2]*r34[0])-(r23[0]*r34[2]); + cross234[2] = (r23[0]*r34[1])-(r23[1]*r34[0]); + + cwnum = (cross321[0]*cross234[0]) + + (cross321[1]*cross234[1]) + (cross321[2]*cross234[2]); + cwnom = r21mag*r34mag*r23mag*r23mag*sin321*sin234; + om1234 = cwnum/cwnom; + cw = om1234; + Etmp += ((1.0-square(om1234))*w21*w34) * (1.0-tspjik)*(1.0-tspijl); + } } } @@ -2338,50 +2349,43 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, rik[2] = x[atomi][2]-x[atomk][2]; rikmag = sqrt(rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2]); lamdajik = 4.0*kronecker(itype,1) * - ((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag)); + ((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag_mod)); wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik); - rjk[0] = rik[0] - rij[0]; - rjk[1] = rik[1] - rij[1]; - rjk[2] = rik[2] - rij[2]; - rjkmag = sqrt(rjk[0] * rjk[0] + rjk[1] * rjk[1] + rjk[2] * rjk[2]); - rijrik = 2 * rijmag * rikmag; - rr = rijmag * rijmag - rikmag * rikmag; - cosjik = (rijmag * rijmag + rikmag * rikmag - rjkmag * rjkmag) / rijrik; + cosjik = (rij[0]*rik[0] + rij[1]*rik[1] + rij[2]*rik[2]) / + (rijmag*rikmag); cosjik = MIN(cosjik,1.0); cosjik = MAX(cosjik,-1.0); - dctdjk = -2 / rijrik; - dctdik = (-rr + rjkmag * rjkmag) / (rijrik * rikmag * rikmag); - // evaluate splines g and derivatives dg + + dcosjikdri[0] = ((rij[0]+rik[0])/(rijmag*rikmag)) - + (cosjik*((rij[0]/(rijmag*rijmag))+(rik[0]/(rikmag*rikmag)))); + dcosjikdri[1] = ((rij[1]+rik[1])/(rijmag*rikmag)) - + (cosjik*((rij[1]/(rijmag*rijmag))+(rik[1]/(rikmag*rikmag)))); + dcosjikdri[2] = ((rij[2]+rik[2])/(rijmag*rikmag)) - + (cosjik*((rij[2]/(rijmag*rijmag))+(rik[2]/(rikmag*rikmag)))); + dcosjikdrk[0] = (-rij[0]/(rijmag*rikmag)) + + (cosjik*(rik[0]/(rikmag*rikmag))); + dcosjikdrk[1] = (-rij[1]/(rijmag*rikmag)) + + (cosjik*(rik[1]/(rikmag*rikmag))); + dcosjikdrk[2] = (-rij[2]/(rijmag*rikmag)) + + (cosjik*(rik[2]/(rikmag*rikmag))); + dcosjikdrj[0] = (-rik[0]/(rijmag*rikmag)) + + (cosjik*(rij[0]/(rijmag*rijmag))); + dcosjikdrj[1] = (-rik[1]/(rijmag*rikmag)) + + (cosjik*(rij[1]/(rijmag*rijmag))); + dcosjikdrj[2] = (-rik[2]/(rijmag*rikmag)) + + (cosjik*(rij[2]/(rijmag*rijmag))); g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN); - tmp2 = VA*.5*(tmp*wik*dgdc*exp(lamdajik)); - - fi[0] = -tmp2 * dctdik * rik[0]; - fi[1] = -tmp2 * dctdik * rik[1]; - fi[2] = -tmp2 * dctdik * rik[2]; - fk[0] = tmp2 * dctdik * rik[0]; - fk[1] = tmp2 * dctdik * rik[1]; - fk[2] = tmp2 * dctdik * rik[2]; - fj[0] = 0; - fj[1] = 0; - fj[2] = 0; - fjk[0] = -tmp2 * dctdjk * rjk[0]; - fjk[1] = -tmp2 * dctdjk * rjk[1]; - fjk[2] = -tmp2 * dctdjk * rjk[2]; - fi[0] += fjk[0]; - fi[1] += fjk[1]; - fi[2] += fjk[2]; - fk[0] -= fjk[0]; - fk[1] -= fjk[1]; - fk[2] -= fjk[2]; - rijmbr = rcmin[itype][jtype] / realrijmag; - fj[0] += rijmbr * (fjk[0] - (realrij[0] * realrij[0] * fjk[0] + realrij[0] * realrij[1] * fjk[1] + realrij[0] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fj[1] += rijmbr * (fjk[1] - (realrij[1] * realrij[0] * fjk[0] + realrij[1] * realrij[1] * fjk[1] + realrij[1] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fj[2] += rijmbr * (fjk[2] - (realrij[2] * realrij[0] * fjk[0] + realrij[2] * realrij[1] * fjk[1] + realrij[2] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fi[0] -= rijmbr * (fjk[0] - (realrij[0] * realrij[0] * fjk[0] + realrij[0] * realrij[1] * fjk[1] + realrij[0] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fi[1] -= rijmbr * (fjk[1] - (realrij[1] * realrij[0] * fjk[0] + realrij[1] * realrij[1] * fjk[1] + realrij[1] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fi[2] -= rijmbr * (fjk[2] - (realrij[2] * realrij[0] * fjk[0] + realrij[2] * realrij[1] * fjk[1] + realrij[2] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); + fj[0] = -tmp2*dcosjikdrj[0]; + fj[1] = -tmp2*dcosjikdrj[1]; + fj[2] = -tmp2*dcosjikdrj[2]; + fi[0] = -tmp2*dcosjikdri[0]; + fi[1] = -tmp2*dcosjikdri[1]; + fi[2] = -tmp2*dcosjikdri[2]; + fk[0] = -tmp2*dcosjikdrk[0]; + fk[1] = -tmp2*dcosjikdrk[1]; + fk[2] = -tmp2*dcosjikdrk[2]; tmp2 = VA*.5*(tmp*wik*g*exp(lamdajik)*4.0*kronecker(itype,1)); fi[0] += tmp2*(rik[0]/rikmag); @@ -2439,6 +2443,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, tmp3 = tmp3pji; dN2[0] = dN2PJI[0]; dN2[1] = dN2PJI[1]; + REBO_neighs = REBO_firstneigh[j]; for (l = 0; l < REBO_numneigh[j]; l++) { atoml = REBO_neighs[l]; @@ -2449,51 +2454,43 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, rjl[2] = x[atomj][2]-x[atoml][2]; rjlmag = sqrt((rjl[0]*rjl[0])+(rjl[1]*rjl[1])+(rjl[2]*rjl[2])); lamdaijl = 4.0*kronecker(jtype,1) * - ((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag)); + ((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag_mod)); wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl); - ril[0] = rij[0] + rjl[0]; - ril[1] = rij[1] + rjl[1]; - ril[2] = rij[2] + rjl[2]; - rilmag = sqrt(ril[0] * ril[0] + ril[1] * ril[1] + ril[2] * ril[2]); - rijrjl = 2 * rijmag * rjlmag; - rr = rijmag * rijmag - rjlmag * rjlmag; - cosijl = (rijmag * rijmag + rjlmag * rjlmag - rilmag * rilmag) / rijrjl; + cosijl = (-1.0*((rij[0]*rjl[0])+(rij[1]*rjl[1])+(rij[2]*rjl[2]))) / + (rijmag*rjlmag); cosijl = MIN(cosijl,1.0); cosijl = MAX(cosijl,-1.0); - dctdil = -2 / rijrjl; - dctdjl = (-rr + rilmag * rilmag) / (rijrjl * rjlmag * rjlmag); + + dcosijldri[0] = (-rjl[0]/(rijmag*rjlmag)) - + (cosijl*rij[0]/(rijmag*rijmag)); + dcosijldri[1] = (-rjl[1]/(rijmag*rjlmag)) - + (cosijl*rij[1]/(rijmag*rijmag)); + dcosijldri[2] = (-rjl[2]/(rijmag*rjlmag)) - + (cosijl*rij[2]/(rijmag*rijmag)); + dcosijldrj[0] = ((-rij[0]+rjl[0])/(rijmag*rjlmag)) + + (cosijl*((rij[0]/square(rijmag))-(rjl[0]/(rjlmag*rjlmag)))); + dcosijldrj[1] = ((-rij[1]+rjl[1])/(rijmag*rjlmag)) + + (cosijl*((rij[1]/square(rijmag))-(rjl[1]/(rjlmag*rjlmag)))); + dcosijldrj[2] = ((-rij[2]+rjl[2])/(rijmag*rjlmag)) + + (cosijl*((rij[2]/square(rijmag))-(rjl[2]/(rjlmag*rjlmag)))); + dcosijldrl[0] = (rij[0]/(rijmag*rjlmag))+(cosijl*rjl[0]/(rjlmag*rjlmag)); + dcosijldrl[1] = (rij[1]/(rijmag*rjlmag))+(cosijl*rjl[1]/(rjlmag*rjlmag)); + dcosijldrl[2] = (rij[2]/(rijmag*rjlmag))+(cosijl*rjl[2]/(rjlmag*rjlmag)); // evaluate splines g and derivatives dg g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wjl*dgdc*exp(lamdaijl)); - fj[0] = -tmp2 * dctdjl * rjl[0]; - fj[1] = -tmp2 * dctdjl * rjl[1]; - fj[2] = -tmp2 * dctdjl * rjl[2]; - fl[0] = tmp2 * dctdjl * rjl[0]; - fl[1] = tmp2 * dctdjl * rjl[1]; - fl[2] = tmp2 * dctdjl * rjl[2]; - fi[0] = 0; - fi[1] = 0; - fi[2] = 0; - fil[0] = -tmp2 * dctdil * ril[0]; - fil[1] = -tmp2 * dctdil * ril[1]; - fil[2] = -tmp2 * dctdil * ril[2]; - fj[0] += fil[0]; - fj[1] += fil[1]; - fj[2] += fil[2]; - fl[0] -= fil[0]; - fl[1] -= fil[1]; - fl[2] -= fil[2]; - rijmbr = rcmin[itype][jtype] / realrijmag; - fi[0] += rijmbr * (fil[0] - (realrij[0] * realrij[0] * fil[0] + realrij[0] * realrij[1] * fil[1] + realrij[0] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fi[1] += rijmbr * (fil[1] - (realrij[1] * realrij[0] * fil[0] + realrij[1] * realrij[1] * fil[1] + realrij[1] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fi[2] += rijmbr * (fil[2] - (realrij[2] * realrij[0] * fil[0] + realrij[2] * realrij[1] * fil[1] + realrij[2] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fj[0] -= rijmbr * (fil[0] - (realrij[0] * realrij[0] * fil[0] + realrij[0] * realrij[1] * fil[1] + realrij[0] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fj[1] -= rijmbr * (fil[1] - (realrij[1] * realrij[0] * fil[0] + realrij[1] * realrij[1] * fil[1] + realrij[1] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fj[2] -= rijmbr * (fil[2] - (realrij[2] * realrij[0] * fil[0] + realrij[2] * realrij[1] * fil[1] + realrij[2] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - - + fi[0] = -tmp2*dcosijldri[0]; + fi[1] = -tmp2*dcosijldri[1]; + fi[2] = -tmp2*dcosijldri[2]; + fj[0] = -tmp2*dcosijldrj[0]; + fj[1] = -tmp2*dcosijldrj[1]; + fj[2] = -tmp2*dcosijldrj[2]; + fl[0] = -tmp2*dcosijldrl[0]; + fl[1] = -tmp2*dcosijldrl[1]; + fl[2] = -tmp2*dcosijldrl[2]; + tmp2 = VA*.5*(tmp*wjl*g*exp(lamdaijl)*4.0*kronecker(jtype,1)); fj[0] += tmp2*(rjl[0]/rjlmag); fj[1] += tmp2*(rjl[1]/rjlmag); @@ -2503,6 +2500,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, fl[2] -= tmp2*(rjl[2]/rjlmag); // coordination forces + // dwik forces tmp2 = VA*.5*(tmp*dwjl*g*exp(lamdaijl))/rjlmag; @@ -2617,8 +2615,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, // piRC forces to J side - REBO_neighs = REBO_firstneigh[j]; - for (l = 0; l < REBO_numneigh[j]; l++) { + REBO_neighs = REBO_firstneigh[atomj]; + for (l = 0; l < REBO_numneigh[atomj]; l++) { atoml = REBO_neighs[l]; if (atoml != atomi) { ltype = map[type[atoml]]; @@ -2688,15 +2686,14 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, dN3[2] = dN3Tij[2]; atom2 = atomi; atom3 = atomj; - r32[0] = -rij[0]; - r32[1] = -rij[1]; - r32[2] = -rij[2]; - r23[0] = rij[0]; - r23[1] = rij[1]; - r23[2] = rij[2]; - r32mag = rijmag; - r23mag = rijmag; - + r32[0] = x[atom3][0]-x[atom2][0]; + r32[1] = x[atom3][1]-x[atom2][1]; + r32[2] = x[atom3][2]-x[atom2][2]; + r32mag = sqrt((r32[0]*r32[0])+(r32[1]*r32[1])+(r32[2]*r32[2])); + r23[0] = -r32[0]; + r23[1] = -r32[1]; + r23[2] = -r32[2]; + r23mag = r32mag; REBO_neighs_i = REBO_firstneigh[i]; for (k = 0; k < REBO_numneigh[i]; k++) { atomk = REBO_neighs_i[k]; @@ -2712,20 +2709,21 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, cos321 = MIN(cos321,1.0); cos321 = MAX(cos321,-1.0); sin321 = sqrt(1.0 - cos321*cos321); - if ((sin321 > TOL) && (r21mag > TOL)) { // XXX was sin321 != 0.0 sink2i = 1.0/(sin321*sin321); rik2i = 1.0/(r21mag*r21mag); rr = (rijmag*rijmag)-(r21mag*r21mag); - rjk[0] = r21[0]-rij[0]; - rjk[1] = r21[1]-rij[1]; - rjk[2] = r21[2]-rij[2]; + rjk[0] = r21[0]-r23[0]; + rjk[1] = r21[1]-r23[1]; + rjk[2] = r21[2]-r23[2]; rjk2 = (rjk[0]*rjk[0])+(rjk[1]*rjk[1])+(rjk[2]*rjk[2]); - rijrik = 2.0*rijmag*r21mag; + rijrik = 2.0*r23mag*r21mag; rik2 = r21mag*r21mag; dctik = (-rr+rjk2)/(rijrik*rik2); + dctij = (rr+rjk2)/(rijrik*r23mag*r23mag); dctjk = -2.0/rijrik; w21 = Sp(r21mag,rcmin[itype][ktype],rcmaxp[itype][ktype],dw21); + rijmag = r32mag; rikmag = r21mag; rij2 = r32mag*r32mag; rik2 = r21mag*r21mag; @@ -2743,8 +2741,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, r34[1] = x[atom3][1]-x[atom4][1]; r34[2] = x[atom3][2]-x[atom4][2]; r34mag = sqrt(r34[0]*r34[0] + r34[1]*r34[1] + r34[2]*r34[2]); - cos234 = -1.0*((rij[0]*r34[0])+(rij[1]*r34[1]) + - (rij[2]*r34[2]))/(rijmag*r34mag); + cos234 = (r32[0]*r34[0] + r32[1]*r34[1] + r32[2]*r34[2]) / + (r32mag*r34mag); cos234 = MIN(cos234,1.0); cos234 = MAX(cos234,-1.0); sin234 = sqrt(1.0 - cos234*cos234); @@ -2762,6 +2760,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, rijrjl = 2.0*r23mag*r34mag; rjl2 = r34mag*r34mag; dctjl = (-rr+ril2)/(rijrjl*rjl2); + dctji = (rr+ril2)/(rijrjl*r23mag*r23mag); dctil = -2.0/rijrjl; rjlmag = r34mag; rjl2 = r34mag*r34mag; @@ -2787,6 +2786,8 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, dt1djk = (-dctjk*sink2i*cos321); dt1djl = (rjl2i)-(dctjl*sinl2i*cos234); dt1dil = (-dctil*sinl2i*cos234); + dt1dij = (2.0/(r23mag*r23mag))-(dctij*sink2i*cos321) - + (dctji*sinl2i*cos234); dt2dik[0] = (-r23[2]*cross234[1])+(r23[1]*cross234[2]); dt2dik[1] = (-r23[0]*cross234[2])+(r23[2]*cross234[0]); @@ -2796,16 +2797,29 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, dt2djl[1] = (-r23[2]*cross321[0])+(r23[0]*cross321[2]); dt2djl[2] = (-r23[0]*cross321[1])+(r23[1]*cross321[0]); + dt2dij[0] = (r21[2]*cross234[1])-(r34[2]*cross321[1]) - + (r21[1]*cross234[2])+(r34[1]*cross321[2]); + dt2dij[1] = (r21[0]*cross234[2])-(r34[0]*cross321[2]) - + (r21[2]*cross234[0])+(r34[2]*cross321[0]); + dt2dij[2] = (r21[1]*cross234[0])-(r34[1]*cross321[0]) - + (r21[0]*cross234[1])+(r34[0]*cross321[1]); + aa = (prefactor*2.0*cw/cwnom)*w21*w34 * (1.0-tspjik)*(1.0-tspijl); aaa2 = -prefactor*(1.0-square(om1234)) * w21*w34; at2 = aa*cwnum; + fcijpc = (-dt1dij*at2)+(aaa2*dtsjik*dctij*(1.0-tspijl)) + + (aaa2*dtsijl*dctji*(1.0-tspjik)); fcikpc = (-dt1dik*at2)+(aaa2*dtsjik*dctik*(1.0-tspijl)); fcjlpc = (-dt1djl*at2)+(aaa2*dtsijl*dctjl*(1.0-tspjik)); fcjkpc = (-dt1djk*at2)+(aaa2*dtsjik*dctjk*(1.0-tspijl)); fcilpc = (-dt1dil*at2)+(aaa2*dtsijl*dctil*(1.0-tspjik)); + F23[0] = (fcijpc*r23[0])+(aa*dt2dij[0]); + F23[1] = (fcijpc*r23[1])+(aa*dt2dij[1]); + F23[2] = (fcijpc*r23[2])+(aa*dt2dij[2]); + F12[0] = (fcikpc*r21[0])+(aa*dt2dik[0]); F12[1] = (fcikpc*r21[1])+(aa*dt2dik[1]); F12[2] = (fcikpc*r21[2])+(aa*dt2dik[2]); @@ -2825,31 +2839,16 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij[3], double rijmag, f1[0] = -F12[0]-F31[0]; f1[1] = -F12[1]-F31[1]; f1[2] = -F12[2]-F31[2]; - f2[0] = F12[0]+F31[0]; - f2[1] = F12[1]+F31[1]; - f2[2] = F12[2]+F31[2]; - f3[0] = F34[0]+F24[0]; - f3[1] = F34[1]+F24[1]; - f3[2] = F34[2]+F24[2]; + f2[0] = F23[0]+F12[0]+F24[0]; + f2[1] = F23[1]+F12[1]+F24[1]; + f2[2] = F23[2]+F12[2]+F24[2]; + f3[0] = -F23[0]+F34[0]+F31[0]; + f3[1] = -F23[1]+F34[1]+F31[1]; + f3[2] = -F23[2]+F34[2]+F31[2]; f4[0] = -F34[0]-F24[0]; f4[1] = -F34[1]-F24[1]; f4[2] = -F34[2]-F24[2]; - rijmbr = rcmin[itype][jtype] / realrijmag; - f2[0] += rijmbr * (F24[0] - (realrij[0] * realrij[0] * F24[0] + realrij[0] * realrij[1] * F24[1] + realrij[0] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f2[1] += rijmbr * (F24[1] - (realrij[1] * realrij[0] * F24[0] + realrij[1] * realrij[1] * F24[1] + realrij[1] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f2[2] += rijmbr * (F24[2] - (realrij[2] * realrij[0] * F24[0] + realrij[2] * realrij[1] * F24[1] + realrij[2] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f3[0] -= rijmbr * (F24[0] - (realrij[0] * realrij[0] * F24[0] + realrij[0] * realrij[1] * F24[1] + realrij[0] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f3[1] -= rijmbr * (F24[1] - (realrij[1] * realrij[0] * F24[0] + realrij[1] * realrij[1] * F24[1] + realrij[1] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f3[2] -= rijmbr * (F24[2] - (realrij[2] * realrij[0] * F24[0] + realrij[2] * realrij[1] * F24[1] + realrij[2] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - - f2[0] -= rijmbr * (F31[0] - (realrij[0] * realrij[0] * F31[0] + realrij[0] * realrij[1] * F31[1] + realrij[0] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f2[1] -= rijmbr * (F31[1] - (realrij[1] * realrij[0] * F31[0] + realrij[1] * realrij[1] * F31[1] + realrij[1] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f2[2] -= rijmbr * (F31[2] - (realrij[2] * realrij[0] * F31[0] + realrij[2] * realrij[1] * F31[1] + realrij[2] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f3[0] += rijmbr * (F31[0] - (realrij[0] * realrij[0] * F31[0] + realrij[0] * realrij[1] * F31[1] + realrij[0] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f3[1] += rijmbr * (F31[1] - (realrij[1] * realrij[0] * F31[0] + realrij[1] * realrij[1] * F31[1] + realrij[1] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f3[2] += rijmbr * (F31[2] - (realrij[2] * realrij[0] * F31[0] + realrij[2] * realrij[1] * F31[1] + realrij[2] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - // coordination forces tmp2 = VA*Tij*((1.0-(om1234*om1234))) * From 7e42af18bc6870b4d8769528228842b27388c97e Mon Sep 17 00:00:00 2001 From: Markus Hoehnerbach Date: Wed, 5 Jul 2017 15:11:58 +0200 Subject: [PATCH 3/9] Feature: AIREBO parametrize cutoff switching In #514 it has been raised that the switching function that ensures a smooth transition to the cutoff is only correct if cutlj = 3.0. This patch gives users an opportunity to configure the switching function together with the cutoff by specifying the start of the transition region. Behaviour in the default case remaing unchanged. This allows users to specify larger cutoffs than 3 (which used to have no effect) and get correct cutoff behaviour for values less then 3. --- doc/src/pair_airebo.txt | 13 ++++++++++--- src/MANYBODY/pair_airebo.cpp | 18 +++++++++++++----- src/MANYBODY/pair_airebo.h | 1 + 3 files changed, 24 insertions(+), 8 deletions(-) diff --git a/doc/src/pair_airebo.txt b/doc/src/pair_airebo.txt index 2e3083c34b..0c03eb3267 100644 --- a/doc/src/pair_airebo.txt +++ b/doc/src/pair_airebo.txt @@ -15,12 +15,13 @@ pair_style rebo/omp command :h3 [Syntax:] -pair_style style cutoff LJ_flag TORSION_flag :pre +pair_style style cutoff LJ_flag TORSION_flag cutoff_min :pre style = {airebo} or {airebo/morse} or {rebo} cutoff = LJ or Morse cutoff (sigma scale factor) (AIREBO and AIREBO-M only) LJ_flag = 0/1 to turn off/on the LJ or Morse term (AIREBO and AIREBO-M only, optional) -TORSION_flag = 0/1 to turn off/on the torsion term (AIREBO and AIREBO-M only, optional) :ul +TORSION_flag = 0/1 to turn off/on the torsion term (AIREBO and AIREBO-M only, optional) +cutoff_min = Start of the transition region of cutoff (sigma scale factor) (AIREBO and AIREBO-M only, optional) :ul [Examples:] @@ -60,7 +61,7 @@ The AIREBO potential consists of three terms: :c,image(Eqs/pair_airebo.jpg) By default, all three terms are included. For the {airebo} style, if -the two optional flag arguments to the pair_style command are +the first two optional flag arguments to the pair_style command are included, the LJ and torsional terms can be turned off. Note that both or neither of the flags must be included. If both of the LJ an torsional terms are turned off, it becomes the 2nd-generation REBO @@ -97,6 +98,12 @@ standard AIREBO potential, sigma_CC = 3.4 Angstroms, so with a scale factor of 3.0 (the argument in pair_style), the resulting E_LJ cutoff would be 10.2 Angstroms. +By default, the longer-ranged interaction is smoothly switched off +between 2.16 and 3.0 sigma. By specifying {cutoff_min} in addition +to {cutoff}, the switching can be configured to take place between +{cutoff_min} and {cutoff}. {cutoff_min} can only be specified if all +optional arguments are given. + The E_TORSION term is an explicit 4-body potential that describes various dihedral angle preferences in hydrocarbon configurations. diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index bc3af2d332..e0cc15fae1 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -68,6 +68,10 @@ PairAIREBO::PairAIREBO(LAMMPS *lmp) : Pair(lmp) nC = nH = NULL; map = NULL; manybody_flag = 1; + + sigwid = 0.84; + sigcut = 3.0; + sigmin = sigcut - sigwid; } /* ---------------------------------------------------------------------- @@ -147,7 +151,8 @@ void PairAIREBO::allocate() void PairAIREBO::settings(int narg, char **arg) { - if (narg != 1 && narg != 3) error->all(FLERR,"Illegal pair_style command"); + if (narg != 1 && narg != 3 && narg != 4) + error->all(FLERR,"Illegal pair_style command"); cutlj = force->numeric(FLERR,arg[0]); @@ -155,6 +160,13 @@ void PairAIREBO::settings(int narg, char **arg) ljflag = force->inumeric(FLERR,arg[1]); torflag = force->inumeric(FLERR,arg[2]); } + if (narg == 4) { + ljflag = force->inumeric(FLERR,arg[1]); + torflag = force->inumeric(FLERR,arg[2]); + sigcut = cutlj; + sigmin = force->numeric(FLERR,arg[3]); + sigwid = sigcut - sigmin; + } // this one parameter for C-C interactions is different in AIREBO vs REBO // see Favata, Micheletti, Ryu, Pugno, Comp Phys Comm (2016) @@ -736,10 +748,6 @@ void PairAIREBO::FLJ(int eflag, int vflag) // compute LJ forces and energy - sigwid = 0.84; - sigcut = 3.0; - sigmin = sigcut - sigwid; - rljmin = sigma[itype][jtype]; rljmax = sigcut * rljmin; rljmin = sigmin * rljmin; diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index aabc3fe17a..d9628f3bb3 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -46,6 +46,7 @@ class PairAIREBO : public Pair { int morseflag; // 1 if Morse instead of LJ for non-bonded double cutlj; // user-specified LJ cutoff + double sigcut,sigwid,sigmin; // corresponding cutoff function double cutljrebosq; // cut for when to compute // REBO neighs of ghost atoms From d2f7f4843a7cc64cb0468f64507d40f203ea0dab Mon Sep 17 00:00:00 2001 From: Markus Hoehnerbach Date: Wed, 5 Jul 2017 15:16:45 +0200 Subject: [PATCH 4/9] AIREBO Fix Credits --- src/MANYBODY/pair_airebo.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index e0cc15fae1..27a5aad21c 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -14,7 +14,8 @@ /* ---------------------------------------------------------------------- Contributing author: Ase Henry (MIT) Bugfixes and optimizations: - Marcel Fallet & Steve Stuart (Clemson), Axel Kohlmeyer (Temple U) + Marcel Fallet & Steve Stuart (Clemson), Axel Kohlmeyer (Temple U), + Markus Hoehnerbach (RWTH Aachen), Github user CF17 AIREBO-M modification to optionally replace LJ with Morse potentials. Thomas C. O'Connor (JHU) 2014 ------------------------------------------------------------------------- */ From ff761d639a78bf40627a26de52724e06fdc375df Mon Sep 17 00:00:00 2001 From: Markus Hoehnerbach Date: Wed, 5 Jul 2017 15:29:40 +0200 Subject: [PATCH 5/9] Sync AIREBO USER-OMP implementation. --- src/MANYBODY/pair_airebo.cpp | 3 +- src/USER-OMP/pair_airebo_omp.cpp | 403 +++++++++++++++---------------- 2 files changed, 200 insertions(+), 206 deletions(-) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 27a5aad21c..a33abfcb63 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -2147,7 +2147,6 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij_mod[3], double rijmag_mo double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3]; double tmp3pij,tmp3pji,Stb,dStb; - double **x = atom->x; int *type = atom->type; @@ -2261,7 +2260,7 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij_mod[3], double rijmag_mo if (itype == 0 && jtype == 0) Tij=TijSpline((NijC+NijH),(NjiC+NjiH),Nijconj,dN3Tij); Etmp = 0.0; - + if (fabs(Tij) > TOL) { atom2 = atomi; atom3 = atomj; diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp index 13df133585..2fd6c93f03 100644 --- a/src/USER-OMP/pair_airebo_omp.cpp +++ b/src/USER-OMP/pair_airebo_omp.cpp @@ -504,10 +504,6 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, // compute LJ forces and energy - sigwid = 0.84; - sigcut = 3.0; - sigmin = sigcut - sigwid; - rljmin = sigma[itype][jtype]; rljmax = sigcut * rljmin; rljmin = sigmin * rljmin; @@ -1853,44 +1849,62 @@ double PairAIREBOOMP::bondorder_thr(int i, int j, double rij[3], double rijmag, /* ---------------------------------------------------------------------- Bij* function -------------------------------------------------------------------------- */ +------------------------------------------------------------------------- -double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag, - double VA, double rij0[3], double rij0mag, +This function calculates S(t_b(b_ij*)) as specified in the AIREBO paper. +To do so, it needs to compute b_ij*, i.e. the bondorder given that the +atoms i and j are placed a ficticious distance rijmag_mod apart. +Now there are two approaches to calculate the resulting forces: +1. Carry through the ficticious distance and corresponding vector + rij_mod, correcting afterwards using the derivative of r/|r|. +2. Perform all the calculations using the real distance, and do not + use a correction, only using rijmag_mod where necessary. +This code opts for (2). Mathematically, the approaches are equivalent +if implemented correctly, since in terms where only the normalized +vector is used, both calculations necessarily lead to the same result +since if f(x) = g(x/|x|) then for x = y/|y| f(x) = g(y/|y|/1). +The actual modified distance is only used in the lamda terms. +Note that these do not contribute to the forces between i and j, since +rijmag_mod is a constant and the corresponding derivatives are +accordingly zero. +This function should be kept in sync with bondorder(), i.e. changes +there probably also need to be performed here. + +*/ + +double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij_mod[3], double rijmag_mod, + double VA, double rij[3], double rijmag, int vflag_atom, ThrData * const thr) { - int k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; - int atomi,atomj,itype,jtype,ktype,ltype,ntype; - double rik[3], rjl[3], rkn[3],rknmag,dNki; + int atomi,atomj,k,n,l,atomk,atoml,atomn,atom1,atom2,atom3,atom4; + int itype,jtype,ktype,ltype,ntype; + double rik[3],rjl[3],rkn[3],rji[3],rki[3],rlj[3],rknmag,dNki,dwjl,bij; double NijC,NijH,NjiC,NjiH,wik,dwik,dwkn,wjl; double rikmag,rjlmag,cosjik,cosijl,g,tmp2,tmp3; - double Etmp,pij,tmp,wij,dwij,NconjtmpI,NconjtmpJ; - double Nki,Nlj,dS,lamdajik,lamdaijl,dgdc,dgdN,pji,Nijconj,piRC; - double dN2[2],dN3[3],dwjl; - double Tij,crosskij[3],crosskijmag; - double crossijl[3],crossijlmag,omkijl; - double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3]; - double bij,tmp3pij,tmp3pji,Stb,dStb; - double r32[3],r32mag,cos321; + double Etmp,pij,tmp,wij,dwij,NconjtmpI,NconjtmpJ,Nki,Nlj,dS; + double lamdajik,lamdaijl,dgdc,dgdN,pji,Nijconj,piRC; + double dcosjikdri[3],dcosijldri[3],dcosjikdrk[3]; + double dN2[2],dN3[3]; + double dcosjikdrj[3],dcosijldrj[3],dcosijldrl[3]; + double Tij; + double r32[3],r32mag,cos321,r43[3],r13[3]; + double dNlj; double om1234,rln[3]; double rlnmag,dwln,r23[3],r23mag,r21[3],r21mag; double w21,dw21,r34[3],r34mag,cos234,w34,dw34; double cross321[3],cross234[3],prefactor,SpN; - double fcikpc,fcjlpc,fcjkpc,fcilpc; - double dt2dik[3],dt2djl[3],aa,aaa2,at2,cw,cwnum,cwnom; + double fcikpc,fcjlpc,fcjkpc,fcilpc,fcijpc; + double dt2dik[3],dt2djl[3],dt2dij[3],aa,aaa1,aaa2,at2,cw,cwnum,cwnom; double sin321,sin234,rr,rijrik,rijrjl,rjk2,rik2,ril2,rjl2; - double dctik,dctjk,dctjl,dctil,rik2i,rjl2i,sink2i,sinl2i; - double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil; - double dNlj; - double PijS,PjiS; + double dctik,dctjk,dctjl,dctij,dctji,dctil,rik2i,rjl2i,sink2i,sinl2i; + double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil,dt1dij; + double F23[3],F12[3],F34[3],F31[3],F24[3],fi[3],fj[3],fk[3],fl[3]; + double f1[3],f2[3],f3[3],f4[4]; + double dcut321,PijS,PjiS; double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp; int *REBO_neighs,*REBO_neighs_i,*REBO_neighs_j,*REBO_neighs_k,*REBO_neighs_l; - double F12[3],F34[3],F31[3],F24[3]; - double fi[3],fj[3],fk[3],fl[3],f1[3],f2[3],f3[3],f4[4]; - double rji[3],rki[3],rlj[3],r13[3],r43[3]; - double realrij[3], realrijmag; - double rjkmag, rilmag, dctdjk, dctdik, dctdil, dctdjl; - double fjk[3], fil[3], rijmbr; + double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3]; + double tmp3pij,tmp3pji,Stb,dStb; const double * const * const x = atom->x; double * const * const f = thr->get_f(); @@ -1900,12 +1914,11 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag atomj = j; itype = map[type[atomi]]; jtype = map[type[atomj]]; - wij = Sp(rij0mag,rcmin[itype][jtype],rcmax[itype][jtype],dwij); + wij = Sp(rijmag,rcmin[itype][jtype],rcmax[itype][jtype],dwij); NijC = nC[atomi]-(wij*kronecker(jtype,0)); NijH = nH[atomi]-(wij*kronecker(jtype,1)); NjiC = nC[atomj]-(wij*kronecker(itype,0)); NjiH = nH[atomj]-(wij*kronecker(itype,1)); - bij = 0.0; tmp = 0.0; tmp2 = 0.0; @@ -1918,12 +1931,6 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag Stb = 0.0; dStb = 0.0; - realrij[0] = x[atomi][0] - x[atomj][0]; - realrij[1] = x[atomi][1] - x[atomj][1]; - realrij[2] = x[atomi][2] - x[atomj][2]; - realrijmag = sqrt(realrij[0] * realrij[0] + realrij[1] * realrij[1] - + realrij[2] * realrij[2]); - REBO_neighs = REBO_firstneigh[i]; for (k = 0; k < REBO_numneigh[i]; k++) { atomk = REBO_neighs[k]; @@ -1934,9 +1941,9 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag rik[2] = x[atomi][2]-x[atomk][2]; rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2])); lamdajik = 4.0*kronecker(itype,1) * - ((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag)); + ((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag_mod)); wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dS); - Nki = nC[atomk]-(wik*kronecker(itype,0)) + + Nki = nC[atomk]-(wik*kronecker(itype,0))+ nH[atomk]-(wik*kronecker(itype,1)); cosjik = ((rij[0]*rik[0])+(rij[1]*rik[1])+(rij[2]*rik[2])) / (rijmag*rikmag); @@ -1959,6 +1966,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag pij = 1.0/sqrt(1.0+Etmp+PijS); tmppij = -.5*pij*pij*pij; tmp3pij = tmp3; + tmp = 0.0; tmp2 = 0.0; tmp3 = 0.0; @@ -1974,7 +1982,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag rjl[2] = x[atomj][2]-x[atoml][2]; rjlmag = sqrt((rjl[0]*rjl[0])+(rjl[1]*rjl[1])+(rjl[2]*rjl[2])); lamdaijl = 4.0*kronecker(jtype,1) * - ((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag)); + ((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag_mod)); wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dS); Nlj = nC[atoml]-(wjl*kronecker(jtype,0))+nH[atoml] - (wjl*kronecker(jtype,1)); @@ -2004,82 +2012,80 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag Nijconj = 1.0+(NconjtmpI*NconjtmpI)+(NconjtmpJ*NconjtmpJ); piRC = piRCSpline(NijC+NijH,NjiC+NjiH,Nijconj,itype,jtype,dN3piRC); + Tij = 0.0; dN3Tij[0] = 0.0; dN3Tij[1] = 0.0; dN3Tij[2] = 0.0; if (itype == 0 && jtype == 0) Tij=TijSpline((NijC+NijH),(NjiC+NjiH),Nijconj,dN3Tij); - Etmp = 0.0; + if (fabs(Tij) > TOL) { + atom2 = atomi; + atom3 = atomj; + r32[0] = x[atom3][0]-x[atom2][0]; + r32[1] = x[atom3][1]-x[atom2][1]; + r32[2] = x[atom3][2]-x[atom2][2]; + r32mag = sqrt((r32[0]*r32[0])+(r32[1]*r32[1])+(r32[2]*r32[2])); + r23[0] = -r32[0]; + r23[1] = -r32[1]; + r23[2] = -r32[2]; + r23mag = r32mag; REBO_neighs_i = REBO_firstneigh[i]; for (k = 0; k < REBO_numneigh[i]; k++) { atomk = REBO_neighs_i[k]; + atom1 = atomk; ktype = map[type[atomk]]; if (atomk != atomj) { - rik[0] = x[atomi][0]-x[atomk][0]; - rik[1] = x[atomi][1]-x[atomk][1]; - rik[2] = x[atomi][2]-x[atomk][2]; - rikmag = sqrt((rik[0]*rik[0])+(rik[1]*rik[1])+(rik[2]*rik[2])); - cos321 = ((rij[0]*rik[0])+(rij[1]*rik[1])+(rij[2]*rik[2])) / - (rijmag*rikmag); + r21[0] = x[atom2][0]-x[atom1][0]; + r21[1] = x[atom2][1]-x[atom1][1]; + r21[2] = x[atom2][2]-x[atom1][2]; + r21mag = sqrt(r21[0]*r21[0] + r21[1]*r21[1] + r21[2]*r21[2]); + cos321 = -1.0*((r21[0]*r32[0])+(r21[1]*r32[1])+(r21[2]*r32[2])) / + (r21mag*r32mag); cos321 = MIN(cos321,1.0); cos321 = MAX(cos321,-1.0); + sin321 = sqrt(1.0 - cos321*cos321); + if ((sin321 > TOL) && (r21mag > TOL)) { // XXX was sin321 != 0.0 + w21 = Sp(r21mag,rcmin[itype][ktype],rcmaxp[itype][ktype],dw21); + tspjik = Sp2(cos321,thmin,thmax,dtsjik); - rjk[0] = rik[0]-rij[0]; - rjk[1] = rik[1]-rij[1]; - rjk[2] = rik[2]-rij[2]; - rjk2 = (rjk[0]*rjk[0])+(rjk[1]*rjk[1])+(rjk[2]*rjk[2]); - rij2 = rijmag*rijmag; - rik2 = rikmag*rikmag; - costmp = 0.5*(rij2+rik2-rjk2)/rijmag/rikmag; - tspjik = Sp2(costmp,thmin,thmax,dtsjik); - - if (sqrt(1.0 - cos321*cos321) > sqrt(TOL)) { - wik = Sp(rikmag,rcmin[itype][ktype],rcmaxp[itype][ktype],dwik); REBO_neighs_j = REBO_firstneigh[j]; for (l = 0; l < REBO_numneigh[j]; l++) { atoml = REBO_neighs_j[l]; + atom4 = atoml; ltype = map[type[atoml]]; if (!(atoml == atomi || atoml == atomk)) { - rjl[0] = x[atomj][0]-x[atoml][0]; - rjl[1] = x[atomj][1]-x[atoml][1]; - rjl[2] = x[atomj][2]-x[atoml][2]; - rjlmag = sqrt(rjl[0]*rjl[0] + rjl[1]*rjl[1] + rjl[2]*rjl[2]); - cos234 = -((rij[0]*rjl[0])+(rij[1]*rjl[1])+(rij[2]*rjl[2])) / - (rijmag*rjlmag); + r34[0] = x[atom3][0]-x[atom4][0]; + r34[1] = x[atom3][1]-x[atom4][1]; + r34[2] = x[atom3][2]-x[atom4][2]; + r34mag = sqrt((r34[0]*r34[0])+(r34[1]*r34[1])+(r34[2]*r34[2])); + cos234 = (r32[0]*r34[0] + r32[1]*r34[1] + r32[2]*r34[2]) / + (r32mag*r34mag); cos234 = MIN(cos234,1.0); cos234 = MAX(cos234,-1.0); + sin234 = sqrt(1.0 - cos234*cos234); - ril[0] = rij[0]+rjl[0]; - ril[1] = rij[1]+rjl[1]; - ril[2] = rij[2]+rjl[2]; - ril2 = (ril[0]*ril[0])+(ril[1]*ril[1])+(ril[2]*ril[2]); - rjl2 = rjlmag*rjlmag; - costmp = 0.5*(rij2+rjl2-ril2)/rijmag/rjlmag; - tspijl = Sp2(costmp,thmin,thmax,dtsijl); + if ((sin234 > TOL) && (r34mag > TOL)) { // XXX was sin234 != 0.0 + w34 = Sp(r34mag,rcmin[jtype][ltype],rcmaxp[jtype][ltype],dw34); + tspijl = Sp2(cos234,thmin,thmax,dtsijl); - if (sqrt(1.0 - cos234*cos234) > sqrt(TOL)) { - wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmaxp[jtype][ltype],dS); - crosskij[0] = (rij[1]*rik[2]-rij[2]*rik[1]); - crosskij[1] = (rij[2]*rik[0]-rij[0]*rik[2]); - crosskij[2] = (rij[0]*rik[1]-rij[1]*rik[0]); - crosskijmag = sqrt(crosskij[0]*crosskij[0] + - crosskij[1]*crosskij[1] + - crosskij[2]*crosskij[2]); - crossijl[0] = (rij[1]*rjl[2]-rij[2]*rjl[1]); - crossijl[1] = (rij[2]*rjl[0]-rij[0]*rjl[2]); - crossijl[2] = (rij[0]*rjl[1]-rij[1]*rjl[0]); - crossijlmag = sqrt(crossijl[0]*crossijl[0] + - crossijl[1]*crossijl[1] + - crossijl[2]*crossijl[2]); - omkijl = -1.0*(((crosskij[0]*crossijl[0]) + - (crosskij[1]*crossijl[1]) + - (crosskij[2]*crossijl[2])) / - (crosskijmag*crossijlmag)); - Etmp += ((1.0-square(omkijl))*wik*wjl) * + cross321[0] = (r32[1]*r21[2])-(r32[2]*r21[1]); + cross321[1] = (r32[2]*r21[0])-(r32[0]*r21[2]); + cross321[2] = (r32[0]*r21[1])-(r32[1]*r21[0]); + cross234[0] = (r23[1]*r34[2])-(r23[2]*r34[1]); + cross234[1] = (r23[2]*r34[0])-(r23[0]*r34[2]); + cross234[2] = (r23[0]*r34[1])-(r23[1]*r34[0]); + + cwnum = (cross321[0]*cross234[0]) + + (cross321[1]*cross234[1]) + (cross321[2]*cross234[2]); + cwnom = r21mag*r34mag*r23mag*r23mag*sin321*sin234; + om1234 = cwnum/cwnom; + cw = om1234; + Etmp += ((1.0-square(om1234))*w21*w34) * (1.0-tspjik)*(1.0-tspijl); + } } } @@ -2111,50 +2117,43 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag rik[2] = x[atomi][2]-x[atomk][2]; rikmag = sqrt(rik[0]*rik[0] + rik[1]*rik[1] + rik[2]*rik[2]); lamdajik = 4.0*kronecker(itype,1) * - ((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag)); + ((rho[ktype][1]-rikmag)-(rho[jtype][1]-rijmag_mod)); wik = Sp(rikmag,rcmin[itype][ktype],rcmax[itype][ktype],dwik); - rjk[0] = rik[0] - rij[0]; - rjk[1] = rik[1] - rij[1]; - rjk[2] = rik[2] - rij[2]; - rjkmag = sqrt(rjk[0] * rjk[0] + rjk[1] * rjk[1] + rjk[2] * rjk[2]); - rijrik = 2 * rijmag * rikmag; - rr = rijmag * rijmag - rikmag * rikmag; - cosjik = (rijmag * rijmag + rikmag * rikmag - rjkmag * rjkmag) / rijrik; + cosjik = (rij[0]*rik[0] + rij[1]*rik[1] + rij[2]*rik[2]) / + (rijmag*rikmag); cosjik = MIN(cosjik,1.0); cosjik = MAX(cosjik,-1.0); - dctdjk = -2 / rijrik; - dctdik = (-rr + rjkmag * rjkmag) / (rijrik * rikmag * rikmag); - // evaluate splines g and derivatives dg + + dcosjikdri[0] = ((rij[0]+rik[0])/(rijmag*rikmag)) - + (cosjik*((rij[0]/(rijmag*rijmag))+(rik[0]/(rikmag*rikmag)))); + dcosjikdri[1] = ((rij[1]+rik[1])/(rijmag*rikmag)) - + (cosjik*((rij[1]/(rijmag*rijmag))+(rik[1]/(rikmag*rikmag)))); + dcosjikdri[2] = ((rij[2]+rik[2])/(rijmag*rikmag)) - + (cosjik*((rij[2]/(rijmag*rijmag))+(rik[2]/(rikmag*rikmag)))); + dcosjikdrk[0] = (-rij[0]/(rijmag*rikmag)) + + (cosjik*(rik[0]/(rikmag*rikmag))); + dcosjikdrk[1] = (-rij[1]/(rijmag*rikmag)) + + (cosjik*(rik[1]/(rikmag*rikmag))); + dcosjikdrk[2] = (-rij[2]/(rijmag*rikmag)) + + (cosjik*(rik[2]/(rikmag*rikmag))); + dcosjikdrj[0] = (-rik[0]/(rijmag*rikmag)) + + (cosjik*(rij[0]/(rijmag*rijmag))); + dcosjikdrj[1] = (-rik[1]/(rijmag*rikmag)) + + (cosjik*(rij[1]/(rijmag*rijmag))); + dcosjikdrj[2] = (-rik[2]/(rijmag*rikmag)) + + (cosjik*(rij[2]/(rijmag*rijmag))); g = gSpline(cosjik,(NijC+NijH),itype,&dgdc,&dgdN); - tmp2 = VA*.5*(tmp*wik*dgdc*exp(lamdajik)); - - fi[0] = -tmp2 * dctdik * rik[0]; - fi[1] = -tmp2 * dctdik * rik[1]; - fi[2] = -tmp2 * dctdik * rik[2]; - fk[0] = tmp2 * dctdik * rik[0]; - fk[1] = tmp2 * dctdik * rik[1]; - fk[2] = tmp2 * dctdik * rik[2]; - fj[0] = 0; - fj[1] = 0; - fj[2] = 0; - fjk[0] = -tmp2 * dctdjk * rjk[0]; - fjk[1] = -tmp2 * dctdjk * rjk[1]; - fjk[2] = -tmp2 * dctdjk * rjk[2]; - fi[0] += fjk[0]; - fi[1] += fjk[1]; - fi[2] += fjk[2]; - fk[0] -= fjk[0]; - fk[1] -= fjk[1]; - fk[2] -= fjk[2]; - rijmbr = rcmin[itype][jtype] / realrijmag; - fj[0] += rijmbr * (fjk[0] - (realrij[0] * realrij[0] * fjk[0] + realrij[0] * realrij[1] * fjk[1] + realrij[0] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fj[1] += rijmbr * (fjk[1] - (realrij[1] * realrij[0] * fjk[0] + realrij[1] * realrij[1] * fjk[1] + realrij[1] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fj[2] += rijmbr * (fjk[2] - (realrij[2] * realrij[0] * fjk[0] + realrij[2] * realrij[1] * fjk[1] + realrij[2] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fi[0] -= rijmbr * (fjk[0] - (realrij[0] * realrij[0] * fjk[0] + realrij[0] * realrij[1] * fjk[1] + realrij[0] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fi[1] -= rijmbr * (fjk[1] - (realrij[1] * realrij[0] * fjk[0] + realrij[1] * realrij[1] * fjk[1] + realrij[1] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); - fi[2] -= rijmbr * (fjk[2] - (realrij[2] * realrij[0] * fjk[0] + realrij[2] * realrij[1] * fjk[1] + realrij[2] * realrij[2] * fjk[2]) / (realrijmag * realrijmag)); + fj[0] = -tmp2*dcosjikdrj[0]; + fj[1] = -tmp2*dcosjikdrj[1]; + fj[2] = -tmp2*dcosjikdrj[2]; + fi[0] = -tmp2*dcosjikdri[0]; + fi[1] = -tmp2*dcosjikdri[1]; + fi[2] = -tmp2*dcosjikdri[2]; + fk[0] = -tmp2*dcosjikdrk[0]; + fk[1] = -tmp2*dcosjikdrk[1]; + fk[2] = -tmp2*dcosjikdrk[2]; tmp2 = VA*.5*(tmp*wik*g*exp(lamdajik)*4.0*kronecker(itype,1)); fi[0] += tmp2*(rik[0]/rikmag); @@ -2212,6 +2211,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag tmp3 = tmp3pji; dN2[0] = dN2PJI[0]; dN2[1] = dN2PJI[1]; + REBO_neighs = REBO_firstneigh[j]; for (l = 0; l < REBO_numneigh[j]; l++) { atoml = REBO_neighs[l]; @@ -2222,50 +2222,43 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag rjl[2] = x[atomj][2]-x[atoml][2]; rjlmag = sqrt((rjl[0]*rjl[0])+(rjl[1]*rjl[1])+(rjl[2]*rjl[2])); lamdaijl = 4.0*kronecker(jtype,1) * - ((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag)); + ((rho[ltype][1]-rjlmag)-(rho[itype][1]-rijmag_mod)); wjl = Sp(rjlmag,rcmin[jtype][ltype],rcmax[jtype][ltype],dwjl); - ril[0] = rij[0] + rjl[0]; - ril[1] = rij[1] + rjl[1]; - ril[2] = rij[2] + rjl[2]; - rilmag = sqrt(ril[0] * ril[0] + ril[1] * ril[1] + ril[2] * ril[2]); - rijrjl = 2 * rijmag * rjlmag; - rr = rijmag * rijmag - rjlmag * rjlmag; - cosijl = (rijmag * rijmag + rjlmag * rjlmag - rilmag * rilmag) / rijrjl; + cosijl = (-1.0*((rij[0]*rjl[0])+(rij[1]*rjl[1])+(rij[2]*rjl[2]))) / + (rijmag*rjlmag); cosijl = MIN(cosijl,1.0); cosijl = MAX(cosijl,-1.0); - dctdil = -2 / rijrjl; - dctdjl = (-rr + rilmag * rilmag) / (rijrjl * rjlmag * rjlmag); + + dcosijldri[0] = (-rjl[0]/(rijmag*rjlmag)) - + (cosijl*rij[0]/(rijmag*rijmag)); + dcosijldri[1] = (-rjl[1]/(rijmag*rjlmag)) - + (cosijl*rij[1]/(rijmag*rijmag)); + dcosijldri[2] = (-rjl[2]/(rijmag*rjlmag)) - + (cosijl*rij[2]/(rijmag*rijmag)); + dcosijldrj[0] = ((-rij[0]+rjl[0])/(rijmag*rjlmag)) + + (cosijl*((rij[0]/square(rijmag))-(rjl[0]/(rjlmag*rjlmag)))); + dcosijldrj[1] = ((-rij[1]+rjl[1])/(rijmag*rjlmag)) + + (cosijl*((rij[1]/square(rijmag))-(rjl[1]/(rjlmag*rjlmag)))); + dcosijldrj[2] = ((-rij[2]+rjl[2])/(rijmag*rjlmag)) + + (cosijl*((rij[2]/square(rijmag))-(rjl[2]/(rjlmag*rjlmag)))); + dcosijldrl[0] = (rij[0]/(rijmag*rjlmag))+(cosijl*rjl[0]/(rjlmag*rjlmag)); + dcosijldrl[1] = (rij[1]/(rijmag*rjlmag))+(cosijl*rjl[1]/(rjlmag*rjlmag)); + dcosijldrl[2] = (rij[2]/(rijmag*rjlmag))+(cosijl*rjl[2]/(rjlmag*rjlmag)); // evaluate splines g and derivatives dg g = gSpline(cosijl,NjiC+NjiH,jtype,&dgdc,&dgdN); tmp2 = VA*.5*(tmp*wjl*dgdc*exp(lamdaijl)); - fj[0] = -tmp2 * dctdjl * rjl[0]; - fj[1] = -tmp2 * dctdjl * rjl[1]; - fj[2] = -tmp2 * dctdjl * rjl[2]; - fl[0] = tmp2 * dctdjl * rjl[0]; - fl[1] = tmp2 * dctdjl * rjl[1]; - fl[2] = tmp2 * dctdjl * rjl[2]; - fi[0] = 0; - fi[1] = 0; - fi[2] = 0; - fil[0] = -tmp2 * dctdil * ril[0]; - fil[1] = -tmp2 * dctdil * ril[1]; - fil[2] = -tmp2 * dctdil * ril[2]; - fj[0] += fil[0]; - fj[1] += fil[1]; - fj[2] += fil[2]; - fl[0] -= fil[0]; - fl[1] -= fil[1]; - fl[2] -= fil[2]; - rijmbr = rcmin[itype][jtype] / realrijmag; - fi[0] += rijmbr * (fil[0] - (realrij[0] * realrij[0] * fil[0] + realrij[0] * realrij[1] * fil[1] + realrij[0] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fi[1] += rijmbr * (fil[1] - (realrij[1] * realrij[0] * fil[0] + realrij[1] * realrij[1] * fil[1] + realrij[1] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fi[2] += rijmbr * (fil[2] - (realrij[2] * realrij[0] * fil[0] + realrij[2] * realrij[1] * fil[1] + realrij[2] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fj[0] -= rijmbr * (fil[0] - (realrij[0] * realrij[0] * fil[0] + realrij[0] * realrij[1] * fil[1] + realrij[0] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fj[1] -= rijmbr * (fil[1] - (realrij[1] * realrij[0] * fil[0] + realrij[1] * realrij[1] * fil[1] + realrij[1] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - fj[2] -= rijmbr * (fil[2] - (realrij[2] * realrij[0] * fil[0] + realrij[2] * realrij[1] * fil[1] + realrij[2] * realrij[2] * fil[2]) / (realrijmag * realrijmag)); - + fi[0] = -tmp2*dcosijldri[0]; + fi[1] = -tmp2*dcosijldri[1]; + fi[2] = -tmp2*dcosijldri[2]; + fj[0] = -tmp2*dcosijldrj[0]; + fj[1] = -tmp2*dcosijldrj[1]; + fj[2] = -tmp2*dcosijldrj[2]; + fl[0] = -tmp2*dcosijldrl[0]; + fl[1] = -tmp2*dcosijldrl[1]; + fl[2] = -tmp2*dcosijldrl[2]; + tmp2 = VA*.5*(tmp*wjl*g*exp(lamdaijl)*4.0*kronecker(jtype,1)); fj[0] += tmp2*(rjl[0]/rjlmag); fj[1] += tmp2*(rjl[1]/rjlmag); @@ -2275,6 +2268,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag fl[2] -= tmp2*(rjl[2]/rjlmag); // coordination forces + // dwik forces tmp2 = VA*.5*(tmp*dwjl*g*exp(lamdaijl))/rjlmag; @@ -2389,8 +2383,8 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag // piRC forces to J side - REBO_neighs = REBO_firstneigh[j]; - for (l = 0; l < REBO_numneigh[j]; l++) { + REBO_neighs = REBO_firstneigh[atomj]; + for (l = 0; l < REBO_numneigh[atomj]; l++) { atoml = REBO_neighs[l]; if (atoml != atomi) { ltype = map[type[atoml]]; @@ -2460,15 +2454,14 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag dN3[2] = dN3Tij[2]; atom2 = atomi; atom3 = atomj; - r32[0] = -rij[0]; - r32[1] = -rij[1]; - r32[2] = -rij[2]; - r23[0] = rij[0]; - r23[1] = rij[1]; - r23[2] = rij[2]; - r32mag = rijmag; - r23mag = rijmag; - + r32[0] = x[atom3][0]-x[atom2][0]; + r32[1] = x[atom3][1]-x[atom2][1]; + r32[2] = x[atom3][2]-x[atom2][2]; + r32mag = sqrt((r32[0]*r32[0])+(r32[1]*r32[1])+(r32[2]*r32[2])); + r23[0] = -r32[0]; + r23[1] = -r32[1]; + r23[2] = -r32[2]; + r23mag = r32mag; REBO_neighs_i = REBO_firstneigh[i]; for (k = 0; k < REBO_numneigh[i]; k++) { atomk = REBO_neighs_i[k]; @@ -2484,20 +2477,21 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag cos321 = MIN(cos321,1.0); cos321 = MAX(cos321,-1.0); sin321 = sqrt(1.0 - cos321*cos321); - if ((sin321 > TOL) && (r21mag > TOL)) { // XXX was sin321 != 0.0 sink2i = 1.0/(sin321*sin321); rik2i = 1.0/(r21mag*r21mag); rr = (rijmag*rijmag)-(r21mag*r21mag); - rjk[0] = r21[0]-rij[0]; - rjk[1] = r21[1]-rij[1]; - rjk[2] = r21[2]-rij[2]; + rjk[0] = r21[0]-r23[0]; + rjk[1] = r21[1]-r23[1]; + rjk[2] = r21[2]-r23[2]; rjk2 = (rjk[0]*rjk[0])+(rjk[1]*rjk[1])+(rjk[2]*rjk[2]); - rijrik = 2.0*rijmag*r21mag; + rijrik = 2.0*r23mag*r21mag; rik2 = r21mag*r21mag; dctik = (-rr+rjk2)/(rijrik*rik2); + dctij = (rr+rjk2)/(rijrik*r23mag*r23mag); dctjk = -2.0/rijrik; w21 = Sp(r21mag,rcmin[itype][ktype],rcmaxp[itype][ktype],dw21); + rijmag = r32mag; rikmag = r21mag; rij2 = r32mag*r32mag; rik2 = r21mag*r21mag; @@ -2515,8 +2509,8 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag r34[1] = x[atom3][1]-x[atom4][1]; r34[2] = x[atom3][2]-x[atom4][2]; r34mag = sqrt(r34[0]*r34[0] + r34[1]*r34[1] + r34[2]*r34[2]); - cos234 = -1.0*((rij[0]*r34[0])+(rij[1]*r34[1]) + - (rij[2]*r34[2]))/(rijmag*r34mag); + cos234 = (r32[0]*r34[0] + r32[1]*r34[1] + r32[2]*r34[2]) / + (r32mag*r34mag); cos234 = MIN(cos234,1.0); cos234 = MAX(cos234,-1.0); sin234 = sqrt(1.0 - cos234*cos234); @@ -2534,6 +2528,7 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag rijrjl = 2.0*r23mag*r34mag; rjl2 = r34mag*r34mag; dctjl = (-rr+ril2)/(rijrjl*rjl2); + dctji = (rr+ril2)/(rijrjl*r23mag*r23mag); dctil = -2.0/rijrjl; rjlmag = r34mag; rjl2 = r34mag*r34mag; @@ -2559,6 +2554,8 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag dt1djk = (-dctjk*sink2i*cos321); dt1djl = (rjl2i)-(dctjl*sinl2i*cos234); dt1dil = (-dctil*sinl2i*cos234); + dt1dij = (2.0/(r23mag*r23mag))-(dctij*sink2i*cos321) - + (dctji*sinl2i*cos234); dt2dik[0] = (-r23[2]*cross234[1])+(r23[1]*cross234[2]); dt2dik[1] = (-r23[0]*cross234[2])+(r23[2]*cross234[0]); @@ -2568,16 +2565,29 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag dt2djl[1] = (-r23[2]*cross321[0])+(r23[0]*cross321[2]); dt2djl[2] = (-r23[0]*cross321[1])+(r23[1]*cross321[0]); + dt2dij[0] = (r21[2]*cross234[1])-(r34[2]*cross321[1]) - + (r21[1]*cross234[2])+(r34[1]*cross321[2]); + dt2dij[1] = (r21[0]*cross234[2])-(r34[0]*cross321[2]) - + (r21[2]*cross234[0])+(r34[2]*cross321[0]); + dt2dij[2] = (r21[1]*cross234[0])-(r34[1]*cross321[0]) - + (r21[0]*cross234[1])+(r34[0]*cross321[1]); + aa = (prefactor*2.0*cw/cwnom)*w21*w34 * (1.0-tspjik)*(1.0-tspijl); aaa2 = -prefactor*(1.0-square(om1234)) * w21*w34; at2 = aa*cwnum; + fcijpc = (-dt1dij*at2)+(aaa2*dtsjik*dctij*(1.0-tspijl)) + + (aaa2*dtsijl*dctji*(1.0-tspjik)); fcikpc = (-dt1dik*at2)+(aaa2*dtsjik*dctik*(1.0-tspijl)); fcjlpc = (-dt1djl*at2)+(aaa2*dtsijl*dctjl*(1.0-tspjik)); fcjkpc = (-dt1djk*at2)+(aaa2*dtsjik*dctjk*(1.0-tspijl)); fcilpc = (-dt1dil*at2)+(aaa2*dtsijl*dctil*(1.0-tspjik)); + F23[0] = (fcijpc*r23[0])+(aa*dt2dij[0]); + F23[1] = (fcijpc*r23[1])+(aa*dt2dij[1]); + F23[2] = (fcijpc*r23[2])+(aa*dt2dij[2]); + F12[0] = (fcikpc*r21[0])+(aa*dt2dik[0]); F12[1] = (fcikpc*r21[1])+(aa*dt2dik[1]); F12[2] = (fcikpc*r21[2])+(aa*dt2dik[2]); @@ -2597,31 +2607,16 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij[3], double rijmag f1[0] = -F12[0]-F31[0]; f1[1] = -F12[1]-F31[1]; f1[2] = -F12[2]-F31[2]; - f2[0] = F12[0]+F31[0]; - f2[1] = F12[1]+F31[1]; - f2[2] = F12[2]+F31[2]; - f3[0] = F34[0]+F24[0]; - f3[1] = F34[1]+F24[1]; - f3[2] = F34[2]+F24[2]; + f2[0] = F23[0]+F12[0]+F24[0]; + f2[1] = F23[1]+F12[1]+F24[1]; + f2[2] = F23[2]+F12[2]+F24[2]; + f3[0] = -F23[0]+F34[0]+F31[0]; + f3[1] = -F23[1]+F34[1]+F31[1]; + f3[2] = -F23[2]+F34[2]+F31[2]; f4[0] = -F34[0]-F24[0]; f4[1] = -F34[1]-F24[1]; f4[2] = -F34[2]-F24[2]; - rijmbr = rcmin[itype][jtype] / realrijmag; - f2[0] += rijmbr * (F24[0] - (realrij[0] * realrij[0] * F24[0] + realrij[0] * realrij[1] * F24[1] + realrij[0] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f2[1] += rijmbr * (F24[1] - (realrij[1] * realrij[0] * F24[0] + realrij[1] * realrij[1] * F24[1] + realrij[1] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f2[2] += rijmbr * (F24[2] - (realrij[2] * realrij[0] * F24[0] + realrij[2] * realrij[1] * F24[1] + realrij[2] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f3[0] -= rijmbr * (F24[0] - (realrij[0] * realrij[0] * F24[0] + realrij[0] * realrij[1] * F24[1] + realrij[0] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f3[1] -= rijmbr * (F24[1] - (realrij[1] * realrij[0] * F24[0] + realrij[1] * realrij[1] * F24[1] + realrij[1] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - f3[2] -= rijmbr * (F24[2] - (realrij[2] * realrij[0] * F24[0] + realrij[2] * realrij[1] * F24[1] + realrij[2] * realrij[2] * F24[2]) / (realrijmag * realrijmag)); - - f2[0] -= rijmbr * (F31[0] - (realrij[0] * realrij[0] * F31[0] + realrij[0] * realrij[1] * F31[1] + realrij[0] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f2[1] -= rijmbr * (F31[1] - (realrij[1] * realrij[0] * F31[0] + realrij[1] * realrij[1] * F31[1] + realrij[1] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f2[2] -= rijmbr * (F31[2] - (realrij[2] * realrij[0] * F31[0] + realrij[2] * realrij[1] * F31[1] + realrij[2] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f3[0] += rijmbr * (F31[0] - (realrij[0] * realrij[0] * F31[0] + realrij[0] * realrij[1] * F31[1] + realrij[0] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f3[1] += rijmbr * (F31[1] - (realrij[1] * realrij[0] * F31[0] + realrij[1] * realrij[1] * F31[1] + realrij[1] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - f3[2] += rijmbr * (F31[2] - (realrij[2] * realrij[0] * F31[0] + realrij[2] * realrij[1] * F31[1] + realrij[2] * realrij[2] * F31[2]) / (realrijmag * realrijmag)); - // coordination forces tmp2 = VA*Tij*((1.0-(om1234*om1234))) * From 1bf1cb150faefd84910a89be70b6823ff3a34e5c Mon Sep 17 00:00:00 2001 From: Markus Hoehnerbach Date: Wed, 5 Jul 2017 18:26:32 +0200 Subject: [PATCH 6/9] Updated credits --- src/MANYBODY/pair_airebo.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index a33abfcb63..5c022cedb2 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -15,7 +15,7 @@ Contributing author: Ase Henry (MIT) Bugfixes and optimizations: Marcel Fallet & Steve Stuart (Clemson), Axel Kohlmeyer (Temple U), - Markus Hoehnerbach (RWTH Aachen), Github user CF17 + Markus Hoehnerbach (RWTH Aachen), Cyril Falvo (Universite Paris Sud) AIREBO-M modification to optionally replace LJ with Morse potentials. Thomas C. O'Connor (JHU) 2014 ------------------------------------------------------------------------- */ From 894e0c3cf5064ab4f524d6cf0a25e564a65a6c1b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 5 Jul 2017 16:19:24 -0400 Subject: [PATCH 7/9] simplify parsing of optional arguments --- src/MANYBODY/pair_airebo.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 5c022cedb2..26800a75d3 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -157,13 +157,11 @@ void PairAIREBO::settings(int narg, char **arg) cutlj = force->numeric(FLERR,arg[0]); - if (narg == 3) { + if (narg >= 3) { ljflag = force->inumeric(FLERR,arg[1]); torflag = force->inumeric(FLERR,arg[2]); } if (narg == 4) { - ljflag = force->inumeric(FLERR,arg[1]); - torflag = force->inumeric(FLERR,arg[2]); sigcut = cutlj; sigmin = force->numeric(FLERR,arg[3]); sigwid = sigcut - sigmin; From 8c3f6947ad480fd6cf22809331af916b99c76146 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 5 Jul 2017 16:19:59 -0400 Subject: [PATCH 8/9] remove unused variables to silence compiler warnings --- src/MANYBODY/pair_airebo.cpp | 8 +++----- src/USER-OMP/pair_airebo_omp.cpp | 7 +++---- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 26800a75d3..1b248f701c 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -534,7 +534,7 @@ void PairAIREBO::FLJ(int eflag, int vflag) double delij[3],rijsq,delik[3],rik,deljk[3]; double rkj,wkj,dC,VLJ,dVLJ,VA,Str,dStr,Stb; double vdw,slw,dvdw,dslw,drij,swidth,tee,tee2; - double rljmin,rljmax,sigcut,sigmin,sigwid; + double rljmin,rljmax,sigcut,sigmin; double delkm[3],rkm,deljm[3],rmj,wmj,r2inv,r6inv,scale,delscale[3]; int *ilist,*jlist,*numneigh,**firstneigh; int *REBO_neighs_i,*REBO_neighs_k; @@ -552,8 +552,6 @@ void PairAIREBO::FLJ(int eflag, int vflag) rljmax = 0.0; sigcut = 0.0; sigmin = 0.0; - sigwid = 0.0; - double **x = atom->x; double **f = atom->f; @@ -2133,13 +2131,13 @@ double PairAIREBO::bondorderLJ(int i, int j, double rij_mod[3], double rijmag_mo double w21,dw21,r34[3],r34mag,cos234,w34,dw34; double cross321[3],cross234[3],prefactor,SpN; double fcikpc,fcjlpc,fcjkpc,fcilpc,fcijpc; - double dt2dik[3],dt2djl[3],dt2dij[3],aa,aaa1,aaa2,at2,cw,cwnum,cwnom; + double dt2dik[3],dt2djl[3],dt2dij[3],aa,aaa2,at2,cw,cwnum,cwnom; double sin321,sin234,rr,rijrik,rijrjl,rjk2,rik2,ril2,rjl2; double dctik,dctjk,dctjl,dctij,dctji,dctil,rik2i,rjl2i,sink2i,sinl2i; double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil,dt1dij; double F23[3],F12[3],F34[3],F31[3],F24[3],fi[3],fj[3],fk[3],fl[3]; double f1[3],f2[3],f3[3],f4[4]; - double dcut321,PijS,PjiS; + double PijS,PjiS; double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp; int *REBO_neighs,*REBO_neighs_i,*REBO_neighs_j,*REBO_neighs_k,*REBO_neighs_l; double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3]; diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp index 2fd6c93f03..bc83ff8e4b 100644 --- a/src/USER-OMP/pair_airebo_omp.cpp +++ b/src/USER-OMP/pair_airebo_omp.cpp @@ -294,7 +294,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, double delij[3],rijsq,delik[3],rik,deljk[3]; double rkj,wkj,dC,VLJ,dVLJ,VA,Str,dStr,Stb; double vdw,slw,dvdw,dslw,drij,swidth,tee,tee2; - double rljmin,rljmax,sigcut,sigmin,sigwid; + double rljmin,rljmax,sigcut,sigmin; double delkm[3],rkm,deljm[3],rmj,wmj,r2inv,r6inv,scale,delscale[3]; int *ilist,*jlist,*numneigh,**firstneigh; int *REBO_neighs_i,*REBO_neighs_k; @@ -312,7 +312,6 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, rljmax = 0.0; sigcut = 0.0; sigmin = 0.0; - sigwid = 0.0; const double * const * const x = atom->x; double * const * const f = thr->get_f(); @@ -1894,13 +1893,13 @@ double PairAIREBOOMP::bondorderLJ_thr(int i, int j, double rij_mod[3], double ri double w21,dw21,r34[3],r34mag,cos234,w34,dw34; double cross321[3],cross234[3],prefactor,SpN; double fcikpc,fcjlpc,fcjkpc,fcilpc,fcijpc; - double dt2dik[3],dt2djl[3],dt2dij[3],aa,aaa1,aaa2,at2,cw,cwnum,cwnom; + double dt2dik[3],dt2djl[3],dt2dij[3],aa,aaa2,at2,cw,cwnum,cwnom; double sin321,sin234,rr,rijrik,rijrjl,rjk2,rik2,ril2,rjl2; double dctik,dctjk,dctjl,dctij,dctji,dctil,rik2i,rjl2i,sink2i,sinl2i; double rjk[3],ril[3],dt1dik,dt1djk,dt1djl,dt1dil,dt1dij; double F23[3],F12[3],F34[3],F31[3],F24[3],fi[3],fj[3],fk[3],fl[3]; double f1[3],f2[3],f3[3],f4[4]; - double dcut321,PijS,PjiS; + double PijS,PjiS; double rij2,tspjik,dtsjik,tspijl,dtsijl,costmp; int *REBO_neighs,*REBO_neighs_i,*REBO_neighs_j,*REBO_neighs_k,*REBO_neighs_l; double tmppij,tmppji,dN2PIJ[2],dN2PJI[2],dN3piRC[3],dN3Tij[3]; From e493b6a648b268a86654824d09536bd846e618d1 Mon Sep 17 00:00:00 2001 From: Markus Hoehnerbach Date: Wed, 5 Jul 2017 22:52:29 +0200 Subject: [PATCH 9/9] Fix sigcut class variable actually used --- src/MANYBODY/pair_airebo.cpp | 4 +--- src/USER-OMP/pair_airebo_omp.cpp | 4 +--- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 1b248f701c..67e1e5262a 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -534,7 +534,7 @@ void PairAIREBO::FLJ(int eflag, int vflag) double delij[3],rijsq,delik[3],rik,deljk[3]; double rkj,wkj,dC,VLJ,dVLJ,VA,Str,dStr,Stb; double vdw,slw,dvdw,dslw,drij,swidth,tee,tee2; - double rljmin,rljmax,sigcut,sigmin; + double rljmin,rljmax; double delkm[3],rkm,deljm[3],rmj,wmj,r2inv,r6inv,scale,delscale[3]; int *ilist,*jlist,*numneigh,**firstneigh; int *REBO_neighs_i,*REBO_neighs_k; @@ -550,8 +550,6 @@ void PairAIREBO::FLJ(int eflag, int vflag) evdwl = 0.0; rljmin = 0.0; rljmax = 0.0; - sigcut = 0.0; - sigmin = 0.0; double **x = atom->x; double **f = atom->f; diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp index bc83ff8e4b..73529847f2 100644 --- a/src/USER-OMP/pair_airebo_omp.cpp +++ b/src/USER-OMP/pair_airebo_omp.cpp @@ -294,7 +294,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, double delij[3],rijsq,delik[3],rik,deljk[3]; double rkj,wkj,dC,VLJ,dVLJ,VA,Str,dStr,Stb; double vdw,slw,dvdw,dslw,drij,swidth,tee,tee2; - double rljmin,rljmax,sigcut,sigmin; + double rljmin,rljmax; double delkm[3],rkm,deljm[3],rmj,wmj,r2inv,r6inv,scale,delscale[3]; int *ilist,*jlist,*numneigh,**firstneigh; int *REBO_neighs_i,*REBO_neighs_k; @@ -310,8 +310,6 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag, evdwl = 0.0; rljmin = 0.0; rljmax = 0.0; - sigcut = 0.0; - sigmin = 0.0; const double * const * const x = atom->x; double * const * const f = thr->get_f();