diff --git a/src/REPLICA/verlet_split.cpp b/src/REPLICA/verlet_split.cpp new file mode 100644 index 0000000000..e84cbc81ec --- /dev/null +++ b/src/REPLICA/verlet_split.cpp @@ -0,0 +1,508 @@ +/* ------------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Yuxing Peng and Chris Knight (U Chicago) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "verlet_split.h" +#include "universe.h" +#include "neighbor.h" +#include "domain.h" +#include "comm.h" +#include "atom.h" +#include "atom_vec.h" +#include "force.h" +#include "pair.h" +#include "bond.h" +#include "angle.h" +#include "dihedral.h" +#include "improper.h" +#include "kspace.h" +#include "output.h" +#include "update.h" +#include "modify.h" +#include "timer.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +VerletSplit::VerletSplit(LAMMPS *lmp, int narg, char **arg) : + Verlet(lmp, narg, arg) +{ + // error checks on partitions + + if (universe->nworlds != 2) + error->universe_all(FLERR,"Verlet/split requires 2 partitions"); + if (universe->procs_per_world[0] % universe->procs_per_world[1]) + error->universe_all(FLERR,"Verlet/split requires Rspace partition " + "size be multiple of Kspace partition size"); + + // master = 1 for Rspace procs, 0 for Kspace procs + + if (universe->iworld == 0) master = 1; + else master = 0; + + ratio = universe->procs_per_world[0] / universe->procs_per_world[1]; + + // Kspace root proc broadcasts info about Kspace proc layout to Rspace procs + + int kspace_procgrid[3]; + + if (universe->me == universe->root_proc[1]) { + kspace_procgrid[0] = comm->procgrid[0]; + kspace_procgrid[1] = comm->procgrid[1]; + kspace_procgrid[2] = comm->procgrid[2]; + } + MPI_Bcast(kspace_procgrid,3,MPI_INT,universe->root_proc[1],universe->uworld); + + int ***kspace_grid2proc; + memory->create(kspace_grid2proc,kspace_procgrid[0], + kspace_procgrid[1],kspace_procgrid[2], + "verlet/split:kspace_grid2proc"); + + if (universe->me == universe->root_proc[1]) { + for (int i = 0; i < comm->procgrid[0]; i++) + for (int j = 0; j < comm->procgrid[1]; j++) + for (int k = 0; k < comm->procgrid[2]; k++) + kspace_grid2proc[i][j][k] = comm->grid2proc[i][j][k]; + } + MPI_Bcast(&kspace_grid2proc[0][0][0], + kspace_procgrid[0]*kspace_procgrid[1]*kspace_procgrid[2],MPI_INT, + universe->root_proc[1],universe->uworld); + + // Rspace partition must be multiple of Kspace partition in each dim + // so atoms of one Kspace proc coincide with atoms of several Rspace procs + + if (master) { + int flag = 0; + if (comm->procgrid[0] % kspace_procgrid[0]) flag = 1; + if (comm->procgrid[1] % kspace_procgrid[1]) flag = 1; + if (comm->procgrid[2] % kspace_procgrid[2]) flag = 1; + if (flag) + error->one(FLERR, + "Verlet/split requires Rspace partition layout be " + "multiple of Kspace partition layout in each dim"); + } + + // block = 1 Kspace proc with set of Rspace procs it overlays + // me_block = 0 for Kspace proc + // me_block = 1 to ratio for Rspace procs + // block = MPI communicator for that set of procs + + int iblock,key; + + if (!master) { + iblock = comm->me; + key = 0; + } else { + int kpx = comm->myloc[0] / (comm->procgrid[0]/kspace_procgrid[0]); + int kpy = comm->myloc[1] / (comm->procgrid[1]/kspace_procgrid[1]); + int kpz = comm->myloc[2] / (comm->procgrid[2]/kspace_procgrid[2]); + iblock = kspace_grid2proc[kpx][kpy][kpz]; + key = 1; + } + + MPI_Comm_split(universe->uworld,iblock,key,&block); + MPI_Comm_rank(block,&me_block); + + // output block groupings to universe screen/logfile + // bmap is ordered by block and then by proc within block + + int *bmap = new int[universe->nprocs]; + for (int i = 0; i < universe->nprocs; i++) bmap[i] = -1; + bmap[iblock*(ratio+1)+me_block] = universe->me; + + int *bmapall = new int[universe->nprocs]; + MPI_Allreduce(bmap,bmapall,universe->nprocs,MPI_INT,MPI_MAX,universe->uworld); + + if (universe->me == 0) { + if (universe->uscreen) { + fprintf(universe->uscreen,"Rspace/Kspace procs in each block:\n"); + int m = 0; + for (int i = 0; i < universe->nprocs/(ratio+1); i++) { + fprintf(universe->uscreen," block %d:",i); + int kspace_proc = bmapall[m]; + for (int j = 1; j <= ratio; j++) + fprintf(universe->uscreen," %d",bmapall[m+j]); + fprintf(universe->uscreen," %d\n",kspace_proc); + /* + kspace_proc = bmapall[m]; + for (int j = 1; j <= ratio; j++) + fprintf(universe->uscreen," %d", + universe->proc2original[bmapall[m+j]]); + fprintf(universe->uscreen," %d\n",universe->proc2original[kspace_proc]); + */ + m += ratio + 1; + } + } + if (universe->ulogfile) { + fprintf(universe->ulogfile,"Rspace/Kspace procs in each block:\n"); + int m = 0; + for (int i = 0; i < universe->nprocs/(ratio+1); i++) { + fprintf(universe->ulogfile," block %d:",i); + int kspace_proc = bmapall[m++]; + for (int j = 1; j <= ratio; j++) + fprintf(universe->ulogfile," %d",bmapall[m++]); + fprintf(universe->ulogfile," %d\n",kspace_proc); + } + } + } + + memory->destroy(kspace_grid2proc); + delete [] bmap; + delete [] bmapall; + + // size/disp = vectors for MPI gather/scatter within block + + qsize = new int[ratio+1]; + qdisp = new int[ratio+1]; + xsize = new int[ratio+1]; + xdisp = new int[ratio+1]; + + // f_kspace = Rspace copy of Kspace forces + // allocate dummy version for Kspace partition + + maxatom = 0; + f_kspace = NULL; + if (!master) memory->create(f_kspace,1,1,"verlet/split:f_kspace"); +} + +/* ---------------------------------------------------------------------- */ + +VerletSplit::~VerletSplit() +{ + delete [] qsize; + delete [] qdisp; + delete [] xsize; + delete [] xdisp; + memory->destroy(f_kspace); + MPI_Comm_free(&block); +} + +/* ---------------------------------------------------------------------- + initialization before run +------------------------------------------------------------------------- */ + +void VerletSplit::init() +{ + if (!force->kspace && comm->me == 0) + error->warning(FLERR,"No Kspace calculation with verlet/split"); + + if (force->kspace_match("tip4p",0)) tip4p_flag = 1; + else tip4p_flag = 0; + + Verlet::init(); +} + +/* ---------------------------------------------------------------------- + run for N steps + master partition does everything but Kspace + servant partition does just Kspace + communicate back and forth every step: + atom coords from master -> servant + kspace forces from servant -> master + also box bounds from master -> servant if necessary +------------------------------------------------------------------------- */ + +void VerletSplit::run(int n) +{ + int nflag,ntimestep,sortflag; + + // sync both partitions before start timer + + MPI_Barrier(universe->uworld); + timer->init(); + timer->barrier_start(TIME_LOOP); + + // setup initial Rspace <-> Kspace comm params + + rk_setup(); + + // flags for timestepping iterations + + int n_post_integrate = modify->n_post_integrate; + int n_pre_exchange = modify->n_pre_exchange; + int n_pre_neighbor = modify->n_pre_neighbor; + int n_pre_force = modify->n_pre_force; + int n_post_force = modify->n_post_force; + int n_end_of_step = modify->n_end_of_step; + + if (atom->sortfreq > 0) sortflag = 1; + else sortflag = 0; + + for (int i = 0; i < n; i++) { + + ntimestep = ++update->ntimestep; + ev_set(ntimestep); + + // initial time integration + + if (master) { + modify->initial_integrate(vflag); + if (n_post_integrate) modify->post_integrate(); + } + + // regular communication vs neighbor list rebuild + + if (master) nflag = neighbor->decide(); + MPI_Bcast(&nflag,1,MPI_INT,1,block); + + if (master) { + if (nflag == 0) { + timer->stamp(); + comm->forward_comm(); + timer->stamp(TIME_COMM); + } else { + if (n_pre_exchange) modify->pre_exchange(); + if (triclinic) domain->x2lamda(atom->nlocal); + domain->pbc(); + if (domain->box_change) { + domain->reset_box(); + comm->setup(); + if (neighbor->style) neighbor->setup_bins(); + } + timer->stamp(); + comm->exchange(); + if (sortflag && ntimestep >= atom->nextsort) atom->sort(); + comm->borders(); + if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + timer->stamp(TIME_COMM); + if (n_pre_neighbor) modify->pre_neighbor(); + neighbor->build(); + timer->stamp(TIME_NEIGHBOR); + } + } + + // if reneighboring occurred, re-setup Rspace <-> Kspace comm params + // comm Rspace atom coords to Kspace procs + + if (nflag) rk_setup(); + r2k_comm(); + + // force computations + + force_clear(); + + if (master) { + if (n_pre_force) modify->pre_force(vflag); + + timer->stamp(); + if (force->pair) { + force->pair->compute(eflag,vflag); + timer->stamp(TIME_PAIR); + } + + if (atom->molecular) { + if (force->bond) force->bond->compute(eflag,vflag); + if (force->angle) force->angle->compute(eflag,vflag); + if (force->dihedral) force->dihedral->compute(eflag,vflag); + if (force->improper) force->improper->compute(eflag,vflag); + timer->stamp(TIME_BOND); + } + + if (force->newton) { + comm->reverse_comm(); + timer->stamp(TIME_COMM); + } + + } else { + if (force->kspace) { + timer->stamp(); + 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 + + k2r_comm(); + + // force modifications, final time integration, diagnostics + // all output + + if (master) { + if (n_post_force) modify->post_force(vflag); + modify->final_integrate(); + if (n_end_of_step) modify->end_of_step(); + + if (ntimestep == output->next) { + timer->stamp(); + output->write(ntimestep); + timer->stamp(TIME_OUTPUT); + } + } + } +} + +/* ---------------------------------------------------------------------- + setup params for Rspace <-> Kspace communication + called initially and after every reneighbor + also communcicate atom charges from Rspace to KSpace since static +------------------------------------------------------------------------- */ + +void VerletSplit::rk_setup() +{ + // grow f_kspace array on master procs if necessary + + if (master) { + if (atom->nlocal > maxatom) { + memory->destroy(f_kspace); + maxatom = atom->nmax; + memory->create(f_kspace,maxatom,3,"verlet/split:f_kspace"); + } + } + + // qsize = # of atoms owned by each master proc in block + + int n = 0; + if (master) n = atom->nlocal; + MPI_Gather(&n,1,MPI_INT,qsize,1,MPI_INT,0,block); + + // setup qdisp, xsize, xdisp based on qsize + // only needed by Kspace proc + // set Kspace nlocal to sum of Rspace nlocals + // insure Kspace atom arrays are large enough + + if (!master) { + qsize[0] = qdisp[0] = xsize[0] = xdisp[0] = 0; + for (int i = 1; i <= ratio; i++) { + qdisp[i] = qdisp[i-1]+qsize[i-1]; + xsize[i] = 3*qsize[i]; + xdisp[i] = xdisp[i-1]+xsize[i-1]; + } + + atom->nlocal = qdisp[ratio] + qsize[ratio]; + while (atom->nmax <= atom->nlocal) atom->avec->grow(0); + atom->nghost = 0; + } + + // one-time gather of Rspace atom charges to Kspace proc + + MPI_Gatherv(atom->q,n,MPI_DOUBLE,atom->q,qsize,qdisp,MPI_DOUBLE,0,block); + + // 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); + MPI_Gatherv(atom->tag,n,MPI_INT,atom->tag,qsize,qdisp,MPI_INT,0,block); + if (!master) { + if (triclinic) domain->x2lamda(atom->nlocal); + if (domain->box_change) comm->setup(); + timer->stamp(); + atom->map_clear(); + comm->borders(); + if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + timer->stamp(TIME_COMM); + } + } +} + +/* ---------------------------------------------------------------------- + communicate Rspace atom coords to Kspace + also eflag,vflag and box bounds if needed +------------------------------------------------------------------------- */ + +void VerletSplit::r2k_comm() +{ + MPI_Status status; + + int n = 0; + if (master) n = atom->nlocal; + MPI_Gatherv(atom->x[0],n*3,MPI_DOUBLE,atom->x[0],xsize,xdisp, + MPI_DOUBLE,0,block); + + // send eflag,vflag from Rspace to Kspace + + if (me_block == 1) { + int flags[2]; + flags[0] = eflag; flags[1] = vflag; + MPI_Send(flags,2,MPI_INT,0,0,block); + } else if (!master) { + int flags[2]; + MPI_Recv(flags,2,MPI_DOUBLE,1,0,block,&status); + eflag = flags[0]; vflag = flags[1]; + } + + // send box bounds from Rspace to Kspace if simulation box is dynamic + + if (domain->box_change) { + if (me_block == 1) { + MPI_Send(domain->boxlo,3,MPI_DOUBLE,0,0,block); + MPI_Send(domain->boxhi,3,MPI_DOUBLE,0,0,block); + } else if (!master) { + MPI_Recv(domain->boxlo,3,MPI_DOUBLE,1,0,block,&status); + MPI_Recv(domain->boxhi,3,MPI_DOUBLE,1,0,block,&status); + domain->set_global_box(); + domain->set_local_box(); + 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); + } +} + +/* ---------------------------------------------------------------------- + communicate and sum Kspace atom forces back to Rspace +------------------------------------------------------------------------- */ + +void VerletSplit::k2r_comm() +{ + if (eflag) MPI_Bcast(&force->kspace->energy,1,MPI_DOUBLE,0,block); + if (vflag) MPI_Bcast(force->kspace->virial,6,MPI_DOUBLE,0,block); + + int n = 0; + if (master) n = atom->nlocal; + MPI_Scatterv(atom->f[0],xsize,xdisp,MPI_DOUBLE, + f_kspace[0],n*3,MPI_DOUBLE,0,block); + + if (master) { + double **f = atom->f; + int nlocal = atom->nlocal; + for (int i = 0; i < nlocal; i++) { + f[i][0] += f_kspace[i][0]; + f[i][1] += f_kspace[i][1]; + f[i][2] += f_kspace[i][2]; + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of Kspace force array on master procs +------------------------------------------------------------------------- */ + +bigint VerletSplit::memory_usage() +{ + bigint bytes = maxatom*3 * sizeof(double); + return bytes; +} diff --git a/src/REPLICA/verlet_split.h b/src/REPLICA/verlet_split.h new file mode 100644 index 0000000000..080a6057eb --- /dev/null +++ b/src/REPLICA/verlet_split.h @@ -0,0 +1,54 @@ +/* ------------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef INTEGRATE_CLASS + +IntegrateStyle(verlet/split,VerletSplit) + +#else + +#ifndef LMP_VERLET_SPLIT_H +#define LMP_VERLET_SPLIT_H + +#include "verlet.h" + +namespace LAMMPS_NS { + +class VerletSplit : public Verlet { + public: + VerletSplit(class LAMMPS *, int, char **); + ~VerletSplit(); + void init(); + void run(int); + bigint memory_usage(); + + private: + int master; // 1 if an Rspace proc, 0 if Kspace + int me_block; // proc ID within Rspace/Kspace block + int ratio; // ratio of Rspace procs to Kspace procs + int *qsize,*qdisp,*xsize,*xdisp; // MPI gather/scatter params for block comm + MPI_Comm block; // communicator within one block + int tip4p_flag; // 1 if PPPM/tip4p so do extra comm + + double **f_kspace; // copy of Kspace forces on Rspace procs + int maxatom; + + void rk_setup(); + void r2k_comm(); + void k2r_comm(); +}; + +} + +#endif +#endif