diff --git a/tools/eff/bohr2ang.py b/tools/eff/bohr2ang.py
new file mode 100644
index 0000000000..e258579aaf
--- /dev/null
+++ b/tools/eff/bohr2ang.py
@@ -0,0 +1,74 @@
+Info="""
+Module name: bohr2ang.py 
+
+Author: (c) Andres Jaramillo-Botero
+California Institute of Technology
+ajaramil@caltech.edu
+Project: pEFF
+Version: August 2009
+
+Usage: python bohr2ang.py
+>>Name of data file (bohr): [datafile]
+
+Results:
+creates a datafile with extension .ang in real units
+
+"""
+import os
+
+currdir=os.getcwd()
+datafile=raw_input("Name of data file (bohr): ")
+bohr2ang=0.529177249
+bperatu2angperfs=0.512396271120794
+f=open(currdir+'/'+datafile,'r')
+w=open(currdir+'/'+datafile+'.ang','w')
+lines=f.readlines()
+atom_flag=False
+vel_flag=False
+for line in lines:
+  if line.find("xlo") > 0:
+    parse=line.split()
+    w.write("%f %f xlo xhi\n"%(float(parse[0])*bohr2ang,float(parse[1])*bohr2ang))
+  elif line.find("ylo") > 0:
+    parse=line.split()
+    w.write("%f %f ylo yhi\n"%(float(parse[0])*bohr2ang,float(parse[1])*bohr2ang))
+  elif line.find("zlo") > 0:
+    parse=line.split()
+    w.write("%f %f zlo zhi\n"%(float(parse[0])*bohr2ang,float(parse[1])*bohr2ang))
+  elif line.find("xy") >= 0:
+    parse=line.split()
+    w.write("%f %f %f xy xz yz\n"%(float(parse[0])*bohr2ang,float(parse[1])*bohr2ang,float(parse[2])*bohr2ang))
+  elif atom_flag and line.strip():
+    parse=line.split()
+    id=parse[0]
+    type=parse[1]
+    q=parse[2]
+    spin=parse[3]
+    eradius=float(parse[4])*bohr2ang
+    x=float(parse[5])*bohr2ang
+    y=float(parse[6])*bohr2ang
+    z=float(parse[7])*bohr2ang
+    rest=" ".join(parse[8:])
+    w.write("%s %s %s %s %f %f %f %f %s\n"%(id,type,q,spin,eradius,x,y,z,rest))
+  elif line.find("Atoms") >= 0:
+    w.write(line)
+    atom_flag=True
+    continue
+  elif vel_flag and line != "\n":
+    parse=line.split()
+    id=parse[0]
+    vx=float(parse[1])*bperatu2angperfs
+    vy=float(parse[2])*bperatu2angperfs
+    vz=float(parse[3])*bperatu2angperfs
+    erv=float(parse[4])*bperatu2angperfs
+    w.write("%s %f %f %f\n"%(id,vx,vy,vz,erv))
+  elif line.find("Velocities") >= 0:
+    w.write(line)
+    atom_flag=False
+    vel_flag=True
+    continue
+  else:
+    w.write(line)
+
+f.close()
+w.close()
diff --git a/tools/eff/restart2data-pEFF.cpp b/tools/eff/restart2data-pEFF.cpp
new file mode 100644
index 0000000000..b5106add2e
--- /dev/null
+++ b/tools/eff/restart2data-pEFF.cpp
@@ -0,0 +1,3520 @@
+/* -----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   www.cs.sandia.gov/~sjplimp/lammps.html
+   Steve Plimpton, sjplimp@sandia.gov, Sandia National Laboratories
+
+   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.
+------------------------------------------------------------------------ */
+
+// Convert a LAMMPS binary restart file into an ASCII text data file
+//
+// Syntax: restart2data restart-file data-file (input-file)
+//         input-file is optional
+//         if specified it will contain LAMMPS input script commands
+//           for mass and force field info
+//         only a few force field styles support this option
+//
+// this serial code must be compiled on a platform that can read the binary
+//   restart file since binary formats are not compatible across all platforms
+// restart-file can have a '%' character to indicate a multiproc restart
+//   file as written by LAMMPS
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+#define MAX_GROUP 32
+
+// same as write_restart.cpp
+
+enum{VERSION,UNITS,NTIMESTEP,DIMENSION,NPROCS,PROCGRID_0,PROCGRID_1,PROCGRID_2,
+       NEWTON_PAIR,NEWTON_BOND,XPERIODIC,YPERIODIC,ZPERIODIC,
+       BOUNDARY_00,BOUNDARY_01,BOUNDARY_10,BOUNDARY_11,BOUNDARY_20,BOUNDARY_21,
+       ATOM_STYLE,NATOMS,NTYPES,
+       NBONDS,NBONDTYPES,BOND_PER_ATOM,
+       NANGLES,NANGLETYPES,ANGLE_PER_ATOM,
+       NDIHEDRALS,NDIHEDRALTYPES,DIHEDRAL_PER_ATOM,
+       NIMPROPERS,NIMPROPERTYPES,IMPROPER_PER_ATOM,
+       BOXLO_0,BOXHI_0,BOXLO_1,BOXHI_1,BOXLO_2,BOXHI_2,
+       SPECIAL_LJ_1,SPECIAL_LJ_2,SPECIAL_LJ_3,
+       SPECIAL_COUL_1,SPECIAL_COUL_2,SPECIAL_COUL_3,
+       XY,XZ,YZ};
+enum{MASS,SHAPE,DIPOLE};
+enum{PAIR,BOND,ANGLE,DIHEDRAL,IMPROPER};
+
+static const char * const cg_type_list[] = 
+  {"none", "lj9_6", "lj12_4", "lj12_6"};
+
+// ---------------------------------------------------------------------
+// Data class to hold problem
+// ---------------------------------------------------------------------
+
+class Data {
+ public:
+
+  // global settings
+
+  char *version;
+  int ntimestep;
+  int nprocs;
+  char *unit_style;
+  int dimension;
+  int px,py,pz;
+  int newton_pair,newton_bond;
+  int xperiodic,yperiodic,zperiodic;
+  int boundary[3][2];
+
+  char *atom_style;
+  int style_angle,style_atomic,style_bond,style_charge,style_dipole;
+  int style_dpd,style_ellipsoid,style_full,style_granular;
+  int style_hybrid,style_molecular,style_peri,style_electron;
+
+  int natoms,nbonds,nangles,ndihedrals,nimpropers;
+  int ntypes,nbondtypes,nangletypes,ndihedraltypes,nimpropertypes;
+  int bond_per_atom,angle_per_atom,dihedral_per_atom,improper_per_atom;
+  int triclinic;
+
+  double xlo,xhi,ylo,yhi,zlo,zhi,xy,xz,yz;
+  double special_lj[4],special_coul[4];
+
+  double cut_lj_global,cut_coul_global,kappa;
+  int offset_flag,mix_flag;
+
+  // force fields
+
+  char *pair_style,*bond_style,*angle_style,*dihedral_style,*improper_style;
+
+  double *pair_born_A,*pair_born_rho,*pair_born_sigma;
+  double *pair_born_C,*pair_born_D;
+  double *pair_buck_A,*pair_buck_rho,*pair_buck_C;
+  double *pair_colloid_A12,*pair_colloid_sigma;
+  double *pair_colloid_d1,*pair_colloid_d2;
+  double *pair_dipole_epsilon,*pair_dipole_sigma;
+  double *pair_dpd_a0,*pair_dpd_gamma;
+  double *pair_charmm_epsilon,*pair_charmm_sigma;
+  double *pair_charmm_eps14,*pair_charmm_sigma14;
+  double *pair_class2_epsilon,*pair_class2_sigma;
+  double *pair_gb_epsilon,*pair_gb_sigma;
+  double *pair_gb_epsa,*pair_gb_epsb,*pair_gb_epsc;
+  double *pair_lj_epsilon,*pair_lj_sigma;
+  double **pair_cg_epsilon,**pair_cg_sigma;
+  int **pair_cg_cmm_type, **pair_setflag;
+  double **pair_cut_coul, **pair_cut_lj, **pair_cut_eff;
+  double *pair_ljexpand_epsilon,*pair_ljexpand_sigma,*pair_ljexpand_shift;
+  double *pair_ljgromacs_epsilon,*pair_ljgromacs_sigma;
+  double *pair_ljsmooth_epsilon,*pair_ljsmooth_sigma;
+  double *pair_morse_d0,*pair_morse_alpha,*pair_morse_r0;
+  double *pair_soft_start,*pair_soft_stop;
+  double *pair_yukawa_A;
+
+  double *bond_class2_r0,*bond_class2_k2,*bond_class2_k3,*bond_class2_k4;
+  double *bond_fene_k,*bond_fene_r0,*bond_fene_epsilon,*bond_fene_sigma;
+  double *bond_feneexpand_k,*bond_feneexpand_r0;
+  double *bond_feneexpand_epsilon,*bond_feneexpand_sigma;
+  double *bond_feneexpand_shift;
+  double *bond_harmonic_k,*bond_harmonic_r0;
+  double *bond_morse_d0,*bond_morse_alpha,*bond_morse_r0;
+  double *bond_nonlinear_epsilon,*bond_nonlinear_r0,*bond_nonlinear_lamda;
+  double *bond_quartic_k,*bond_quartic_b1,*bond_quartic_b2;
+  double *bond_quartic_rc,*bond_quartic_u0;
+
+  double *angle_charmm_k,*angle_charmm_theta0;
+  double *angle_charmm_k_ub,*angle_charmm_r_ub;
+  double *angle_class2_theta0;
+  double *angle_class2_k2,*angle_class2_k3,*angle_class2_k4;
+  double *angle_class2_bb_k,*angle_class2_bb_r1,*angle_class2_bb_r2;
+  double *angle_class2_ba_k1,*angle_class2_ba_k2;
+  double *angle_class2_ba_r1,*angle_class2_ba_r2;
+  double *angle_cosine_k;
+  double *angle_cosine_squared_k,*angle_cosine_squared_theta0;
+  double *angle_harmonic_k,*angle_harmonic_theta0;
+  double *angle_cg_cmm_epsilon,*angle_cg_cmm_sigma;
+  int *angle_cg_cmm_type;
+
+  double *dihedral_charmm_k,*dihedral_charmm_weight;
+  int *dihedral_charmm_multiplicity,*dihedral_charmm_sign;
+  double *dihedral_class2_k1,*dihedral_class2_k2,*dihedral_class2_k3;
+  double *dihedral_class2_phi1,*dihedral_class2_phi2,*dihedral_class2_phi3;
+  double *dihedral_class2_mbt_f1,*dihedral_class2_mbt_f2;
+  double *dihedral_class2_mbt_f3,*dihedral_class2_mbt_r0;
+  double *dihedral_class2_ebt_f1_1,*dihedral_class2_ebt_f2_1;
+  double *dihedral_class2_ebt_f3_1,*dihedral_class2_ebt_r0_1;
+  double *dihedral_class2_ebt_f1_2,*dihedral_class2_ebt_f2_2;
+  double *dihedral_class2_ebt_f3_2,*dihedral_class2_ebt_r0_2;
+  double *dihedral_class2_at_f1_1,*dihedral_class2_at_f2_1;
+  double *dihedral_class2_at_f3_1,*dihedral_class2_at_theta0_1;
+  double *dihedral_class2_at_f1_2,*dihedral_class2_at_f2_2;
+  double *dihedral_class2_at_f3_2,*dihedral_class2_at_theta0_2;
+  double *dihedral_class2_aat_k;
+  double *dihedral_class2_aat_theta0_1,*dihedral_class2_aat_theta0_2;
+  double *dihedral_class2_bb13_k;
+  double *dihedral_class2_bb13_r10,*dihedral_class2_bb13_r30;
+  double *dihedral_harmonic_k;
+  int *dihedral_harmonic_multiplicity,*dihedral_harmonic_sign;
+  double *dihedral_helix_aphi,*dihedral_helix_bphi,*dihedral_helix_cphi;
+  double *dihedral_multi_a1,*dihedral_multi_a2,*dihedral_multi_a3;
+  double *dihedral_multi_a4,*dihedral_multi_a5;
+  double *dihedral_opls_k1,*dihedral_opls_k2;
+  double *dihedral_opls_k3,*dihedral_opls_k4;
+
+  double *improper_class2_k0,*improper_class2_chi0;
+  double *improper_class2_aa_k1,*improper_class2_aa_k2,*improper_class2_aa_k3;
+  double *improper_class2_aa_theta0_1,*improper_class2_aa_theta0_2;
+  double *improper_class2_aa_theta0_3;
+  double *improper_cvff_k;
+  int *improper_cvff_sign,*improper_cvff_multiplicity;
+  double *improper_harmonic_k,*improper_harmonic_chi;
+
+  // atom quantities
+
+  int iatoms,ibonds,iangles,idihedrals,iimpropers;
+
+  double *mass,*shape,*dipole;
+  double *x,*y,*z,*vx,*vy,*vz;
+  double *omegax,*omegay,*omegaz;
+  int *tag,*type,*mask,*image;
+  int *molecule;
+  int *spin;
+  double *ervel, *erforce;
+  double *q,*mux,*muy,*muz,*radius,*density,*vfrac,*rmass;
+  double *s0,*x0x,*x0y,*x0z;
+  double *quatw,*quati,*quatj,*quatk,*angmomx,*angmomy,*angmomz;
+  int *bond_type,*angle_type,*dihedral_type,*improper_type;
+  int *bond_atom1,*bond_atom2;
+  int *angle_atom1,*angle_atom2,*angle_atom3;
+  int *dihedral_atom1,*dihedral_atom2,*dihedral_atom3,*dihedral_atom4;
+  int *improper_atom1,*improper_atom2,*improper_atom3,*improper_atom4;
+  
+  // functions
+
+  Data();
+  void stats();
+  void write(FILE *fp, FILE *fp2=NULL);
+
+  void write_atom_angle(FILE *, int, int, int, int);
+  void write_atom_atomic(FILE *, int, int, int, int);
+  void write_atom_bond(FILE *, int, int, int, int);
+  void write_atom_charge(FILE *, int, int, int, int);
+  void write_atom_dipole(FILE *, int, int, int, int);
+  void write_atom_dpd(FILE *, int, int, int, int);
+  void write_atom_ellipsoid(FILE *, int, int, int, int);
+  void write_atom_full(FILE *, int, int, int, int);
+  void write_atom_granular(FILE *, int, int, int, int);
+  void write_atom_molecular(FILE *, int, int, int, int);
+  void write_atom_peri(FILE *, int, int, int, int);
+  void write_atom_electron(FILE *, int, int, int, int);
+
+  void write_atom_angle_extra(FILE *, int);
+  void write_atom_atomic_extra(FILE *, int);
+  void write_atom_bond_extra(FILE *, int);
+  void write_atom_charge_extra(FILE *, int);
+  void write_atom_dipole_extra(FILE *, int);
+  void write_atom_dpd_extra(FILE *, int);
+  void write_atom_ellipsoid_extra(FILE *, int);
+  void write_atom_full_extra(FILE *, int);
+  void write_atom_granular_extra(FILE *, int);
+  void write_atom_molecular_extra(FILE *, int);
+  void write_atom_peri_extra(FILE *, int);
+  void write_atom_electron_extra(FILE *, int);
+
+  void write_vel_angle(FILE *, int);
+  void write_vel_atomic(FILE *, int);
+  void write_vel_bond(FILE *, int);
+  void write_vel_charge(FILE *, int);
+  void write_vel_dipole(FILE *, int);
+  void write_vel_dpd(FILE *, int);
+  void write_vel_ellipsoid(FILE *, int);
+  void write_vel_full(FILE *, int);
+  void write_vel_granular(FILE *, int);
+  void write_vel_molecular(FILE *, int);
+  void write_vel_peri(FILE *, int);
+  void write_vel_electron(FILE *, int);
+
+  void write_vel_angle_extra(FILE *, int);
+  void write_vel_atomic_extra(FILE *, int);
+  void write_vel_bond_extra(FILE *, int);
+  void write_vel_charge_extra(FILE *, int);
+  void write_vel_dipole_extra(FILE *, int);
+  void write_vel_dpd_extra(FILE *, int);
+  void write_vel_ellipsoid_extra(FILE *, int);
+  void write_vel_full_extra(FILE *, int);
+  void write_vel_granular_extra(FILE *, int);
+  void write_vel_molecular_extra(FILE *, int);
+  void write_vel_peri_extra(FILE *, int);
+  void write_vel_electron_extra(FILE *, int);
+};
+
+// ---------------------------------------------------------------------
+// function prototypes
+// ---------------------------------------------------------------------
+
+void header(FILE *, Data &);
+void set_style(char *, Data &, int);
+void groups(FILE *);
+void type_arrays(FILE *, Data &);
+void force_fields(FILE *, Data &);
+void modify(FILE *);
+void pair(FILE *fp, Data &data, char *style, int flag);
+void bond(FILE *fp, Data &data);
+void angle(FILE *fp, Data &data);
+void dihedral(FILE *fp, Data &data);
+void improper(FILE *fp, Data &data);
+int atom(double *, Data &data);
+
+void allocate_angle(Data &data);
+void allocate_atomic(Data &data);
+void allocate_bond(Data &data);
+void allocate_charge(Data &data);
+void allocate_dipole(Data &data);
+void allocate_dpd(Data &data);
+void allocate_ellipsoid(Data &data);
+void allocate_full(Data &data);
+void allocate_granular(Data &data);
+void allocate_molecular(Data &data);
+void allocate_peri(Data &data);
+void allocate_electron(Data &data);
+
+int atom_angle(double *, Data &, int);
+int atom_atomic(double *, Data &, int);
+int atom_bond(double *, Data &, int);
+int atom_charge(double *, Data &, int);
+int atom_dipole(double *, Data &, int);
+int atom_dpd(double *, Data &, int);
+int atom_ellipsoid(double *, Data &, int);
+int atom_full(double *, Data &, int);
+int atom_granular(double *, Data &, int);
+int atom_molecular(double *, Data &, int);
+int atom_peri(double *, Data &, int);
+int atom_electron(double *, Data &, int);
+
+int read_int(FILE *fp);
+double read_double(FILE *fp);
+char *read_char(FILE *fp);
+
+// ---------------------------------------------------------------------
+// main program
+// ---------------------------------------------------------------------
+
+int main (int argc, char **argv)
+{
+  // syntax error check
+
+    if ((argc != 3) && (argc !=4)) {
+    printf("Syntax: restart2data restart-file data-file (input-file)\n");
+    return 1;
+  }
+
+  // if restart file contains '%', file = filename with % replaced by "base"
+  // else file = single file
+
+  int multiproc;
+  char *file,*ptr;
+
+  if (ptr = strchr(argv[1],'%')) {
+    multiproc = 1;
+    file = new char[strlen(argv[1]) + 16];
+    *ptr = '\0';
+    sprintf(file,"%s%s%s",argv[1],"base",ptr+1);
+  } else {
+    multiproc = 0;
+    file = argv[1];
+  }
+
+  // open single restart file or base file for multiproc case
+
+  printf("Reading restart file ...\n");
+  FILE *fp = fopen(file,"rb");
+  if (fp == NULL) {
+    printf("ERROR: Cannot open restart file %s\n",file);
+    return 1;
+  }
+
+  // read beginning of restart file
+
+  Data data;
+
+  header(fp,data);
+  groups(fp);
+  type_arrays(fp,data);
+  force_fields(fp,data);
+  modify(fp);
+
+  // read atoms from single or multiple restart files
+
+  double *buf = NULL;
+  int n,m;
+  int maxbuf = 0;
+  data.iatoms = data.ibonds = data.iangles = 
+    data.idihedrals = data.iimpropers = 0;
+
+  for (int iproc = 0; iproc < data.nprocs; iproc++) {
+    if (multiproc) {
+      fclose(fp);
+      sprintf(file,"%s%d%s",argv[1],iproc,ptr+1);
+      fp = fopen(file,"rb");
+      if (fp == NULL) {
+        printf("ERROR: Cannot open restart file %s\n",file);
+        return 1;
+      }
+    }
+    n = read_int(fp);
+
+    if (n > maxbuf) {
+      maxbuf = n;
+      delete [] buf;
+      buf = new double[maxbuf];
+    }
+
+    fread(buf,sizeof(double),n,fp);
+
+    m = 0;
+    while (m < n) m += atom(&buf[m],data);
+  }
+
+  fclose(fp);
+
+  // print out stats
+
+  data.stats();
+
+  // write out data file and no input file
+
+  if (argc == 3) {
+    printf("Writing data file ...\n");
+    fp = fopen(argv[2],"w");
+    if (fp == NULL) {
+      printf("ERROR: Cannot open data file %s\n",argv[2]);
+      return 1;
+    }
+    data.write(fp);
+    fclose(fp);
+
+  // write out data file and input file
+
+  } else {
+    printf("Writing data file ...\n");
+    fp = fopen(argv[2],"w");
+    if (fp == NULL) {
+      printf("ERROR: Cannot open data file %s\n",argv[2]);
+      return 1;
+    }
+    printf("Writing input file ...\n");
+    FILE *fp2 = fopen(argv[3],"w");
+    if (fp2 == NULL) {
+      printf("ERROR: Cannot open input file %s\n",argv[3]);
+      return 1;
+    }
+
+    data.write(fp,fp2);
+    fclose(fp);
+    fclose(fp2);
+  }
+  
+  return 0;
+}
+
+// ---------------------------------------------------------------------
+// read header of restart file
+// ---------------------------------------------------------------------
+
+void header(FILE *fp, Data &data)
+{
+  char *version = "7 Jul 2009";
+
+  data.triclinic = 0;
+
+  int flag;
+  flag = read_int(fp);
+
+  while (flag >= 0) {
+
+    if (flag == VERSION) {
+      data.version = read_char(fp);
+      if (strcmp(version,data.version) != 0) {
+	char *str = "Restart file version does not match restart2data version";
+	printf("WARNING %s\n",str);
+	printf("  restart2data version = %s\n",version);
+      }
+    }
+    else if (flag == UNITS) data.unit_style = read_char(fp);
+    else if (flag == NTIMESTEP) data.ntimestep = read_int(fp);
+    else if (flag == DIMENSION) data.dimension = read_int(fp);
+    else if (flag == NPROCS) data.nprocs = read_int(fp);
+    else if (flag == PROCGRID_0) data.px = read_int(fp);
+    else if (flag == PROCGRID_1) data.py = read_int(fp);
+    else if (flag == PROCGRID_2) data.pz = read_int(fp);
+    else if (flag == NEWTON_PAIR) data.newton_pair = read_int(fp);
+    else if (flag == NEWTON_BOND) data.newton_bond = read_int(fp);
+    else if (flag == XPERIODIC) data.xperiodic = read_int(fp);
+    else if (flag == YPERIODIC) data.yperiodic = read_int(fp);
+    else if (flag == ZPERIODIC) data.zperiodic = read_int(fp);
+    else if (flag == BOUNDARY_00) data.boundary[0][0] = read_int(fp);
+    else if (flag == BOUNDARY_01) data.boundary[0][1] = read_int(fp);
+    else if (flag == BOUNDARY_10) data.boundary[1][0] = read_int(fp);
+    else if (flag == BOUNDARY_11) data.boundary[1][1] = read_int(fp);
+    else if (flag == BOUNDARY_20) data.boundary[2][0] = read_int(fp);
+    else if (flag == BOUNDARY_21) data.boundary[2][1] = read_int(fp);
+
+    // if atom_style = hybrid:
+    //   set data_style_hybrid to # of sub-styles
+    //   read additional sub-class arguments
+    //   set sub-styles to 1 to N
+
+    else if (flag == ATOM_STYLE) {
+      data.style_angle = data.style_atomic = data.style_bond = 
+	data.style_charge = data.style_dipole =	data.style_dpd =
+	data.style_ellipsoid = data.style_full = data.style_granular =
+	data.style_hybrid = data.style_molecular = data.style_peri = data.style_electron = 0;
+
+      data.atom_style = read_char(fp);
+      set_style(data.atom_style,data,1);
+
+      if (strcmp(data.atom_style,"hybrid") == 0) {
+	int nwords = read_int(fp);
+	set_style(data.atom_style,data,nwords);
+	char *substyle;
+	for (int i = 1; i <= nwords; i++) {
+	  substyle = read_char(fp);
+	  set_style(substyle,data,i);
+	}
+      }
+    }
+
+    else if (flag == NATOMS) data.natoms = static_cast<int> (read_double(fp));
+    else if (flag == NTYPES) data.ntypes = read_int(fp);
+    else if (flag == NBONDS) data.nbonds = read_int(fp);
+    else if (flag == NBONDTYPES) data.nbondtypes = read_int(fp);
+    else if (flag == BOND_PER_ATOM) data.bond_per_atom = read_int(fp);
+    else if (flag == NANGLES) data.nangles = read_int(fp);
+    else if (flag == NANGLETYPES) data.nangletypes = read_int(fp);
+    else if (flag == ANGLE_PER_ATOM) data.angle_per_atom = read_int(fp);
+    else if (flag == NDIHEDRALS) data.ndihedrals = read_int(fp);
+    else if (flag == NDIHEDRALTYPES) data.ndihedraltypes = read_int(fp);
+    else if (flag == DIHEDRAL_PER_ATOM) data.dihedral_per_atom = read_int(fp);
+    else if (flag == NIMPROPERS) data.nimpropers = read_int(fp);
+    else if (flag == NIMPROPERTYPES) data.nimpropertypes = read_int(fp);
+    else if (flag == IMPROPER_PER_ATOM) data.improper_per_atom = read_int(fp);
+    else if (flag == BOXLO_0) data.xlo = read_double(fp);
+    else if (flag == BOXHI_0) data.xhi = read_double(fp);
+    else if (flag == BOXLO_1) data.ylo = read_double(fp);
+    else if (flag == BOXHI_1) data.yhi = read_double(fp);
+    else if (flag == BOXLO_2) data.zlo = read_double(fp);
+    else if (flag == BOXHI_2) data.zhi = read_double(fp);
+    else if (flag == SPECIAL_LJ_1) data.special_lj[1] = read_double(fp);
+    else if (flag == SPECIAL_LJ_2) data.special_lj[2] = read_double(fp);
+    else if (flag == SPECIAL_LJ_3) data.special_lj[3] = read_double(fp);
+    else if (flag == SPECIAL_COUL_1) data.special_coul[1] = read_double(fp);
+    else if (flag == SPECIAL_COUL_2) data.special_coul[2] = read_double(fp);
+    else if (flag == SPECIAL_COUL_3) data.special_coul[3] = read_double(fp);
+    else if (flag == XY) {
+      data.triclinic = 1;
+      data.xy = read_double(fp);
+    } else if (flag == XZ) {
+      data.triclinic = 1;
+      data.xz = read_double(fp);
+    } else if (flag == YZ) {
+      data.triclinic = 1;
+      data.yz = read_double(fp);
+    } else {
+      printf("ERROR: Invalid flag in header section of restart file %d\n",
+	     flag);
+      exit(1);
+    }
+
+    flag = read_int(fp);
+  }
+}
+
+// ---------------------------------------------------------------------
+// set atom style to flag
+// ---------------------------------------------------------------------
+
+void set_style(char *name, Data &data, int flag)
+{
+  if (strcmp(name,"angle") == 0) data.style_angle = flag;
+  else if (strcmp(name,"atomic") == 0) data.style_atomic = flag;
+  else if (strcmp(name,"bond") == 0) data.style_bond = flag;
+  else if (strcmp(name,"charge") == 0) data.style_charge = flag;
+  else if (strcmp(name,"dipole") == 0) data.style_dipole = flag;
+  else if (strcmp(name,"dpd") == 0) data.style_dpd = flag;
+  else if (strcmp(name,"ellipsoid") == 0) data.style_ellipsoid = flag;
+  else if (strcmp(name,"full") == 0) data.style_full = flag;
+  else if (strcmp(name,"granular") == 0) data.style_granular = flag;
+  else if (strcmp(name,"hybrid") == 0) data.style_hybrid = flag;
+  else if (strcmp(name,"molecular") == 0) data.style_molecular = flag;
+  else if (strcmp(name,"peri") == 0) data.style_peri = flag;
+  else if (strcmp(name,"electron") == 0) data.style_electron = flag;
+  else {
+    printf("ERROR: Unknown atom style %s\n",name);
+    exit(1);
+  }
+}
+
+// ---------------------------------------------------------------------
+// read group info from restart file, just ignore it
+// ---------------------------------------------------------------------
+
+void groups(FILE *fp)
+{
+  int ngroup = read_int(fp);
+
+  int n;
+  char *name;
+
+  // use count to not change restart format with deleted groups
+  // remove this on next major release
+
+  int count = 0;
+  for (int i = 0; i < MAX_GROUP; i++) {
+    name = read_char(fp);
+    delete [] name;
+    count++;
+    if (count == ngroup) break;
+  }
+}
+
+// ---------------------------------------------------------------------
+// read type arrays from restart file
+// ---------------------------------------------------------------------
+
+void type_arrays(FILE *fp, Data &data)
+{
+  data.mass = NULL;
+  data.shape = NULL;
+  data.dipole = NULL;
+
+  int flag;
+  flag = read_int(fp);
+
+  while (flag >= 0) {
+
+    if (flag == MASS) {
+      data.mass = new double[data.ntypes+1];
+      fread(&data.mass[1],sizeof(double),data.ntypes,fp);
+    } else if (flag == SHAPE) {
+      data.shape = new double[3*(data.ntypes+1)];
+      fread(&data.shape[3],sizeof(double),3*data.ntypes,fp);
+    } else if (flag == DIPOLE) {
+      data.dipole = new double[data.ntypes+1];
+      fread(&data.dipole[1],sizeof(double),data.ntypes,fp);
+    } else {
+      printf("ERROR: Invalid flag in type arrays section of restart file %d\n",
+	     flag);
+      exit(1);
+    }
+
+    flag = read_int(fp);
+  }
+}
+
+// ---------------------------------------------------------------------
+// read force-field info from restart file
+// ---------------------------------------------------------------------
+
+void force_fields(FILE *fp, Data &data)
+{
+  data.pair_style = data.bond_style = data.angle_style =
+    data.dihedral_style = data.improper_style = NULL;
+
+  int flag;
+  flag = read_int(fp);
+
+  while (flag >= 0) {
+
+    if (flag == PAIR) {
+      data.pair_style = read_char(fp);
+      pair(fp,data,data.pair_style,1);
+    } else if (flag == BOND) {
+      data.bond_style = read_char(fp);
+      bond(fp,data);
+    } else if (flag == ANGLE) {
+      data.angle_style = read_char(fp);
+      angle(fp,data);
+    } else if (flag == DIHEDRAL) {
+      data.dihedral_style = read_char(fp);
+      dihedral(fp,data);
+    } else if (flag == IMPROPER) {
+      data.improper_style = read_char(fp);
+      improper(fp,data);
+    } else {
+      printf("ERROR: Invalid flag in force fields section of restart file %d\n",
+	     flag);
+      exit(1);
+    }
+
+    flag = read_int(fp);
+  }
+}
+
+// ---------------------------------------------------------------------
+// read fix info from restart file, just ignore it
+// ---------------------------------------------------------------------
+
+void modify(FILE *fp)
+{
+  char *buf;
+  int n;
+
+  // nfix = # of fix entries with state
+
+  int nfix = read_int(fp);
+
+  // read each entry with id string, style string, chunk of data
+
+  for (int i = 0; i < nfix; i++) {
+    buf = read_char(fp); delete [] buf;
+    buf = read_char(fp); delete [] buf;
+    buf = read_char(fp); delete [] buf;
+  }
+
+  // nfix = # of fix entries with peratom info
+
+  int nfix_peratom = read_int(fp);
+
+  // read each entry with id string, style string, maxsize of one atom data
+
+  for (int i = 0; i < nfix_peratom; i++) {
+    buf = read_char(fp); delete [] buf;
+    buf = read_char(fp); delete [] buf;
+    n = read_int(fp);
+  }
+}
+
+// ---------------------------------------------------------------------
+// read atom info from restart file and store in data struct
+// ---------------------------------------------------------------------
+
+int atom(double *buf, Data &data)
+{
+  // allocate per-atom arrays
+
+  if (data.iatoms == 0) {
+
+    // common to all atom styles
+
+    data.x = new double[data.natoms];
+    data.y = new double[data.natoms];
+    data.z = new double[data.natoms];
+    data.tag = new int[data.natoms];
+    data.type = new int[data.natoms];
+    data.mask = new int[data.natoms];
+    data.image = new int[data.natoms];
+    data.vx = new double[data.natoms];
+    data.vy = new double[data.natoms];
+    data.vz = new double[data.natoms];
+
+    // style-specific arrays
+    // don't worry about re-allocating if style = hybrid
+
+    if (data.style_angle) allocate_angle(data);
+    if (data.style_atomic) allocate_atomic(data);
+    if (data.style_bond) allocate_bond(data);
+    if (data.style_charge) allocate_charge(data);
+    if (data.style_dipole) allocate_dipole(data);
+    if (data.style_dpd) allocate_dpd(data);
+    if (data.style_ellipsoid) allocate_ellipsoid(data);
+    if (data.style_full) allocate_full(data);
+    if (data.style_granular) allocate_granular(data);
+    if (data.style_molecular) allocate_molecular(data);
+    if (data.style_peri) allocate_peri(data);
+    if (data.style_electron) allocate_electron(data);
+  }
+
+  // read atom quantities from buf
+  // if hybrid, loop over all sub-styles in order listed
+  // if hybrid, loop index k will match style setting to insure correct order
+
+  int nloop = 1;
+  if (data.style_hybrid) nloop = data.style_hybrid;
+
+  int iatoms = data.iatoms;
+  int m = 0;
+  for (int k = 1; k <= nloop; k++) {
+    if (k == data.style_angle) m += atom_angle(&buf[m],data,iatoms);
+    if (k == data.style_atomic) m += atom_atomic(&buf[m],data,iatoms);
+    if (k == data.style_bond) m += atom_bond(&buf[m],data,iatoms);
+    if (k == data.style_charge) m += atom_charge(&buf[m],data,iatoms);
+    if (k == data.style_dipole) m += atom_dipole(&buf[m],data,iatoms);
+    if (k == data.style_dpd) m += atom_dpd(&buf[m],data,iatoms);
+    if (k == data.style_ellipsoid) m += atom_ellipsoid(&buf[m],data,iatoms);
+    if (k == data.style_full) m += atom_full(&buf[m],data,iatoms);
+    if (k == data.style_granular) m += atom_granular(&buf[m],data,iatoms);
+    if (k == data.style_molecular) m += atom_molecular(&buf[m],data,iatoms);
+    if (k == data.style_peri) m += atom_peri(&buf[m],data,iatoms);
+    if (k == data.style_electron) m += atom_electron(&buf[m],data,iatoms);
+  }
+
+  data.iatoms++;
+  m = static_cast<int> (buf[0]);
+  return m;
+}
+
+// ---------------------------------------------------------------------
+// read one atom's info from buffer
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+int atom_angle(double *buf, Data &data, int iatoms)
+{
+  int type,atom1,atom2,atom3;
+
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.molecule[iatoms] = static_cast<int> (buf[m++]);
+
+  int n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] < atom1) {
+      data.bond_type[data.ibonds] = type;
+      data.bond_atom1[data.ibonds] = data.tag[iatoms];
+      data.bond_atom2[data.ibonds] = atom1;
+      data.ibonds++;
+    }
+  }
+
+  n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    atom2 = static_cast<int> (buf[m++]);
+    atom3 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] == atom2) {
+      data.angle_type[data.iangles] = type;
+      data.angle_atom1[data.iangles] = atom1;
+      data.angle_atom2[data.iangles] = atom2;
+      data.angle_atom3[data.iangles] = atom3;
+	  data.iangles++;
+    }
+  }
+
+  return m;
+}
+
+int atom_atomic(double *buf, Data &data, int iatoms)
+{
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  return m;
+}
+
+int atom_bond(double *buf, Data &data, int iatoms)
+{
+  int type,atom1;
+
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.molecule[iatoms] = static_cast<int> (buf[m++]);
+
+  int n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] < atom1) {
+      data.bond_type[data.ibonds] = type;
+      data.bond_atom1[data.ibonds] = data.tag[iatoms];
+      data.bond_atom2[data.ibonds] = atom1;
+      data.ibonds++;
+    }
+  }
+
+  return m;
+}
+
+int atom_charge(double *buf, Data &data, int iatoms)
+{
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.q[iatoms] = buf[m++];
+
+  return m;
+}
+
+int atom_dipole(double *buf, Data &data, int iatoms)
+{
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.q[iatoms] = buf[m++];
+  data.mux[iatoms] = buf[m++];
+  data.muy[iatoms] = buf[m++];
+  data.muz[iatoms] = buf[m++];
+  
+  return m;
+}
+
+int atom_dpd(double *buf, Data &data, int iatoms)
+{
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  return m;
+}
+
+int atom_ellipsoid(double *buf, Data &data, int iatoms)
+{
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.quatw[iatoms] = buf[m++];
+  data.quati[iatoms] = buf[m++];
+  data.quatj[iatoms] = buf[m++];
+  data.quatk[iatoms] = buf[m++];
+  data.angmomx[iatoms] = buf[m++];
+  data.angmomy[iatoms] = buf[m++];
+  data.angmomz[iatoms] = buf[m++];
+  
+  return m;
+}
+
+int atom_granular(double *buf, Data &data, int iatoms)
+{
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.radius[iatoms] = buf[m++];
+  data.density[iatoms] = buf[m++];
+  data.omegax[iatoms] = buf[m++];
+  data.omegay[iatoms] = buf[m++];
+  data.omegaz[iatoms] = buf[m++];
+
+  return m;
+}
+
+int atom_full(double *buf, Data &data, int iatoms)
+{
+  int type,atom1,atom2,atom3,atom4;
+
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.q[iatoms] = buf[m++];
+  data.molecule[iatoms] = static_cast<int> (buf[m++]);
+
+  int n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] < atom1) {
+      data.bond_type[data.ibonds] = type;
+      data.bond_atom1[data.ibonds] = data.tag[iatoms];
+      data.bond_atom2[data.ibonds] = atom1;
+      data.ibonds++;
+    }
+  }
+
+  n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    atom2 = static_cast<int> (buf[m++]);
+    atom3 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] == atom2) {
+      data.angle_type[data.iangles] = type;
+      data.angle_atom1[data.iangles] = atom1;
+      data.angle_atom2[data.iangles] = atom2;
+      data.angle_atom3[data.iangles] = atom3;
+      data.iangles++;
+    }
+  }
+
+  n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    atom2 = static_cast<int> (buf[m++]);
+    atom3 = static_cast<int> (buf[m++]);
+    atom4 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] == atom2) {
+      data.dihedral_type[data.idihedrals] = type;
+      data.dihedral_atom1[data.idihedrals] = atom1;
+      data.dihedral_atom2[data.idihedrals] = atom2;
+      data.dihedral_atom3[data.idihedrals] = atom3;
+      data.dihedral_atom4[data.idihedrals] = atom4;
+      data.idihedrals++;
+    }
+  }
+    
+  n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    atom2 = static_cast<int> (buf[m++]);
+    atom3 = static_cast<int> (buf[m++]);
+    atom4 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] == atom2) {
+      data.improper_type[data.iimpropers] = type;
+      data.improper_atom1[data.iimpropers] = atom1;
+      data.improper_atom2[data.iimpropers] = atom2;
+      data.improper_atom3[data.iimpropers] = atom3;
+      data.improper_atom4[data.iimpropers] = atom4;
+      data.iimpropers++;
+    }
+  }
+
+  return m;
+}
+
+int atom_molecular(double *buf, Data &data, int iatoms)
+{
+  int type,atom1,atom2,atom3,atom4;
+
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.molecule[iatoms] = static_cast<int> (buf[m++]);
+
+  int n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] < atom1) {
+      data.bond_type[data.ibonds] = type;
+      data.bond_atom1[data.ibonds] = data.tag[iatoms];
+      data.bond_atom2[data.ibonds] = atom1;
+      data.ibonds++;
+    }
+  }
+
+  n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    atom2 = static_cast<int> (buf[m++]);
+    atom3 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] == atom2) {
+      data.angle_type[data.iangles] = type;
+      data.angle_atom1[data.iangles] = atom1;
+      data.angle_atom2[data.iangles] = atom2;
+      data.angle_atom3[data.iangles] = atom3;
+      data.iangles++;
+    }
+  }
+
+  n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    atom2 = static_cast<int> (buf[m++]);
+    atom3 = static_cast<int> (buf[m++]);
+    atom4 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] == atom2) {
+      data.dihedral_type[data.idihedrals] = type;
+      data.dihedral_atom1[data.idihedrals] = atom1;
+      data.dihedral_atom2[data.idihedrals] = atom2;
+      data.dihedral_atom3[data.idihedrals] = atom3;
+      data.dihedral_atom4[data.idihedrals] = atom4;
+      data.idihedrals++;
+    }
+  }
+
+  n = static_cast<int> (buf[m++]);
+  for (int k = 0; k < n; k++) {
+    type = static_cast<int> (buf[m++]);
+    atom1 = static_cast<int> (buf[m++]);
+    atom2 = static_cast<int> (buf[m++]);
+    atom3 = static_cast<int> (buf[m++]);
+    atom4 = static_cast<int> (buf[m++]);
+    if (data.newton_bond || data.tag[iatoms] == atom2) {
+      data.improper_type[data.iimpropers] = type;
+      data.improper_atom1[data.iimpropers] = atom1;
+      data.improper_atom2[data.iimpropers] = atom2;
+      data.improper_atom3[data.iimpropers] = atom3;
+      data.improper_atom4[data.iimpropers] = atom4;
+      data.iimpropers++;
+    }
+  }
+
+  return m;
+}
+
+int atom_peri(double *buf, Data &data, int iatoms)
+{
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.vfrac[iatoms] = buf[m++];
+  data.density[iatoms] = buf[m++];
+  data.rmass[iatoms] = buf[m++];
+  data.s0[iatoms] = buf[m++];
+  data.x0x[iatoms] = buf[m++];
+  data.x0y[iatoms] = buf[m++];
+  data.x0z[iatoms] = buf[m++];
+
+  return m;
+}
+
+int atom_electron(double *buf, Data &data, int iatoms)
+{
+  int m = 1;
+  data.x[iatoms] = buf[m++];
+  data.y[iatoms] = buf[m++];
+  data.z[iatoms] = buf[m++];
+
+  data.tag[iatoms] = static_cast<int> (buf[m++]);
+  data.type[iatoms] = static_cast<int> (buf[m++]);
+  data.mask[iatoms] = static_cast<int> (buf[m++]);
+  data.image[iatoms] = static_cast<int> (buf[m++]);
+
+  data.vx[iatoms] = buf[m++];
+  data.vy[iatoms] = buf[m++];
+  data.vz[iatoms] = buf[m++];
+
+  data.q[iatoms] = buf[m++];
+  data.spin[iatoms] = static_cast<int> (buf[m++]);
+  data.radius[iatoms] = buf[m++];
+  data.ervel[iatoms] = buf[m++];
+
+  return m;
+}
+
+
+// ---------------------------------------------------------------------
+// per-atom memory allocation routines
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void allocate_angle(Data &data)
+{
+  data.molecule = new int[data.natoms];
+  data.bond_type = new int[data.nbonds];
+  data.bond_atom1 = new int[data.nbonds];
+  data.bond_atom2 = new int[data.nbonds];
+  data.angle_type = new int[data.nangles];
+  data.angle_atom1 = new int[data.nangles];
+  data.angle_atom2 = new int[data.nangles];
+  data.angle_atom3 = new int[data.nangles];
+}
+
+void allocate_atomic(Data &data) {}
+
+void allocate_bond(Data &data)
+{
+  data.molecule = new int[data.natoms];
+  data.bond_type = new int[data.nbonds];
+  data.bond_atom1 = new int[data.nbonds];
+  data.bond_atom2 = new int[data.nbonds];
+}
+
+void allocate_charge(Data &data)
+{
+  data.q = new double[data.natoms];
+}
+
+void allocate_dipole(Data &data)
+{
+  data.q = new double[data.natoms];
+  data.mux = new double[data.natoms];
+  data.muy = new double[data.natoms];
+  data.muz = new double[data.natoms];
+}
+
+void allocate_dpd(Data &data) {}
+
+void allocate_full(Data &data)
+{
+  data.q = new double[data.natoms];
+  data.molecule = new int[data.natoms];
+  data.bond_type = new int[data.nbonds];
+  data.bond_atom1 = new int[data.nbonds];
+  data.bond_atom2 = new int[data.nbonds];
+  data.angle_type = new int[data.nangles];
+  data.angle_atom1 = new int[data.nangles];
+  data.angle_atom2 = new int[data.nangles];
+  data.angle_atom3 = new int[data.nangles];
+  data.dihedral_type = new int[data.ndihedrals];
+  data.dihedral_atom1 = new int[data.ndihedrals];
+  data.dihedral_atom2 = new int[data.ndihedrals];
+  data.dihedral_atom3 = new int[data.ndihedrals];
+  data.dihedral_atom4 = new int[data.ndihedrals];
+  data.improper_type = new int[data.nimpropers];
+  data.improper_atom1 = new int[data.nimpropers];
+  data.improper_atom2 = new int[data.nimpropers];
+  data.improper_atom3 = new int[data.nimpropers];
+  data.improper_atom4 = new int[data.nimpropers];
+}
+
+void allocate_ellipsoid(Data &data)
+{
+  data.quatw = new double[data.natoms];
+  data.quati = new double[data.natoms];
+  data.quatj = new double[data.natoms];
+  data.quatk = new double[data.natoms];
+  data.angmomx = new double[data.natoms];
+  data.angmomy = new double[data.natoms];
+  data.angmomz = new double[data.natoms];
+}
+
+void allocate_granular(Data &data)
+{
+  data.radius = new double[data.natoms];
+  data.density = new double[data.natoms];
+  data.omegax = new double[data.natoms];
+  data.omegay = new double[data.natoms];
+  data.omegaz = new double[data.natoms];
+}
+
+void allocate_molecular(Data &data)
+{
+  data.molecule = new int[data.natoms];
+  data.bond_type = new int[data.nbonds];
+  data.bond_atom1 = new int[data.nbonds];
+  data.bond_atom2 = new int[data.nbonds];
+  data.angle_type = new int[data.nangles];
+  data.angle_atom1 = new int[data.nangles];
+  data.angle_atom2 = new int[data.nangles];
+  data.angle_atom3 = new int[data.nangles];
+  data.dihedral_type = new int[data.ndihedrals];
+  data.dihedral_atom1 = new int[data.ndihedrals];
+  data.dihedral_atom2 = new int[data.ndihedrals];
+  data.dihedral_atom3 = new int[data.ndihedrals];
+  data.dihedral_atom4 = new int[data.ndihedrals];
+  data.improper_type = new int[data.nimpropers];
+  data.improper_atom1 = new int[data.nimpropers];
+  data.improper_atom2 = new int[data.nimpropers];
+  data.improper_atom3 = new int[data.nimpropers];
+  data.improper_atom4 = new int[data.nimpropers];
+}
+
+void allocate_peri(Data &data)
+{
+  data.vfrac = new double[data.natoms];
+  data.density = new double[data.natoms];
+  data.rmass = new double[data.natoms];
+  data.s0 = new double[data.natoms];
+  data.x0x = new double[data.natoms];
+  data.x0y = new double[data.natoms];
+  data.x0z = new double[data.natoms];
+}
+
+void allocate_electron(Data &data)
+{
+  data.q = new double[data.natoms];
+  data.spin = new int[data.natoms];
+  data.radius = new double[data.natoms];
+  data.ervel = new double[data.natoms];
+}
+
+
+// ---------------------------------------------------------------------
+// pair coeffs
+// one section for each pair style
+// flag = 1, read all coeff info and allocation arrays
+// flag = 0, just read global settings (when called recursively by hybrid)
+// ---------------------------------------------------------------------
+
+void pair(FILE *fp, Data &data, char *style, int flag)
+{
+  int i,j,m;
+  int itmp;
+  double rtmp;
+
+  if (strcmp(style,"none") == 0) {
+
+  } else if (strcmp(style,"airebo") == 0) {
+
+  } else if (strcmp(style,"born/coul/long") == 0) {
+    double cut_lj_global = read_double(fp);
+    double cut_coul = read_double(fp);
+    int offset_flag = read_int(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_born_A = new double[data.ntypes+1];
+    data.pair_born_rho = new double[data.ntypes+1];
+    data.pair_born_sigma = new double[data.ntypes+1];
+    data.pair_born_C = new double[data.ntypes+1];
+    data.pair_born_D = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_born_A[i] = read_double(fp);
+	    data.pair_born_rho[i] = read_double(fp);
+	    data.pair_born_sigma[i] = read_double(fp);
+	    data.pair_born_C[i] = read_double(fp);
+	    data.pair_born_D[i] = read_double(fp);
+	    double cut_lj = read_double(fp);
+	  } else {
+	    double born_A = read_double(fp);
+	    double born_rho = read_double(fp);
+	    double born_sigma = read_double(fp);
+	    double born_C = read_double(fp);
+	    double born_D = read_double(fp);
+	    double cut_lj = read_double(fp);
+	  }
+	}
+      }
+
+  } else if ((strcmp(style,"buck") == 0)  ||
+	   (strcmp(style,"buck/coul/cut") == 0) ||
+	   (strcmp(style,"buck/coul/long") == 0) ||
+	   (strcmp(style,"buck/coul") == 0)) {
+
+    if (strcmp(style,"buck") == 0) {
+      m = 0;
+      double cut_lj_global = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"buck/coul/cut") == 0) {
+      m = 1;
+      double cut_lj_global = read_double(fp);
+      double cut_lj_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"buck/coul/long") == 0) {
+      m = 0;
+      double cut_lj_global = read_double(fp);
+      double cut_lj_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"buck/coul") == 0) {
+      m = 0;
+      double cut_buck_global = read_double(fp);
+      double cut_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+      int ewald_order = read_int(fp);
+    }
+
+    if (!flag) return;
+
+    data.pair_buck_A = new double[data.ntypes+1];
+    data.pair_buck_rho = new double[data.ntypes+1];
+    data.pair_buck_C = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_buck_A[i] = read_double(fp);
+	    data.pair_buck_rho[i] = read_double(fp);
+	    data.pair_buck_C[i] = read_double(fp);
+	    double cut_lj = read_double(fp);
+	    if (m) double cut_coul = read_double(fp);
+	  } else {
+	    double buck_A = read_double(fp);
+	    double buck_rho = read_double(fp);
+	    double buck_C = read_double(fp);
+	    double cut_lj = read_double(fp);
+	    if (m) double cut_coul = read_double(fp);
+	  }
+	}
+      }
+
+  } else if (strcmp(style,"colloid") == 0) {
+
+    double cut_global = read_double(fp);
+    int offset_flag = read_int(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_colloid_A12 = new double[data.ntypes+1];
+    data.pair_colloid_sigma = new double[data.ntypes+1];
+    data.pair_colloid_d1 = new double[data.ntypes+1];
+    data.pair_colloid_d2 = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_colloid_A12[i] = read_double(fp);
+	    data.pair_colloid_sigma[i] = read_double(fp);
+	    data.pair_colloid_d1[i] = read_double(fp);
+	    data.pair_colloid_d2[i] = read_double(fp);
+	    double cut_lj = read_double(fp);
+	  } else {
+	    double colloid_A12 = read_double(fp);
+	    double colloid_sigma = read_double(fp);
+	    double colloid_d1 = read_double(fp);
+	    double colloid_d2 = read_double(fp);
+	    double cut_lj = read_double(fp);
+	  }
+	}
+      }
+    
+  } else if ((strcmp(style,"coul/cut") == 0) ||
+	     (strcmp(style,"coul/debye") == 0) ||
+	     (strcmp(style,"coul/long") == 0)) {
+
+    if (strcmp(style,"coul/cut") == 0) {
+      double cut_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"coul/debye") == 0) {
+      m = 1;
+      double cut_coul = read_double(fp);
+      double kappa = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"coul/long") == 0) {
+      double cut_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    }
+
+    if (!flag) return;
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    double cut_coul = read_double(fp);
+	  } else {
+	    double cut_coul = read_double(fp);
+	  }
+	}
+      }
+
+  } else if (strcmp(style,"eff/cut") == 0) {
+
+    if (strcmp(style,"eff/cut") == 0) {
+      double cut_eff = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    }
+
+    if (!flag) return;
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+        itmp = read_int(fp);
+        if (i == j && itmp == 0) {
+          printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+          exit(1);
+        }
+        if (itmp) {
+          if (i == j) {
+            double cut_eff = read_double(fp);
+          } else {
+            double cut_eff = read_double(fp);
+          }
+        }
+      }
+
+  }  else if (strcmp(style,"dipole/cut") == 0) {
+
+    double cut_lj_global = read_double(fp);
+    double cut_coul_global = read_double(fp);
+    int offset_flag = read_int(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_dipole_epsilon = new double[data.ntypes+1];
+    data.pair_dipole_sigma = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_dipole_epsilon[i] = read_double(fp);
+	    data.pair_dipole_sigma[i] = read_double(fp);
+	    double cut_lj = read_double(fp);
+	    double cut_coul = read_double(fp);
+	  } else {
+	    double dipole_epsilon = read_double(fp);
+	    double dipole_sigma = read_double(fp);
+	    double cut_lj = read_double(fp);
+	    double cut_coul = read_double(fp);
+	  }
+	}
+      }
+
+  } else if (strcmp(style,"dpd") == 0) {
+
+    double temperature = read_double(fp);
+    double cut_global = read_double(fp);
+    int seed = read_int(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_dpd_a0 = new double[data.ntypes+1];
+    data.pair_dpd_gamma = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_dpd_a0[i] = read_double(fp);
+	    data.pair_dpd_gamma[i] = read_double(fp);
+	    double cut = read_double(fp);
+	  } else {
+	    double dpd_a0 = read_double(fp);
+	    double dpd_gamma = read_double(fp);
+	    double cut = read_double(fp);
+	  }
+	}
+      }
+
+  } else if (strcmp(style,"eam") == 0) {
+  } else if (strcmp(style,"eam/opt") == 0) {
+  } else if (strcmp(style,"eam/alloy") == 0) {
+  } else if (strcmp(style,"eam/alloy/opt") == 0) {
+  } else if (strcmp(style,"eam/fs") == 0) {
+  } else if (strcmp(style,"eam/fs/opt") == 0) {
+
+  } else if (strcmp(style,"gayberne") == 0) {
+
+    double gamma = read_double(fp);
+    double upsilon = read_double(fp);
+    double mu = read_double(fp);
+    double cut_global = read_double(fp);
+    int offset_flag = read_int(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_gb_epsilon = new double[data.ntypes+1];
+    data.pair_gb_sigma = new double[data.ntypes+1];
+    data.pair_gb_epsa = new double[data.ntypes+1];
+    data.pair_gb_epsb = new double[data.ntypes+1];
+    data.pair_gb_epsc = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++) {
+      itmp = read_int(fp);
+      if (itmp) {
+	data.pair_gb_epsa[i] = read_double(fp);
+	data.pair_gb_epsb[i] = read_double(fp);
+	data.pair_gb_epsc[i] = read_double(fp);
+	data.pair_gb_epsa[i] = pow(data.pair_gb_epsa[i],-mu);
+	data.pair_gb_epsb[i] = pow(data.pair_gb_epsb[i],-mu);
+	data.pair_gb_epsc[i] = pow(data.pair_gb_epsc[i],-mu);
+      }
+      
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_gb_epsilon[i] = read_double(fp);
+	    data.pair_gb_sigma[i] = read_double(fp);
+	    double cut = read_double(fp);
+	  } else {
+	    double gb_epsilon = read_double(fp);
+	    double gb_sigma = read_double(fp);
+	    double cut = read_double(fp);
+	  }
+	}
+      }
+    }
+
+  } else if ((strcmp(style,"gran/hooke") == 0) ||
+	   (strcmp(style,"gran/hooke/history") == 0) ||
+	   (strcmp(style,"gran/hertz/history") == 0)) {
+
+    double kn = read_double(fp);
+    double kt = read_double(fp);
+    double gamman = read_double(fp);
+    double gammat = read_double(fp);
+    double xmu = read_double(fp);
+    int dampflag = read_int(fp);
+
+  } else if ((strcmp(style,"lj/charmm/coul/charmm") == 0) ||
+	   (strcmp(style,"lj/charmm/coul/charmm/implicit") == 0) ||
+	   (strcmp(style,"lj/charmm/coul/long") == 0) ||
+	   (strcmp(style,"lj/charmm/coul/long/opt") == 0)) {
+
+    if (strcmp(style,"lj/charmm/coul/charmm") == 0) {
+      double cut_lj_inner = read_double(fp);
+      double cut_lj = read_double(fp);
+      double cut_coul_inner = read_double(fp);
+      double cut_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"lj/charmm/coul/charmm/implicit") == 0) {
+      double cut_lj_inner = read_double(fp);
+      double cut_lj = read_double(fp);
+      double cut_coul_inner = read_double(fp);
+      double cut_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if ((strcmp(style,"lj/charmm/coul/long") == 0) ||
+	       (strcmp(style,"lj/charmm/coul/long/opt") == 0)) {
+      double cut_lj_inner = read_double(fp);
+      double cut_lj = read_double(fp);
+      double cut_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    }
+
+    if (!flag) return;
+
+    data.pair_charmm_epsilon = new double[data.ntypes+1];
+    data.pair_charmm_sigma = new double[data.ntypes+1];
+    data.pair_charmm_eps14 = new double[data.ntypes+1];
+    data.pair_charmm_sigma14 = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_charmm_epsilon[i] = read_double(fp);
+	    data.pair_charmm_sigma[i] = read_double(fp);
+	    data.pair_charmm_eps14[i] = read_double(fp);
+	    data.pair_charmm_sigma14[i] = read_double(fp);
+	  } else {
+	    double charmm_epsilon = read_double(fp);
+	    double charmm_sigma = read_double(fp);
+	    double charmm_eps14 = read_double(fp);
+	    double charmm_sigma14 = read_double(fp);
+	  }
+	}
+      }
+
+  } else if ((strcmp(style,"lj/class2") == 0) ||
+	   (strcmp(style,"lj/class2/coul/cut") == 0) ||
+	   (strcmp(style,"lj/class2/coul/long") == 0)) {
+
+    if (strcmp(style,"lj/class2") == 0) {
+      m = 0;
+      double cut_lj_global = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"lj/class2/coul/cut") == 0) {
+      m = 1;
+      double cut_lj_global = read_double(fp);
+      double cut_lj_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"lj/class2/coul/long") == 0) {
+      m = 0;
+      double cut_lj_global = read_double(fp);
+      double cut_lj_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    }
+
+    if (!flag) return;
+
+    data.pair_class2_epsilon = new double[data.ntypes+1];
+    data.pair_class2_sigma = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_class2_epsilon[i] = read_double(fp);
+	    data.pair_class2_sigma[i] = read_double(fp);
+	    double cut_lj = read_double(fp);
+	    if (m) double cut_coul = read_double(fp);
+	  } else {
+	    double class2_epsilon = read_double(fp);
+	    double class2_sigma = read_double(fp);
+	    double cut_lj = read_double(fp);
+	    if (m) double cut_coul = read_double(fp);
+	  }
+	}
+      }
+
+  } else if ((strcmp(style,"lj/cut") == 0) ||
+	   (strcmp(style,"lj/cut/opt") == 0) ||
+	   (strcmp(style,"lj/cut/coul/cut") == 0) ||
+	   (strcmp(style,"lj/cut/coul/debye") == 0) ||
+	   (strcmp(style,"lj/cut/coul/long") == 0) ||
+	   (strcmp(style,"lj/cut/coul/long/tip4p") == 0) ||
+	   (strcmp(style,"lj/coul") == 0)) {
+
+    if ((strcmp(style,"lj/cut") == 0) || (strcmp(style,"lj/cut/opt") == 0)) {
+      m = 0;
+      double cut_lj_global = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"lj/cut/coul/cut") == 0) {
+      m = 1;
+      double cut_lj_global = read_double(fp);
+      double cut_lj_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"lj/cut/coul/debye") == 0) {
+      m = 1;
+      double cut_lj_global = read_double(fp);
+      double cut_lj_coul = read_double(fp);
+      double kappa = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"lj/cut/coul/long") == 0) {
+      m = 0;
+      double cut_lj_global = read_double(fp);
+      double cut_lj_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"lj/cut/coul/long/tip4p") == 0) {
+      m = 0;
+      int typeO = read_int(fp);
+      int typeH = read_int(fp);
+      int typeB = read_int(fp);
+      int typeA = read_int(fp);
+      double qdist = read_double(fp);
+      double cut_lj_global = read_double(fp);
+      double cut_lj_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"lj/coul") == 0) {
+      m = 0;
+      double cut_lj_global = read_double(fp);
+      double cut_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+      int ewald_order = read_int(fp);
+    }
+
+    if (!flag) return;
+
+    data.pair_lj_epsilon = new double[data.ntypes+1];
+    data.pair_lj_sigma = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_lj_epsilon[i] = read_double(fp);
+	    data.pair_lj_sigma[i] = read_double(fp);
+	    double cut_lj = read_double(fp);
+	    if (m) double cut_coul = read_double(fp);
+	  } else {
+	    double lj_epsilon = read_double(fp);
+	    double lj_sigma = read_double(fp);
+	    double cut_lj = read_double(fp);
+	    if (m) double cut_coul = read_double(fp);
+	  }
+	}
+      }
+
+  } else if (strcmp(style,"lj/expand") == 0) {
+
+    double cut_global = read_double(fp);
+    int offset_flag = read_int(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_ljexpand_epsilon = new double[data.ntypes+1];
+    data.pair_ljexpand_sigma = new double[data.ntypes+1];
+    data.pair_ljexpand_shift = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_ljexpand_epsilon[i] = read_double(fp);
+	    data.pair_ljexpand_sigma[i] = read_double(fp);
+	    data.pair_ljexpand_shift[i] = read_double(fp);
+	    double cut_lj = read_double(fp);
+	  } else {
+	    double ljexpand_epsilon = read_double(fp);
+	    double ljexpand_sigma = read_double(fp);
+	    double ljexpand_shift = read_double(fp);
+	    double cut_lj = read_double(fp);
+	  }
+	}
+      }
+
+  } else if ((strcmp(style,"lj/gromacs") == 0) ||
+	     (strcmp(style,"lj/gromacs/coul/gromacs") == 0)) {
+
+    if (strcmp(style,"lj/gromacs") == 0) {
+      m = 1;
+      double cut_inner_global = read_double(fp);
+      double cut_global = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    } else if (strcmp(style,"lj/gromacs/coul/gromacs") == 0) {
+      m = 0;
+      double cut_lj_inner_global = read_double(fp);
+      double cut_lj = read_double(fp);
+      double cut_coul_inner_global = read_double(fp);
+      double cut_coul = read_double(fp);
+      int offset_flag = read_int(fp);
+      int mix_flag = read_int(fp);
+    }
+
+    if (!flag) return;
+
+    data.pair_ljgromacs_epsilon = new double[data.ntypes+1];
+    data.pair_ljgromacs_sigma = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_ljgromacs_epsilon[i] = read_double(fp);
+	    data.pair_ljgromacs_sigma[i] = read_double(fp);
+	    if (m) {
+	      double cut_inner = read_double(fp);
+	      double cut = read_double(fp);
+	    }
+	    } else {
+	    double ljgromacs_epsilon = read_double(fp);
+	    double ljgromacs_sigma = read_double(fp);
+	    if (m) {
+	      double cut_inner = read_double(fp);
+	      double cut = read_double(fp);
+	    }
+	  }
+	}
+      }
+
+  } else if (strcmp(style,"lj/smooth") == 0) {
+
+    double cut_inner_global = read_double(fp);
+    double cut_global = read_double(fp);
+    int offset_flag = read_int(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_ljsmooth_epsilon = new double[data.ntypes+1];
+    data.pair_ljsmooth_sigma = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_ljsmooth_epsilon[i] = read_double(fp);
+	    data.pair_ljsmooth_sigma[i] = read_double(fp);
+	    double cut_inner = read_double(fp);
+	    double cut = read_double(fp);
+	  } else {
+	    double ljsmooth_epsilon = read_double(fp);
+	    double ljsmooth_sigma = read_double(fp);
+	    double cut_inner = read_double(fp);
+	    double cut = read_double(fp);
+	  }
+	}
+      }
+
+  } else if (strcmp(style,"meam") == 0) {
+
+  } else if ((strcmp(style,"morse") == 0) ||
+	     (strcmp(style,"morse/opt") == 0)) {
+
+    double cut_global = read_double(fp);
+    int offset_flag = read_int(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_morse_d0 = new double[data.ntypes+1];
+    data.pair_morse_alpha = new double[data.ntypes+1];
+    data.pair_morse_r0 = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_morse_d0[i] = read_double(fp);
+	    data.pair_morse_alpha[i] = read_double(fp);
+	    data.pair_morse_r0[i] = read_double(fp);
+	    double cut = read_double(fp);
+	  } else {
+	    double morse_d0 = read_double(fp);
+	    double morse_alpha = read_double(fp);
+	    double morse_r0 = read_double(fp);
+	    double cut = read_double(fp);
+	  }
+	}
+      }
+
+  } else if (strcmp(style,"reax") == 0) {
+
+  } else if (strcmp(style,"soft") == 0) {
+
+    double cut_global = read_double(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_soft_start = new double[data.ntypes+1];
+    data.pair_soft_stop = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_soft_start[i] = read_double(fp);
+	    data.pair_soft_stop[i] = read_double(fp);
+	    double cut = read_double(fp);
+	  } else {
+	    double soft_start = read_double(fp);
+	    double soft_stop = read_double(fp);
+	    double cut = read_double(fp);
+	  }
+	}
+      }
+
+  } else if (strcmp(style,"sw") == 0) {
+
+  } else if (strcmp(style,"table") == 0) {
+
+    int tabstyle = read_int(fp);
+    int n = read_int(fp);
+
+  } else if (strcmp(style,"tersoff") == 0) {
+
+  } else if (strcmp(style,"yukawa") == 0) {
+
+    double kappa = read_double(fp);
+    double cut_global = read_double(fp);
+    int offset_flag = read_int(fp);
+    int mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    data.pair_yukawa_A = new double[data.ntypes+1];
+
+    for (i = 1; i <= data.ntypes; i++)
+      for (j = i; j <= data.ntypes; j++) {
+	itmp = read_int(fp);
+	if (i == j && itmp == 0) {
+	  printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+	  exit(1);
+	}
+	if (itmp) {
+	  if (i == j) {
+	    data.pair_yukawa_A[i] = read_double(fp);
+	    double cut = read_double(fp);
+	  } else {
+	    double yukawa_A = read_double(fp);
+	    double cut = read_double(fp);
+	  }
+	}
+      }
+
+  } else if ((strcmp(style,"cg/cmm") == 0) ||
+             (strcmp(style,"cg/cmm/coul/cut") == 0) ||
+             (strcmp(style,"cg/cmm/coul/long") == 0) ) {
+    m = 0;
+    data.cut_lj_global = read_double(fp);
+    data.cut_coul_global = read_double(fp);
+    data.kappa = read_double(fp);
+    data.offset_flag = read_int(fp);
+    data.mix_flag = read_int(fp);
+
+    if (!flag) return;
+
+    const int numtyp=data.ntypes+1;
+    data.pair_cg_cmm_type = new int*[numtyp];
+    data.pair_setflag = new int*[numtyp];
+    data.pair_cg_epsilon = new double*[numtyp];
+    data.pair_cg_sigma = new double*[numtyp];
+    data.pair_cut_lj = new double*[numtyp];
+    if  ((strcmp(style,"cg/cmm/coul/cut") == 0) ||
+             (strcmp(style,"cg/cmm/coul/long") == 0) ) {
+      data.pair_cut_coul = new double*[numtyp];
+      m=1;
+    } else {
+      data.pair_cut_coul = NULL;
+      m=0;
+    }
+    
+    for (i = 1; i <= data.ntypes; i++) {
+      data.pair_cg_cmm_type[i] = new int[numtyp];
+      data.pair_setflag[i] = new int[numtyp];
+      data.pair_cg_epsilon[i] = new double[numtyp];
+      data.pair_cg_sigma[i] = new double[numtyp];
+      data.pair_cut_lj[i] = new double[numtyp];
+      if ((strcmp(style,"cg/cmm/coul/cut") == 0) ||
+          (strcmp(style,"cg/cmm/coul/long") == 0) ) {
+        data.pair_cut_coul[i] = new double[numtyp];
+      }
+
+      for (j = i; j <= data.ntypes; j++) {
+        itmp = read_int(fp);
+        data.pair_setflag[i][j] = itmp;        
+        if (i == j && itmp == 0) {
+          printf("ERROR: Pair coeff %d,%d is not in restart file\n",i,j);
+          exit(1);
+        }
+        if (itmp) {
+          data.pair_cg_cmm_type[i][j] = read_int(fp);
+          data.pair_cg_epsilon[i][j] = read_double(fp);
+          data.pair_cg_sigma[i][j] = read_double(fp);
+          data.pair_cut_lj[i][j] = read_double(fp);
+          if (m) {
+	    data.pair_cut_lj[i][j] = read_double(fp);
+	    data.pair_cut_coul[i][j] = read_double(fp);
+	  }
+        } 
+      }
+    }
+
+  } else if ((strcmp(style,"hybrid") == 0) ||
+	     (strcmp(style,"hybrid/overlay") == 0)) {
+
+    // for each substyle of hybrid,
+    //   read its settings by calling pair() recursively with flag = 0
+    //   so that coeff array allocation is skipped
+
+    int nstyles = read_int(fp);
+    for (int i = 0; i < nstyles; i++) {
+      char *substyle = read_char(fp);
+      pair(fp,data,substyle,0);
+    }
+
+  } else {
+    printf("ERROR: Unknown pair style %s\n",style);
+    exit(1);
+  }
+}
+
+// ---------------------------------------------------------------------
+// bond coeffs
+// one section for each bond style
+// ---------------------------------------------------------------------
+
+void bond(FILE *fp, Data &data)
+{
+  if (strcmp(data.bond_style,"none") == 0) {
+
+  } else if (strcmp(data.bond_style,"class2") == 0) {
+
+    data.bond_class2_r0 = new double[data.nbondtypes+1];
+    data.bond_class2_k2 = new double[data.nbondtypes+1];
+    data.bond_class2_k3 = new double[data.nbondtypes+1];
+    data.bond_class2_k4 = new double[data.nbondtypes+1];
+    fread(&data.bond_class2_r0[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_class2_k2[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_class2_k3[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_class2_k4[1],sizeof(double),data.nbondtypes,fp);
+
+  } else if (strcmp(data.bond_style,"fene") == 0) {
+
+    data.bond_fene_k = new double[data.nbondtypes+1];
+    data.bond_fene_r0 = new double[data.nbondtypes+1];
+    data.bond_fene_epsilon = new double[data.nbondtypes+1];
+    data.bond_fene_sigma = new double[data.nbondtypes+1];
+    fread(&data.bond_fene_k[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_fene_r0[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_fene_epsilon[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_fene_sigma[1],sizeof(double),data.nbondtypes,fp);
+
+  } else if (strcmp(data.bond_style,"fene/expand") == 0) {
+
+    data.bond_feneexpand_k = new double[data.nbondtypes+1];
+    data.bond_feneexpand_r0 = new double[data.nbondtypes+1];
+    data.bond_feneexpand_epsilon = new double[data.nbondtypes+1];
+    data.bond_feneexpand_sigma = new double[data.nbondtypes+1];
+    data.bond_feneexpand_shift = new double[data.nbondtypes+1];
+    fread(&data.bond_feneexpand_k[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_feneexpand_r0[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_feneexpand_epsilon[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_feneexpand_sigma[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_feneexpand_shift[1],sizeof(double),data.nbondtypes,fp);
+
+  } else if (strcmp(data.bond_style,"harmonic") == 0) {
+
+    data.bond_harmonic_k = new double[data.nbondtypes+1];
+    data.bond_harmonic_r0 = new double[data.nbondtypes+1];
+    fread(&data.bond_harmonic_k[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_harmonic_r0[1],sizeof(double),data.nbondtypes,fp);
+
+  } else if (strcmp(data.bond_style,"morse") == 0) {
+
+    data.bond_morse_d0 = new double[data.nbondtypes+1];
+    data.bond_morse_alpha = new double[data.nbondtypes+1];
+    data.bond_morse_r0 = new double[data.nbondtypes+1];
+    fread(&data.bond_morse_d0[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_morse_alpha[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_morse_r0[1],sizeof(double),data.nbondtypes,fp);
+
+  } else if (strcmp(data.bond_style,"nonlinear") == 0) {
+
+    data.bond_nonlinear_epsilon = new double[data.nbondtypes+1];
+    data.bond_nonlinear_r0 = new double[data.nbondtypes+1];
+    data.bond_nonlinear_lamda = new double[data.nbondtypes+1];
+    fread(&data.bond_nonlinear_epsilon[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_nonlinear_r0[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_nonlinear_lamda[1],sizeof(double),data.nbondtypes,fp);
+
+  } else if (strcmp(data.bond_style,"quartic") == 0) {
+
+    data.bond_quartic_k = new double[data.nbondtypes+1];
+    data.bond_quartic_b1 = new double[data.nbondtypes+1];
+    data.bond_quartic_b2 = new double[data.nbondtypes+1];
+    data.bond_quartic_rc = new double[data.nbondtypes+1];
+    data.bond_quartic_u0 = new double[data.nbondtypes+1];
+    fread(&data.bond_quartic_k[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_quartic_b1[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_quartic_b2[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_quartic_rc[1],sizeof(double),data.nbondtypes,fp);
+    fread(&data.bond_quartic_u0[1],sizeof(double),data.nbondtypes,fp);
+
+  } else if (strcmp(data.bond_style,"hybrid") == 0) {
+
+    int nstyles = read_int(fp);
+    for (int i = 0; i < nstyles; i++)
+      char *substyle = read_char(fp);
+
+  } else {
+    printf("ERROR: Unknown bond style %s\n",data.bond_style);
+    exit(1);
+  }
+}
+
+// ---------------------------------------------------------------------
+// angle coeffs
+// one section for each angle style
+// ---------------------------------------------------------------------
+
+void angle(FILE *fp, Data &data)
+{
+  if (strcmp(data.angle_style,"none") == 0) {
+
+  } else if (strcmp(data.angle_style,"charmm") == 0) {
+
+    data.angle_charmm_k = new double[data.nangletypes+1];
+    data.angle_charmm_theta0 = new double[data.nangletypes+1];
+    data.angle_charmm_k_ub = new double[data.nangletypes+1];
+    data.angle_charmm_r_ub = new double[data.nangletypes+1];
+    fread(&data.angle_charmm_k[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_charmm_theta0[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_charmm_k_ub[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_charmm_r_ub[1],sizeof(double),data.nangletypes,fp);
+
+  } else if (strcmp(data.angle_style,"class2") == 0) {
+
+    data.angle_class2_theta0 = new double[data.nangletypes+1];
+    data.angle_class2_k2 = new double[data.nangletypes+1];
+    data.angle_class2_k3 = new double[data.nangletypes+1];
+    data.angle_class2_k4 = new double[data.nangletypes+1];
+
+    data.angle_class2_bb_k = new double[data.nangletypes+1];
+    data.angle_class2_bb_r1 = new double[data.nangletypes+1];
+    data.angle_class2_bb_r2 = new double[data.nangletypes+1];
+
+    data.angle_class2_ba_k1 = new double[data.nangletypes+1];
+    data.angle_class2_ba_k2 = new double[data.nangletypes+1];
+    data.angle_class2_ba_r1 = new double[data.nangletypes+1];
+    data.angle_class2_ba_r2 = new double[data.nangletypes+1];
+
+    fread(&data.angle_class2_theta0[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_class2_k2[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_class2_k3[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_class2_k4[1],sizeof(double),data.nangletypes,fp);
+
+    fread(&data.angle_class2_bb_k[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_class2_bb_r1[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_class2_bb_r2[1],sizeof(double),data.nangletypes,fp);
+
+    fread(&data.angle_class2_ba_k1[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_class2_ba_k2[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_class2_ba_r1[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_class2_ba_r2[1],sizeof(double),data.nangletypes,fp);
+
+  } else if (strcmp(data.angle_style,"cosine") == 0) {
+
+    data.angle_cosine_k = new double[data.nangletypes+1];
+    fread(&data.angle_cosine_k[1],sizeof(double),data.nangletypes,fp);
+
+  } else if ((strcmp(data.angle_style,"cosine/squared") == 0) ||
+             (strcmp(data.angle_style,"cosine/delta") == 0)) {
+
+    data.angle_cosine_squared_k = new double[data.nangletypes+1];
+    data.angle_cosine_squared_theta0 = new double[data.nangletypes+1];
+    fread(&data.angle_cosine_squared_k[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_cosine_squared_theta0[1],
+	  sizeof(double),data.nangletypes,fp);
+
+  } else if (strcmp(data.angle_style,"harmonic") == 0) {
+
+    data.angle_harmonic_k = new double[data.nangletypes+1];
+    data.angle_harmonic_theta0 = new double[data.nangletypes+1];
+    fread(&data.angle_harmonic_k[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_harmonic_theta0[1],sizeof(double),data.nangletypes,fp);
+
+  } else if (strcmp(data.angle_style,"cg/cmm") == 0) {
+
+    data.angle_harmonic_k = new double[data.nangletypes+1];
+    data.angle_harmonic_theta0 = new double[data.nangletypes+1];
+    data.angle_cg_cmm_epsilon = new double[data.nangletypes+1];
+    data.angle_cg_cmm_sigma = new double[data.nangletypes+1];
+    double *angle_cg_cmm_rcut = new double[data.nangletypes+1];
+    data.angle_cg_cmm_type = new int[data.nangletypes+1];
+    
+    fread(&data.angle_harmonic_k[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_harmonic_theta0[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_cg_cmm_epsilon[1],sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_cg_cmm_sigma[1],sizeof(double),data.nangletypes,fp);
+    fread(angle_cg_cmm_rcut,sizeof(double),data.nangletypes,fp);
+    fread(&data.angle_cg_cmm_type[1],sizeof(int),data.nangletypes,fp);
+
+  } else if (strcmp(data.angle_style,"hybrid") == 0) {
+
+    int nstyles = read_int(fp);
+    for (int i = 0; i < nstyles; i++)
+      char *substyle = read_char(fp);
+
+  } else {
+    printf("ERROR: Unknown angle style %s\n",data.angle_style);
+    exit(1);
+  }
+}
+
+// ---------------------------------------------------------------------
+// dihedral coeffs
+// one section for each dihedral style
+// ---------------------------------------------------------------------
+
+void dihedral(FILE *fp, Data &data)
+{
+  if (strcmp(data.dihedral_style,"none") == 0) {
+
+  } else if (strcmp(data.dihedral_style,"charmm") == 0) {
+
+    data.dihedral_charmm_k = new double[data.ndihedraltypes+1];
+    data.dihedral_charmm_multiplicity = new int[data.ndihedraltypes+1];
+    data.dihedral_charmm_sign = new int[data.ndihedraltypes+1];
+    data.dihedral_charmm_weight = new double[data.ndihedraltypes+1];
+    fread(&data.dihedral_charmm_k[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_charmm_multiplicity[1],sizeof(int),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_charmm_sign[1],sizeof(int),data.ndihedraltypes,fp);
+    fread(&data.dihedral_charmm_weight[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+
+  } else if (strcmp(data.dihedral_style,"class2") == 0) {
+
+    data.dihedral_class2_k1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_k2 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_k3 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_phi1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_phi2 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_phi3 = new double[data.ndihedraltypes+1];
+
+    data.dihedral_class2_mbt_f1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_mbt_f2 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_mbt_f3 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_mbt_r0 = new double[data.ndihedraltypes+1];
+
+    data.dihedral_class2_ebt_f1_1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_ebt_f2_1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_ebt_f3_1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_ebt_r0_1 = new double[data.ndihedraltypes+1];
+
+    data.dihedral_class2_ebt_f1_2 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_ebt_f2_2 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_ebt_f3_2 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_ebt_r0_2 = new double[data.ndihedraltypes+1];
+
+    data.dihedral_class2_at_f1_1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_at_f2_1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_at_f3_1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_at_theta0_1 = new double[data.ndihedraltypes+1];
+
+    data.dihedral_class2_at_f1_2 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_at_f2_2 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_at_f3_2 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_at_theta0_2 = new double[data.ndihedraltypes+1];
+
+    data.dihedral_class2_aat_k = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_aat_theta0_1 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_aat_theta0_2 = new double[data.ndihedraltypes+1];
+
+    data.dihedral_class2_bb13_k = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_bb13_r10 = new double[data.ndihedraltypes+1];
+    data.dihedral_class2_bb13_r30 = new double[data.ndihedraltypes+1];
+
+    fread(&data.dihedral_class2_k1[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_k2[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_k3[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_phi1[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_phi2[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_phi3[1],sizeof(double),data.ndihedraltypes,fp);
+
+    fread(&data.dihedral_class2_mbt_f1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_mbt_f2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_mbt_f3[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_mbt_r0[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+
+    fread(&data.dihedral_class2_ebt_f1_1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_ebt_f2_1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_ebt_f3_1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_ebt_r0_1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+
+    fread(&data.dihedral_class2_ebt_f1_2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_ebt_f2_2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_ebt_f3_2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_ebt_r0_2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+
+    fread(&data.dihedral_class2_at_f1_1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_at_f2_1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_at_f3_1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_at_theta0_1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+
+    fread(&data.dihedral_class2_at_f1_2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_at_f2_2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_at_f3_2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_at_theta0_2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+
+    fread(&data.dihedral_class2_aat_k[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_aat_theta0_1[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_aat_theta0_2[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+
+    fread(&data.dihedral_class2_bb13_k[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_bb13_r10[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_class2_bb13_r30[1],sizeof(double),
+	  data.ndihedraltypes,fp);
+
+  } else if (strcmp(data.dihedral_style,"harmonic") == 0) {
+
+    data.dihedral_harmonic_k = new double[data.ndihedraltypes+1];
+    data.dihedral_harmonic_multiplicity = new int[data.ndihedraltypes+1];
+    data.dihedral_harmonic_sign = new int[data.ndihedraltypes+1];
+    fread(&data.dihedral_harmonic_k[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_harmonic_multiplicity[1],sizeof(int),
+	  data.ndihedraltypes,fp);
+    fread(&data.dihedral_harmonic_sign[1],sizeof(int),data.ndihedraltypes,fp);
+
+  } else if (strcmp(data.dihedral_style,"helix") == 0) {
+
+    data.dihedral_helix_aphi = new double[data.ndihedraltypes+1];
+    data.dihedral_helix_bphi = new double[data.ndihedraltypes+1];
+    data.dihedral_helix_cphi = new double[data.ndihedraltypes+1];
+    fread(&data.dihedral_helix_aphi[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_helix_bphi[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_helix_cphi[1],sizeof(double),data.ndihedraltypes,fp);
+
+  } else if (strcmp(data.dihedral_style,"multi/harmonic") == 0) {
+
+    data.dihedral_multi_a1 = new double[data.ndihedraltypes+1];
+    data.dihedral_multi_a2 = new double[data.ndihedraltypes+1];
+    data.dihedral_multi_a3 = new double[data.ndihedraltypes+1];
+    data.dihedral_multi_a4 = new double[data.ndihedraltypes+1];
+    data.dihedral_multi_a5 = new double[data.ndihedraltypes+1];
+    fread(&data.dihedral_multi_a1[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_multi_a2[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_multi_a3[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_multi_a4[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_multi_a5[1],sizeof(double),data.ndihedraltypes,fp);
+
+  } else if (strcmp(data.dihedral_style,"opls") == 0) {
+
+    data.dihedral_opls_k1 = new double[data.ndihedraltypes+1];
+    data.dihedral_opls_k2 = new double[data.ndihedraltypes+1];
+    data.dihedral_opls_k3 = new double[data.ndihedraltypes+1];
+    data.dihedral_opls_k4 = new double[data.ndihedraltypes+1];
+    fread(&data.dihedral_opls_k1[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_opls_k2[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_opls_k3[1],sizeof(double),data.ndihedraltypes,fp);
+    fread(&data.dihedral_opls_k4[1],sizeof(double),data.ndihedraltypes,fp);
+
+  } else if (strcmp(data.dihedral_style,"hybrid") == 0) {
+
+    int nstyles = read_int(fp);
+    for (int i = 0; i < nstyles; i++)
+      char *substyle = read_char(fp);
+
+  } else {
+    printf("ERROR: Unknown dihedral style %s\n",data.dihedral_style);
+    exit(1);
+  }
+}
+
+// ---------------------------------------------------------------------
+// improper coeffs
+// one section for each improper style
+// ---------------------------------------------------------------------
+
+void improper(FILE *fp, Data &data)
+{
+  if (strcmp(data.improper_style,"none") == 0) {
+
+  } else if (strcmp(data.improper_style,"class2") == 0) {
+
+    data.improper_class2_k0 = new double[data.nimpropertypes+1];
+    data.improper_class2_chi0 = new double[data.nimpropertypes+1];
+
+    data.improper_class2_aa_k1 = new double[data.nimpropertypes+1];
+    data.improper_class2_aa_k2 = new double[data.nimpropertypes+1];
+    data.improper_class2_aa_k3 = new double[data.nimpropertypes+1];
+    data.improper_class2_aa_theta0_1 = new double[data.nimpropertypes+1];
+    data.improper_class2_aa_theta0_2 = new double[data.nimpropertypes+1];
+    data.improper_class2_aa_theta0_3 = new double[data.nimpropertypes+1];
+
+    fread(&data.improper_class2_k0[1],sizeof(double),
+	  data.nimpropertypes,fp);
+    fread(&data.improper_class2_chi0[1],sizeof(double),
+	  data.nimpropertypes,fp);
+
+    fread(&data.improper_class2_aa_k1[1],sizeof(double),
+	  data.nimpropertypes,fp);
+    fread(&data.improper_class2_aa_k2[1],sizeof(double),
+	  data.nimpropertypes,fp);
+    fread(&data.improper_class2_aa_k3[1],sizeof(double),
+	  data.nimpropertypes,fp);
+    fread(&data.improper_class2_aa_theta0_1[1],sizeof(double),
+	  data.nimpropertypes,fp);
+    fread(&data.improper_class2_aa_theta0_2[1],sizeof(double),
+	  data.nimpropertypes,fp);
+    fread(&data.improper_class2_aa_theta0_3[1],sizeof(double),
+	  data.nimpropertypes,fp);
+
+  } else if (strcmp(data.improper_style,"cvff") == 0) {
+
+    data.improper_cvff_k = new double[data.nimpropertypes+1];
+    data.improper_cvff_sign = new int[data.nimpropertypes+1];
+    data.improper_cvff_multiplicity = new int[data.nimpropertypes+1];
+    fread(&data.improper_cvff_k[1],sizeof(double),data.nimpropertypes,fp);
+    fread(&data.improper_cvff_sign[1],sizeof(int),data.nimpropertypes,fp);
+    fread(&data.improper_cvff_multiplicity[1],sizeof(int),
+	  data.nimpropertypes,fp);
+
+  } else if (strcmp(data.improper_style,"harmonic") == 0) {
+
+    data.improper_harmonic_k = new double[data.nimpropertypes+1];
+    data.improper_harmonic_chi = new double[data.nimpropertypes+1];
+    fread(&data.improper_harmonic_k[1],sizeof(double),data.nimpropertypes,fp);
+    fread(&data.improper_harmonic_chi[1],sizeof(double),
+	  data.nimpropertypes,fp);
+
+  } else if (strcmp(data.improper_style,"hybrid") == 0) {
+
+    int nstyles = read_int(fp);
+    for (int i = 0; i < nstyles; i++)
+      char *substyle = read_char(fp);
+
+  } else {
+    printf("ERROR: Unknown improper style %s\n",data.improper_style);
+    exit(1);
+  }
+}
+
+// ---------------------------------------------------------------------
+// initialize Data
+// ---------------------------------------------------------------------
+
+Data::Data() {}
+
+// ---------------------------------------------------------------------
+// print out stats on problem
+// ---------------------------------------------------------------------
+
+void Data::stats()
+{
+  printf("  Restart file version = %s\n",version);
+  printf("  Ntimestep = %d\n",ntimestep);
+  printf("  Nprocs = %d\n",nprocs);
+  printf("  Natoms = %d\n",natoms);
+  printf("  Nbonds = %d\n",nbonds);
+  printf("  Nangles = %d\n",nangles);
+  printf("  Ndihedrals = %d\n",ndihedrals);
+  printf("  Nimpropers = %d\n",nimpropers);
+  printf("  Unit style = %s\n",unit_style);
+  printf("  Atom style = %s\n",atom_style);
+  printf("  Pair style = %s\n",pair_style);
+  printf("  Bond style = %s\n",bond_style);
+  printf("  Angle style = %s\n",angle_style);
+  printf("  Dihedral style = %s\n",dihedral_style);
+  printf("  Improper style = %s\n",improper_style);
+  printf("  Xlo xhi = %g %g\n",xlo,xhi);
+  printf("  Ylo yhi = %g %g\n",ylo,yhi);
+  printf("  Zlo zhi = %g %g\n",zlo,zhi);
+  if (triclinic) printf("  Xy xz yz = %g %g %g\n",xy,xz,yz);
+  printf("  Periodicity = %d %d %d\n",xperiodic,yperiodic,zperiodic);
+  printf("  Boundary = %d %d, %d %d, %d %d\n",boundary[0][0],boundary[0][1],
+	 boundary[1][0],boundary[1][1],boundary[2][0],boundary[2][1]);
+}
+
+// ---------------------------------------------------------------------
+// write the data file
+// ---------------------------------------------------------------------
+
+void Data::write(FILE *fp, FILE *fp2)
+{
+  fprintf(fp,"LAMMPS data file from restart file: timestep = %d, procs = %d\n",
+	  ntimestep,nprocs);
+  
+  fprintf(fp,"\n");
+
+  fprintf(fp,"%d atoms\n",natoms);
+  if (nbonds) fprintf(fp,"%d bonds\n",nbonds);
+  if (nangles) fprintf(fp,"%d angles\n",nangles);
+  if (ndihedrals) fprintf(fp,"%d dihedrals\n",ndihedrals);
+  if (nimpropers) fprintf(fp,"%d impropers\n",nimpropers);
+
+  fprintf(fp,"\n");
+
+  fprintf(fp,"%d atom types\n",ntypes);
+  if (nbondtypes) fprintf(fp,"%d bond types\n",nbondtypes);
+  if (nangletypes) fprintf(fp,"%d angle types\n",nangletypes);
+  if (ndihedraltypes) fprintf(fp,"%d dihedral types\n",ndihedraltypes);
+  if (nimpropertypes) fprintf(fp,"%d improper types\n",nimpropertypes);
+
+  fprintf(fp,"\n");
+  
+  fprintf(fp,"%-1.16e %-1.16e xlo xhi\n",xlo,xhi);
+  fprintf(fp,"%-1.16e %-1.16e ylo yhi\n",ylo,yhi);
+  fprintf(fp,"%-1.16e %-1.16e zlo zhi\n",zlo,zhi);
+  if (triclinic) fprintf(fp,"%-1.16e %-1.16e %-1.16e xy xz yz\n",xy,xz,yz);
+
+  // write ff styles to input file
+
+  if (fp2) {
+    fprintf(fp2,"# LAMMPS input file from restart file: timestep = %d, procs = %d\n\n",
+	  ntimestep,nprocs);
+    if (pair_style) fprintf(fp2,"pair_style %s\n",pair_style);
+    if (bond_style) fprintf(fp2,"bond_style %s\n",bond_style);
+    if (angle_style) fprintf(fp2,"angle_style %s\n",angle_style);
+    if (dihedral_style) fprintf(fp2,"dihedral_style %s\n",dihedral_style);
+    if (improper_style) fprintf(fp2,"improper_style %s\n",improper_style);
+    fprintf(fp2,"special_bonds %g %g %g %g %g %g\n",
+            special_lj[1],special_lj[2],special_lj[3],
+            special_lj[1],special_coul[2],special_coul[3]);
+    fprintf(fp2,"\n");
+  }
+
+  // mass to either data file or input file
+
+  if (mass) {
+    if (fp2) {
+      fprintf(fp2,"\n");
+      for (int i = 1; i <= ntypes; i++) fprintf(fp2,"mass %d %g\n",i,mass[i]);
+      fprintf(fp2,"\n");
+    } else {
+      fprintf(fp,"\nMasses\n\n");
+      for (int i = 1; i <= ntypes; i++) fprintf(fp,"%d %g\n",i,mass[i]);
+    }
+  }
+
+  // shape and dipole to data file
+  // convert shape from radius to diameter
+
+  if (shape) {
+    fprintf(fp,"\nShapes\n\n");
+    for (int i = 1; i <= ntypes; i++)
+      fprintf(fp,"%d %g %g %g\n",i,
+	      2.0*shape[3*i+0],2.0*shape[3*i+1],2.0*shape[3*i+2]);
+  }
+
+  if (dipole) {
+    fprintf(fp,"\nDipoles\n\n");
+    for (int i = 1; i <= ntypes; i++) fprintf(fp,"%d %g\n",i,dipole[i]);
+  }
+
+  // pair coeffs to data file
+
+  if (pair_style && fp2 == NULL) {
+    if ((strcmp(pair_style,"none") != 0) && 
+	(strcmp(pair_style,"airebo") != 0) &&
+	(strcmp(pair_style,"coul/cut") != 0) &&
+	(strcmp(pair_style,"coul/debye") != 0) &&
+        (strcmp(pair_style,"eff/cut") != 0) &&
+	(strcmp(pair_style,"coul/long") != 0) &&
+	(strcmp(pair_style,"eam") != 0) &&
+	(strcmp(pair_style,"eam/opt") != 0) &&
+	(strcmp(pair_style,"eam/alloy") != 0) &&
+	(strcmp(pair_style,"eam/alloy/opt") != 0) &&
+	(strcmp(pair_style,"eam/fs") != 0) &&
+	(strcmp(pair_style,"eam/fs/opt") != 0) &&
+	(strcmp(pair_style,"gran/history") != 0) &&
+	(strcmp(pair_style,"gran/no_history") != 0) &&
+	(strcmp(pair_style,"gran/hertzian") != 0) &&
+	(strcmp(pair_style,"meam") != 0) &&
+	(strcmp(pair_style,"reax") != 0) &&
+	(strcmp(pair_style,"sw") != 0) &&
+	(strcmp(pair_style,"table") != 0) &&
+	(strcmp(pair_style,"tersoff") != 0) &&
+	(strcmp(pair_style,"hybrid") != 0) &&
+	(strcmp(pair_style,"hybrid/overlay") != 0))
+      fprintf(fp,"\nPair Coeffs\n\n");
+    
+    if (strcmp(pair_style,"born/coul/long") == 0) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g %g %g %g\n",i,
+		pair_born_A[i],pair_born_rho[i],pair_born_sigma[i],
+		pair_born_C[i],pair_born_D[i]);
+
+    } else if ((strcmp(pair_style,"buck") == 0) || 
+	(strcmp(pair_style,"buck/coul/cut") == 0) ||
+	(strcmp(pair_style,"buck/coul/long") == 0) ||
+	(strcmp(pair_style,"buck/long") == 0)) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g %g\n",i,
+		pair_buck_A[i],pair_buck_rho[i],pair_buck_C[i]);
+      
+    } else if (strcmp(pair_style,"colloid") == 0) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g %g %g\n",i,
+		pair_colloid_A12[i],pair_colloid_sigma[i],
+		pair_colloid_d2[i],pair_colloid_d2[i]);
+
+    } else if (strcmp(pair_style,"dipole/cut") == 0) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		pair_dipole_epsilon[i],pair_dipole_sigma[i]);
+      
+    } else if (strcmp(pair_style,"dpd") == 0) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		pair_dpd_a0[i],pair_dpd_gamma[i]);
+      
+    } else if (strcmp(pair_style,"gayberne") == 0) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g %g %g %g %g %g %g\n",i,
+		pair_gb_epsilon[i],pair_gb_sigma[i],
+		pair_gb_epsa[i],pair_gb_epsb[i],pair_gb_epsc[i],
+		pair_gb_epsa[i],pair_gb_epsb[i],pair_gb_epsc[i]);
+
+    } else if ((strcmp(pair_style,"lj/charmm/coul/charmm") == 0) ||
+	       (strcmp(pair_style,"lj/charmm/coul/charmm/implicit") == 0) ||
+	       (strcmp(pair_style,"lj/charmm/coul/long") == 0) ||
+	       (strcmp(pair_style,"lj/charmm/coul/long") == 0)) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g %g %g\n",i,
+		pair_charmm_epsilon[i],pair_charmm_sigma[i],
+		pair_charmm_eps14[i],pair_charmm_sigma14[i]);
+      
+    } else if ((strcmp(pair_style,"lj/class2") == 0) || 
+	       (strcmp(pair_style,"lj/class2/coul/cut") == 0) ||
+	       (strcmp(pair_style,"lj/class2/coul/long") == 0)) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		pair_class2_epsilon[i],pair_class2_sigma[i]);
+      
+    } else if ((strcmp(pair_style,"lj/cut") == 0) || 
+	       (strcmp(pair_style,"lj/cut/opt") == 0) ||
+	       (strcmp(pair_style,"lj/cut/coul/cut") == 0) ||
+	       (strcmp(pair_style,"lj/cut/coul/debye") == 0) ||
+	       (strcmp(pair_style,"lj/cut/coul/long") == 0) ||
+	       (strcmp(pair_style,"lj/cut/coul/long/tip4p") == 0) |
+	       (strcmp(pair_style,"lj/coul") == 0)) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		pair_lj_epsilon[i],pair_lj_sigma[i]);
+
+    } else if (strcmp(pair_style,"lj/expand") == 0) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g %g\n",i,
+		pair_ljexpand_epsilon[i],pair_ljexpand_sigma[i],
+		pair_ljexpand_shift[i]);
+
+    } else if ((strcmp(pair_style,"lj/gromacs") == 0) ||
+	       (strcmp(pair_style,"lj/gromacs/coul/gromacs") == 0)) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		pair_ljgromacs_epsilon[i],pair_ljgromacs_sigma[i]);
+
+    } else if (strcmp(pair_style,"lj/smooth") == 0) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		pair_ljsmooth_epsilon[i],pair_ljsmooth_sigma[i]);
+      
+    } else if ((strcmp(pair_style,"morse") == 0) ||
+	       (strcmp(pair_style,"morse/opt") == 0)) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g %g\n",i,
+		pair_morse_d0[i],pair_morse_alpha[i],pair_morse_r0[i]);
+      
+    } else if (strcmp(pair_style,"soft") == 0) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		pair_soft_start[i],pair_soft_stop[i]);
+      
+    } else if (strcmp(pair_style,"yukawa") == 0) {
+      for (int i = 1; i <= ntypes; i++)
+	fprintf(fp,"%d %g\n",i,
+		pair_yukawa_A[i]);
+
+    } else if ((strcmp(pair_style,"cg/cmm") == 0) || 
+               (strcmp(pair_style,"cg/cmm/coul/cut") == 0) ||
+               (strcmp(pair_style,"cg/cmm/coul/long") == 0)) {
+      printf("ERROR: Cannot write pair_style %s to data file\n",
+	     pair_style);
+      exit(1);
+    }
+  }
+
+  // pair coeffs to input file
+  // only supported styles = cg/cmm
+
+  if (pair_style && fp2) {
+    if ((strcmp(pair_style,"cg/cmm") == 0) || 
+	(strcmp(pair_style,"cg/cmm/coul/cut") == 0) ||
+	(strcmp(pair_style,"cg/cmm/coul/long") == 0)) {
+      for (int i = 1; i <= ntypes; i++) {
+	for (int j = i; j <= ntypes; j++) {
+	  fprintf(fp2,"pair_coeff %d %d %s %g %g\n",i,j,
+		  cg_type_list[pair_cg_cmm_type[i][j]],
+		  pair_cg_epsilon[i][j],pair_cg_sigma[i][j]);
+	}
+      }
+
+    } else {
+      printf("ERROR: Cannot write pair_style %s to input file\n",
+	     pair_style);
+      exit(1);
+    }
+  }
+
+  // bond coeffs to data file
+
+  if (bond_style && fp2 == NULL) {
+    if ((strcmp(bond_style,"none") != 0) && 
+	(strcmp(bond_style,"hybrid") != 0))
+      fprintf(fp,"\nBond Coeffs\n\n");
+
+    if (strcmp(bond_style,"class2") == 0) {
+      for (int i = 1; i <= nbondtypes; i++)
+	fprintf(fp,"%d %g %g %g %g\n",i,
+		bond_class2_r0[i],bond_class2_k2[i],
+		bond_class2_k3[i],bond_class2_k4[i]);
+
+    } else if (strcmp(bond_style,"fene") == 0) {
+      for (int i = 1; i <= nbondtypes; i++)
+	fprintf(fp,"%d %g %g %g %g\n",i,
+		bond_fene_k[i],bond_fene_r0[i],
+		bond_fene_epsilon[i],bond_fene_sigma[i]);
+
+    } else if (strcmp(bond_style,"fene/expand") == 0) {
+      for (int i = 1; i <= nbondtypes; i++)
+	fprintf(fp,"%d %g %g %g %g %g\n",i,
+		bond_feneexpand_k[i],bond_feneexpand_r0[i],
+		bond_feneexpand_epsilon[i],bond_feneexpand_sigma[i],
+		bond_feneexpand_shift[i]);
+
+    } else if (strcmp(bond_style,"harmonic") == 0) {
+      for (int i = 1; i <= nbondtypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		bond_harmonic_k[i],bond_harmonic_r0[i]);
+
+    } else if (strcmp(bond_style,"morse") == 0) {
+      for (int i = 1; i <= nbondtypes; i++)
+	fprintf(fp,"%d %g %g %g\n",i,
+		bond_morse_d0[i],bond_morse_alpha[i],bond_morse_r0[i]);
+
+    } else if (strcmp(bond_style,"nonlinear") == 0) {
+      for (int i = 1; i <= nbondtypes; i++)
+	fprintf(fp,"%d %g %g %g\n",i,
+		bond_nonlinear_epsilon[i],bond_nonlinear_r0[i],
+		bond_nonlinear_lamda[i]);
+
+    } else if (strcmp(bond_style,"quartic") == 0) {
+      for (int i = 1; i <= nbondtypes; i++)
+        fprintf(fp,"%d %g %g %g %g %g\n",i,
+                bond_quartic_k[i],bond_quartic_b1[i],bond_quartic_b2[i],
+                bond_quartic_rc[i],bond_quartic_u0[i]);
+    }
+  }
+
+  // bond coeffs to input file
+  // only supported styles = harmonic
+
+  if (bond_style && fp2) {
+    if (strcmp(bond_style,"harmonic") == 0) {
+      for (int i = 1; i <= nbondtypes; i++)
+	fprintf(fp2,"bond_coeff  %d %g %g\n",i,
+		bond_harmonic_k[i],bond_harmonic_r0[i]);
+
+    } else {
+      printf("ERROR: Cannot write bond_style %s to input file\n",
+	     bond_style);
+      exit(1);
+    }
+  }
+
+  // angle coeffs to data file
+
+  if (angle_style && fp2 == NULL) {
+    double PI = 3.1415926;           // convert back to degrees
+
+    if ((strcmp(angle_style,"none") != 0) && 
+	(strcmp(angle_style,"hybrid") != 0))
+      fprintf(fp,"\nAngle Coeffs\n\n");
+    
+    if (strcmp(angle_style,"charmm") == 0) {
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp,"%d %g %g %g %g\n",i,
+		angle_charmm_k[i],angle_charmm_theta0[i]/PI*180.0,
+		angle_charmm_k_ub[i],angle_charmm_r_ub[i]);
+
+    } else if (strcmp(angle_style,"class2") == 0) {
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp,"%d %g %g %g %g\n",i,
+		angle_class2_theta0[i]/PI*180.0,angle_class2_k2[i], 
+		angle_class2_k3[i],angle_class2_k4[i]);
+      
+      fprintf(fp,"\nBondBond Coeffs\n\n");
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp,"%d %g %g %g\n",i,
+		angle_class2_bb_k[i],
+		angle_class2_bb_r1[i],angle_class2_bb_r2[i]);
+      
+      fprintf(fp,"\nBondAngle Coeffs\n\n");
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp,"%d %g %g %g %g\n",i,
+		angle_class2_ba_k1[i],angle_class2_ba_k2[i],
+		angle_class2_ba_r1[i],angle_class2_ba_r2[i]);
+
+    } else if (strcmp(angle_style,"cosine") == 0) {
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp,"%d %g\n",i,angle_cosine_k[i]);
+      
+    } else if ((strcmp(angle_style,"cosine/squared") == 0) ||
+               (strcmp(angle_style,"cosine/delta") == 0)) {
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		angle_cosine_squared_k[i],
+		angle_cosine_squared_theta0[i]/PI*180.0);
+      
+    } else if (strcmp(angle_style,"harmonic") == 0) {
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0);
+      
+    } else if (strcmp(angle_style,"cg/cmm") == 0) {
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp,"%d %g %g %s %g %g\n",i,
+		angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0,
+		cg_type_list[angle_cg_cmm_type[i]],angle_cg_cmm_epsilon[i],
+		angle_cg_cmm_sigma[i]);
+    }
+  }
+
+  // angle coeffs to input file
+  // only supported styles = cosine/squared, harmonic, cg/cmm
+
+  if (angle_style && fp2) {
+    double PI = 3.1415926;           // convert back to degrees
+
+    if ((strcmp(angle_style,"cosine/squared") == 0) ||
+        (strcmp(angle_style,"cosine/delta") == 0)) {
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp2,"angle_coeffs  %d %g %g\n",i,
+		angle_cosine_squared_k[i],
+		angle_cosine_squared_theta0[i]/PI*180.0);
+      
+    } else if (strcmp(angle_style,"harmonic") == 0) {
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp2,"angle_coeffs  %d %g %g\n",i,
+		angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0);
+      
+    } else if (strcmp(angle_style,"cg/cmm") == 0) {
+      for (int i = 1; i <= nangletypes; i++)
+	fprintf(fp2,"angle_coeffs  %d %g %g %s %g %g\n",i,
+		angle_harmonic_k[i],angle_harmonic_theta0[i]/PI*180.0,
+		cg_type_list[angle_cg_cmm_type[i]],angle_cg_cmm_epsilon[i],
+		angle_cg_cmm_sigma[i]);
+      
+    } else {
+      printf("ERROR: Cannot write angle_style %s to input file\n",
+	     angle_style);
+      exit(1);
+    }
+  }
+
+  if (dihedral_style) {
+    double PI = 3.1415926;           // convert back to degrees
+
+    if ((strcmp(dihedral_style,"none") != 0) && 
+	(strcmp(dihedral_style,"hybrid") != 0))
+      fprintf(fp,"\nDihedral Coeffs\n\n");
+
+    if (strcmp(dihedral_style,"charmm") == 0) {
+      for (int i = 1; i <= ndihedraltypes; i++)
+	fprintf(fp,"%d %g %d %d %g\n",i,
+		dihedral_charmm_k[i],dihedral_charmm_multiplicity[i],
+		dihedral_charmm_sign[i],dihedral_charmm_weight[i]);
+
+    } else if (strcmp(dihedral_style,"class2") == 0) {
+      for (int i = 1; i <= ndihedraltypes; i++)
+	fprintf(fp,"%d %g %g %g %g %g %g\n",i,
+		dihedral_class2_k1[i],
+		dihedral_class2_phi1[i]/PI*180.0,
+		dihedral_class2_k2[i],
+		dihedral_class2_phi2[i]/PI*180.0,
+		dihedral_class2_k3[i],
+		dihedral_class2_phi3[i]/PI*180.0);
+
+      fprintf(fp,"\nMiddleBondTorsion Coeffs\n\n");
+      for (int i = 1; i <= ndihedraltypes; i++)
+	fprintf(fp,"%d %g %g %g %g\n",i,
+		dihedral_class2_mbt_f1[i],dihedral_class2_mbt_f2[i],
+		dihedral_class2_mbt_f3[i],dihedral_class2_mbt_r0[i]);
+
+      fprintf(fp,"\nEndBondTorsion Coeffs\n\n");
+      for (int i = 1; i <= ndihedraltypes; i++)
+	fprintf(fp,"%d %g %g %g %g %g %g %g %g\n",i,
+		dihedral_class2_ebt_f1_1[i],dihedral_class2_ebt_f2_1[i],
+		dihedral_class2_ebt_f3_1[i],
+		dihedral_class2_ebt_f1_2[i],dihedral_class2_ebt_f2_2[i],
+		dihedral_class2_ebt_f3_2[i],
+		dihedral_class2_ebt_r0_1[i],
+		dihedral_class2_ebt_r0_2[i]);
+
+      fprintf(fp,"\nAngleTorsion Coeffs\n\n");
+      for (int i = 1; i <= ndihedraltypes; i++)
+	fprintf(fp,"%d %g %g %g %g %g %g %g %g\n",i,
+		dihedral_class2_at_f1_1[i],dihedral_class2_at_f2_1[i],
+		dihedral_class2_at_f3_1[i],
+		dihedral_class2_at_f1_2[i],dihedral_class2_at_f2_2[i],
+		dihedral_class2_at_f3_2[i],
+		dihedral_class2_at_theta0_1[i]/PI*180.0,
+		dihedral_class2_at_theta0_2[i]/PI*180.0);
+
+      fprintf(fp,"\nAngleAngleTorsion Coeffs\n\n");
+      for (int i = 1; i <= ndihedraltypes; i++)
+	fprintf(fp,"%d %g %g %g\n",i,
+		dihedral_class2_aat_k[i],
+		dihedral_class2_aat_theta0_1[i]/PI*180.0,
+		dihedral_class2_aat_theta0_2[i]/PI*180.0);
+
+      fprintf(fp,"\nBondBond13 Coeffs\n\n");
+      for (int i = 1; i <= ndihedraltypes; i++)
+	fprintf(fp,"%d %g %g %g\n",i,
+		dihedral_class2_bb13_k[i],
+		dihedral_class2_bb13_r10[i],dihedral_class2_bb13_r30[i]);
+
+    } else if (strcmp(dihedral_style,"harmonic") == 0) {
+      for (int i = 1; i <= ndihedraltypes; i++)
+	fprintf(fp,"%d %g %d %d\n",i,
+		dihedral_harmonic_k[i],dihedral_harmonic_multiplicity[i],
+		dihedral_harmonic_sign[i]);
+
+    } else if (strcmp(dihedral_style,"helix") == 0) {
+      for (int i = 1; i <= ndihedraltypes; i++) 
+	fprintf(fp,"%d %g %g %g\n",i,dihedral_helix_aphi[i],
+		dihedral_helix_bphi[i],dihedral_helix_cphi[i]);
+
+    } else if (strcmp(dihedral_style,"multi/harmonic") == 0) {
+      for (int i = 1; i <= ndihedraltypes; i++)
+	fprintf(fp,"%d %g %g %g %g %g\n",i,
+		dihedral_multi_a1[i],dihedral_multi_a2[i],
+		dihedral_multi_a3[i],dihedral_multi_a4[i],
+		dihedral_multi_a5[i]);
+
+    } else if (strcmp(dihedral_style,"opls") == 0) {
+      for (int i = 1; i <= ndihedraltypes; i++)        // restore factor of 2
+	fprintf(fp,"%d %g %g %g %g\n",i,
+		2.0*dihedral_opls_k1[i],2.0*dihedral_opls_k2[i],
+		2.0*dihedral_opls_k3[i],2.0*dihedral_opls_k4[i]);
+    }
+  }
+
+  if (improper_style) {
+    double PI = 3.1415926;           // convert back to degrees
+
+    if ((strcmp(improper_style,"none") != 0) && 
+	(strcmp(improper_style,"hybrid") != 0))
+      fprintf(fp,"\nImproper Coeffs\n\n");
+
+    if (strcmp(improper_style,"class2") == 0) {
+      for (int i = 1; i <= nimpropertypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		improper_class2_k0[i],improper_class2_chi0[i]/PI*180.0);
+
+      fprintf(fp,"\nAngleAngle Coeffs\n\n");
+      for (int i = 1; i <= nimpropertypes; i++)
+	fprintf(fp,"%d %g %g %g %g %g %g\n",i,
+		improper_class2_aa_k1[i],improper_class2_aa_k2[i],
+		improper_class2_aa_k3[i],
+		improper_class2_aa_theta0_1[i]/PI*180.0,
+		improper_class2_aa_theta0_2[i]/PI*180.0,
+		improper_class2_aa_theta0_3[i]/PI*180.0);
+
+    } else if (strcmp(improper_style,"cvff") == 0) {
+      for (int i = 1; i <= nimpropertypes; i++)
+	fprintf(fp,"%d %g %d %d\n",i,
+		improper_cvff_k[i],improper_cvff_sign[i],
+		improper_cvff_multiplicity[i]);
+
+    } else if (strcmp(improper_style,"harmonic") == 0) {
+      for (int i = 1; i <= nimpropertypes; i++)
+	fprintf(fp,"%d %g %g\n",i,
+		improper_harmonic_k[i],improper_harmonic_chi[i]/PI*180.0);
+    }
+  }
+
+  if (natoms) {
+    fprintf(fp,"\nAtoms\n\n");
+
+    int ix,iy,iz;
+    for (int i = 0; i < natoms; i++) {
+
+      ix = (image[i] & 1023) - 512;
+      iy = (image[i] >> 10 & 1023) - 512;
+      iz = (image[i] >> 20) - 512;
+
+      if (style_hybrid == 0) {
+	if (style_angle) write_atom_angle(fp,i,ix,iy,iz);
+	if (style_atomic) write_atom_atomic(fp,i,ix,iy,iz);
+	if (style_bond) write_atom_bond(fp,i,ix,iy,iz);
+	if (style_charge) write_atom_charge(fp,i,ix,iy,iz);
+	if (style_dipole) write_atom_dipole(fp,i,ix,iy,iz);
+	if (style_dpd) write_atom_dpd(fp,i,ix,iy,iz);
+	if (style_ellipsoid) write_atom_ellipsoid(fp,i,ix,iy,iz);
+	if (style_full) write_atom_full(fp,i,ix,iy,iz);
+	if (style_granular) write_atom_granular(fp,i,ix,iy,iz);
+	if (style_molecular) write_atom_molecular(fp,i,ix,iy,iz);
+	if (style_peri) write_atom_peri(fp,i,ix,iy,iz);
+        if (style_electron) write_atom_electron(fp,i,ix,iy,iz);
+	fprintf(fp,"\n");
+
+      } else {
+	fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e",
+		tag[i],type[i],x[i],y[i],z[i]);
+	for (int k = 1; k <= style_hybrid; k++) {
+	  if (k == style_angle) write_atom_angle_extra(fp,i);
+	  if (k == style_atomic) write_atom_atomic_extra(fp,i);
+	  if (k == style_bond) write_atom_bond_extra(fp,i);
+	  if (k == style_charge) write_atom_charge_extra(fp,i);
+	  if (k == style_dipole) write_atom_dipole_extra(fp,i);
+	  if (k == style_dpd) write_atom_dpd_extra(fp,i);
+	  if (k == style_ellipsoid) write_atom_ellipsoid_extra(fp,i);
+	  if (k == style_full) write_atom_full_extra(fp,i);
+	  if (k == style_granular) write_atom_granular_extra(fp,i);
+	  if (k == style_molecular) write_atom_molecular_extra(fp,i);
+	  if (k == style_peri) write_atom_peri_extra(fp,i);
+          if (k == style_electron) write_atom_electron_extra(fp,i);
+	}
+	fprintf(fp," %d %d %d\n",ix,iy,iz);
+      }
+    }
+  }
+
+  if (natoms) {
+    fprintf(fp,"\nVelocities\n\n");
+    for (int i = 0; i < natoms; i++)
+
+      if (style_hybrid == 0) {
+	if (style_angle) write_vel_angle(fp,i);
+	if (style_atomic) write_vel_atomic(fp,i);
+	if (style_bond) write_vel_bond(fp,i);
+	if (style_charge) write_vel_charge(fp,i);
+	if (style_dipole) write_vel_dipole(fp,i);
+	if (style_dpd) write_vel_dpd(fp,i);
+	if (style_ellipsoid) write_vel_ellipsoid(fp,i);
+	if (style_full) write_vel_full(fp,i);
+	if (style_granular) write_vel_granular(fp,i);
+	if (style_molecular) write_vel_molecular(fp,i);
+	if (style_peri) write_vel_peri(fp,i);
+        if (style_electron) write_vel_electron(fp,i);
+	fprintf(fp,"\n");
+
+      } else {
+	fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+	for (int k = 1; k <= style_hybrid; k++) {
+	  if (k == style_angle) write_vel_angle_extra(fp,i);
+	  if (k == style_atomic) write_vel_atomic_extra(fp,i);
+	  if (k == style_bond) write_vel_bond_extra(fp,i);
+	  if (k == style_charge) write_vel_charge_extra(fp,i);
+	  if (k == style_dipole) write_vel_dipole_extra(fp,i);
+	  if (k == style_dpd) write_vel_dpd_extra(fp,i);
+	  if (k == style_ellipsoid) write_vel_ellipsoid_extra(fp,i);
+	  if (k == style_full) write_vel_full_extra(fp,i);
+	  if (k == style_granular) write_vel_granular_extra(fp,i);
+	  if (k == style_molecular) write_vel_molecular_extra(fp,i);
+	  if (k == style_peri) write_vel_peri_extra(fp,i);
+          if (k == style_electron) write_vel_electron_extra(fp,i);
+	}
+	fprintf(fp,"\n");
+      }
+  }
+
+  if (nbonds) {
+    fprintf(fp,"\nBonds\n\n");
+    for (int i = 0; i < nbonds; i++)
+      fprintf(fp,"%d %d %d %d\n",
+	      i+1,bond_type[i],bond_atom1[i],bond_atom2[i]);
+  }
+
+  if (nangles) {
+    fprintf(fp,"\nAngles\n\n");
+    for (int i = 0; i < nangles; i++)
+      fprintf(fp,"%d %d %d %d %d\n",
+	      i+1,angle_type[i],angle_atom1[i],angle_atom2[i],angle_atom3[i]);
+  }
+
+  if (ndihedrals) {
+    fprintf(fp,"\nDihedrals\n\n");
+    for (int i = 0; i < ndihedrals; i++)
+      fprintf(fp,"%d %d %d %d %d %d\n",
+	      i+1,dihedral_type[i],dihedral_atom1[i],dihedral_atom2[i],
+	      dihedral_atom3[i],dihedral_atom4[i]);
+  }
+
+  if (nimpropers) {
+    fprintf(fp,"\nImpropers\n\n");
+    for (int i = 0; i < nimpropers; i++)
+      fprintf(fp,"%d %d %d %d %d %d\n",
+	      i+1,improper_type[i],improper_atom1[i],improper_atom2[i],
+	      improper_atom3[i],improper_atom4[i]);
+  }
+}
+
+// ---------------------------------------------------------------------
+// per-atom write routines
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void Data::write_atom_angle(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],molecule[i],type[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_atomic(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],type[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_bond(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],molecule[i],type[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_charge(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],type[i],q[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_dipole(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],type[i],q[i],x[i],y[i],z[i],mux[i],muy[i],muz[i],ix,iy,iz);
+}
+
+void Data::write_atom_dpd(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],type[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_ellipsoid(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],type[i],x[i],y[i],z[i],
+	  quatw[i],quati[i],quatj[i],quatk[i],ix,iy,iz);
+}
+
+void Data::write_atom_full(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],molecule[i],type[i],q[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_granular(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],type[i],2.0*radius[i],density[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_molecular(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %d %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],molecule[i],type[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_peri(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+	  tag[i],type[i],vfrac[i],density[i],rmass[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+void Data::write_atom_electron(FILE *fp, int i, int ix, int iy, int iz)
+{
+  fprintf(fp,"%d %d %-1.16e %d %-1.16e %-1.16e %-1.16e %-1.16e %d %d %d",
+          tag[i],type[i],q[i],spin[i],radius[i],x[i],y[i],z[i],ix,iy,iz);
+}
+
+// ---------------------------------------------------------------------
+// per-atom write routines of extra quantities unique to style
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void Data::write_atom_angle_extra(FILE *fp, int i)
+{
+  fprintf(fp," %d",molecule[i]);
+}
+
+void Data::write_atom_atomic_extra(FILE *fp, int i) {}
+
+void Data::write_atom_bond_extra(FILE *fp, int i)
+{
+  fprintf(fp," %d",molecule[i]);
+}
+
+void Data::write_atom_charge_extra(FILE *fp, int i)
+{
+  fprintf(fp," %-1.16e",q[i]);
+}
+
+void Data::write_atom_dipole_extra(FILE *fp, int i)
+{
+  fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e",q[i],mux[i],muy[i],muz[i]);
+}
+
+void Data::write_atom_dpd_extra(FILE *fp, int i) {}
+
+void Data::write_atom_ellipsoid_extra(FILE *fp, int i)
+{
+  fprintf(fp," %-1.16e %-1.16e %-1.16e %-1.16e",quatw[i],quati[i],quatj[i],quatk[i]);
+}
+
+void Data::write_atom_full_extra(FILE *fp, int i)
+{
+  fprintf(fp," %d %-1.16e",molecule[i],q[i]);
+}
+
+void Data::write_atom_granular_extra(FILE *fp, int i)
+{
+  fprintf(fp," %-1.16e %-1.16e",2.0*radius[i],density[i]);
+}
+
+void Data::write_atom_molecular_extra(FILE *fp, int i)
+{
+  fprintf(fp," %d",molecule[i]);
+}
+
+void Data::write_atom_peri_extra(FILE *fp, int i)
+{
+  fprintf(fp," %-1.16e %-1.16e %-1.16e",vfrac[i],density[i],rmass[i]);
+}
+
+void Data::write_atom_electron_extra(FILE *fp, int i)
+{
+  fprintf(fp," %-1.16e %d %-1.16e",q[i],spin[i],radius[i]);
+}
+
+// ---------------------------------------------------------------------
+// per-atom velocity write routines
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void Data::write_vel_angle(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_atomic(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_bond(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_charge(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_dipole(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_dpd(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_ellipsoid(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",
+	  tag[i],vx[i],vy[i],vz[i],angmomx[i],angmomy[i],angmomz[i]);
+}
+
+void Data::write_vel_full(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_granular(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e %-1.16e",
+	  tag[i],vx[i],vy[i],vz[i],omegax[i],omegay[i],omegaz[i]);
+}
+
+void Data::write_vel_molecular(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_peri(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i]);
+}
+
+void Data::write_vel_electron(FILE *fp, int i)
+{
+  fprintf(fp,"%d %-1.16e %-1.16e %-1.16e %-1.16e",tag[i],vx[i],vy[i],vz[i],ervel[i]);
+}
+
+// ---------------------------------------------------------------------
+// per-atom velocity write routines of extra quantities unique to style
+// one routine per atom style
+// ---------------------------------------------------------------------
+
+void Data::write_vel_angle_extra(FILE *fp, int i) {}
+void Data::write_vel_atomic_extra(FILE *fp, int i) {}
+void Data::write_vel_bond_extra(FILE *fp, int i) {}
+void Data::write_vel_charge_extra(FILE *fp, int i) {}
+void Data::write_vel_dipole_extra(FILE *fp, int i) {}
+void Data::write_vel_dpd_extra(FILE *fp, int i) {}
+
+void Data::write_vel_ellipsoid_extra(FILE *fp, int i)
+{
+  fprintf(fp," %-1.16e %-1.16e %-1.16e",angmomx[i],angmomy[i],angmomz[i]);
+}
+
+void Data::write_vel_full_extra(FILE *fp, int i) {}
+
+void Data::write_vel_granular_extra(FILE *fp, int i)
+{
+  fprintf(fp," %-1.16e %-1.16e %-1.16e",omegax[i],omegay[i],omegaz[i]);
+}
+
+void Data::write_vel_molecular_extra(FILE *fp, int i) {}
+void Data::write_vel_peri_extra(FILE *fp, int i) {}
+void Data::write_vel_electron_extra(FILE *fp, int i) 
+{
+  fprintf(fp," %-1.16e",ervel[i]);
+}
+// ---------------------------------------------------------------------
+// binary reads from restart file
+// ---------------------------------------------------------------------
+
+int read_int(FILE *fp)
+{
+  int value;
+  fread(&value,sizeof(int),1,fp);
+  return value;
+}
+
+double read_double(FILE *fp)
+{
+  double value;
+  fread(&value,sizeof(double),1,fp);
+  return value;
+}
+
+char *read_char(FILE *fp)
+{
+  int n;
+  fread(&n,sizeof(int),1,fp);
+  if (n == 0) return NULL;
+  char *value = new char[n];
+  fread(value,sizeof(char),n,fp);
+  return value;
+}
diff --git a/tools/eff/structure_generators/Li-solid-angstrom.pl b/tools/eff/structure_generators/Li-solid-angstrom.pl
new file mode 100644
index 0000000000..a0d7bdb533
--- /dev/null
+++ b/tools/eff/structure_generators/Li-solid-angstrom.pl
@@ -0,0 +1,92 @@
+#!/usr/bin/perl
+
+# Usage: Li-solid.pl <nx> <ny> <nz> > data.Li
+#
+# Generates a lithium solid LAMMPS data input file
+# with FCC nuclear positions and interstitial electron positions.
+
+$nx = shift(@ARGV);
+$ny = shift(@ARGV);
+$nz = shift(@ARGV);
+
+$L = 4.419688383648;  # eFF optimized (8.32 expt)
+
+# This part changes for different lattices
+@xunit = (0, 0.2645886245, 0, 0.2645886245, 0.2645886245, 0, 0, 0.2645886245);
+@yunit = (0, 0.2645886245, 0.2645886245, 0, 0, 0.2645886245, 0, 0.2645886245);
+@zunit = (0, 0, 0.2645886245, 0.2645886245, 0, 0, 0.2645886245, 0.2645886245);
+
+$r_elec = 0.3704240743;
+$r_elec2 = 1.587531747;
+
+$Lx = $L;
+$Ly = $L;
+$Lz = $L;
+
+$idx = 0;
+for ($x = 0; $x < $nx; $x++)
+{
+ for ($y = 0; $y < $ny; $y++)
+ {
+   for ($z = 0; $z < $nz; $z++)
+   {
+     for ($i = 0; $i <= $#xunit; $i++)
+     {
+       $xnuc[$idx] = $x * $Lx + $xunit[$i] * $L + 0.2645886245 * $L;
+       $ynuc[$idx] = $y * $Ly + $yunit[$i] * $L + 0.2645886245 * $L;
+       $znuc[$idx] = $z * $Lz + $zunit[$i] * $L + 0.2645886245 * $L;
+       $idx++;
+     }
+   }
+ }
+} 
+
+$numnuc = $idx;
+
+# Print length of supercell
+
+printf("Created by AJB\n\n");
+printf("%d atoms\n",2*$numnuc);
+printf("2 atom types\n\n");
+printf("%f %f xlo xhi\n", 0, $Lx * $nx);
+printf("%f %f ylo yhi\n", 0, $Ly * $ny);
+printf("%f %f zlo zhi\n\n", 0, $Lz * $nz);
+printf("Masses\n\n");
+printf("1 6.941000\n");
+printf("2 1.000000\n\n");
+printf("Atoms\n\n");
+
+$j = 0;
+# Print out the nuclei
+for ($i = 0; $i < $numnuc; $i += 8)
+{
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1,1,3.0, 0, 0.0,$xnuc[$i], $ynuc[$i], $znuc[$i]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1,1,3.0, 0, 0.0,$xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1,1,3.0, 0, 0.0,$xnuc[$i+2], $ynuc[$i+2], $znuc[$i+2]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1,1,3.0, 0, 0.0,$xnuc[$i+3], $ynuc[$i+3], $znuc[$i+3]);
+}
+
+# Print out the core electrons
+for ($i = 0; $i < $numnuc; $i += 8)
+{
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_elec,$xnuc[$i  ], $ynuc[$i  ], $znuc[$i  ]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_elec,$xnuc[$i  ], $ynuc[$i  ], $znuc[$i  ]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_elec,$xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_elec,$xnuc[$i+1], $ynuc[$i+1], $znuc[$i+1]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_elec,$xnuc[$i+2], $ynuc[$i+2], $znuc[$i+2]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_elec,$xnuc[$i+2], $ynuc[$i+2], $znuc[$i+2]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, 1, $r_elec,$xnuc[$i+3], $ynuc[$i+3], $znuc[$i+3]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -1, $r_elec,$xnuc[$i+3], $ynuc[$i+3], $znuc[$i+3]);
+}
+
+# Print out the valence electrons
+for ($i = 0; $i < $numnuc; $i += 8)
+{
+ if (rand() < .5) {$spin = 1;} else {$spin = -1;}  
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, $spin, $r_elec2,$xnuc[$i+4], $ynuc[$i+4], $znuc[$i+4]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -$spin, $r_elec2,$xnuc[$i+5], $ynuc[$i+5], $znuc[$i+5]);
+ if (rand() < .5) {$spin = 1;} else {$spin = -1;}  
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, $spin, $r_elec2,$xnuc[$i+6], $ynuc[$i+6], $znuc[$i+6]);
+ printf("%i %i %f %i %f %f %f %f\n", $j+=1, 2, 0.0, -$spin, $r_elec2,$xnuc[$i+7], $ynuc[$i+7], $znuc[$i+7]);
+}
+