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

Merge branch 'master' into fix-ttm-grid

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-08-27 19:56:42 -04:00
parent e93ae9ba12 4d1a6cbd73
commit a66f3dfbdd
129 changed files with 4023 additions and 3006 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
src/INTERLAYER/pair_ilp_graphene_hbn.cpp
View File

@ -47,6 +47,7 @@ using namespace InterLayer;
#define PGDELTA 1
static const char cite_ilp[] =
"ilp/graphene/hbn potential doi:10.1021/acs.nanolett.8b02848\n"
"@Article{Ouyang2018\n"
" author = {W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod},\n"
" title = {Nanoserpents: Graphene Nanoribbon Motion on Two-Dimensional Hexagonal Materials},\n"
@ -482,7 +483,7 @@ void PairILPGrapheneHBN::calc_FRep(int eflag, int /* vflag */)
double dprodnorm1[3] = {0.0, 0.0, 0.0};
double fp1[3] = {0.0, 0.0, 0.0};
double fprod1[3] = {0.0, 0.0, 0.0};
double delkj[3] = {0.0, 0.0, 0.0};
double delki[3] = {0.0, 0.0, 0.0};
double fk[3] = {0.0, 0.0, 0.0};
inum = list->inum;
@ -588,12 +589,12 @@ void PairILPGrapheneHBN::calc_FRep(int eflag, int /* vflag */)
f[k][0] += fk[0];
f[k][1] += fk[1];
f[k][2] += fk[2];
delkj[0] = x[k][0] - x[j][0];
delkj[1] = x[k][1] - x[j][1];
delkj[2] = x[k][2] - x[j][2];
delki[0] = x[k][0] - x[i][0];
delki[1] = x[k][1] - x[i][1];
delki[2] = x[k][2] - x[i][2];
if (evflag)
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delkj[0],
delkj[1], delkj[2]);
ev_tally_xyz(k, i, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
delki[1], delki[2]);
}
if (eflag) pvector[1] += evdwl = Tap * Vilp;

13
src/INTERLAYER/pair_kolmogorov_crespi_full.cpp
View File

@ -47,6 +47,7 @@ using namespace InterLayer;
#define PGDELTA 1
static const char cite_kc[] =
"kolmogorov/crespi/full potential doi:10.1021/acs.nanolett.8b02848\n"
"@Article{Ouyang2018\n"
" author = {W. Ouyang, D. Mandelli, M. Urbakh, and O. Hod},\n"
" title = {Nanoserpents: Graphene Nanoribbon Motion on Two-Dimensional Hexagonal Materials},\n"
@ -478,7 +479,7 @@ void PairKolmogorovCrespiFull::calc_FRep(int eflag, int /* vflag */)
double dprodnorm1[3] = {0.0, 0.0, 0.0};
double fp1[3] = {0.0, 0.0, 0.0};
double fprod1[3] = {0.0, 0.0, 0.0};
double delkj[3] = {0.0, 0.0, 0.0};
double delki[3] = {0.0, 0.0, 0.0};
double fk[3] = {0.0, 0.0, 0.0};
inum = list->inum;
@ -585,12 +586,12 @@ void PairKolmogorovCrespiFull::calc_FRep(int eflag, int /* vflag */)
f[k][0] += fk[0];
f[k][1] += fk[1];
f[k][2] += fk[2];
delkj[0] = x[k][0] - x[j][0];
delkj[1] = x[k][1] - x[j][1];
delkj[2] = x[k][2] - x[j][2];
delki[0] = x[k][0] - x[i][0];
delki[1] = x[k][1] - x[i][1];
delki[2] = x[k][2] - x[i][2];
if (evflag)
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delkj[0],
delkj[1], delkj[2]);
ev_tally_xyz(k, j, nlocal, newton_pair, 0.0, 0.0, fk[0], fk[1], fk[2], delki[0],
delki[1], delki[2]);
}
if (eflag) {

26
src/INTERLAYER/pair_kolmogorov_crespi_z.cpp
View File

@ -28,6 +28,7 @@
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "potential_file_reader.h"
#include "tokenizer.h"
@ -158,7 +159,18 @@ void PairKolmogorovCrespiZ::compute(int eflag, int vflag)
if (eflag) { evdwl = -p.A * p.z06 / r6 + exp1 * sumCff - offset[itype][jtype]; }
if (evflag) {
ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0, fsum, fsum, fpair, delx, dely, delz);
ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
if (vflag_either) {
double fi[3],fj[3];
fi[0] = delx * fpair1;
fi[1] = dely * fpair1;
fi[2] = 0;
fj[0] = -delx * fpair1;
fj[1] = -dely * fpair1;
fj[2] = 0;
v_tally2_newton(i,fi,x[i]);
v_tally2_newton(j,fj,x[j]);
}
}
}
}
@ -228,6 +240,18 @@ void PairKolmogorovCrespiZ::coeff(int narg, char **arg)
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairKolmogorovCrespiZ::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style kolmogorov/crespi/z requires newton pair on");
neighbor->request(this,instance_me);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */

1
src/INTERLAYER/pair_kolmogorov_crespi_z.h
View File

@ -32,6 +32,7 @@ class PairKolmogorovCrespiZ : public Pair {
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
static constexpr int NPARAMS_PER_LINE = 11;

43
src/INTERLAYER/pair_lebedeva_z.cpp
View File

@ -30,6 +30,7 @@
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "potential_file_reader.h"
#include "tokenizer.h"
@ -137,16 +138,16 @@ void PairLebedevaZ::compute(int eflag, int vflag)
// derivatives
fpair = -6.0*p.A*p.z06/r8+p.B*p.alpha*exp1/r; // used for x,y,z
der = p.D1+2*p.D2*rhosq-p.lambda1*sumD; // used for x,y
fxy = fpair - 2*p.C*exp2*exp3*der;
fz = fpair + 2*p.C*p.lambda2*sumD*exp2*exp3;
fxy = 2*p.C*exp2*exp3*der;
fz = 2*p.C*p.lambda2*sumD*exp2*exp3;
f[i][0] += delx*fxy;
f[i][1] += dely*fxy;
f[i][2] += delz*fz;
f[i][0] += delx*(fpair-fxy);
f[i][1] += dely*(fpair-fxy);
f[i][2] += delz*(fpair+fz);
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fxy;
f[j][1] -= dely*fxy;
f[j][2] -= delz*fz;
f[j][0] -= delx*(fpair-fxy);
f[j][1] -= dely*(fpair-fxy);
f[j][2] -= delz*(fpair+fz);
}
if (eflag) {
@ -154,8 +155,18 @@ void PairLebedevaZ::compute(int eflag, int vflag)
}
if (evflag) {
ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0,
-fxy,-fxy,-fz,delx,dely,delz);
ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
if (vflag_either) {
double fi[3],fj[3];
fi[0] = -delx * fxy;
fi[1] = -dely * fxy;
fi[2] = delz * fz;
fj[0] = delx * fxy;
fj[1] = dely * fxy;
fj[2] = -delz * fz;
v_tally2_newton(i,fi,x[i]);
v_tally2_newton(j,fj,x[j]);
}
}
}
}
@ -226,6 +237,18 @@ void PairLebedevaZ::coeff(int narg, char **arg)
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairLebedevaZ::init_style()
{
if (force->newton_pair == 0)
error->all(FLERR,"Pair style lebedeva/z requires newton pair on");
neighbor->request(this,instance_me);
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */

1
src/INTERLAYER/pair_lebedeva_z.h
View File

@ -32,6 +32,7 @@ class PairLebedevaZ : public Pair {
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
static constexpr int NPARAMS_PER_LINE = 12;

18
src/KOKKOS/pair_snap_kokkos.h
View File

@ -44,7 +44,8 @@ struct TagPairSNAPComputeForce{};
struct TagPairSNAPComputeNeigh{};
struct TagPairSNAPComputeCayleyKlein{};
struct TagPairSNAPPreUi{};
struct TagPairSNAPComputeUi{};
struct TagPairSNAPComputeUiSmall{}; // more parallelism, more divergence
struct TagPairSNAPComputeUiLarge{}; // less parallelism, no divergence
struct TagPairSNAPTransformUi{}; // re-order ulisttot from SoA to AoSoA, zero ylist
struct TagPairSNAPComputeZi{};
struct TagPairSNAPBeta{};
@ -53,7 +54,9 @@ struct TagPairSNAPTransformBi{}; // re-order blist from AoSoA to AoS
struct TagPairSNAPComputeYi{};
struct TagPairSNAPComputeYiWithZlist{};
template<int dir>
struct TagPairSNAPComputeFusedDeidrj{};
struct TagPairSNAPComputeFusedDeidrjSmall{}; // more parallelism, more divergence
template<int dir>
struct TagPairSNAPComputeFusedDeidrjLarge{}; // less parallelism, no divergence
// CPU backend only
struct TagPairSNAPComputeNeighCPU{};
@ -143,7 +146,10 @@ public:
void operator() (TagPairSNAPPreUi,const int iatom_mod, const int j, const int iatom_div) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPComputeUi,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPComputeUi>::member_type& team) const;
void operator() (TagPairSNAPComputeUiSmall,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPComputeUiSmall>::member_type& team) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPComputeUiLarge,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPComputeUiLarge>::member_type& team) const;
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPTransformUi,const int iatom_mod, const int j, const int iatom_div) const;
@ -168,7 +174,11 @@ public:
template<int dir>
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPComputeFusedDeidrj<dir>,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPComputeFusedDeidrj<dir> >::member_type& team) const;
void operator() (TagPairSNAPComputeFusedDeidrjSmall<dir>,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPComputeFusedDeidrjSmall<dir> >::member_type& team) const;
template<int dir>
KOKKOS_INLINE_FUNCTION
void operator() (TagPairSNAPComputeFusedDeidrjLarge<dir>,const typename Kokkos::TeamPolicy<DeviceType, TagPairSNAPComputeFusedDeidrjLarge<dir> >::member_type& team) const;
// CPU backend only
KOKKOS_INLINE_FUNCTION

162
src/KOKKOS/pair_snap_kokkos_impl.h
View File

