From 8c372f6c410dde1df5e4210c08dbbb877cb5de48 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 1 Dec 2011 17:21:19 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7266
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/GPU/pppm_gpu.cpp | 8 +-
src/PERI/fix_peri_neigh.h | 4 +-
src/USER-OMP/fix_peri_neigh_omp.h | 3 +-
src/USER-OMP/pair_peri_lps_omp.cpp | 1 -
src/USER-OMP/pair_peri_pmb_omp.cpp | 1 -
src/angle.cpp | 2 +
src/angle.h | 2 +-
src/bond_hybrid.h | 7 +-
src/comm.cpp | 4 +-
src/input.cpp | 36 +++--
src/main.cpp | 9 ++
src/memory.cpp | 5 +
src/min.cpp | 43 ++----
src/neigh_list.cpp | 5 +-
src/neigh_list.h | 2 +-
src/neighbor.cpp | 231 ++++++++++++++++++++---------
src/neighbor.h | 36 +++++
src/verlet.cpp | 43 ++----
18 files changed, 272 insertions(+), 170 deletions(-)
diff --git a/src/GPU/pppm_gpu.cpp b/src/GPU/pppm_gpu.cpp
index 0fa91a2518..13af9fcfbe 100644
--- a/src/GPU/pppm_gpu.cpp
+++ b/src/GPU/pppm_gpu.cpp
@@ -179,8 +179,8 @@ void PPPMGPU::compute(int eflag, int vflag)
// calculate the force on my particles
- FFT_SCALAR qqrd2e_scale=qqrd2e*scale;
- PPPM_GPU_API(interp)(qqrd2e_scale);
+ FFT_SCALAR qscale = force->qqrd2e * scale;
+ PPPM_GPU_API(interp)(qscale);
// sum energy across procs and add in volume-dependent term
@@ -192,7 +192,7 @@ void PPPMGPU::compute(int eflag, int vflag)
energy *= 0.5*volume;
energy -= g_ewald*qsqsum/1.772453851 +
MY_PI2*qsum*qsum / (g_ewald*g_ewald*volume);
- energy *= qqrd2e*scale;
+ energy *= qscale;
}
// sum virial across procs
@@ -200,7 +200,7 @@ void PPPMGPU::compute(int eflag, int vflag)
if (vflag) {
double virial_all[6];
MPI_Allreduce(virial,virial_all,6,MPI_DOUBLE,MPI_SUM,world);
- for (i = 0; i < 6; i++) virial[i] = 0.5*qqrd2e*scale*volume*virial_all[i];
+ for (i = 0; i < 6; i++) virial[i] = 0.5*qscale*volume*virial_all[i];
}
// 2d slab correction
diff --git a/src/PERI/fix_peri_neigh.h b/src/PERI/fix_peri_neigh.h
index 8a441f22cb..c823403a46 100644
--- a/src/PERI/fix_peri_neigh.h
+++ b/src/PERI/fix_peri_neigh.h
@@ -33,7 +33,7 @@ class FixPeriNeigh : public Fix {
public:
FixPeriNeigh(class LAMMPS *,int, char **);
- ~FixPeriNeigh();
+ virtual ~FixPeriNeigh();
int setmask();
void init();
void init_list(int, class NeighList *);
@@ -55,7 +55,7 @@ class FixPeriNeigh : public Fix {
void unpack_comm(int, int, double *);
- private:
+ protected:
int first; // flag for first time initialization
int maxpartner; // max # of peridynamic neighs for any atom
int *npartner; // # of neighbors for each atom
diff --git a/src/USER-OMP/fix_peri_neigh_omp.h b/src/USER-OMP/fix_peri_neigh_omp.h
index 605bbae5f4..11a383a2fe 100644
--- a/src/USER-OMP/fix_peri_neigh_omp.h
+++ b/src/USER-OMP/fix_peri_neigh_omp.h
@@ -27,7 +27,8 @@ namespace LAMMPS_NS {
class FixPeriNeighOMP : public FixPeriNeigh {
public:
- FixPeriNeighOMP(class LAMMPS *lmp, int narg, char **argv) : FixPeriNeigh(lmp,narg,argv) {};
+ FixPeriNeighOMP(class LAMMPS *lmp, int narg, char **argv) :
+ FixPeriNeigh(lmp,narg,argv) {};
virtual void init();
};
diff --git a/src/USER-OMP/pair_peri_lps_omp.cpp b/src/USER-OMP/pair_peri_lps_omp.cpp
index e052271e4f..a4479180bf 100644
--- a/src/USER-OMP/pair_peri_lps_omp.cpp
+++ b/src/USER-OMP/pair_peri_lps_omp.cpp
@@ -37,7 +37,6 @@ PairPeriLPSOMP::PairPeriLPSOMP(LAMMPS *lmp) :
PairPeriLPS(lmp), ThrOMP(lmp, THR_PAIR)
{
respa_enable = 0;
- fix_name = "PERI_NEIGH_OMP";
}
/* ---------------------------------------------------------------------- */
diff --git a/src/USER-OMP/pair_peri_pmb_omp.cpp b/src/USER-OMP/pair_peri_pmb_omp.cpp
index 96e991bab6..de6387ae6b 100644
--- a/src/USER-OMP/pair_peri_pmb_omp.cpp
+++ b/src/USER-OMP/pair_peri_pmb_omp.cpp
@@ -35,7 +35,6 @@ PairPeriPMBOMP::PairPeriPMBOMP(LAMMPS *lmp) :
PairPeriPMB(lmp), ThrOMP(lmp, THR_PAIR)
{
respa_enable = 0;
- fix_name = "PERI_NEIGH_OMP";
}
/* ---------------------------------------------------------------------- */
diff --git a/src/angle.cpp b/src/angle.cpp
index 54a3a5d051..709ace1826 100644
--- a/src/angle.cpp
+++ b/src/angle.cpp
@@ -53,6 +53,8 @@ void Angle::init()
if (!allocated) error->all(FLERR,"Angle coeffs are not set");
for (int i = 1; i <= atom->nangletypes; i++)
if (setflag[i] == 0) error->all(FLERR,"All angle coeffs are not set");
+
+ init_style();
}
/* ----------------------------------------------------------------------
diff --git a/src/angle.h b/src/angle.h
index 7aa1fe85e2..95a35a5a5a 100644
--- a/src/angle.h
+++ b/src/angle.h
@@ -31,10 +31,10 @@ class Angle : protected Pointers {
Angle(class LAMMPS *);
virtual ~Angle();
virtual void init();
- virtual void init_style() {}
virtual void compute(int, int) = 0;
virtual void settings(int, char **) {}
virtual void coeff(int, char **) = 0;
+ virtual void init_style() {};
virtual double equilibrium_angle(int) = 0;
virtual void write_restart(FILE *) = 0;
virtual void read_restart(FILE *) = 0;
diff --git a/src/bond_hybrid.h b/src/bond_hybrid.h
index acdee75063..6c607ec6c4 100644
--- a/src/bond_hybrid.h
+++ b/src/bond_hybrid.h
@@ -29,6 +29,10 @@ class BondHybrid : public Bond {
friend class Force;
public:
+ int nstyles; // # of different bond styles
+ Bond **styles; // class list for each Bond style
+ char **keywords; // keyword for each Bond style
+
BondHybrid(class LAMMPS *);
~BondHybrid();
void compute(int, int);
@@ -42,9 +46,6 @@ class BondHybrid : public Bond {
double memory_usage();
private:
- int nstyles; // # of different bond styles
- Bond **styles; // class list for each Bond style
- char **keywords; // keyword for each Bond style
int *map; // which style each bond type points to
int *nbondlist; // # of bonds in sub-style bondlists
diff --git a/src/comm.cpp b/src/comm.cpp
index e28baab37f..1651e0dfd7 100644
--- a/src/comm.cpp
+++ b/src/comm.cpp
@@ -174,9 +174,9 @@ void Comm::set_procs()
if (domain->triclinic) domain->set_lamda_box();
if (me == 0) {
- if (screen) fprintf(screen," %d by %d by %d processor grid\n",
+ if (screen) fprintf(screen," %d by %d by %d MPI processor grid\n",
procgrid[0],procgrid[1],procgrid[2]);
- if (logfile) fprintf(logfile," %d by %d by %d processor grid\n",
+ if (logfile) fprintf(logfile," %d by %d by %d MPI processor grid\n",
procgrid[0],procgrid[1],procgrid[2]);
}
}
diff --git a/src/input.cpp b/src/input.cpp
index 457faf4ac8..3af9062d05 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -509,6 +509,14 @@ void Input::clear()
if (narg > 0) error->all(FLERR,"Illegal clear command");
lmp->destroy();
lmp->create();
+
+ // use of "omp" suffix implies using "package omp"
+ // re-create the effect of that package command which creates a fix
+
+ if (lmp->suffix && lmp->suffix_enable) {
+ if (strcmp(lmp->suffix,"omp") == 0)
+ lmp->input->one("package omp *");
+ }
}
/* ---------------------------------------------------------------------- */
@@ -1141,25 +1149,15 @@ void Input::package()
delete [] fixarg;
} else if (strcmp(arg[0],"omp") == 0) {
-
-#ifdef _OPENMP
- if (narg != 2) error->all(FLERR,"Illegal package command");
- comm->nthreads = atoi(arg[1]);
- if (comm->nthreads < 1) error->all(FLERR,"Illegal package command");
-
- omp_set_num_threads(comm->nthreads);
- if (me == 0) {
- if (screen)
- fprintf(screen," reset %d OpenMP thread(s) per MPI task\n",
- comm->nthreads);
- if (logfile)
- fprintf(logfile," reset %d OpenMP thread(s) per MPI task\n",
- comm->nthreads);
- }
-#else
- error->all(FLERR,"Cannot use package omp command with no OpenMP support");
-#endif
-
+ char **fixarg = new char*[2+narg];
+ fixarg[0] = "package_omp";
+ fixarg[1] = "all";
+ fixarg[2] = "OMP";
+ for (int i = 1; i < narg; i++) fixarg[i+2] = arg[i];
+ modify->allow_early_fix = 1;
+ modify->add_fix(2+narg,fixarg,NULL);
+ modify->allow_early_fix = 0;
+ delete [] fixarg;
} else error->all(FLERR,"Illegal package command");
}
diff --git a/src/main.cpp b/src/main.cpp
index a3764872ec..0776e7e977 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -14,6 +14,7 @@
#include "mpi.h"
#include "lammps.h"
#include "input.h"
+#include "string.h"
using namespace LAMMPS_NS;
@@ -26,6 +27,14 @@ int main(int argc, char **argv)
MPI_Init(&argc,&argv);
LAMMPS *lammps = new LAMMPS(argc,argv,MPI_COMM_WORLD);
+
+ // using /omp suffix implies running the "package omp"
+ // command with default settings.
+ if (lammps->suffix && lammps->suffix_enable) {
+ if (strcmp(lammps->suffix,"omp") == 0)
+ lammps->input->one("package omp *");
+ }
+
lammps->input->file();
delete lammps;
diff --git a/src/memory.cpp b/src/memory.cpp
index 11e15f10e5..c98f40fd3f 100644
--- a/src/memory.cpp
+++ b/src/memory.cpp
@@ -32,7 +32,12 @@ void *Memory::smalloc(bigint nbytes, const char *name)
{
if (nbytes == 0) return NULL;
+#if defined(LAMMPS_MEMALIGN)
+ void *ptr;
+ posix_memalign(&ptr, LAMMPS_MEMALIGN, nbytes);
+#else
void *ptr = malloc(nbytes);
+#endif
if (ptr == NULL) {
char str[128];
sprintf(str,"Failed to allocate " BIGINT_FORMAT " bytes for array %s",
diff --git a/src/min.cpp b/src/min.cpp
index f380a103d5..6322b96f77 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -135,6 +135,11 @@ void Min::init()
ev_setup();
+ // detect if fix omp is present for clearing force arrays
+
+ int ifix = modify->find_fix("package_omp");
+ if (ifix >= 0) external_force_clear = 1;
+
// set flags for what arrays to clear in force_clear()
// need to clear additionals arrays if they exist
@@ -513,6 +518,8 @@ double Min::energy_force(int resetflag)
void Min::force_clear()
{
+ if (external_force_clear) return;
+
int i;
if (external_force_clear) return;
@@ -523,38 +530,14 @@ void Min::force_clear()
int nall;
if (force->newton) nall = atom->nlocal + atom->nghost;
else nall = atom->nlocal;
- int ntot = nall * comm->nthreads;
- double **f = atom->f;
- for (i = 0; i < ntot; i++) {
- f[i][0] = 0.0;
- f[i][1] = 0.0;
- f[i][2] = 0.0;
- }
+ size_t nbytes = sizeof(double) * nall;
- if (torqueflag) {
- double **torque = atom->torque;
- for (i = 0; i < nall; i++) {
- torque[i][0] = 0.0;
- torque[i][1] = 0.0;
- torque[i][2] = 0.0;
- }
- }
-
- if (erforceflag) {
- double *erforce = atom->erforce;
- for (i = 0; i < nall; i++) erforce[i] = 0.0;
- }
-
- if (e_flag) {
- double *de = atom->de;
- for (i = 0; i < nall; i++) de[i] = 0.0;
- }
-
- if (rho_flag) {
- double *drho = atom->drho;
- for (i = 0; i < nall; i++) drho[i] = 0.0;
- }
+ memset(&(atom->f[0][0]),0,3*nbytes);
+ if (torqueflag) memset(&(atom->torque[0][0]),0,3*nbytes);
+ if (erforceflag) memset(&(atom->erforce[0]), 0, nbytes);
+ if (e_flag) memset(&(atom->de[0]), 0, nbytes);
+ if (rho_flag) memset(&(atom->drho[0]), 0, nbytes);
}
/* ----------------------------------------------------------------------
diff --git a/src/neigh_list.cpp b/src/neigh_list.cpp
index 13dc4a217f..963ee03cf5 100644
--- a/src/neigh_list.cpp
+++ b/src/neigh_list.cpp
@@ -181,10 +181,10 @@ void NeighList::stencil_allocate(int smax, int style)
add PGDELTA pages to neighbor list
------------------------------------------------------------------------- */
-int **NeighList::add_pages()
+int **NeighList::add_pages(int howmany)
{
int toppage = maxpage;
- maxpage += PGDELTA;
+ maxpage += howmany*PGDELTA;
pages = (int **)
memory->srealloc(pages,maxpage*sizeof(int *),"neighlist:pages");
@@ -249,6 +249,7 @@ void NeighList::print_attributes()
printf("\n");
printf(" %d = occasional\n",rq->occasional);
printf(" %d = dnum\n",rq->dnum);
+ printf(" %d = omp\n",rq->omp);
printf(" %d = ghost\n",rq->ghost);
printf(" %d = cudable\n",rq->cudable);
printf(" %d = omp\n",rq->omp);
diff --git a/src/neigh_list.h b/src/neigh_list.h
index 63675503cf..a8857811a3 100644
--- a/src/neigh_list.h
+++ b/src/neigh_list.h
@@ -79,7 +79,7 @@ class NeighList : protected Pointers {
~NeighList();
void grow(int); // grow maxlocal
void stencil_allocate(int, int); // allocate stencil arrays
- int **add_pages(); // add pages to neigh list
+ int **add_pages(int howmany=1); // add pages to neigh list
void copy_skip_info(int *, int **); // copy skip info from a neigh request
void print_attributes(); // debug routine
int get_maxlocal() {return maxatoms;}
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index 24e4285b50..46b0a1e0c8 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -808,85 +808,170 @@ void Neighbor::choose_build(int index, NeighRequest *rq)
{
PairPtr pb = NULL;
- if (rq->copy) pb = &Neighbor::copy_from;
+ if (rq->omp == 0) {
+
+ if (rq->copy) pb = &Neighbor::copy_from;
+
+ else if (rq->skip) {
+ if (rq->gran && lists[index]->listgranhistory)
+ pb = &Neighbor::skip_from_granular;
+ else if (rq->respaouter) pb = &Neighbor::skip_from_respa;
+ else pb = &Neighbor::skip_from;
+
+ } else if (rq->half_from_full) {
+ if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton;
+ else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton;
+
+ } else if (rq->half) {
+ if (style == NSQ) {
+ if (rq->newton == 0) {
+ if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton;
+ else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton;
+ } else if (rq->newton == 1) {
+ pb = &Neighbor::half_nsq_newton;
+ } else if (rq->newton == 2) {
+ pb = &Neighbor::half_nsq_no_newton;
+ }
+ } else if (style == BIN) {
+ if (rq->newton == 0) {
+ if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton;
+ else if (triclinic == 0) pb = &Neighbor::half_bin_newton;
+ else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri;
+ } else if (rq->newton == 1) {
+ if (triclinic == 0) pb = &Neighbor::half_bin_newton;
+ else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri;
+ } else if (rq->newton == 2) pb = &Neighbor::half_bin_no_newton;
+ } else if (style == MULTI) {
+ if (rq->newton == 0) {
+ if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton;
+ else if (triclinic == 0) pb = &Neighbor::half_multi_newton;
+ else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri;
+ } else if (rq->newton == 1) {
+ if (triclinic == 0) pb = &Neighbor::half_multi_newton;
+ else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri;
+ } else if (rq->newton == 2) pb = &Neighbor::half_multi_no_newton;
+ }
- else if (rq->skip) {
- if (rq->gran && lists[index]->listgranhistory)
- pb = &Neighbor::skip_from_granular;
- else if (rq->respaouter) pb = &Neighbor::skip_from_respa;
- else pb = &Neighbor::skip_from;
+ } else if (rq->full) {
+ if (style == NSQ) {
+ if (rq->ghost == 0) pb = &Neighbor::full_nsq;
+ else if (includegroup)
+ error->all(FLERR,"Neighbor include group not allowed with ghost neighbors");
+ else if (rq->ghost == 1) pb = &Neighbor::full_nsq_ghost;
+ } else if (style == BIN) {
+ if (rq->ghost == 0) pb = &Neighbor::full_bin;
+ else if (includegroup)
+ error->all(FLERR,"Neighbor include group not allowed with ghost neighbors");
+ else if (rq->ghost == 1) pb = &Neighbor::full_bin_ghost;
+ } else if (style == MULTI) {
+ if (rq->ghost == 0) pb = &Neighbor::full_multi;
+ else error->all(FLERR,"Neighbor multi not yet enabled for ghost neighbors");
+ }
- } else if (rq->half_from_full) {
- if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton;
- else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton;
+ } else if (rq->gran) {
+ if (style == NSQ) {
+ if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton;
+ else if (newton_pair == 1) pb = &Neighbor::granular_nsq_newton;
+ } else if (style == BIN) {
+ if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton;
+ else if (triclinic == 0) pb = &Neighbor::granular_bin_newton;
+ else if (triclinic == 1) pb = &Neighbor::granular_bin_newton_tri;
+ } else if (style == MULTI)
+ error->all(FLERR,"Neighbor multi not yet enabled for granular");
+
+ } else if (rq->respaouter) {
+ if (style == NSQ) {
+ if (newton_pair == 0) pb = &Neighbor::respa_nsq_no_newton;
+ else if (newton_pair == 1) pb = &Neighbor::respa_nsq_newton;
+ } else if (style == BIN) {
+ if (newton_pair == 0) pb = &Neighbor::respa_bin_no_newton;
+ else if (triclinic == 0) pb = &Neighbor::respa_bin_newton;
+ else if (triclinic == 1) pb = &Neighbor::respa_bin_newton_tri;
+ } else if (style == MULTI)
+ error->all(FLERR,"Neighbor multi not yet enabled for rRESPA");
+ }
+ } else {
- } else if (rq->half) {
- if (style == NSQ) {
- if (rq->newton == 0) {
- if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton;
- else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton;
- } else if (rq->newton == 1) {
- pb = &Neighbor::half_nsq_newton;
- } else if (rq->newton == 2) {
- pb = &Neighbor::half_nsq_no_newton;
+ if (rq->copy) pb = &Neighbor::copy_from;
+
+ else if (rq->skip) {
+ if (rq->gran && lists[index]->listgranhistory)
+ pb = &Neighbor::skip_from_granular;
+ else if (rq->respaouter) pb = &Neighbor::skip_from_respa;
+ else pb = &Neighbor::skip_from;
+
+ } else if (rq->half_from_full) {
+ if (newton_pair == 0) pb = &Neighbor::half_from_full_no_newton_omp;
+ else if (newton_pair == 1) pb = &Neighbor::half_from_full_newton_omp;
+
+ } else if (rq->half) {
+ if (style == NSQ) {
+ if (rq->newton == 0) {
+ if (newton_pair == 0) pb = &Neighbor::half_nsq_no_newton_omp;
+ else if (newton_pair == 1) pb = &Neighbor::half_nsq_newton_omp;
+ } else if (rq->newton == 1) {
+ pb = &Neighbor::half_nsq_newton_omp;
+ } else if (rq->newton == 2) {
+ pb = &Neighbor::half_nsq_no_newton_omp;
+ }
+ } else if (style == BIN) {
+ if (rq->newton == 0) {
+ if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton_omp;
+ else if (triclinic == 0) pb = &Neighbor::half_bin_newton_omp;
+ else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri_omp;
+ } else if (rq->newton == 1) {
+ if (triclinic == 0) pb = &Neighbor::half_bin_newton_omp;
+ else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri_omp;
+ } else if (rq->newton == 2) pb = &Neighbor::half_bin_no_newton_omp;
+ } else if (style == MULTI) {
+ if (rq->newton == 0) {
+ if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton_omp;
+ else if (triclinic == 0) pb = &Neighbor::half_multi_newton_omp;
+ else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri_omp;
+ } else if (rq->newton == 1) {
+ if (triclinic == 0) pb = &Neighbor::half_multi_newton_omp;
+ else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri_omp;
+ } else if (rq->newton == 2) pb = &Neighbor::half_multi_no_newton_omp;
}
- } else if (style == BIN) {
- if (rq->newton == 0) {
- if (newton_pair == 0) pb = &Neighbor::half_bin_no_newton;
- else if (triclinic == 0) pb = &Neighbor::half_bin_newton;
- else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri;
- } else if (rq->newton == 1) {
- if (triclinic == 0) pb = &Neighbor::half_bin_newton;
- else if (triclinic == 1) pb = &Neighbor::half_bin_newton_tri;
- } else if (rq->newton == 2) pb = &Neighbor::half_bin_no_newton;
- } else if (style == MULTI) {
- if (rq->newton == 0) {
- if (newton_pair == 0) pb = &Neighbor::half_multi_no_newton;
- else if (triclinic == 0) pb = &Neighbor::half_multi_newton;
- else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri;
- } else if (rq->newton == 1) {
- if (triclinic == 0) pb = &Neighbor::half_multi_newton;
- else if (triclinic == 1) pb = &Neighbor::half_multi_newton_tri;
- } else if (rq->newton == 2) pb = &Neighbor::half_multi_no_newton;
- }
- } else if (rq->full) {
- if (style == NSQ) {
- if (rq->ghost == 0) pb = &Neighbor::full_nsq;
- else if (includegroup)
- error->all(FLERR,"Neighbor include group not allowed with ghost neighbors");
- else if (rq->ghost == 1) pb = &Neighbor::full_nsq_ghost;
- } else if (style == BIN) {
- if (rq->ghost == 0) pb = &Neighbor::full_bin;
- else if (includegroup)
- error->all(FLERR,"Neighbor include group not allowed with ghost neighbors");
- else if (rq->ghost == 1) pb = &Neighbor::full_bin_ghost;
- } else if (style == MULTI) {
- if (rq->ghost == 0) pb = &Neighbor::full_multi;
- else error->all(FLERR,"Neighbor multi not yet enabled for ghost neighbors");
- }
+ } else if (rq->full) {
+ if (style == NSQ) {
+ if (rq->ghost == 0) pb = &Neighbor::full_nsq_omp;
+ else if (includegroup)
+ error->all(FLERR,"Neighbor include group not allowed with ghost neighbors");
+ else if (rq->ghost == 1) pb = &Neighbor::full_nsq_ghost_omp;
+ } else if (style == BIN) {
+ if (rq->ghost == 0) pb = &Neighbor::full_bin_omp;
+ else if (includegroup)
+ error->all(FLERR,"Neighbor include group not allowed with ghost neighbors");
+ else if (rq->ghost == 1) pb = &Neighbor::full_bin_ghost_omp;
+ } else if (style == MULTI) {
+ if (rq->ghost == 0) pb = &Neighbor::full_multi_omp;
+ else error->all(FLERR,"Neighbor multi not yet enabled for ghost neighbors");
+ }
- } else if (rq->gran) {
- if (style == NSQ) {
- if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton;
- else if (newton_pair == 1) pb = &Neighbor::granular_nsq_newton;
- } else if (style == BIN) {
- if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton;
- else if (triclinic == 0) pb = &Neighbor::granular_bin_newton;
- else if (triclinic == 1) pb = &Neighbor::granular_bin_newton_tri;
- } else if (style == MULTI)
- error->all(FLERR,"Neighbor multi not yet enabled for granular");
-
- } else if (rq->respaouter) {
- if (style == NSQ) {
- if (newton_pair == 0) pb = &Neighbor::respa_nsq_no_newton;
- else if (newton_pair == 1) pb = &Neighbor::respa_nsq_newton;
- } else if (style == BIN) {
- if (newton_pair == 0) pb = &Neighbor::respa_bin_no_newton;
- else if (triclinic == 0) pb = &Neighbor::respa_bin_newton;
- else if (triclinic == 1) pb = &Neighbor::respa_bin_newton_tri;
- } else if (style == MULTI)
- error->all(FLERR,"Neighbor multi not yet enabled for rRESPA");
+ } else if (rq->gran) {
+ if (style == NSQ) {
+ if (newton_pair == 0) pb = &Neighbor::granular_nsq_no_newton_omp;
+ else if (newton_pair == 1) pb = &Neighbor::granular_nsq_newton_omp;
+ } else if (style == BIN) {
+ if (newton_pair == 0) pb = &Neighbor::granular_bin_no_newton_omp;
+ else if (triclinic == 0) pb = &Neighbor::granular_bin_newton_omp;
+ else if (triclinic == 1) pb = &Neighbor::granular_bin_newton_tri_omp;
+ } else if (style == MULTI)
+ error->all(FLERR,"Neighbor multi not yet enabled for granular");
+
+ } else if (rq->respaouter) {
+ if (style == NSQ) {
+ if (newton_pair == 0) pb = &Neighbor::respa_nsq_no_newton_omp;
+ else if (newton_pair == 1) pb = &Neighbor::respa_nsq_newton_omp;
+ } else if (style == BIN) {
+ if (newton_pair == 0) pb = &Neighbor::respa_bin_no_newton_omp;
+ else if (triclinic == 0) pb = &Neighbor::respa_bin_newton_omp;
+ else if (triclinic == 1) pb = &Neighbor::respa_bin_newton_tri_omp;
+ } else if (style == MULTI)
+ error->all(FLERR,"Neighbor multi not yet enabled for rRESPA");
+ }
}
// general error check
diff --git a/src/neighbor.h b/src/neighbor.h
index 5248c0bd6e..aa2fcbcba1 100644
--- a/src/neighbor.h
+++ b/src/neighbor.h
@@ -75,6 +75,8 @@ class Neighbor : protected Pointers {
void set(int, char **); // set neighbor style and skin distance
void modify_params(int, char**); // modify parameters that control builds
bigint memory_usage();
+
+ inline int exclude_setting() { return exclude; }
protected:
int me,nprocs;
@@ -209,6 +211,40 @@ class Neighbor : protected Pointers {
void respa_bin_newton(class NeighList *);
void respa_bin_newton_tri(class NeighList *);
+ // OpenMP multi-threaded neighbor list build versions
+
+ void half_nsq_no_newton_omp(class NeighList *);
+ void half_nsq_newton_omp(class NeighList *);
+
+ void half_bin_no_newton_omp(class NeighList *);
+ void half_bin_newton_omp(class NeighList *);
+ void half_bin_newton_tri_omp(class NeighList *);
+
+ void half_multi_no_newton_omp(class NeighList *);
+ void half_multi_newton_omp(class NeighList *);
+ void half_multi_newton_tri_omp(class NeighList *);
+
+ void full_nsq_omp(class NeighList *);
+ void full_nsq_ghost_omp(class NeighList *);
+ void full_bin_omp(class NeighList *);
+ void full_bin_ghost_omp(class NeighList *);
+ void full_multi_omp(class NeighList *);
+
+ void half_from_full_no_newton_omp(class NeighList *);
+ void half_from_full_newton_omp(class NeighList *);
+
+ void granular_nsq_no_newton_omp(class NeighList *);
+ void granular_nsq_newton_omp(class NeighList *);
+ void granular_bin_no_newton_omp(class NeighList *);
+ void granular_bin_newton_omp(class NeighList *);
+ void granular_bin_newton_tri_omp(class NeighList *);
+
+ void respa_nsq_no_newton_omp(class NeighList *);
+ void respa_nsq_newton_omp(class NeighList *);
+ void respa_bin_no_newton_omp(class NeighList *);
+ void respa_bin_newton_omp(class NeighList *);
+ void respa_bin_newton_tri_omp(class NeighList *);
+
// pairwise stencil creation functions
typedef void (Neighbor::*StencilPtr)(class NeighList *, int, int, int);
diff --git a/src/verlet.cpp b/src/verlet.cpp
index c6d53e5240..284f9726a4 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -62,6 +62,11 @@ void Verlet::init()
ev_setup();
+ // detect if fix omp is present for clearing force arrays
+
+ int ifix = modify->find_fix("package_omp");
+ if (ifix >= 0) external_force_clear = 1;
+
// set flags for what arrays to clear in force_clear()
// need to clear additionals arrays if they exist
@@ -306,49 +311,27 @@ void Verlet::cleanup()
void Verlet::force_clear()
{
+ if (external_force_clear) return;
int i;
if (external_force_clear) return;
// clear force on all particles
// if either newton flag is set, also include ghosts
+ // when using threads always clear all forces.
if (neighbor->includegroup == 0) {
int nall;
if (force->newton) nall = atom->nlocal + atom->nghost;
else nall = atom->nlocal;
- int ntot = nall * comm->nthreads;
- double **f = atom->f;
- for (i = 0; i < ntot; i++) {
- f[i][0] = 0.0;
- f[i][1] = 0.0;
- f[i][2] = 0.0;
- }
-
- if (torqueflag) {
- double **torque = atom->torque;
- for (i = 0; i < nall; i++) {
- torque[i][0] = 0.0;
- torque[i][1] = 0.0;
- torque[i][2] = 0.0;
- }
- }
-
- if (erforceflag) {
- double *erforce = atom->erforce;
- for (i = 0; i < nall; i++) erforce[i] = 0.0;
- }
+ size_t nbytes = sizeof(double) * nall;
- if (e_flag) {
- double *de = atom->de;
- for (i = 0; i < nall; i++) de[i] = 0.0;
- }
-
- if (rho_flag) {
- double *drho = atom->drho;
- for (i = 0; i < nall; i++) drho[i] = 0.0;
- }
+ memset(&(atom->f[0][0]),0,3*nbytes);
+ if (torqueflag) memset(&(atom->torque[0][0]),0,3*nbytes);
+ if (erforceflag) memset(&(atom->erforce[0]), 0, nbytes);
+ if (e_flag) memset(&(atom->de[0]), 0, nbytes);
+ if (rho_flag) memset(&(atom->drho[0]), 0, nbytes);
// neighbor includegroup flag is set
// clear force only on initial nfirst particles
--
GitLab