Skip to content

Commit

Permalink
Merge pull request #8 from mir-group/multiple-cutoffs
Browse files Browse the repository at this point in the history
Multiple cutoffs
  • Loading branch information
jonpvandermause authored Feb 21, 2021
2 parents f63f7f5 + 44ee981 commit 8f877b7
Show file tree
Hide file tree
Showing 8 changed files with 388 additions and 17 deletions.
2 changes: 1 addition & 1 deletion flare_pp/sparse_gp.py
Original file line number Diff line number Diff line change
Expand Up @@ -259,7 +259,7 @@ def update_db(
custom_range=(),
energy: float = None,
stress: "ndarray" = None,
mode: str = "all",
mode: str = "specific",
sgp: SparseGP = None, # for creating sgp_var
update_qr=True,
):
Expand Down
102 changes: 102 additions & 0 deletions lammps/lammps_descriptor.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,108 @@
#include <cmath>
#include <iostream>

void single_bond_multiple_cutoffs(
double **x, int *type, int jnum, int n_inner, int i, double xtmp,
double ytmp, double ztmp, int *jlist,
std::function<void(std::vector<double> &, std::vector<double> &, double,
int, std::vector<double>)>
basis_function,
std::function<void(std::vector<double> &, double, double,
std::vector<double>)>
cutoff_function,
int n_species, int N, int lmax,
const std::vector<double> &radial_hyps,
const std::vector<double> &cutoff_hyps, Eigen::VectorXd &single_bond_vals,
Eigen::MatrixXd &single_bond_env_dervs,
const Eigen::MatrixXd &cutoff_matrix) {

// Initialize basis vectors and spherical harmonics.
std::vector<double> g = std::vector<double>(N, 0);
std::vector<double> gx = std::vector<double>(N, 0);
std::vector<double> gy = std::vector<double>(N, 0);
std::vector<double> gz = std::vector<double>(N, 0);

int n_harmonics = (lmax + 1) * (lmax + 1);
std::vector<double> h = std::vector<double>(n_harmonics, 0);
std::vector<double> hx = std::vector<double>(n_harmonics, 0);
std::vector<double> hy = std::vector<double>(n_harmonics, 0);
std::vector<double> hz = std::vector<double>(n_harmonics, 0);

// Prepare LAMMPS variables.
int central_species = type[i] - 1;
double delx, dely, delz, rsq, r, bond, bond_x, bond_y, bond_z, g_val, gx_val,
gy_val, gz_val, h_val;
int j, s, descriptor_counter;

// Initialize vectors.
int n_radial = n_species * N;
int n_bond = n_radial * n_harmonics;
single_bond_vals = Eigen::VectorXd::Zero(n_bond);
single_bond_env_dervs = Eigen::MatrixXd::Zero(n_inner * 3, n_bond);

// Initialize radial hyperparameters.
std::vector<double> new_radial_hyps = radial_hyps;

// Loop over neighbors.
int n_count = 0;
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];

delx = x[j][0] - xtmp;
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = delx * delx + dely * dely + delz * delz;
r = sqrt(rsq);

// Retrieve the cutoff.
int s = type[j] - 1;
double cutoff = cutoff_matrix(central_species, s);
double cutforcesq = cutoff * cutoff;

if (rsq < cutforcesq) { // minus a small value to prevent numerial error
// Reset endpoint of the radial basis set.
new_radial_hyps[1] = cutoff;

calculate_radial(g, gx, gy, gz, basis_function, cutoff_function, delx,
dely, delz, r, cutoff, N, new_radial_hyps, cutoff_hyps);
get_Y(h, hx, hy, hz, delx, dely, delz, lmax);

// Store the products and their derivatives.
descriptor_counter = s * N * n_harmonics;

for (int radial_counter = 0; radial_counter < N; radial_counter++) {
// Retrieve radial values.
g_val = g[radial_counter];
gx_val = gx[radial_counter];
gy_val = gy[radial_counter];
gz_val = gz[radial_counter];

for (int angular_counter = 0; angular_counter < n_harmonics;
angular_counter++) {

h_val = h[angular_counter];
bond = g_val * h_val;

// Calculate derivatives with the product rule.
bond_x = gx_val * h_val + g_val * hx[angular_counter];
bond_y = gy_val * h_val + g_val * hy[angular_counter];
bond_z = gz_val * h_val + g_val * hz[angular_counter];

// Update single bond basis arrays.
single_bond_vals(descriptor_counter) += bond;

single_bond_env_dervs(n_count * 3, descriptor_counter) += bond_x;
single_bond_env_dervs(n_count * 3 + 1, descriptor_counter) += bond_y;
single_bond_env_dervs(n_count * 3 + 2, descriptor_counter) += bond_z;

descriptor_counter++;
}
}
n_count++;
}
}
}

