ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added ml-uf3 src files

Browse Source
This commit is contained in:
Ajinkya Hire
2024-03-25 12:07:18 -04:00
parent 810ae3cc45
commit 6dded43b2c
10 changed files with 2479 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1312
src/ML-UF3/pair_uf3.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

88
src/ML-UF3/pair_uf3.h Normal file
View File

@ -0,0 +1,88 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
* Contributing authors: Ajinkya Hire(U of Florida),
* Hendrik Kraß (U of Constance),
* Richard Hennig (U of Florida)
* ---------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(uf3,PairUF3);
// clang-format on
#else
#ifndef LMP_PAIR_UF3_H
#define LMP_PAIR_UF3_H
#include "uf3_pair_bspline.h"
#include "uf3_triplet_bspline.h"
#include "pair.h"
#include <unordered_map>
namespace LAMMPS_NS {
class PairUF3 : public Pair {
public:
PairUF3(class LAMMPS *);
~PairUF3() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
void init_list(int, class NeighList *) override; // needed for ptr to full neigh list
double init_one(int, int) override; // needed for cutoff radius for neighbour list
double single(int, int, int, int, double, double, double, double &) override;
double memory_usage() override;
protected:
void uf3_read_pot_file(char *potf_name);
void uf3_read_pot_file(int i, int j, char *potf_name);
void uf3_read_pot_file(int i, int j, int k, char *potf_name);
int num_of_elements, nbody_flag, n2body_pot_files, n3body_pot_files, tot_pot_files;
int bsplines_created;
int coeff_matrix_dim1, coeff_matrix_dim2, coeff_matrix_dim3, coeff_matrix_elements_len;
bool pot_3b;
int ***setflag_3b, **knot_spacing_type_2b, ***knot_spacing_type_3b;
double **cut, ***cut_3b, **cut_3b_list, ****min_cut_3b;
virtual void allocate();
void create_bsplines();
std::vector<std::vector<std::vector<double>>> n2b_knot, n2b_coeff;
std::vector<std::vector<std::vector<std::vector<std::vector<double>>>>> n3b_knot_matrix;
std::unordered_map<std::string, std::vector<std::vector<std::vector<double>>>> n3b_coeff_matrix;
std::vector<std::vector<uf3_pair_bspline>> UFBS2b;
std::vector<std::vector<std::vector<uf3_triplet_bspline>>> UFBS3b;
int *neighshort, maxshort; // short neighbor list array for 3body interaction
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/

88
src/ML-UF3/uf3_bspline_basis2.cpp Normal file
View File

@ -0,0 +1,88 @@
#include "uf3_bspline_basis2.h"
#include "utils.h"
#include <vector>
using namespace LAMMPS_NS;
// Constructor
// Initializes coefficients and knots
// Requires [knots] to have length 4
uf3_bspline_basis2::uf3_bspline_basis2(LAMMPS *ulmp, const double *knots, double coefficient)
{
lmp = ulmp;
double c0, c1, c2;
c0 = coefficient *
(pow(knots[0], 2) /
(pow(knots[0], 2) - knots[0] * knots[1] - knots[0] * knots[2] + knots[1] * knots[2]));
c1 = coefficient *
(-2 * knots[0] /
(pow(knots[0], 2) - knots[0] * knots[1] - knots[0] * knots[2] + knots[1] * knots[2]));
c2 = coefficient *
(1 / (pow(knots[0], 2) - knots[0] * knots[1] - knots[0] * knots[2] + knots[1] * knots[2]));
constants.push_back(c0);
constants.push_back(c1);
constants.push_back(c2);
c0 = coefficient *
(-knots[1] * knots[3] /
(pow(knots[1], 2) - knots[1] * knots[2] - knots[1] * knots[3] + knots[2] * knots[3]) -
knots[0] * knots[2] /
(knots[0] * knots[1] - knots[0] * knots[2] - knots[1] * knots[2] + pow(knots[2], 2)));
c1 = coefficient *
(knots[1] /
(pow(knots[1], 2) - knots[1] * knots[2] - knots[1] * knots[3] + knots[2] * knots[3]) +
knots[3] /
(pow(knots[1], 2) - knots[1] * knots[2] - knots[1] * knots[3] + knots[2] * knots[3]) +
knots[0] /
(knots[0] * knots[1] - knots[0] * knots[2] - knots[1] * knots[2] + pow(knots[2], 2)) +
knots[2] /
(knots[0] * knots[1] - knots[0] * knots[2] - knots[1] * knots[2] + pow(knots[2], 2)));
c2 = coefficient *
(-1 / (pow(knots[1], 2) - knots[1] * knots[2] - knots[1] * knots[3] + knots[2] * knots[3]) -
1 / (knots[0] * knots[1] - knots[0] * knots[2] - knots[1] * knots[2] + pow(knots[2], 2)));
constants.push_back(c0);
constants.push_back(c1);
constants.push_back(c2);
c0 = coefficient *
(pow(knots[3], 2) /
(knots[1] * knots[2] - knots[1] * knots[3] - knots[2] * knots[3] + pow(knots[3], 2)));
c1 = coefficient *
(-2 * knots[3] /
(knots[1] * knots[2] - knots[1] * knots[3] - knots[2] * knots[3] + pow(knots[3], 2)));
c2 = coefficient *
(1 / (knots[1] * knots[2] - knots[1] * knots[3] - knots[2] * knots[3] + pow(knots[3], 2)));
constants.push_back(c0);
constants.push_back(c1);
constants.push_back(c2);
}
uf3_bspline_basis2::~uf3_bspline_basis2() {}
// Evaluate outer-left part of spline
double uf3_bspline_basis2::eval0(double rsq, double r)
{
return rsq * constants[2] + r * constants[1] + constants[0];
}
// Evaluate center-left part of spline
double uf3_bspline_basis2::eval1(double rsq, double r)
{
return rsq * constants[5] + r * constants[4] + constants[3];
}
// Evaluate center-right part of spline
double uf3_bspline_basis2::eval2(double rsq, double r)
{
return rsq * constants[8] + r * constants[7] + constants[6];
}
double uf3_bspline_basis2::memory_usage()
{
double bytes = 0;
bytes += (double)constants.size()*sizeof(double);
return bytes;
}

32
src/ML-UF3/uf3_bspline_basis2.h Normal file
View File

@ -0,0 +1,32 @@
//De Boor's algorithm @
//https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/de-Boor.html
//For values outside the domain,
//extrapoltaes the left(right) hand side piece of the curve
//Only works for bspline degree upto 3 becuase of definiation of P
//
#include "pointers.h"
#include <vector>
#ifndef UF3_BSPLINE_BASIS2_H
#define UF3_BSPLINE_BASIS2_H
namespace LAMMPS_NS {
class uf3_bspline_basis2 {
private:
LAMMPS *lmp;
std::vector<double> constants;
public:
uf3_bspline_basis2(LAMMPS *ulmp, const double *knots, double coefficient);
~uf3_bspline_basis2();
double eval0(double, double);
double eval1(double, double);
double eval2(double, double);
double memory_usage();
};
} // namespace LAMMPS_NS
#endif

324
src/ML-UF3/uf3_bspline_basis3.cpp Normal file
View File

@ -0,0 +1,324 @@
#include "uf3_bspline_basis3.h"
#include "utils.h"
#include <vector>
using namespace LAMMPS_NS;
// Constructor
// Initializes coefficients and knots
// [knots] needs to have length 4
uf3_bspline_basis3::uf3_bspline_basis3(LAMMPS *ulmp, const double *knots, double coefficient)
{
lmp = ulmp;
double c0, c1, c2, c3;
c0 = coefficient *
(-pow(knots[0], 3) /
(-pow(knots[0], 3) + pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
pow(knots[0], 2) * knots[3] - knots[0] * knots[1] * knots[2] -
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[3]));
c1 = coefficient *
(3 * pow(knots[0], 2) /
(-pow(knots[0], 3) + pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
pow(knots[0], 2) * knots[3] - knots[0] * knots[1] * knots[2] -
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[3]));
c2 = coefficient *
(-3 * knots[0] /
(-pow(knots[0], 3) + pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
pow(knots[0], 2) * knots[3] - knots[0] * knots[1] * knots[2] -
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[3]));
c3 = coefficient *
(1 /
(-pow(knots[0], 3) + pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
pow(knots[0], 2) * knots[3] - knots[0] * knots[1] * knots[2] -
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[3]));
constants.push_back(c0);
constants.push_back(c1);
constants.push_back(c2);
constants.push_back(c3);
c0 = coefficient *
(pow(knots[1], 2) * knots[4] /
(-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
knots[2] * knots[3] * knots[4]) +
pow(knots[0], 2) * knots[2] /
(-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) +
knots[0] * knots[1] * knots[3] /
(-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)));
c1 = coefficient *
(-pow(knots[1], 2) /
(-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
knots[2] * knots[3] * knots[4]) -
2 * knots[1] * knots[4] /
(-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
knots[2] * knots[3] * knots[4]) -
pow(knots[0], 2) /
(-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) -
2 * knots[0] * knots[2] /
(-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) -
knots[0] * knots[1] /
(-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)) -
knots[0] * knots[3] /
(-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)) -
knots[1] * knots[3] /
(-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)));
c2 = coefficient *
(2 * knots[1] /
(-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
knots[2] * knots[3] * knots[4]) +
knots[4] /
(-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
knots[2] * knots[3] * knots[4]) +
2 * knots[0] /
(-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) +
knots[2] /
(-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) +
knots[0] /
(-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)) +
knots[1] /
(-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)) +
knots[3] /
(-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)));
c3 = coefficient *
(-1 /
(-pow(knots[1], 3) + pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
pow(knots[1], 2) * knots[4] - knots[1] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
knots[2] * knots[3] * knots[4]) -
1 /
(-pow(knots[0], 2) * knots[1] + pow(knots[0], 2) * knots[2] +
knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] -
knots[0] * pow(knots[2], 2) - knots[0] * knots[2] * knots[3] -
knots[1] * knots[2] * knots[3] + pow(knots[2], 2) * knots[3]) -
1 /
(-knots[0] * pow(knots[1], 2) + knots[0] * knots[1] * knots[2] +
knots[0] * knots[1] * knots[3] - knots[0] * knots[2] * knots[3] +
pow(knots[1], 2) * knots[3] - knots[1] * knots[2] * knots[3] -
knots[1] * pow(knots[3], 2) + knots[2] * pow(knots[3], 2)));
constants.push_back(c0);
constants.push_back(c1);
constants.push_back(c2);
constants.push_back(c3);
c0 = coefficient *
(-knots[0] * pow(knots[3], 2) /
(-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) -
knots[1] * knots[3] * knots[4] /
(-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) -
knots[2] * pow(knots[4], 2) /
(-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)));
c1 = coefficient *
(2 * knots[0] * knots[3] /
(-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) +
pow(knots[3], 2) /
(-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) +
knots[1] * knots[3] /
(-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) +
knots[1] * knots[4] /
(-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) +
knots[3] * knots[4] /
(-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) +
2 * knots[2] * knots[4] /
(-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)) +
pow(knots[4], 2) /
(-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)));
c2 = coefficient *
(-knots[0] /
(-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) -
2 * knots[3] /
(-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) -
knots[1] /
(-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) -
knots[3] /
(-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) -
knots[4] /
(-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) -
knots[2] /
(-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)) -
2 * knots[4] /
(-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)));
c3 = coefficient *
(1 /
(-knots[0] * knots[1] * knots[2] + knots[0] * knots[1] * knots[3] +
knots[0] * knots[2] * knots[3] - knots[0] * pow(knots[3], 2) +
knots[1] * knots[2] * knots[3] - knots[1] * pow(knots[3], 2) -
knots[2] * pow(knots[3], 2) + pow(knots[3], 3)) +
1 /
(-pow(knots[1], 2) * knots[2] + pow(knots[1], 2) * knots[3] +
knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] -
knots[1] * pow(knots[3], 2) - knots[1] * knots[3] * knots[4] -
knots[2] * knots[3] * knots[4] + pow(knots[3], 2) * knots[4]) +
1 /
(-knots[1] * pow(knots[2], 2) + knots[1] * knots[2] * knots[3] +
knots[1] * knots[2] * knots[4] - knots[1] * knots[3] * knots[4] +
pow(knots[2], 2) * knots[4] - knots[2] * knots[3] * knots[4] -
knots[2] * pow(knots[4], 2) + knots[3] * pow(knots[4], 2)));
constants.push_back(c0);
constants.push_back(c1);
constants.push_back(c2);
constants.push_back(c3);
c0 = coefficient *
(pow(knots[4], 3) /
(-knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] +
knots[1] * knots[3] * knots[4] - knots[1] * pow(knots[4], 2) +
knots[2] * knots[3] * knots[4] - knots[2] * pow(knots[4], 2) - knots[3] * pow(knots[4], 2) +
pow(knots[4], 3)));
c1 = coefficient *
(-3 * pow(knots[4], 2) /
(-knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] +
knots[1] * knots[3] * knots[4] - knots[1] * pow(knots[4], 2) +
knots[2] * knots[3] * knots[4] - knots[2] * pow(knots[4], 2) - knots[3] * pow(knots[4], 2) +
pow(knots[4], 3)));
c2 = coefficient *
(3 * knots[4] /
(-knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] +
knots[1] * knots[3] * knots[4] - knots[1] * pow(knots[4], 2) +
knots[2] * knots[3] * knots[4] - knots[2] * pow(knots[4], 2) - knots[3] * pow(knots[4], 2) +
pow(knots[4], 3)));
c3 = coefficient *
(-1 /
(-knots[1] * knots[2] * knots[3] + knots[1] * knots[2] * knots[4] +
knots[1] * knots[3] * knots[4] - knots[1] * pow(knots[4], 2) +
knots[2] * knots[3] * knots[4] - knots[2] * pow(knots[4], 2) - knots[3] * pow(knots[4], 2) +
pow(knots[4], 3)));
constants.push_back(c0);
constants.push_back(c1);
constants.push_back(c2);
constants.push_back(c3);
}
uf3_bspline_basis3::~uf3_bspline_basis3() {}
// Evaluate outer-left part of spline
double uf3_bspline_basis3::eval0(double rth, double rsq, double r)
{
return rth * constants[3] + rsq * constants[2] + r * constants[1] + constants[0];
}
// Evaluate center-left part of spline
double uf3_bspline_basis3::eval1(double rth, double rsq, double r)
{
return rth * constants[7] + rsq * constants[6] + r * constants[5] + constants[4];
}
// Evaluate center-right part of spline
double uf3_bspline_basis3::eval2(double rth, double rsq, double r)
{
return rth * constants[11] + rsq * constants[10] + r * constants[9] + constants[8];
}
// Evaluate outer-right part of spline
double uf3_bspline_basis3::eval3(double rth, double rsq, double r)
{
return rth * constants[15] + rsq * constants[14] + r * constants[13] + constants[12];
}
double uf3_bspline_basis3::memory_usage()
{
double bytes = 0;
bytes += (double)constants.size()*sizeof(double);
return bytes;
}

