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

src/STUBS/mpi.cpp

@@ -14,6 +14,7 @@
 /* Single-processor "stub" versions of MPI routines */
 
 #include "stdlib.h"
+#include "string.h"
 #include "stdio.h"
 #include <sys/time.h>
 #include "mpi.h"
@@ -26,27 +27,46 @@ void mpi_copy_double(void *, void *, int);
 void mpi_copy_char(void *, void *, int);
 void mpi_copy_byte(void *, void *, int);
 
+/* lo-level data structure */
+
+struct {
+  double value;
+  int proc;
+} double_int;
+
+/* ---------------------------------------------------------------------- */
 /* MPI Functions */
+/* ---------------------------------------------------------------------- */
 
 void MPI_Init(int *argc, char ***argv) {}
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Comm_rank(MPI_Comm comm, int *me)
 {
   *me = 0;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Comm_size(MPI_Comm comm, int *nprocs)
 {
   *nprocs = 1;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Abort(MPI_Comm comm, int errorcode)
 {
   exit(1);
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Finalize() {}
 
+/* ---------------------------------------------------------------------- */
+
 double MPI_Wtime()
 {
   double time;
@@ -57,46 +77,62 @@ double MPI_Wtime()
   return time;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Send(void *buf, int count, MPI_Datatype datatype,
 	      int dest, int tag, MPI_Comm comm)
 {
   printf("MPI Stub WARNING: Should not send message to self\n");
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Rsend(void *buf, int count, MPI_Datatype datatype,
 	       int dest, int tag, MPI_Comm comm)
 {
   printf("MPI Stub WARNING: Should not rsend message to self\n");
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Recv(void *buf, int count, MPI_Datatype datatype,
 	      int source, int tag, MPI_Comm comm, MPI_Status *status)
 {
   printf("MPI Stub WARNING: Should not recv message from self\n");
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Irecv(void *buf, int count, MPI_Datatype datatype,
 	       int source, int tag, MPI_Comm comm, MPI_Request *request)
 {
   printf("MPI Stub WARNING: Should not recv message from self\n");
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Wait(MPI_Request *request, MPI_Status *status)
 {
   printf("MPI Stub WARNING: Should not wait on message from self\n");
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Waitall(int n, MPI_Request *request, MPI_Status *status)
 {
   printf("MPI Stub WARNING: Should not wait on message from self\n");
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Waitany(int count, MPI_Request *request, int *index, 
 		 MPI_Status *status)
 {
   printf("MPI Stub WARNING: Should not wait on message from self\n");
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Sendrecv(void *sbuf, int scount, MPI_Datatype sdatatype,
 		  int dest, int stag, void *rbuf, int rcount,
 		  MPI_Datatype rdatatype, int source, int rtag,
@@ -105,29 +141,41 @@ void MPI_Sendrecv(void *sbuf, int scount, MPI_Datatype sdatatype,
   printf("MPI Stub WARNING: Should not send message to self\n");
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Get_count(MPI_Status *status, MPI_Datatype datatype, int *count)
 {
   printf("MPI Stub WARNING: Should not get count of message to self\n");
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Comm_split(MPI_Comm comm, int color, int key, MPI_Comm *comm_out)
 {
   *comm_out = comm;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Comm_dup(MPI_Comm comm, MPI_Comm *comm_out)
 {
   *comm_out = comm;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Comm_free(MPI_Comm *comm) { }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Cart_create(MPI_Comm comm_old, int ndims, int *dims, int *periods,
 		     int reorder, MPI_Comm *comm_cart)
 {
   *comm_cart = comm_old;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Cart_get(MPI_Comm comm, int maxdims, int *dims, int *periods,
 		  int *coords)
 {
@@ -136,169 +184,135 @@ void MPI_Cart_get(MPI_Comm comm, int maxdims, int *dims, int *periods,
   coords[0] = coords[1] = coords[2] = 0;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Cart_shift(MPI_Comm comm, int direction, int displ,
 		    int *source, int *dest)
 {
   *source = *dest = 0;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Cart_rank(MPI_Comm comm, int *coords, int *rank)
 {
   *rank = 0;
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Barrier(MPI_Comm comm) {}
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Bcast(void *buf, int count, MPI_Datatype datatype,
 	       int root, MPI_Comm comm) {}
 
+/* ---------------------------------------------------------------------- */
+
 /* copy values from data1 to data2 */
 
 void MPI_Allreduce(void *sendbuf, void *recvbuf, int count,
 		   MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
 {
-  if (datatype == MPI_INT)
+  int n;
-    mpi_copy_int(sendbuf,recvbuf,count);
+  if (datatype == MPI_INT) n = count*sizeof(int);
-  else if (datatype == MPI_FLOAT)
+  else if (datatype == MPI_FLOAT) n = count*sizeof(float);
-    mpi_copy_float(sendbuf,recvbuf,count);
+  else if (datatype == MPI_DOUBLE) n = count*sizeof(double);
-  else if (datatype == MPI_DOUBLE)
+  else if (datatype == MPI_CHAR) n = count*sizeof(char);
-    mpi_copy_double(sendbuf,recvbuf,count);
+  else if (datatype == MPI_BYTE) n = count*sizeof(char);
-  else if (datatype == MPI_CHAR)
+  else if (datatype == MPI_DOUBLE_INT) n = count*sizeof(double_int);
-    mpi_copy_char(sendbuf,recvbuf,count);
+
-  else if (datatype == MPI_BYTE)
+  memcpy(recvbuf,sendbuf,n);
-    mpi_copy_byte(sendbuf,recvbuf,count);
 }
 
+/* ---------------------------------------------------------------------- */
+
 void MPI_Scan(void *sendbuf, void *recvbuf, int count,
 	      MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
 {
-  if (datatype == MPI_INT)
+  int n;
-    mpi_copy_int(sendbuf,recvbuf,count);
+  if (datatype == MPI_INT) n = count*sizeof(int);
-  else if (datatype == MPI_FLOAT)
+  else if (datatype == MPI_FLOAT) n = count*sizeof(float);
-    mpi_copy_float(sendbuf,recvbuf,count);
+  else if (datatype == MPI_DOUBLE) n = count*sizeof(double);
-  else if (datatype == MPI_DOUBLE)
+  else if (datatype == MPI_CHAR) n = count*sizeof(char);
-    mpi_copy_double(sendbuf,recvbuf,count);
+  else if (datatype == MPI_BYTE) n = count*sizeof(char);
-  else if (datatype == MPI_CHAR)
+  else if (datatype == MPI_DOUBLE_INT) n = count*sizeof(double_int);
-    mpi_copy_char(sendbuf,recvbuf,count);
+
-  else if (datatype == MPI_BYTE)
+  memcpy(recvbuf,sendbuf,n);
-    mpi_copy_byte(sendbuf,recvbuf,count);
 }
 
+/* ---------------------------------------------------------------------- */
+
 /* copy values from data1 to data2 */
 
 void MPI_Allgather(void *sendbuf, int sendcount, MPI_Datatype sendtype,
 		   void *recvbuf, int recvcount, MPI_Datatype recvtype,
 		   MPI_Comm comm)
 {
-  if (sendtype == MPI_INT)
+  int n;
-    mpi_copy_int(sendbuf,recvbuf,sendcount);
+  if (sendtype == MPI_INT) n = sendcount*sizeof(int);
-  else if (sendtype == MPI_FLOAT)
+  else if (sendtype == MPI_FLOAT) n = sendcount*sizeof(float);
-    mpi_copy_float(sendbuf,recvbuf,sendcount);
+  else if (sendtype == MPI_DOUBLE) n = sendcount*sizeof(double);
-  else if (sendtype == MPI_DOUBLE)
+  else if (sendtype == MPI_CHAR) n = sendcount*sizeof(char);
-    mpi_copy_double(sendbuf,recvbuf,sendcount);
+  else if (sendtype == MPI_BYTE) n = sendcount*sizeof(char);
-  else if (sendtype == MPI_CHAR)
+  else if (sendtype == MPI_DOUBLE_INT) n = sendcount*sizeof(double_int);
-    mpi_copy_char(sendbuf,recvbuf,sendcount);
+
-  else if (sendtype == MPI_BYTE)
+  memcpy(recvbuf,sendbuf,n);
-    mpi_copy_byte(sendbuf,recvbuf,sendcount);
 }
 
+/* ---------------------------------------------------------------------- */
+
 /* copy values from data1 to data2 */
 
 void MPI_Allgatherv(void *sendbuf, int sendcount, MPI_Datatype sendtype,
 		    void *recvbuf, int *recvcounts, int *displs,
 		    MPI_Datatype recvtype, MPI_Comm comm)
 {
-  if (sendtype == MPI_INT)
+  int n;
-    mpi_copy_int(sendbuf,recvbuf,sendcount);
+  if (sendtype == MPI_INT) n = sendcount*sizeof(int);
-  else if (sendtype == MPI_FLOAT)
+  else if (sendtype == MPI_FLOAT) n = sendcount*sizeof(float);
-    mpi_copy_float(sendbuf,recvbuf,sendcount);
+  else if (sendtype == MPI_DOUBLE) n = sendcount*sizeof(double);
-  else if (sendtype == MPI_DOUBLE)
+  else if (sendtype == MPI_CHAR) n = sendcount*sizeof(char);
-    mpi_copy_double(sendbuf,recvbuf,sendcount);
+  else if (sendtype == MPI_BYTE) n = sendcount*sizeof(char);
-  else if (sendtype == MPI_CHAR)
+  else if (sendtype == MPI_DOUBLE_INT) n = sendcount*sizeof(double_int);
-    mpi_copy_char(sendbuf,recvbuf,sendcount);
+
-  else if (sendtype == MPI_BYTE)
+  memcpy(recvbuf,sendbuf,n);
-    mpi_copy_byte(sendbuf,recvbuf,sendcount);
 }
 
+/* ---------------------------------------------------------------------- */
+
 /* copy values from data1 to data2 */
 
 void MPI_Reduce_scatter(void *sendbuf, void *recvbuf, int *recvcounts,
 			MPI_Datatype datatype, MPI_Op op, MPI_Comm comm)
 {
-  if (datatype == MPI_INT)
+  int n;
-    mpi_copy_int(sendbuf,recvbuf,*recvcounts);
+  if (datatype == MPI_INT) n = *recvcounts*sizeof(int);
-  else if (datatype == MPI_FLOAT)
+  else if (datatype == MPI_FLOAT) n = *recvcounts*sizeof(float);
-    mpi_copy_float(sendbuf,recvbuf,*recvcounts);
+  else if (datatype == MPI_DOUBLE) n = *recvcounts*sizeof(double);
-  else if (datatype == MPI_DOUBLE)
+  else if (datatype == MPI_CHAR) n = *recvcounts*sizeof(char);
-    mpi_copy_double(sendbuf,recvbuf,*recvcounts);
+  else if (datatype == MPI_BYTE) n = *recvcounts*sizeof(char);
-  else if (datatype == MPI_CHAR)
+  else if (datatype == MPI_DOUBLE_INT) n = *recvcounts*sizeof(double_int);
-    mpi_copy_char(sendbuf,recvbuf,*recvcounts);
+
-  else if (datatype == MPI_BYTE)
+  memcpy(recvbuf,sendbuf,n);
-    mpi_copy_byte(sendbuf,recvbuf,*recvcounts);
 }
 
+/* ---------------------------------------------------------------------- */
+
 /* copy values from data1 to data2 */
 
 void MPI_Gather(void *sendbuf, int sendcount, MPI_Datatype sendtype,
 		   void *recvbuf, int recvcount, MPI_Datatype recvtype,
 		   int root, MPI_Comm comm)
 {
-  if (sendtype == MPI_INT)
+  int n;
-    mpi_copy_int(sendbuf,recvbuf,sendcount);
+  if (sendtype == MPI_INT) n = sendcount*sizeof(int);
-  else if (sendtype == MPI_FLOAT)
+  else if (sendtype == MPI_FLOAT) n = sendcount*sizeof(float);
-    mpi_copy_float(sendbuf,recvbuf,sendcount);
+  else if (sendtype == MPI_DOUBLE) n = sendcount*sizeof(double);
-  else if (sendtype == MPI_DOUBLE)
+  else if (sendtype == MPI_CHAR) n = sendcount*sizeof(char);
-    mpi_copy_double(sendbuf,recvbuf,sendcount);
+  else if (sendtype == MPI_BYTE) n = sendcount*sizeof(char);
-  else if (sendtype == MPI_CHAR)
+  else if (sendtype == MPI_DOUBLE_INT) n = sendcount*sizeof(double_int);
-    mpi_copy_char(sendbuf,recvbuf,sendcount);
+
-  else if (sendtype == MPI_BYTE)
+  memcpy(recvbuf,sendbuf,n);
-    mpi_copy_byte(sendbuf,recvbuf,sendcount);
-}
-
-/* ------------------------------------------------------------------------ */
-/* Added routines for data copying */
-
-void mpi_copy_int(void *data1, void *data2, int n)
-{
-  int i;
-  int *pdata1 = (int *) data1;
-  int *pdata2 = (int *) data2;
-
-  for (i = 0; i < n; i++) pdata2[i] = pdata1[i];
-}
-
-void mpi_copy_float(void *data1, void *data2, int n)
-{
-  int i;
-  float *pdata1 = (float *) data1;
-  float *pdata2 = (float *) data2;
-
-  for (i = 0; i < n; i++) pdata2[i] = pdata1[i];
-}
-
-void mpi_copy_double(void *data1, void *data2, int n)
-{
-  int i;
-  double *pdata1 = (double *) data1;
-  double *pdata2 = (double *) data2;
-
-  for (i = 0; i < n; i++) pdata2[i] = pdata1[i];
-}
-
-void mpi_copy_char(void *data1, void *data2, int n)
-{
-  int i;
-  char *pdata1 = (char *) data1;
-  char *pdata2 = (char *) data2;
-
-  for (i = 0; i < n; i++) pdata2[i] = pdata1[i];
-}
-
-void mpi_copy_byte(void *data1, void *data2, int n)
-{
-  int i;
-  char *pdata1 = (char *) data1;
-  char *pdata2 = (char *) data2;
-
-  for (i = 0; i < n; i++) pdata2[i] = pdata1[i];
 }

src/STUBS/mpi.h

@@ -23,10 +23,13 @@
 #define MPI_DOUBLE 3
 #define MPI_CHAR 4
 #define MPI_BYTE 5
+#define MPI_DOUBLE_INT 6
 
 #define MPI_SUM 1
 #define MPI_MAX 2
 #define MPI_MIN 3
+#define MPI_MAXLOC 4
+#define MPI_MINLOC 5
 
 #define MPI_ANY_SOURCE -1
 

src/fix_viscosity.cpp

@@ -0,0 +1,208 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#include "mpi.h"
+#include "string.h"
+#include "stdlib.h"
+#include "fix_viscosity.h"
+#include "atom.h"
+#include "domain.h"
+#include "error.h"
+
+#include "update.h"
+
+using namespace LAMMPS_NS;
+
+#define BIG 1.0e20
+
+/* ---------------------------------------------------------------------- */
+
+FixViscosity::FixViscosity(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg != 7) error->all("Illegal fix viscosity command");
+
+  MPI_Comm_rank(world,&me);
+
+  nevery = atoi(arg[3]);
+  if (nevery <= 0) error->all("Illegal fix viscosity command");
+
+  scalar_flag = 1;
+  scalar_vector_freq = nevery;
+  extensive = 0;
+
+  if (strcmp(arg[4],"x") == 0) vdim = 0;
+  else if (strcmp(arg[4],"y") == 0) vdim = 1;
+  else if (strcmp(arg[4],"z") == 0) vdim = 2;
+  else error->all("Illegal fix viscosity command");
+
+  if (strcmp(arg[5],"x") == 0) pdim = 0;
+  else if (strcmp(arg[5],"y") == 0) pdim = 1;
+  else if (strcmp(arg[5],"z") == 0) pdim = 2;
+  else error->all("Illegal fix viscosity command");
+
+  nbin = atoi(arg[6]);
+  if (nbin < 3) error->all("Illegal fix viscosity command");
+
+  flux = 0.0;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixViscosity::setmask()
+{
+  int mask = 0;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixViscosity::init()
+{
+  // set bounds of 2 slabs in pdim
+  // only necessary for static box, else re-computed in end_of_step()
+  // lo bin is always bottom bin
+  // if nbin even, hi bin is just below half height
+  // if nbin odd, hi bin straddles half height
+
+  if (domain->box_change == 0) {
+    prd = domain->prd[pdim];
+    boxlo = domain->boxlo[pdim];
+    boxhi = domain->boxhi[pdim];
+    double binsize = (boxhi-boxlo) / nbin;
+    slablo_lo = boxlo;
+    slablo_hi = boxlo + binsize;
+    slabhi_lo = boxlo + ((nbin-1)/2)*binsize;
+    slabhi_hi = boxlo + ((nbin-1)/2 + 1)*binsize;
+  }
+
+  periodicity = domain->periodicity[pdim];
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixViscosity::end_of_step()
+{
+  int i;
+  double p,coord;
+  MPI_Status status;
+  struct {
+    double value;
+    int proc;
+  } mine[2],all[2];
+
+  // if box changes, recompute bounds of 2 slabs in pdim
+
+  if (domain->box_change) {
+    prd = domain->prd[pdim];
+    boxlo = domain->boxlo[pdim];
+    boxhi = domain->boxhi[pdim];
+    double binsize = (boxhi-boxlo) / nbin;
+    slablo_lo = boxlo;
+    slablo_hi = boxlo + binsize;
+    slabhi_lo = boxlo + ((nbin-1)/2)*binsize;
+    slabhi_hi = boxlo + ((nbin-1)/2 + 1)*binsize;
+  }
+
+  // ipos,ineg = my 2 atoms with most pos/neg momenta in bottom/top slabs
+  // map atom back into periodic box if necessary
+
+  double **x = atom->x;
+  double **v = atom->v;
+  double *mass = atom->mass;
+  double *rmass = atom->rmass;
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  int ipos = -1;
+  int ineg = -1;
+  double pmin = BIG;
+  double pmax = -BIG;
+
+  for (i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) {
+      if (mass) p = mass[type[i]] * v[i][vdim];
+      else p = rmass[i] * v[i][vdim];
+      coord = x[i][pdim];
+      if (coord < boxlo && periodicity) coord += prd;
+      else if (coord >= boxhi && periodicity) coord -= prd;
+      if (coord >= slablo_lo && coord < slablo_hi) {
+	if (p > pmax) {
+	  pmax = p;
+	  ipos = i;
+	}
+      }
+      if (coord >= slabhi_lo && coord < slabhi_hi) {
+	if (p < pmin) {
+	  pmin = p;
+	  ineg = i;
+	}
+      }
+    }
+
+  // find 2 global atoms with most pos/neg momenta in bottom/top slabs
+  // MAXLOC also communicates which procs own them
+
+  mine[0].value = pmax;
+  mine[0].proc = me;
+  mine[1].value = -pmin;
+  mine[1].proc = me;
+
+  MPI_Allreduce(mine,all,2,MPI_DOUBLE_INT,MPI_MAXLOC,world);
+
+  // exchange momenta between the 2 particles
+  // if I own both particles just swap, else point2point comm of mass,vel
+  
+  double sbuf[2],rbuf[2];
+
+  if (me == all[0].proc && me == all[1].proc) {
+    sbuf[0] = v[ineg][vdim];
+    if (mass) sbuf[1] = mass[type[ineg]];
+    else sbuf[1] = rmass[ineg];
+    rbuf[0] = v[ipos][vdim];
+    if (mass) rbuf[1] = mass[type[ipos]];
+    else rbuf[1] = rmass[ipos];
+    v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1];
+    v[ipos][vdim] = sbuf[0] * sbuf[1]/rbuf[1];
+
+  } else if (me == all[0].proc) {
+    sbuf[0] = v[ipos][vdim];
+    if (mass) sbuf[1] = mass[type[ipos]];
+    else sbuf[1] = rmass[ipos];
+    MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[1].proc,0,
+		 rbuf,2,MPI_DOUBLE,all[1].proc,0,world,&status);
+    v[ipos][vdim] = rbuf[0] * rbuf[1]/sbuf[1];
+
+  } else if (me == all[1].proc) {
+    sbuf[0] = v[ineg][vdim];
+    if (mass) sbuf[1] = mass[type[ineg]];
+    else sbuf[1] = rmass[ineg];
+    MPI_Sendrecv(sbuf,2,MPI_DOUBLE,all[0].proc,0,
+		 rbuf,2,MPI_DOUBLE,all[0].proc,0,world,&status);
+    v[ineg][vdim] = rbuf[0] * rbuf[1]/sbuf[1];
+  }
+
+  // tally momentum flux
+  // sign of all[1].value was flipped for MPI_Allreduce
+
+  flux += all[0].value + all[1].value;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixViscosity::compute_scalar()
+{
+  return flux;
+}

src/fix_viscosity.h

@@ -0,0 +1,40 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_VISCOSITY_H
+#define FIX_VISCOSITY_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixViscosity : public Fix {
+ public:
+  FixViscosity(class LAMMPS *, int, char **);
+  ~FixViscosity() {}
+  int setmask();
+  void init();
+  void end_of_step();
+  double compute_scalar();
+
+ private:
+  int me;
+  int vdim,pdim,nbin,periodicity;
+  double prd,boxlo,boxhi;
+  double slablo_lo,slablo_hi,slabhi_lo,slabhi_hi;
+  double flux;
+};
+
+}
+
+#endif

src/lammps.cpp

@@ -306,7 +306,7 @@ void LAMMPS::init()
                          //   atom deletes extra array
                          //   used by fix shear_history::unpack_restart()
                          //   when force->pair->gran_history creates fix ??
-  modify->init();        // modify must come after update, force, atom
+  modify->init();        // modify must come after update, force, atom, domain
   neighbor->init();      // neighbor must come after force, modify
   output->init();        // output must come after domain, force, modify
   comm->init();          // comm must come after force, modify

src/style.h

@@ -181,6 +181,7 @@ DumpStyle(xyz,DumpXYZ)
 #include "fix_spring_self.h"
 #include "fix_temp_rescale.h"
 #include "fix_tmd.h"
+#include "fix_viscosity.h"
 #include "fix_viscous.h"
 #include "fix_wall_lj126.h"
 #include "fix_wall_lj93.h"
@@ -232,6 +233,7 @@ FixStyle(spring/rg,FixSpringRG)
 FixStyle(spring/self,FixSpringSelf)
 FixStyle(temp/rescale,FixTempRescale)
 FixStyle(tmd,FixTMD)
+FixStyle(viscosity,FixViscosity)
 FixStyle(viscous,FixViscous)
 FixStyle(wall/lj126,FixWallLJ126)
 FixStyle(wall/lj93,FixWallLJ93)

src/thermo.cpp

@@ -228,7 +228,8 @@ void Thermo::init()
     if (format_user[i]) ptr = format_user[i];
     else if (vtype[i] == INT && format_int_user) ptr = format_int_user;
     else if (vtype[i] == INT && lineflag == ONELINE) ptr = format_int_one_def;
-    else if (vtype[i] == INT && lineflag == MULTILINE) ptr = format_int_multi_def;
+    else if (vtype[i] == INT && lineflag == MULTILINE) 
+      ptr = format_int_multi_def;
     else if (vtype[i] == FLOAT && format_float_user) ptr = format_float_user;
     else if (lineflag == ONELINE) ptr = format_g_def;
     else if (lineflag == MULTILINE) ptr = format_f_def;