void single_bond(
double **x, int *type, int jnum, int n_inner, int i, double xtmp,
double ytmp, double ztmp, int *jlist,
Expand Down
15 changes: 15 additions & 0 deletions lammps/lammps_descriptor.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,21 @@ void single_bond(
const std::vector<double> &cutoff_hyps, Eigen::VectorXd &single_bond_vals,
Eigen::MatrixXd &single_bond_env_dervs);

void single_bond_multiple_cutoffs(
double **x, int *type, int jnum, int n_inner, int i, double xtmp,
double ytmp, double ztmp, int *jlist,
std::function<void(std::vector<double> &, std::vector<double> &, double,
int, std::vector<double>)>
basis_function,
std::function<void(std::vector<double> &, double, double,
std::vector<double>)>
cutoff_function,
int n_species, int N, int lmax,
const std::vector<double> &radial_hyps,
const std::vector<double> &cutoff_hyps, Eigen::VectorXd &single_bond_vals,
Eigen::MatrixXd &single_bond_env_dervs,
const Eigen::MatrixXd &cutoff_matrix);

void B2_descriptor(Eigen::VectorXd &B2_vals, Eigen::MatrixXd &B2_env_dervs,
double &norm_squared, Eigen::VectorXd &B2_env_dot,
const Eigen::VectorXd &single_bond_vals,
Expand Down
40 changes: 32 additions & 8 deletions lammps/pair_flare.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -92,19 +92,23 @@ void PairFLARE::compute(int eflag, int vflag) {
n_inner = 0;
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
int s = type[j] - 1;
double cutoff_val = cutoff_matrix(itype-1, s);

delx = x[j][0] - xtmp;
dely = x[j][1] - ytmp;
delz = x[j][2] - ztmp;
rsq = delx * delx + dely * dely + delz * delz;
if (rsq < (cutoff * cutoff))
if (rsq < (cutoff_val * cutoff_val))
n_inner++;
}

// Compute covariant descriptors.
single_bond(x, type, jnum, n_inner, i, xtmp, ytmp, ztmp, jlist,
basis_function, cutoff_function, cutoff, n_species, n_max,
l_max, radial_hyps, cutoff_hyps, single_bond_vals,
single_bond_env_dervs);
single_bond_multiple_cutoffs(x, type, jnum, n_inner, i, xtmp, ytmp, ztmp,
jlist, basis_function, cutoff_function,
n_species, n_max, l_max, radial_hyps,
cutoff_hyps, single_bond_vals,
single_bond_env_dervs, cutoff_matrix);

// Compute invariant descriptors.
B2_descriptor(B2_vals, B2_env_dervs, B2_norm_squared, B2_env_dot,
Expand All @@ -125,12 +129,14 @@ void PairFLARE::compute(int eflag, int vflag) {
n_count = 0;
for (int jj = 0; jj < jnum; jj++) {
j = jlist[jj];
int s = type[j] - 1;
double cutoff_val = cutoff_matrix(itype-1, s);
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;

if (rsq < (cutoff * cutoff)) {
if (rsq < (cutoff_val * cutoff_val)) {
double fx = -partial_forces(n_count * 3);
double fy = -partial_forces(n_count * 3 + 1);
double fz = -partial_forces(n_count * 3 + 2);
Expand Down Expand Up @@ -262,8 +268,6 @@ void PairFLARE::read_file(char *filename) {
fgets(line, MAXLINE, fptr);
sscanf(line, "%s", cutoff_string); // Cutoff function
cutoff_string_length = strlen(cutoff_string);
fgets(line, MAXLINE, fptr);
sscanf(line, "%lg", &cutoff); // Cutoff
}

MPI_Bcast(&n_species, 1, MPI_INT, 0, world);
Expand All @@ -276,6 +280,26 @@ void PairFLARE::read_file(char *filename) {
MPI_Bcast(radial_string, radial_string_length + 1, MPI_CHAR, 0, world);
MPI_Bcast(cutoff_string, cutoff_string_length + 1, MPI_CHAR, 0, world);

// Parse the cutoffs.
int n_cutoffs = n_species * n_species;
memory->create(cutoffs, n_cutoffs, "pair:cutoffs");
if (me == 0)
grab(fptr, n_cutoffs, cutoffs);
MPI_Bcast(cutoffs, n_cutoffs, MPI_DOUBLE, 0, world);

// Fill in the cutoff matrix.
cutoff = -1;
cutoff_matrix = Eigen::MatrixXd::Zero(n_species, n_species);
int cutoff_count = 0;
for (int i = 0; i < n_species; i++){
for (int j = 0; j < n_species; j++){
double cutoff_val = cutoffs[cutoff_count];
cutoff_matrix(i, j) = cutoff_val;
if (cutoff_val > cutoff) cutoff = cutoff_val;
cutoff_count ++;
}
}

// Set number of descriptors.
int n_radial = n_max * n_species;
n_descriptors = (n_radial * (n_radial + 1) / 2) * (l_max + 1);
Expand Down
4 changes: 2 additions & 2 deletions lammps/pair_flare.h
Original file line number Diff line number Diff line change
Expand Up @@ -40,8 +40,8 @@ class PairFLARE : public Pair {
std::vector<double> radial_hyps, cutoff_hyps;

double cutoff;
double *beta;
Eigen::MatrixXd beta_matrix;
double *beta, *cutoffs;
Eigen::MatrixXd beta_matrix, cutoff_matrix;
std::vector<Eigen::MatrixXd> beta_matrices;

virtual void allocate();
Expand Down
5 changes: 5 additions & 0 deletions src/ace_binding.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -100,10 +100,15 @@ PYBIND11_MODULE(_C_flare, m) {
.def(py::init<const std::string &, const std::string &,
const std::vector<double> &, const std::vector<double> &,
const std::vector<int> &>())
.def(py::init<const std::string &, const std::string &,
const std::vector<double> &, const std::vector<double> &,
const std::vector<int> &,
const Eigen::MatrixXd &>())
.def_readonly("radial_basis", &B2::radial_basis)
.def_readonly("cutoff_function", &B2::cutoff_function)
.def_readonly("radial_hyps", &B2::radial_hyps)
.def_readonly("cutoff_hyps", &B2::cutoff_hyps)
.def_readonly("cutoffs", &B2::cutoffs)
.def_readonly("descriptor_settings", &B2::descriptor_settings);

py::class_<B2_Simple, Descriptor>(m, "B2_Simple")
Expand Down
Loading

0 comments on commit 8f877b7

Please sign in to comment.