33
src/ML-UF3/uf3_bspline_basis3.h Normal file
View File

@ -0,0 +1,33 @@
//De Boor's algorithm @
//https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/de-Boor.html
//For values outside the domain,
//extrapoltaes the left(right) hand side piece of the curve
//Only works for bspline degree upto 3 becuase of definiation of P
//
#include "pointers.h"
#include <vector>
#ifndef UF3_BSPLINE_BASIS3_H
#define UF3_BSPLINE_BASIS3_H
namespace LAMMPS_NS {
class uf3_bspline_basis3 {
private:
LAMMPS *lmp;
std::vector<double> constants;
public:
uf3_bspline_basis3(LAMMPS *ulmp, const double *knots, double coefficient);
~uf3_bspline_basis3();
double eval0(double, double, double);
double eval1(double, double, double);
double eval2(double, double, double);
double eval3(double, double, double);
double memory_usage();
};
} // namespace LAMMPS_NS
#endif

144
src/ML-UF3/uf3_pair_bspline.cpp Normal file
View File

@ -0,0 +1,144 @@
#include "uf3_pair_bspline.h"
#include "uf3_bspline_basis2.h"
#include "uf3_bspline_basis3.h"
#include "utils.h"
#include "error.h"
#include <vector>
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<double> &uknot_vect,
const std::vector<double> &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]));
}
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;
}

