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

author information only in cpp, kokkos version needs protected instead of private variables

Browse Source
This commit is contained in:
alphataubio
2024-08-07 00:12:11 -04:00
parent 1166531594
commit ce4e01fb78
8 changed files with 5491 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4576
src/KOKKOS/fix_bond_react_kokkos.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

238
src/KOKKOS/fix_bond_react_kokkos.h Normal file
View File

@ -0,0 +1,238 @@
/* -*- c++ -*- ----------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing Author: Jacob Gissinger (jgissing@stevens.edu)
KOKKOS version (2024/08): Mitch Murphy (alphataubio@gmail.com)
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(bond/react/kk,FixBondReactKokkos<LMPDeviceType>);
FixStyle(bond/react/kk/device,FixBondReactKokkos<LMPDeviceType>);
FixStyle(bond/react/kk/host,FixBondReactKokkos<LMPHostType>);
// clang-format on
#else
// clang-format off
#ifndef LMP_FIX_BOND_REACT_KOKKOS_H
#define LMP_FIX_BOND_REACT_KOKKOS_H
#include "fix_bond_react.h"
#include "kokkos_type.h"
#include <map>
#include <set>
namespace LAMMPS_NS {
template<class DeviceType>
class FixBondReactKokkos : public FixBondReact {
public:
FixBondReactKokkos(class LAMMPS *, int, char **);
~FixBondReactKokkos() override;
//int setmask() override;
void post_constructor() override;
void init() override;
void init_list(int, class NeighList *) override;
void post_integrate() override;
//void post_integrate_respa(int, int) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
double compute_vector(int) override;
//double memory_usage() override;
private:
int *nevery;
FILE *fp;
int *iatomtype, *jatomtype;
int *seed;
double **cutsq, *fraction;
int *max_rxn, *nlocalskips, *nghostlyskips;
int **rate_limit;
int **store_rxn_count;
int *stabilize_steps_flag;
int *custom_charges_fragid;
int *rescale_charges_flag; // if nonzero, indicates number of atoms whose charges are updated
double *mol_total_charge; // sum of charges of post-reaction atoms whose charges are updated
int *create_atoms_flag;
int *modify_create_fragid;
double *overlapsq;
int *molecule_keyword;
int *nconstraints;
char **constraintstr;
std::vector<std::string> rxnfunclist; // lists current special rxn function
std::vector<int> peratomflag; // 1 if special rxn function uses per-atom variable (vs. per-bond)
int **var_flag, **var_id; // for keyword values with variable inputs
int *groupbits;
char **rxn_name; // name of reaction
int *reaction_count;
int *reaction_count_total;
tagint *partner, *finalpartner;
double **distsq;
int *nattempt;
tagint ***attempt;
class Molecule *onemol; // pre-reacted molecule template
class Molecule *twomol; // post-reacted molecule template
Fix *fix1; // nve/limit used to relax reaction sites
Fix *fix2; // properties/atom used to indicate 1) relaxing atoms
// 2) to which 'react' atom belongs
Fix *fix3; // property/atom used for system-wide thermostat
class RanMars **random; // random number for 'prob' keyword
class RanMars **rrhandom; // random number for Arrhenius constraint
class NeighList *list;
class ResetAtomsMol *reset_mol_ids; // class for resetting mol IDs
int *reacted_mol, *unreacted_mol;
int *limit_duration; // indicates how long to relax
char *nve_limit_xmax; // indicates max distance allowed to move when relaxing
char *id_fix1; // id of internally created fix nve/limit
char *id_fix2; // id of internally created fix per-atom properties
char *id_fix3; // id of internally created 'stabilization group' per-atom property fix
char *statted_id; // name of 'stabilization group' per-atom property
char *master_group; // group containing relaxing atoms from all fix rxns
char *exclude_group; // group for system-wide thermostat
void superimpose_algorithm(); // main function of the superimpose algorithm
int *ibonding, *jbonding;
int *closeneigh; // indicates if bonding atoms of a rxn are 1-2, 1-3, or 1-4 neighbors
int nedge, nequivalent, ndelete, ncreate, nchiral; // # edge, equivalent atoms in mapping file
int attempted_rxn; // there was an attempt!
int *local_rxn_count;
int *ghostly_rxn_count;
int avail_guesses; // num of restore points available
int *guess_branch; // used when there is more than two choices when guessing
int **restore_pt; // contains info about restore points
tagint **restore; // contains info about restore points
int *pioneer_count; // counts pioneers
int **edge; // atoms in molecule templates with incorrect valences
int ***equivalences; // relation between pre- and post-reacted templates
int ***reverse_equiv; // re-ordered equivalences
int **landlocked_atoms; // all atoms at least three bonds away from edge atoms
int **custom_charges; // atoms whose charge should be updated
int **delete_atoms; // atoms in pre-reacted templates to delete
int **create_atoms; // atoms in post-reacted templates to create
int ***chiral_atoms; // pre-react chiral atoms. 1) flag 2) orientation 3-4) ordered atom types
int **nxspecial, **onemol_nxspecial, **twomol_nxspecial; // full number of 1-4 neighbors
tagint **xspecial, **onemol_xspecial, **twomol_xspecial; // full 1-4 neighbor list
int pion, neigh, trace; // important indices for various loops. required for restore points
int lcl_inst; // reaction instance
tagint **glove; // 1st colmn: pre-reacted template, 2nd colmn: global IDs
// for all mega_gloves: first row is the ID of bond/react
// 'cuff' leaves room for additional values carried around
int cuff; // default = 1, w/ rescale_charges_flag = 2
double **my_mega_glove; // local + ghostly reaction instances
double **local_mega_glove; // consolidation of local reaction instances
double **ghostly_mega_glove; // consolidation of nonlocal reaction instances
double **global_mega_glove; // consolidation (inter-processor) of gloves
// containing nonlocal atoms
int *localsendlist; // indicates ghosts of other procs
int my_num_mega; // local + ghostly reaction instances (on this proc)
int local_num_mega; // num of local reaction instances
int ghostly_num_mega; // num of ghostly reaction instances
int global_megasize; // num of reaction instances in global_mega_glove
int *pioneers; // during Superimpose Algorithm, atoms which have been assigned,
// but whose first neighbors haven't
int glove_counter; // used to determine when to terminate Superimpose Algorithm
void read_variable_keyword(const char *, int, int);
void read_map_file(int);
void EdgeIDs(char *, int);
void Equivalences(char *, int);
void DeleteAtoms(char *, int);
void CreateAtoms(char *, int);
void CustomCharges(int, int);
void ChiralCenters(char *, int);
void ReadConstraints(char *, int);
void readID(char *, int, int, int);
void make_a_guess();
void neighbor_loop();
void check_a_neighbor();
void crosscheck_the_neighbor();
void inner_crosscheck_loop();
int ring_check();
int check_constraints();
void get_IDcoords(int, int, double *);
double get_temperature(tagint **, int, int);
double get_totalcharge();
void customvarnames(); // get per-atom variables names used by custom constraint
void get_customvars(); // evaluate local values for variables names used by custom constraint
double custom_constraint(const std::string &); // evaulate expression for custom constraint
double rxnfunction(const std::string &, const std::string &,
const std::string &); // eval rxn_sum and rxn_ave
void get_atoms2bond(int);
int get_chirality(double[12]); // get handedness given an ordered set of coordinates
void open(char *);
void readline(char *);
void parse_keyword(int, char *, char *);
void far_partner();
void close_partner();
void get_molxspecials();
void find_landlocked_atoms(int);
void glove_ghostcheck();
void ghost_glovecast();
void update_everything();
int insert_atoms(tagint **, int);
void unlimit_bond(); // removes atoms from stabilization, and other post-reaction every-step operations
void dedup_mega_gloves(int); //dedup global mega_glove
void write_restart(FILE *) override;
void restart(char *buf) override;
// store restart data
struct Set {
int nreacts;
char rxn_name[MAXNAME];
int reaction_count_total;
int max_rate_limit_steps;
};
Set *set;
struct Constraint {
int type;
int id[MAXCONIDS];
int idtype[MAXCONIDS];
double par[MAXCONPAR];
std::string str;
};
int ncustomvars;
std::vector<std::string> customvarstrs;
int nvvec;
double **vvec; // per-atom vector to store custom constraint atom-style variable values
class Compute *cperbond; // pointer to 'compute bond/local' used by custom constraint ('rxnbond' function)
std::map<std::set<tagint>, int> atoms2bond; // maps atom pair to index of local bond array
std::vector<std::vector<Constraint>> constraints;
// DEBUG
void print_bb();
};
} // namespace LAMMPS_NS
#endif
#endif

35
src/KOKKOS/fix_cmap_kokkos.cpp
View File

@ -35,6 +35,16 @@
using namespace LAMMPS_NS;
using namespace MathConst;
static constexpr int LISTDELTA = 10000;
static constexpr double LB_FACTOR = 1.5;
static constexpr int CMAPMAX = 6; // max # of CMAP terms stored by one atom
static constexpr int CMAPDIM = 24; // grid map dimension is 24 x 24
static constexpr double CMAPXMIN = -360.0;
static constexpr double CMAPXMIN2 = -180.0;
static constexpr double CMAPDX = 15.0; // 360/CMAPDIM
/* ---------------------------------------------------------------------- */
template<class DeviceType>
@ -246,6 +256,9 @@ void FixCMAPKokkos<DeviceType>::post_force(int vflag)
Kokkos::parallel_for(ncrosstermlist, *this);
copymode = 0;
atomKK->modified(execution_space,F_MASK);
std::cerr << fmt::format("*** post_force ncrosstermlist {} vflag {} ecmap {}\n",ncrosstermlist,vflag,ecmap);
}
@ -423,14 +436,19 @@ void FixCMAPKokkos<DeviceType>::operator()(const int n) const
// sum up cmap energy contributions
/* FIXME: needed for compute_scalar()
double ecmapKK = 0.0;
// FIXME: needed for compute_scalar()
double engfraction = 0.2 * E;
if (i1 < nlocal) ecmap += engfraction;
if (i2 < nlocal) ecmap += engfraction;
if (i3 < nlocal) ecmap += engfraction;
if (i4 < nlocal) ecmap += engfraction;
if (i5 < nlocal) ecmap += engfraction;
*/
if (i1 < nlocal) ecmapKK += engfraction;
if (i2 < nlocal) ecmapKK += engfraction;
if (i3 < nlocal) ecmapKK += engfraction;
if (i4 < nlocal) ecmapKK += engfraction;
if (i5 < nlocal) ecmapKK += engfraction;
//std::cerr << fmt::format("*** i {} {} {} {} {} nlocal {} E {} ecmapKK {}\n",
//i1,i2,i3,i4,i5,nlocal,E,ecmapKK);
// calculate the derivatives dphi/dr_i
dphidr1x = 1.0*r32/a1sq*a1x;
@ -705,7 +723,8 @@ void FixCMAPKokkos<DeviceType>::bc_interpol(double x1, double x2, int low1, int
double &E, double &dEdPhi, double &dEdPsi ) const
{
// FUSE bc_coeff() and bc_interpol() inline functions for kokkos version
// FUSE bc_coeff() and bc_interpol() inline functions for
// KOKKOS version to avoid passing cij[][] array back and forth
// calculate the bicubic interpolation coefficients c_ij

439
src/KOKKOS/superpose3d_kokkos.h Normal file
View File

@ -0,0 +1,439 @@
/* -*- c++ -*- ----------------------------------------------------------
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. (Some of the code in this file is also
available using a more premissive license. See below for details.)
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Andrew Jewett (Scripps Research)
Availability: https://github.com/jewettaij/superpose3d_cpp (MIT license)
------------------------------------------------------------------------- */
/// @file superpose3d.hpp
/// @brief Calculate the optimal rotation, translation and scale needed to
/// optimally fit two different point clouds containing n points.
/// @author Andrew Jewett
/// @license MIT
#ifndef LMP_SUPERPOSE3D_H
#define LMP_SUPERPOSE3D_H
#include "math_eigen_impl.h" //functions to calculate eigenvalues and eigenvectors
// -----------------------------------------------------------
// ------------------------ INTERFACE ------------------------
// -----------------------------------------------------------
/// @brief Superpose3d is a class with only one important member function
/// Superpose(). It is useful for calculating the optimal
/// superposition (rotations, translations, and scale transformations)
/// between two point clouds of the same size.
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray = Scalar const *>
class Superpose3D {
private:
size_t N; //number of points in the point clouds
Scalar *aWeights; //weights applied to points when computing RMSD
MathEigen::Jacobi<double, double *, double **> eigen_calc; // calc eigenvectors
Scalar **aaXf_shifted; //preallocated space for fixed point cloud (Nx3 array)
Scalar **aaXm_shifted; //preallocated space for mobile point cloud (Nx3 array)
public:
// The following data members store the rotation, translation and scale
// after optimal superposition
Scalar **R; //!< store optimal rotation here (this is a 3x3 array).
Scalar T[3]; //!< store optimal translation here
Scalar c; //!< store optimal scale (typically 1 unless requested by the user)
Scalar q[4]; //!< quaternion corresponding to the rotation stored in R.
// The first entry of q is cos(θ/2). The remaining 3 entries
// of q are the axis of rotation (with length sin(θ/2)).
// (Note: This is not the same as "p" from Diamond's 1988 paper.)
Superpose3D(size_t N = 0); //!< N=number of points in both point clouds
Superpose3D(size_t N, //!< N = number of points in both point clouds
ConstArray aWeights); //!< weight per point for computing RMSD
~Superpose3D();
/// @brief specify he number of points in both point clouds
void SetNumPoints(size_t N);
/// @brief return the number of points in both point clouds
size_t GetNumPoints() { return N; }
/// @brief specify the weight applied to each point when computing RMSD
void SetWeights(ConstArray aWeights);
/// @brief Use rigid-body transformations (rotations, translations, and
/// optionally scale transformations) to superimpose two point clouds.
///
/// @details
/// This function takes two lists of xyz coordinates (of the same length) and
/// attempts to superimpose them using rotations, translations, and
/// (optionally) scale transformations. These transformations are applied to
/// to the coordinates in the "aaXm_orig" array (the "mobile" point cloud)
/// in order to minimize the root-mean-squared-distance (RMSD) between the
/// corresponding points in each cloud, where RMSD is defined as:
///
/// @verbatim
/// sqrt((Σ_n w[n]*Σ_i |X[n][i] - (Σ_j c*R[i][j]*x[n][j]+T[i])|^2)/(Σ_n w[n]))
/// @endverbatim
///
/// In this formula, the "X_i" and "x_i" are coordinates of the ith fixed and
/// mobile point clouds (represented by "aaXf" and "aaXm" in the code below)
/// and "w_i" are optional weights (represented by "aWeights" in the code).
/// This function implements a more general variant of the method from:
/// @verbatim
/// R. Diamond, (1988) "A Note on the Rotational Superposition Problem",
/// Acta Cryst. A44, pp. 211-216
/// @endverbatim
///
/// @note:
/// This code has been augmented with a new feature. The version in the
/// original paper only considers rotation and translation and does not allow
/// coordinates of either cloud to be rescaled (multiplied by a scalar).
/// To enable the ability to rescale the coordinates, set allow_rescale=true.
/// (By default, this feature is disabled.)
///
/// @returns
/// The RMSD between the 2 pointclouds after optimal rotation, translation
/// (and scaling if requested) was applied to the "mobile" point cloud.
/// After this function is called, the optimal rotation, translation,
/// and scale (if requested) will be stored in the "R", "T", and "c"
/// public data members.
Scalar Superpose(ConstArrayOfCoords aaXf, //!< coords for the "frozen" object
ConstArrayOfCoords aaXm, //!< coords for the "mobile" object
bool allow_rescale = false //!< rescale mobile object? (c≠1?)
);
// C++ boilerplate: copy and move constructor, swap, and assignment operator
Superpose3D(const Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> &source);
Superpose3D(Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> &&other);
void swap(Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> &other);
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> &
operator=(Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> source);
private:
// memory management:
void Alloc(size_t N);
void Init();
void Dealloc();
}; // class Superpose3D
// -------------- IMPLEMENTATION --------------
template <typename Scalar> static inline Scalar SQR(Scalar x)
{
return x * x;
}
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
Scalar Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::Superpose(
ConstArrayOfCoords aaXf, // coords for the "frozen" object
ConstArrayOfCoords aaXm, // coords for the "mobile" object
bool allow_rescale) // rescale mobile object? (c!=1?)
{
// Find the center of mass of each object:
Scalar aCenter_f[3] = {0.0, 0.0, 0.0};
Scalar aCenter_m[3] = {0.0, 0.0, 0.0};
Scalar sum_weights = 0.0;
for (size_t n = 0; n < N; n++) {
Scalar weight = aWeights[n];
for (int d = 0; d < 3; d++) {
aCenter_f[d] += aaXf[n][d] * weight;
aCenter_m[d] += aaXm[n][d] * weight;
}
sum_weights += weight;
}
//assert(sum_weights != 0.0);
for (int d = 0; d < 3; d++) {
aCenter_f[d] /= sum_weights;
aCenter_m[d] /= sum_weights;
}
//Subtract the centers-of-mass from the original coordinates for each object
for (size_t n = 0; n < N; n++) {
for (int d = 0; d < 3; d++) {
// shift the coordinates so that the new center of mass is at the origin
aaXf_shifted[n][d] = aaXf[n][d] - aCenter_f[d];
aaXm_shifted[n][d] = aaXm[n][d] - aCenter_m[d];
}
}
// Calculate the "M" array from the Diamond paper (equation 16)
Scalar M[3][3];
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) M[i][j] = 0.0;
for (size_t n = 0; n < N; n++) {
Scalar weight = aWeights[n];
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) { M[i][j] += weight * aaXm_shifted[n][i] * aaXf_shifted[n][j]; }
}
}
// Calculate Q (equation 17)
Scalar traceM = 0.0;
for (int i = 0; i < 3; i++) traceM += M[i][i];
Scalar Q[3][3];
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
Q[i][j] = M[i][j] + M[j][i];
if (i == j) Q[i][j] -= 2.0 * traceM;
}
}
// Calculate V (equation 18)
Scalar V[3];
V[0] = M[1][2] - M[2][1];
V[1] = M[2][0] - M[0][2];
V[2] = M[0][1] - M[1][0];
// Calculate "P" (equation 22)
// First we must allocate space for the P matrix. It's not safe to declare:
// Scalar P[4][4];
// ...because most matrix solvers expect arrays in pointer-to-pointer format.
// (a different format). Below I create a fixed size matrix P in this format.
Scalar _PF[4 * 4]; // Contiguous 1D array for storing contents of the 2D P array
Scalar *P[4]; // This version of P has has ** (pointer-to-pointer) format.
for (int i = 0; i < 4; i++) // We must make sure that
P[i] = &(_PF[4 * i]); // P[i] points to the appropriate location in memory
// Now fill the P array
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) P[i][j] = Q[i][j];
P[0][3] = V[0];
P[3][0] = V[0];
P[1][3] = V[1];
P[3][1] = V[1];
P[2][3] = V[2];
P[3][2] = V[2];
P[3][3] = 0.0;
// The vector "p" contains the optimal rotation (backwards quaternion format)
Scalar p[4] = {0.0, 0.0, 0.0, 1.0}; // default value
Scalar pPp = 0.0; // = p^T * P * p (zero by default)
Scalar rmsd = 0.0; // default value
bool singular = N < 2; // (it doesn't make sense to rotate a single point)
if (!singular) {
// Calculate the principal eigenvalue and eigenvector of matrix P.
// Store the principal eigenvector in "p"
// The vector "p" will contain the optimal rotation (in quaternion format)
Scalar Evl[4]; // Store the eigenvalues of P here.
Scalar *Evc[4]; // Store the eigevectors here. This version has ** format.
Scalar _Evc[4 * 4]; // Contiguous 1D array for storing contents of "Evc" array
for (int i = 0; i < 4; i++) // We must make sure that
Evc[i] = &(_Evc[4 * i]); // Evc[i] points to the correct location in memory
eigen_calc.Diagonalize(P, Evl, Evc);
// Note: The eigenvalues are sorted in decreasing order by default.
pPp = Evl[0]; // = the maximum eigenvalue of P
for (int i = 0; i < 4; i++)
p[i] = Evc[0][i]; //copy eigenvector corresponding to this eigenvalue to p
} //if (! singular)
// Now normalize p
Scalar pnorm = 0.0;
for (int i = 0; i < 4; i++) pnorm += p[i] * p[i];
pnorm = sqrt(pnorm);
for (int i = 0; i < 4; i++) p[i] /= pnorm;
// Finally, calculate the rotation matrix corresponding to "p"
// (convert a quaternion into a 3x3 rotation matrix)
R[0][0] = (p[0] * p[0]) - (p[1] * p[1]) - (p[2] * p[2]) + (p[3] * p[3]);
R[1][1] = -(p[0] * p[0]) + (p[1] * p[1]) - (p[2] * p[2]) + (p[3] * p[3]);
R[2][2] = -(p[0] * p[0]) - (p[1] * p[1]) + (p[2] * p[2]) + (p[3] * p[3]);
R[0][1] = 2 * (p[0] * p[1] - p[2] * p[3]);
R[1][0] = 2 * (p[0] * p[1] + p[2] * p[3]);
R[1][2] = 2 * (p[1] * p[2] - p[0] * p[3]);
R[2][1] = 2 * (p[1] * p[2] + p[0] * p[3]);
R[0][2] = 2 * (p[0] * p[2] + p[1] * p[3]);
R[2][0] = 2 * (p[0] * p[2] - p[1] * p[3]);
q[0] = p[3]; // Note: The "p" variable is not a quaternion in the
q[1] = p[0]; // conventional sense because its elements
q[2] = p[1]; // are in the wrong order. I correct for that here.
q[3] = p[2]; // "q" is the quaternion correspond to rotation R.
// Optional: Decide the scale factor, c
c = 1.0; // by default, don't rescale the coordinates
if ((allow_rescale) && (!singular)) {
Scalar Waxaixai = 0.0;
Scalar WaxaiXai = 0.0;
for (size_t a = 0; a < N; a++) {
Scalar weight = aWeights[a];
for (int i = 0; i < 3; i++) {
Waxaixai += weight * aaXm_shifted[a][i] * aaXm_shifted[a][i];
WaxaiXai += weight * aaXm_shifted[a][i] * aaXf_shifted[a][i];
}
}
c = (WaxaiXai + pPp) / Waxaixai;
} // if (allow_rescale)
// Finally compute the RMSD between the two coordinate sets:
// First compute E0 from equation 24 of the paper
Scalar E0 = 0.0;
for (size_t n = 0; n < N; n++) {
Scalar weight = aWeights[n];
for (int d = 0; d < 3; d++)
// (remember to include the scale factor "c" that we inserted)
E0 += weight * (SQR(aaXf_shifted[n][d] - c * aaXm_shifted[n][d]));
}
Scalar sum_sqr_dist = E0 - c * 2.0 * pPp;
if (sum_sqr_dist < 0.0) //(edge case due to rounding error)
sum_sqr_dist = 0.0;
if (!singular) rmsd = sqrt(sum_sqr_dist / sum_weights);
// Lastly, calculate the translational offset.
// If c!=1, this is slightly more complicated than it seems. Recall that:
//RMSD=sqrt((Sum_i w_i * |X_i - Sum_j(c*R_ij*x_j + T_i))|^2) / (Sum_j w_j))
// =sqrt((Sum_i w_i * |X_i - x_i')|^2) / (Sum_j w_j))
// where
// x_i' = Sum_j(c*R_ij*x_j) + T_i
// = Xcm_i + c*R_ij*(x_j - xcm_j)
// and Xcm and xcm = center_of_mass for the frozen and mobile point clouds
//
// Hence:
// T_i = Xcm_i - Sum_j c*R_ij*xcm_j
// In the code, Xcm_i is represented by "aCenter_f[i]"
// and xcm_j is represented by "aCenter_m[j]"
for (int i = 0; i < 3; i++) {
T[i] = aCenter_f[i];
for (int j = 0; j < 3; j++) { T[i] -= c * R[i][j] * aCenter_m[j]; }
}
return rmsd;
} //Superpose3D::Superpose(aaXf, aaXm, allow_rescale)
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
void Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::SetNumPoints(size_t N)
{
Dealloc();
Alloc(N);
}
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
void Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::SetWeights(ConstArray aWeights)
{
for (size_t i = 0; i < N; i++) this->aWeights[i] = aWeights[i];
}
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::Superpose3D(size_t N) : eigen_calc(4)
{
Init();
Alloc(N);
}
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::Superpose3D(size_t N, ConstArray aWeights) :
eigen_calc(4)
{
Init();
Alloc(N);
SetWeights(aWeights);
}
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::~Superpose3D()
{
Dealloc();
}
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
void Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::Init()
{
R = nullptr;
aWeights = nullptr;
aaXf_shifted = nullptr;
aaXm_shifted = nullptr;
}
// memory management:
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
void Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::Alloc(size_t N)
{
this->N = N;
aWeights = new Scalar[N];
for (size_t i = 0; i < N; i++) aWeights[i] = 1.0 / N;
MathEigen::Alloc2D(3, 3, &R);
MathEigen::Alloc2D(N, 3, &aaXf_shifted);
MathEigen::Alloc2D(N, 3, &aaXm_shifted);
}
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
void Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::Dealloc()
{
if (R) MathEigen::Dealloc2D(&R);
if (aWeights) delete[] aWeights;
if (aaXf_shifted) MathEigen::Dealloc2D(&aaXf_shifted);
if (aaXm_shifted) MathEigen::Dealloc2D(&aaXm_shifted);
}
// memory management: copy and move constructor, swap, and assignment operator:
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::Superpose3D(
const Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> &source) :
eigen_calc(4)
{
Init();
Alloc(source.N);
//assert(N == source.N);
for (int i = 0; i < N; i++) {
std::copy(source.aaXf_shifted[i], source.aaXf_shifted[i] + 3, aaXf_shifted[i]);
std::copy(source.aaXm_shifted[i], source.aaXm_shifted[i] + 3, aaXm_shifted[i]);
}
}
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
void Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::swap(
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> &other)
{
std::swap(N, other.N);
std::swap(R, other.R);
std::swap(aaXf_shifted, other.aaXf_shifted);
std::swap(aaXm_shifted, other.aaXm_shifted);
}
// Move constructor (C++11)
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::Superpose3D(
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> &&other)
{
Init();
swap(*this, other);
}
// Using the "copy-swap" idiom for the assignment operator
template <typename Scalar, typename ConstArrayOfCoords, typename ConstArray>
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> &
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray>::operator=(
Superpose3D<Scalar, ConstArrayOfCoords, ConstArray> source)
{
this->swap(source);
return *this;
}
#endif //#ifndef LMP_SUPERPOSE3D_H

