diff --git a/.gitignore b/.gitignore
index 1ce39c035..d4e6e5dcd 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,40 +1,42 @@
 *~
 /src/Makefile.package
 /src/Makefile.package.settings
 /src/lmp_mpi
 /src/style_*.h
 *.o
 *.so
 *.cu_o
 *.ptx
 *_ptx.h
 *.a
 *.d
 *.x
 *.exe
 *.dll
 *.pyc
 
 Obj_*
 log.lammps
 log.cite
 *.bz2
 *.gz
 *.tar
 .*.swp
 *.orig
 *.rej
 .vagrant
 .DS_Store
 .DS_Store?
 ._*
 .Spotlight-V100
 .Trashes
 ehthumbs.db
 Thumbs.db
 *.mod
 /src/pair_meam.cpp
 /src/pair_meam.h
 /lib/meam/Makefile.lammps
 /doc/html
 /src/TAGS
+/src/pair_neuronet.cpp
+/src/pair_neuronet.h
diff --git a/src/NEURONET/pair_neuronet.cpp b/src/NEURONET/pair_neuronet.cpp
index 5b22a04dc..8e478b5c9 100644
--- a/src/NEURONET/pair_neuronet.cpp
+++ b/src/NEURONET/pair_neuronet.cpp
@@ -1,440 +1,440 @@
 /* ----------------------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
 This file is written by Till Junge <till.junge@epfl.ch>
 ------------------------------------------------------------------------- */
 
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 #include "pair_neuronet.h"
 #include "atom.h"
 #include "comm.h"
 #include "force.h"
 #include "neighbor.h"
 #include "neigh_list.h"
 #include "neigh_request.h"
 #include "update.h"
 #include "integrate.h"
 //#include "respa.h"
 #include "math_const.h"
 #include "memory.h"
 #include "error.h"
 #include <fstream>
 
 
 extern "C" {
   void __nn_MOD_force_nn(int * namax, int * natm, int* tag, double* ra,
                  int * nnmax, double * force, double * vatom, double *h,
                  double *hi, double * tcom, int * nb, int * nbmax,
                  int * lsb, int * nex, int * lsrc, int * myparity, int * nn,
                  double * sv, double * rc, int *lspr, int * mpi_world,
                  int * myid, double * epi, double * epot, bool * vflag_atom);
   void __nn_MOD_read_params(int * myid, int * mpi_world, double * rcin,
                        double * rc3);
   double __nn_MOD_dsigmoid(double * x);
 }
 
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
 PairNeuroNet::PairNeuroNet(LAMMPS *lmp) : Pair(lmp)
 {
   respa_enable = 1;
   writedata = 1;
 
 
   this->ncnst_type[0]= 2;   // Gaussian
   this->ncnst_type[1]= 1;   // cosine
   this->ncnst_type[2]= 1;   // polynomial
   this->ncnst_type[3]= 2;   // Morse
   this->ncnst_type[100]= 1; // angular
 
   for (int i = 0; i < 100; i++) {
     this->ncomb_type[i]= 2;    // pair
     this->ncomb_type[100+i] = 3;  // triplet
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 PairNeuroNet::~PairNeuroNet()
 {
   if (allocated) {
-    delete[] this->itype;
-    delete[] this->cnst;
-    delete[] this->iaddr2;
-    delete[] this->iaddr3;
-    delete[] this->wgt11;
-    delete[] this->wgt12;
-    delete[] this->wgt21;
-    delete[] this->wgt22;
-    delete[] this->wgt23;
-
   }
 }
 
 /* ---------------------------------------------------------------------- */
 
 void PairNeuroNet::compute(int eflag, int vflag)
 {
-  int i,j,ii,jj,inum,jnum,itype,jtype;
+  int i,j,ii,jj,inum,jnum,jtype;
   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
   double rsq,r2inv,r6inv,forcelj,factor_lj;
   int *ilist,*jlist,*numneigh,**firstneigh;
 
   evdwl = 0.0;
   if (eflag || vflag) ev_setup(eflag,vflag);
   else evflag = vflag_fdotr = 0;
 
   double **x = atom->x;
   double **f = atom->f;
   int *type = atom->type;
   int nlocal = atom->nlocal;
   int nghost = atom->nghost;
   int nmax = atom->nmax;
   double *special_lj = force->special_lj;
   int newton_pair = force->newton_pair;
 
   inum = list->inum; //number of atoms that have a neighbor list
   ilist = list->ilist; //indices of those atoms
   numneigh = list->numneigh; // number of their neighbors
   firstneigh = list->firstneigh; // pointer to first neighbor index (indices contiguous
 
   // loop over neighbors of my atoms
 
   int nnmax = 0;
   for (int ii = 0; ii < inum; ++ii) {
     i = ilist[ii];
     nnmax = (nnmax > numneigh[i]) ? nnmax : numneigh[i];
   }
   // ryo neighlist is matrix
   int fort_neighbor[(nnmax+1) * nmax] = {0};
   int fort_type[nmax] = {0};
   double fort_x [3*nmax] = {0};
   double fort_f [3*nmax] = {0};
   double fort_vatom [6*nmax] = {0};
   double fort_eatom [nmax] = {0};
 
   for (int ii = 0; ii < inum; ++ii) {
     i = ilist[ii];
     for (int dir = 0; dir < 3; ++dir) {
       fort_x[dir+3*ii] = x[i][dir];
     }
     int *neigh_ptr = firstneigh[i];
     int nb_neigh = numneigh[i];
     fort_neighbor[0 + (nnmax+1)*ii] = nb_neigh;
     for (int j = 0; j < nb_neigh; j++) {
       fort_neighbor[j+1 + (nnmax+1)*ii] = neigh_ptr[j];
     }
     fort_type[ii] = type[i];
   }
 
   double h[] = {1,0,0,
                 0,1,0,
                 0,0,1};
   double dummy = 0;
   int idummy = 0;
   double delta_eng_vdwl;
   bool fort_vflag_atom = bool(this->vflag_atom);
 
   __nn_MOD_force_nn(&nmax, &inum, fort_type, fort_x, &nnmax, fort_f,
             fort_vatom, h, h, &dummy, &nghost, &idummy, &idummy, &idummy,
             &idummy, &idummy, &idummy, &dummy, &rcin, fort_neighbor, &idummy,
             &idummy, fort_eatom, &delta_eng_vdwl, &fort_vflag_atom);
 
   this->eng_vdwl += delta_eng_vdwl;
   // rewrite results into lammps form
 
   for (int ii = 0; ii < inum; ++ii) {
     i = ilist[ii];
     eatom[i] += fort_eatom[ii];
     for (int dir = 0; dir < 3; ++dir) {
       f[i][dir] += fort_f[dir+3*ii];
     }
     for (int dir = 0; dir < 6; ++dir) {
       vatom[i][dir] += fort_vatom[dir+6*ii];
     }
   }
 
   if (vflag_fdotr) virial_fdotr_compute();
 }
 
 /* ---------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
    global settings
 ------------------------------------------------------------------------- */
 
 void PairNeuroNet::settings(int narg, char **arg) // read_params
 {
   if (narg != 0) error->all(FLERR,"Illegal pair_style command");
 }
 
 /* ----------------------------------------------------------------------
    set coeffs for one or more type pairs
 ------------------------------------------------------------------------- */
 
 void PairNeuroNet::coeff(int narg, char **arg)
 {
   if (narg < 5)
     error->all(FLERR,"Not ennough arguments, should be * * constants_input_file params_input_file ");
 
   // insure I,J args are * *
 
   if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
     error->all(FLERR,"Incorrect args for pair coefficients");
 
   // Read constants input file
   std::ifstream const_file(arg[2]);
   //TODO: check that file exist etc
   int dummy;
   const_file >> this->nl >> this->nsp >> dummy;
   // dummy contains number of input nodes
   this->nhl.push_back(dummy);
   for (int i = 0; i < nl; ++i) {
     // read the number of nodes per hidden layer
     const_file >> dummy;
     this->nhl.push_back(dummy);
   }
 
   // last layer always has one node
   nhl.push_back(1);
 
   // allocate arrays for consts
-  this->itype = new int[this->nhl[0]];
-  this->cnst = new double[this->max_ncnst*this->nhl[0]];
-  this->iaddr2 = new int[2*this->nsp*this->nsp];
-  this->iaddr3 = new int[2*this->nsp*this->nsp*this->nsp];
+  this->itype.resize(this->nhl[0]);
+  this->cnst.resize(this->max_ncnst*this->nhl[0]);
+  this->iaddr2.resize(2*this->nsp*this->nsp);
+  this->iaddr3.resize(2*this->nsp*this->nsp*this->nsp);
   for (int i = 0; i < 2*this->nsp*this->nsp; i++) {
     this-> iaddr2[i] = -1;
     for (int j = 0; j < this->nsp; j++) {
       this->iaddr3[j+i*nsp] = -1;
     }
   }
 
   //read in consts
 
   int icmb[3] = {0}; // three is hardcoded maximum n in n-body potential
 
   int iap = 0;
   int jap = 0;
   int kap = 0;
   int nsf1 = 0, nsf2 = 0;
    
   for (int i = 0; i < this->nhl[0]; i++) {
     const_file >> this->itype[i];
     for (int j = 0; i < this->ncomb_type[itype[i]]; i++) {
       const_file >> icmb[j];
     }
 
     for (int j = 0; i < this->ncnst_type[itype[i]]; i++) {
       const_file >> this->cnst[j+i*nhl[0]];
     }
 
     // distinguish n-body potentials
     
     if (itype[i] < 100) {
       if ((icmb[0] != iap) || (icmb[1] != jap)) {
         this->iaddr2[0 + this->nsp*(icmb[0] + nsp*icmb[1])] = i+1;
         this->iaddr2[0 + this->nsp*(icmb[1] + nsp*icmb[0])] = i+1;
       }
       this->iaddr2[1 + this->nsp*(icmb[0] + nsp*icmb[1])] = i+1;
       this->iaddr2[1 + this->nsp*(icmb[1] + nsp*icmb[0])] = i+1;
       nsf1++;
       iap = icmb[0];
       jap = icmb[1];
     }
     else if (itype[i] < 200){
       if ((icmb[0] != iap) || (icmb[1] != jap) || (icmb[2] != kap)) {
         this->iaddr3[0 + this->nsp*(icmb[0] + nsp*(icmb[1] + nsp*icmb[2]))] = i+1;
         this->iaddr3[0 + this->nsp*(icmb[0] + nsp*(icmb[2] + nsp*icmb[1]))] = i+1;
       }
       this->iaddr3[1 + this->nsp*(icmb[0] + nsp*(icmb[1] + nsp*icmb[2]))] = i+1;
       this->iaddr3[1 + this->nsp*(icmb[0] + nsp*(icmb[2] + nsp*icmb[1]))] = i+1;
       nsf2++;
       iap = icmb[0];
       jap = icmb[1];
       kap = icmb[2];
     }
     else {
       error->all(FLERR,"Unknown itype");
     }
   }
 
   if (this-> nhl[0] != nsf1 + nsf2) {
     error->all(FLERR, "[Error] nsf.ne.nsf1+nsf2 !!!");
   }
 
   // read in weights
   double nwgt[nl+1];
   for (int i = 0 ; i < nl+1; i++) {
     nwgt[i] = nhl[i] * nhl[i+1];
   }
 
   std::ifstream weights_file(arg[3]);
   // TODO check file existence etc
   int ncoeff;
   weights_file >> ncoeff >> this->rcin >> this->rc3;
   if (rc3 > rcin) {
     rc3 = rcin;
     // TODO add warning as in ryo_force_NN.F90:680
   }
   int nc = 0;
   for (int i = 0; i < nl+1; i++) {
     nc += nwgt[i];
   }
   if( ncoeff != nc ) {
     error->all(FLERR, "[Error] num of parameters is not correct !!!");
   }
 
   if (nl == 1) {
-    this->wgt11 = new double[nhl[0] * nhl[1]];
-    this->wgt12 = new double[nhl[1]];
-    this->wgt21 = new double[1];
-    this->wgt22 = new double[1];
-    this->wgt23 = new double[1];
+    this->wgt11.resize(nhl[0] * nhl[1]);
+    this->wgt12.resize(nhl[1]);
 
     for (int ihl0 = 0; ihl0 < nhl[0]; ihl0++) {
       for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
         weights_file >> this->wgt11[ihl0 + ihl1*nhl[0]];
       }
     }
     for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
       weights_file >> this->wgt12[ihl1];
     }
 
   } else if (nl == 2) {
-    this->wgt11 = new double[1];
-    this->wgt12 = new double[1];
-    this->wgt21 = new double[nhl[0] * nhl[1]];
-    this->wgt22 = new double[nhl[1] * nhl[2]];
-    this->wgt23 = new double[nhl[2]];
+    this->wgt21.resize(nhl[0] * nhl[1]);
+    this->wgt22.resize(nhl[1] * nhl[2]);
+    this->wgt23.resize(nhl[2]);
 
     for (int ihl0 = 0; ihl0 < nhl[0]; ihl0++) {
       for (int ihl1 = 0; ihl1 < nhl[1]; ihl1++) {
         weights_file >> this->wgt21[ihl0 + ihl1*nhl[0]];
       }
     }
     for (int ihl0 = 0; ihl0 < nhl[1]; ihl0++) {
       for (int ihl1 = 0; ihl1 < nhl[2]; ihl1++) {
         weights_file >> this->wgt22[ihl0 + ihl1*nhl[0]];
       }
     }
     for (int ihl1 = 0; ihl1 < nhl[2]; ihl1++) {
       weights_file >> this->wgt23[ihl1];
     }
   }
 
   // allocate arrays for weights
   this->allocated = true;
 }
 
 /* ----------------------------------------------------------------------
    init specific to this pair style
 ------------------------------------------------------------------------- */
 
 void PairNeuroNet::init_style()
 {
   // Don't understand this yet, but seems to be what lammps wants here
   int irequest = neighbor->request(this,instance_me);
 }
 
 /* ----------------------------------------------------------------------
    init for one type pair i,j and corresponding j,i
 ------------------------------------------------------------------------- */
 
 double PairNeuroNet::init_one(int i, int j)
 {
   return this->rcin;
 }
 
 /* ----------------------------------------------------------------------
    proc 0 writes to restart file
 ------------------------------------------------------------------------- */
