From 24f7d34ec3b90ab6e369c40d183535b6bfc7d179 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Fri, 22 Jun 2007 17:23:42 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@644
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 tools/Makefile    |   7 +-
 tools/replicate.c | 861 ----------------------------------------------
 2 files changed, 2 insertions(+), 866 deletions(-)
 delete mode 100644 tools/replicate.c

diff --git a/tools/Makefile b/tools/Makefile
index 6e4c93e925..a9ec59bf66 100644
--- a/tools/Makefile
+++ b/tools/Makefile
@@ -5,14 +5,11 @@
 #
 
 all:
-	$(MAKE) binary2txt replicate restart2data chain micelle2d data2xmovie
+	$(MAKE) binary2txt restart2data chain micelle2d data2xmovie
 
 binary2txt:	binary2txt.o
 	g++ -g binary2txt.o -o binary2txt
 
-replicate:	replicate.o
-	g++ -g replicate.o -o replicate
-
 restart2data:	restart2data.o
 	g++ -g restart2data.o -o restart2data
 
@@ -29,7 +26,7 @@ thermo_extract:	thermo_extract.o
 	gcc -g thermo_extract.o -o thermo_extract
 
 clean:
-	rm binary2txt replicate restart2data chain micelle2d data2xmovie
+	rm binary2txt restart2data chain micelle2d data2xmovie
 	rm thermo_extract
 	rm *.o
 
