git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7303 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/KSPACE/pppm.cpp

@@ -115,7 +115,8 @@ void PPPM::init()
 
   if (domain->triclinic)
     error->all(FLERR,"Cannot (yet) use PPPM with triclinic box");
-  if (domain->dimension == 2) error->all(FLERR,"Cannot use PPPM with 2d simulation");
+  if (domain->dimension == 2) error->all(FLERR,
+					 "Cannot use PPPM with 2d simulation");
 
   if (!atom->q_flag) error->all(FLERR,"Kspace style requires atom attribute q");
 
@@ -186,8 +187,8 @@ void PPPM::init()
 
   // if we have a /proxy pppm version check if the pair style is compatible
 
-  if ( (strcmp(force->kspace_style,"pppm/proxy") == 0) ||
+  if ((strcmp(force->kspace_style,"pppm/proxy") == 0) ||
-       (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0) ) {
+      (strcmp(force->kspace_style,"pppm/tip4p/proxy") == 0) ) {
     if (force->pair == NULL)
       error->all(FLERR,"KSpace style is incompatible with Pair style");
     if (strstr(force->pair_style,"pppm/") == NULL )
@@ -1507,7 +1508,8 @@ void PPPM::particle_map()
 
     if (nx+nlower < nxlo_out || nx+nupper > nxhi_out ||
 	ny+nlower < nylo_out || ny+nupper > nyhi_out ||
-	nz+nlower < nzlo_out || nz+nupper > nzhi_out) flag = 1;
+	nz+nlower < nzlo_out || nz+nupper > nzhi_out)
+      flag = 1;
   }
 
   if (flag) error->one(FLERR,"Out of range atoms - cannot compute PPPM");
@@ -1527,7 +1529,8 @@ void PPPM::make_rho()
 
   // clear 3d density array
 
-  memset(&(density_brick[nzlo_out][nylo_out][nxlo_out]),0,ngrid*sizeof(FFT_SCALAR));
+  memset(&(density_brick[nzlo_out][nylo_out][nxlo_out]),0,
+	 ngrid*sizeof(FFT_SCALAR));
 
   // loop over my charges, add their contribution to nearby grid points
   // (nx,ny,nz) = global coords of grid pt to "lower left" of charge
@@ -1784,6 +1787,7 @@ void PPPM::procs2grid2d(int nprocs, int nx, int ny, int *px, int *py)
    charge assignment into rho1d
    dx,dy,dz = distance of particle from "lower left" grid point 
 ------------------------------------------------------------------------- */
+
 void PPPM::compute_rho1d(const FFT_SCALAR &dx, const FFT_SCALAR &dy,
 			 const FFT_SCALAR &dz)
 {

src/KSPACE/pppm.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

src/PERI/pair_peri_lps.cpp

@@ -41,7 +41,7 @@ using namespace LAMMPS_NS;
 PairPeriLPS::PairPeriLPS(LAMMPS *lmp) : Pair(lmp)			
 {
   for (int i = 0; i < 6; i++) virial[i] = 0.0;
-  no_virial_fdotr_compute=1;
+  no_virial_fdotr_compute = 1;
 
   ifix_peri = -1;
 
@@ -189,7 +189,8 @@ void PairPeriLPS::compute(int eflag, int vflag)
         }
 
         if (eflag) evdwl = 0.5*rk*dr;
-	if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair*vfrac[i],delx,dely,delz);
+	if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,
+			     fpair*vfrac[i],delx,dely,delz);
       }
     }
   }
@@ -303,7 +304,8 @@ void PairPeriLPS::compute(int eflag, int vflag)
       if (eflag) evdwl = 0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * 
 		   omega_plus*(deviatoric_extension * deviatoric_extension) *
 		   vfrac[j] * vfrac_scale;
-      if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0,0.5*fbond*vfrac[i],delx,dely,delz);
+      if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0,
+			   0.5*fbond*vfrac[i],delx,dely,delz);
 
       // find stretch in bond I-J and break if necessary
       // use s0 from previous timestep