+template<typename T>
+void vector_write(const std::vector<T> & vec, FILE *fp) {
+  size_t siz = vec.size();
+  fwrite(&siz, sizeof(size_t), 1, fp);
+  fwrite(&vec[0], sizeof(T), siz, fp);
+}
+template<typename T>
+void vector_read(std::vector<T> & vec, FILE *fp) {
+  size_t siz;
+  fread(&siz, sizeof(size_t), 1, fp);
+  vec.resize(siz);
+  fwrite(&vec[0], sizeof(T), siz, fp);
+}
+template<typename T>
+void scalar_write(const T & scalar, FILE *fp) {
+  fwrite(&scalar, sizeof(T), 1, fp);
+}
+template<typename T>
+void scalar_read(T & scalar, FILE *fp) {
+  fread(&scalar, sizeof(T), 1, fp);
+}
 
 void PairNeuroNet::write_restart(FILE *fp)
 {
   write_restart_settings(fp);
-
-  //int i,j;
-  //for (i = 1; i <= atom->ntypes; i++)
-  //  for (j = i; j <= atom->ntypes; j++) {
-  //    fwrite(&setflag[i][j],sizeof(int),1,fp);
-  //    if (setflag[i][j]) {
-  //      fwrite(&epsilon[i][j],sizeof(double),1,fp);
-  //      fwrite(&sigma[i][j],sizeof(double),1,fp);
-  //      fwrite(&cut[i][j],sizeof(double),1,fp);
-  //    }
-  //  }
+  vector_write(this->wgt11, fp);
+  vector_write(this->wgt12, fp);
+  vector_write(this->wgt21, fp);
+  vector_write(this->wgt22, fp);
+  vector_write(this->wgt23, fp);
+  vector_write(this->nhl, fp);
+  vector_write(this->itype, fp);
+  vector_write(this->iaddr2, fp);
+  vector_write(this->iaddr3, fp);
+  vector_write(this->cnst, fp);
+  scalar_write(this->rcin, fp);
+  scalar_write(this->rc3, fp);
 }
 
 /* ----------------------------------------------------------------------
    proc 0 reads from restart file, bcasts
 ------------------------------------------------------------------------- */
 
 void PairNeuroNet::read_restart(FILE *fp)
 {
   read_restart_settings(fp);
-  //allocate();
-  //
-  //int i,j;
-  //int me = comm->me;
-  //for (i = 1; i <= atom->ntypes; i++)
-  //  for (j = i; j <= atom->ntypes; j++) {
-  //    if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp);
-  //    MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world);
-  //    if (setflag[i][j]) {
-  //      if (me == 0) {
-  //        fread(&epsilon[i][j],sizeof(double),1,fp);
-  //        fread(&sigma[i][j],sizeof(double),1,fp);
-  //        fread(&cut[i][j],sizeof(double),1,fp);
-  //      }
-  //      MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world);
-  //      MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world);
-  //      MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world);
-  //    }
-  //  }
+  vector_read(this->wgt11, fp);
+  vector_read(this->wgt12, fp);
+  vector_read(this->wgt21, fp);
+  vector_read(this->wgt22, fp);
+  vector_read(this->wgt23, fp);
+  vector_read(this->nhl, fp);
+  vector_read(this->itype, fp);
+  vector_read(this->iaddr2, fp);
+  vector_read(this->iaddr3, fp);
+  vector_read(this->cnst, fp);
+  scalar_read(this->rcin, fp);
+  scalar_read(this->rc3, fp);
 }
 
 /* ----------------------------------------------------------------------
    proc 0 writes to restart file
 ------------------------------------------------------------------------- */
 
 void PairNeuroNet::write_restart_settings(FILE *fp)
 {
   // think this is supposed to be empty for us
 }
 
 /* ----------------------------------------------------------------------
    proc 0 reads from restart file, bcasts
 ------------------------------------------------------------------------- */
 
 void PairNeuroNet::read_restart_settings(FILE *fp)
 {
   int me = comm->me;
   if (me == 0) {
     // think this is supposed to be empty for us
   }
   //MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world);
   //MPI_Bcast(&offset_flag,1,MPI_INT,0,world);
   //MPI_Bcast(&mix_flag,1,MPI_INT,0,world);
   //MPI_Bcast(&tail_flag,1,MPI_INT,0,world);
 }
 
 const int PairNeuroNet::nlmax = 2;
 const int PairNeuroNet::max_ncnst = 2;
