diff --git a/src/KSPACE/pppm.cpp b/src/KSPACE/pppm.cpp index 34b17533ea..89054a7685 100644 --- a/src/KSPACE/pppm.cpp +++ b/src/KSPACE/pppm.cpp @@ -22,6 +22,7 @@ #include "stdlib.h" #include "math.h" #include "pppm.h" +#include "math_const.h" #include "atom.h" #include "comm.h" #include "neighbor.h" @@ -35,8 +36,6 @@ #include "memory.h" #include "error.h" -#include "math_const.h" - using namespace LAMMPS_NS; using namespace MathConst; @@ -1050,7 +1049,8 @@ double PPPM::rms(double h, double prd, bigint natoms, compute difference in real-space and kspace RMS precision ------------------------------------------------------------------------- */ -double PPPM::diffpr(double h_x, double h_y, double h_z, double q2, double **acons) +double PPPM::diffpr(double h_x, double h_y, double h_z, double q2, + double **acons) { double lprx,lpry,lprz,kspace_prec,real_prec; double xprd = domain->xprd; @@ -1738,6 +1738,7 @@ void PPPM::fieldforce() } // convert E-field to force + const double qfactor = force->qqrd2e * scale * q[i]; f[i][0] += qfactor*ekx; f[i][1] += qfactor*eky; diff --git a/src/KSPACE/pppm_tip4p.cpp b/src/KSPACE/pppm_tip4p.cpp index 16a0f7a04c..bd6e880a3b 100644 --- a/src/KSPACE/pppm_tip4p.cpp +++ b/src/KSPACE/pppm_tip4p.cpp @@ -40,6 +40,18 @@ using namespace LAMMPS_NS; PPPMTIP4P::PPPMTIP4P(LAMMPS *lmp, int narg, char **arg) : PPPM(lmp, narg, arg) {} +/* ---------------------------------------------------------------------- */ + +void PPPMTIP4P::init() +{ + // TIP4P PPPM requires newton on, b/c it computes forces on ghost atoms + + if (force->newton == 0) + error->all(FLERR,"Kspace style pppm/tip4p requires newton on"); + + PPPM::init(); +} + /* ---------------------------------------------------------------------- find center grid pt for each of my particles check that full stencil for the particle will fit in my 3d brick @@ -206,13 +218,14 @@ void PPPMTIP4P::fieldforce() } // convert E-field to force + const double qfactor = force->qqrd2e * scale * q[i]; if (type[i] != typeO) { f[i][0] += qfactor*ekx; f[i][1] += qfactor*eky; f[i][2] += qfactor*ekz; - } else { + } else { fx = qfactor * ekx; fy = qfactor * eky; fz = qfactor * ekz; diff --git a/src/KSPACE/pppm_tip4p.h b/src/KSPACE/pppm_tip4p.h index 764db68edd..35dd747ad4 100644 --- a/src/KSPACE/pppm_tip4p.h +++ b/src/KSPACE/pppm_tip4p.h @@ -28,6 +28,7 @@ class PPPMTIP4P : public PPPM { public: PPPMTIP4P(class LAMMPS *, int, char **); virtual ~PPPMTIP4P () {}; + void init(); protected: virtual void particle_map(); diff --git a/src/MANYBODY/pair_airebo.h b/src/MANYBODY/pair_airebo.h index 5ef5d0ee14..2f78fc4511 100644 --- a/src/MANYBODY/pair_airebo.h +++ b/src/MANYBODY/pair_airebo.h @@ -163,6 +163,7 @@ class PairAIREBO : public Pair { }; /* kronecker delta function returning a double */ + inline double kronecker(const int a, const int b) const { return (a == b) ? 1.0 : 0.0; }; diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp index c8747aa811..f6719fa3a5 100644 --- a/src/REPLICA/neb.cpp +++ b/src/REPLICA/neb.cpp @@ -215,6 +215,7 @@ void NEB::run() // break induced if converged // damped dynamic min styles insure all replicas converge together + timer->init(); timer->barrier_start(TIME_LOOP); while (update->minimize->niter < n1steps) { @@ -284,7 +285,9 @@ void NEB::run() // break induced if converged // damped dynamic min styles insure all replicas converge together + timer->init(); timer->barrier_start(TIME_LOOP); + while (update->minimize->niter < n2steps) { update->minimize->run(nevery); print_status(); diff --git a/src/REPLICA/prd.cpp b/src/REPLICA/prd.cpp index 909d7e35ec..ff20afe735 100644 --- a/src/REPLICA/prd.cpp +++ b/src/REPLICA/prd.cpp @@ -261,8 +261,11 @@ void PRD::command(int narg, char **arg) quench(); ncoincident = 0; share_event(0,0); + + timer->init(); timer->barrier_start(TIME_LOOP); time_start = timer->array[TIME_LOOP]; + log_event(); // do full init/setup since are starting all replicas after event @@ -343,6 +346,7 @@ void PRD::command(int narg, char **arg) update->integrate->setup(); timer->barrier_start(TIME_LOOP); + if (t_corr > 0) replicate(ireplica); if (temp_flag == 0) { if (ireplica == universe->iworld) @@ -350,6 +354,7 @@ void PRD::command(int narg, char **arg) MPI_Bcast(&temp_dephase,1,MPI_DOUBLE,universe->root_proc[ireplica], universe->uworld); } + timer->barrier_stop(TIME_LOOP); time_comm += timer->array[TIME_LOOP]; diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp index 3f4ba9c160..50e352ebd5 100644 --- a/src/REPLICA/tad.cpp +++ b/src/REPLICA/tad.cpp @@ -250,8 +250,11 @@ void TAD::command(int narg, char **arg) // This should work with if uncommented, but does not // if (universe->iworld == 0) { + fix_event->store_state(); quench(); + + timer->init(); timer->barrier_start(TIME_LOOP); time_start = timer->array[TIME_LOOP]; fix_event->store_event_tad(update->ntimestep); @@ -731,7 +734,7 @@ void TAD::perform_neb(int ievent) memory->destroy(buf_init); memory->destroy(buf_final); - // Run NEB + // run NEB int beginstep_hold = update->beginstep; int endstep_hold = update->endstep; @@ -763,7 +766,7 @@ void TAD::perform_neb(int ievent) universe->uscreen = uscreen_lammps; } - // Extract barrier energy from NEB + // extract barrier energy from NEB if (universe->iworld == 0) fix_event_list[ievent]->ebarrier = neb->ebf; @@ -786,7 +789,7 @@ void TAD::perform_neb(int ievent) delete [] args; - // Clean up + // clean up modify->delete_fix("neb"); delete neb; diff --git a/src/REPLICA/temper.cpp b/src/REPLICA/temper.cpp index d95fb3d89b..cd6c3e8e80 100644 --- a/src/REPLICA/temper.cpp +++ b/src/REPLICA/temper.cpp @@ -68,7 +68,8 @@ void Temper::command(int narg, char **arg) error->all(FLERR,"Must have more than one processor partition to temper"); if (domain->box_exist == 0) error->all(FLERR,"Temper command before simulation box is defined"); - if (narg != 6 && narg != 7) error->universe_all(FLERR,"Illegal temper command"); + if (narg != 6 && narg != 7) + error->universe_all(FLERR,"Illegal temper command"); int nsteps = atoi(arg[0]); nevery = atoi(arg[1]); @@ -87,7 +88,8 @@ void Temper::command(int narg, char **arg) // swap frequency must evenly divide total # of timesteps - if (nevery == 0) error->universe_all(FLERR,"Invalid frequency in temper command"); + if (nevery == 0) + error->universe_all(FLERR,"Invalid frequency in temper command"); nswaps = nsteps/nevery; if (nswaps*nevery != nsteps) error->universe_all(FLERR,"Non integer # of swaps in temper command"); @@ -203,6 +205,7 @@ void Temper::command(int narg, char **arg) print_status(); } + timer->init(); timer->barrier_start(TIME_LOOP); for (int iswap = 0; iswap < nswaps; iswap++) { diff --git a/src/comm.cpp b/src/comm.cpp index 8309382b78..a1ca10cdf7 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -50,6 +50,7 @@ using namespace LAMMPS_NS; #define BIG 1.0e20 enum{SINGLE,MULTI}; +enum{NONE,MULTIPLE}; /* ---------------------------------------------------------------------- setup MPI and allocate buffer space @@ -70,6 +71,7 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp) cutghostuser = 0.0; ghost_velocity = 0; recv_from_partition = send_to_partition = -1; + other_partition_style = NONE; // use of OpenMP threads // query OpenMP for number of threads/process set by user at run-time @@ -159,11 +161,23 @@ void Comm::set_processors(int narg, char **arg) int irecv = atoi(arg[iarg+2]); if (isend < 1 || isend > universe->nworlds || irecv < 1 || irecv > universe->nworlds || isend == irecv) - error->all(FLERR,"Invalid partition in processors part command"); - if (isend-1 == universe->iworld) send_to_partition = irecv-1; - if (irecv-1 == universe->iworld) recv_from_partition = isend-1; + error->all(FLERR,"Invalid partitions in processors part command"); + if (isend-1 == universe->iworld) { + if (send_to_partition >= 0) + error->all(FLERR, + "Sending partition in processors part command " + "is already a sender"); + send_to_partition = irecv-1; + } + if (irecv-1 == universe->iworld) { + if (recv_from_partition >= 0) + error->all(FLERR, + "Receiving partition in processors part command " + "is already a receiver"); + recv_from_partition = isend-1; + } if (strcmp(arg[iarg+3],"multiple") == 0) { - if (universe->iworld == irecv-1) other_partition_style = 0; + if (universe->iworld == irecv-1) other_partition_style = MULTIPLE; } else error->all(FLERR,"Illegal processors command"); iarg += 4; } else error->all(FLERR,"Illegal processors command"); @@ -1229,8 +1243,9 @@ void Comm::procs2box() while (ipx <= nprocs) { valid = 1; if (user_procgrid[0] && ipx != user_procgrid[0]) valid = 0; - if (other_procgrid[0]) { - if (other_partition_style == 0 && other_procgrid[0] % ipx) valid = 0; + if (other_partition_style) { + if (other_partition_style == MULTIPLE && other_procgrid[0] % ipx) + valid = 0; } if (nprocs % ipx) valid = 0; if (!valid) { @@ -1242,8 +1257,9 @@ void Comm::procs2box() while (ipy <= nprocs/ipx) { valid = 1; if (user_procgrid[1] && ipy != user_procgrid[1]) valid = 0; - if (other_procgrid[1]) { - if (other_partition_style == 0 && other_procgrid[1] % ipy) valid = 0; + if (other_partition_style) { + if (other_partition_style == MULTIPLE && other_procgrid[1] % ipy) + valid = 0; } if ((nprocs/ipx) % ipy) valid = 0; if (!valid) { @@ -1254,8 +1270,9 @@ void Comm::procs2box() ipz = nprocs/ipx/ipy; valid = 1; if (user_procgrid[2] && ipz != user_procgrid[2]) valid = 0; - if (other_procgrid[2]) { - if (other_partition_style == 0 && other_procgrid[2] % ipz) valid = 0; + if (other_partition_style) { + if (other_partition_style == MULTIPLE && other_procgrid[2] % ipz) + valid = 0; } if (domain->dimension == 2 && ipz != 1) valid = 0; if (!valid) { @@ -1276,7 +1293,8 @@ void Comm::procs2box() ipx++; } - if (procgrid[0] == 0) error->one(FLERR,"Could not layout grid of processors"); + if (procgrid[0] == 0) + error->all(FLERR,"Could not layout grid of processors",1); } /* ---------------------------------------------------------------------- diff --git a/src/error.cpp b/src/error.cpp index c9297eb75d..b7836b45c1 100644 --- a/src/error.cpp +++ b/src/error.cpp @@ -27,6 +27,7 @@ Error::Error(LAMMPS *lmp) : Pointers(lmp) {} /* ---------------------------------------------------------------------- called by all procs in universe close all output, screen, and log files in world and universe + no abort, so insure all procs in universe call, else will hang ------------------------------------------------------------------------- */ void Error::universe_all(const char *file, int line, const char *str) @@ -53,6 +54,7 @@ void Error::universe_all(const char *file, int line, const char *str) /* ---------------------------------------------------------------------- called by one proc in universe + forces abort of entire universe if any proc in universe calls ------------------------------------------------------------------------- */ void Error::universe_one(const char *file, int line, const char *str) @@ -66,9 +68,13 @@ void Error::universe_one(const char *file, int line, const char *str) /* ---------------------------------------------------------------------- called by all procs in one world close all output, screen, and log files in world + insure all procs in world call, else will hang + if abort = 0 (default): + if only one world in universe calls, universe will hang + if abort = 1: force abort of entire universe if any world in universe calls ------------------------------------------------------------------------- */ -void Error::all(const char *file, int line, const char *str) +void Error::all(const char *file, int line, const char *str, int abort) { MPI_Barrier(world); @@ -84,6 +90,7 @@ void Error::all(const char *file, int line, const char *str) if (screen && screen != stdout) fclose(screen); if (logfile) fclose(logfile); + if (abort) MPI_Abort(world,1); MPI_Finalize(); exit(1); } @@ -92,6 +99,7 @@ void Error::all(const char *file, int line, const char *str) called by one proc in world write to world screen only if non-NULL on this proc always write to universe screen + forces abort of entire world (and universe) if any proc in world calls ------------------------------------------------------------------------- */ void Error::one(const char *file, int line, const char *str) @@ -132,6 +140,7 @@ void Error::message(const char *file, int line, char *str, int logflag) /* ---------------------------------------------------------------------- called by all procs in one world close all output, screen, and log files in world + no abort, so insure all procs in world call, else will hang ------------------------------------------------------------------------- */ void Error::done() diff --git a/src/error.h b/src/error.h index 395581dee3..3ab51866bf 100644 --- a/src/error.h +++ b/src/error.h @@ -25,7 +25,7 @@ class Error : protected Pointers { void universe_all(const char *, int, const char *); void universe_one(const char *, int, const char *); - void all(const char *, int, const char *); + void all(const char *, int, const char *, int = 0); void one(const char *, int, const char *); void warning(const char *, int, const char *, int = 1); void message(const char *, int, char *, int = 1); diff --git a/src/lammps.cpp b/src/lammps.cpp index c845599412..146c030fa9 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -482,6 +482,7 @@ void LAMMPS::post_create() /* ---------------------------------------------------------------------- initialize top-level classes + do not initialize Timer class, other classes like Run() do that explicitly ------------------------------------------------------------------------- */ void LAMMPS::init() @@ -500,7 +501,6 @@ void LAMMPS::init() neighbor->init(); // neighbor must come after force, modify comm->init(); // comm must come after force, modify, neighbor, atom output->init(); // output must come after domain, force, modify - timer->init(); } /* ---------------------------------------------------------------------- diff --git a/src/minimize.cpp b/src/minimize.cpp index f93fbf5385..8f8f9d6546 100644 --- a/src/minimize.cpp +++ b/src/minimize.cpp @@ -53,6 +53,7 @@ void Minimize::command(int narg, char **arg) lmp->init(); update->minimize->setup(); + timer->init(); timer->barrier_start(TIME_LOOP); update->minimize->run(update->nsteps); timer->barrier_stop(TIME_LOOP); diff --git a/src/run.cpp b/src/run.cpp index 1f413072a0..001320f6fd 100644 --- a/src/run.cpp +++ b/src/run.cpp @@ -168,11 +168,9 @@ void Run::command(int narg, char **arg) if (preflag || update->first_update == 0) { lmp->init(); update->integrate->setup(); - } else { - timer->init(); - output->setup(0); - } + } else output->setup(0); + timer->init(); timer->barrier_start(TIME_LOOP); update->integrate->run(nsteps); timer->barrier_stop(TIME_LOOP); @@ -208,11 +206,9 @@ void Run::command(int narg, char **arg) if (preflag || iter == 0) { lmp->init(); update->integrate->setup(); - } else { - timer->init(); - output->setup(0); - } + } else output->setup(0); + timer->init(); timer->barrier_start(TIME_LOOP); update->integrate->run(nsteps); timer->barrier_stop(TIME_LOOP); diff --git a/src/verlet_split.cpp b/src/verlet_split.cpp index 39d7185800..22a4564e20 100644 --- a/src/verlet_split.cpp +++ b/src/verlet_split.cpp @@ -215,11 +215,14 @@ void VerletSplit::run(int n) { int nflag,ntimestep,sortflag; - // reset LOOP timer due to imbalance in setup() on 2 partitions - // setup initial Rspace <-> Kspace comm params + // sync both partitions before start timer MPI_Barrier(universe->uworld); - timer->array[TIME_LOOP] = MPI_Wtime(); + timer->init(); + timer->barrier_start(TIME_LOOP); + + // setup initial Rspace <-> Kspace comm params + rk_setup(); // flags for timestepping iterations @@ -315,6 +318,13 @@ void VerletSplit::run(int n) force->kspace->compute(eflag,vflag); timer->stamp(TIME_KSPACE); } + + // TIP4P PPPM puts forces on ghost atoms, so must reverse_comm() + + if (tip4p_flag && force->newton) { + comm->reverse_comm(); + timer->stamp(TIME_COMM); + } } // comm and sum Kspace forces back to Rspace procs @@ -384,8 +394,10 @@ void VerletSplit::rk_setup() MPI_Gatherv(atom->q,n,MPI_DOUBLE,atom->q,qsize,qdisp,MPI_DOUBLE,0,block); - // for TIP4P also need to send type and tag - // KSpace proc needs to acquire ghost atoms and map all its atoms + // for TIP4P also need to send atom type and tag + // KSpace procs need to acquire ghost atoms and map all their atoms + // map_clear() call is in lieu of comm->exchange() which performs map_clear + // borders() call acquires ghost atoms and maps them if (tip4p_flag) { MPI_Gatherv(atom->type,n,MPI_INT,atom->type,qsize,qdisp,MPI_INT,0,block); @@ -394,7 +406,7 @@ void VerletSplit::rk_setup() if (triclinic) domain->x2lamda(atom->nlocal); if (domain->box_change) comm->setup(); timer->stamp(); - atom->map_clear(); // in lieu of comm->exchange() + atom->map_clear(); comm->borders(); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); timer->stamp(TIME_COMM); @@ -442,6 +454,14 @@ void VerletSplit::r2k_comm() force->kspace->setup(); } } + + // for TIP4P, Kspace partition needs to update its ghost atoms + + if (tip4p_flag && !master) { + timer->stamp(); + comm->forward_comm(); + timer->stamp(TIME_COMM); + } } /* ----------------------------------------------------------------------