@ -341,18 +341,32 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::compute(int eflag_in,
// ComputeUi w/vector parallelism, shared memory, direct atomicAdd into ulisttot
{
// team_size_compute_ui is defined in `pair_snap_kokkos.h`
// scratch size: 32 atoms * (twojmax+1) cached values, no double buffer
const int tile_size = vector_length * (twojmax + 1);
const int scratch_size = scratch_size_helper<complex>(team_size_compute_ui * tile_size);
// total number of teams needed: (natoms / 32) * (max_neighs) * ("bend" locations)
const int n_teams = chunk_size_div * max_neighs * (twojmax + 1);
const int n_teams_div = (n_teams + team_size_compute_ui - 1) / team_size_compute_ui;
if (chunk_size < parallel_thresh)
{
// Version with parallelism over j_bend
SnapAoSoATeamPolicy<DeviceType, team_size_compute_ui, TagPairSNAPComputeUi> policy_ui(n_teams_div, team_size_compute_ui, vector_length);
policy_ui = policy_ui.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeUi",policy_ui,*this);
// total number of teams needed: (natoms / 32) * (max_neighs) * ("bend" locations)
const int n_teams = chunk_size_div * max_neighs * (twojmax + 1);
const int n_teams_div = (n_teams + team_size_compute_ui - 1) / team_size_compute_ui;
SnapAoSoATeamPolicy<DeviceType, team_size_compute_ui, TagPairSNAPComputeUiSmall> policy_ui(n_teams_div, team_size_compute_ui, vector_length);
policy_ui = policy_ui.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeUiSmall",policy_ui,*this);
} else {
// Version w/out parallelism over j_bend
// total number of teams needed: (natoms / 32) * (max_neighs)
const int n_teams = chunk_size_div * max_neighs;
const int n_teams_div = (n_teams + team_size_compute_ui - 1) / team_size_compute_ui;
SnapAoSoATeamPolicy<DeviceType, team_size_compute_ui, TagPairSNAPComputeUiLarge> policy_ui(n_teams_div, team_size_compute_ui, vector_length);
policy_ui = policy_ui.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeUiLarge",policy_ui,*this);
}
}
//TransformUi: un-"fold" ulisttot, zero ylist
@ -412,25 +426,51 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::compute(int eflag_in,
const int tile_size = vector_length * (twojmax + 1);
const int scratch_size = scratch_size_helper<complex>(2 * team_size_compute_fused_deidrj * tile_size);
// total number of teams needed: (natoms / 32) * (max_neighs) * ("bend" locations)
const int n_teams = chunk_size_div * max_neighs * (twojmax + 1);
const int n_teams_div = (n_teams + team_size_compute_fused_deidrj - 1) / team_size_compute_fused_deidrj;
if (chunk_size < parallel_thresh)
{
// Version with parallelism over j_bend
// x direction
SnapAoSoATeamPolicy<DeviceType, team_size_compute_fused_deidrj, TagPairSNAPComputeFusedDeidrj<0> > policy_fused_deidrj_x(n_teams_div,team_size_compute_fused_deidrj,vector_length);
policy_fused_deidrj_x = policy_fused_deidrj_x.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeFusedDeidrj<0>",policy_fused_deidrj_x,*this);
// total number of teams needed: (natoms / 32) * (max_neighs) * ("bend" locations)
const int n_teams = chunk_size_div * max_neighs * (twojmax + 1);
const int n_teams_div = (n_teams + team_size_compute_fused_deidrj - 1) / team_size_compute_fused_deidrj;
// y direction
SnapAoSoATeamPolicy<DeviceType, team_size_compute_fused_deidrj, TagPairSNAPComputeFusedDeidrj<1> > policy_fused_deidrj_y(n_teams_div,team_size_compute_fused_deidrj,vector_length);
policy_fused_deidrj_y = policy_fused_deidrj_y.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeFusedDeidrj<1>",policy_fused_deidrj_y,*this);
// x direction
SnapAoSoATeamPolicy<DeviceType, team_size_compute_fused_deidrj, TagPairSNAPComputeFusedDeidrjSmall<0> > policy_fused_deidrj_x(n_teams_div,team_size_compute_fused_deidrj,vector_length);
policy_fused_deidrj_x = policy_fused_deidrj_x.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeFusedDeidrjSmall<0>",policy_fused_deidrj_x,*this);
// z direction
SnapAoSoATeamPolicy<DeviceType, team_size_compute_fused_deidrj, TagPairSNAPComputeFusedDeidrj<2> > policy_fused_deidrj_z(n_teams_div,team_size_compute_fused_deidrj,vector_length);
policy_fused_deidrj_z = policy_fused_deidrj_z.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeFusedDeidrj<2>",policy_fused_deidrj_z,*this);
// y direction
SnapAoSoATeamPolicy<DeviceType, team_size_compute_fused_deidrj, TagPairSNAPComputeFusedDeidrjSmall<1> > policy_fused_deidrj_y(n_teams_div,team_size_compute_fused_deidrj,vector_length);
policy_fused_deidrj_y = policy_fused_deidrj_y.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeFusedDeidrjSmall<1>",policy_fused_deidrj_y,*this);
// z direction
SnapAoSoATeamPolicy<DeviceType, team_size_compute_fused_deidrj, TagPairSNAPComputeFusedDeidrjSmall<2> > policy_fused_deidrj_z(n_teams_div,team_size_compute_fused_deidrj,vector_length);
policy_fused_deidrj_z = policy_fused_deidrj_z.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeFusedDeidrjSmall<2>",policy_fused_deidrj_z,*this);
} else {
// Version w/out parallelism over j_bend
// total number of teams needed: (natoms / 32) * (max_neighs)
const int n_teams = chunk_size_div * max_neighs;
const int n_teams_div = (n_teams + team_size_compute_fused_deidrj - 1) / team_size_compute_fused_deidrj;
// x direction
SnapAoSoATeamPolicy<DeviceType, team_size_compute_fused_deidrj, TagPairSNAPComputeFusedDeidrjLarge<0> > policy_fused_deidrj_x(n_teams_div,team_size_compute_fused_deidrj,vector_length);
policy_fused_deidrj_x = policy_fused_deidrj_x.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeFusedDeidrjLarge<0>",policy_fused_deidrj_x,*this);
// y direction
SnapAoSoATeamPolicy<DeviceType, team_size_compute_fused_deidrj, TagPairSNAPComputeFusedDeidrjLarge<1> > policy_fused_deidrj_y(n_teams_div,team_size_compute_fused_deidrj,vector_length);
policy_fused_deidrj_y = policy_fused_deidrj_y.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeFusedDeidrjLarge<1>",policy_fused_deidrj_y,*this);
// z direction
SnapAoSoATeamPolicy<DeviceType, team_size_compute_fused_deidrj, TagPairSNAPComputeFusedDeidrjLarge<2> > policy_fused_deidrj_z(n_teams_div,team_size_compute_fused_deidrj,vector_length);
policy_fused_deidrj_z = policy_fused_deidrj_z.set_scratch_size(0, Kokkos::PerTeam(scratch_size));
Kokkos::parallel_for("ComputeFusedDeidrjLarge<2>",policy_fused_deidrj_z,*this);
}
}
#endif // LMP_KOKKOS_GPU
@ -603,13 +643,13 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSN
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
const auto idxb = icoeff % idxb_max;
const auto idx_chem = icoeff / idxb_max;
auto bveci = my_sna.blist(idxb, idx_chem, ii);
real_type bveci = my_sna.blist(ii, idx_chem, idxb);
d_beta_pack(iatom_mod,icoeff,iatom_div) += d_coeffi[k]*bveci;
k++;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
const auto jdxb = jcoeff % idxb_max;
const auto jdx_chem = jcoeff / idxb_max;
real_type bvecj = my_sna.blist(jdxb, jdx_chem, ii);
real_type bvecj = my_sna.blist(ii, jdx_chem, jdxb);
d_beta_pack(iatom_mod,icoeff,iatom_div) += d_coeffi[k]*bvecj;
d_beta_pack(iatom_mod,jcoeff,iatom_div) += d_coeffi[k]*bveci;
k++;
@ -736,7 +776,7 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSN
template<class DeviceType, typename real_type, int vector_length>
KOKKOS_INLINE_FUNCTION
void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSNAPComputeUi,const typename Kokkos::TeamPolicy<DeviceType,TagPairSNAPComputeUi>::member_type& team) const {
void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSNAPComputeUiSmall,const typename Kokkos::TeamPolicy<DeviceType,TagPairSNAPComputeUiSmall>::member_type& team) const {
SNAKokkos<DeviceType, real_type, vector_length> my_sna = snaKK;
// extract flattened atom_div / neighbor number / bend location
@ -756,11 +796,37 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSN
const int ninside = d_ninside(ii);
if (jj >= ninside) return;
my_sna.compute_ui(team,iatom_mod, jbend, jj, iatom_div);
my_sna.compute_ui_small(team, iatom_mod, jbend, jj, iatom_div);
});
}
template<class DeviceType, typename real_type, int vector_length>
KOKKOS_INLINE_FUNCTION
void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSNAPComputeUiLarge,const typename Kokkos::TeamPolicy<DeviceType,TagPairSNAPComputeUiLarge>::member_type& team) const {
SNAKokkos<DeviceType, real_type, vector_length> my_sna = snaKK;
// extract flattened atom_div / neighbor number / bend location
int flattened_idx = team.team_rank() + team.league_rank() * team_size_compute_ui;
// extract neighbor index, iatom_div
int iatom_div = flattened_idx / max_neighs; // removed "const" to work around GCC 7 bug
int jj = flattened_idx - iatom_div * max_neighs;
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team, vector_length),
[&] (const int iatom_mod) {
const int ii = iatom_mod + vector_length * iatom_div;
if (ii >= chunk_size) return;
const int ninside = d_ninside(ii);
if (jj >= ninside) return;
my_sna.compute_ui_large(team,iatom_mod, jj, iatom_div);
});
}
template<class DeviceType, typename real_type, int vector_length>
KOKKOS_INLINE_FUNCTION
void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSNAPTransformUi,const int iatom_mod, const int idxu, const int iatom_div) const {
@ -861,9 +927,9 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSN
for (int itriple = 0; itriple < ntriples; itriple++) {
const auto blocal = my_sna.blist_pack(iatom_mod, idxb, itriple, iatom_div);
const real_type blocal = my_sna.blist_pack(iatom_mod, idxb, itriple, iatom_div);
my_sna.blist(idxb, itriple, iatom) = blocal;
my_sna.blist(iatom, itriple, idxb) = blocal;
}
}
@ -871,7 +937,7 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSN
template<class DeviceType, typename real_type, int vector_length>
template<int dir>
KOKKOS_INLINE_FUNCTION
void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSNAPComputeFusedDeidrj<dir>,const typename Kokkos::TeamPolicy<DeviceType,TagPairSNAPComputeFusedDeidrj<dir> >::member_type& team) const {
void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSNAPComputeFusedDeidrjSmall<dir>,const typename Kokkos::TeamPolicy<DeviceType,TagPairSNAPComputeFusedDeidrjSmall<dir> >::member_type& team) const {
SNAKokkos<DeviceType, real_type, vector_length> my_sna = snaKK;
// extract flattened atom_div / neighbor number / bend location
@ -891,12 +957,38 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSN
const int ninside = d_ninside(ii);
if (jj >= ninside) return;
my_sna.template compute_fused_deidrj<dir>(team, iatom_mod, jbend, jj, iatom_div);
my_sna.template compute_fused_deidrj_small<dir>(team, iatom_mod, jbend, jj, iatom_div);
});
}
template<class DeviceType, typename real_type, int vector_length>
template<int dir>
KOKKOS_INLINE_FUNCTION
void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSNAPComputeFusedDeidrjLarge<dir>,const typename Kokkos::TeamPolicy<DeviceType,TagPairSNAPComputeFusedDeidrjLarge<dir> >::member_type& team) const {
SNAKokkos<DeviceType, real_type, vector_length> my_sna = snaKK;
// extract flattened atom_div / neighbor number / bend location
int flattened_idx = team.team_rank() + team.league_rank() * team_size_compute_fused_deidrj;
// extract neighbor index, iatom_div
int iatom_div = flattened_idx / max_neighs; // removed "const" to work around GCC 7 bug
int jj = flattened_idx - max_neighs * iatom_div;
Kokkos::parallel_for(Kokkos::ThreadVectorRange(team, vector_length),
[&] (const int iatom_mod) {
const int ii = iatom_mod + vector_length * iatom_div;
if (ii >= chunk_size) return;
const int ninside = d_ninside(ii);
if (jj >= ninside) return;
my_sna.template compute_fused_deidrj_large<dir>(team, iatom_mod, jj, iatom_div);
});
}
/* ----------------------------------------------------------------------
Begin routines that are unique to the CPU codepath. These do not take
advantage of AoSoA data layouts, but that could be a good point of
@ -925,13 +1017,13 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSN
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
const auto idxb = icoeff % idxb_max;
const auto idx_chem = icoeff / idxb_max;
auto bveci = my_sna.blist(idxb,idx_chem,ii);
real_type bveci = my_sna.blist(ii,idx_chem,idxb);
d_beta(icoeff,ii) += d_coeffi[k]*bveci;
k++;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
const auto jdxb = jcoeff % idxb_max;
const auto jdx_chem = jcoeff / idxb_max;
auto bvecj = my_sna.blist(jdxb,jdx_chem,ii);
real_type bvecj = my_sna.blist(ii,jdx_chem,jdxb);
d_beta(icoeff,ii) += d_coeffi[k]*bvecj;
d_beta(jcoeff,ii) += d_coeffi[k]*bveci;
k++;
@ -1221,7 +1313,7 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSN
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
const auto idxb = icoeff % idxb_max;
const auto idx_chem = icoeff / idxb_max;
evdwl += d_coeffi[icoeff+1]*my_sna.blist(idxb,idx_chem,ii);
evdwl += d_coeffi[icoeff+1]*my_sna.blist(ii,idx_chem,idxb);
}
// quadratic contributions
@ -1230,12 +1322,12 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::operator() (TagPairSN
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
const auto idxb = icoeff % idxb_max;
const auto idx_chem = icoeff / idxb_max;
auto bveci = my_sna.blist(idxb,idx_chem,ii);
real_type bveci = my_sna.blist(ii,idx_chem,idxb);
evdwl += 0.5*d_coeffi[k++]*bveci*bveci;
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++) {
auto jdxb = jcoeff % idxb_max;
auto jdx_chem = jcoeff / idxb_max;
auto bvecj = my_sna.blist(jdxb,jdx_chem,ii);
auto bvecj = my_sna.blist(ii,jdx_chem,jdxb);
evdwl += d_coeffi[k++]*bveci*bvecj;
}
}

65
src/KOKKOS/sna_kokkos.h
View File

@ -45,12 +45,12 @@ struct WignerWrapper {
{ ; }
KOKKOS_INLINE_FUNCTION
complex get(const int& ma) {
complex get(const int& ma) const {
return complex(buffer[offset + 2 * vector_length * ma], buffer[offset + vector_length + 2 * vector_length * ma]);
}
KOKKOS_INLINE_FUNCTION
void set(const int& ma, const complex& store) {
void set(const int& ma, const complex& store) const {
buffer[offset + 2 * vector_length * ma] = store.re;
buffer[offset + vector_length + 2 * vector_length * ma] = store.im;
}
@ -122,8 +122,14 @@ inline
void compute_cayley_klein(const int&, const int&, const int&);
KOKKOS_INLINE_FUNCTION
void pre_ui(const int&, const int&, const int&, const int&); // ForceSNAP
// version of the code with parallelism over j_bend
KOKKOS_INLINE_FUNCTION
void compute_ui(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, const int, const int, const int, const int); // ForceSNAP
void compute_ui_small(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, const int, const int, const int, const int); // ForceSNAP
// version of the code without parallelism over j_bend
KOKKOS_INLINE_FUNCTION
void compute_ui_large(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, const int, const int, const int); // ForceSNAP
KOKKOS_INLINE_FUNCTION
void compute_zi(const int&, const int&, const int&); // ForceSNAP
KOKKOS_INLINE_FUNCTION
@ -135,6 +141,35 @@ inline
KOKKOS_INLINE_FUNCTION
void compute_bi(const int&, const int&, const int&); // ForceSNAP
// functions for derivatives, GPU only
// version of the code with parallelism over j_bend
template<int dir>
KOKKOS_INLINE_FUNCTION
void compute_fused_deidrj_small(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, const int, const int, const int, const int); //ForceSNAP
// version of the code without parallelism over j_bend
template<int dir>
KOKKOS_INLINE_FUNCTION
void compute_fused_deidrj_large(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, const int, const int, const int); //ForceSNAP
// core "evaluation" functions that get plugged into "compute" functions
// plugged into compute_ui_small, compute_ui_large
KOKKOS_FORCEINLINE_FUNCTION
void evaluate_ui_jbend(const WignerWrapper<real_type, vector_length>&, const complex&, const complex&, const real_type&, const int&,
const int&, const int&, const int&);
// plugged into compute_zi, compute_yi
KOKKOS_FORCEINLINE_FUNCTION
complex evaluate_zi(const int&, const int&, const int&, const int&, const int&, const int&, const int&, const int&, const int&,
const int&, const int&, const int&, const int&, const real_type*);
// plugged into compute_yi, compute_yi_with_zlist
KOKKOS_FORCEINLINE_FUNCTION
real_type evaluate_beta_scaled(const int&, const int&, const int&, const int&, const int&, const int&, const int&, const int&,
const Kokkos::View<real_type***, Kokkos::LayoutLeft, DeviceType> &);
// plugged into compute_fused_deidrj_small, compute_fused_deidrj_large
KOKKOS_FORCEINLINE_FUNCTION
real_type evaluate_duidrj_jbend(const WignerWrapper<real_type, vector_length>&, const complex&, const complex&, const real_type&,
const WignerWrapper<real_type, vector_length>&, const complex&, const complex&, const real_type&,
const int&, const int&, const int&, const int&);
// functions for bispectrum coefficients, CPU only
KOKKOS_INLINE_FUNCTION
void pre_ui_cpu(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team,const int&,const int&); // ForceSNAP
@ -148,11 +183,6 @@ inline
KOKKOS_INLINE_FUNCTION
void compute_bi_cpu(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, int); // ForceSNAP
// functions for derivatives, GPU only
template<int dir>
KOKKOS_INLINE_FUNCTION
void compute_fused_deidrj(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, const int, const int, const int, const int); //ForceSNAP
// functions for derivatives, CPU only
KOKKOS_INLINE_FUNCTION
void compute_duidrj_cpu(const typename Kokkos::TeamPolicy<DeviceType>::member_type& team, int, int); //ForceSNAP
@ -168,23 +198,6 @@ inline
KOKKOS_INLINE_FUNCTION
void compute_s_dsfac(const real_type, const real_type, real_type&, real_type&); // compute_cayley_klein
static KOKKOS_FORCEINLINE_FUNCTION
void sincos_wrapper(double x, double* sin_, double *cos_) {
#ifdef __SYCL_DEVICE_ONLY__
*sin_ = sycl::sincos(x, cos_);
#else
sincos(x, sin_, cos_);
#endif
}
static KOKKOS_FORCEINLINE_FUNCTION
void sincos_wrapper(float x, float* sin_, float *cos_) {
#ifdef __SYCL_DEVICE_ONLY__
*sin_ = sycl::sincos(x, cos_);
#else
sincosf(x, sin_, cos_);
#endif
}
#ifdef TIMING_INFO
double* timers;
timespec starttime, endtime;
@ -207,7 +220,7 @@ inline
int twojmax, diagonalstyle;
t_sna_3d_ll blist;
t_sna_3d blist;
t_sna_3c_ll ulisttot;
t_sna_3c_ll ulisttot_full; // un-folded ulisttot, cpu only
t_sna_3c_ll zlist;

855
src/KOKKOS/sna_kokkos_impl.h
View File

File diff suppressed because it is too large Load Diff

1
src/MC/fix_atom_swap.cpp
View File

@ -516,7 +516,6 @@ double FixAtomSwap::energy_full()
if (force->kspace) force->kspace->compute(eflag,vflag);
if (modify->n_post_force) modify->post_force(vflag);
if (modify->n_end_of_step) modify->end_of_step();
update->eflag_global = update->ntimestep;
double total_energy = c_pe->compute_scalar();

2
src/MC/fix_charge_regulation.cpp
View File

@ -1108,7 +1108,7 @@ double FixChargeRegulation::energy_full() {
if (modify->n_pre_reverse) modify->pre_reverse(eflag,vflag);
if (modify->n_post_force) modify->post_force(vflag);
if (modify->n_end_of_step) modify->end_of_step();
update->eflag_global = update->ntimestep;
double total_energy = c_pe->compute_scalar();
return total_energy;

3
src/MC/fix_gcmc.cpp
View File

@ -2319,10 +2319,9 @@ double FixGCMC::energy_full()
if (modify->n_pre_reverse) modify->pre_reverse(eflag,vflag);
if (modify->n_post_force) modify->post_force(vflag);
if (modify->n_end_of_step) modify->end_of_step();
// NOTE: all fixes with energy_global_flag set and which
// operate at pre_force() or post_force() or end_of_step()
// operate at pre_force() or post_force()
// and which user has enabled via fix_modify energy yes,
// will contribute to total MC energy via pe->compute_scalar()

3
src/MC/fix_widom.cpp
View File

@ -1053,10 +1053,9 @@ double FixWidom::energy_full()
if (modify->n_pre_reverse) modify->pre_reverse(eflag,vflag);
if (modify->n_pre_force) modify->pre_force(vflag);
if (modify->n_end_of_step) modify->end_of_step();
// NOTE: all fixes with energy_global_flag set and which
// operate at pre_force() or post_force() or end_of_step()
// operate at pre_force() or post_force()
// and which user has enabled via fix_modify energy yes,
// will contribute to total MC energy via pe->compute_scalar()

5
src/ML-SNAP/pair_snap.cpp
View File

@ -628,7 +628,8 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;
chunksize = 4096;
chunksize = 32768;
parallel_thresh = 8192;
// open SNAP parameter file on proc 0
@ -696,6 +697,8 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
wselfallflag = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
else if (keywd == "chunksize")
chunksize = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
else if (keywd == "parallelthresh")
parallel_thresh = utils::inumeric(FLERR,keyval.c_str(),false,lmp);
else
error->all(FLERR,"Unknown parameter '{}' in SNAP "
"parameter file", keywd);

2
src/ML-SNAP/pair_snap.h
View File

@ -59,7 +59,7 @@ class PairSNAP : public Pair {
double **scale; // for thermodynamic integration
int twojmax, switchflag, bzeroflag, bnormflag;
int chemflag, wselfallflag;
int chunksize;
int chunksize,parallel_thresh;
double rfac0, rmin0, wj1, wj2;
int rcutfacflag, twojmaxflag; // flags for required parameters
int beta_max; // length of beta

2
src/MOFFF/README
View File

@ -20,7 +20,7 @@ charges (dsf and long-range treatment of charges)
out-of-plane angle
See the file doc/drude_tutorial.html for getting started.
See the doc pages for "pair_style buck6d/coul/gauss", "anlge_style class2",
See the doc pages for "pair_style buck6d/coul/gauss", "angle_style class2",
"angle_style cosine/buck6d", and "improper_style inversion/harmonic"
commands to get started. Also see the above mentioned website and
literature for further documentation about the force field.

64
src/OPENMP/domain_omp.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -16,13 +15,16 @@
Contributing author : Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "omp_compat.h"
#include "accelerator_omp.h"
#include "atom.h"
#include "error.h"
#include "omp_compat.h"
using namespace LAMMPS_NS;
typedef struct { double x,y,z; } dbl3_t;
typedef struct {
double x, y, z;
} dbl3_t;
/* ----------------------------------------------------------------------
enforce PBC and modify box image flags for each atom
@ -37,20 +39,35 @@ typedef struct { double x,y,z; } dbl3_t;
void DomainOMP::pbc()
{
dbl3_t * _noalias const x = (dbl3_t *)&atom->x[0][0];
dbl3_t * _noalias const v = (dbl3_t *)&atom->v[0][0];
const double * _noalias const lo = (triclinic == 0) ? boxlo : boxlo_lamda;
const double * _noalias const hi = (triclinic == 0) ? boxhi : boxhi_lamda;
const double * _noalias const period = (triclinic == 0) ? prd : prd_lamda;
const int * _noalias const mask = atom->mask;
imageint * _noalias const image = atom->image;
const int nlocal = atom->nlocal;
if (!nlocal) return;
// verify owned atoms have valid numerical coords
// may not if computed pairwise force between 2 atoms at same location
const double *_noalias const coord = &atom->x[0][0];
const int n3 = 3 * nlocal;
int flag = 0;
#if defined(_OPENMP) // clang-format off
#pragma omp parallel for LMP_DEFAULT_NONE schedule(static) reduction(+:flag)
#endif // clang-format on
for (int i = 0; i < n3; i++)
if (!std::isfinite(coord[i])) flag = 1;
if (flag) error->one(FLERR, "Non-numeric atom coords - simulation unstable");
dbl3_t *_noalias const x = (dbl3_t *) &atom->x[0][0];
dbl3_t *_noalias const v = (dbl3_t *) &atom->v[0][0];
const double *_noalias const lo = (triclinic == 0) ? boxlo : boxlo_lamda;
const double *_noalias const hi = (triclinic == 0) ? boxhi : boxhi_lamda;
const double *_noalias const period = (triclinic == 0) ? prd : prd_lamda;
const int *_noalias const mask = atom->mask;
imageint *_noalias const image = atom->image;
#if defined(_OPENMP)
#pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
#endif
for (int i = 0; i < nlocal; i++) {
imageint idim,otherdims;
imageint idim, otherdims;
if (xperiodic) {
if (x[i].x < lo[0]) {
@ -64,7 +81,7 @@ void DomainOMP::pbc()
}
if (x[i].x >= hi[0]) {
x[i].x -= period[0];
x[i].x = MAX(x[i].x,lo[0]);
x[i].x = MAX(x[i].x, lo[0]);
if (deform_vremap && mask[i] & deform_groupbit) v[i].x -= h_rate[0];
idim = image[i] & IMGMASK;
otherdims = image[i] ^ idim;
@ -89,7 +106,7 @@ void DomainOMP::pbc()
}
if (x[i].y >= hi[1]) {
x[i].y -= period[1];
x[i].y = MAX(x[i].y,lo[1]);
x[i].y = MAX(x[i].y, lo[1]);
if (deform_vremap && mask[i] & deform_groupbit) {
v[i].x -= h_rate[5];
v[i].y -= h_rate[1];
@ -118,7 +135,7 @@ void DomainOMP::pbc()
}
if (x[i].z >= hi[2]) {
x[i].z -= period[2];
x[i].z = MAX(x[i].z,lo[2]);
x[i].z = MAX(x[i].z, lo[2]);
if (deform_vremap && mask[i] & deform_groupbit) {
v[i].x -= h_rate[4];
v[i].y -= h_rate[3];
@ -141,16 +158,17 @@ void DomainOMP::pbc()
void DomainOMP::lamda2x(int n)
{
dbl3_t * _noalias const x = (dbl3_t *)&atom->x[0][0];
const int num = n;
if (!n) return;
dbl3_t *_noalias const x = (dbl3_t *) &atom->x[0][0];
#if defined(_OPENMP)
#pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
#endif
for (int i = 0; i < num; i++) {
x[i].x = h[0]*x[i].x + h[5]*x[i].y + h[4]*x[i].z + boxlo[0];
x[i].y = h[1]*x[i].y + h[3]*x[i].z + boxlo[1];
x[i].z = h[2]*x[i].z + boxlo[2];
x[i].x = h[0] * x[i].x + h[5] * x[i].y + h[4] * x[i].z + boxlo[0];
x[i].y = h[1] * x[i].y + h[3] * x[i].z + boxlo[1];
x[i].z = h[2] * x[i].z + boxlo[2];
}
}
@ -161,8 +179,9 @@ void DomainOMP::lamda2x(int n)
void DomainOMP::x2lamda(int n)
{
dbl3_t * _noalias const x = (dbl3_t *)&atom->x[0][0];
const int num = n;
if (!n) return;
dbl3_t *_noalias const x = (dbl3_t *) &atom->x[0][0];
#if defined(_OPENMP)
#pragma omp parallel for LMP_DEFAULT_NONE schedule(static)
@ -172,9 +191,8 @@ void DomainOMP::x2lamda(int n)
double delta1 = x[i].y - boxlo[1];
double delta2 = x[i].z - boxlo[2];
x[i].x = h_inv[0]*delta0 + h_inv[5]*delta1 + h_inv[4]*delta2;
x[i].y = h_inv[1]*delta1 + h_inv[3]*delta2;
x[i].z = h_inv[2]*delta2;
x[i].x = h_inv[0] * delta0 + h_inv[5] * delta1 + h_inv[4] * delta2;
x[i].y = h_inv[1] * delta1 + h_inv[3] * delta2;
x[i].z = h_inv[2] * delta2;
}
}

187
src/REACTION/fix_bond_react.cpp
View File

@ -140,6 +140,7 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
rxnfunclist[0] = "rxnsum";
rxnfunclist[1] = "rxnave";
nvvec = 0;
ncustomvars = 0;
vvec = nullptr;
nxspecial = nullptr;
@ -232,7 +233,6 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
memory->create(molecule_keyword,nreacts,"bond/react:molecule_keyword");
memory->create(nconstraints,nreacts,"bond/react:nconstraints");
memory->create(constraintstr,nreacts,MAXLINE,"bond/react:constraintstr");
memory->create(constraints,0,nreacts,"bond/react:constraints");
memory->create(var_flag,NUMVARVALS,nreacts,"bond/react:var_flag");
memory->create(var_id,NUMVARVALS,nreacts,"bond/react:var_id");
memory->create(iatomtype,nreacts,"bond/react:iatomtype");
@ -488,6 +488,9 @@ FixBondReact::FixBondReact(LAMMPS *lmp, int narg, char **arg) :
if (custom_charges_fragid[i] >= 0) CustomCharges(custom_charges_fragid[i],i);
}
// get the names of per-atom variables needed by 'rxn' functions of custom constraint
customvarnames();
// initialize Marsaglia RNG with processor-unique seed (Arrhenius prob)
rrhandom = new RanMars*[narrhenius];
@ -617,7 +620,6 @@ FixBondReact::~FixBondReact()
memory->destroy(stabilize_steps_flag);
memory->destroy(custom_charges_fragid);
memory->destroy(molecule_keyword);
memory->destroy(constraints);
memory->destroy(nconstraints);
memory->destroy(constraintstr);
memory->destroy(create_atoms_flag);
@ -1020,6 +1022,11 @@ void FixBondReact::post_integrate()
return;
}
// evaluate custom constraint variable values here and forward_comm
get_customvars();
commflag = 1;
comm->forward_comm_fix(this,ncustomvars);
// run through the superimpose algorithm
// this checks if simulation topology matches unreacted mol template
superimpose_algorithm();
@ -2086,6 +2093,100 @@ double FixBondReact::get_temperature(tagint **myglove, int row_offset, int col)
return t;
}
/* ----------------------------------------------------------------------
get per-atom variable names used by custom constraint
------------------------------------------------------------------------- */
void FixBondReact::customvarnames()
{
int pos,pos1,pos2,pos3,prev3;
std::string varstr,argstr,varid;
// search all constraints' varstr for special 'rxn' functions
// add variable names to customvarstrs
// add values to customvars
for (rxnID = 0; rxnID < nreacts; rxnID++) {
for (int i = 0; i < nconstraints[rxnID]; i++) {
if (constraints[i][rxnID].type == CUSTOM) {
varstr = constraints[i][rxnID].str;
prev3 = -1;
while (true) {
// find next reaction special function occurrence
pos1 = INT_MAX;
for (int i = 0; i < nrxnfunction; i++) {
pos = varstr.find(rxnfunclist[i],prev3+1);
if (pos == std::string::npos) continue;
if (pos < pos1) pos1 = pos;
}
if (pos1 == INT_MAX) break;
pos2 = varstr.find("(",pos1);
pos3 = varstr.find(")",pos2);
if (pos2 == std::string::npos || pos3 == std::string::npos)
error->all(FLERR,"Bond/react: Illegal rxn function syntax\n");
prev3 = pos3;
argstr = varstr.substr(pos2+1,pos3-pos2-1);
argstr.erase(remove_if(argstr.begin(), argstr.end(), isspace), argstr.end()); // remove whitespace
pos2 = argstr.find(",");
if (pos2 != std::string::npos) varid = argstr.substr(0,pos2);
else varid = argstr;
// check if we already know about this variable
int varidflag = 0;
for (int j = 0; j < ncustomvars; j++) {
if (customvarstrs[j] == varid) {
varidflag = 1;
break;
}
}
if (!varidflag) {
customvarstrs.resize(ncustomvars+1);
customvarstrs[ncustomvars++] = varid;
}
}
}
}
}
}
/* ----------------------------------------------------------------------
evaluate per-atom variables needed for custom constraint
------------------------------------------------------------------------- */
void FixBondReact::get_customvars()
{
double *tempvvec;
std::string varid;
int nall = atom->nlocal + atom->nghost;
memory->create(tempvvec,nall,"bond/react:tempvvec");
if (vvec == nullptr) {
memory->create(vvec,nall,ncustomvars,"bond/react:vvec");
nvvec = nall;
}
if (nvvec < nall) {
memory->grow(vvec,nall,ncustomvars,"bond/react:vvec");
nvvec = nall;
}
for (int i = 0; i < ncustomvars; i++) {
varid = customvarstrs[i];
if (varid.substr(0,2) != "v_") error->all(FLERR,"Bond/react: Reaction special function variable "
"name should begin with 'v_'");
varid = varid.substr(2);
int ivar = input->variable->find(varid.c_str());
if (ivar < 0)
error->all(FLERR,"Bond/react: Reaction special function variable "
"name does not exist");
if (!input->variable->atomstyle(ivar))
error->all(FLERR,"Bond/react: Reaction special function must "
"reference an atom-style variable");
input->variable->compute_atom(ivar,igroup,tempvvec,1,0);
for (int j = 0; j < nall; j++) vvec[j][i] = tempvvec[j];
}
memory->destroy(tempvvec);
}
/* ----------------------------------------------------------------------
evaulate expression for variable constraint
------------------------------------------------------------------------- */
@ -2120,6 +2221,7 @@ double FixBondReact::custom_constraint(std::string varstr)
evlstr.push_back(varstr.substr(prev3+1,pos1-(prev3+1)));
prev3 = pos3;
argstr = varstr.substr(pos2+1,pos3-pos2-1);
argstr.erase(remove_if(argstr.begin(), argstr.end(), isspace), argstr.end()); // remove whitespace
pos2 = argstr.find(",");
if (pos2 != std::string::npos) {
varid = argstr.substr(0,pos2);
@ -2145,25 +2247,19 @@ currently two 'rxn' functions: rxnsum and rxnave
double FixBondReact::rxnfunction(std::string rxnfunc, std::string varid,
std::string fragid)
{
if (varid.substr(0,2) != "v_") error->one(FLERR,"Bond/react: Reaction special function variable "
"name should begin with 'v_'");
varid = varid.substr(2);
int ivar = input->variable->find(varid.c_str());
int ivar = -1;
for (int i = 0; i < ncustomvars; i++) {
if (varid == customvarstrs[i]) {
ivar = i;
break;
}
}
// variable name should always be found, at this point
// however, let's double check for completeness
if (ivar < 0)
error->one(FLERR,"Bond/react: Reaction special function variable "
"name does not exist");
if (!input->variable->atomstyle(ivar))
error->one(FLERR,"Bond/react: Reaction special function must "
"reference an atom-style variable");
if (vvec == nullptr) {
memory->create(vvec,atom->nlocal,"bond/react:vvec");
nvvec = atom->nlocal;
}
if (nvvec < atom->nlocal) {
memory->grow(vvec,atom->nlocal,"bond/react:vvec");
nvvec = atom->nlocal;
}
input->variable->compute_atom(ivar,igroup,vvec,1,0);
int ifrag = -1;
if (fragid != "all") {
ifrag = onemol->findfragment(fragid.c_str());
@ -2178,14 +2274,14 @@ double FixBondReact::rxnfunction(std::string rxnfunc, std::string varid,
if (fragid == "all") {
for (int i = 0; i < onemol->natoms; i++) {
iatom = atom->map(glove[i][1]);
sumvvec += vvec[iatom];
sumvvec += vvec[iatom][ivar];
}
nsum = onemol->natoms;
} else {
for (int i = 0; i < onemol->natoms; i++) {
if (onemol->fragmentmask[ifrag][i]) {
iatom = atom->map(glove[i][1]);
sumvvec += vvec[iatom];
sumvvec += vvec[iatom][ivar];
nsum++;
}
}
@ -2830,6 +2926,19 @@ void FixBondReact::update_everything()
rxnID = local_mega_glove[0][i];
// reactions already shuffled from dedup procedure, so can skip first N
if (iskip[rxnID]++ < nlocalskips[rxnID]) continue;
// atoms inserted here for serial MPI_STUBS build only
if (create_atoms_flag[rxnID] == 1) {
onemol = atom->molecules[unreacted_mol[rxnID]];
twomol = atom->molecules[reacted_mol[rxnID]];
if (insert_atoms(local_mega_glove,i)) {
inserted_atoms_flag = 1;
} else { // create aborted
reaction_count_total[rxnID]--;
continue;
}
}
for (int j = 0; j < max_natoms+1; j++)
update_mega_glove[j][update_num_mega] = local_mega_glove[j][i];
update_num_mega++;
@ -2842,7 +2951,7 @@ void FixBondReact::update_everything()
// we can insert atoms here, now that reactions are finalized
// can't do it any earlier, due to skipped reactions (max_rxn)
// reactions that create atoms are always treated as 'global'
// for MPI build, reactions that create atoms are always treated as 'global'
if (create_atoms_flag[rxnID] == 1) {
onemol = atom->molecules[unreacted_mol[rxnID]];
twomol = atom->molecules[reacted_mol[rxnID]];
@ -2858,17 +2967,18 @@ void FixBondReact::update_everything()
update_mega_glove[j][update_num_mega] = global_mega_glove[j][i];
update_num_mega++;
}
// if inserted atoms and global map exists, reset map now instead
// of waiting for comm since other pre-exchange fixes may use it
// invoke map_init() b/c atom count has grown
// do this once after all atom insertions
if (inserted_atoms_flag == 1 && atom->map_style != Atom::MAP_NONE) {
atom->map_init();
atom->map_set();
}
}
delete [] iskip;
// if inserted atoms and global map exists, reset map now instead
// of waiting for comm since other pre-exchange fixes may use it
// invoke map_init() b/c atom count has grown
// do this once after all atom insertions
if (inserted_atoms_flag == 1 && atom->map_style != Atom::MAP_NONE) {
atom->map_init();
atom->map_set();
}
// mark to-delete atoms
nlocal = atom->nlocal;
if (nlocal > nmark) {
@ -3668,7 +3778,7 @@ void FixBondReact::read(int myrxn)
else if (strstr(line,"constraints")) {
sscanf(line,"%d",&nconstraints[myrxn]);
if (maxnconstraints < nconstraints[myrxn]) maxnconstraints = nconstraints[myrxn];
memory->grow(constraints,maxnconstraints,nreacts,"bond/react:constraints");
constraints.resize(maxnconstraints, std::vector<Constraint>(nreacts));
} else break;
}
@ -4027,6 +4137,15 @@ int FixBondReact::pack_forward_comm(int n, int *list, double *buf,
m = 0;
if (commflag == 1) {
for (i = 0; i < n; i++) {
j = list[i];
for (k = 0; k < ncustomvars; k++)
buf[m++] = vvec[j][k];
}
return m;
}
if (commflag == 2) {
for (i = 0; i < n; i++) {
j = list[i];
@ -4051,12 +4170,16 @@ int FixBondReact::pack_forward_comm(int n, int *list, double *buf,
void FixBondReact::unpack_forward_comm(int n, int first, double *buf)
{
int i,j,m,ns,last;
int i,j,k,m,ns,last;
m = 0;
last = first + n;
if (commflag == 2) {
if (commflag == 1) {
for (i = first; i < last; i++)
for (k = 0; k < ncustomvars; k++)
vvec[i][k] = buf[m++];
} else if (commflag == 2) {
for (i = first; i < last; i++)
partner[i] = (tagint) ubuf(buf[m++]).i;
} else {

20
src/REACTION/fix_bond_react.h
View File

@ -139,8 +139,6 @@ class FixBondReact : public Fix {
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 nvvec;
double *vvec; // per-atom vector to store variable constraint atom-style variable values
int **nxspecial, **onemol_nxspecial, **twomol_nxspecial; // full number of 1-4 neighbors
tagint **xspecial, **onemol_xspecial, **twomol_xspecial; // full 1-4 neighbor list
@ -151,14 +149,12 @@ class FixBondReact : public Fix {
// for all mega_gloves and global_mega_glove: first row is the ID of bond/react
tagint **local_mega_glove; // consolidation local of reaction instances
tagint **ghostly_mega_glove; // consolidation nonlocal of reaction instances
tagint *
*global_mega_glove; // consolidation (inter-processor) of gloves containing nonlocal atoms
tagint **global_mega_glove; // consolidation (inter-processor) of gloves containing nonlocal atoms
int *localsendlist; // indicates ghosts of other procs
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 *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(int);
@ -180,8 +176,10 @@ class FixBondReact : public Fix {
int check_constraints();
void get_IDcoords(int, int, double *);
double get_temperature(tagint **, int, int);
double custom_constraint(std::string);
double rxnfunction(std::string, std::string, std::string);
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(std::string); // evaulate expression for custom constraint
double rxnfunction(std::string, std::string, std::string); // eval rxn_sum and rxn_ave
int get_chirality(double[12]); // get handedness given an ordered set of coordinates
void open(char *);
@ -216,7 +214,11 @@ class FixBondReact : public Fix {
double par[MAXCONPAR];
std::string str;
};
Constraint **constraints;
int ncustomvars;
std::vector<std::string> customvarstrs;
int nvvec;
double **vvec; // per-atom vector to store variable constraint atom-style variable values
std::vector<std::vector<Constraint>> constraints;
// DEBUG

5
src/atom.cpp
View File

@ -2471,7 +2471,7 @@ This function is called, e.g. from :doc:`fix property/atom <fix_property_atom>`.
*/
int Atom::add_custom(const char *name, int flag, int cols)
{
int index;
int index = -1;
if ((flag == 0) && (cols == 0)) {
index = nivector;
@ -2511,7 +2511,8 @@ int Atom::add_custom(const char *name, int flag, int cols)
dcols = (int *) memory->srealloc(dcols,ndarray*sizeof(int),"atom:dcols");
dcols[index] = cols;
}
if (index < 0)
error->all(FLERR,"Invalid call to Atom::add_custom()");
return index;
}

5
src/domain.cpp
View File

@ -528,10 +528,11 @@ void Domain::reset_box()
void Domain::pbc()
{
int nlocal = atom->nlocal;
if (!nlocal) return;
int i;
imageint idim,otherdims;
double *lo,*hi,*period;
int nlocal = atom->nlocal;
double **x = atom->x;
double **v = atom->v;
int *mask = atom->mask;
@ -542,7 +543,7 @@ void Domain::pbc()
double *coord;
int n3 = 3*nlocal;
coord = &x[0][0]; // note: x is always initialized to at least one element.
coord = &x[0][0];
int flag = 0;
for (i = 0; i < n3; i++)
if (!std::isfinite(*coord++)) flag = 1;

2
src/dump.cpp
View File

@ -100,7 +100,7 @@ Dump::Dump(LAMMPS *lmp, int /*narg*/, char **arg) : Pointers(lmp)
maxsbuf = 0;
sbuf = nullptr;
maxpbc = 0;
maxpbc = -1;
xpbc = vpbc = nullptr;
imagepbc = nullptr;

230
src/fix_property_atom.cpp
View File

@ -31,7 +31,7 @@ enum{MOLECULE,CHARGE,RMASS,IVEC,DVEC,IARRAY,DARRAY};
FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
nvalue(0), style(nullptr), index(nullptr), astyle(nullptr)
nvalue(0), styles(nullptr), index(nullptr), astyle(nullptr)
{
if (narg < 4) error->all(FLERR,"Illegal fix property/atom command");
@ -40,7 +40,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
int iarg = 3;
nvalue = narg-iarg;
style = new int[nvalue];
styles = new int[nvalue];
cols = new int[nvalue];
index = new int[nvalue];
@ -58,7 +58,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
"already has molecule attribute");
if (molecule_flag)
error->all(FLERR,"Fix property/atom cannot specify mol twice");
style[nvalue] = MOLECULE;
styles[nvalue] = MOLECULE;
cols[nvalue] = 0;
atom->molecule_flag = molecule_flag = 1;
values_peratom++;
@ -69,7 +69,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Fix property/atom q when atom_style already has charge attribute");
if (q_flag)
error->all(FLERR,"Fix property/atom cannot specify q twice");
style[nvalue] = CHARGE;
styles[nvalue] = CHARGE;
cols[nvalue] = 0;
atom->q_flag = q_flag = 1;
values_peratom++;
@ -80,7 +80,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Fix property/atom rmass when atom_style already has rmass attribute");
if (rmass_flag)
error->all(FLERR,"Fix property/atom cannot specify rmass twice");
style[nvalue] = RMASS;
styles[nvalue] = RMASS;
cols[nvalue] = 0;
atom->rmass_flag = rmass_flag = 1;
values_peratom++;
@ -90,7 +90,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
// custom atom vector
} else if (utils::strmatch(arg[iarg],"^i_")) {
style[nvalue] = IVEC;
styles[nvalue] = IVEC;
int flag,ncols;
index[nvalue] = atom->find_custom(&arg[iarg][2],flag,ncols);
if (index[nvalue] >= 0)
@ -102,7 +102,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
iarg++;
} else if (utils::strmatch(arg[iarg],"^d_")) {
style[nvalue] = DVEC;
styles[nvalue] = DVEC;
int flag,ncols;
index[nvalue] = atom->find_custom(&arg[iarg][2],flag,ncols);
if (index[nvalue] >= 0)
@ -129,10 +129,10 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
if (arg[iarg][0] == 'i') {
which = 0;
style[nvalue] = IARRAY;
styles[nvalue] = IARRAY;
} else {
which = 1;
style[nvalue] = DARRAY;
styles[nvalue] = DARRAY;
}
index[nvalue] = atom->add_custom(&arg[iarg][3],which,ncols);
cols[nvalue] = ncols;
@ -165,9 +165,9 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
if (border == 0) {
int flag = 0;
for (int i = 0; i < nvalue; i++)
if (style[i] == MOLECULE
|| style[i] == CHARGE
|| style[i] == RMASS) flag = 1;
if (styles[i] == MOLECULE
|| styles[i] == CHARGE
|| styles[i] == RMASS) flag = 1;
if (flag && comm->me == 0)
error->warning(FLERR,"Fix property/atom mol or charge or rmass "
"w/out ghost communication");
@ -201,30 +201,30 @@ FixPropertyAtom::~FixPropertyAtom()
// set ptrs to a null pointer, so they no longer exist for Atom class
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE) {
if (styles[nv] == MOLECULE) {
atom->molecule_flag = 0;
memory->destroy(atom->molecule);
atom->molecule = nullptr;
} else if (style[nv] == CHARGE) {
} else if (styles[nv] == CHARGE) {
atom->q_flag = 0;
memory->destroy(atom->q);
atom->q = nullptr;
} else if (style[nv] == RMASS) {
} else if (styles[nv] == RMASS) {
atom->rmass_flag = 0;
memory->destroy(atom->rmass);
atom->rmass = nullptr;
} else if (style[nv] == IVEC) {
} else if (styles[nv] == IVEC) {
atom->remove_custom(index[nv],0,cols[nv]);
} else if (style[nv] == DVEC) {
} else if (styles[nv] == DVEC) {
atom->remove_custom(index[nv],1,cols[nv]);
} else if (style[nv] == IARRAY) {
} else if (styles[nv] == IARRAY) {
atom->remove_custom(index[nv],0,cols[nv]);
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
atom->remove_custom(index[nv],1,cols[nv]);
}
}
delete [] style;
delete [] styles;
delete [] cols;
delete [] index;
delete [] astyle;
@ -291,21 +291,21 @@ void FixPropertyAtom::read_data_section(char *keyword, int n, char *buf, tagint
if ((m = atom->map(itag)) >= 0) {
for (j = 0; j < nvalue; j++) {
if (style[j] == MOLECULE) {
if (styles[j] == MOLECULE) {
atom->molecule[m] = values.next_tagint();
} else if (style[j] == CHARGE) {
} else if (styles[j] == CHARGE) {
atom->q[m] = values.next_double();
} else if (style[j] == RMASS) {
} else if (styles[j] == RMASS) {
atom->rmass[m] = values.next_double();
} else if (style[j] == IVEC) {
} else if (styles[j] == IVEC) {
atom->ivector[index[j]][m] = values.next_int();
} else if (style[j] == DVEC) {
} else if (styles[j] == DVEC) {
atom->dvector[index[j]][m] = values.next_double();
} else if (style[j] == IARRAY) {
} else if (styles[j] == IARRAY) {
ncol = cols[j];
for (k = 0; k < ncol; k++)
atom->iarray[index[j]][m][k] = values.next_int();
} else if (style[j] == DARRAY) {
} else if (styles[j] == DARRAY) {
ncol = cols[j];
for (k = 0; k < ncol; k++)
atom->darray[index[j]][m][k] = values.next_double();
@ -365,34 +365,34 @@ void FixPropertyAtom::write_data_section_pack(int /*mth*/, double **buf)
int icol = 1;
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE) {
if (styles[nv] == MOLECULE) {
tagint *molecule = atom->molecule;
for (i = 0; i < nlocal; i++) buf[i][icol] = ubuf(molecule[i]).d;
icol++;
} else if (style[nv] == CHARGE) {
} else if (styles[nv] == CHARGE) {
double *q = atom->q;
for (i = 0; i < nlocal; i++) buf[i][icol] = q[i];
icol++;
} else if (style[nv] == RMASS) {
} else if (styles[nv] == RMASS) {
double *rmass = atom->rmass;
for (i = 0; i < nlocal; i++) buf[i][icol] = rmass[i];
icol++;
} else if (style[nv] == IVEC) {
} else if (styles[nv] == IVEC) {
int *ivec = atom->ivector[index[nv]];
for (i = 0; i < nlocal; i++) buf[i][icol] = ubuf(ivec[i]).d;
icol++;
} else if (style[nv] == DVEC) {
} else if (styles[nv] == DVEC) {
double *dvec = atom->dvector[index[nv]];
for (i = 0; i < nlocal; i++) buf[i][icol] = dvec[i];
icol++;
} else if (style[nv] == IARRAY) {
} else if (styles[nv] == IARRAY) {
int **iarray = atom->iarray[index[nv]];
ncol = cols[nv];
for (i = 0; i < nlocal; i++)
for (k = 0; k < ncol; k++)
buf[i][icol+k] = ubuf(iarray[i][k]).d;
icol += ncol;
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
double **darray = atom->darray[index[nv]];
ncol = cols[nv];
for (i = 0; i < nlocal; i++)
@ -412,19 +412,19 @@ void FixPropertyAtom::write_data_section_pack(int /*mth*/, double **buf)
void FixPropertyAtom::write_data_section_keyword(int /*mth*/, FILE *fp)
{
if (nvalue == 1 && style[0] == MOLECULE) fprintf(fp,"\nMolecules\n\n");
else if (nvalue == 1 && style[0] == CHARGE) fprintf(fp,"\nCharges\n\n");
if (nvalue == 1 && styles[0] == MOLECULE) fprintf(fp,"\nMolecules\n\n");
else if (nvalue == 1 && styles[0] == CHARGE) fprintf(fp,"\nCharges\n\n");
else {
fprintf(fp,"\n%s #",id);
// write column hint as comment
for (int i = 0; i < nvalue; ++i) {
if (style[i] == MOLECULE) fputs(" mol",fp);
else if (style[i] == CHARGE) fputs(" q",fp);
else if (style[i] == RMASS) fputs(" rmass",fp);
else if (style[i] == IVEC) fprintf(fp," i_%s", atom->ivname[index[i]]);
else if (style[i] == DVEC) fprintf(fp, " d_%s", atom->dvname[index[i]]);
else if (style[i] == IARRAY) fprintf(fp, " i_%s", atom->ianame[index[i]]);
else if (style[i] == DARRAY) fprintf(fp, " d_%s", atom->daname[index[i]]);
if (styles[i] == MOLECULE) fputs(" mol",fp);
else if (styles[i] == CHARGE) fputs(" q",fp);
else if (styles[i] == RMASS) fputs(" rmass",fp);
else if (styles[i] == IVEC) fprintf(fp," i_%s", atom->ivname[index[i]]);
else if (styles[i] == DVEC) fprintf(fp, " d_%s", atom->dvname[index[i]]);
else if (styles[i] == IARRAY) fprintf(fp, " i_%s", atom->ianame[index[i]]);
else if (styles[i] == DARRAY) fprintf(fp, " d_%s", atom->daname[index[i]]);
}
fputs("\n\n",fp);
}
@ -446,22 +446,22 @@ void FixPropertyAtom::write_data_section(int /*mth*/, FILE *fp,
fprintf(fp,TAGINT_FORMAT,(tagint) ubuf(buf[i][0]).i);
icol = 1;
for (nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE)
if (styles[nv] == MOLECULE)
fprintf(fp," " TAGINT_FORMAT,(tagint) ubuf(buf[i][icol++]).i);
else if (style[nv] == CHARGE)
else if (styles[nv] == CHARGE)
fprintf(fp," %g",buf[i][icol++]);
else if (style[nv] == RMASS)
else if (styles[nv] == RMASS)
fprintf(fp," %g",buf[i][icol++]);
else if (style[nv] == IVEC)
else if (styles[nv] == IVEC)
fprintf(fp," %d",(int) ubuf(buf[i][icol++]).i);
else if (style[nv] == DVEC)
else if (styles[nv] == DVEC)
fprintf(fp," %g",buf[i][icol++]);
else if (style[nv] == IARRAY) {
else if (styles[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
fprintf(fp," %d",(int) ubuf(buf[i][icol+k]).i);
icol += ncol;
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
fprintf(fp," %g",buf[i][icol+k]);
@ -480,13 +480,13 @@ double FixPropertyAtom::memory_usage()
{
double bytes = 0.0;
for (int m = 0; m < nvalue; m++) {
if (style[m] == MOLECULE) bytes = atom->nmax * sizeof(tagint);
else if (style[m] == CHARGE) bytes = atom->nmax * sizeof(double);
else if (style[m] == RMASS) bytes = atom->nmax * sizeof(double);
else if (style[m] == IVEC) bytes = atom->nmax * sizeof(int);
else if (style[m] == DVEC) bytes = atom->nmax * sizeof(double);
else if (style[m] == IARRAY) bytes = atom->nmax * cols[m] * sizeof(int);
else if (style[m] == DARRAY) bytes = atom->nmax * cols[m] * sizeof(double);
if (styles[m] == MOLECULE) bytes = atom->nmax * sizeof(tagint);
else if (styles[m] == CHARGE) bytes = atom->nmax * sizeof(double);
else if (styles[m] == RMASS) bytes = atom->nmax * sizeof(double);
else if (styles[m] == IVEC) bytes = atom->nmax * sizeof(int);
else if (styles[m] == DVEC) bytes = atom->nmax * sizeof(double);
else if (styles[m] == IARRAY) bytes = (size_t) atom->nmax * cols[m] * sizeof(int);
else if (styles[m] == DARRAY) bytes = (size_t) atom->nmax * cols[m] * sizeof(double);
}
return bytes;
}
@ -501,33 +501,33 @@ double FixPropertyAtom::memory_usage()
void FixPropertyAtom::grow_arrays(int nmax)
{
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE) {
if (styles[nv] == MOLECULE) {
memory->grow(atom->molecule,nmax,"atom:molecule");
size_t nbytes = (nmax-nmax_old) * sizeof(tagint);
memset(&atom->molecule[nmax_old],0,nbytes);
} else if (style[nv] == CHARGE) {
} else if (styles[nv] == CHARGE) {
memory->grow(atom->q,nmax,"atom:q");
size_t nbytes = (nmax-nmax_old) * sizeof(double);
memset(&atom->q[nmax_old],0,nbytes);
} else if (style[nv] == RMASS) {
} else if (styles[nv] == RMASS) {
memory->grow(atom->rmass,nmax,"atom:rmass");
size_t nbytes = (nmax-nmax_old) * sizeof(double);
memset(&atom->rmass[nmax_old],0,nbytes);
} else if (style[nv] == IVEC) {
} else if (styles[nv] == IVEC) {
memory->grow(atom->ivector[index[nv]],nmax,"atom:ivector");
size_t nbytes = (nmax-nmax_old) * sizeof(int);
memset(&atom->ivector[index[nv]][nmax_old],0,nbytes);
} else if (style[nv] == DVEC) {
} else if (styles[nv] == DVEC) {
memory->grow(atom->dvector[index[nv]],nmax,"atom:dvector");
size_t nbytes = (nmax-nmax_old) * sizeof(double);
memset(&atom->dvector[index[nv]][nmax_old],0,nbytes);
} else if (style[nv] == IARRAY) {
} else if (styles[nv] == IARRAY) {
memory->grow(atom->iarray[index[nv]],nmax,cols[nv],"atom:iarray");
size_t nbytes = (nmax-nmax_old) * cols[nv] * sizeof(int);
size_t nbytes = (size_t) (nmax-nmax_old) * cols[nv] * sizeof(int);
if (nbytes) memset(&atom->iarray[index[nv]][nmax_old][0],0,nbytes);
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
memory->grow(atom->darray[index[nv]],nmax,cols[nv],"atom:darray");
size_t nbytes = (nmax-nmax_old) * cols[nv] * sizeof(double);
size_t nbytes = (size_t) (nmax-nmax_old) * cols[nv] * sizeof(double);
if (nbytes) memset(&atom->darray[index[nv]][nmax_old][0],0,nbytes);
}
}
@ -544,21 +544,21 @@ void FixPropertyAtom::copy_arrays(int i, int j, int /*delflag*/)
int k,ncol;
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE)
if (styles[nv] == MOLECULE)
atom->molecule[j] = atom->molecule[i];
else if (style[nv] == CHARGE)
else if (styles[nv] == CHARGE)
atom->q[j] = atom->q[i];
else if (style[nv] == RMASS)
else if (styles[nv] == RMASS)
atom->rmass[j] = atom->rmass[i];
else if (style[nv] == IVEC)
else if (styles[nv] == IVEC)
atom->ivector[index[nv]][j] = atom->ivector[index[nv]][i];
else if (style[nv] == DVEC)
else if (styles[nv] == DVEC)
atom->dvector[index[nv]][j] = atom->dvector[index[nv]][i];
else if (style[nv] == IARRAY) {
else if (styles[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->iarray[index[nv]][j][k] = atom->iarray[index[nv]][i][k];
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->darray[index[nv]][j][k] = atom->darray[index[nv]][i][k];
@ -576,37 +576,37 @@ int FixPropertyAtom::pack_border(int n, int *list, double *buf)
int m = 0;
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE) {
if (styles[nv] == MOLECULE) {
tagint *molecule = atom->molecule;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = ubuf(molecule[j]).d;
}
} else if (style[nv] == CHARGE) {
} else if (styles[nv] == CHARGE) {
double *q = atom->q;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = q[j];
}
} else if (style[nv] == RMASS) {
} else if (styles[nv] == RMASS) {
double *rmass = atom->rmass;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = rmass[j];
}
} else if (style[nv] == IVEC) {
} else if (styles[nv] == IVEC) {
int *ivector = atom->ivector[index[nv]];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = ubuf(ivector[j]).d;
}
} else if (style[nv] == DVEC) {
} else if (styles[nv] == DVEC) {
double *dvector = atom->dvector[index[nv]];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = dvector[j];
}
} else if (style[nv] == IARRAY) {
} else if (styles[nv] == IARRAY) {
int **iarray = atom->iarray[index[nv]];
ncol = cols[nv];
for (i = 0; i < n; i++) {
@ -614,7 +614,7 @@ int FixPropertyAtom::pack_border(int n, int *list, double *buf)
for (k = 0; k < ncol; k++)
buf[m++] = ubuf(iarray[j][k]).d;
}
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
double **darray = atom->darray[index[nv]];
ncol = cols[nv];
for (i = 0; i < n; i++) {
@ -638,39 +638,39 @@ int FixPropertyAtom::unpack_border(int n, int first, double *buf)
int m = 0;
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE) {
if (styles[nv] == MOLECULE) {
tagint *molecule = atom->molecule;
last = first + n;
for (i = first; i < last; i++)
molecule[i] = (tagint) ubuf(buf[m++]).i;
} else if (style[nv] == CHARGE) {
} else if (styles[nv] == CHARGE) {
double *q = atom->q;
last = first + n;
for (i = first; i < last; i++)
q[i] = buf[m++];
} else if (style[nv] == RMASS) {
} else if (styles[nv] == RMASS) {
double *rmass = atom->rmass;
last = first + n;
for (i = first; i < last; i++)
rmass[i] = buf[m++];
} else if (style[nv] == IVEC) {
} else if (styles[nv] == IVEC) {
int *ivector = atom->ivector[index[nv]];
last = first + n;
for (i = first; i < last; i++)
ivector[i] = (int) ubuf(buf[m++]).i;
} else if (style[nv] == DVEC) {
} else if (styles[nv] == DVEC) {
double *dvector = atom->dvector[index[nv]];
last = first + n;
for (i = first; i < last; i++)
dvector[i] = buf[m++];
} else if (style[nv] == IARRAY) {
} else if (styles[nv] == IARRAY) {
int **iarray = atom->iarray[index[nv]];
ncol = cols[nv];
last = first + n;
for (i = first; i < last; i++)
for (k = 0; k < ncol; k++)
iarray[i][k] = (int) ubuf(buf[m++]).i;
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
double **darray = atom->darray[index[nv]];
ncol = cols[nv];
last = first + n;
@ -693,16 +693,16 @@ int FixPropertyAtom::pack_exchange(int i, double *buf)
int m = 0;
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d;
else if (style[nv] == CHARGE) buf[m++] = atom->q[i];
else if (style[nv] == RMASS) buf[m++] = atom->rmass[i];
else if (style[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
else if (style[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i];
else if (style[nv] == IARRAY) {
if (styles[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d;
else if (styles[nv] == CHARGE) buf[m++] = atom->q[i];
else if (styles[nv] == RMASS) buf[m++] = atom->rmass[i];
else if (styles[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
else if (styles[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i];
else if (styles[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d;
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
buf[m++] = atom->darray[index[nv]][i][k];
@ -722,21 +722,21 @@ int FixPropertyAtom::unpack_exchange(int nlocal, double *buf)
int m = 0;
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE)
if (styles[nv] == MOLECULE)
atom->molecule[nlocal] = (tagint) ubuf(buf[m++]).i;
else if (style[nv] == CHARGE)
else if (styles[nv] == CHARGE)
atom->q[nlocal] = buf[m++];
else if (style[nv] == RMASS)
else if (styles[nv] == RMASS)
atom->rmass[nlocal] = buf[m++];
else if (style[nv] == IVEC)
else if (styles[nv] == IVEC)
atom->ivector[index[nv]][nlocal] = (int) ubuf(buf[m++]).i;
else if (style[nv] == DVEC)
else if (styles[nv] == DVEC)
atom->dvector[index[nv]][nlocal] = buf[m++];
else if (style[nv] == IARRAY) {
else if (styles[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->iarray[index[nv]][nlocal][k] = (int) ubuf(buf[m++]).i;
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->darray[index[nv]][nlocal][k] = buf[m++];
@ -760,16 +760,16 @@ int FixPropertyAtom::pack_restart(int i, double *buf)
int m = 1;
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d;
else if (style[nv] == CHARGE) buf[m++] = atom->q[i];
else if (style[nv] == RMASS) buf[m++] = atom->rmass[i];
else if (style[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
else if (style[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i];
else if (style[nv] == IARRAY) {
if (styles[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d;
else if (styles[nv] == CHARGE) buf[m++] = atom->q[i];
else if (styles[nv] == RMASS) buf[m++] = atom->rmass[i];
else if (styles[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d;
else if (styles[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i];
else if (styles[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d;
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
buf[m++] = atom->darray[index[nv]][i][k];
@ -796,21 +796,21 @@ void FixPropertyAtom::unpack_restart(int nlocal, int nth)
m++;
for (int nv = 0; nv < nvalue; nv++) {
if (style[nv] == MOLECULE)
if (styles[nv] == MOLECULE)
atom->molecule[nlocal] = (tagint) ubuf(extra[nlocal][m++]).i;
else if (style[nv] == CHARGE)
else if (styles[nv] == CHARGE)
atom->q[nlocal] = extra[nlocal][m++];
else if (style[nv] == RMASS)
else if (styles[nv] == RMASS)
atom->rmass[nlocal] = extra[nlocal][m++];
else if (style[nv] == IVEC)
else if (styles[nv] == IVEC)
atom->ivector[index[nv]][nlocal] = (int) ubuf(extra[nlocal][m++]).i;
else if (style[nv] == DVEC)
else if (styles[nv] == DVEC)
atom->dvector[index[nv]][nlocal] = extra[nlocal][m++];
else if (style[nv] == IARRAY) {
else if (styles[nv] == IARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->iarray[index[nv]][nlocal][k] = (int) ubuf(extra[nlocal][m++]).i;
} else if (style[nv] == DARRAY) {
} else if (styles[nv] == DARRAY) {
ncol = cols[nv];
for (k = 0; k < ncol; k++)
atom->darray[index[nv]][nlocal][k] = extra[nlocal][m++];

2
src/fix_property_atom.h
View File

@ -53,7 +53,7 @@ class FixPropertyAtom : public Fix {
protected:
int nvalue, border;
int molecule_flag, q_flag, rmass_flag; // flags for specific fields
int *style; // style of each value, see enum
int *styles; // style of each value, see enum
int *index; // indices into atom custom data structs
int *cols; // columns per value, for arrays
char *astyle; // atom style at instantiation

6
src/info.cpp
View File

@ -1356,8 +1356,10 @@ std::string Info::get_openmp_info()
// Supported OpenMP version corresponds to the release date of the
// specifications as posted at https://www.openmp.org/specifications/
#if _OPENMP > 201811
return "OpenMP newer than version 5.0";
#if _OPENMP > 202011
return "OpenMP newer than version 5.1";
#elif _OPENMP == 202011
return "OpenMP 5.1";
#elif _OPENMP == 201811
return "OpenMP 5.0";
#elif _OPENMP == 201611

2
src/library.cpp
View File

@ -2712,7 +2712,7 @@ Below is a brief C code demonstrating accessing this collected bond information.
#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <stdint.h>
#include "library.h"
int main(int argc, char **argv)

2
src/library.h
View File

@ -35,7 +35,7 @@
#endif
#if defined(LAMMPS_BIGBIG) || defined(LAMMPS_SMALLBIG)
#include <inttypes.h> /* for int64_t */
#include <stdint.h> /* for int64_t */
#endif
/** Data type constants for extracting data from atoms, computes and fixes

110
src/memory.h
View File

@ -27,7 +27,7 @@ class Memory : protected Pointers {
void sfree(void *);
void fail(const char *);
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create/grow/destroy vecs and multidim arrays with contiguous memory blocks
only use with primitive data types, e.g. 1d vec of ints, 2d array of doubles
fail() prevents use with pointers,
@ -36,12 +36,15 @@ class Memory : protected Pointers {
for these other cases, use smalloc/srealloc/sfree directly
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 1d array
------------------------------------------------------------------------- */
template <typename TYPE> TYPE *create(TYPE *&array, int n, const char *name)
{
// POSSIBLE future change
//if (n <= 0) return nullptr;
bigint nbytes = ((bigint) sizeof(TYPE)) * n;
array = (TYPE *) smalloc(nbytes, name);
return array;
@ -53,12 +56,18 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
grow or shrink 1d array
------------------------------------------------------------------------- */
template <typename TYPE> TYPE *grow(TYPE *&array, int n, const char *name)
{
// POSSIBLE future change
//if (n <= 0) {
// destroy(array);
// return nullptr;
// }
if (array == nullptr) return create(array, n, name);
bigint nbytes = ((bigint) sizeof(TYPE)) * n;
@ -72,7 +81,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
destroy a 1d array
------------------------------------------------------------------------- */
@ -82,13 +91,16 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 1d array with index from nlo to nhi inclusive
cannot grow it
------------------------------------------------------------------------- */
template <typename TYPE> TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name)
{
// POSSIBLE future change
// if (nlo > nhi) return nullptr;
bigint nbytes = ((bigint) sizeof(TYPE)) * (nhi - nlo + 1);
array = (TYPE *) smalloc(nbytes, name);
array -= nlo;
@ -102,7 +114,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
destroy a 1d array with index offset
------------------------------------------------------------------------- */
@ -112,12 +124,15 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 2d array
------------------------------------------------------------------------- */
template <typename TYPE> TYPE **create(TYPE **&array, int n1, int n2, const char *name)
{
// POSSIBLE future change
//if (n1 <= 0 || n2 <= 0) return nullptr;
bigint nbytes = ((bigint) sizeof(TYPE)) * n1 * n2;
TYPE *data = (TYPE *) smalloc(nbytes, name);
nbytes = ((bigint) sizeof(TYPE *)) * n1;
@ -138,13 +153,19 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
grow or shrink 1st dim of a 2d array
last dim must stay the same
------------------------------------------------------------------------- */
template <typename TYPE> TYPE **grow(TYPE **&array, int n1, int n2, const char *name)
{
// POSSIBLE future change
//if (n1 <= 0 || n2 <= 0) {
// destroy(array);
// return nullptr;
// }
if (array == nullptr) return create(array, n1, n2, name);
bigint nbytes = ((bigint) sizeof(TYPE)) * n1 * n2;
@ -167,7 +188,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
destroy a 2d array
------------------------------------------------------------------------- */
@ -179,12 +200,15 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 2d array with a ragged 2nd dimension
------------------------------------------------------------------------- */
template <typename TYPE> TYPE **create_ragged(TYPE **&array, int n1, int *n2, const char *name)
{
// POSSIBLE future change
//if (n1 <= 0) return nullptr;
bigint n2sum = 0;
for (int i = 0; i < n1; i++) n2sum += n2[i];
@ -208,7 +232,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 2d array with 2nd index from n2lo to n2hi inclusive
cannot grow it
------------------------------------------------------------------------- */
@ -216,6 +240,9 @@ class Memory : protected Pointers {
template <typename TYPE>
TYPE **create2d_offset(TYPE **&array, int n1, int n2lo, int n2hi, const char *name)
{
// POSSIBLE future change
//if (n1 <= 0 || n2lo > n2hi) return nullptr;
int n2 = n2hi - n2lo + 1;
create(array, n1, n2, name);
for (int i = 0; i < n1; i++) array[i] -= n2lo;
@ -230,7 +257,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
destroy a 2d array with 2nd index offset
------------------------------------------------------------------------- */
@ -242,12 +269,15 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 3d array
------------------------------------------------------------------------- */
template <typename TYPE> TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name)
{
// POSSIBLE future change
//if (n1 <= 0 || n2 <= 0 || n3 <= 0) return nullptr;
bigint nbytes = ((bigint) sizeof(TYPE)) * n1 * n2 * n3;
TYPE *data = (TYPE *) smalloc(nbytes, name);
nbytes = ((bigint) sizeof(TYPE *)) * n1 * n2;
@ -276,13 +306,19 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
grow or shrink 1st dim of a 3d array
last 2 dims must stay the same
------------------------------------------------------------------------- */
template <typename TYPE> TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name)
{
// POSSIBLE future change
//if (n1 <= 0 || n2 <= 0 || n3 <= 0) {
// destroy(array);
// return nullptr;
//};
if (array == nullptr) return create(array, n1, n2, n3, name);
bigint nbytes = ((bigint) sizeof(TYPE)) * n1 * n2 * n3;
@ -313,7 +349,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
destroy a 3d array
------------------------------------------------------------------------- */
@ -326,7 +362,7 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 3d array with 1st index from n1lo to n1hi inclusive
cannot grow it
------------------------------------------------------------------------- */
@ -334,6 +370,8 @@ class Memory : protected Pointers {
template <typename TYPE>
TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, int n2, int n3, const char *name)
{
if (n1lo > n1hi || n2 <= 0 || n3 <= 0) return nullptr;
int n1 = n1hi - n1lo + 1;
create(array, n1, n2, n3, name);
array -= n1lo;
@ -348,7 +386,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
free a 3d array with 1st index offset
------------------------------------------------------------------------- */
@ -361,7 +399,7 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 3d array with
1st index from n1lo to n1hi inclusive,
2nd index from n2lo to n2hi inclusive,
@ -373,6 +411,8 @@ class Memory : protected Pointers {
TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, int n2lo, int n2hi, int n3lo,
int n3hi, const char *name)
{
if (n1lo > n1hi || n2lo > n2hi || n3lo > n3hi) return nullptr;
int n1 = n1hi - n1lo + 1;
int n2 = n2hi - n2lo + 1;
int n3 = n3hi - n3lo + 1;
@ -393,7 +433,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
free a 3d array with all 3 indices offset
------------------------------------------------------------------------- */
@ -407,13 +447,16 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 4d array
------------------------------------------------------------------------- */
template <typename TYPE>
TYPE ****create(TYPE ****&array, int n1, int n2, int n3, int n4, const char *name)
{
// POSSIBLE future change
//if (n1 <= 0 || n2 <= 0 || n3 <= 0 || n4 <= 0) return nullptr;
bigint nbytes = ((bigint) sizeof(TYPE)) * n1 * n2 * n3 * n4;
TYPE *data = (TYPE *) smalloc(nbytes, name);
nbytes = ((bigint) sizeof(TYPE *)) * n1 * n2 * n3;
@ -451,7 +494,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
grow or shrink 1st dim of a 4d array
last 3 dims must stay the same
------------------------------------------------------------------------- */
@ -459,6 +502,12 @@ class Memory : protected Pointers {
template <typename TYPE>
TYPE ****grow(TYPE ****&array, int n1, int n2, int n3, int n4, const char *name)
{
// POSSIBLE future change
//if (n1 <= 0 || n2 <= 0 || n3 <= 0 || n4 <= 0) {
// destroy(array);
// return nullptr;
// }
if (array == nullptr) return create(array, n1, n2, n3, n4, name);
bigint nbytes = ((bigint) sizeof(TYPE)) * n1 * n2 * n3 * n4;
@ -498,7 +547,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
destroy a 4d array
------------------------------------------------------------------------- */
@ -512,7 +561,7 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 4d array with indices
2nd index from n2lo to n2hi inclusive
3rd index from n3lo to n3hi inclusive
@ -524,6 +573,8 @@ class Memory : protected Pointers {
TYPE ****create4d_offset(TYPE ****&array, int n1, int n2lo, int n2hi, int n3lo, int n3hi,
int n4lo, int n4hi, const char *name)
{
if (n1 <= 0 || n2lo > n2hi || n3lo > n3hi || n4lo > n4hi) return nullptr;
int n2 = n2hi - n2lo + 1;
int n3 = n3hi - n3lo + 1;
int n4 = n4hi - n4lo + 1;
@ -545,8 +596,8 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
free a 4d array with indices 2,3, and 4 offset
/* ----------------------------------------------------------------------
free a 4d array with indices 2,3,4 offset
------------------------------------------------------------------------- */
template <typename TYPE>
@ -560,13 +611,16 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create a 5d array
------------------------------------------------------------------------- */
template <typename TYPE>
TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4, int n5, const char *name)
{
// POSSIBLE future change
//if (n1 <= 0 || n2 <= 0 || n3 <= 0 || n4 <= 0 || n5 <= 0) return nullptr;
bigint nbytes = ((bigint) sizeof(TYPE)) * n1 * n2 * n3 * n4 * n5;
TYPE *data = (TYPE *) smalloc(nbytes, name);
nbytes = ((bigint) sizeof(TYPE *)) * n1 * n2 * n3 * n4;
@ -611,7 +665,7 @@ class Memory : protected Pointers {
return nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
destroy a 5d array
------------------------------------------------------------------------- */
@ -626,7 +680,7 @@ class Memory : protected Pointers {
array = nullptr;
}
/* ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
memory usage of arrays, including pointers
------------------------------------------------------------------------- */

6
src/molecule.cpp
View File

@ -148,13 +148,15 @@ Molecule::Molecule(LAMMPS *lmp, int narg, char **arg, int &index) :
if (me == 0)
utils::logmesg(lmp,"Read molecule template {}:\n {} molecules\n"
" {} fragments\n"
" {} atoms with max type {}\n"
" {} bonds with max type {}\n"
" {} angles with max type {}\n"
" {} dihedrals with max type {}\n"
" {} impropers with max type {}\n", id,nmolecules,
natoms,ntypes,nbonds,nbondtypes,nangles,nangletypes,
ndihedrals,ndihedraltypes,nimpropers,nimpropertypes);
nfragments,natoms,ntypes,nbonds,nbondtypes,nangles,
nangletypes,ndihedrals,ndihedraltypes,nimpropers,
nimpropertypes);
}
/* ---------------------------------------------------------------------- */

605
src/pair_table.cpp
View File

File diff suppressed because it is too large Load Diff