55
src/ML-UF3/uf3_pair_bspline.h Normal file
View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
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.
------------------------------------------------------------------------- */
// De Boor's algorithm @
// https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/de-Boor.html
// For values outside the domain, it exhibits undefined behavior.
// Uses fixed B-Spline degree 3.
#include "pointers.h"
#include "uf3_bspline_basis2.h"
#include "uf3_bspline_basis3.h"
#include <vector>
#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<double> knot_vect, dnknot_vect;
std::vector<double> coeff_vect, dncoeff_vect;
std::vector<uf3_bspline_basis3> bspline_bases;
std::vector<uf3_bspline_basis2> 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<double> &uknot_vect,
const std::vector<double> &ucoeff_vect,
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

347
src/ML-UF3/uf3_triplet_bspline.cpp Normal file
View File

@ -0,0 +1,347 @@
#include "uf3_triplet_bspline.h"
#include "error.h"
#include <iostream>
#include <vector>
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<std::vector<double>> &uknot_matrix,
const std::vector<std::vector<std::vector<double>>> &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<std::vector<double>> dncoeff_vect2;
for (int j = 0; j < coeff_matrix[0].size(); j++) {
std::vector<double> 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<std::vector<double>> 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<double> 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<std::vector<double>> 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<double> 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<std::vector<double>> dnknot_matrix;
for (int i = 0; i < knot_matrix.size(); i++) {
std::vector<double> 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<double> 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<double> 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<double> 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;
}
}
}
}
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<coeff_matrix.size(); i++)
for (int j=0; j<coeff_matrix[i].size(); j++)
bytes += (double)coeff_matrix[i][j].size()*sizeof(double);
for (int i=0; i<dncoeff_matrix_ij.size(); i++)
for (int j=0; j<dncoeff_matrix_ij[i].size(); j++)
bytes += (double)dncoeff_matrix_ij[i][j].size()*sizeof(double);
for (int i=0; i<dncoeff_matrix_ik.size(); i++)
for (int j=0; j<dncoeff_matrix_ik[i].size(); j++)
bytes += (double)dncoeff_matrix_ik[i][j].size()*sizeof(double);
for (int i=0; i<dncoeff_matrix_jk.size(); i++)
for (int j=0; j<dncoeff_matrix_jk[i].size(); j++)
bytes += (double)dncoeff_matrix_jk[i][j].size()*sizeof(double);
bytes += (double)knot_matrix[0].size()*sizeof(double);
bytes += (double)knot_matrix[1].size()*sizeof(double);
bytes += (double)knot_matrix[2].size()*sizeof(double);
for (int i=0; i < bsplines_ij.size(); i++)
bytes += (double)bsplines_ij[i].memory_usage();
for (int i=0; i < bsplines_ik.size(); i++)
bytes += (double)bsplines_ik[i].memory_usage();
for (int i=0; i < bsplines_jk.size(); i++)
bytes += (double)bsplines_jk[i].memory_usage();
for (int i=0; i < dnbsplines_ij.size(); i++)
bytes += (double)dnbsplines_ij[i].memory_usage();
for (int i=0; i < dnbsplines_ik.size(); i++)
bytes += (double)dnbsplines_ik[i].memory_usage();
for (int i=0; i < dnbsplines_jk.size(); i++)
bytes += (double)dnbsplines_jk[i].memory_usage();
bytes += (double)4*sizeof(double); //ret_val
return bytes;
}

56
src/ML-UF3/uf3_triplet_bspline.h Normal file
View File

@ -0,0 +1,56 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/ Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
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 "uf3_pair_bspline.h"
#include <vector>
#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<std::vector<std::vector<double>>> coeff_matrix, dncoeff_matrix_ij, dncoeff_matrix_ik,
dncoeff_matrix_jk;
std::vector<std::vector<double>> knot_matrix;
std::vector<uf3_bspline_basis3> bsplines_ij, bsplines_ik, bsplines_jk;
std::vector<uf3_bspline_basis2> dnbsplines_ij, dnbsplines_ik, dnbsplines_jk;
int get_starting_index_uniform(const std::vector<double>, int, double, double);
int get_starting_index_nonuniform(const std::vector<double>, int, double, double);
int (uf3_triplet_bspline::*get_starting_index)(const std::vector<double>, 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<double>, int, double);
public:
//Dummy Constructor
uf3_triplet_bspline();
uf3_triplet_bspline(LAMMPS *ulmp, const std::vector<std::vector<double>> &uknot_matrix,
const std::vector<std::vector<std::vector<double>>> &ucoeff_matrix,
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