31
src/MOLECULE/fix_cmap.cpp
View File

@ -53,9 +53,14 @@ using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
//static constexpr int LISTDELTA = 10000;
//static constexpr double LB_FACTOR = 1.5;
static constexpr int LISTDELTA = 10000;
static constexpr double LB_FACTOR = 1.5;
static constexpr int CMAPMAX = 6; // max # of CMAP terms stored by one atom
static constexpr int CMAPDIM = 24; // grid map dimension is 24 x 24
static constexpr double CMAPXMIN = -360.0;
static constexpr double CMAPXMIN2 = -180.0;
static constexpr double CMAPDX = 15.0; // 360/CMAPDIM
/* ---------------------------------------------------------------------- */
@ -70,6 +75,8 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) :
{
if (narg != 4) error->all(FLERR,"Illegal fix cmap command");
std::cerr << "*** FixCMAP constructor\n";
restart_global = 1;
restart_peratom = 1;
energy_global_flag = energy_peratom_flag = 1;
@ -305,7 +312,7 @@ void FixCMAP::pre_reverse(int eflag, int /*vflag*/)
void FixCMAP::post_force(int vflag)
{
int n,i1,i2,i3,i4,i5,type,nlist;
int i1,i2,i3,i4,i5,type,nlist;
int li1, li2, mli1,mli2,mli11,mli21,t1,li3,li4,mli3,mli4,mli31,mli41;
int list[5];
// vectors needed to calculate the cross-term dihedral angles
@ -338,11 +345,11 @@ void FixCMAP::post_force(int vflag)
double **f = atom->f;
int nlocal = atom->nlocal;
ecmap = 0.0;
//if( ncrosstermlist>0 ) ecmap = 0.0;
int eflag = eflag_caller;
ev_init(eflag,vflag);
for (n = 0; n < ncrosstermlist; n++) {
for (int n = 0; n < ncrosstermlist; n++) {
i1 = crosstermlist[n][0];
i2 = crosstermlist[n][1];
i3 = crosstermlist[n][2];
@ -489,6 +496,9 @@ void FixCMAP::post_force(int vflag)
if (i4 < nlocal) ecmap += engfraction;
if (i5 < nlocal) ecmap += engfraction;
//std::cerr << fmt::format("*** i {} {} {} {} {} nlocal {} E {} ecmap {}\n",
//i1,i2,i3,i4,i5,nlocal,E,ecmap);
// calculate the derivatives dphi/dr_i
dphidr1x = 1.0*r32/a1sq*a1x;
@ -598,8 +608,13 @@ void FixCMAP::post_force(int vflag)
ev_tally(nlist,list,5.0,E,vcmap);
//ev_tally(5,list,nlocal,newton_bond,E,vcmap);
}
std::cerr << fmt::format("*** n {} ecmap {}\n",n,ecmap);
}
std::cerr << fmt::format("*** post_force eflag {} eflag_caller {} evflag {} thermo_energy {} ncrosstermlist {} vflag {} ecmap {}\n",eflag,eflag_caller,evflag,thermo_energy,ncrosstermlist,vflag,ecmap);
}
/* ---------------------------------------------------------------------- */
@ -607,6 +622,7 @@ void FixCMAP::post_force(int vflag)
void FixCMAP::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == ilevel_respa) post_force(vflag);
std::cerr << fmt::format("*** post_force_respa ecmap {}\n",ecmap);
}
/* ---------------------------------------------------------------------- */
@ -614,6 +630,7 @@ void FixCMAP::post_force_respa(int vflag, int ilevel, int /*iloop*/)
void FixCMAP::min_post_force(int vflag)
{
post_force(vflag);
std::cerr << fmt::format("*** min_post_force vflag {} ecmap {}\n",vflag,ecmap);
}
/* ----------------------------------------------------------------------
@ -623,8 +640,10 @@ void FixCMAP::min_post_force(int vflag)
double FixCMAP::compute_scalar()
{
double all;
MPI_Allreduce(&ecmap,&all,1,MPI_DOUBLE,MPI_SUM,world);
utils::logmesg(lmp, "compute_scalar = {}\n", all);
utils::logmesg(lmp, "compute_scalar: ecmap {} all {}\n", ecmap, all);
return all;
}

11
src/MOLECULE/fix_cmap.h
View File

@ -24,14 +24,6 @@ FixStyle(cmap,FixCMAP);
namespace LAMMPS_NS {
#define CMAPMAX 6 // max # of CMAP terms stored by one atom
#define CMAPDIM 24 // grid map dimension is 24 x 24
#define CMAPXMIN -360.0
#define CMAPXMIN2 -180.0
#define CMAPDX 15.0 // 360/CMAPDIM
#define LB_FACTOR 1.5
#define LISTDELTA 10000
class FixCMAP : public Fix {
public:
FixCMAP(class LAMMPS *, int, char **);
@ -73,6 +65,8 @@ class FixCMAP : public Fix {
double memory_usage() override;
double ecmap;
protected:
int eflag_caller;
int ctype, ilevel_respa;
@ -88,7 +82,6 @@ class FixCMAP : public Fix {
tagint **crossterm_atom1, **crossterm_atom2, **crossterm_atom3;
tagint **crossterm_atom4, **crossterm_atom5;
double ecmap;
double *g_axis;
// CMAP grid points obtained from external file

6
src/REACTION/fix_bond_react.h
View File

@ -11,10 +11,6 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing Author: Jacob Gissinger (jgissing@stevens.edu)
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(bond/react,FixBondReact);
@ -54,7 +50,7 @@ class FixBondReact : public Fix {
double compute_vector(int) override;
double memory_usage() override;
private:
protected:
int newton_bond;
int nreacts;
int *nevery;

184
unittest/force-styles/test_fix_timestep.cpp
View File

@ -335,10 +335,14 @@ TEST(FixTimestep, plain)
restart_lammps(lmp, test_config, false, false);
if (!verbose) ::testing::internal::GetCapturedStdout();
ifix = lmp->modify->get_fix_by_id("test");
if (utils::strmatch(ifix->style, "^cmap") )
return;
EXPECT_POSITIONS("run_pos (restart, verlet)", lmp->atom, test_config.run_pos, epsilon);
EXPECT_VELOCITIES("run_vel (restart, verlet)", lmp->atom, test_config.run_vel, epsilon);
ifix = lmp->modify->get_fix_by_id("test");
if (!ifix) {
FAIL() << "ERROR: no fix defined with fix ID 'test'\n";
} else {
@ -832,3 +836,181 @@ TEST(FixTimestep, omp)
cleanup_lammps(lmp, test_config);
if (!verbose) ::testing::internal::GetCapturedStdout();
};
TEST(FixTimestep, kokkos_omp)
{
if (!LAMMPS::is_installed_pkg("KOKKOS")) GTEST_SKIP();
if (test_config.skip_tests.count(test_info_->name())) GTEST_SKIP();
if (!Info::has_accelerator_feature("KOKKOS", "api", "openmp")) GTEST_SKIP();
LAMMPS::argv args = {"FixTimestep", "-log", "none", "-echo", "screen", "-nocite",
"-k", "on", "t", "4", "-sf", "kk"};
::testing::internal::CaptureStdout();
LAMMPS *lmp = init_lammps(args, test_config);
std::string output = ::testing::internal::GetCapturedStdout();
if (verbose) std::cout << output;
if (!lmp) {
std::cerr << "One or more prerequisite styles with /kk suffix\n"
"are not available in this LAMMPS configuration:\n";
for (auto &prerequisite : test_config.prerequisites) {
std::cerr << prerequisite.first << "_style " << prerequisite.second << "\n";
}
GTEST_SKIP();
}
EXPECT_THAT(output, StartsWith("LAMMPS ("));
EXPECT_THAT(output, HasSubstr("Loop time"));
// abort if running in parallel and not all atoms are local
const int nlocal = lmp->atom->nlocal;
ASSERT_EQ(lmp->atom->natoms, nlocal);
// relax error a bit for KOKKOS package
double epsilon = 5.0 * test_config.epsilon;
// relax test precision when using pppm and single precision FFTs
#if defined(FFT_SINGLE)
if (lmp->force->kspace && utils::strmatch(lmp->force->kspace_style, "^pppm")) epsilon *= 2.0e8;
#endif
ErrorStats stats;
EXPECT_POSITIONS("run_pos (normal run, verlet)", lmp->atom, test_config.run_pos, epsilon);
EXPECT_VELOCITIES("run_vel (normal run, verlet)", lmp->atom, test_config.run_vel, epsilon);
int ifix = lmp->modify->find_fix("test");
if (ifix < 0) {
FAIL() << "ERROR: no fix defined with fix ID 'test'\n";
} else {
Fix *fix = lmp->modify->fix[ifix];
if (fix->thermo_virial) {
EXPECT_STRESS("run_stress (normal run, verlet)", fix->virial, test_config.run_stress,
epsilon);
}
stats.reset();
// global scalar
if (fix->scalar_flag) {
double value = fix->compute_scalar();
EXPECT_FP_LE_WITH_EPS(test_config.global_scalar, value, epsilon);
}
// global vector
if (fix->vector_flag) {
int num = fix->size_vector;
EXPECT_EQ(num, test_config.global_vector.size());
for (int i = 0; i < num; ++i)
EXPECT_FP_LE_WITH_EPS(test_config.global_vector[i], fix->compute_vector(i),
epsilon);
}
// check t_target for thermostats
int dim = -1;
double *ptr = (double *)fix->extract("t_target", dim);
if ((ptr != nullptr) && (dim == 0)) {
int ivar = lmp->input->variable->find("t_target");
if (ivar >= 0) {
double t_ref = atof(lmp->input->variable->retrieve("t_target"));
double t_target = *ptr;
EXPECT_FP_LE_WITH_EPS(t_target, t_ref, epsilon);
}
}
if (print_stats && stats.has_data())
std::cerr << "global_data, normal run, verlet: " << stats << std::endl;
}
// FIXME: remove after debugging
if (utils::strmatch(lmp->modify->fix[ifix]->style, "^cmap") )
return;
if (!verbose) ::testing::internal::CaptureStdout();
restart_lammps(lmp, test_config, false, false);
if (!verbose) ::testing::internal::GetCapturedStdout();
EXPECT_POSITIONS("run_pos (restart, verlet)", lmp->atom, test_config.run_pos, epsilon);
EXPECT_VELOCITIES("run_vel (restart, verlet)", lmp->atom, test_config.run_vel, epsilon);
ifix = lmp->modify->find_fix("test");
if (ifix < 0) {
FAIL() << "ERROR: no fix defined with fix ID 'test'\n";
} else {
Fix *fix = lmp->modify->fix[ifix];
if (fix->thermo_virial) {
EXPECT_STRESS("run_stress (restart, verlet)", fix->virial, test_config.run_stress,
epsilon);
}
stats.reset();
// global scalar
if (fix->scalar_flag) {
double value = fix->compute_scalar();
EXPECT_FP_LE_WITH_EPS(test_config.global_scalar, value, epsilon);
}
// global vector
if (fix->vector_flag) {
int num = fix->size_vector;
EXPECT_EQ(num, test_config.global_vector.size());
for (int i = 0; i < num; ++i)
EXPECT_FP_LE_WITH_EPS(test_config.global_vector[i], fix->compute_vector(i),
epsilon);
}
if (print_stats && stats.has_data())
std::cerr << "global_data, restart, verlet: " << stats << std::endl;
}
if (lmp->atom->rmass == nullptr) {
if (!verbose) ::testing::internal::CaptureStdout();
restart_lammps(lmp, test_config, true, false);
if (!verbose) ::testing::internal::GetCapturedStdout();
EXPECT_POSITIONS("run_pos (rmass, verlet)", lmp->atom, test_config.run_pos, epsilon);
EXPECT_VELOCITIES("run_vel (rmass, verlet)", lmp->atom, test_config.run_vel, epsilon);
ifix = lmp->modify->find_fix("test");
if (ifix < 0) {
FAIL() << "ERROR: no fix defined with fix ID 'test'\n";
} else {
Fix *fix = lmp->modify->fix[ifix];
if (fix->thermo_virial) {
EXPECT_STRESS("run_stress (rmass, verlet)", fix->virial, test_config.run_stress,
epsilon);
}
stats.reset();
// global scalar
if (fix->scalar_flag) {
double value = fix->compute_scalar();
EXPECT_FP_LE_WITH_EPS(test_config.global_scalar, value, epsilon);
}
// global vector
if (fix->vector_flag) {
int num = fix->size_vector;
EXPECT_EQ(num, test_config.global_vector.size());
for (int i = 0; i < num; ++i)
EXPECT_FP_LE_WITH_EPS(test_config.global_vector[i], fix->compute_vector(i),
epsilon);
}
if (print_stats && stats.has_data())
std::cerr << "global_data, rmass, verlet: " << stats << std::endl;
}
}
// skip RESPA tests for KOKKOS
if (!verbose) ::testing::internal::CaptureStdout();
cleanup_lammps(lmp, test_config);
if (!verbose) ::testing::internal::GetCapturedStdout();
};