diff --git a/src/ASPHERE/pair_line_lj.cpp b/src/ASPHERE/pair_line_lj.cpp
index ca50168be5ebfff17e18b44d349bc560f7288837..2f3f602ba78b9390b272f4f621011ae10fa6bf84 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 b30140acea1ddb8818392bfc5cc233295dd6db3d..4b60c80eafd71c153a49dd15c77f71956e7f8c2f 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 13e57fa64ee3ecc8d7be3ab38c80eb0792ba533c..9fb782b681f1138f262f390c4794229567a32884 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 0ecf19aee5397bb1448968b512d7b6369139f985..dfcdd819ef501fb45fb3ddbeabe64bcc84105caf 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 a8cf2bf7218f2d9821a770f3f3117bd929334df5..fb0028cce0d585906156c95d4ab11c2138b3bf35 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 cbea3cfe7e4f93956ec569ed5c593820ad3fa437..f24ecce524b999a899fa1a043ee3c49f3e8a1481 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 2a5a42c6482dd383837bc65a4b15605462cdd5ef..3c201e923fbccbcf11f7c3fe9d6ee01efd769ed3 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 9a86e9cd8645de80ad85c8ecd91d6273b545092d..cbee3ac8105b4ddf5defc8572f0236e1a0e4a517 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 a44498670117faaba9d3112f68d2c526b6a982ef..eec1391bfc14ec6d6b294aa8c437a517d4a0497f 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 03a1c3a62ac762924c8de59c90ab55ae5e49d86c..fe0ce40f315ed6dd11efbe53f54c5d8fb631c602 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 ae260e478ee0054ddfc2fe6bd10a3e0dfe89390b..343511c8609e333f3e0b4d3858a1777e86fb48a8 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 edd030b55b44ed33e117c3ddc33a2822654071b0..b6cff3d932da7c5abee6d6362fbb78ae4732dc52 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 95482cfbfee201e01b757d1711f15f3b4524d4dc..cd22f987f733ec5308c7f0a555686174f6829059 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 bc3cfc98dd2e0e5e9531c80283cb818440d3dc5c..ece15c68397ffa7465d61c3c77ac2bb8a60bac12 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 221778e6ee7092fd299d663e3b033a3065a5a386..e2d017a6e5f7258b7bbec2be1a5ad5f02b47cfc0 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 5f57715f9f89c95ebcfe64e2989bdafb0d818d7a..daffb5ff2dd10f3d995ca8a6911407fd3a03a8bc 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 0f6b0a6a9db24ff2a35c515be731e877cfa9c00d..952c1b73a5cedcdb6916b8ea0c5b7fc628ee244b 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 81d23867b66663c0113c06818c8628af11615c65..affd5f1ec39dbf94ecf6c2ebb5d5af2634beb15f 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 1f4aabe456d541d35c1df7443f4599e4f85bc2a1..6b73cf50f1fe68ce557f66fe1f253abca85c7450 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 a6848cc9a505c608b69d743ea41fab1873106c03..52399422fcd5e60dd406f066f167b86a87e5785d 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 4fb901ebe5f5fe026fd1afeaaecc8aacbf62f0f6..fa0a56ab89ec5b5a5f8a7c0693710363fb9ff1bc 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 f91a84f7fae41ac8838e4fbd347dd6facfa5263d..7c13d926fb4681f62e8742ef8f6c34b031aaa5d1 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 23ce4874ca337ec7a9233c968981b5e6a58f5425..b4044083494c69724bdb09d4b3e70cc41c448acd 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 25ca928e3190d4a8aba70f526445a03417d99067..6fb262c332fe7576a27c49e6e5181e12a0d83953 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 a86dfdba1420c303c74444df9383edf6c0ad568f..89f28c4bf363e290b68c0140df9bb89350862e17 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 a4364cc8bb397ef57ff274682b7ced3a7fd17df1..87eefd6a0dcaf90a92f3dcb7c40a77d21b16b7fe 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 34fbf265cc1919f03ddb06216c412bf654e012ac..1dcf3653e15e7a4e1b2d08ee83a39444b58be041 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 1390d2519385a5433f15256c4ba98b9bf7b2550a..8d715e2e6158d081c84128177599128b7eea2cfe 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 67ae294dec646f69c35f84bd46d82921c0277879..628d2c237bf1069f1a44fbc68f626747a65626e3 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 f0d289f0552e11053331239babc15a00aa245e4f..810693262f863877670ecff42d8da579cb48d65e 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 ccef917d3a15ef73eceb5c186ec60aebf0691c2a..86de363a15552009312a81a0c7769503cf6e7a60 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 a14e222d60f483c947075f37c1952846313010ba..b3d5ef81f6df0c0cbfab857ea54ee4b512f031a2 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 cef7583743950469c825f7aa24a10e9bfd1247aa..483979826eb7b7d8ea8e9c31355c9bb2d7f8dbdb 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 24db63e5795f86375ae17d9cd5131f8dc5bc95b7..aee802db21693e62744288d6315d9dde1d45a271 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 2d610ea2e0937da6837f3026229a4d65ab9648b8..3b66583715d28cd9b76c05d9f947cad8bce42ee0 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 ce44ea3045e01b0689f2d813ef59547495bbaabc..d8b7ea17411a9b5fdb42bc8916422998309813a1 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 cd6a4960d93228ce4368d966f188c11bd8487c0a..031cfec5635df8fb4e4372394f0b5608ff80c6c6 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 c51ece7c44bfd2ad67967b41378090061ab43a1e..50883dc0346551b68354855c990edd6745e4bac3 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 91bb723323b804dbac22ee307ec5caf2e24ff881..b3378aea3f8dc0e1bed060be5e67056e32ad105f 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 7bb63ca64b6d277f53b2819b0db4e555debd0c96..6e957d286c413ae9d4a9f16fce15dc90e64415cf 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 58ba679948f3d976008545a0570c64f4c681199a..54a3a5d051bf5655fc27787fc35cac1ff16d4418 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 7f25f3828972be7eaeeb3f0263b89cda98b52251..46ee5e1dc63709e941746b0261910968c27cdc26 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 ebf8b7d7a8b40dfeac0d44a06b50503ba8841df8..c49cd6b6ad4d88bdfb9269538bbfa4c0ce9af53b 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 03aed7aa65694ca157765fcc646690a927f2c9ca..0f3fbe161a3cf4a85f1d05f26c4e4e9615e5a04c 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 92310df6532fc9ab6d12b2e7c1028a5ee3d3f311..aa168ff45cc6eeaa99075b347dd60f042c0f346e 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 551d84b51fff497762eeff9b3f591adb2c77a451..29f9862fb7cd5ea5528f2b6b9e4450dd05a1ddb0 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 974df81b1425e0f26acfb9fb2b71aa26a156c2e0..f1a62cba163d1121b64e333a569ea2907d59a844 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 a1a0045fa94c274bd0b9b855d9a718ea382a5d74..3140ea6a5487fc97799a24b7d49e63bd7dec0971 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 c07808ca028ab0674944735e3fe3d237c6c7d02c..1e04ab4ab9b6b3aa8f4c9d5c907c2af3f0a9edfc 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 67b387330926c2cdc44364e689d9d8b346d5d759..5a26d16aa0178268458a74bc176fbeb4285656b5 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 654b81a3a8888a4087e8a90f1019213b57f8acfa..a7089070d222af385a7ca0214bae3edec078cb15 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 8f5ab4d6e3b5aed7846fd92023f66d519bbe3be0..9a1f47e5ce717f624d91f88b9c5d833837c110b6 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 d8310b2475f37132762a7ec54c417f1a02f52e71..d766015e34f0f154ecd8aaf6ea667068d628da99 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 2251335b4373c6e2486c50785d05c02689bb9946..457faf4ac800e8ce39bdcc7ce0187aea5416cdbd 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 80bf89f676e611ec991d850e94350753988a8e1f..6de4df26c3fa3abf7990c2259e32d754071732dc 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 8f8942c80c3613731b0ad2fadbc010b221879b93..3c75176a105b22975b19e2d09273f1da4269a713 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 3ad2193f11db40c588e58e58da7c96fad480e5b1..f77bacfcdc120370546e6d57aaca143e7b40ab55 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 cf7e8a6112b20e597e5cd564301081cd79b8f988..5cbc2aee28b5a6f80a0d598a01cd276021317e4b 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 bf83065261f530d18abffc98ad34a86911e7345f..0df114edaf8000fa9aa91d15eeeda46599947fd5 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 f367fa7b93783ceef199060f2699aafb52fd248a..47956059d8157077541cbc2c32e260844bf0dbb5 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 b4b7fd0aeb91a5aee67c94722b1c1b79e304acf5..f35acdb4e39811bf925d255ffd5f85c35689e18f 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 dbd8f0d4a00de429dc4829b3d2bb260971dd0b68..738bd5d7e166b56e263fa797711a992970c7a8dc 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 82c9f0b267c0d531eb4959b092ad65540a81b34d..2146009a99039baee432133a0a235eac5fa1b216 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 afff5966c04b0f2755ca79cfc1bf804e3b6b46c1..e3dcc915a7bc0018864d9a0d148bdce8142a9125 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 58f8ad94dd99df7475119dec529c566f20ee031d..24e4285b5050065f1bfa63ce5e7d5d23065518d7 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 c2fa635d63404e76e6591f7ff8faa4fde13469e9..5248c0bd6e862febfc3095d4e83ba30852e4aec8 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 c83907831080092548cd32b1d6b893e0f59a41fe..c921fdebdf596286698db92f6239a088f2d9b126 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 7e396bbc5045e1c159a376180a425545586645f7..09800c7d3cbdcc1ce86f05991c33c86becdde6a3 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 c283b6454cb5abb23cb8dd30db2571ac2bb33f62..3e54b58424616ce5a12c867823b2ec7fe0637348 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 eeb758b8e73908650786505cabb023c7244d372d..d68fff887fe139e1576e1e7181269d2e3592dcbd 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 7e38367400d2cbc4ab94d32ca354c03c81323484..c6d53e52406ab6a2e5ecff1352d20ff3c7c596ec 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