diff --git a/lib/gpu/lal_table.cpp b/lib/gpu/lal_table.cpp index 6e0a3964f0..5b3d934e53 100644 --- a/lib/gpu/lal_table.cpp +++ b/lib/gpu/lal_table.cpp @@ -122,7 +122,7 @@ int TableT::init(const int ntypes, host_write[ix*lj_types+iy].w = (numtyp)0.0; } ucl_copy(coeff2,host_write,false); - + // Allocate tablength arrays UCL_H_Vec host_write2(_ntables*_tablength,*(this->ucl_device), UCL_WRITE_OPTIMIZED); @@ -307,29 +307,29 @@ void TableT::loop(const bool _eflag, const bool _vflag) { &this->ans->dev_engv.begin(), &eflag, &vflag, &ainum, &nbor_pitch, &this->_threads_per_atom, &_tablength); - } + } } else { if (_tabstyle == LOOKUP) { this->k_pair.set_size(GX,BX); this->k_pair.run(&this->atom->dev_x.begin(), &tabindex.begin(), - &coeff2.begin(), &coeff3.begin(), &coeff4.begin(), &cutsq.begin(), - &_lj_types, &sp_lj.begin(), &this->nbor->dev_nbor.begin(), + &coeff2.begin(), &coeff3.begin(), &coeff4.begin(), &_lj_types, + &cutsq.begin(), &sp_lj.begin(), &this->nbor->dev_nbor.begin(), &this->_nbor_data->begin(), &this->ans->dev_ans.begin(), &this->ans->dev_engv.begin(), &eflag, &vflag, &ainum, &nbor_pitch, &this->_threads_per_atom, &_tablength); } else if (_tabstyle == LINEAR) { this->k_pair_linear.set_size(GX,BX); this->k_pair_linear.run(&this->atom->dev_x.begin(), &tabindex.begin(), - &coeff2.begin(), &coeff3.begin(), &coeff4.begin(), &cutsq.begin(), - &_lj_types, &sp_lj.begin(), &this->nbor->dev_nbor.begin(), + &coeff2.begin(), &coeff3.begin(), &coeff4.begin(), &_lj_types, + &cutsq.begin(), &sp_lj.begin(), &this->nbor->dev_nbor.begin(), &this->_nbor_data->begin(), &this->ans->dev_ans.begin(), &this->ans->dev_engv.begin(), &eflag, &vflag, &ainum, &nbor_pitch, &this->_threads_per_atom, &_tablength); } else if (_tabstyle == SPLINE) { this->k_pair_spline.set_size(GX,BX); this->k_pair_spline.run(&this->atom->dev_x.begin(), &tabindex.begin(), - &coeff2.begin(), &coeff3.begin(), &coeff4.begin(), &cutsq.begin(), - &_lj_types, &sp_lj.begin(), &this->nbor->dev_nbor.begin(), + &coeff2.begin(), &coeff3.begin(), &coeff4.begin(), &_lj_types, + &cutsq.begin(), &sp_lj.begin(), &this->nbor->dev_nbor.begin(), &this->_nbor_data->begin(), &this->ans->dev_ans.begin(), &this->ans->dev_engv.begin(), &eflag, &vflag, &ainum, &nbor_pitch, &this->_threads_per_atom, &_tablength); @@ -337,8 +337,8 @@ void TableT::loop(const bool _eflag, const bool _vflag) { this->k_pair_bitmap.set_size(GX,BX); this->k_pair_bitmap.run(&this->atom->dev_x.begin(), &tabindex.begin(), &nshiftbits.begin(), &nmask.begin(), - &coeff2.begin(), &coeff3.begin(), &coeff4.begin(), &cutsq.begin(), - &_lj_types, &sp_lj.begin(), &this->nbor->dev_nbor.begin(), + &coeff2.begin(), &coeff3.begin(), &coeff4.begin(), &_lj_types, + &cutsq.begin(), &sp_lj.begin(), &this->nbor->dev_nbor.begin(), &this->_nbor_data->begin(), &this->ans->dev_ans.begin(), &this->ans->dev_engv.begin(), &eflag, &vflag, &ainum, &nbor_pitch, &this->_threads_per_atom, &_tablength); diff --git a/src/GPU/pair_buck_coul_long_gpu.cpp b/src/GPU/pair_buck_coul_long_gpu.cpp index a5f6401117..fdca85cf1f 100644 --- a/src/GPU/pair_buck_coul_long_gpu.cpp +++ b/src/GPU/pair_buck_coul_long_gpu.cpp @@ -248,7 +248,7 @@ void PairBuckCoulLongGPU::cpu_compute(int start, int inum, int eflag, if (rsq < cutsq[itype][jtype]) { r2inv = 1.0/rsq; - + r = sqrt(rsq); if (rsq < cut_coulsq) { grij = g_ewald * r; expm2 = exp(-grij*grij); diff --git a/src/STUBS/mpi.c b/src/STUBS/mpi.c index ea383ef722..7f3dd3f8ef 100644 --- a/src/STUBS/mpi.c +++ b/src/STUBS/mpi.c @@ -125,6 +125,15 @@ int MPI_Send(void *buf, int count, MPI_Datatype datatype, /* ---------------------------------------------------------------------- */ +int MPI_Isend(void *buf, int count, MPI_Datatype datatype, + int source, int tag, MPI_Comm comm, MPI_Request *request) +{ + printf("MPI Stub WARNING: Should not send message to self\n"); + return 0; +} + +/* ---------------------------------------------------------------------- */ + int MPI_Rsend(void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm) { diff --git a/src/STUBS/mpi.h b/src/STUBS/mpi.h index 7f2c5d931c..090fcb1554 100644 --- a/src/STUBS/mpi.h +++ b/src/STUBS/mpi.h @@ -73,6 +73,8 @@ int MPI_Type_size(int, int *); int MPI_Send(void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm); +int MPI_Isend(void *buf, int count, MPI_Datatype datatype, + int source, int tag, MPI_Comm comm, MPI_Request *request); int MPI_Rsend(void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm); int MPI_Recv(void *buf, int count, MPI_Datatype datatype, diff --git a/src/USER-MISC/dihedral_table.h b/src/USER-MISC/dihedral_table.h index ff6a1d307e..369a6b480b 100644 --- a/src/USER-MISC/dihedral_table.h +++ b/src/USER-MISC/dihedral_table.h @@ -36,15 +36,15 @@ namespace LAMMPS_NS { class DihedralTable : public Dihedral { public: DihedralTable(class LAMMPS *); - ~DihedralTable(); - void compute(int, int); + virtual ~DihedralTable(); + virtual void compute(int, int); void settings(int, char **); void coeff(int, char **); void write_restart(FILE *); void read_restart(FILE *); double single(int type, int i1, int i2, int i3, int i4); - private: + protected: int tabstyle,tablength; // double *phi0; <- equilibrium angles not supported diff --git a/src/USER-OMP/bond_fene_expand_omp.cpp b/src/USER-OMP/bond_fene_expand_omp.cpp index f77e11d1df..4e6b7d521a 100644 --- a/src/USER-OMP/bond_fene_expand_omp.cpp +++ b/src/USER-OMP/bond_fene_expand_omp.cpp @@ -81,6 +81,7 @@ void BondFENEExpandOMP::eval(int nfrom, int nto, ThrData * const thr) double * const * const f = thr->get_f(); const int * const * const bondlist = neighbor->bondlist; const int nlocal = atom->nlocal; + const int tid = thr->get_tid(); for (n = nfrom; n < nto; n++) { i1 = bondlist[n][0]; @@ -105,10 +106,14 @@ void BondFENEExpandOMP::eval(int nfrom, int nto, ThrData * const thr) if (rlogarg < 0.1) { char str[128]; + sprintf(str,"FENE bond too long: " BIGINT_FORMAT " %d %d %g", update->ntimestep,atom->tag[i1],atom->tag[i2],sqrt(rsq)); error->warning(FLERR,str,0); - if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond"); + + if (check_error_thr((rlogarg <= -3.0),tid,FLERR,"Bad FENE bond")) + return; + rlogarg = 0.1; } diff --git a/src/USER-OMP/bond_fene_omp.cpp b/src/USER-OMP/bond_fene_omp.cpp index 2758ed9e9e..65c8ff62cf 100644 --- a/src/USER-OMP/bond_fene_omp.cpp +++ b/src/USER-OMP/bond_fene_omp.cpp @@ -80,6 +80,7 @@ void BondFENEOMP::eval(int nfrom, int nto, ThrData * const thr) double * const * const f = thr->get_f(); const int * const * const bondlist = neighbor->bondlist; const int nlocal = atom->nlocal; + const int tid = thr->get_tid(); for (n = nfrom; n < nto; n++) { i1 = bondlist[n][0]; @@ -101,10 +102,14 @@ void BondFENEOMP::eval(int nfrom, int nto, ThrData * const thr) if (rlogarg < 0.1) { char str[128]; + sprintf(str,"FENE bond too long: " BIGINT_FORMAT " %d %d %g", update->ntimestep,atom->tag[i1],atom->tag[i2],sqrt(rsq)); error->warning(FLERR,str,0); - if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond"); + + if (check_error_thr((rlogarg <= -3.0),tid,FLERR,"Bad FENE bond")) + return; + rlogarg = 0.1; } diff --git a/src/USER-OMP/neigh_derive_omp.cpp b/src/USER-OMP/neigh_derive_omp.cpp index 3b8fd79a20..5f42dd2b50 100644 --- a/src/USER-OMP/neigh_derive_omp.cpp +++ b/src/USER-OMP/neigh_derive_omp.cpp @@ -86,8 +86,8 @@ void Neighbor::half_from_full_no_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = inum_full; @@ -175,8 +175,8 @@ void Neighbor::half_from_full_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = inum_full; diff --git a/src/USER-OMP/neigh_full_omp.cpp b/src/USER-OMP/neigh_full_omp.cpp index 5f3ebd4d96..dcfaad9a49 100644 --- a/src/USER-OMP/neigh_full_omp.cpp +++ b/src/USER-OMP/neigh_full_omp.cpp @@ -105,8 +105,8 @@ void Neighbor::full_nsq_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -213,8 +213,8 @@ void Neighbor::full_nsq_ghost_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -314,8 +314,8 @@ void Neighbor::full_bin_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -324,7 +324,7 @@ void Neighbor::full_bin_omp(NeighList *list) /* ---------------------------------------------------------------------- binned neighbor list construction for all neighbors - include neighbors of ghost atoms (no "special neighbors" for ghosts) + include neighbors of ghost atoms, but no "special neighbors" for ghosts every neighbor pair appears in list of both atoms i and j ------------------------------------------------------------------------- */ @@ -446,8 +446,8 @@ void Neighbor::full_bin_ghost_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -555,8 +555,8 @@ void Neighbor::full_multi_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; diff --git a/src/USER-OMP/neigh_gran_omp.cpp b/src/USER-OMP/neigh_gran_omp.cpp index eca58ddccd..8a03da8e64 100644 --- a/src/USER-OMP/neigh_gran_omp.cpp +++ b/src/USER-OMP/neigh_gran_omp.cpp @@ -163,8 +163,8 @@ void Neighbor::granular_nsq_no_newton_omp(NeighList *list) firstshear[i] = shearptr; } npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -266,8 +266,8 @@ void Neighbor::granular_nsq_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -427,8 +427,8 @@ void Neighbor::granular_bin_no_newton_omp(NeighList *list) firstshear[i] = shearptr; } npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -544,8 +544,8 @@ void Neighbor::granular_bin_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -649,8 +649,8 @@ void Neighbor::granular_bin_newton_tri_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; diff --git a/src/USER-OMP/neigh_half_bin_omp.cpp b/src/USER-OMP/neigh_half_bin_omp.cpp index 3df11d38ac..3ff15b11da 100644 --- a/src/USER-OMP/neigh_half_bin_omp.cpp +++ b/src/USER-OMP/neigh_half_bin_omp.cpp @@ -118,8 +118,8 @@ void Neighbor::half_bin_no_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -245,8 +245,8 @@ void Neighbor::half_bin_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -356,8 +356,8 @@ void Neighbor::half_bin_newton_tri_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; diff --git a/src/USER-OMP/neigh_half_multi_omp.cpp b/src/USER-OMP/neigh_half_multi_omp.cpp index dcc2661310..ca74276139 100644 --- a/src/USER-OMP/neigh_half_multi_omp.cpp +++ b/src/USER-OMP/neigh_half_multi_omp.cpp @@ -126,8 +126,8 @@ void Neighbor::half_multi_no_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -263,8 +263,8 @@ void Neighbor::half_multi_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -384,8 +384,8 @@ void Neighbor::half_multi_newton_tri_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; diff --git a/src/USER-OMP/neigh_half_nsq_omp.cpp b/src/USER-OMP/neigh_half_nsq_omp.cpp index 0d389b45f6..c557fe7d84 100644 --- a/src/USER-OMP/neigh_half_nsq_omp.cpp +++ b/src/USER-OMP/neigh_half_nsq_omp.cpp @@ -105,8 +105,8 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; @@ -214,8 +214,8 @@ void Neighbor::half_nsq_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } NEIGH_OMP_CLOSE; list->inum = nlocal; diff --git a/src/USER-OMP/neigh_respa_omp.cpp b/src/USER-OMP/neigh_respa_omp.cpp index f8487244a5..e0817385d6 100644 --- a/src/USER-OMP/neigh_respa_omp.cpp +++ b/src/USER-OMP/neigh_respa_omp.cpp @@ -52,7 +52,7 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,n,itype,jtype,which,n_inner,n_middle; + int i,j,n,itype,jtype,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; @@ -91,6 +91,8 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list) int npage_middle = tid; int npnt_middle = 0; + int which = 0; + for (i = ifrom; i < ito; i++) { #if defined(_OPENMP) @@ -168,23 +170,23 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); ilist_inner[i] = i; firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; - if (npnt_inner >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_inner > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); if (respamiddle) { ilist_middle[i] = i; firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; - if (npnt_middle >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_middle > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } } NEIGH_OMP_CLOSE; @@ -225,7 +227,7 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,n,itype,jtype,itag,jtag,which,n_inner,n_middle; + int i,j,n,itype,jtype,itag,jtag,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; @@ -264,6 +266,8 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list) int npage_middle = tid; int npnt_middle = 0; + int which = 0; + for (i = ifrom; i < ito; i++) { #if defined(_OPENMP) @@ -359,23 +363,23 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); ilist_inner[i] = i; firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; - if (npnt_inner >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_inner > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); if (respamiddle) { ilist_middle[i] = i; firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; - if (npnt_middle >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_middle > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } } NEIGH_OMP_CLOSE; @@ -419,7 +423,7 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle; + int i,j,k,n,itype,jtype,ibin,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; @@ -459,6 +463,8 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list) int npage_middle = tid; int npnt_middle = 0; + int which = 0; + for (i = ifrom; i < ito; i++) { #if defined(_OPENMP) @@ -545,23 +551,23 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); ilist_inner[i] = i; firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; - if (npnt_inner >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_inner > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); if (respamiddle) { ilist_middle[i] = i; firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; - if (npnt_middle >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_middle > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } } NEIGH_OMP_CLOSE; @@ -604,7 +610,7 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle; + int i,j,k,n,itype,jtype,ibin,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; @@ -644,6 +650,8 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) int npage_middle = tid; int npnt_middle = 0; + int which = 0; + for (i = ifrom; i < ito; i++) { #if defined(_OPENMP) @@ -766,23 +774,23 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); ilist_inner[i] = i; firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; - if (npnt_inner >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_inner > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); if (respamiddle) { ilist_middle[i] = i; firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; - if (npnt_middle >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_middle > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } } NEIGH_OMP_CLOSE; @@ -825,7 +833,7 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which,n_inner,n_middle; + int i,j,k,n,itype,jtype,ibin,n_inner,n_middle; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; @@ -865,6 +873,8 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list) int npage_middle = tid; int npnt_middle = 0; + int which = 0; + for (i = ifrom; i < ito; i++) { #if defined(_OPENMP) @@ -959,23 +969,23 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list) firstneigh[i] = neighptr; numneigh[i] = n; npnt += n; - if (n > oneatom || npnt >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); ilist_inner[i] = i; firstneigh_inner[i] = neighptr_inner; numneigh_inner[i] = n_inner; npnt_inner += n_inner; - if (npnt_inner >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_inner > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); if (respamiddle) { ilist_middle[i] = i; firstneigh_middle[i] = neighptr_middle; numneigh_middle[i] = n_middle; npnt_middle += n_middle; - if (npnt_middle >= pgsize) - error->one(FLERR,"Neighbor list overflow, boost neigh_modify one or page"); + if (n_middle > oneatom) + error->one(FLERR,"Neighbor list overflow, boost neigh_modify one"); } } NEIGH_OMP_CLOSE; diff --git a/src/USER-OMP/neighbor_omp.h b/src/USER-OMP/neighbor_omp.h index 8cdd61d5b1..138ef96d5b 100644 --- a/src/USER-OMP/neighbor_omp.h +++ b/src/USER-OMP/neighbor_omp.h @@ -1,4 +1,4 @@ -/* ---------------------------------------------------------------------- +/* -*- c++ -*- ---------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov @@ -35,9 +35,8 @@ namespace LAMMPS_NS { const int tid = omp_get_thread_num(); \ const int idelta = 1 + num/nthreads; \ const int ifrom = tid*idelta; \ - int ito = ifrom + idelta; \ - if (ito > num) \ - ito = num + const int ito = ((ifrom + idelta) > num) \ + ? num : (ifrom+idelta); \ #define NEIGH_OMP_CLOSE } diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp index bb8b38a5b1..0ff25d16bd 100644 --- a/src/USER-OMP/pair_airebo_omp.cpp +++ b/src/USER-OMP/pair_airebo_omp.cpp @@ -2653,7 +2653,6 @@ void PairAIREBOOMP::TORSION_thr(int ifrom, int ito, void PairAIREBOOMP::REBO_neigh_thr() { const int nlocal = atom->nlocal; - const int nall = nlocal + atom->nghost; const int nthreads = comm->nthreads; if (atom->nmax > maxlocal) { @@ -2697,9 +2696,7 @@ void PairAIREBOOMP::REBO_neigh_thr() const int iidelta = 1 + allnum/nthreads; const int iifrom = tid*iidelta; - int iito = iifrom + iidelta; - if (iito > allnum) - iito = allnum; + const int iito = ((iifrom+iidelta)>allnum) ? allnum : (iifrom+iidelta); // store all REBO neighs of owned and ghost atoms // scan full neighbor list of I @@ -2757,6 +2754,7 @@ void PairAIREBOOMP::REBO_neigh_thr() } } + /* ---------------------------------------------------------------------- */ double PairAIREBOOMP::memory_usage() diff --git a/src/USER-OMP/pair_colloid_omp.cpp b/src/USER-OMP/pair_colloid_omp.cpp index 7bfe1c04de..5cb47c664c 100644 --- a/src/USER-OMP/pair_colloid_omp.cpp +++ b/src/USER-OMP/pair_colloid_omp.cpp @@ -85,6 +85,7 @@ void PairColloidOMP::eval(int iifrom, int iito, ThrData * const thr) double * const * const f = thr->get_f(); const int * const type = atom->type; const int nlocal = atom->nlocal; + const int tid = thr->get_tid(); const double * const special_lj = force->special_lj; double fxtmp,fytmp,fztmp; @@ -146,7 +147,11 @@ void PairColloidOMP::eval(int iifrom, int iito, ThrData * const thr) evdwl = 2.0/9.0*fR * (1.0-(K[1]*(K[1]*(K[1]/3.0+3.0*K[2])+4.2*K[4])+K[2]*K[4]) * sigma6[itype][jtype]/K[6]) - offset[itype][jtype]; - if (rsq <= K[1]) error->one(FLERR,"Overlapping small/large in pair colloid"); + + if (check_error_thr((rsq <= K[1]),tid,FLERR, + "Overlapping small/large in pair colloid")) + return; + break; case LARGE_LARGE: diff --git a/src/USER-OMP/pair_comb_omp.cpp b/src/USER-OMP/pair_comb_omp.cpp index 80607b7225..797b9cc8fd 100644 --- a/src/USER-OMP/pair_comb_omp.cpp +++ b/src/USER-OMP/pair_comb_omp.cpp @@ -22,12 +22,10 @@ #include "neighbor.h" #include "neigh_list.h" -#if defined(_OPENMP) -#include -#endif - using namespace LAMMPS_NS; +#define MAXNEIGH 24 + /* ---------------------------------------------------------------------- */ PairCombOMP::PairCombOMP(LAMMPS *lmp) : @@ -48,17 +46,8 @@ void PairCombOMP::compute(int eflag, int vflag) const int nthreads = comm->nthreads; const int inum = list->inum; - // grow coordination array if necessary - - if (atom->nmax > nmax) { - memory->destroy(NCo); - memory->destroy(bbij); - nmax = atom->nmax; - memory->create(NCo,nmax,"pair:NCo"); - memory->create(bbij,nmax,nmax,"pair:bbij"); - } - // Build short range neighbor list + Short_neigh_thr(); #if defined(_OPENMP) @@ -102,7 +91,7 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) double yaself; double potal,fac11,fac11e; double vionij,fvionij,sr1,sr2,sr3,Eov,Fov; - int sht_jnum, *sht_jlist; + int sht_jnum, *sht_jlist, nj; evdwl = ecoul = 0.0; @@ -139,7 +128,8 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) fxtmp = fytmp = fztmp = 0.0; iq = q[i]; - NCo[i] = 0; + NCo[i] = 0; + nj = 0; iparam_i = elem2param[itype][itype][itype]; // self energy, only on i atom @@ -242,7 +232,6 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) int numcoor = 0; for (jj = 0; jj < sht_jnum; jj++) { j = sht_jlist[jj]; - j &= NEIGHMASK; jtype = map[type[j]]; iparam_ij = elem2param[itype][jtype][jtype]; @@ -264,7 +253,6 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) for (jj = 0; jj < sht_jnum; jj++) { j = sht_jlist[jj]; - j &= NEIGHMASK; jtype = map[type[j]]; iparam_ij = elem2param[itype][jtype][jtype]; @@ -279,6 +267,7 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) rsq1 = vec3_dot(delr1,delr1); if (rsq1 > params[iparam_ij].cutsq) continue; + nj ++; // accumulate bondorder zeta for each i-j interaction via loop over k @@ -289,7 +278,6 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) for (kk = 0; kk < sht_jnum; kk++) { k = sht_jlist[kk]; if (j == k) continue; - k &= NEIGHMASK; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; @@ -309,7 +297,7 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) if (cuo_flag1 && cuo_flag2) cuo_flag = 1; else cuo_flag = 0; - force_zeta(¶ms[iparam_ij],EFLAG,i,j,rsq1,zeta_ij, + force_zeta(¶ms[iparam_ij],EFLAG,i,nj,rsq1,zeta_ij, iq,jq,fpair,prefactor,evdwl); // over-coordination correction for HfO2 @@ -334,7 +322,6 @@ void PairCombOMP::eval(int iifrom, int iito, ThrData * const thr) for (kk = 0; kk < sht_jnum; kk++) { k = sht_jlist[kk]; if (j == k) continue; - k &= NEIGHMASK; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; @@ -429,6 +416,7 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup) const int i = ilist[ii]; const int itag = tag[i]; + int nj = 0; if (mask[i] & groupbit) { fqi = fqj = fqij = fqji = fqjj = 0.0; // should not be needed. @@ -511,10 +499,11 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup) const int iparam_ij = elem2param[itype][jtype][jtype]; if (rsq1 > params[iparam_ij].cutsq) continue; + nj ++; // charge force in Aij and Bij - qfo_short(¶ms[iparam_ij],i,j,rsq1,iq,jq,fqij,fqjj); + qfo_short(¶ms[iparam_ij],i,nj,rsq1,iq,jq,fqij,fqjj); fqi += fqij; #if defined(_OPENMP) #pragma omp atomic @@ -522,7 +511,7 @@ double PairCombOMP::yasu_char(double *qf_fix, int &igroup) qf[j] += fqjj; } -#if defined(_OPENMP) && 0 +#if defined(_OPENMP) #pragma omp atomic #endif qf[i] += fqi; @@ -559,12 +548,13 @@ void PairCombOMP::Short_neigh_thr() { if (atom->nmax > nmax) { - nmax = int(1.0 * atom->nmax); - memory->sfree(sht_num); + nmax = atom->nmax; memory->sfree(sht_first); - memory->create(sht_num,nmax,"pair:sht_num"); sht_first = (int **) memory->smalloc(nmax*sizeof(int *), "pair:sht_first"); + memory->grow(sht_num,nmax,"pair:sht_num"); + memory->grow(NCo,nmax,"pair:NCo"); + memory->grow(bbij,nmax,nmax,"pair:bbij"); } const int nthreads = comm->nthreads; @@ -602,46 +592,47 @@ void PairCombOMP::Short_neigh_thr() int npage = tid; npntj = 0; - for (ii = iifrom; ii < iito; ii++) { - i = ilist[ii]; - itype = type[i]; + if (iifrom < inum) { + for (ii = iifrom; ii < iito; ii++) { + i = ilist[ii]; + itype = type[i]; #if defined(_OPENMP) #pragma omp critical #endif - if(pgsize - npntj < oneatom) { - npntj = 0; - npage++; - if (npage == maxpage) add_pages(nthreads); - } + if(pgsize - npntj < oneatom) { + npntj = 0; + npage += nthreads; + if (npage >= maxpage) add_pages(nthreads); + } - neighptrj = &pages[npage][npntj]; - nj = 0; + neighptrj = &pages[npage][npntj]; + nj = 0; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; - jlist = firstneigh[i]; - jnum = numneigh[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - jtype = type[j]; - iparam_ij = elem2param[itype][jtype][jtype]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + jtype = type[j]; - delrj[0] = xtmp - x[j][0]; - delrj[1] = ytmp - x[j][1]; - delrj[2] = ztmp - x[j][2]; - rsq = vec3_dot(delrj,delrj); + delrj[0] = xtmp - x[j][0]; + delrj[1] = ytmp - x[j][1]; + delrj[2] = ztmp - x[j][2]; + rsq = vec3_dot(delrj,delrj); - if (rsq > cutmin) continue; - neighptrj[nj++] = j; + if (rsq > cutmin) continue; + neighptrj[nj++] = j; + } + sht_first[i] = neighptrj; + sht_num[i] = nj; + npntj += nj; } - sht_first[i] = neighptrj; - sht_num[i] = nj; - npntj += nj; } } } diff --git a/src/USER-OMP/pair_lj_cut_coul_long_tip4p_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_long_tip4p_omp.cpp index 78f35709a2..c647fbb067 100644 --- a/src/USER-OMP/pair_lj_cut_coul_long_tip4p_omp.cpp +++ b/src/USER-OMP/pair_lj_cut_coul_long_tip4p_omp.cpp @@ -151,6 +151,7 @@ void PairLJCutCoulLongTIP4POMP::eval(int iifrom, int iito, ThrData * const thr) const double * const q = atom->q; const int * const type = atom->type; const int nlocal = atom->nlocal; + const int tid = thr->get_tid(); const double * const special_coul = force->special_coul; const double * const special_lj = force->special_lj; const double qqrd2e = force->qqrd2e; @@ -221,8 +222,11 @@ void PairLJCutCoulLongTIP4POMP::eval(int iifrom, int iito, ThrData * const thr) x2 = mpos[j]; jH1 = h1idx[j]; jH2 = h2idx[j]; - if (jtype == typeO && ( jH1 < 0 || jH2 < 0 )) - error->one(FLERR,"TIP4P hydrogen is missing"); + + if (check_error_thr((jtype == typeO && ( jH1 < 0 || jH2 < 0 )), + tid, FLERR,"TIP4P hydrogen is missing")) + return; + delx = x1[0] - x2[0]; dely = x1[1] - x2[1]; delz = x1[2] - x2[2]; diff --git a/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp b/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp index a4aa811075..ec58456887 100644 --- a/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp +++ b/src/USER-OMP/pair_lj_cut_coul_pppm_tip4p_omp.cpp @@ -179,6 +179,7 @@ void PairLJCutCoulPPPMTIP4POMP::eval(int iifrom, int iito, ThrData * const thr) const double * const q = atom->q; const int * const type = atom->type; const int nlocal = atom->nlocal; + const int tid = thr->get_tid(); const double * const special_coul = force->special_coul; const double * const special_lj = force->special_lj; const double qqrd2e = force->qqrd2e; @@ -249,8 +250,11 @@ void PairLJCutCoulPPPMTIP4POMP::eval(int iifrom, int iito, ThrData * const thr) x2 = mpos[j]; jH1 = h1idx[j]; jH2 = h2idx[j]; - if (jtype == typeO && ( jH1 < 0 || jH2 < 0 )) - error->one(FLERR,"TIP4P hydrogen is missing"); + + if (check_error_thr((jtype == typeO && ( jH1 < 0 || jH2 < 0 )), + tid, FLERR,"TIP4P hydrogen is missing")) + return; + delx = x1[0] - x2[0]; dely = x1[1] - x2[1]; delz = x1[2] - x2[2]; diff --git a/src/USER-OMP/pair_peri_lps_omp.cpp b/src/USER-OMP/pair_peri_lps_omp.cpp index a4479180bf..a23f202455 100644 --- a/src/USER-OMP/pair_peri_lps_omp.cpp +++ b/src/USER-OMP/pair_peri_lps_omp.cpp @@ -213,16 +213,15 @@ void PairPeriLPSOMP::eval(int iifrom, int iito, ThrData * const thr) // each thread works on a fixed chunk of atoms. const int idelta = 1 + nlocal/comm->nthreads; iifrom = thr->get_tid()*idelta; - iito = iifrom + idelta; - if (iito > nlocal) - iito = nlocal; + iito = ((iifrom + idelta) > nlocal) ? nlocal : (iifrom + idelta); #else iifrom = 0; iito = nlocal; #endif // Compute the dilatation on each particle - compute_dilatation_thr(iifrom, iito); + if (iifrom < nlocal) + compute_dilatation_thr(iifrom, iito); // wait until all threads are done before communication sync_threads(); @@ -350,7 +349,8 @@ void PairPeriLPSOMP::eval(int iifrom, int iito, ThrData * const thr) sync_threads(); // store new s0 (in parallel) - for (i = iifrom; i < iito; i++) s0[i] = s0_new[i]; + if (iifrom < nlocal) + for (i = iifrom; i < iito; i++) s0[i] = s0_new[i]; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/pair_peri_pmb_omp.cpp b/src/USER-OMP/pair_peri_pmb_omp.cpp index de6387ae6b..2173bf9e8f 100644 --- a/src/USER-OMP/pair_peri_pmb_omp.cpp +++ b/src/USER-OMP/pair_peri_pmb_omp.cpp @@ -203,9 +203,7 @@ void PairPeriPMBOMP::eval(int iifrom, int iito, ThrData * const thr) // each thread works on a fixed chunk of atoms. const int idelta = 1 + nlocal/comm->nthreads; iifrom = thr->get_tid()*idelta; - iito = iifrom + idelta; - if (iito > nlocal) - iito = nlocal; + iito = ((iifrom + idelta) > nlocal) ? nlocal : (iifrom + idelta); #else iifrom = 0; iito = nlocal; @@ -234,8 +232,8 @@ void PairPeriPMBOMP::eval(int iifrom, int iito, ThrData * const thr) // check if lost a partner without first breaking bond if (j < 0) { - partner[i][jj] = 0; - continue; + partner[i][jj] = 0; + continue; } // compute force density, add to PD equation of motion @@ -257,7 +255,7 @@ void PairPeriPMBOMP::eval(int iifrom, int iito, ThrData * const thr) // scale vfrac[j] if particle j near the horizon if ((fabs(r0[i][jj] - delta)) <= half_lc) - vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); else vfrac_scale = 1.0; @@ -285,9 +283,9 @@ void PairPeriPMBOMP::eval(int iifrom, int iito, ThrData * const thr) // update s0 for next timestep if (first) - s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch); + s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch); else - s0_new[i] = MAX(s0_new[i],s00[itype][jtype] - (alpha[itype][jtype] * stretch)); + s0_new[i] = MAX(s0_new[i],s00[itype][jtype] - (alpha[itype][jtype] * stretch)); first = false; } @@ -296,7 +294,8 @@ void PairPeriPMBOMP::eval(int iifrom, int iito, ThrData * const thr) sync_threads(); // store new s0 (in parallel) - for (i = iifrom; i < iito; i++) s0[i] = s0_new[i]; + if (iifrom < nlocal) + for (i = iifrom; i < iito; i++) s0[i] = s0_new[i]; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/pair_table_omp.cpp b/src/USER-OMP/pair_table_omp.cpp index e8d63e590d..a20538f2c0 100644 --- a/src/USER-OMP/pair_table_omp.cpp +++ b/src/USER-OMP/pair_table_omp.cpp @@ -88,6 +88,7 @@ void PairTableOMP::eval(int iifrom, int iito, ThrData * const thr) double * const * const f = thr->get_f(); const int * const type = atom->type; const int nlocal = atom->nlocal; + const int tid = thr->get_tid(); const double * const special_lj = force->special_lj; double fxtmp,fytmp,fztmp; @@ -121,25 +122,36 @@ void PairTableOMP::eval(int iifrom, int iito, ThrData * const thr) if (rsq < cutsq[itype][jtype]) { tb = &tables[tabindex[itype][jtype]]; - if (rsq < tb->innersq) - error->one(FLERR,"Pair distance < table inner cutoff"); + + if (check_error_thr((rsq < tb->innersq),tid, + FLERR,"Pair distance < table inner cutoff")) + return; if (tabstyle == LOOKUP) { itable = static_cast ((rsq - tb->innersq) * tb->invdelta); - if (itable >= tlm1) - error->one(FLERR,"Pair distance > table outer cutoff"); + + if (check_error_thr((itable >= tlm1),tid, + FLERR,"Pair distance > table outer cutoff")) + return; + fpair = factor_lj * tb->f[itable]; } else if (tabstyle == LINEAR) { itable = static_cast ((rsq - tb->innersq) * tb->invdelta); - if (itable >= tlm1) - error->one(FLERR,"Pair distance > table outer cutoff"); + + if (check_error_thr((itable >= tlm1),tid, + FLERR,"Pair distance > table outer cutoff")) + return; + fraction = (rsq - tb->rsq[itable]) * tb->invdelta; value = tb->f[itable] + fraction*tb->df[itable]; fpair = factor_lj * value; } else if (tabstyle == SPLINE) { itable = static_cast ((rsq - tb->innersq) * tb->invdelta); - if (itable >= tlm1) - error->one(FLERR,"Pair distance > table outer cutoff"); + + if (check_error_thr((itable >= tlm1),tid, + FLERR,"Pair distance > table outer cutoff")) + return; + b = (rsq - tb->rsq[itable]) * tb->invdelta; a = 1.0 - b; value = a * tb->f[itable] + b * tb->f[itable+1] + diff --git a/src/USER-OMP/pair_tersoff_table_omp.cpp b/src/USER-OMP/pair_tersoff_table_omp.cpp index 6c9a680f09..20c47ed598 100644 --- a/src/USER-OMP/pair_tersoff_table_omp.cpp +++ b/src/USER-OMP/pair_tersoff_table_omp.cpp @@ -16,6 +16,7 @@ #include "pair_tersoff_table_omp.h" #include "atom.h" #include "comm.h" +#include "error.h" #include "force.h" #include "neighbor.h" #include "neigh_list.h" @@ -89,6 +90,7 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr) double * const * const f = thr->get_f(); const int * const type = atom->type; const int nlocal = atom->nlocal; + const int tid = thr->get_tid(); inum = list->inum; ilist = list->ilist; @@ -108,6 +110,12 @@ void PairTersoffTableOMP::eval(int iifrom, int iito, ThrData * const thr) jlist = firstneigh[i]; jnum = numneigh[i]; + if (check_error_thr((jnum > leadingDimensionInteractionList), tid, + FLERR,"Too many neighbors for interaction list.\n" + "Check your system or increase 'leadingDimension" + "InteractionList'")) + return; + // Pre-calculate gteta and cutoff function for (int neighbor_j = 0; neighbor_j < jnum; neighbor_j++) { diff --git a/src/USER-OMP/thr_omp.cpp b/src/USER-OMP/thr_omp.cpp index 5ae7cc87da..8ff0c9b6cf 100644 --- a/src/USER-OMP/thr_omp.cpp +++ b/src/USER-OMP/thr_omp.cpp @@ -42,7 +42,8 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -ThrOMP::ThrOMP(LAMMPS *ptr, int style) : lmp(ptr), fix(NULL), thr_style(style) +ThrOMP::ThrOMP(LAMMPS *ptr, int style) : lmp(ptr), fix(NULL), + thr_style(style), thr_error(0) { // register fix omp with this class int ifix = lmp->modify->find_fix("package_omp"); @@ -66,6 +67,7 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom, double **vatom, ThrData *thr) { const int tid = thr->get_tid(); + if (tid == 0) thr_error = 0; if (thr_style & THR_PAIR) { if (eflag & 2) { @@ -171,6 +173,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, double **f = lmp->atom->f; double **x = lmp->atom->x; + if (evflag) + sync_threads(); + switch (thr_style) { case THR_PAIR: { @@ -201,11 +206,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } } if (eflag & 2) { - sync_threads(); data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid); } if (vflag & 4) { - sync_threads(); data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid); } } @@ -248,11 +251,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } } if (eflag & 2) { - sync_threads(); data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid); } if (vflag & 4) { - sync_threads(); data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid); } } @@ -272,11 +273,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, bond->virial[i] += thr->virial_bond[i]; } if (eflag & 2) { - sync_threads(); data_reduce_thr(&(bond->eatom[0]), nall, nthreads, 1, tid); } if (vflag & 4) { - sync_threads(); data_reduce_thr(&(bond->vatom[0][0]), nall, nthreads, 6, tid); } } @@ -295,11 +294,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, angle->virial[i] += thr->virial_angle[i]; } if (eflag & 2) { - sync_threads(); data_reduce_thr(&(angle->eatom[0]), nall, nthreads, 1, tid); } if (vflag & 4) { - sync_threads(); data_reduce_thr(&(angle->vatom[0][0]), nall, nthreads, 6, tid); } } @@ -318,11 +315,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, dihedral->virial[i] += thr->virial_dihed[i]; } if (eflag & 2) { - sync_threads(); data_reduce_thr(&(dihedral->eatom[0]), nall, nthreads, 1, tid); } if (vflag & 4) { - sync_threads(); data_reduce_thr(&(dihedral->vatom[0][0]), nall, nthreads, 6, tid); } } @@ -351,12 +346,10 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } } if (eflag & 2) { - sync_threads(); data_reduce_thr(&(dihedral->eatom[0]), nall, nthreads, 1, tid); data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid); } if (vflag & 4) { - sync_threads(); data_reduce_thr(&(dihedral->vatom[0][0]), nall, nthreads, 6, tid); data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid); } @@ -376,11 +369,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, improper->virial[i] += thr->virial_imprp[i]; } if (eflag & 2) { - sync_threads(); data_reduce_thr(&(improper->eatom[0]), nall, nthreads, 1, tid); } if (vflag & 4) { - sync_threads(); data_reduce_thr(&(improper->vatom[0][0]), nall, nthreads, 6, tid); } } @@ -401,11 +392,9 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, kspace->virial[i] += thr->virial_kspce[i]; } if (eflag & 2) { - sync_threads(); data_reduce_thr(&(kspace->eatom[0]), nall, nthreads, 1, tid); } if (vflag & 4) { - sync_threads(); data_reduce_thr(&(kspace->vatom[0][0]), nall, nthreads, 6, tid); } } @@ -418,7 +407,6 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } if (style == fix->last_omp_style) { - sync_threads(); data_reduce_thr(&(f[0][0]), nall, nthreads, 3, tid); if (lmp->atom->torque) data_reduce_thr(&(lmp->atom->torque[0][0]), nall, nthreads, 3, tid); diff --git a/src/USER-OMP/thr_omp.h b/src/USER-OMP/thr_omp.h index 32f7045124..72fc330cb4 100644 --- a/src/USER-OMP/thr_omp.h +++ b/src/USER-OMP/thr_omp.h @@ -19,6 +19,7 @@ #define LMP_THR_OMP_H #include "pointers.h" +#include "error.h" #include "fix_omp.h" #include "thr_data.h" @@ -40,6 +41,7 @@ class ThrOMP { FixOMP *fix; // pointer to fix_omp; const int thr_style; + int thr_error; public: ThrOMP(LAMMPS *, int); @@ -68,7 +70,34 @@ class ThrOMP { const int, const int, const int); // reduce per thread data as needed - void reduce_thr(void * const style, const int eflag, const int vflag, ThrData * const thr, const int nproxy=0); + void reduce_thr(void * const style, const int eflag, const int vflag, + ThrData * const thr, const int nproxy=0); + + // thread safe variant error abort support. + // signals an error condition in any thread by making + // thr_error > 0, if condition "cond" is true. + // will abort from thread 0 if thr_error is > 0 + // otherwise return true. + // returns false if no error found on any thread. + // use return value to jump/return to end of threaded region. + + bool check_error_thr(const bool cond, const int tid, const char *fname, + const int line, const char *errmsg) { + if (cond) { +#if defined(_OPENMP) +#pragma omp atomic + ++thr_error; +#endif + if (tid > 0) return true; + else lmp->error->one(fname,line,errmsg); + } else { + if (thr_error > 0) { + if (tid == 0) lmp->error->one(fname,line,errmsg); + else return true; + } else return false; + } + return false; + }; protected: diff --git a/src/memory.cpp b/src/memory.cpp index c98f40fd3f..4819303df2 100644 --- a/src/memory.cpp +++ b/src/memory.cpp @@ -34,7 +34,8 @@ void *Memory::smalloc(bigint nbytes, const char *name) #if defined(LAMMPS_MEMALIGN) void *ptr; - posix_memalign(&ptr, LAMMPS_MEMALIGN, nbytes); + int retval = posix_memalign(&ptr, LAMMPS_MEMALIGN, nbytes); + if (retval) ptr = NULL; #else void *ptr = malloc(nbytes); #endif diff --git a/tools/restart2data.cpp b/tools/restart2data.cpp index 6882251bcc..3654605add 100644 --- a/tools/restart2data.cpp +++ b/tools/restart2data.cpp @@ -475,7 +475,7 @@ int main (int narg, char **arg) void header(FILE *fp, Data &data) { - const char *version = "5 Jan 2012"; + const char *version = "31 Jan 2012"; data.triclinic = 0; @@ -1506,6 +1506,8 @@ void pair(FILE *fp, Data &data, char *style, int flag) } } + } else if (strcmp(style,"comb") == 0) { + } else if (strcmp(style,"coul/diel") == 0) { m = 1; double cut_coul = read_double(fp); @@ -2534,6 +2536,16 @@ void angle(FILE *fp, Data &data) } else if ((strcmp(data.angle_style,"cg/cmm") == 0) || (strcmp(data.angle_style,"sdk") == 0)) { + data.angle_harmonic_k = new double[data.nangletypes+1]; + data.angle_harmonic_theta0 = new double[data.nangletypes+1]; + data.angle_cg_cmm_epsilon = new double[data.nangletypes+1]; + + fread(&data.angle_harmonic_k[1],sizeof(double),data.nangletypes,fp); + fread(&data.angle_harmonic_theta0[1],sizeof(double),data.nangletypes,fp); + fread(&data.angle_cg_cmm_epsilon[1],sizeof(double),data.nangletypes,fp); + + } else if (strcmp(data.angle_style,"cg/cmm/old") == 0) { + data.angle_harmonic_k = new double[data.nangletypes+1]; data.angle_harmonic_theta0 = new double[data.nangletypes+1]; data.angle_cg_cmm_epsilon = new double[data.nangletypes+1]; @@ -3046,6 +3058,7 @@ void Data::write(FILE *fp, FILE *fp2) (strcmp(pair_style,"adp") != 0) && (strcmp(pair_style,"airebo") != 0) && (strcmp(pair_style,"brownian") != 0) && + (strcmp(pair_style,"comb") != 0) && (strcmp(pair_style,"coul/cut") != 0) && (strcmp(pair_style,"coul/debye") != 0) && (strcmp(pair_style,"coul/diel") != 0) && @@ -3181,7 +3194,7 @@ void Data::write(FILE *fp, FILE *fp2) (strcmp(pair_style,"cg/cmm/coul/long") == 0) || (strcmp(pair_style,"lj/sdk") == 0) || (strcmp(pair_style,"lj/sdk/coul/long") == 0)) { - printf("ERROR: Cannot write pair_style %s to data file\n", + printf("ERROR: Cannot write pair_style %s to data file alone. please provide an input file, too.\n", pair_style); exit(1); } @@ -3351,7 +3364,13 @@ void Data::write(FILE *fp, FILE *fp2) fprintf(fp,"%d %g %g\n",i, angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0); - } else if (strcmp(angle_style,"cg/cmm") == 0) { + } else if ((strcmp(angle_style,"cg/cmm") == 0) || + (strcmp(angle_style,"sdk") == 0)) { + for (int i = 1; i <= nangletypes; i++) + fprintf(fp,"%d %g %g %g\n",i,angle_harmonic_k[i], + angle_harmonic_theta0[i]/PI*180.0,angle_cg_cmm_epsilon[i]); + + } else if (strcmp(angle_style,"cg/cmm/old") == 0) { for (int i = 1; i <= nangletypes; i++) fprintf(fp,"%d %g %g %s %g %g\n",i, angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0, @@ -3376,7 +3395,13 @@ void Data::write(FILE *fp, FILE *fp2) fprintf(fp2,"angle_coeffs %d %g %g\n",i, angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0); - } else if (strcmp(angle_style,"cg/cmm") == 0) { + } else if ((strcmp(angle_style,"cg/cmm") == 0) || + (strcmp(angle_style,"sdk") == 0)) { + for (int i = 1; i <= nangletypes; i++) + fprintf(fp2,"angle_coeffs %d %g %g %g\n",i,angle_harmonic_k[i], + angle_harmonic_theta0[i]/PI*180.0,angle_cg_cmm_epsilon[i]); + + } else if (strcmp(angle_style,"cg/cmm/old") == 0) { for (int i = 1; i <= nangletypes; i++) fprintf(fp2,"angle_coeffs %d %g %g %s %g %g\n",i, angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0,