From c6bbdcfbf2af5d717f244495cd3cb3e2b35286e5 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Mon, 24 Oct 2011 20:28:02 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7184
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/ASPHERE/pair_line_lj.cpp | 1 +
src/ASPHERE/pair_line_lj.h | 4 +-
src/ASPHERE/pair_resquared.cpp | 2 +
src/ASPHERE/pair_tri_lj.cpp | 3 +-
src/ASPHERE/pair_tri_lj.h | 4 +-
src/CLASS2/angle_class2.h | 6 +--
src/CLASS2/bond_class2.h | 6 +--
src/CLASS2/improper_class2.h | 6 +--
src/GPU/fix_gpu.cpp | 9 ++--
src/GRANULAR/fix_pour.cpp | 8 +++-
src/GRANULAR/fix_wall_gran.h | 8 ++--
src/KSPACE/ewald.cpp | 2 +-
src/KSPACE/ewald.h | 10 ++--
src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp | 3 +-
src/KSPACE/pppm.cpp | 32 +++++++++----
src/KSPACE/pppm_cg.cpp | 2 +-
src/MANYBODY/pair_adp.cpp | 3 +-
src/MOLECULE/angle_charmm.h | 6 +--
src/MOLECULE/angle_cosine.h | 6 +--
src/MOLECULE/angle_cosine_delta.h | 2 +-
src/MOLECULE/angle_cosine_periodic.h | 6 +--
src/MOLECULE/angle_harmonic.h | 6 +--
src/MOLECULE/angle_table.h | 6 +--
src/MOLECULE/bond_fene.h | 6 +--
src/MOLECULE/bond_fene_expand.h | 6 +--
src/MOLECULE/bond_harmonic.h | 6 +--
src/MOLECULE/bond_morse.h | 6 +--
src/MOLECULE/bond_nonlinear.h | 6 +--
src/MOLECULE/bond_quartic.h | 6 +--
src/MOLECULE/bond_table.h | 6 +--
src/MOLECULE/improper_cvff.h | 6 +--
src/MOLECULE/improper_harmonic.h | 6 +--
src/MOLECULE/improper_umbrella.h | 6 +--
src/OPT/pair_lj_cut_coul_long_tip4p_opt.cpp | 2 +-
src/SHOCK/fix_wall_piston.cpp | 6 ++-
src/USER-MISC/angle_cosine_shift.h | 4 +-
src/USER-MISC/angle_cosine_shift_exp.cpp | 2 +-
src/USER-MISC/angle_cosine_shift_exp.h | 6 +--
src/USER-MISC/bond_harmonic_shift.h | 6 +--
src/USER-MISC/bond_harmonic_shift_cut.h | 6 +--
src/angle.cpp | 9 ++--
src/angle.h | 1 +
src/atom_vec_line.cpp | 11 +++--
src/atom_vec_tri.cpp | 11 +++--
src/bond.cpp | 9 ++--
src/bond.h | 1 +
src/dihedral.cpp | 10 ++--
src/domain.cpp | 2 +-
src/domain.h | 19 ++++----
src/force.cpp | 51 ++++++++++++++++++---
src/force.h | 3 +-
src/improper.cpp | 9 ++--
src/improper.h | 1 +
src/input.cpp | 5 +-
src/integrate.cpp | 1 +
src/integrate.h | 1 +
src/kspace.h | 1 +
src/math_const.h | 2 +
src/min.cpp | 11 +++--
src/min.h | 1 +
src/neigh_full.cpp | 32 ++++++-------
src/neigh_list.cpp | 9 ++--
src/neigh_request.cpp | 6 +++
src/neigh_request.h | 4 ++
src/neighbor.cpp | 5 +-
src/neighbor.h | 2 +-
src/pair.cpp | 8 ++--
src/pair.h | 1 +
src/respa.cpp | 2 +
src/thermo.h | 6 +--
src/verlet.cpp | 2 +
71 files changed, 279 insertions(+), 190 deletions(-)
diff --git a/src/ASPHERE/pair_line_lj.cpp b/src/ASPHERE/pair_line_lj.cpp
index ca50168be5..2f3f602ba7 100644
--- a/src/ASPHERE/pair_line_lj.cpp
+++ b/src/ASPHERE/pair_line_lj.cpp
@@ -96,6 +96,7 @@ void PairLineLJ::compute(int eflag, int vflag)
// grow discrete list if necessary and initialize
if (nall > nmax) {
+ nmax = nall;
memory->destroy(dnum);
memory->destroy(dfirst);
memory->create(dnum,nall,"pair:dnum");
diff --git a/src/ASPHERE/pair_line_lj.h b/src/ASPHERE/pair_line_lj.h
index b30140acea..4b60c80eaf 100644
--- a/src/ASPHERE/pair_line_lj.h
+++ b/src/ASPHERE/pair_line_lj.h
@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
class PairLineLJ : public Pair {
public:
PairLineLJ(class LAMMPS *);
- ~PairLineLJ();
- void compute(int, int);
+ virtual ~PairLineLJ();
+ virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
diff --git a/src/ASPHERE/pair_resquared.cpp b/src/ASPHERE/pair_resquared.cpp
index 13e57fa64e..9fb782b681 100755
--- a/src/ASPHERE/pair_resquared.cpp
+++ b/src/ASPHERE/pair_resquared.cpp
@@ -133,6 +133,8 @@ void PairRESquared::compute(int eflag, int vflag)
// compute if less than cutoff
if (rsq < cutsq[itype][jtype]) {
+ fforce[0] = fforce[1] = fforce[2] = 0.0;
+
switch (form[itype][jtype]) {
case SPHERE_SPHERE:
diff --git a/src/ASPHERE/pair_tri_lj.cpp b/src/ASPHERE/pair_tri_lj.cpp
index 0ecf19aee5..dfcdd819ef 100644
--- a/src/ASPHERE/pair_tri_lj.cpp
+++ b/src/ASPHERE/pair_tri_lj.cpp
@@ -100,6 +100,7 @@ void PairTriLJ::compute(int eflag, int vflag)
// grow discrete list if necessary and initialize
if (nall > nmax) {
+ nmax = nall;
memory->destroy(dnum);
memory->destroy(dfirst);
memory->create(dnum,nall,"pair:dnum");
@@ -118,7 +119,7 @@ void PairTriLJ::compute(int eflag, int vflag)
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
-
+
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
diff --git a/src/ASPHERE/pair_tri_lj.h b/src/ASPHERE/pair_tri_lj.h
index a8cf2bf721..fb0028cce0 100644
--- a/src/ASPHERE/pair_tri_lj.h
+++ b/src/ASPHERE/pair_tri_lj.h
@@ -27,8 +27,8 @@ namespace LAMMPS_NS {
class PairTriLJ : public Pair {
public:
PairTriLJ(class LAMMPS *);
- ~PairTriLJ();
- void compute(int, int);
+ virtual ~PairTriLJ();
+ virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double init_one(int, int);
diff --git a/src/CLASS2/angle_class2.h b/src/CLASS2/angle_class2.h
index cbea3cfe7e..f24ecce524 100644
--- a/src/CLASS2/angle_class2.h
+++ b/src/CLASS2/angle_class2.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class AngleClass2 : public Angle {
public:
AngleClass2(class LAMMPS *);
- ~AngleClass2();
- void compute(int, int);
+ virtual ~AngleClass2();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_angle(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, int, int, int);
- private:
+ protected:
double *theta0,*k2,*k3,*k4;
double *bb_k,*bb_r1,*bb_r2;
double *ba_k1,*ba_k2,*ba_r1,*ba_r2;
diff --git a/src/CLASS2/bond_class2.h b/src/CLASS2/bond_class2.h
index 2a5a42c648..3c201e923f 100644
--- a/src/CLASS2/bond_class2.h
+++ b/src/CLASS2/bond_class2.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class BondClass2 : public Bond {
public:
BondClass2(class LAMMPS *);
- ~BondClass2();
- void compute(int, int);
+ virtual ~BondClass2();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
double *r0,*k2,*k3,*k4;
void allocate();
diff --git a/src/CLASS2/improper_class2.h b/src/CLASS2/improper_class2.h
index 9a86e9cd86..cbee3ac810 100644
--- a/src/CLASS2/improper_class2.h
+++ b/src/CLASS2/improper_class2.h
@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
class ImproperClass2 : public Improper {
public:
ImproperClass2(class LAMMPS *);
- ~ImproperClass2();
- void compute(int, int);
+ virtual ~ImproperClass2();
+ virtual void compute(int, int);
void coeff(int, char **);
void write_restart(FILE *);
void read_restart(FILE *);
- private:
+ protected:
double *k0,*chi0;
double *aa_k1,*aa_k2,*aa_k3,*aa_theta0_1,*aa_theta0_2,*aa_theta0_3;
int *setflag_i,*setflag_aa;
diff --git a/src/GPU/fix_gpu.cpp b/src/GPU/fix_gpu.cpp
index a444986701..eec1391bfc 100644
--- a/src/GPU/fix_gpu.cpp
+++ b/src/GPU/fix_gpu.cpp
@@ -43,7 +43,8 @@ extern double lmp_gpu_forces(double **f, double **tor, double *eatom,
FixGPU::FixGPU(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
- if (lmp->cuda) error->all(FLERR,"Cannot use fix GPU with USER-CUDA mode enabled");
+ if (lmp->cuda)
+ error->all(FLERR,"Cannot use fix GPU with USER-CUDA mode enabled");
if (narg < 7) error->all(FLERR,"Illegal fix GPU command");
if (strcmp(arg[1],"all") != 0) error->all(FLERR,"Illegal fix GPU command");
@@ -113,10 +114,8 @@ int FixGPU::setmask()
void FixGPU::init()
{
- // Can only have 1 gpu fix that must be the first fix for a run
- if ((void*)modify->fix[0] != (void*)this)
- error->all(FLERR,"GPU is not the first fix for this run");
- // Hybrid cannot be used with force/neigh option
+ // hybrid cannot be used with force/neigh option
+
if (_gpu_mode == GPU_NEIGH)
if (force->pair_match("hybrid",1) != NULL ||
force->pair_match("hybrid/overlay",1) != NULL)
diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp
index 03a1c3a62a..fe0ce40f31 100644
--- a/src/GRANULAR/fix_pour.cpp
+++ b/src/GRANULAR/fix_pour.cpp
@@ -168,8 +168,10 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
// assume grav = -magnitude at this point, enforce in init()
int ifix;
- for (ifix = 0; ifix < modify->nfix; ifix++)
+ for (ifix = 0; ifix < modify->nfix; ifix++) {
if (strcmp(modify->fix[ifix]->style,"gravity") == 0) break;
+ if (strcmp(modify->fix[ifix]->style,"gravity/omp") == 0) break;
+ }
if (ifix == modify->nfix)
error->all(FLERR,"No fix gravity defined for fix pour");
grav = - ((FixGravity *) modify->fix[ifix])->magnitude * force->ftm2v;
@@ -267,8 +269,10 @@ void FixPour::init()
// else insertion cannot work
int ifix;
- for (ifix = 0; ifix < modify->nfix; ifix++)
+ for (ifix = 0; ifix < modify->nfix; ifix++) {
if (strcmp(modify->fix[ifix]->style,"gravity") == 0) break;
+ if (strcmp(modify->fix[ifix]->style,"gravity/omp") == 0) break;
+ }
if (ifix == modify->nfix)
error->all(FLERR,"No fix gravity defined for fix pour");
diff --git a/src/GRANULAR/fix_wall_gran.h b/src/GRANULAR/fix_wall_gran.h
index ae260e478e..343511c860 100644
--- a/src/GRANULAR/fix_wall_gran.h
+++ b/src/GRANULAR/fix_wall_gran.h
@@ -27,12 +27,12 @@ namespace LAMMPS_NS {
class FixWallGran : public Fix {
public:
FixWallGran(class LAMMPS *, int, char **);
- ~FixWallGran();
+ virtual ~FixWallGran();
int setmask();
void init();
void setup(int);
- void post_force(int);
- void post_force_respa(int, int, int);
+ virtual void post_force(int);
+ virtual void post_force_respa(int, int, int);
double memory_usage();
void grow_arrays(int);
@@ -46,7 +46,7 @@ class FixWallGran : public Fix {
int maxsize_restart();
void reset_dt();
- private:
+ protected:
int wallstyle,pairstyle,wiggle,wshear,axis;
double kn,kt,gamman,gammat,xmu;
double lo,hi,cylradius;
diff --git a/src/KSPACE/ewald.cpp b/src/KSPACE/ewald.cpp
index edd030b55b..b6cff3d932 100644
--- a/src/KSPACE/ewald.cpp
+++ b/src/KSPACE/ewald.cpp
@@ -282,7 +282,7 @@ void Ewald::compute(int eflag, int vflag)
for (k = 0; k < kcount; k++)
energy += ug[k] * (sfacrl_all[k]*sfacrl_all[k] +
sfacim_all[k]*sfacim_all[k]);
- energy -= g_ewald*qsqsum/1.772453851 +
+ energy -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qqrd2e*scale;
}
diff --git a/src/KSPACE/ewald.h b/src/KSPACE/ewald.h
index 95482cfbfe..cd22f987f7 100644
--- a/src/KSPACE/ewald.h
+++ b/src/KSPACE/ewald.h
@@ -27,13 +27,13 @@ namespace LAMMPS_NS {
class Ewald : public KSpace {
public:
Ewald(class LAMMPS *, int, char **);
- ~Ewald();
+ virtual ~Ewald();
void init();
void setup();
- void compute(int, int);
+ virtual void compute(int, int);
double memory_usage();
- private:
+ protected:
double precision;
int kcount,kmax,kmax3d,kmax_created;
double qqrd2e;
@@ -48,9 +48,9 @@ class Ewald : public KSpace {
double *sfacrl,*sfacim,*sfacrl_all,*sfacim_all;
double ***cs,***sn;
- void eik_dot_r();
+ virtual void eik_dot_r();
void coeffs();
- void allocate();
+ virtual void allocate();
void deallocate();
void slabcorr(int);
};
diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
index bc3cfc98dd..ece15c6839 100644
--- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
+++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
@@ -400,7 +400,8 @@ void PairLJCutCoulLongTIP4P::init_style()
error->all(FLERR,"Pair style lj/cut/coul/long/tip4p requires newton pair on");
if (!atom->q_flag)
error->all(FLERR,"Pair style lj/cut/coul/long/tip4p requires atom attribute q");
- if (strcmp(force->kspace_style,"pppm/tip4p") != 0)
+ if ( (strcmp(force->kspace_style,"pppm/tip4p") != 0) &&
+ (strcmp(force->kspace_style,"pppm/tip4p/proxy") != 0) )
error->all(FLERR,"Pair style is incompatible with KSpace style");
if (force->bond == NULL)
error->all(FLERR,"Must use a bond style with TIP4P potential");
diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp
index 221778e6ee..e2d017a6e5 100644
--- a/src/KSPACE/pppm.cpp
+++ b/src/KSPACE/pppm.cpp
@@ -155,7 +155,8 @@ void PPPM::init()
qdist = 0.0;
- if (strcmp(force->kspace_style,"pppm/tip4p") == 0) {
+ if ( (strcmp(force->kspace_style,"pppm/tip4p") == 0) ||
+ (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0) ) {
if (force->pair == NULL)
error->all(FLERR,"KSpace style is incompatible with Pair style");
double *p_qdist = (double *) force->pair->extract("qdist",itmp);
@@ -184,6 +185,16 @@ void PPPM::init()
alpha = qdist / (cos(0.5*theta) * blen);
}
+ // if we have a /proxy pppm version check if the pair style is compatible
+
+ if ( (strcmp(force->kspace_style,"pppm/proxy") == 0) ||
+ (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0) ) {
+ if (force->pair == NULL)
+ error->all(FLERR,"KSpace style is incompatible with Pair style");
+ if (strstr(force->pair_style,"pppm/") == NULL )
+ error->all(FLERR,"KSpace style is incompatible with Pair style");
+ }
+
// compute qsum & qsqsum and warn if not charge-neutral
qsum = qsqsum = 0.0;
@@ -709,7 +720,7 @@ void PPPM::compute(int eflag, int vflag)
energy = energy_all;
energy *= 0.5*volume;
- energy -= g_ewald*qsqsum/1.772453851 +
+ energy -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qqrd2e*scale;
}
@@ -1515,8 +1526,7 @@ void PPPM::make_rho()
// clear 3d density array
- FFT_SCALAR *vec = &density_brick[nzlo_out][nylo_out][nxlo_out];
- for (i = 0; i < ngrid; i++) vec[i] = ZEROF;
+ memset(&(density_brick[nzlo_out][nylo_out][nxlo_out]),0,ngrid*sizeof(FFT_SCALAR));
// loop over my charges, add their contribution to nearby grid points
// (nx,ny,nz) = global coords of grid pt to "lower left" of charge
@@ -1777,17 +1787,19 @@ void PPPM::compute_rho1d(const FFT_SCALAR &dx, const FFT_SCALAR &dy,
const FFT_SCALAR &dz)
{
int k,l;
+ FFT_SCALAR r1,r2,r3;
for (k = (1-order)/2; k <= order/2; k++) {
- rho1d[0][k] = ZEROF;
- rho1d[1][k] = ZEROF;
- rho1d[2][k] = ZEROF;
+ r1 = r2 = r3 = ZEROF;
for (l = order-1; l >= 0; l--) {
- rho1d[0][k] = rho_coeff[l][k] + rho1d[0][k]*dx;
- rho1d[1][k] = rho_coeff[l][k] + rho1d[1][k]*dy;
- rho1d[2][k] = rho_coeff[l][k] + rho1d[2][k]*dz;
+ r1 = rho_coeff[l][k] + r1*dx;
+ r2 = rho_coeff[l][k] + r2*dy;
+ r3 = rho_coeff[l][k] + r3*dz;
}
+ rho1d[0][k] = r1;
+ rho1d[1][k] = r2;
+ rho1d[2][k] = r3;
}
}
diff --git a/src/KSPACE/pppm_cg.cpp b/src/KSPACE/pppm_cg.cpp
index 5f57715f9f..daffb5ff2d 100644
--- a/src/KSPACE/pppm_cg.cpp
+++ b/src/KSPACE/pppm_cg.cpp
@@ -178,7 +178,7 @@ void PPPMCG::compute(int eflag, int vflag)
energy = energy_all;
energy *= 0.5*volume;
- energy -= g_ewald*qsqsum/1.772453851 +
+ energy -= g_ewald*qsqsum/MY_PIS +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
energy *= qqrd2e*scale;
}
diff --git a/src/MANYBODY/pair_adp.cpp b/src/MANYBODY/pair_adp.cpp
index 0f6b0a6a9d..952c1b73a5 100644
--- a/src/MANYBODY/pair_adp.cpp
+++ b/src/MANYBODY/pair_adp.cpp
@@ -1022,8 +1022,7 @@ void PairADP::unpack_reverse_comm(int n, int *list, double *buf)
double PairADP::memory_usage()
{
- double bytes = maxeatom * sizeof(double);
- bytes += maxvatom*6 * sizeof(double);
+ double bytes = Pair::memory_usage();
bytes += 21 * nmax * sizeof(double);
return bytes;
}
diff --git a/src/MOLECULE/angle_charmm.h b/src/MOLECULE/angle_charmm.h
index 81d23867b6..affd5f1ec3 100644
--- a/src/MOLECULE/angle_charmm.h
+++ b/src/MOLECULE/angle_charmm.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class AngleCharmm : public Angle {
public:
AngleCharmm(class LAMMPS *);
- ~AngleCharmm();
- void compute(int, int);
+ virtual ~AngleCharmm();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_angle(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, int, int, int);
- private:
+ protected:
double *k,*theta0,*k_ub,*r_ub;
void allocate();
diff --git a/src/MOLECULE/angle_cosine.h b/src/MOLECULE/angle_cosine.h
index 1f4aabe456..6b73cf50f1 100644
--- a/src/MOLECULE/angle_cosine.h
+++ b/src/MOLECULE/angle_cosine.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class AngleCosine : public Angle {
public:
AngleCosine(class LAMMPS *);
- ~AngleCosine();
- void compute(int, int);
+ virtual ~AngleCosine();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_angle(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, int, int, int);
- private:
+ protected:
double *k;
void allocate();
diff --git a/src/MOLECULE/angle_cosine_delta.h b/src/MOLECULE/angle_cosine_delta.h
index a6848cc9a5..52399422fc 100644
--- a/src/MOLECULE/angle_cosine_delta.h
+++ b/src/MOLECULE/angle_cosine_delta.h
@@ -28,7 +28,7 @@ namespace LAMMPS_NS {
class AngleCosineDelta : public AngleCosineSquared {
public:
AngleCosineDelta(class LAMMPS *);
- void compute(int, int);
+ virtual void compute(int, int);
double single(int, int, int, int);
};
diff --git a/src/MOLECULE/angle_cosine_periodic.h b/src/MOLECULE/angle_cosine_periodic.h
index 4fb901ebe5..fa0a56ab89 100644
--- a/src/MOLECULE/angle_cosine_periodic.h
+++ b/src/MOLECULE/angle_cosine_periodic.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class AngleCosinePeriodic : public Angle {
public:
AngleCosinePeriodic(class LAMMPS *);
- ~AngleCosinePeriodic();
- void compute(int, int);
+ virtual ~AngleCosinePeriodic();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_angle(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, int, int, int);
- private:
+ protected:
double *k;
int *multiplicity,*b;
diff --git a/src/MOLECULE/angle_harmonic.h b/src/MOLECULE/angle_harmonic.h
index f91a84f7fa..7c13d926fb 100644
--- a/src/MOLECULE/angle_harmonic.h
+++ b/src/MOLECULE/angle_harmonic.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class AngleHarmonic : public Angle {
public:
AngleHarmonic(class LAMMPS *);
- ~AngleHarmonic();
- void compute(int, int);
+ virtual ~AngleHarmonic();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_angle(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, int, int, int);
- private:
+ protected:
double *k,*theta0;
void allocate();
diff --git a/src/MOLECULE/angle_table.h b/src/MOLECULE/angle_table.h
index 23ce4874ca..b404408349 100644
--- a/src/MOLECULE/angle_table.h
+++ b/src/MOLECULE/angle_table.h
@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
class AngleTable : public Angle {
public:
AngleTable(class LAMMPS *);
- ~AngleTable();
- void compute(int, int);
+ virtual ~AngleTable();
+ virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double equilibrium_angle(int);
@@ -37,7 +37,7 @@ class AngleTable : public Angle {
void read_restart(FILE *);
double single(int, int, int, int);
- private:
+ protected:
int tabstyle,tablength;
double *theta0;
diff --git a/src/MOLECULE/bond_fene.h b/src/MOLECULE/bond_fene.h
index 25ca928e31..6fb262c332 100644
--- a/src/MOLECULE/bond_fene.h
+++ b/src/MOLECULE/bond_fene.h
@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
class BondFENE : public Bond {
public:
BondFENE(class LAMMPS *);
- ~BondFENE();
- void compute(int, int);
+ virtual ~BondFENE();
+ virtual void compute(int, int);
void coeff(int, char **);
void init_style();
double equilibrium_distance(int);
@@ -37,7 +37,7 @@ class BondFENE : public Bond {
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
double TWO_1_3;
double *k,*r0,*epsilon,*sigma;
diff --git a/src/MOLECULE/bond_fene_expand.h b/src/MOLECULE/bond_fene_expand.h
index a86dfdba14..89f28c4bf3 100644
--- a/src/MOLECULE/bond_fene_expand.h
+++ b/src/MOLECULE/bond_fene_expand.h
@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
class BondFENEExpand : public Bond {
public:
BondFENEExpand(class LAMMPS *);
- ~BondFENEExpand();
- void compute(int, int);
+ virtual ~BondFENEExpand();
+ virtual void compute(int, int);
void coeff(int, char **);
void init_style();
double equilibrium_distance(int);
@@ -37,7 +37,7 @@ class BondFENEExpand : public Bond {
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
double TWO_1_3;
double *k,*r0,*epsilon,*sigma,*shift;
diff --git a/src/MOLECULE/bond_harmonic.h b/src/MOLECULE/bond_harmonic.h
index a4364cc8bb..87eefd6a0d 100644
--- a/src/MOLECULE/bond_harmonic.h
+++ b/src/MOLECULE/bond_harmonic.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class BondHarmonic : public Bond {
public:
BondHarmonic(class LAMMPS *);
- ~BondHarmonic();
- void compute(int, int);
+ virtual ~BondHarmonic();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
double *k,*r0;
void allocate();
diff --git a/src/MOLECULE/bond_morse.h b/src/MOLECULE/bond_morse.h
index 34fbf265cc..1dcf3653e1 100644
--- a/src/MOLECULE/bond_morse.h
+++ b/src/MOLECULE/bond_morse.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class BondMorse : public Bond {
public:
BondMorse(class LAMMPS *);
- ~BondMorse();
- void compute(int, int);
+ virtual ~BondMorse();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
double *d0,*alpha,*r0;
void allocate();
diff --git a/src/MOLECULE/bond_nonlinear.h b/src/MOLECULE/bond_nonlinear.h
index 1390d25193..8d715e2e61 100644
--- a/src/MOLECULE/bond_nonlinear.h
+++ b/src/MOLECULE/bond_nonlinear.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class BondNonlinear : public Bond {
public:
BondNonlinear(class LAMMPS *);
- ~BondNonlinear();
- void compute(int, int);
+ virtual ~BondNonlinear();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
double *epsilon,*r0,*lamda;
void allocate();
diff --git a/src/MOLECULE/bond_quartic.h b/src/MOLECULE/bond_quartic.h
index 67ae294dec..628d2c237b 100644
--- a/src/MOLECULE/bond_quartic.h
+++ b/src/MOLECULE/bond_quartic.h
@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
class BondQuartic : public Bond {
public:
BondQuartic(class LAMMPS *);
- ~BondQuartic();
- void compute(int, int);
+ virtual ~BondQuartic();
+ virtual void compute(int, int);
void coeff(int, char **);
void init_style();
double equilibrium_distance(int);
@@ -37,7 +37,7 @@ class BondQuartic : public Bond {
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
double TWO_1_3;
double *k,*b1,*b2,*rc,*u0;
diff --git a/src/MOLECULE/bond_table.h b/src/MOLECULE/bond_table.h
index f0d289f055..810693262f 100644
--- a/src/MOLECULE/bond_table.h
+++ b/src/MOLECULE/bond_table.h
@@ -28,8 +28,8 @@ namespace LAMMPS_NS {
class BondTable : public Bond {
public:
BondTable(class LAMMPS *);
- ~BondTable();
- void compute(int, int);
+ virtual ~BondTable();
+ virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
double equilibrium_distance(int);
@@ -37,7 +37,7 @@ class BondTable : public Bond {
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
int tabstyle,tablength;
double *r0;
diff --git a/src/MOLECULE/improper_cvff.h b/src/MOLECULE/improper_cvff.h
index ccef917d3a..86de363a15 100644
--- a/src/MOLECULE/improper_cvff.h
+++ b/src/MOLECULE/improper_cvff.h
@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
class ImproperCvff : public Improper {
public:
ImproperCvff(class LAMMPS *);
- ~ImproperCvff();
- void compute(int, int);
+ virtual ~ImproperCvff();
+ virtual void compute(int, int);
void coeff(int, char **);
void write_restart(FILE *);
void read_restart(FILE *);
- private:
+ protected:
double *k;
int *sign,*multiplicity;
diff --git a/src/MOLECULE/improper_harmonic.h b/src/MOLECULE/improper_harmonic.h
index a14e222d60..b3d5ef81f6 100644
--- a/src/MOLECULE/improper_harmonic.h
+++ b/src/MOLECULE/improper_harmonic.h
@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
class ImproperHarmonic : public Improper {
public:
ImproperHarmonic(class LAMMPS *);
- ~ImproperHarmonic();
- void compute(int, int);
+ virtual ~ImproperHarmonic();
+ virtual void compute(int, int);
void coeff(int, char **);
void write_restart(FILE *);
void read_restart(FILE *);
- private:
+ protected:
double *k,*chi;
void allocate();
diff --git a/src/MOLECULE/improper_umbrella.h b/src/MOLECULE/improper_umbrella.h
index cef7583743..483979826e 100644
--- a/src/MOLECULE/improper_umbrella.h
+++ b/src/MOLECULE/improper_umbrella.h
@@ -28,13 +28,13 @@ namespace LAMMPS_NS {
class ImproperUmbrella : public Improper {
public:
ImproperUmbrella(class LAMMPS *);
- ~ImproperUmbrella();
- void compute(int, int);
+ virtual ~ImproperUmbrella();
+ virtual void compute(int, int);
void coeff(int, char **);
void write_restart(FILE *);
void read_restart(FILE *);
- private:
+ protected:
double *kw, *w0, *C;
void allocate();
diff --git a/src/OPT/pair_lj_cut_coul_long_tip4p_opt.cpp b/src/OPT/pair_lj_cut_coul_long_tip4p_opt.cpp
index 24db63e579..aee802db21 100644
--- a/src/OPT/pair_lj_cut_coul_long_tip4p_opt.cpp
+++ b/src/OPT/pair_lj_cut_coul_long_tip4p_opt.cpp
@@ -162,7 +162,7 @@ void PairLJCutCoulLongTIP4POpt::eval()
double r,r2inv,r6inv,forcecoul,forcelj,cforce;
double factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc,ddotf;
- double xiM[3],xjM[3],v[6],xH1[3],xH2[3];
+ double v[6],xH1[3],xH2[3];
double fdx,fdy,fdz,f1x,f1y,f1z,fOx,fOy,fOz,fHx,fHy,fHz;
double *x1,*x2;
int *ilist,*jlist,*numneigh,**firstneigh;
diff --git a/src/SHOCK/fix_wall_piston.cpp b/src/SHOCK/fix_wall_piston.cpp
index 2d610ea2e0..3b66583715 100644
--- a/src/SHOCK/fix_wall_piston.cpp
+++ b/src/SHOCK/fix_wall_piston.cpp
@@ -24,8 +24,10 @@
#include "random_mars.h"
#include "force.h"
#include "comm.h"
+#include "math_const.h"
using namespace LAMMPS_NS;
+using namespace MathConst;
/* ---------------------------------------------------------------------- */
@@ -203,7 +205,7 @@ void FixWallPiston::post_integrate()
}
else if (rampNL1flag) {
paccelz = maxvz / tott;
- angfreq = 6.283185 / (0.5 * tott);
+ angfreq = MY_2PI / (0.5 * tott);
if (zloflag) {
zlo = z0 + paccelz * (0.5*tt + 1.0/(angfreq*angfreq) - 1.0/(angfreq*angfreq)*cos(angfreq*t));
@@ -213,7 +215,7 @@ void FixWallPiston::post_integrate()
}
else if (rampNL2flag) {
paccelz = maxvz / tott;
- angfreq = 18.84956 / tott;
+ angfreq = 3.0*MY_2PI / tott;
if (zloflag) {
zlo = z0 + paccelz * (0.5*tt + 4.0/(3.0*angfreq*angfreq)*(1.0-cos(angfreq*t)) + 1.0/(6.0*angfreq*angfreq)*(1.0-cos(2.0*angfreq*t)));
diff --git a/src/USER-MISC/angle_cosine_shift.h b/src/USER-MISC/angle_cosine_shift.h
index ce44ea3045..d8b7ea1741 100644
--- a/src/USER-MISC/angle_cosine_shift.h
+++ b/src/USER-MISC/angle_cosine_shift.h
@@ -28,7 +28,7 @@ namespace LAMMPS_NS {
class AngleCosineShift : public Angle {
public:
AngleCosineShift(class LAMMPS *);
- ~AngleCosineShift();
+ virtual ~AngleCosineShift();
virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_angle(int);
@@ -36,7 +36,7 @@ class AngleCosineShift : public Angle {
void read_restart(FILE *);
double single(int, int, int, int);
- private:
+ protected:
double *k;
double *a;
double *theta;
diff --git a/src/USER-MISC/angle_cosine_shift_exp.cpp b/src/USER-MISC/angle_cosine_shift_exp.cpp
index cd6a4960d9..031cfec563 100644
--- a/src/USER-MISC/angle_cosine_shift_exp.cpp
+++ b/src/USER-MISC/angle_cosine_shift_exp.cpp
@@ -57,7 +57,7 @@ void AngleCosineShiftExp::compute(int eflag, int vflag)
int i1,i2,i3,n,type;
double delx1,dely1,delz1,delx2,dely2,delz2;
double eangle,f1[3],f3[3],ff;
- double rsq1,rsq2,r1,r2,c,s,cc,ss,a11,a12,a22;
+ double rsq1,rsq2,r1,r2,c,s,a11,a12,a22;
double exp2,aa,uumin,cccpsss,cssmscc;
eangle = 0.0;
diff --git a/src/USER-MISC/angle_cosine_shift_exp.h b/src/USER-MISC/angle_cosine_shift_exp.h
index c51ece7c44..50883dc034 100644
--- a/src/USER-MISC/angle_cosine_shift_exp.h
+++ b/src/USER-MISC/angle_cosine_shift_exp.h
@@ -26,15 +26,15 @@ namespace LAMMPS_NS {
class AngleCosineShiftExp : public Angle {
public:
AngleCosineShiftExp(class LAMMPS *);
- ~AngleCosineShiftExp();
- void compute(int, int);
+ virtual ~AngleCosineShiftExp();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_angle(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, int, int, int);
- private:
+ protected:
bool *doExpansion;
double *umin,*a,*opt1;
double *theta0;
diff --git a/src/USER-MISC/bond_harmonic_shift.h b/src/USER-MISC/bond_harmonic_shift.h
index 91bb723323..b3378aea3f 100644
--- a/src/USER-MISC/bond_harmonic_shift.h
+++ b/src/USER-MISC/bond_harmonic_shift.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class BondHarmonicShift : public Bond {
public:
BondHarmonicShift(class LAMMPS *);
- ~BondHarmonicShift();
- void compute(int, int);
+ virtual ~BondHarmonicShift();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
double *k,*r0,*r1;
void allocate();
diff --git a/src/USER-MISC/bond_harmonic_shift_cut.h b/src/USER-MISC/bond_harmonic_shift_cut.h
index 7bb63ca64b..6e957d286c 100644
--- a/src/USER-MISC/bond_harmonic_shift_cut.h
+++ b/src/USER-MISC/bond_harmonic_shift_cut.h
@@ -28,15 +28,15 @@ namespace LAMMPS_NS {
class BondHarmonicShiftCut : public Bond {
public:
BondHarmonicShiftCut(class LAMMPS *);
- ~BondHarmonicShiftCut();
- void compute(int, int);
+ virtual ~BondHarmonicShiftCut();
+ virtual void compute(int, int);
void coeff(int, char **);
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
double single(int, double, int, int);
- private:
+ protected:
double *k,*r0,*r1;
void allocate();
diff --git a/src/angle.cpp b/src/angle.cpp
index 58ba679948..54a3a5d051 100644
--- a/src/angle.cpp
+++ b/src/angle.cpp
@@ -14,6 +14,7 @@
#include "math.h"
#include "angle.h"
#include "atom.h"
+#include "comm.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
@@ -78,12 +79,12 @@ void Angle::ev_setup(int eflag, int vflag)
if (eflag_atom && atom->nmax > maxeatom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
- memory->create(eatom,maxeatom,"bond:eatom");
+ memory->create(eatom,comm->nthreads*maxeatom,"bond:eatom");
}
if (vflag_atom && atom->nmax > maxvatom) {
maxvatom = atom->nmax;
memory->destroy(vatom);
- memory->create(vatom,maxvatom,6,"bond:vatom");
+ memory->create(vatom,comm->nthreads*maxvatom,6,"bond:vatom");
}
// zero accumulators
@@ -216,7 +217,7 @@ void Angle::ev_tally(int i, int j, int k, int nlocal, int newton_bond,
double Angle::memory_usage()
{
- double bytes = maxeatom * sizeof(double);
- bytes += maxvatom*6 * sizeof(double);
+ double bytes = comm->nthreads*maxeatom * sizeof(double);
+ bytes += comm->nthreads*maxvatom*6 * sizeof(double);
return bytes;
}
diff --git a/src/angle.h b/src/angle.h
index 7f25f38289..46ee5e1dc6 100644
--- a/src/angle.h
+++ b/src/angle.h
@@ -20,6 +20,7 @@
namespace LAMMPS_NS {
class Angle : protected Pointers {
+ friend class ThrOMP;
public:
int allocated;
int *setflag;
diff --git a/src/atom_vec_line.cpp b/src/atom_vec_line.cpp
index ebf8b7d7a8..c49cd6b6ad 100644
--- a/src/atom_vec_line.cpp
+++ b/src/atom_vec_line.cpp
@@ -17,6 +17,7 @@
#include "string.h"
#include "atom_vec_line.h"
#include "atom.h"
+#include "comm.h"
#include "domain.h"
#include "modify.h"
#include "force.h"
@@ -83,6 +84,8 @@ void AtomVecLine::grow(int n)
if (n == 0) nmax += DELTA;
else nmax = n;
atom->nmax = nmax;
+ if (nmax < 0 || nmax > MAXSMALLINT)
+ error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");
type = memory->grow(atom->type,nmax,"atom:type");
@@ -90,12 +93,12 @@ void AtomVecLine::grow(int n)
image = memory->grow(atom->image,nmax,"atom:image");
x = memory->grow(atom->x,nmax,3,"atom:x");
v = memory->grow(atom->v,nmax,3,"atom:v");
- f = memory->grow(atom->f,nmax,3,"atom:f");
+ f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
molecule = memory->grow(atom->molecule,nmax,"atom:molecule");
rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
omega = memory->grow(atom->omega,nmax,3,"atom:omega");
- torque = memory->grow(atom->torque,nmax,3,"atom:torque");
+ torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque");
line = memory->grow(atom->line,nmax,"atom:line");
if (atom->nextra_grow)
@@ -1115,12 +1118,12 @@ bigint AtomVecLine::memory_usage()
if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
- if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
+ if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
if (atom->memcheck("omega")) bytes += memory->usage(omega,nmax,3);
- if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3);
+ if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax*comm->nthreads,3);
if (atom->memcheck("line")) bytes += memory->usage(line,nmax);
bytes += nmax_bonus*sizeof(Bonus);
diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp
index 03aed7aa65..0f3fbe161a 100644
--- a/src/atom_vec_tri.cpp
+++ b/src/atom_vec_tri.cpp
@@ -18,6 +18,7 @@
#include "atom_vec_tri.h"
#include "math_extra.h"
#include "atom.h"
+#include "comm.h"
#include "domain.h"
#include "modify.h"
#include "force.h"
@@ -84,6 +85,8 @@ void AtomVecTri::grow(int n)
if (n == 0) nmax += DELTA;
else nmax = n;
atom->nmax = nmax;
+ if (nmax < 0 || nmax > MAXSMALLINT)
+ error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");
type = memory->grow(atom->type,nmax,"atom:type");
@@ -91,12 +94,12 @@ void AtomVecTri::grow(int n)
image = memory->grow(atom->image,nmax,"atom:image");
x = memory->grow(atom->x,nmax,3,"atom:x");
v = memory->grow(atom->v,nmax,3,"atom:v");
- f = memory->grow(atom->f,nmax,3,"atom:f");
+ f = memory->grow(atom->f,nmax*comm->nthreads,3,"atom:f");
molecule = memory->grow(atom->molecule,nmax,"atom:molecule");
rmass = memory->grow(atom->rmass,nmax,"atom:rmass");
angmom = memory->grow(atom->angmom,nmax,3,"atom:angmom");
- torque = memory->grow(atom->torque,nmax,3,"atom:torque");
+ torque = memory->grow(atom->torque,nmax*comm->nthreads,3,"atom:torque");
tri = memory->grow(atom->tri,nmax,"atom:tri");
if (atom->nextra_grow)
@@ -1547,12 +1550,12 @@ bigint AtomVecTri::memory_usage()
if (atom->memcheck("image")) bytes += memory->usage(image,nmax);
if (atom->memcheck("x")) bytes += memory->usage(x,nmax,3);
if (atom->memcheck("v")) bytes += memory->usage(v,nmax,3);
- if (atom->memcheck("f")) bytes += memory->usage(f,nmax,3);
+ if (atom->memcheck("f")) bytes += memory->usage(f,nmax*comm->nthreads,3);
if (atom->memcheck("molecule")) bytes += memory->usage(molecule,nmax);
if (atom->memcheck("rmass")) bytes += memory->usage(rmass,nmax);
if (atom->memcheck("angmom")) bytes += memory->usage(angmom,nmax,3);
- if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax,3);
+ if (atom->memcheck("torque")) bytes += memory->usage(torque,nmax*comm->nthreads,3);
if (atom->memcheck("tri")) bytes += memory->usage(tri,nmax);
bytes += nmax_bonus*sizeof(Bonus);
diff --git a/src/bond.cpp b/src/bond.cpp
index 92310df653..aa168ff45c 100644
--- a/src/bond.cpp
+++ b/src/bond.cpp
@@ -14,6 +14,7 @@
#include "string.h"
#include "bond.h"
#include "atom.h"
+#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
@@ -80,12 +81,12 @@ void Bond::ev_setup(int eflag, int vflag)
if (eflag_atom && atom->nmax > maxeatom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
- memory->create(eatom,maxeatom,"bond:eatom");
+ memory->create(eatom,comm->nthreads*maxeatom,"bond:eatom");
}
if (vflag_atom && atom->nmax > maxvatom) {
maxvatom = atom->nmax;
memory->destroy(vatom);
- memory->create(vatom,maxvatom,6,"bond:vatom");
+ memory->create(vatom,comm->nthreads*maxvatom,6,"bond:vatom");
}
// zero accumulators
@@ -198,7 +199,7 @@ void Bond::ev_tally(int i, int j, int nlocal, int newton_bond,
double Bond::memory_usage()
{
- double bytes = maxeatom * sizeof(double);
- bytes += maxvatom*6 * sizeof(double);
+ double bytes = comm->nthreads*maxeatom * sizeof(double);
+ bytes += comm->nthreads*maxvatom*6 * sizeof(double);
return bytes;
}
diff --git a/src/bond.h b/src/bond.h
index 551d84b51f..29f9862fb7 100644
--- a/src/bond.h
+++ b/src/bond.h
@@ -20,6 +20,7 @@
namespace LAMMPS_NS {
class Bond : protected Pointers {
+ friend class ThrOMP;
public:
int allocated;
int *setflag;
diff --git a/src/dihedral.cpp b/src/dihedral.cpp
index 974df81b14..f1a62cba16 100644
--- a/src/dihedral.cpp
+++ b/src/dihedral.cpp
@@ -14,7 +14,7 @@
#include "math.h"
#include "dihedral.h"
#include "atom.h"
-#include "atom.h"
+#include "comm.h"
#include "force.h"
#include "pair.h"
#include "memory.h"
@@ -82,12 +82,12 @@ void Dihedral::ev_setup(int eflag, int vflag)
if (eflag_atom && atom->nmax > maxeatom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
- memory->create(eatom,maxeatom,"bond:eatom");
+ memory->create(eatom,comm->nthreads*maxeatom,"bond:eatom");
}
if (vflag_atom && atom->nmax > maxvatom) {
maxvatom = atom->nmax;
memory->destroy(vatom);
- memory->create(vatom,maxvatom,6,"bond:vatom");
+ memory->create(vatom,comm->nthreads*maxvatom,6,"bond:vatom");
}
// zero accumulators
@@ -242,7 +242,7 @@ void Dihedral::ev_tally(int i1, int i2, int i3, int i4,
double Dihedral::memory_usage()
{
- double bytes = maxeatom * sizeof(double);
- bytes += maxvatom*6 * sizeof(double);
+ double bytes = comm->nthreads*maxeatom * sizeof(double);
+ bytes += comm->nthreads*maxvatom*6 * sizeof(double);
return bytes;
}
diff --git a/src/domain.cpp b/src/domain.cpp
index a1a0045fa9..3140ea6a54 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -603,7 +603,7 @@ void Domain::minimum_image(double *delta)
for triclinic, also add/subtract tilt factors in other dims as needed
------------------------------------------------------------------------- */
-void Domain::closest_image(double *xi, double *xj, double *xjimage)
+void Domain::closest_image(const double * const xi, const double * const xj, double * const xjimage)
{
double dx,dy,dz;
diff --git a/src/domain.h b/src/domain.h
index c07808ca02..1e04ab4ab9 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -93,16 +93,15 @@ class Domain : protected Pointers {
virtual void set_local_box();
virtual void reset_box();
virtual void pbc();
- virtual void remap(double *, int &);
- virtual void remap(double *);
- virtual void remap_near(double *, double *);
- virtual void unmap(double *, int);
- virtual void unmap(double *, int, double *);
- virtual int minimum_image_check(double, double, double);
- virtual void minimum_image(double &, double &, double &);
- virtual void minimum_image(double *);
- virtual void closest_image(double *, double *, double *);
-
+ void remap(double *, int &);
+ void remap(double *);
+ void remap_near(double *, double *);
+ void unmap(double *, int);
+ void unmap(double *, int, double *);
+ int minimum_image_check(double, double, double);
+ void minimum_image(double &, double &, double &);
+ void minimum_image(double *);
+ void closest_image(const double * const, const double * const, double * const);
void set_lattice(int, char **);
void add_region(int, char **);
void delete_region(int, char **);
diff --git a/src/force.cpp b/src/force.cpp
index 67b3873309..5a26d16aa0 100644
--- a/src/force.cpp
+++ b/src/force.cpp
@@ -478,24 +478,61 @@ Improper *Force::new_improper(const char *style, const char *suffix, int &sflag)
new kspace style
------------------------------------------------------------------------- */
-void Force::create_kspace(int narg, char **arg)
+void Force::create_kspace(int narg, char **arg, const char *suffix)
{
delete [] kspace_style;
if (kspace) delete kspace;
- if (strcmp(arg[0],"none") == 0) kspace = NULL;
+ int sflag;
+ kspace = new_kspace(narg,arg,suffix,sflag);
+
+ if (sflag) {
+ char estyle[256];
+ sprintf(estyle,"%s/%s",arg[0],suffix);
+ int n = strlen(estyle) + 1;
+ kspace_style = new char[n];
+ strcpy(kspace_style,estyle);
+ } else {
+ int n = strlen(arg[0]) + 1;
+ kspace_style = new char[n];
+ strcpy(kspace_style,arg[0]);
+ }
+}
+
+/* ----------------------------------------------------------------------
+ generate a kspace class
+------------------------------------------------------------------------- */
+
+KSpace *Force::new_kspace(int narg, char **arg, const char *suffix, int &sflag)
+{
+ if (suffix && lmp->suffix_enable) {
+ sflag = 1;
+ char estyle[256];
+ sprintf(estyle,"%s/%s",arg[0],suffix);
+
+ if (0) return NULL;
#define KSPACE_CLASS
#define KSpaceStyle(key,Class) \
- else if (strcmp(arg[0],#key) == 0) kspace = new Class(lmp,narg-1,&arg[1]);
+ else if (strcmp(estyle,#key) == 0) return new Class(lmp,narg-1,&arg[1]);
#include "style_kspace.h"
+#undef KSpaceStyle
#undef KSPACE_CLASS
- else error->all(FLERR,"Invalid kspace style");
+ }
- int n = strlen(arg[0]) + 1;
- kspace_style = new char[n];
- strcpy(kspace_style,arg[0]);
+ sflag = 0;
+
+ if (strcmp(arg[0],"none") == 0) return NULL;
+
+#define KSPACE_CLASS
+#define KSpaceStyle(key,Class) \
+ else if (strcmp(arg[0],#key) == 0) return new Class(lmp,narg-1,&arg[1]);
+#include "style_kspace.h"
+#undef KSPACE_CLASS
+
+ else error->all(FLERR,"Invalid kspace style");
+ return NULL;
}
/* ----------------------------------------------------------------------
diff --git a/src/force.h b/src/force.h
index 654b81a3a8..a7089070d2 100644
--- a/src/force.h
+++ b/src/force.h
@@ -86,7 +86,8 @@ class Force : protected Pointers {
void create_improper(const char *, const char *suffix = NULL);
class Improper *new_improper(const char *, const char *, int &);
- void create_kspace(int, char **);
+ void create_kspace(int, char **, const char *suffix = NULL);
+ class KSpace *new_kspace(int, char **, const char *, int &);
void set_special(int, char **);
void bounds(char *, int, int &, int &);
diff --git a/src/improper.cpp b/src/improper.cpp
index 8f5ab4d6e3..9a1f47e5ce 100644
--- a/src/improper.cpp
+++ b/src/improper.cpp
@@ -14,6 +14,7 @@
#include "math.h"
#include "improper.h"
#include "atom.h"
+#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
@@ -76,12 +77,12 @@ void Improper::ev_setup(int eflag, int vflag)
if (eflag_atom && atom->nmax > maxeatom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
- memory->create(eatom,maxeatom,"bond:eatom");
+ memory->create(eatom,comm->nthreads*maxeatom,"bond:eatom");
}
if (vflag_atom && atom->nmax > maxvatom) {
maxvatom = atom->nmax;
memory->destroy(vatom);
- memory->create(vatom,maxvatom,6,"bond:vatom");
+ memory->create(vatom,comm->nthreads*maxvatom,6,"bond:vatom");
}
// zero accumulators
@@ -236,7 +237,7 @@ void Improper::ev_tally(int i1, int i2, int i3, int i4,
double Improper::memory_usage()
{
- double bytes = maxeatom * sizeof(double);
- bytes += maxvatom*6 * sizeof(double);
+ double bytes = comm->nthreads*maxeatom * sizeof(double);
+ bytes += comm->nthreads*maxvatom*6 * sizeof(double);
return bytes;
}
diff --git a/src/improper.h b/src/improper.h
index d8310b2475..d766015e34 100644
--- a/src/improper.h
+++ b/src/improper.h
@@ -20,6 +20,7 @@
namespace LAMMPS_NS {
class Improper : protected Pointers {
+ friend class ThrOMP;
public:
int allocated;
int *setflag;
diff --git a/src/input.cpp b/src/input.cpp
index 2251335b43..457faf4ac8 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -1022,7 +1022,8 @@ void Input::improper_style()
void Input::kspace_modify()
{
- if (force->kspace == NULL) error->all(FLERR,"KSpace style has not yet been set");
+ if (force->kspace == NULL)
+ error->all(FLERR,"KSpace style has not yet been set");
force->kspace->modify_params(narg,arg);
}
@@ -1030,7 +1031,7 @@ void Input::kspace_modify()
void Input::kspace_style()
{
- force->create_kspace(narg,arg);
+ force->create_kspace(narg,arg,lmp->suffix);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/integrate.cpp b/src/integrate.cpp
index 80bf89f676..6de4df26c3 100644
--- a/src/integrate.cpp
+++ b/src/integrate.cpp
@@ -26,6 +26,7 @@ Integrate::Integrate(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
{
elist_global = elist_atom = NULL;
vlist_global = vlist_atom = NULL;
+ external_force_clear = 0;
}
/* ---------------------------------------------------------------------- */
diff --git a/src/integrate.h b/src/integrate.h
index 8f8942c80c..3c75176a10 100644
--- a/src/integrate.h
+++ b/src/integrate.h
@@ -33,6 +33,7 @@ class Integrate : protected Pointers {
protected:
int eflag,vflag; // flags for energy/virial computation
int virial_style; // compute virial explicitly or implicitly
+ int external_force_clear; // clear forces locally or externally
int nelist_global,nelist_atom; // # of PE,virial computes to check
int nvlist_global,nvlist_atom;
diff --git a/src/kspace.h b/src/kspace.h
index 3ad2193f11..f77bacfcdc 100644
--- a/src/kspace.h
+++ b/src/kspace.h
@@ -19,6 +19,7 @@
namespace LAMMPS_NS {
class KSpace : protected Pointers {
+ friend class ThrOMP;
public:
double energy;
double virial[6];
diff --git a/src/math_const.h b/src/math_const.h
index cf7e8a6112..5cbc2aee28 100644
--- a/src/math_const.h
+++ b/src/math_const.h
@@ -19,6 +19,8 @@ namespace LAMMPS_NS {
namespace MathConst {
static const double THIRD = 1.0/3.0;
static const double MY_PI = 3.14159265358979323846; // pi
+ static const double MY_2PI = 6.28318530717958647692; // 2pi
+ static const double MY_3PI = 9.42477796076937971538; // 3pi
static const double MY_PI2 = 1.57079632679489661923; // pi/2
static const double MY_PI4 = 0.78539816339744830962; // pi/4
static const double MY_PIS = 1.77245385090551602729; // sqrt(pi)
diff --git a/src/min.cpp b/src/min.cpp
index bf83065261..0df114edaf 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -157,7 +157,8 @@ void Min::init()
if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
if (comm->me == 0)
- error->warning(FLERR,"Resetting reneighboring criteria during minimization");
+ error->warning(FLERR,
+ "Resetting reneighboring criteria during minimization");
}
neighbor->every = 1;
@@ -227,9 +228,11 @@ void Min::setup()
// remove these restriction eventually
if (nextra_global && searchflag == 0)
- error->all(FLERR,"Cannot use a damped dynamics min style with fix box/relax");
+ error->all(FLERR,
+ "Cannot use a damped dynamics min style with fix box/relax");
if (nextra_atom && searchflag == 0)
- error->all(FLERR,"Cannot use a damped dynamics min style with per-atom DOF");
+ error->all(FLERR,
+ "Cannot use a damped dynamics min style with per-atom DOF");
// atoms may have migrated in comm->exchange()
@@ -510,6 +513,8 @@ void Min::force_clear()
{
int i;
+ if (external_force_clear) return;
+
// clear global force array
// nall includes ghosts only if either newton flag is set
diff --git a/src/min.h b/src/min.h
index f367fa7b93..47956059d8 100644
--- a/src/min.h
+++ b/src/min.h
@@ -49,6 +49,7 @@ class Min : protected Pointers {
protected:
int eflag,vflag; // flags for energy/virial computation
int virial_style; // compute virial explicitly or implicitly
+ int external_force_clear; // clear forces locally or externally
double dmax; // max dist to move any atom in one step
int linestyle; // 0 = backtrack, 1 = quadratic
diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp
index b4b7fd0aeb..f35acdb4e3 100644
--- a/src/neigh_full.cpp
+++ b/src/neigh_full.cpp
@@ -108,7 +108,7 @@ void Neighbor::full_nsq(NeighList *list)
/* ----------------------------------------------------------------------
N^2 search for all neighbors
- include neighbors of ghost atoms
+ include neighbors of ghost atoms, but no "special neighbors" for ghosts
every neighbor pair appears in list of both atoms i and j
------------------------------------------------------------------------- */
@@ -187,12 +187,9 @@ void Neighbor::full_nsq_ghost(NeighList *list)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
- if (rsq <= cutneighghostsq[itype][jtype]) {
- if (molecular) {
- which = find_special(special[i],nspecial[i],tag[j]);
- if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
- } else neighptr[n++] = j;
- }
+
+ if (rsq <= cutneighghostsq[itype][jtype])
+ neighptr[n++] = j;
}
}
@@ -201,7 +198,8 @@ void Neighbor::full_nsq_ghost(NeighList *list)
numneigh[i] = n;
npnt += n;
if (n > oneatom || npnt >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ error->one(FLERR,
+ "Neighbor list overflow, boost neigh_modify one or page");
}
list->inum = atom->nlocal;
@@ -295,7 +293,8 @@ void Neighbor::full_bin(NeighList *list)
numneigh[i] = n;
npnt += n;
if (n > oneatom || npnt >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ error->one(FLERR,
+ "Neighbor list overflow, boost neigh_modify one or page");
}
list->inum = inum;
@@ -304,7 +303,7 @@ void Neighbor::full_bin(NeighList *list)
/* ----------------------------------------------------------------------
binned neighbor list construction for all neighbors
- include neighbors of ghost atoms
+ include neighbors of ghost atoms, but no "special neighbors" for ghosts
every neighbor pair appears in list of both atoms i and j
------------------------------------------------------------------------- */
@@ -407,13 +406,9 @@ void Neighbor::full_bin_ghost(NeighList *list)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
-
- if (rsq <= cutneighghostsq[itype][jtype]) {
- if (molecular) {
- which = find_special(special[i],nspecial[i],tag[j]);
- if (which >= 0) neighptr[n++] = j ^ (which << SBBITS);
- } else neighptr[n++] = j;
- }
+
+ if (rsq <= cutneighghostsq[itype][jtype])
+ neighptr[n++] = j;
}
}
}
@@ -423,7 +418,8 @@ void Neighbor::full_bin_ghost(NeighList *list)
numneigh[i] = n;
npnt += n;
if (n > oneatom || npnt >= pgsize)
- error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page");
+ error->one(FLERR,
+ "Neighbor list overflow, boost neigh_modify one or page");
}
list->inum = atom->nlocal;
diff --git a/src/neigh_list.cpp b/src/neigh_list.cpp
index dbd8f0d4a0..738bd5d7e1 100644
--- a/src/neigh_list.cpp
+++ b/src/neigh_list.cpp
@@ -15,17 +15,12 @@
#include "neigh_list.h"
#include "atom.h"
#include "comm.h"
+#include "update.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
-
-
-#include "update.h"
-
-
-
using namespace LAMMPS_NS;
#define PGDELTA 1
@@ -255,6 +250,8 @@ void NeighList::print_attributes()
printf(" %d = occasional\n",rq->occasional);
printf(" %d = dnum\n",rq->dnum);
printf(" %d = ghost\n",rq->ghost);
+ printf(" %d = cudable\n",rq->cudable);
+ printf(" %d = omp\n",rq->omp);
printf(" %d = copy\n",rq->copy);
printf(" %d = skip\n",rq->skip);
printf(" %d = otherlist\n",rq->otherlist);
diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp
index 82c9f0b267..2146009a99 100644
--- a/src/neigh_request.cpp
+++ b/src/neigh_request.cpp
@@ -44,6 +44,8 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
// default is encode special bond flags
// default is no auxiliary floating point values
// default is no neighbors of ghosts
+ // default is no CUDA neighbor list build
+ // default is no multi-threaded neighbor list build
occasional = 0;
newton = 0;
@@ -51,6 +53,7 @@ NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
dnum = 0;
ghost = 0;
cudable = 0;
+ omp = 0;
// default is no copy or skip
@@ -101,6 +104,7 @@ int NeighRequest::identical(NeighRequest *other)
if (dnum != other->dnum) same = 0;
if (ghost != other->ghost) same = 0;
if (cudable != other->cudable) same = 0;
+ if (omp != other->omp) same = 0;
if (copy != other->copy) same = 0;
if (same_skip(other) == 0) same = 0;
@@ -129,6 +133,7 @@ int NeighRequest::same_kind(NeighRequest *other)
if (newton != other->newton) same = 0;
if (ghost != other->ghost) same = 0;
if (cudable != other->cudable) same = 0;
+ if (omp != other->omp) same = 0;
return same;
}
@@ -178,4 +183,5 @@ void NeighRequest::copy_request(NeighRequest *other)
dnum = other->dnum;
ghost = other->ghost;
cudable = other->cudable;
+ omp = other->omp;
}
diff --git a/src/neigh_request.h b/src/neigh_request.h
index afff5966c0..e3dcc915a7 100644
--- a/src/neigh_request.h
+++ b/src/neigh_request.h
@@ -76,6 +76,10 @@ class NeighRequest : protected Pointers {
int cudable;
+ // 1 if using multi-threaded neighbor list build
+
+ int omp;
+
// set by neighbor and pair_hybrid after all requests are made
// these settings do not change kind value
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index 58f8ad94dd..24e4285b50 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -1698,9 +1698,8 @@ int Neighbor::coord2bin(double *x, int &ix, int &iy, int &iz)
return 1 if should be excluded, 0 if included
------------------------------------------------------------------------- */
-int Neighbor::exclusion(int i, int j, int itype, int jtype,
- int *mask, int *molecule)
-{
+int Neighbor::exclusion(int i, int j, int itype, int jtype,
+ int *mask, int *molecule) const {
int m;
if (nex_type && ex_type[itype][jtype]) return 1;
diff --git a/src/neighbor.h b/src/neighbor.h
index c2fa635d63..5248c0bd6e 100644
--- a/src/neighbor.h
+++ b/src/neighbor.h
@@ -164,7 +164,7 @@ class Neighbor : protected Pointers {
int coord2bin(double *); // mapping atom coord to a bin
int coord2bin(double *, int &, int &, int&); // ditto
- int exclusion(int, int, int, int, int *, int *); // test for pair exclusion
+ int exclusion(int, int, int, int, int *, int *) const; // test for pair exclusion
virtual void choose_build(int, class NeighRequest *);
void choose_stencil(int, class NeighRequest *);
diff --git a/src/pair.cpp b/src/pair.cpp
index c839078310..c921fdebdf 100644
--- a/src/pair.cpp
+++ b/src/pair.cpp
@@ -281,12 +281,12 @@ void Pair::ev_setup(int eflag, int vflag)
if (eflag_atom && atom->nmax > maxeatom) {
maxeatom = atom->nmax;
memory->destroy(eatom);
- memory->create(eatom,maxeatom,"pair:eatom");
+ memory->create(eatom,comm->nthreads*maxeatom,"pair:eatom");
}
if (vflag_atom && atom->nmax > maxvatom) {
maxvatom = atom->nmax;
memory->destroy(vatom);
- memory->create(vatom,maxvatom,6,"pair:vatom");
+ memory->create(vatom,comm->nthreads*maxvatom,6,"pair:vatom");
}
// zero accumulators
@@ -1151,7 +1151,7 @@ void Pair::init_bitmap(double inner, double outer, int ntablebits,
double Pair::memory_usage()
{
- double bytes = maxeatom * sizeof(double);
- bytes += maxvatom*6 * sizeof(double);
+ double bytes = comm->nthreads*maxeatom * sizeof(double);
+ bytes += comm->nthreads*maxvatom*6 * sizeof(double);
return bytes;
}
diff --git a/src/pair.h b/src/pair.h
index 7e396bbc50..09800c7d3c 100644
--- a/src/pair.h
+++ b/src/pair.h
@@ -20,6 +20,7 @@ namespace LAMMPS_NS {
class Pair : protected Pointers {
friend class BondQuartic;
+ friend class BondQuarticOMP;
friend class DihedralCharmm;
friend class DihedralCharmmOMP;
friend class FixGPU;
diff --git a/src/respa.cpp b/src/respa.cpp
index c283b6454c..3e54b58424 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -599,6 +599,8 @@ void Respa::force_clear(int newtonflag)
{
int i;
+ if (external_force_clear) return;
+
// clear global force array
// nall includes ghosts only if newton flag is set
diff --git a/src/thermo.h b/src/thermo.h
index eeb758b8e7..d68fff887f 100644
--- a/src/thermo.h
+++ b/src/thermo.h
@@ -38,13 +38,13 @@ class Thermo : protected Pointers {
int evaluate_keyword(char *, double *);
private:
- int me;
-
char *line;
- int nfield,nfield_initial;
char **keyword;
int *vtype;
+ int nfield,nfield_initial;
+ int me;
+
char **format,**format_user;
char *format_float_one_def,*format_float_multi_def;
char *format_int_one_def,*format_int_multi_def;
diff --git a/src/verlet.cpp b/src/verlet.cpp
index 7e38367400..c6d53e5240 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -308,6 +308,8 @@ void Verlet::force_clear()
{
int i;
+ if (external_force_clear) return;
+
// clear force on all particles
// if either newton flag is set, also include ghosts
--
GitLab