@@ -316,7 +318,8 @@ void PairPeriLPS::compute(int eflag, int vflag)
       if (first)
          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;
     }

src/comm.cpp

@@ -22,6 +22,7 @@
 #include "stdio.h"
 #include "stdlib.h"
 #include "comm.h"
+#include "universe.h"
 #include "atom.h"
 #include "atom_vec.h"
 #include "force.h"
@@ -68,6 +69,7 @@ Comm::Comm(LAMMPS *lmp) : Pointers(lmp)
   cutghostmulti = NULL;
   cutghostuser = 0.0;
   ghost_velocity = 0;
+  recv_from_partition = send_to_partition = -1;
 
   // use of OpenMP threads
   // query OpenMP for number of threads/process set by user at run-time
@@ -126,14 +128,70 @@ Comm::~Comm()
   memory->destroy(buf_recv);
 }
 
+/* ----------------------------------------------------------------------
+   set dimensions of processor grid
+   invoked from input script by processors command
+------------------------------------------------------------------------- */
+
+void Comm::set_processors(int narg, char **arg)
+{
+  if (narg < 3) error->all(FLERR,"Illegal processors command");
+
+  if (strcmp(arg[0],"*") == 0) user_procgrid[0] = 0;
+  else user_procgrid[0] = atoi(arg[0]);
+  if (strcmp(arg[1],"*") == 0) user_procgrid[1] = 0;
+  else user_procgrid[1] = atoi(arg[1]);
+  if (strcmp(arg[2],"*") == 0) user_procgrid[2] = 0;
+  else user_procgrid[2] = atoi(arg[2]);
+
+  if (user_procgrid[0] < 0 || user_procgrid[1] < 0 || user_procgrid[2] < 0) 
+    error->all(FLERR,"Illegal processors command");
+
+  int iarg = 3;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"part") == 0) {
+      if (iarg+4 > narg) error->all(FLERR,"Illegal processors command");
+      if (universe->nworlds == 1)
+	error->all(FLERR,
+		   "Cannot use processors part command "
+		   "without using partitions");
+      int isend = atoi(arg[iarg+1]);
+      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;
+      if (strcmp(arg[iarg+3],"multiple") == 0) {
+	if (universe->iworld == irecv-1) other_partition_style = 0;
+      } else error->all(FLERR,"Illegal processors command");
+      iarg += 4;
+    } else error->all(FLERR,"Illegal processors command");
+  }
+}
+
 /* ----------------------------------------------------------------------
    setup 3d grid of procs based on box size
 ------------------------------------------------------------------------- */
 
-void Comm::set_procs()
+void Comm::set_proc_grid()
 {
+  // recv proc layout of another partition if my layout depends on it
+
+  if (recv_from_partition >= 0) {
+    MPI_Status status;
+    if (me == 0) MPI_Recv(other_procgrid,3,MPI_INT,
+			  universe->root_proc[recv_from_partition],0,
+			  universe->uworld,&status);
+    MPI_Bcast(other_procgrid,3,MPI_INT,0,world);
+  }
+
+  // create layout of procs mapped to simulation box
+
   procs2box();
 
+  // error check
+
   if (procgrid[0]*procgrid[1]*procgrid[2] != nprocs)
     error->all(FLERR,"Bad grid of processors");
   if (domain->dimension == 2 && procgrid[2] != 1)
@@ -179,6 +237,14 @@ void Comm::set_procs()
     if (logfile) fprintf(logfile,"  %d by %d by %d MPI processor grid\n",
 			 procgrid[0],procgrid[1],procgrid[2]);
   }