diff --git a/tools/replicate.c b/tools/replicate.c
deleted file mode 100644
index ec7f89a9e6..0000000000
--- a/tools/replicate.c
+++ /dev/null
@@ -1,861 +0,0 @@
-/* replicate tool
-
-   read LAMMPS data file as input
-   create new LAMMPS data file as output
-
-   new file can be translated, dilated, mapped to periodic box
-     and/or replicated to create a larger system
-   create new molecule tags as needed
-
-   Syntax:  replicate [options] < infile > outfile
-
-     Options:
-
-     -style atom_style
-
-          use the LAMMPS atom style that corresponds to this file
-	    e.g. atomic or bond or angle or full or granular
-          will be used for reading and writing files
-	  this option must be specified
-	  for atom_style = hybrid, append subsequent styles
-	    e.g. -style hybrid bond charge
-
-     -repeat nx ny nz
-
-          use the input data as a unit cell and replicate the unit
-	    cell by (nx,ny,nz) in the +x, +y, +z directions
-	  this creates new atoms, bonds, angles, etc
-	  if not specified, nx = ny = nz = 1
-
-     -inbox xlo xhi ylo yhi zlo zhi
-
-          use these values for the bounding box of the input system
-	    and for the unit cell used for replicating purposes
-	  if not specified use the box bounds read in from infile
-
-     -outbox xlo xhi ylo yhi zlo zhi
-
-          use these values for the bounding box of the output system
-	  if not specified the output bounding box is computed from the input
-	    bounding box by holding the "lo" values constant and creating 
-	    new "hi" values by multiplying by the replicating factors
-
-     -unmap
-
-          unmap all input atom positions using input image flags
-            and inbox in a periodic sense
-	  this enables molecular bonds to be replicated
-
-     -map
-
-          map all final atom positions into the output bounding box
-	    in a periodic sense
-	  works no matter how far away from the box the atoms are
-	  this operation is independent of unmap option
-
-     -imageout
-
-          image flags will be written to outfile, else they won't be
-	  if unmapping/mapping is done, image flags are modified accordingly
-
-   Notes:
-      
-      works with "Velocities" entry by giving same velocity to
-        all replicated images of an atom, which is probably not ideal
-
-      can perform system dilation/contraction by using "-repeat 1 1 1"
-        and specifying a slightly larger/smaller output box than input
-        box.
-
-      can add/delete image flags by using -imageout
-*/
-
-#include <string.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-
-struct box {
-  double xlo,xhi,xsize;
-  double ylo,yhi,ysize;
-  double zlo,zhi,zsize;
-};
-
-#define MAXLINE 1000
-
-main(int argc, char *argv[])
-
-{
-  char line[MAXLINE];        /* strings for reading/parsing input file */
-  char linetoken[MAXLINE];
-  char *token,*ptr;             
-                             /* numbers of various quantities as read in */
-
-  int natoms,nbonds,nangles,ndihedrals,nimpropers;
-  int ntypes,nbondtypes,nangletypes,ndihedtypes,nimprotypes;
-
-  int nmolecules;            /* total # of mols specified by mol tags */
-
-                             /* vectors of various quantities as read in */
-
-  int *tag,*molecule,*type;
-  double *q,*radius,*density,*x,*y,*z,*mux,*muy,*muz;
-  int *imagex,*imagey,*imagez;
-  int *vtag;
-  double *vx,*vy,*vz,*vphix,*vphiy,*vphiz;
-  int *btype,*atype,*dtype,*itype;
-  int *bond1,*bond2,*angle1,*angle2,*angle3;
-  int *dihed1,*dihed2,*dihed3,*dihed4,*impro1,*impro2,*impro3,*impro4;
-
-  int atomstyle;             /* atom styles and settings */
-  int style_angle,style_atomic,style_bond,style_charge,style_dipole;
-  int style_dpd,style_full,style_granular,style_molecular;
-  int charge,molecular,dipole;
-
-  int nx,ny,nz;              /* replication factors */
-  int inbox,outbox;          /* flags for whether in/out boxes were explicitly
-				specified */
-  int unmap;                 /* 0 = no unmapping, 1 = unmap from input box */
-  int map;                   /* 0 = no mapping, 1 = map to output box */
-  int imageout;              /* 0/1 = no/yes output of image flags */
-  int imageflag;             /* 0/1 = no/yes if image flags exist in input */
-
-  struct box in,out;         /* input and output bounding boxes */
-  int ntotal;                /* total replication factor = nx*ny*nz */
-  double newx,newy,newz;     /* replicated atom position */
-  int size_atom_valid;       /* values that should be on atom line */
-  int size_atom_actual;      /* values that are on atom line */
-
-  int i,j,k,m,n,tmp,del,tag_offset,mol_offset;   /* temp variables */
-  char **values;
-
-  char *gettoken(char *, int);    /* function defs */
-  int count_words(char *); 
-  void copyline(int); 
-  void scale(struct box, int, int, int, 
-	     struct box, double *, double *, double *);
-  void pbc(struct box, double *, double *, double *, int *, int *, int *);
-
-  /* default input values */
-
-  atomstyle = 0;
-  style_angle = style_atomic = style_bond = style_charge = style_dipole = 
-    style_dpd = style_full = style_granular = style_molecular = 0;
-  nx = ny = nz = 1;
-  inbox = outbox = 0;
-  unmap = 0;
-  map = 0;
-  imageout = 0;
-  
-  /* read input options from command line
-     should do more error checking for missing args */
-
-  i = 1;
-  while (i < argc) {
-
-    if (!strcmp(argv[i],"-style")) {
-      atomstyle = 1;
-      if (strcmp(argv[i+1],"angle") == 0) style_angle = 1;
-      else if (strcmp(argv[i+1],"atomic") == 0) style_atomic = 1;
-      else if (strcmp(argv[i+1],"bond") == 0) style_bond = 1;
-      else if (strcmp(argv[i+1],"charge") == 0) style_charge = 1;
-      else if (strcmp(argv[i+1],"dipole") == 0) style_dipole = 1;
-      else if (strcmp(argv[i+1],"dpd") == 0) style_dpd = 1;
-      else if (strcmp(argv[i+1],"full") == 0) style_full = 1;
-      else if (strcmp(argv[i+1],"granular") == 0) style_granular = 1;
-      else if (strcmp(argv[i+1],"molecular") == 0) style_molecular = 1;
-      else if (strcmp(argv[i+1],"hybrid") == 0) {
-	n = i + 2;
-	while (n < argc && argv[n][0] != '-') {
-	  if (strcmp(argv[n],"angle") == 0) style_angle = 1;
-	  else if (strcmp(argv[n],"atomic") == 0) style_atomic = 1;
-	  else if (strcmp(argv[n],"bond") == 0) style_bond = 1;
-	  else if (strcmp(argv[n],"charge") == 0) style_charge = 1;
-	  else if (strcmp(argv[n],"dipole") == 0) style_dipole = 1;
-	  else if (strcmp(argv[n],"dpd") == 0) style_dpd = 1;
-	  else if (strcmp(argv[n],"full") == 0) style_full = 1;
-	  else if (strcmp(argv[n],"granular") == 0) style_granular = 1;
-	  else if (strcmp(argv[n],"molecular") == 0) style_molecular = 1;
-	  else {
-	    fprintf(stderr,"Error with command-line arg style\n");
-	    exit(1);
-	  }
-	  n++;
-	}
-	if (i == n - 2) {
-	  fprintf(stderr,"Error with command-line arg style\n");
-	  exit(1);
-	}
-	i = n - 2;
-      } else {
-	fprintf(stderr,"Error with command-line arg style\n");
-	exit(1);
-      }
-      i += 2;
-    } else if (!strcmp(argv[i],"-repeat")) {
-      sscanf(argv[i+1],"%d",&nx);
-      sscanf(argv[i+2],"%d",&ny);
-      sscanf(argv[i+3],"%d",&nz);
-      i += 4;
-    } else if (!strcmp(argv[i],"-inbox")) {
-      inbox = 1;
-      sscanf(argv[i+1],"%lg",&in.xlo);
-      sscanf(argv[i+2],"%lg",&in.xhi);
-      sscanf(argv[i+3],"%lg",&in.ylo);
-      sscanf(argv[i+4],"%lg",&in.yhi);
-      sscanf(argv[i+5],"%lg",&in.zlo);
-      sscanf(argv[i+6],"%lg",&in.zhi);
-      i += 7;
-    } else if (!strcmp(argv[i],"-outbox")) {
-      outbox = 1;
-      sscanf(argv[i+1],"%lg",&out.xlo);
-      sscanf(argv[i+2],"%lg",&out.xhi);
-      sscanf(argv[i+3],"%lg",&out.ylo);
-      sscanf(argv[i+4],"%lg",&out.yhi);
-      sscanf(argv[i+5],"%lg",&out.zlo);
-      sscanf(argv[i+6],"%lg",&out.zhi);
-      i += 7;
-    } else if (!strcmp(argv[i],"-unmap")) {
-      unmap = 1;
-      i += 1;
-    } else if (!strcmp(argv[i],"-map")) {
-      map = 1;
-      i += 1;
-    } else if (!strcmp(argv[i],"-imageout")) {
-      imageout = 1;
-      i += 1;
-    } else {
-      fprintf(stderr,"Syntax error: replicate [options] < infile > outfile\n");
-      exit(1);
-    }
-  }
-
-  /* set low-level info from style flags */
-
-  charge = 0;
-  if (style_charge || style_full || style_dipole) charge = 1;
-
-  molecular = 0;
-  if (style_angle || style_bond || style_full || style_molecular)
-    molecular = 1;
-
-  dipole = 0;
-  if (style_dipole) dipole = 1;
-
-  size_atom_valid = 5;
-  if (charge) size_atom_valid += 1;
-  if (style_dipole) size_atom_valid += 3;
-  if (style_granular) size_atom_valid += 2;
-  if (molecular) size_atom_valid += 1;
-
-  /* error checks */
-
-  if (atomstyle == 0) {
-    fprintf(stderr,"ERROR: must use -style to set atom style\n");
-    exit(1);
-  }
-
-  /* ntotal = total replication factor */
-
-  ntotal = nx*ny*nz;
-
-  /* defaults */
-
-  natoms = nbonds = nangles = ndihedrals = nimpropers = 0;
-  ntypes = nbondtypes = nangletypes = ndihedtypes = nimprotypes = 0;
-
-  /* read/write header */
-
-  copyline(1);   // echo 1st line
-
-  while (1) {
-
-    fgets(line,MAXLINE,stdin);
-
-    // trim anything from '#' onward
-    // if line is blank, continue
-
-    if (ptr = strchr(line,'#')) *ptr = '\0';
-    if (strspn(line," \t\n\r") == strlen(line)) {
-      printf("\n");
-      continue;
-    }
-
-    if (strstr(line,"atoms")) {
-      sscanf(line,"%d",&natoms);
-      printf("%d atoms\n",natoms*ntotal);
-    } else if (strstr(line,"bonds")) {
-      sscanf(line,"%d",&nbonds);
-      printf("%d bonds\n",nbonds*ntotal);
-    } else if (strstr(line,"angles")) {
-      sscanf(line,"%d",&nangles);
-      printf("%d angles\n",nangles*ntotal);
-    } else if (strstr(line,"dihedrals")) {
-      sscanf(line,"%d",&ndihedrals);
-      printf("%d dihedrals\n",ndihedrals*ntotal);
-    } else if (strstr(line,"impropers")) {
-      sscanf(line,"%d",&nimpropers);
-      printf("%d impropers\n",nimpropers*ntotal);
-    } else if (strstr(line,"atom types")) {
-      sscanf(line,"%d",&ntypes);
-      printf("%d atom types\n",ntypes);
-    } else if (strstr(line,"bond types")) {
-      sscanf(line,"%d",&nbondtypes);
-      printf("%d bond types\n",nbondtypes);
-    } else if (strstr(line,"angle types")) {
-      sscanf(line,"%d",&nangletypes);
-      printf("%d angle types\n",nangletypes);
-    } else if (strstr(line,"dihedral types")) {
-      sscanf(line,"%d",&ndihedtypes);
-      printf("%d dihedral types\n",ndihedtypes);
-    } else if (strstr(line,"improper types")) {
-      sscanf(line,"%d",&nimprotypes);
-      printf("%d improper types\n",nimprotypes);
-    } else if (strstr(line,"xlo xhi")) {
-      if (!inbox) sscanf(line,"%lg %lg",&in.xlo,&in.xhi);
-      in.xsize = in.xhi - in.xlo;
-      if (!outbox) {
-	out.xlo = in.xlo;
-	if (nx > 1) out.xhi = out.xlo + (in.xhi-in.xlo)*nx;
-	else out.xhi = in.xhi;
-      }
-      out.xsize = out.xhi - out.xlo;
-      printf("%g %g xlo xhi\n",out.xlo,out.xhi);
-    } else if (strstr(line,"ylo yhi")) {
-      if (!inbox) sscanf(line,"%lg %lg",&in.ylo,&in.yhi);
-      in.ysize = in.yhi - in.ylo;
-      if (!outbox) {
-	out.ylo = in.ylo;
-	if (ny > 1) out.yhi = out.ylo + (in.yhi-in.ylo)*ny;
-	else out.yhi = in.yhi;
-      }
-      out.ysize = out.yhi - out.ylo;
-      printf("%g %g ylo yhi\n",out.ylo,out.yhi);
-    } else if (strstr(line,"zlo zhi")) {
-      if (!inbox) sscanf(line,"%lg %lg",&in.zlo,&in.zhi);
-      in.zsize = in.zhi - in.zlo;
-      if (!outbox) {
-	out.zlo = in.zlo;
-	if (nz > 1) out.zhi = out.zlo + (in.zhi-in.zlo)*nz;
-	else out.zhi = in.zhi;
-      }
-      out.zsize = out.zhi - out.zlo;
-      printf("%g %g zlo zhi\n",out.zlo,out.zhi);
-    } else break;
-  }
-
-  /* malloc space for input atoms and topology */
-
-  tag = (int *) malloc(natoms*sizeof(int));
-  type = (int *) malloc(natoms*sizeof(int));
-  x = (double *) malloc(natoms*sizeof(double));
-  y = (double *) malloc(natoms*sizeof(double));
-  z = (double *) malloc(natoms*sizeof(double));
-  imagex = (int *) malloc(natoms*sizeof(int));
-  imagey = (int *) malloc(natoms*sizeof(int));
-  imagez = (int *) malloc(natoms*sizeof(int));
-
-  vtag = (int *) malloc(natoms*sizeof(int));
-  vx = (double *) malloc(natoms*sizeof(double));
-  vy = (double *) malloc(natoms*sizeof(double));
-  vz = (double *) malloc(natoms*sizeof(double));
-
-  if (tag == NULL || type == NULL ||
-      x == NULL || y == NULL || z == NULL ||
-      imagex == NULL || imagey == NULL || imagez == NULL ||
-      vtag == NULL || vx == NULL || vy == NULL || vz == NULL) {
-    fprintf(stderr,"Error in atom malloc - no space\n");
-    exit(1);
-  }
-
-  if (charge) {
-    q = (double *) malloc(natoms*sizeof(double));
-    if (q == NULL) {
-      fprintf(stderr,"Error in atom malloc - no space\n");
-      exit(1);
-    }
-  }
-
-  if (molecular) {
-    molecule = (int *) malloc(natoms*sizeof(int));
-    if (molecule == NULL) {
-      fprintf(stderr,"Error in atom malloc - no space\n");
-      exit(1);
-    }
-  }
-
-  if (dipole) {
-    mux = (double *) malloc(natoms*sizeof(double));
-    muy = (double *) malloc(natoms*sizeof(double));
-    muz = (double *) malloc(natoms*sizeof(double));
-    if (mux == NULL || muy == NULL || muz == NULL) {
-      fprintf(stderr,"Error in atom malloc - no space\n");
-      exit(1);
-    }
-  }
-
-  if (style_granular) {
-    radius = (double *) malloc(natoms*sizeof(double));
-    density = (double *) malloc(natoms*sizeof(double));
-    vphix = (double *) malloc(natoms*sizeof(double));
-    vphiy = (double *) malloc(natoms*sizeof(double));
-    vphiz = (double *) malloc(natoms*sizeof(double));
-    if (radius == NULL || density == NULL || 
-	vphix == NULL || vphiy == NULL || vphiz == NULL) {
-      fprintf(stderr,"Error in atom malloc - no space\n");
-      exit(1);
-    }
-  }
-
-  if (nbonds) {
-    btype = (int *) malloc(nbonds*sizeof(int));
-    bond1 = (int *) malloc(nbonds*sizeof(int));
-    bond2 = (int *) malloc(nbonds*sizeof(int));
-    
-    if (btype == NULL || bond1 == NULL || bond2 == NULL) {
-      fprintf(stderr,"Error in bond malloc - no space\n");
-      exit(1);
-    }
-  }
-
-  if (nangles) {
-    atype = (int *) malloc(nangles*sizeof(int));
-    angle1 = (int *) malloc(nangles*sizeof(int));
-    angle2 = (int *) malloc(nangles*sizeof(int));
-    angle3 = (int *) malloc(nangles*sizeof(int));
-    
-    if (atype == NULL || angle1 == NULL || angle2 == NULL || angle3 == NULL) {
-      fprintf(stderr,"Error in angle malloc - no space\n");
-      exit(1);
-    }
-  }
-
-  if (ndihedrals) {
-    dtype = (int *) malloc(ndihedrals*sizeof(int));
-    dihed1 = (int *) malloc(ndihedrals*sizeof(int));
-    dihed2 = (int *) malloc(ndihedrals*sizeof(int));
-    dihed3 = (int *) malloc(ndihedrals*sizeof(int));
-    dihed4 = (int *) malloc(ndihedrals*sizeof(int));
-    
-    if (dtype == NULL || dihed1 == NULL || dihed2 == NULL || 
-	dihed3 == NULL || dihed4 == NULL) {
-      fprintf(stderr,"Error in dihedral malloc - no space\n");
-      exit(1);
-    }
-  }
-
-  if (nimpropers) {
-    itype = (int *) malloc(nimpropers*sizeof(int));
-    impro1 = (int *) malloc(nimpropers*sizeof(int));
-    impro2 = (int *) malloc(nimpropers*sizeof(int));
-    impro3 = (int *) malloc(nimpropers*sizeof(int));
-    impro4 = (int *) malloc(nimpropers*sizeof(int));
-    
-    if (itype == NULL || impro1 == NULL || impro2 == NULL || 
-	impro3 == NULL || impro4 == NULL) {
-      fprintf(stderr,"Error in improper malloc - no space\n");
-      exit(1);
-    }
-  }
-
-  /* read identifier strings one by one in free-form part of data file */
-
-  token = gettoken(line,1);
-
-  while (token) {
-
-    /* read atoms */
-
-    if (!strcmp(token,"Atoms")) {
-
-      for (i = 0; i < natoms; i++) {
-	fgets(line,MAXLINE,stdin);
-
-	if (i == 0) {
-	  size_atom_actual = count_words(line);
-	  if (size_atom_actual == size_atom_valid) imageflag = 0;
-	  else if (size_atom_actual == size_atom_valid + 3) imageflag = 1;
-	  else {
-	    fprintf(stderr,"Error in atom lines\n");
-	    exit(1);
-	  }
-	  if (unmap == 1 && imageflag == 0) {
-	    fprintf(stderr,
-		    "If unmap is requested, file must have image flags\n");
-	    exit(1);
-	  }
-	  values = (char **) malloc(size_atom_actual*sizeof(char *));
-	}
-
-	values[0] = strtok(line," \t\n");
-	for (m = 1; m < size_atom_actual; m++)
-	  values[m] = strtok(NULL," \t\n");
-
-	m = 0;
-	tag[i] = atoi(values[m++]);
-	if (molecular) molecule[i] = atoi(values[m++]);
-	type[i] = atoi(values[m++]);
-	if (charge) q[i] = atof(values[m++]);
-	if (style_granular) {
-	  radius[i] = 0.5 * atof(values[m++]);
-	  density[i] = atof(values[m++]);
-	}
-	x[i] = atof(values[m++]);
-	y[i] = atof(values[m++]);
-	z[i] = atof(values[m++]);
-	if (style_dipole) {
-	  mux[i] = atof(values[m++]);
-	  muy[i] = atof(values[m++]);
-	  muz[i] = atof(values[m++]);
-	}
-	if (imageflag) {
-	  imagex[i] = atoi(values[m++]);
-	  imagey[i] = atoi(values[m++]);
-	  imagez[i] = atoi(values[m++]);
-	} else imagex[i] = imagey[i] = imagez[i] = 0;
-      }
-      
-      /* unmap atoms if requested */
-
-      if (unmap) {
-	for (i = 0; i < natoms; i++) {
-	  x[i] = x[i] + imagex[i]*in.xsize;
-	  y[i] = y[i] + imagey[i]*in.ysize;
-	  z[i] = z[i] + imagez[i]*in.zsize;
-	  imagex[i] = imagey[i] = imagez[i] = 0;
-	}
-      }
-
-      /* find number of molecules */
-
-      nmolecules = 0;
-      if (molecular)
-	for (i = 0; i < natoms; i++)
-	  if (molecule[i] > nmolecules) nmolecules = molecule[i];
-
-      /* replicate set of N atoms as many times as requested
-	 generate new tags, mol number, coords, image flags as needed */
-
-      for (k = 0; k < nz; k++) {
-	for (j = 0; j < ny; j++) {
-	  for (i = 0; i < nx; i++) {
-	    for (m = 0; m < natoms; m++) {
-	      newx = x[m] + i*in.xsize;
-	      newy = y[m] + j*in.ysize;
-	      newz = z[m] + k*in.zsize;
-	      if (outbox) scale(in,nx,ny,nz,out,&newx,&newy,&newz);
-	      if (map) pbc(out,&newx,&newy,&newz,
-			   &imagex[m],&imagey[m],&imagez[m]);
-              tag_offset = k*ny*nx*natoms + j*nx*natoms + i*natoms;
-              mol_offset = k*ny*nx*nmolecules + j*nx*nmolecules + i*nmolecules;
-
-	      line[0] = '\0';
-	      sprintf(&line[strlen(line)],"%d",tag[m]+tag_offset);
-	      if (molecular) sprintf(&line[strlen(line)]," %d",
-				     molecule[m]+mol_offset);
-	      sprintf(&line[strlen(line)]," %d",type[m]);
-	      if (charge) sprintf(&line[strlen(line)]," %g",q[m]);
-	      if (style_granular) sprintf(&line[strlen(line)]," %g %g",
-					  radius[m],density[m]);
-	      sprintf(&line[strlen(line)]," %g %g %g",newx,newy,newz);
-	      if (style_dipole) sprintf(&line[strlen(line)]," %g %g %g",
-					mux[m],muy[m],muz[m]);
-	      if (imageout) sprintf(&line[strlen(line)]," %d %d %d",
-					imagex[m],imagey[m],imagez[m]);
-	      printf("%s\n",line);
-	    }
-	  }
-	}
-      }
-
-    }
-
-    /* read velocities and replicate */
-    
-    else if (!strcmp(token,"Velocities")) {
-
-      for (i = 0; i < natoms; i++) {
-	fgets(line,MAXLINE,stdin);
-	if (!style_granular) 
-	  sscanf(line,"%d %lg %lg %lg",&vtag[i],&vx[i],&vy[i],&vz[i]);
-	else 
-	  sscanf(line,"%d %lg %lg %lg %lg %lg %lg",&vtag[i],
-		 &vx[i],&vy[i],&vz[i],&vphix[i],&vphiy[i],&vphiz[i]);
-      }
-      
-      for (k = 0; k < nz; k++) {
-	for (j = 0; j < ny; j++) {
-	  for (i = 0; i < nx; i++) {
-	    for (m = 0; m < natoms; m++) {
-              tag_offset = k*ny*nx*natoms + j*nx*natoms + i*natoms;
-	      if (!style_granular) 
-		printf("%d %g %g %g\n",vtag[m]+tag_offset,vx[m],vy[m],vz[m]);
-	      else 
-		printf("%d %g %g %g %g %g %g\n",vtag[m]+tag_offset,
-		       vx[m],vy[m],vz[m],vphix[m],vphiy[m],vphiz[m]);
-	    }
-	  }
-	}
-      }
-
-    }
-
-    /* read bonds and replicate */
-
-    else if (!strcmp(token,"Bonds")) {
-
-      for (i = 0; i < nbonds; i++) {
-	fgets(line,MAXLINE,stdin);
-	sscanf(line,"%d %d %d %d",&n,&btype[i],&bond1[i],&bond2[i]);
-      }
-    
-      n = 0;
-      for (k = 0; k < nz; k++)
-	for (j = 0; j < ny; j++)
-	  for (i = 0; i < nx; i++)
-	    for (m = 0; m < nbonds; m++) {
-	      n++;
-	      del = k*ny*nx*natoms + j*nx*natoms + i*natoms;
-	      printf("%d %d %d %d\n",n,btype[m],bond1[m]+del,bond2[m]+del);
-	    }
-    }
-
-    /* read angles and replicate */
-
-    else if (!strcmp(token,"Angles")) {
-
-      for (i = 0; i < nangles; i++) {
-	fgets(line,MAXLINE,stdin);
-	sscanf(line,"%d %d %d %d %d",&n,&atype[i],
-	       &angle1[i],&angle2[i],&angle3[i]);
-      }
-    
-      n = 0;
-      for (k = 0; k < nz; k++)
-	for (j = 0; j < ny; j++)
-	  for (i = 0; i < nx; i++)
-	    for (m = 0; m < nangles; m++) {
-	      n++;
-	      del = k*ny*nx*natoms + j*nx*natoms + i*natoms;
-	      printf("%d %d %d %d %d\n",n,atype[m],
-		     angle1[m]+del,angle2[m]+del,angle3[m]+del);
-	    }
-    }
-    
-    /* read dihedrals and replicate */
-    
-    else if (!strcmp(token,"Dihedrals")) {
-
-      for (i = 0; i < ndihedrals; i++) {
-	fgets(line,MAXLINE,stdin);
-	sscanf(line,"%d %d %d %d %d %d",&n,&dtype[i],&dihed1[i],&dihed2[i],
-	       &dihed3[i],&dihed4[i]);
-      }
-
-      n = 0;
-      for (k = 0; k < nz; k++)
-	for (j = 0; j < ny; j++)
-	  for (i = 0; i < nx; i++)
-	    for (m = 0; m < ndihedrals; m++) {
-	      n++;
-	      del = k*ny*nx*natoms + j*nx*natoms + i*natoms;
-	      printf("%d %d %d %d %d %d\n",n,dtype[m],
-		     dihed1[m]+del,dihed2[m]+del,dihed3[m]+del,dihed4[m]+del);
-	    }
-    }
-
-    /* read impropers and replicate */
-
-    else if (!strcmp(token,"Impropers")) {
-
-      for (i = 0; i < nimpropers; i++) {
-	fgets(line,MAXLINE,stdin);
-	sscanf(line,"%d %d %d %d %d %d",&n,&itype[i],&impro1[i],&impro2[i],
-	     &impro3[i],&impro4[i]);
-      }
-
-      n = 0;
-      for (k = 0; k < nz; k++)
-	for (j = 0; j < ny; j++)
-	  for (i = 0; i < nx; i++)
-	    for (m = 0; m < nimpropers; m++) {
-	      n++;
-	      del = k*ny*nx*natoms + j*nx*natoms + i*natoms;
-	      printf("%d %d %d %d %d %d\n",n,itype[m],
-		     impro1[m]+del,impro2[m]+del,impro3[m]+del,impro4[m]+del);
-	    }
-    }
-    
-    /* non-replicated sections
-       just copy proper number of lines from in to out */
-
-    else if (!strcmp(token,"Masses"))
-      copyline(ntypes);
-    else if (!strcmp(token,"Dipoles"))
-      copyline(ntypes);
-    else if (!strcmp(token,"Pair Coeffs"))
-      copyline(ntypes);
-    else if (!strcmp(token,"Bond Coeffs"))
-      copyline(nbondtypes);
-    else if (!strcmp(token,"Angle Coeffs"))
-      copyline(nangletypes);
-    else if (!strcmp(token,"Dihedral Coeffs"))
-      copyline(ndihedtypes);
-    else if (!strcmp(token,"Improper Coeffs"))
-      copyline(nimprotypes);
-    else if (!strcmp(token,"BondBond Coeffs"))
-      copyline(nangletypes);
-    else if (!strcmp(token,"BondAngle Coeffs"))
-      copyline(nangletypes);
-    else if (!strcmp(token,"MiddleBondTorsion Coeffs"))
-      copyline(ndihedtypes);
-    else if (!strcmp(token,"EndBondTorsion Coeffs"))
-      copyline(ndihedtypes);
-    else if (!strcmp(token,"AngleTorsion Coeffs"))
-      copyline(ndihedtypes);
-    else if (!strcmp(token,"AngleAngleTorsion Coeffs"))
-      copyline(ndihedtypes);
-    else if (!strcmp(token,"BondBond13 Coeffs"))
-      copyline(ndihedtypes);
-    else if (!strcmp(token,"AngleAngle Coeffs"))
-      copyline(nimprotypes);
-    else {
-      fprintf(stderr,
-	      "Error in input data file - unknown identifier %s\n",token);
-      exit(1);
-    }
-
-    token = gettoken(line,0);
-  }
-  
-}
-
-/* ------------------------------------------------------------------- */
-
-/* return a LAMMPS keyword from data file
-   if first is 1, non-blank line with token is in line
-   else read until find a non-blank line
-   keyword is all text on line w/out leading & trailing white space
-   read one additional line after non-blank line (assumed blank)
-   return ptr to keyword
-   return NULL if end of file */
-
-char *gettoken(char *line, int first)
-{
-  char buffer[MAXLINE];
-
-  // read upto non-blank line plus 1 following line
-  // eof is set to 1 if any read hits end-of-file
-
-  int eof = 0;
-
-  if (!first) if (fgets(line,MAXLINE,stdin) == NULL) eof = 1;
-  while (eof == 0 && strspn(line," \t\n\r") == strlen(line))
-    if (fgets(line,MAXLINE,stdin) == NULL) eof = 1;
-  if (fgets(buffer,MAXLINE,stdin) == NULL) eof = 1;
-
-  // if eof, return NULL
-
-  if (eof) return NULL;
-
-  // bracket non-whitespace portion of line
-
-  int start = strspn(line," \t\n\r");
-  int stop = strlen(line) - 1;
-  while (line[stop] == ' ' || line[stop] == '\t' 
-	 || line[stop] == '\n' || line[stop] == '\r') stop--;
-  line[stop+1] = '\0';
-
-  // print keyword into output file
-
-  if (first) printf("%s\n\n",&line[start]);
-  else printf("\n%s\n\n",&line[start]);
-
-  // return ptr to keyword
-
-  return &line[start];
-}
-
-/* count and return words in a single line
-   make copy of line before using strtok so as not to change line
-   trim anything from '#' onward */
-
-int count_words(char *line)
-{
-  int n = strlen(line) + 1;
-  char *copy = (char *) malloc(n*sizeof(char));
-  strcpy(copy,line);
-
-  char *ptr;
-  if (ptr = strchr(copy,'#')) *ptr = '\0';
-
-  if (strtok(copy," \t\n") == NULL) return 0;
-  n = 1;
-  while (strtok(NULL," \t\n")) n++;
-
-  free(copy);
-  return n;
-}
-
-/* copy n lines from stdin to stdout */
-
-void copyline(int n)
-{
-  char line[1000];
-
-  while (n) {
-    fgets(line,MAXLINE,stdin);
-    printf("%s",line);
-    n--;
-  }
-}
-
-/* point (x,y,z) is somewhere in (or near) the replicated box "in"
-   rescale (x,y,z) so it is in the same relative location in box "out"
-   do this by 1st converting to 0-1 coord system in "in" space */
-   
-void scale(struct box in, int nx, int ny, int nz,
-	   struct box out, double *x, double *y, double *z)
-{
-  double rel;
-
-  rel = (*x-in.xlo) / (nx*in.xsize);
-  *x = out.xlo + rel*out.xsize;
-  rel = (*y-in.ylo) / (ny*in.ysize);
-  *y = out.ylo + rel*out.ysize;
-  rel = (*z-in.zlo) / (nz*in.zsize);
-  *z = out.zlo + rel*out.zsize;
-}
-
-/* remap the point (x,y,z) so it is guaranteed to be inside the box 
-   using standard periodic imaging, update image flags */
-
-void pbc(struct box box, double *x, double *y, double *z, 
-	 int *imagex, int *imagey, int *imagez)
-{
-  while (*x < box.xlo) {
-    *x += box.xsize;
-    *imagex -= 1;
-  }
-  while (*x >= box.xhi) {
-    *x -= box.xsize;
-    *imagex += 1;
-  }
-  while (*y < box.ylo) {
-    *y += box.ysize;
-    *imagey -= 1;
-  }
-  while (*y >= box.yhi) {
-    *y -= box.ysize;
-    *imagey += 1;
-  }
-  while (*z < box.zlo) {
-    *z += box.zsize;
-    *imagez -= 1;
-  }
-  while (*z >= box.zhi) {
-    *z -= box.zsize;
-    *imagez += 1;
-  }
-}
-