diff --git a/src/NEURONET/pair_neuronet.h b/src/NEURONET/pair_neuronet.h
index e05901cb0..64f3da2ca 100644
--- a/src/NEURONET/pair_neuronet.h
+++ b/src/NEURONET/pair_neuronet.h
@@ -1,87 +1,87 @@
 /* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov
 
    Copyright (2003) Sandia Corporation.  Under the terms of Contract
    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
    certain rights in this software.  This software is distributed under
    the GNU General Public License.
 
    See the README file in the top-level LAMMPS directory.
 
 This file is written by Till Junge <till.junge@epfl.ch>
 ------------------------------------------------------------------------- */
 
 #ifdef PAIR_CLASS
 
 PairStyle(neuronet,PairNeuroNet)
 
 #else
 
 #ifndef LMP_PAIR_NEURO_NET_H
 #define LMP_PAIR_NEURO_NET_H
 
 #include "pair.h"
 #include <vector>
 
 
 namespace LAMMPS_NS {
 
 class PairNeuroNet : public Pair {
  public:
   PairNeuroNet(class LAMMPS *);
   virtual ~PairNeuroNet();
   virtual void compute(int, int);
   void settings(int, char **);
   void coeff(int, char **);
   void init_style();
   double init_one(int, int);
   void write_restart(FILE *);
   void read_restart(FILE *);
   void write_restart_settings(FILE *);
   void read_restart_settings(FILE *);
 
  protected:
 
   // weight parameters (indices: nb of hidden layers, layer)
   // weights for one hidden layer
-  double * wgt11; // row is hidden id, column is input id
-  double * wgt12;
+  std::vector<double> wgt11; // row is hidden id, column is input id
+  std::vector<double> wgt12;
   // weights for two hidden layers
-  double * wgt21;
-  double * wgt22;
-  double * wgt23;
+  std::vector<double> wgt21;
+  std::vector<double> wgt22;
+  std::vector<double> wgt23;
 
   // size values
   int nl; //number of layers
   int nsp; // number of species
   std::vector<int>  nhl; // number of nodes per hidden layer
 
-  int * itype; // array of symmetry function ids
+  std::vector<int> itype; // array of symmetry function ids
   // contains beginnings and ends of symmetry function id to consider for each species pair: iaddr2(i, j, k) i: 1 for start, 2 for end; j, k pair of species 
-  int * iaddr2; // for 2-body potential
-  int * iaddr3; // for 3-body potential
+  std::vector<int> iaddr2; // for 2-body potential
+  std::vector<int> iaddr3; // for 3-body potential
   // parameters of the symmetry functions
   // eta   = cnst(1, isf) // isf is symmetry function id
   // alpha = cnst(2, isf)
-  double * cnst;
+  std::vector<double> cnst;
 
   // cutoff radii for 2 and 3 body potentials
   double rcin, rc3;
 
   // nb of parameters per symmetry function type
   int ncnst_type[200];
   // n in n-body potential per symmetry function type
   int ncomb_type[200];
 
   static const int nlmax;
   static const int max_ncnst;
 
 };
 
 }
 
 #endif
 #endif