+
+  // send my proc layout to another partition if requested
+
+  if (send_to_partition >= 0) {
+    if (me == 0) MPI_Send(procgrid,3,MPI_INT,
+			  universe->root_proc[send_to_partition],0,
+			  universe->uworld);
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -1149,8 +1215,12 @@ void Comm::procs2box()
   double bestsurf = 2.0 * (area[0]+area[1]+area[2]);
 
   // loop thru all possible factorizations of nprocs
-  // only consider valid cases that match procgrid settings
+  // choose factorization with least surface area
   // surf = surface area of a proc sub-domain
+  // only consider factorizations that match procgrid & other_procgrid settings
+  // if set, only consider factorizations that match other_procgrid settings
+
+  procgrid[0] = procgrid[1] = procgrid[2] = 0;
 
   int ipx,ipy,ipz,valid;
   double surf;
@@ -1159,6 +1229,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 (nprocs % ipx) valid = 0;
     if (!valid) {
       ipx++;
@@ -1169,6 +1242,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 ((nprocs/ipx) % ipy) valid = 0;
       if (!valid) {
 	ipy++;
@@ -1178,6 +1254,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 (domain->dimension == 2 && ipz != 1) valid = 0;
       if (!valid) {
 	ipy++;
@@ -1196,6 +1275,8 @@ void Comm::procs2box()
 
     ipx++;
   }
+
+  if (procgrid[0] == 0) error->one(FLERR,"Could not layout grid of processors");
 }
 
 /* ----------------------------------------------------------------------

src/comm.h

@@ -29,18 +29,24 @@ class Comm : protected Pointers {
   double cutghost[3];               // cutoffs used for acquiring ghost atoms
   double cutghostuser;              // user-specified ghost cutoff
   int ***grid2proc;                 // which proc owns i,j,k loc in 3d grid
+  int recv_from_partition;          // recv proc layout from this partition
+  int send_to_partition;            // send my proc layout to this partition
+                                    // -1 if no recv or send
+  int other_partition_style;        // 0 = recv layout dims must be multiple of
+                                    //     my layout dims
   int nthreads;                     // OpenMP threads per MPI process
 
   Comm(class LAMMPS *);
   virtual ~Comm();
 
   virtual void init();
-  virtual void set_procs();                 // setup 3d grid of procs
+  void set_processors(int, char **);          // set procs from input script
-  virtual void setup();                     // setup 3d communication pattern
+  virtual void set_proc_grid();               // setup 3d grid of procs
-  virtual void forward_comm(int dummy = 0); // forward communication of atom coords
+  virtual void setup();                       // setup 3d comm pattern
-  virtual void reverse_comm();              // reverse communication of forces
+  virtual void forward_comm(int dummy = 0);   // forward comm of atom coords
-  virtual void exchange();                  // move atoms to new procs
+  virtual void reverse_comm();                // reverse comm of forces
-  virtual void borders();                   // setup list of atoms to communicate
+  virtual void exchange();                    // move atoms to new procs
+  virtual void borders();                     // setup list of atoms to comm
 
   virtual void forward_comm_pair(class Pair *);    // forward comm from a Pair
   virtual void reverse_comm_pair(class Pair *);    // reverse comm from a Pair
@@ -77,6 +83,8 @@ class Comm : protected Pointers {
   int map_style;                    // non-0 if global->local mapping is done
   int bordergroup;                  // only communicate this group in borders
 
+  int other_procgrid[3];            // proc layout
+
   int *firstrecv;                   // where to put 1st recv atom in each swap
   int **sendlist;                   // list of atoms to send in each swap
   int *maxsendlist;                 // max size of send list for each swap

src/create_box.cpp

@@ -110,6 +110,6 @@ void CreateBox::command(int narg, char **arg)
   domain->print_box("Created ");
   domain->set_initial_box();
   domain->set_global_box();
-  comm->set_procs();
+  comm->set_proc_grid();
   domain->set_local_box();
 }

src/domain.cpp

@@ -603,7 +603,9 @@ void Domain::minimum_image(double *delta)
    for triclinic, also add/subtract tilt factors in other dims as needed
 ------------------------------------------------------------------------- */
 
-void Domain::closest_image(const double * const xi, const double * const xj, double * const xjimage)
+void Domain::closest_image(const double * const xi, 
+			   const double * const xj,
+			   double * const xjimage)
 {
   double dx,dy,dz;
 

src/force.cpp

@@ -534,6 +534,20 @@ KSpace *Force::new_kspace(int narg, char **arg, const char *suffix, int &sflag)
   return NULL;
 }
 
+/* ----------------------------------------------------------------------
+   return ptr to Kspace class if matches word
+   if exact, then style name must be exact match to word
+   if not exact, style name must contain word
+   return NULL if no match
+------------------------------------------------------------------------- */
+
+KSpace *Force::kspace_match(const char *word, int exact)
+{
+  if (exact && strcmp(kspace_style,word) == 0) return kspace;
+  else if (!exact && strstr(kspace_style,word)) return kspace;
+  return NULL;
+}
+
 /* ----------------------------------------------------------------------
    set special bond values 
 ------------------------------------------------------------------------- */
@@ -627,7 +641,7 @@ void Force::set_special(int narg, char **arg)
 
 /* ----------------------------------------------------------------------
    compute bounds implied by numeric str with a possible wildcard asterik
-   nmax = upper bound
+   1 = lower bound, nmax = upper bound
    5 possibilities:
      (1) i = i to i, (2) * = 1 to nmax,
      (3) i* = i to nmax, (4) *j = 1 to j, (5) i*j = i to j
@@ -639,21 +653,22 @@ void Force::bounds(char *str, int nmax, int &nlo, int &nhi)
   char *ptr = strchr(str,'*');
 
   if (ptr == NULL) {
-    nlo = MAX(atoi(str),1);
+    nlo = nhi = atoi(str);
-    nhi = MIN(atoi(str),nmax);
   } else if (strlen(str) == 1) {
     nlo = 1;
     nhi = nmax;
   } else if (ptr == str) {
     nlo = 1;
-    nhi = MIN(atoi(ptr+1),nmax);
+    nhi = atoi(ptr+1);
   } else if (strlen(ptr+1) == 0) {
-    nlo = MAX(atoi(str),1);
+    nlo = atoi(str);
     nhi = nmax;
   } else {
-    nlo = MAX(atoi(str),1);
+    nlo = atoi(str);
-    nhi = MIN(atoi(ptr+1),nmax);
+    nhi = atoi(ptr+1);
   }
+
+  if (nlo < 1 || nhi > nmax) error->all(FLERR,"Numeric index is out of bounds");
 }
 
 /* ----------------------------------------------------------------------

src/force.h

@@ -88,6 +88,7 @@ class Force : protected Pointers {
 
   void create_kspace(int, char **, const char *suffix = NULL);
   class KSpace *new_kspace(int, char **, const char *, int &);
+  class KSpace *kspace_match(const char *, int);
 
   void set_special(int, char **);
   void bounds(char *, int, int &, int &);

src/input.cpp

@@ -418,6 +418,7 @@ int Input::execute_command()
   else if (!strcmp(command,"label")) label();
   else if (!strcmp(command,"log")) log();
   else if (!strcmp(command,"next")) next_command();
+  else if (!strcmp(command,"partition")) partition();
   else if (!strcmp(command,"print")) print();
   else if (!strcmp(command,"shell")) shell();
   else if (!strcmp(command,"variable")) variable_command();
@@ -580,7 +581,7 @@ void Input::ifthenelse()
       ncommands++;
     }
     
-    for (int i = 0; i < ncommands; i++) input->one(commands[i]);
+    for (int i = 0; i < ncommands; i++) one(commands[i]);
     
     for (int i = 0; i < ncommands; i++) delete [] commands[i];
     delete [] commands;
@@ -636,7 +637,7 @@ void Input::ifthenelse()
     
     // execute the list of commands
     
-    for (int i = 0; i < ncommands; i++) input->one(commands[i]);
+    for (int i = 0; i < ncommands; i++) one(commands[i]);
     
     // clean up
 
@@ -742,6 +743,40 @@ void Input::next_command()
 
 /* ---------------------------------------------------------------------- */
 
+void Input::partition()
+{
+  if (narg < 3) error->all(FLERR,"Illegal partition command");
+
+  int yesflag;
+  if (strcmp(arg[0],"yes") == 0) yesflag = 1;
+  else if (strcmp(arg[0],"no") == 0) yesflag = 0;
+  else error->all(FLERR,"Illegal partition command");
+
+  int ilo,ihi;
+  force->bounds(arg[1],universe->nworlds,ilo,ihi);
+
+  // copy original line to copy, since will use strtok() on it
+  // ptr = start of 4th word
+
+  strcpy(copy,line);
+  copy[strlen(copy)-1] = '\0';
+  char *ptr = strtok(copy," \t\n\r\f");
+  ptr = strtok(NULL," \t\n\r\f");
+  ptr = strtok(NULL," \t\n\r\f");
+  ptr += strlen(ptr) + 1;
+  ptr += strspn(ptr," \t\n\r\f");
+
+  // execute the remaining command line on requested partitions
+
+  if (yesflag) {
+    if (universe->iworld+1 >= ilo && universe->iworld+1 <= ihi) one(ptr);
+  } else {
+    if (universe->iworld+1 < ilo || universe->iworld+1 > ihi) one(ptr);
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void Input::print()
 {
   if (narg != 1) error->all(FLERR,"Illegal print command");
@@ -1204,19 +1239,9 @@ void Input::pair_write()
 
 void Input::processors()
 {
-  if (narg != 3) error->all(FLERR,"Illegal processors command");
   if (domain->box_exist)
     error->all(FLERR,"Processors command after simulation box is defined");
-
+  comm->set_processors(narg,arg);
-  if (strcmp(arg[0],"*") == 0) comm->user_procgrid[0] = 0;
-  else comm->user_procgrid[0] = atoi(arg[0]);
-  if (strcmp(arg[1],"*") == 0) comm->user_procgrid[1] = 0;
-  else comm->user_procgrid[1] = atoi(arg[1]);
-  if (strcmp(arg[2],"*") == 0) comm->user_procgrid[2] = 0;
-  else comm->user_procgrid[2] = atoi(arg[2]);
-
-  if (comm->user_procgrid[0] < 0 || comm->user_procgrid[1] < 0 ||
-      comm->user_procgrid[2] < 0) error->all(FLERR,"Illegal processors command");
 }
 
 /* ---------------------------------------------------------------------- */

src/input.h

@@ -57,6 +57,7 @@ class Input : protected Pointers {
   void label();
   void log();
   void next_command();
+  void partition();
   void print();
   void shell();
   void variable_command();

src/read_data.cpp

@@ -137,7 +137,7 @@ void ReadData::command(int narg, char **arg)
   domain->print_box("  ");
   domain->set_initial_box();
   domain->set_global_box();
-  comm->set_procs();
+  comm->set_proc_grid();
   domain->set_local_box();
 
   // customize for new sections

src/read_restart.cpp

@@ -135,7 +135,7 @@ void ReadRestart::command(int narg, char **arg)
   domain->print_box("  ");
   domain->set_initial_box();
   domain->set_global_box();
-  comm->set_procs();
+  comm->set_proc_grid();
   domain->set_local_box();
 
   // read groups, ntype-length arrays, force field, fix info from file

src/replicate.cpp

@@ -193,7 +193,7 @@ void Replicate::command(int narg, char **arg)
   domain->print_box("  ");
   domain->set_initial_box();
   domain->set_global_box();
-  comm->set_procs();
+  comm->set_proc_grid();
   domain->set_local_box();
 
   // copy type arrays to new atom class

src/verlet_split.cpp

@@ -0,0 +1,480 @@
+/* -------------------------------------------------------------------------
+   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++]);
+	fprintf(universe->uscreen," %d\n",kspace_proc);
+      }
+    }
+    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;
+
+  // reset LOOP timer due to imbalance in setup() on 2 partitions
+  // setup initial Rspace <-> Kspace comm params
+
+  MPI_Barrier(universe->uworld);
+  timer->array[TIME_LOOP] = MPI_Wtime();
+  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);
+      }
+    }
+
+    // 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 type and tag
+  // KSpace proc needs to acquire ghost atoms and map all its atoms
+
+  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();          // in lieu of comm->exchange()
+      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();
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   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;
+}

src/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