diff --git a/src/ML-UF3/uf3_pair_bspline.cpp b/src/ML-UF3/uf3_pair_bspline.cpp deleted file mode 100644 index 92068d117d..0000000000 --- a/src/ML-UF3/uf3_pair_bspline.cpp +++ /dev/null @@ -1,213 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - LAMMPS development team: developers@lammps.org - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "uf3_pair_bspline.h" - -#include "uf3_bspline_basis2.h" -#include "uf3_bspline_basis3.h" - -#include "utils.h" -#include "error.h" -#include - -using namespace LAMMPS_NS; - -// Dummy constructor -uf3_pair_bspline::uf3_pair_bspline() {} - -// Constructor -// Passing vectors by reference -uf3_pair_bspline::uf3_pair_bspline(LAMMPS *ulmp, const std::vector &uknot_vect, - const std::vector &ucoeff_vect, - const int &uknot_spacing_type) -{ - lmp = ulmp; - knot_vect = uknot_vect; - coeff_vect = ucoeff_vect; - - knot_spacing_type = uknot_spacing_type; - if (knot_spacing_type==0){ - knot_spacing = knot_vect[4]-knot_vect[3]; - get_starting_index=&uf3_pair_bspline::get_starting_index_uniform; - } - else if (knot_spacing_type==1){ - knot_spacing = 0; - get_starting_index=&uf3_pair_bspline::get_starting_index_nonuniform; - } - - else - lmp->error->all(FLERR, "UF3: Expected either '0'(uniform-knots) or \n\ - '1'(non-uniform knots)"); - - knot_vect_size = uknot_vect.size(); - coeff_vect_size = ucoeff_vect.size(); - - // Initialize B-Spline Basis Functions - for (int i = 0; i < knot_vect.size() - 4; i++) - bspline_bases.push_back(uf3_bspline_basis3(lmp, &knot_vect[i], coeff_vect[i])); - - // Initialize Coefficients and Knots for Derivatives - // The last coefficient needs to be droped - for (int i = 0; i < coeff_vect_size - 1; i++) { - double dntemp4 = 3 / (knot_vect[i + 4] - knot_vect[i + 1]); - dncoeff_vect.push_back((coeff_vect[i + 1] - coeff_vect[i]) * dntemp4); - } - //What we have is a clamped bspline -->i.e value of the bspline curve at the - //knots with multiplicity equal to the degree of bspline is equal to the coefficient - // - //Therefore for the derivative bspline the very first and last knot needs to be droped - //to change their multiplicity from 4 (necessary condition for clamped cubic bspline) - //to 3 (necessary condition for clamped quadratic bspline) - // - //Also if the coeff vector size of decreases by 1 for the derivative bspline - //knots size needs to go down by 2 as ==> knots = coefficient + degree + 1 - for (int i = 1; i < knot_vect_size - 1; i++) dnknot_vect.push_back(knot_vect[i]); - - // Initialize B-Spline Derivative Basis Functions - for (int i = 0; i < dnknot_vect.size() - 3; i++) - dnbspline_bases.push_back(uf3_bspline_basis2(lmp, &dnknot_vect[i], dncoeff_vect[i])); -} - -// Constructor -// Passing arrays -uf3_pair_bspline::uf3_pair_bspline(LAMMPS *ulmp, const double* uknot_array, - const int uknot_array_size, - const double* ucoeff_array, - const int ucoeff_array_size, - const int uknot_spacing_type) -{ - lmp = ulmp; - - knot_vect = std::vector (uknot_array, uknot_array + uknot_array_size); - coeff_vect = std::vector (ucoeff_array, ucoeff_array + ucoeff_array_size); - - knot_spacing_type = uknot_spacing_type; - if (knot_spacing_type==0){ - knot_spacing = knot_vect[4]-knot_vect[3]; - get_starting_index=&uf3_pair_bspline::get_starting_index_uniform; - } - else if (knot_spacing_type==1){ - knot_spacing = 0; - get_starting_index=&uf3_pair_bspline::get_starting_index_nonuniform; - } - - else - lmp->error->all(FLERR, "UF3: Expected either '0'(uniform-knots) or \n\ - '1'(non-uniform knots)"); - - knot_vect_size = uknot_array_size; - coeff_vect_size = ucoeff_array_size; - - // Initialize B-Spline Basis Functions - for (int i = 0; i < knot_vect.size() - 4; i++) - bspline_bases.push_back(uf3_bspline_basis3(lmp, &knot_vect[i], coeff_vect[i])); - - // Initialize Coefficients and Knots for Derivatives - // The last coefficient needs to be droped - for (int i = 0; i < coeff_vect_size - 1; i++) { - double dntemp4 = 3 / (knot_vect[i + 4] - knot_vect[i + 1]); - dncoeff_vect.push_back((coeff_vect[i + 1] - coeff_vect[i]) * dntemp4); - } - //What we have is a clamped bspline -->i.e value of the bspline curve at the - //knots with multiplicity equal to the degree of bspline is equal to the coefficient - // - //Therefore for the derivative bspline the very first and last knot needs to be droped - //to change their multiplicity from 4 (necessary condition for clamped cubic bspline) - //to 3 (necessary condition for clamped quadratic bspline) - // - //Also if the coeff vector size of decreases by 1 for the derivative bspline - //knots size needs to go down by 2 as ==> knots = coefficient + degree + 1 - for (int i = 1; i < knot_vect_size - 1; i++) dnknot_vect.push_back(knot_vect[i]); - - // Initialize B-Spline Derivative Basis Functions - for (int i = 0; i < dnknot_vect.size() - 3; i++) - dnbspline_bases.push_back(uf3_bspline_basis2(lmp, &dnknot_vect[i], dncoeff_vect[i])); -} - -uf3_pair_bspline::~uf3_pair_bspline() {} - -int uf3_pair_bspline::get_starting_index_uniform(double r) -{ - return 3+(int)((r-knot_vect[0])/knot_spacing); -} - -int uf3_pair_bspline::get_starting_index_nonuniform(double r) -{ - if (knot_vect.front() <= r && r < knot_vect.back()) { - //Determine the interval for value_rij - for (int i = 3; i < knot_vect_size - 1; ++i) { - if (knot_vect[i] <= r && r < knot_vect[i + 1]) { - return i; - } - } - } -} - -double *uf3_pair_bspline::eval(double r) -{ - - // Find knot starting position - - int start_index=(this->*get_starting_index)(r); - /*if (knot_vect.front() <= r && r < knot_vect.back()) { - //Determine the interval for value_rij - for (int i = 3; i < knot_vect_size - 1; ++i) { - if (knot_vect[i] <= r && r < knot_vect[i + 1]) { - start_index = i; - break; - } - } - }*/ - - int knot_affect_start = start_index - 3; - - double rsq = r * r; - double rth = rsq * r; - - // Calculate energy - - ret_val[0] = bspline_bases[knot_affect_start + 3].eval0(rth, rsq, r); - ret_val[0] += bspline_bases[knot_affect_start + 2].eval1(rth, rsq, r); - ret_val[0] += bspline_bases[knot_affect_start + 1].eval2(rth, rsq, r); - ret_val[0] += bspline_bases[knot_affect_start].eval3(rth, rsq, r); - - // Calculate force - - ret_val[1] = dnbspline_bases[knot_affect_start + 2].eval0(rsq, r); - ret_val[1] += dnbspline_bases[knot_affect_start + 1].eval1(rsq, r); - ret_val[1] += dnbspline_bases[knot_affect_start].eval2(rsq, r); - - return ret_val; -} - -double uf3_pair_bspline::memory_usage() -{ - double bytes = 0; - - bytes += (double)2*sizeof(int); //knot_vect_size, - //coeff_vect_size - bytes += (double)knot_vect.size()*sizeof(double); //knot_vect - bytes += (double)dnknot_vect.size()*sizeof(double); //dnknot_vect - bytes += (double)coeff_vect.size()*sizeof(double); //coeff_vect - bytes += (double)dncoeff_vect.size()*sizeof(double); //dncoeff_vect - - for (int i = 0; i < knot_vect.size() - 4; i++) - bytes += (double)bspline_bases[i].memory_usage(); //bspline_basis3 - - for (int i = 0; i < dnknot_vect.size() - 3; i++) - bytes += (double)dnbspline_bases[i].memory_usage(); //bspline_basis2 - - bytes += (double)2*sizeof(double); //ret_val - - return bytes; -} diff --git a/src/ML-UF3/uf3_pair_bspline.h b/src/ML-UF3/uf3_pair_bspline.h deleted file mode 100644 index 3f4fe696ef..0000000000 --- a/src/ML-UF3/uf3_pair_bspline.h +++ /dev/null @@ -1,59 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - LAMMPS development team: developers@lammps.org - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "pointers.h" - -#include "uf3_bspline_basis2.h" -#include "uf3_bspline_basis3.h" - -#include - -#ifndef UF3_PAIR_BSPLINE_H -#define UF3_PAIR_BSPLINE_H - -namespace LAMMPS_NS { - -class uf3_pair_bspline { - private: - int knot_vect_size, coeff_vect_size; - std::vector knot_vect, dnknot_vect; - std::vector coeff_vect, dncoeff_vect; - std::vector bspline_bases; - std::vector dnbspline_bases; - int get_starting_index_uniform(double), get_starting_index_nonuniform(double); - int (uf3_pair_bspline::*get_starting_index)(double); - //double knot_spacing=0; - LAMMPS *lmp; - - public: - // dummy constructor - uf3_pair_bspline(); - uf3_pair_bspline(LAMMPS *ulmp, const std::vector &uknot_vect, - const std::vector &ucoeff_vect, - const int &uknot_spacing_type); - - uf3_pair_bspline(LAMMPS *ulmp, const double* uknot_array, - const int uknot_array_size, - const double* ucoeff_array, - const int ucoeff_array_size, - const int uknot_spacing_type); - - ~uf3_pair_bspline(); - int knot_spacing_type; - double knot_spacing=0; - double ret_val[2]; - double *eval(double value_rij); - double memory_usage(); -}; -} // namespace LAMMPS_NS -#endif diff --git a/src/ML-UF3/uf3_triplet_bspline.cpp b/src/ML-UF3/uf3_triplet_bspline.cpp deleted file mode 100644 index 0527599965..0000000000 --- a/src/ML-UF3/uf3_triplet_bspline.cpp +++ /dev/null @@ -1,508 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - LAMMPS development team: developers@lammps.org - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#include "uf3_triplet_bspline.h" -#include "error.h" - -#include - -using namespace LAMMPS_NS; - -// Dummy constructor -uf3_triplet_bspline::uf3_triplet_bspline(){}; - -// Construct a new 3D B-Spline -uf3_triplet_bspline::uf3_triplet_bspline( - LAMMPS *ulmp, const std::vector> &uknot_matrix, - const std::vector>> &ucoeff_matrix, - const int &uknot_spacing_type) -{ - lmp = ulmp; - knot_matrix = uknot_matrix; - coeff_matrix = ucoeff_matrix; - - knot_spacing_type = uknot_spacing_type; - if (knot_spacing_type==0){ - knot_spacing_ij = knot_matrix[2][4]-knot_matrix[2][3]; - knot_spacing_ik = knot_matrix[1][4]-knot_matrix[1][3]; - knot_spacing_jk = knot_matrix[0][4]-knot_matrix[0][3]; - get_starting_index=&uf3_triplet_bspline::get_starting_index_uniform; - } - else if (knot_spacing_type==1){ - knot_spacing_ij = 0; - knot_spacing_ik = 0; - knot_spacing_jk = 0; - get_starting_index=&uf3_triplet_bspline::get_starting_index_nonuniform; - } - - else - lmp->error->all(FLERR, "UF3: Expected either '0'(uniform-knots) or \n\ - '1'(non-uniform knots)"); - - knot_vect_size_ij = knot_matrix[2].size(); - knot_vect_size_ik = knot_matrix[1].size(); - knot_vect_size_jk = knot_matrix[0].size(); - - int resolution_ij = knot_vect_size_ij - 4; - int resolution_ik = knot_vect_size_ik - 4; - int resolution_jk = knot_vect_size_jk - 4; - - // Cache Spline Basis Functions - for (int l = 0; l < resolution_ij; l++) { - bsplines_ij.push_back(uf3_bspline_basis3(lmp, &knot_matrix[2][l], 1)); - } - - for (int l = 0; l < resolution_ik; l++) { - // Reuse jk Basis if Knots match - if (knot_matrix[1][l] == knot_matrix[2][l] && knot_matrix[1][l + 1] == knot_matrix[2][l + 1] && - knot_matrix[1][l + 2] == knot_matrix[2][l + 2] && - knot_matrix[1][l + 3] == knot_matrix[2][l + 3]) - bsplines_ik.push_back(bsplines_ij[l]); - else - bsplines_ik.push_back(uf3_bspline_basis3(lmp, &knot_matrix[1][l], 1)); - } - - for (int l = 0; l < resolution_jk; l++) { - bsplines_jk.push_back(uf3_bspline_basis3(lmp, &knot_matrix[0][l], 1)); - } - - // Initialize Coefficients for Derivatives - for (int i = 0; i < coeff_matrix.size(); i++) { - std::vector> dncoeff_vect2; - for (int j = 0; j < coeff_matrix[0].size(); j++) { - std::vector dncoeff_vect; - for (int k = 0; k < coeff_matrix[0][0].size() - 1; k++) { - double dntemp4 = 3 / (knot_matrix[0][k + 4] - knot_matrix[0][k + 1]); - dncoeff_vect.push_back((coeff_matrix[i][j][k + 1] - coeff_matrix[i][j][k]) * dntemp4); - } - dncoeff_vect2.push_back(dncoeff_vect); - } - dncoeff_matrix_jk.push_back(dncoeff_vect2); - } - - for (int i = 0; i < coeff_matrix.size(); i++) { - std::vector> dncoeff_vect2; - for (int j = 0; j < coeff_matrix[0].size() - 1; j++) { - double dntemp4 = 3 / (knot_matrix[1][j + 4] - knot_matrix[1][j + 1]); - std::vector dncoeff_vect; - for (int k = 0; k < coeff_matrix[0][0].size(); k++) { - dncoeff_vect.push_back((coeff_matrix[i][j + 1][k] - coeff_matrix[i][j][k]) * dntemp4); - } - dncoeff_vect2.push_back(dncoeff_vect); - } - dncoeff_matrix_ik.push_back(dncoeff_vect2); - } - - for (int i = 0; i < coeff_matrix.size() - 1; i++) { - std::vector> dncoeff_vect2; - double dntemp4 = 3 / (knot_matrix[2][i + 4] - knot_matrix[2][i + 1]); - for (int j = 0; j < coeff_matrix[0].size(); j++) { - std::vector dncoeff_vect; - for (int k = 0; k < coeff_matrix[0][0].size(); k++) { - dncoeff_vect.push_back((coeff_matrix[i + 1][j][k] - coeff_matrix[i][j][k]) * dntemp4); - } - dncoeff_vect2.push_back(dncoeff_vect); - } - dncoeff_matrix_ij.push_back(dncoeff_vect2); - } - - std::vector> dnknot_matrix; - for (int i = 0; i < knot_matrix.size(); i++) { - std::vector dnknot_vect; - for (int j = 1; j < knot_matrix[0].size() - 1; j++) { - dnknot_vect.push_back(knot_matrix[i][j]); - } - dnknot_matrix.push_back(dnknot_vect); - } - - // Cache Derivative Spline Basis Functions - for (int l = 0; l < resolution_ij - 1; l++) { - dnbsplines_ij.push_back(uf3_bspline_basis2(lmp, &dnknot_matrix[2][l], 1)); - } - - for (int l = 0; l < resolution_ik - 1; l++) { - // Reuse jk Basis if Knots match - if (dnknot_matrix[1][l] == dnknot_matrix[2][l] && - dnknot_matrix[1][l + 1] == dnknot_matrix[2][l + 1] && - dnknot_matrix[1][l + 2] == dnknot_matrix[2][l + 2]) - dnbsplines_ik.push_back(dnbsplines_ij[l]); - else - dnbsplines_ik.push_back(uf3_bspline_basis2(lmp, &dnknot_matrix[1][l], 1)); - } - - for (int l = 0; l < resolution_jk - 1; l++) { - dnbsplines_jk.push_back(uf3_bspline_basis2(lmp, &dnknot_matrix[0][l], 1)); - } -} - -// Construct a new 3D B-Spline from arrays -uf3_triplet_bspline::uf3_triplet_bspline( - LAMMPS *ulmp, double **uknot_array, const int *uknot_array_size, - double ***ucoeff_array, const int *ucoeff_array_size, - const int &uknot_spacing_type) -{ - lmp = ulmp; - - knot_matrix.resize(3); - //utils::logmesg(lmp, "knot_matrix dim = {} {} {}\nknots = ",uknot_array_size[0], - // uknot_array_size[1], uknot_array_size[2]); - for (int i = 0; i < 3; i++) { - knot_matrix[i].resize(uknot_array_size[i]); - //utils::logmesg(lmp, "{}= ",i); - for (int j = 0; j < uknot_array_size[i]; j++) { - //utils::logmesg(lmp, "{} ", uknot_array[i][j]); - knot_matrix[i][j] = uknot_array[i][j]; - } - //utils::logmesg(lmp,"\n"); - } - - coeff_matrix.resize(ucoeff_array_size[0]); - for (int i = 0; i < ucoeff_array_size[0]; i++) { - coeff_matrix[i].resize(ucoeff_array_size[1]); - for (int j = 0; j < ucoeff_array_size[1]; j++) { - coeff_matrix[i][j].resize(ucoeff_array_size[2]); - for (int k = 0; k < ucoeff_array_size[2]; k++){ - coeff_matrix[i][j][k] = ucoeff_array[i][j][k]; - } - } - } - - knot_spacing_type = uknot_spacing_type; - if (knot_spacing_type==0){ - knot_spacing_ij = knot_matrix[2][4]-knot_matrix[2][3]; - knot_spacing_ik = knot_matrix[1][4]-knot_matrix[1][3]; - knot_spacing_jk = knot_matrix[0][4]-knot_matrix[0][3]; - get_starting_index=&uf3_triplet_bspline::get_starting_index_uniform; - } - else if (knot_spacing_type==1){ - knot_spacing_ij = 0; - knot_spacing_ik = 0; - knot_spacing_jk = 0; - get_starting_index=&uf3_triplet_bspline::get_starting_index_nonuniform; - } - - else - lmp->error->all(FLERR, "UF3: Expected either '0'(uniform-knots) or \n\ - '1'(non-uniform knots)"); - - knot_vect_size_ij = knot_matrix[2].size(); - knot_vect_size_ik = knot_matrix[1].size(); - knot_vect_size_jk = knot_matrix[0].size(); - - int resolution_ij = knot_vect_size_ij - 4; - int resolution_ik = knot_vect_size_ik - 4; - int resolution_jk = knot_vect_size_jk - 4; - - // Cache Spline Basis Functions - for (int l = 0; l < resolution_ij; l++) { - bsplines_ij.push_back(uf3_bspline_basis3(lmp, &knot_matrix[2][l], 1)); - } - - for (int l = 0; l < resolution_ik; l++) { - // Reuse jk Basis if Knots match - if (knot_matrix[1][l] == knot_matrix[2][l] && knot_matrix[1][l + 1] == knot_matrix[2][l + 1] && - knot_matrix[1][l + 2] == knot_matrix[2][l + 2] && - knot_matrix[1][l + 3] == knot_matrix[2][l + 3]) - bsplines_ik.push_back(bsplines_ij[l]); - else - bsplines_ik.push_back(uf3_bspline_basis3(lmp, &knot_matrix[1][l], 1)); - } - - for (int l = 0; l < resolution_jk; l++) { - bsplines_jk.push_back(uf3_bspline_basis3(lmp, &knot_matrix[0][l], 1)); - } - - // Initialize Coefficients for Derivatives - for (int i = 0; i < coeff_matrix.size(); i++) { - std::vector> dncoeff_vect2; - for (int j = 0; j < coeff_matrix[0].size(); j++) { - std::vector dncoeff_vect; - for (int k = 0; k < coeff_matrix[0][0].size() - 1; k++) { - double dntemp4 = 3 / (knot_matrix[0][k + 4] - knot_matrix[0][k + 1]); - dncoeff_vect.push_back((coeff_matrix[i][j][k + 1] - coeff_matrix[i][j][k]) * dntemp4); - } - dncoeff_vect2.push_back(dncoeff_vect); - } - dncoeff_matrix_jk.push_back(dncoeff_vect2); - } - - for (int i = 0; i < coeff_matrix.size(); i++) { - std::vector> dncoeff_vect2; - for (int j = 0; j < coeff_matrix[0].size() - 1; j++) { - double dntemp4 = 3 / (knot_matrix[1][j + 4] - knot_matrix[1][j + 1]); - std::vector dncoeff_vect; - for (int k = 0; k < coeff_matrix[0][0].size(); k++) { - dncoeff_vect.push_back((coeff_matrix[i][j + 1][k] - coeff_matrix[i][j][k]) * dntemp4); - } - dncoeff_vect2.push_back(dncoeff_vect); - } - dncoeff_matrix_ik.push_back(dncoeff_vect2); - } - - for (int i = 0; i < coeff_matrix.size() - 1; i++) { - std::vector> dncoeff_vect2; - double dntemp4 = 3 / (knot_matrix[2][i + 4] - knot_matrix[2][i + 1]); - for (int j = 0; j < coeff_matrix[0].size(); j++) { - std::vector dncoeff_vect; - for (int k = 0; k < coeff_matrix[0][0].size(); k++) { - dncoeff_vect.push_back((coeff_matrix[i + 1][j][k] - coeff_matrix[i][j][k]) * dntemp4); - } - dncoeff_vect2.push_back(dncoeff_vect); - } - dncoeff_matrix_ij.push_back(dncoeff_vect2); - } - - std::vector> dnknot_matrix; - for (int i = 0; i < knot_matrix.size(); i++) { - std::vector dnknot_vect; - for (int j = 1; j < knot_matrix[0].size() - 1; j++) { - dnknot_vect.push_back(knot_matrix[i][j]); - } - dnknot_matrix.push_back(dnknot_vect); - } - - // Cache Derivative Spline Basis Functions - for (int l = 0; l < resolution_ij - 1; l++) { - dnbsplines_ij.push_back(uf3_bspline_basis2(lmp, &dnknot_matrix[2][l], 1)); - } - - for (int l = 0; l < resolution_ik - 1; l++) { - // Reuse jk Basis if Knots match - if (dnknot_matrix[1][l] == dnknot_matrix[2][l] && - dnknot_matrix[1][l + 1] == dnknot_matrix[2][l + 1] && - dnknot_matrix[1][l + 2] == dnknot_matrix[2][l + 2]) - dnbsplines_ik.push_back(dnbsplines_ij[l]); - else - dnbsplines_ik.push_back(uf3_bspline_basis2(lmp, &dnknot_matrix[1][l], 1)); - } - - for (int l = 0; l < resolution_jk - 1; l++) { - dnbsplines_jk.push_back(uf3_bspline_basis2(lmp, &dnknot_matrix[0][l], 1)); - } -} - - -// Destructor -uf3_triplet_bspline::~uf3_triplet_bspline() {} - -// Evaluate 3D B-Spline value -double *uf3_triplet_bspline::eval(double value_rij, double value_rik, double value_rjk) -{ - - // Find starting knots - - //int iknot_ij = starting_knot(knot_matrix[2], knot_vect_size_ij, value_rij) - 3; - //int iknot_ik = starting_knot(knot_matrix[1], knot_vect_size_ik, value_rik) - 3; - //int iknot_jk = starting_knot(knot_matrix[0], knot_vect_size_jk, value_rjk) - 3; - int iknot_ij = (this->*get_starting_index)(knot_matrix[2], knot_vect_size_ij, value_rij,knot_spacing_ij) - 3; - int iknot_ik = (this->*get_starting_index)(knot_matrix[1], knot_vect_size_ik, value_rik,knot_spacing_ik) - 3; - int iknot_jk = (this->*get_starting_index)(knot_matrix[0], knot_vect_size_jk, value_rjk,knot_spacing_jk) - 3; - - double rsq_ij = value_rij * value_rij; - double rsq_ik = value_rik * value_rik; - double rsq_jk = value_rjk * value_rjk; - double rth_ij = rsq_ij * value_rij; - double rth_ik = rsq_ik * value_rik; - double rth_jk = rsq_jk * value_rjk; - - // Calculate energies - - double basis_ij[4]; - basis_ij[0] = bsplines_ij[iknot_ij].eval3(rth_ij, rsq_ij, value_rij); - basis_ij[1] = bsplines_ij[iknot_ij + 1].eval2(rth_ij, rsq_ij, value_rij); - basis_ij[2] = bsplines_ij[iknot_ij + 2].eval1(rth_ij, rsq_ij, value_rij); - basis_ij[3] = bsplines_ij[iknot_ij + 3].eval0(rth_ij, rsq_ij, value_rij); - - double basis_ik[4]; - basis_ik[0] = bsplines_ik[iknot_ik].eval3(rth_ik, rsq_ik, value_rik); - basis_ik[1] = bsplines_ik[iknot_ik + 1].eval2(rth_ik, rsq_ik, value_rik); - basis_ik[2] = bsplines_ik[iknot_ik + 2].eval1(rth_ik, rsq_ik, value_rik); - basis_ik[3] = bsplines_ik[iknot_ik + 3].eval0(rth_ik, rsq_ik, value_rik); - - double basis_jk[4]; - basis_jk[0] = bsplines_jk[iknot_jk].eval3(rth_jk, rsq_jk, value_rjk); - basis_jk[1] = bsplines_jk[iknot_jk + 1].eval2(rth_jk, rsq_jk, value_rjk); - basis_jk[2] = bsplines_jk[iknot_jk + 2].eval1(rth_jk, rsq_jk, value_rjk); - basis_jk[3] = bsplines_jk[iknot_jk + 3].eval0(rth_jk, rsq_jk, value_rjk); - - ret_val[0] = 0; - ret_val[1] = 0; - ret_val[2] = 0; - ret_val[3] = 0; - - for (int i = 0; i < 4; i++) { - const double basis_iji = basis_ij[i]; // prevent repeated access of same memory location - for (int j = 0; j < 4; j++) { - const double factor = basis_iji * basis_ik[j]; // prevent repeated access of same memory location - const double* slice = &coeff_matrix[i + iknot_ij][j + iknot_ik][iknot_jk]; // declare a contigues 1D slice of memory - double tmp[4]; // declare tmp array that holds the 4 tmp values so the can be computed simultaniously in 4 separate registeres. - tmp[0] = slice[0] * basis_jk[0]; - tmp[1] = slice[1] * basis_jk[1]; - tmp[2] = slice[2] * basis_jk[2]; - tmp[3] = slice[3] * basis_jk[3]; - double sum = tmp[0] + tmp[1] + tmp[2] + tmp[3]; - ret_val[0] += factor * sum; // use 1 fused multiply-add (FMA) - } - } - - // Calculate forces - - double dnbasis_ij[4]; - dnbasis_ij[0] = dnbsplines_ij[iknot_ij].eval2(rsq_ij, value_rij); - dnbasis_ij[1] = dnbsplines_ij[iknot_ij + 1].eval1(rsq_ij, value_rij); - dnbasis_ij[2] = dnbsplines_ij[iknot_ij + 2].eval0(rsq_ij, value_rij); - dnbasis_ij[3] = 0; - - double dnbasis_ik[4]; - dnbasis_ik[0] = dnbsplines_ik[iknot_ik].eval2(rsq_ik, value_rik); - dnbasis_ik[1] = dnbsplines_ik[iknot_ik + 1].eval1(rsq_ik, value_rik); - dnbasis_ik[2] = dnbsplines_ik[iknot_ik + 2].eval0(rsq_ik, value_rik); - dnbasis_ik[3] = 0; - - double dnbasis_jk[4]; - dnbasis_jk[0] = dnbsplines_jk[iknot_jk].eval2(rsq_jk, value_rjk); - dnbasis_jk[1] = dnbsplines_jk[iknot_jk + 1].eval1(rsq_jk, value_rjk); - dnbasis_jk[2] = dnbsplines_jk[iknot_jk + 2].eval0(rsq_jk, value_rjk); - dnbasis_jk[3] = 0; - - for (int i = 0; i < 3; i++) { - const double dnbasis_iji = dnbasis_ij[i]; - for (int j = 0; j < 4; j++) { - const double factor = dnbasis_iji * basis_ik[j]; - const double* slice = &dncoeff_matrix_ij[iknot_ij + i][iknot_ik + j][iknot_jk]; - double tmp[4]; - tmp[0] = slice[0] * basis_jk[0]; - tmp[1] = slice[1] * basis_jk[1]; - tmp[2] = slice[2] * basis_jk[2]; - tmp[3] = slice[3] * basis_jk[3]; - double sum = tmp[0] + tmp[1] + tmp[2] + tmp[3]; - ret_val[1] += factor * sum; - } - } - - for (int i = 0; i < 4; i++) { - const double basis_iji = basis_ij[i]; - for (int j = 0; j < 3; j++) { - const double factor = basis_iji * dnbasis_ik[j]; - const double* slice = &dncoeff_matrix_ik[iknot_ij + i][iknot_ik + j][iknot_jk]; - double tmp[4]; - tmp[0] = slice[0] * basis_jk[0]; - tmp[1] = slice[1] * basis_jk[1]; - tmp[2] = slice[2] * basis_jk[2]; - tmp[3] = slice[3] * basis_jk[3]; - double sum = tmp[0] + tmp[1] + tmp[2] + tmp[3]; - ret_val[2] += factor * sum; - } - } - - for (int i = 0; i < 4; i++) { - const double basis_iji = basis_ij[i]; - for (int j = 0; j < 4; j++) { - const double factor = basis_iji * basis_ik[j]; - const double* slice = &dncoeff_matrix_jk[iknot_ij + i][iknot_ik + j][iknot_jk]; - double tmp[3]; - tmp[0] = slice[0] * dnbasis_jk[0]; - tmp[1] = slice[1] * dnbasis_jk[1]; - tmp[2] = slice[2] * dnbasis_jk[2]; - double sum = tmp[0] + tmp[1] + tmp[2]; - ret_val[3] += factor * sum; - } - } - - return ret_val; -} - -// Find starting knot for spline evaluation - -int uf3_triplet_bspline::starting_knot(const std::vector knot_vect, int knot_vect_size, - double r) -{ - if (knot_vect.front() <= r && r < knot_vect.back()) { - for (int i = 3; i < knot_vect_size - 1; i++) { - if (knot_vect[i] <= r && r < knot_vect[i + 1]) return i; - } - } - - return 0; -} - -int uf3_triplet_bspline::get_starting_index_uniform(const std::vector knot_vect, int knot_vect_size, - double r, double knot_spacing) -{ - return 3+(int)((r-knot_vect[0])/knot_spacing); -} - -int uf3_triplet_bspline::get_starting_index_nonuniform(const std::vector knot_vect, int knot_vect_size, - double r, double knot_spacing) -{ - if (knot_vect.front() <= r && r < knot_vect.back()) { - //Determine the interval for value_rij - for (int i = 3; i < knot_vect_size - 1; ++i) { - if (knot_vect[i] <= r && r < knot_vect[i + 1]) { - return i; - } - } - } - return -1; -} - -double uf3_triplet_bspline::memory_usage() -{ - double bytes = 0; - - bytes += (double) 3*sizeof(int); //knot_vect_size_ij, - //knot_vect_size_ik, - //knot_vect_size_jk; - - for (int i=0; i - -#ifndef UF3_TRIPLET_BSPLINE_H -#define UF3_TRIPLET_BSPLINE_H - -namespace LAMMPS_NS { -class uf3_triplet_bspline { - private: - LAMMPS *lmp; - int knot_vect_size_ij, knot_vect_size_ik, knot_vect_size_jk; - std::vector>> coeff_matrix, dncoeff_matrix_ij, dncoeff_matrix_ik, - dncoeff_matrix_jk; - std::vector> knot_matrix; - std::vector bsplines_ij, bsplines_ik, bsplines_jk; - std::vector dnbsplines_ij, dnbsplines_ik, dnbsplines_jk; - int get_starting_index_uniform(const std::vector, int, double, double); - int get_starting_index_nonuniform(const std::vector, int, double, double); - int (uf3_triplet_bspline::*get_starting_index)(const std::vector, int, double, double); - //double knot_spacing_ij=0,knot_spacing_ik=0,knot_spacing_jk=0; - //double _alignvar(, 8) ret_val[4]; // Force memory alignment on 8 byte boundaries - double ret_val[4]; - - int starting_knot(const std::vector, int, double); - - public: - //Dummy Constructor - uf3_triplet_bspline(); - uf3_triplet_bspline(LAMMPS *ulmp, const std::vector> &uknot_matrix, - const std::vector>> &ucoeff_matrix, - const int &uknot_spacing_type); - - uf3_triplet_bspline(LAMMPS *ulmp, double **uknot_array, const int *uknot_array_size, - double ***ucoeff_array, const int *ucoeff_array_size, - const int &uknot_spacing_type); - - ~uf3_triplet_bspline(); - int knot_spacing_type; - double knot_spacing_ij=0,knot_spacing_ik=0,knot_spacing_jk=0; - double *eval(double value_rij, double value_rik, double value_rjk); - - double memory_usage(); -}; -} // namespace LAMMPS_NS -#endif