diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp
index d4e57c2c6cb62eda855b0cc0e9556bae83f3e3bb..2bb1c1da22f1baf9f8bf65034149e6cd984db8a2 100644
--- a/src/MANYBODY/pair_airebo.cpp
+++ b/src/MANYBODY/pair_airebo.cpp
@@ -592,6 +592,8 @@ void PairAIREBO::FLJ(int eflag, int vflag)
       // if outside of 2-path cutoff but inside 4-path cutoff,
       //   best = 0.0, test 3-,4-paths
       // if inside 2-path cutoff, best = wij, only test 3-,4-paths if best < 1
+      npath = testpath = done = 0;
+      best = 0.0;
 
       if (rijsq >= cutljsq[itype][jtype]) continue;
       rij = sqrt(rijsq);
@@ -608,7 +610,6 @@ void PairAIREBO::FLJ(int eflag, int vflag)
         else testpath = 0;
       }
 
-      done = 0;
       if (testpath) {
 
         // test all 3-body paths = I-K-J
@@ -629,7 +630,7 @@ void PairAIREBO::FLJ(int eflag, int vflag)
           if (rsq < rcmaxsq[itype][ktype]) {
             rik = sqrt(rsq);
             wik = Sp(rik,rcmin[itype][ktype],rcmax[itype][ktype],dwik);
-          } else wik = 0.0;
+          } else { dwik = wik = 0.0; rikS = rik = 1.0; }
 
           if (wik > best) {
             deljk[0] = x[j][0] - x[k][0];
@@ -679,7 +680,7 @@ void PairAIREBO::FLJ(int eflag, int vflag)
               if (rsq < rcmaxsq[ktype][mtype]) {
                 rkm = sqrt(rsq);
                 wkm = Sp(rkm,rcmin[ktype][mtype],rcmax[ktype][mtype],dwkm);
-              } else wkm = 0.0;
+              } else { dwkm = wkm = 0.0; rkmS = rkm = 1.0; }
 
               if (wik*wkm > best) {
                 deljm[0] = x[j][0] - x[m][0];
diff --git a/src/MOLECULE/improper_umbrella.cpp b/src/MOLECULE/improper_umbrella.cpp
index 14551e6c5c4f191b6b7964602f5f6e57e27ffe5e..30a62ffb1828343f484ae42c440b92196ca82f5c 100644
--- a/src/MOLECULE/improper_umbrella.cpp
+++ b/src/MOLECULE/improper_umbrella.cpp
@@ -37,7 +37,10 @@ using namespace MathConst;
 
 /* ---------------------------------------------------------------------- */
 
-ImproperUmbrella::ImproperUmbrella(LAMMPS *lmp) : Improper(lmp) {}
+ImproperUmbrella::ImproperUmbrella(LAMMPS *lmp) : Improper(lmp)
+{
+  writedata = 1;
+}
 
 /* ---------------------------------------------------------------------- */
 
diff --git a/src/Makefile b/src/Makefile
index 85b212d539003945110b5d2554a2d8b231362c85..7cf9d1730fe9fb0722419076c1ab421d8b7a2c9e 100755
--- a/src/Makefile
+++ b/src/Makefile
@@ -54,7 +54,7 @@ PACKUSER = user-atc user-awpmd user-cg-cmm user-colvars \
 	   user-vtk
 
 PACKLIB = compress gpu kim kokkos meam poems python reax voronoi \
-	  user-atc user-awpmd user-colvars user-cuda user-h5md user-molfile \
+	  user-atc user-awpmd user-colvars user-h5md user-molfile \
 	  user-qmmm user-quip user-vtk
 
 PACKALL = $(PACKAGE) $(PACKUSER)
diff --git a/src/Purge.list b/src/Purge.list
index 79ad3c0c42cc2c9454f648de0a7c1f3ae6ace926..37d137e838782bf9233f88bbf2aae01d75cfab5d 100644
--- a/src/Purge.list
+++ b/src/Purge.list
@@ -13,6 +13,140 @@ style_kspace.h
 style_minimize.h
 style_pair.h
 style_region.h
+# deleted on 31 May 2016
+fix_ave_spatial_sphere.cpp
+fix_ave_spatial_sphere.h
+atom_vec_angle_cuda.cpp
+atom_vec_angle_cuda.h
+atom_vec_atomic_cuda.cpp
+atom_vec_atomic_cuda.h
+atom_vec_charge_cuda.cpp
+atom_vec_charge_cuda.h
+atom_vec_full_cuda.cpp
+atom_vec_full_cuda.h
+comm_cuda.cpp
+comm_cuda.h
+compute_pe_cuda.cpp
+compute_pe_cuda.h
+compute_pressure_cuda.cpp
+compute_pressure_cuda.h
+compute_temp_cuda.cpp
+compute_temp_cuda.h
+compute_temp_partial_cuda.cpp
+compute_temp_partial_cuda.h
+cuda.cpp
+cuda_data.h
+cuda_modify_flags.h
+cuda_neigh_list.cpp
+cuda_neigh_list.h
+domain_cuda.cpp
+domain_cuda.h
+fft3d_cuda.cpp
+fft3d_cuda.h
+fft3d_wrap_cuda.cpp
+fft3d_wrap_cuda.h
+fix_addforce_cuda.cpp
+fix_addforce_cuda.h
+fix_aveforce_cuda.cpp
+fix_aveforce_cuda.h
+fix_enforce2d_cuda.cpp
+fix_enforce2d_cuda.h
+fix_freeze_cuda.cpp
+fix_freeze_cuda.h
+fix_gravity_cuda.cpp
+fix_gravity_cuda.h
+fix_nh_cuda.cpp
+fix_nh_cuda.h
+fix_npt_cuda.cpp
+fix_npt_cuda.h
+fix_nve_cuda.cpp
+fix_nve_cuda.h
+fix_nvt_cuda.cpp
+fix_nvt_cuda.h
+fix_set_force_cuda.cpp
+fix_set_force_cuda.h
+fix_shake_cuda.cpp
+fix_shake_cuda.h
+fix_temp_berendsen_cuda.cpp
+fix_temp_berendsen_cuda.h
+fix_temp_rescale_cuda.cpp
+fix_temp_rescale_cuda.h
+fix_temp_rescale_limit_cuda.cpp
+fix_temp_rescale_limit_cuda.h
+fix_viscous_cuda.cpp
+fix_viscous_cuda.h
+modify_cuda.cpp
+modify_cuda.h
+neighbor_cuda.cpp
+neighbor_cuda.h
+neigh_full_cuda.cpp
+pair_born_coul_long_cuda.cpp
+pair_born_coul_long_cuda.h
+pair_buck_coul_cut_cuda.cpp
+pair_buck_coul_cut_cuda.h
+pair_buck_coul_long_cuda.cpp
+pair_buck_coul_long_cuda.h
+pair_buck_cuda.cpp
+pair_buck_cuda.h
+pair_eam_alloy_cuda.cpp
+pair_eam_alloy_cuda.h
+pair_eam_cuda.cpp
+pair_eam_cuda.h
+pair_eam_fs_cuda.cpp
+pair_eam_fs_cuda.h
+pair_gran_hooke_cuda.cpp
+pair_gran_hooke_cuda.h
+pair_lj96_cut_cuda.cpp
+pair_lj96_cut_cuda.h
+pair_lj_charmm_coul_charmm_cuda.cpp
+pair_lj_charmm_coul_charmm_cuda.h
+pair_lj_charmm_coul_charmm_implicit_cuda.cpp
+pair_lj_charmm_coul_charmm_implicit_cuda.h
+pair_lj_charmm_coul_long_cuda.cpp
+pair_lj_charmm_coul_long_cuda.h
+pair_lj_class2_coul_cut_cuda.cpp
+pair_lj_class2_coul_cut_cuda.h
+pair_lj_class2_coul_long_cuda.cpp
+pair_lj_class2_coul_long_cuda.h
+pair_lj_class2_cuda.cpp
+pair_lj_class2_cuda.h
+pair_lj_cut_coul_cut_cuda.cpp
+pair_lj_cut_coul_cut_cuda.h
+pair_lj_cut_coul_debye_cuda.cpp
+pair_lj_cut_coul_debye_cuda.h
+pair_lj_cut_coul_long_cuda.cpp
+pair_lj_cut_coul_long_cuda.h
+pair_lj_cut_cuda.cpp
+pair_lj_cut_cuda.h
+pair_lj_cut_experimental_cuda.cpp
+pair_lj_cut_experimental_cuda.h
+pair_lj_expand_cuda.cpp
+pair_lj_expand_cuda.h
+pair_lj_gromacs_coul_gromacs_cuda.cpp
+pair_lj_gromacs_coul_gromacs_cuda.h
+pair_lj_gromacs_cuda.cpp
+pair_lj_gromacs_cuda.h
+pair_lj_sdk_coul_long_cuda.cpp
+pair_lj_sdk_coul_long_cuda.h
+pair_lj_sdk_cuda.cpp
+pair_lj_sdk_cuda.h
+pair_lj_smooth_cuda.cpp
+pair_lj_smooth_cuda.h
+pair_morse_cuda.cpp
+pair_morse_cuda.h
+pair_sw_cuda.cpp
+pair_sw_cuda.h
+pair_tersoff_cuda.cpp
+pair_tersoff_cuda.h
+pair_tersoff_zbl_cuda.cpp
+pair_tersoff_zbl_cuda.h
+pppm_cuda.cpp
+pppm_cuda.h
+pppm_old.cpp
+pppm_old.h
+user_cuda.h
+verlet_cuda.cpp
+verlet_cuda.h
 # deleted on 11 May 2016
 pair_dpd_conservative.cpp
 pair_dpd_conservative.h
diff --git a/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp b/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp
index d2d73de6bdf4fa346804a5a20a4ce48271c2e21d..f485c4b7b4b41b8b19a1f1fb1984b190c45a56aa 100644
--- a/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp
+++ b/src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp
@@ -64,14 +64,16 @@ enum { CONST, EQUAL }; // For treating the variables.
 
 static const char* cite_fix_nve_manifold_rattle =
   "fix nve/manifold/rattle command:\n\n"
-  "@article{paquay-2016,\n"
-  " author = {Paquay, Stefan and Kusters, Remy}, \n"
-  " eprint = {arXiv:1411.3019}, \n"
-  " title = {A method for molecular dynamics on curved surfaces}, \n"
-  " url = {http://arxiv.org/abs/1411.3019}, \n"
-  " year = 2016, \n"
-  " journal = {Arxiv preprint, \\url{http://arxiv.org/abs/1411.3019}}, \n"
-  " note = {To be published in Biophys. J.}, \n"
+  "   author        = {Paquay, Stefan and Kusters, Remy},\n"
+  "   doi           = {10.1016/j.bpj.2016.02.017},\n"
+  "   issn          = {0006-3495},\n"
+  "   journal       = {Biophysical Journal},\n"
+  "   month         = apr,\n"
+  "   number        = {6},\n"
+  "   pages         = {1226--1233},\n"
+  "   title         = {{A Method for Molecular Dynamics on Curved Surfaces}},\n"
+  "   volume        = {110},\n"
+  "   year          = {2016}\n"
   "}\n\n";
 
 
diff --git a/src/USER-MANIFOLD/fix_nvt_manifold_rattle.cpp b/src/USER-MANIFOLD/fix_nvt_manifold_rattle.cpp
index 3801508336299ef6bd7bcebcc81ab920b7de6f9d..c304a8718d54360e2a44bd6ba7c2a3a9f0adc302 100644
--- a/src/USER-MANIFOLD/fix_nvt_manifold_rattle.cpp
+++ b/src/USER-MANIFOLD/fix_nvt_manifold_rattle.cpp
@@ -65,18 +65,18 @@ enum {NOBIAS,BIAS};
 
 
 static const char* cite_fix_nvt_manifold_rattle =
-  "fix nve/manifold/rattle command:\n\n"
-  "@article{paquay-2016,\n"
-  " author = {Paquay, Stefan and Kusters, Remy}, \n"
-  " eprint = {arXiv:1411.3019}, \n"
-  " title = {A method for molecular dynamics on curved surfaces}, \n"
-  " url = {http://arxiv.org/abs/1411.3019}, \n"
-  " year = 2016, \n"
-  " journal = {Arxiv preprint, \\url{http://arxiv.org/abs/1411.3019}}, \n"
-  " note = {To be published in Biophys. J.}, \n"
+  "fix nvt/manifold/rattle command:\n\n"
+  "   author        = {Paquay, Stefan and Kusters, Remy},\n"
+  "   doi           = {10.1016/j.bpj.2016.02.017},\n"
+  "   issn          = {0006-3495},\n"
+  "   journal       = {Biophysical Journal},\n"
+  "   month         = apr,\n"
+  "   number        = {6},\n"
+  "   pages         = {1226--1233},\n"
+  "   title         = {{A Method for Molecular Dynamics on Curved Surfaces}},\n"
+  "   volume        = {110},\n"
+  "   year          = {2016}\n"
   "}\n\n";
-
-
 /* ---------------------------------------------------------------------- */
 
 FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
diff --git a/src/USER-OMP/pair_airebo_omp.cpp b/src/USER-OMP/pair_airebo_omp.cpp
index 35a644007ac9060a93061d9438baf4ba1ffe434f..27089fd6b3ef4ba673b48be9bee265945652dde8 100644
--- a/src/USER-OMP/pair_airebo_omp.cpp
+++ b/src/USER-OMP/pair_airebo_omp.cpp
@@ -2022,6 +2022,8 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
       // if outside of 2-path cutoff but inside 4-path cutoff,
       //   best = 0.0, test 3-,4-paths
       // if inside 2-path cutoff, best = wij, only test 3-,4-paths if best < 1
+      npath = testpath = done = 0;
+      best = 0.0;
 
       if (rijsq >= cutljsq[itype][jtype]) continue;
       rij = sqrt(rijsq);
@@ -2038,7 +2040,6 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
         else testpath = 0;
       }
 
-      done = 0;
       if (testpath) {
 
         // test all 3-body paths = I-K-J
@@ -2059,7 +2060,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
           if (rsq < rcmaxsq[itype][ktype]) {
             rik = sqrt(rsq);
             wik = Sp(rik,rcmin[itype][ktype],rcmax[itype][ktype],dwik);
-          } else wik = 0.0;
+          } else { dwik = wik = 0.0; rikS = rik = 1.0; }
 
           if (wik > best) {
             deljk[0] = x[j][0] - x[k][0];
@@ -2109,7 +2110,7 @@ void PairAIREBOOMP::FLJ_thr(int ifrom, int ito, int evflag, int eflag,
               if (rsq < rcmaxsq[ktype][mtype]) {
                 rkm = sqrt(rsq);
                 wkm = Sp(rkm,rcmin[ktype][mtype],rcmax[ktype][mtype],dwkm);
-              } else wkm = 0.0;
+              } else { dwkm = wkm = 0.0; rkmS = rkm = 1.0; }
 
               if (wik*wkm > best) {
                 deljm[0] = x[j][0] - x[m][0];
diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp
index 352a9468b33358b37a24e8384dd21afd7e5eabe8..59b39d45345cf320c27298a81a7f394b1ae40d03 100644
--- a/src/USER-REAXC/pair_reax_c.cpp
+++ b/src/USER-REAXC/pair_reax_c.cpp
@@ -247,7 +247,7 @@ void PairReaxC::settings(int narg, char **arg)
       system->safezone = force->numeric(FLERR,arg[iarg+1]);
       if (system->safezone < 0.0)
 	error->all(FLERR,"Illegal pair_style reax/c safezone command");
-      system->saferzone = system->safezone + 0.2;
+      system->saferzone = system->safezone*1.2;
       iarg += 2;
     } else if (strcmp(arg[iarg],"mincap") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal pair_style reax/c command");
diff --git a/src/VORONOI/compute_voronoi_atom.cpp b/src/VORONOI/compute_voronoi_atom.cpp
index f4f48fd5e8f40f35f9a80c8ff9ee73c8c8a323d0..5239a27ecc6b2f4e29601a309c683443db98b237 100644
--- a/src/VORONOI/compute_voronoi_atom.cpp
+++ b/src/VORONOI/compute_voronoi_atom.cpp
@@ -61,6 +61,7 @@ ComputeVoronoi::ComputeVoronoi(LAMMPS *lmp, int narg, char **arg) :
   con_mono = NULL;
   con_poly = NULL;
   tags = NULL;
+  oldmaxtag = 0;
   occvec = sendocc = lroot = lnext = NULL;
   faces = NULL;
 
@@ -177,6 +178,8 @@ ComputeVoronoi::~ComputeVoronoi()
 
 void ComputeVoronoi::init()
 {
+  if (occupation && (atom->tag_enable == 0))
+    error->all(FLERR,"Compute voronoi/atom occupation requires atom IDs");
 }
 
 /* ----------------------------------------------------------------------
@@ -213,11 +216,14 @@ void ComputeVoronoi::compute_peratom()
       lnext = NULL;
       lmax = 0;
 
-      // build the occupation buffer
+      // build the occupation buffer.
+      // NOTE: we cannot make it of length oldnatoms, as tags may not be
+      // consecutive at this point, e.g. due to deleted or lost atoms.
       oldnatoms = atom->natoms;
-      memory->create(occvec,oldnatoms,"voronoi/atom:occvec");
+      oldmaxtag = atom->map_tag_max;
+      memory->create(occvec,oldmaxtag,"voronoi/atom:occvec");
 #ifdef NOTINPLACE
-      memory->create(sendocc,oldnatoms,"voronoi/atom:sendocc");
+      memory->create(sendocc,oldmaxtag,"voronoi/atom:sendocc");
 #endif
     }
 
@@ -440,7 +446,14 @@ void ComputeVoronoi::checkOccupation()
   // cherry pick currently owned atoms
   for (i=0; i<nlocal; i++) {
     // set the new atom count in the atom's first frame voronoi cell
-    voro[i][0] = occvec[atom->tag[i]-1];
+    // but take into account that new atoms might have been added to
+    // the system, so we can only look up occupancy for tags that are
+    // smaller or equal to the recorded largest tag.
+    tagint mytag = atom->tag[i];
+    if (mytag > oldmaxtag)
+      voro[i][0] = 0;
+    else
+      voro[i][0] = occvec[mytag-1];
   }
 }
 
diff --git a/src/VORONOI/compute_voronoi_atom.h b/src/VORONOI/compute_voronoi_atom.h
index 7d23a611bdcf83bc2e82cc4247742a497400533e..8fd187cea348e55d508e6ab61ad46868c364a920 100644
--- a/src/VORONOI/compute_voronoi_atom.h
+++ b/src/VORONOI/compute_voronoi_atom.h
@@ -55,7 +55,7 @@ class ComputeVoronoi : public Compute {
   enum { VOROSURF_NONE, VOROSURF_ALL, VOROSURF_GROUP } surface;
   bool onlyGroup, occupation;
 
-  tagint *tags;
+  tagint *tags, oldmaxtag;
   int *occvec, *sendocc, *lroot, *lnext, lmax, oldnatoms, oldnall;
   int faces_flag, nfaces, nfacesmax;
   double **faces;
diff --git a/src/compute_com.cpp b/src/compute_com.cpp
index e086210322481cacf70c0fb240b47748b50e63e1..27606cbc6f147a8af6b311d17e1605bfe02d2f31 100644
--- a/src/compute_com.cpp
+++ b/src/compute_com.cpp
@@ -51,6 +51,7 @@ void ComputeCOM::init()
 void ComputeCOM::compute_vector()
 {
   invoked_vector = update->ntimestep;
+  if (group->dynamic[igroup]) masstotal = group->mass(igroup);
 
   group->xcm(igroup,masstotal,vector);
 }
diff --git a/src/compute_gyration.cpp b/src/compute_gyration.cpp
index fab413e670a3409bd6ca8b59497f5c1236d927c2..0c1c58a4fd5134ae831107b4118f9d47e998663d 100644
--- a/src/compute_gyration.cpp
+++ b/src/compute_gyration.cpp
@@ -57,6 +57,7 @@ double ComputeGyration::compute_scalar()
   invoked_scalar = update->ntimestep;
 
   double xcm[3];
+  if (group->dynamic[igroup]) masstotal = group->mass(igroup);
   group->xcm(igroup,masstotal,xcm);
   scalar = group->gyration(igroup,masstotal,xcm);
   return scalar;
diff --git a/src/compute_msd.cpp b/src/compute_msd.cpp
index 32d93f933426d0cb063261eac2ea0adf59f3d661..7f005d27593a7e75d98e6c169442ff99c2c5ddd6 100644
--- a/src/compute_msd.cpp
+++ b/src/compute_msd.cpp
@@ -34,6 +34,7 @@ ComputeMSD::ComputeMSD(LAMMPS *lmp, int narg, char **arg) :
   size_vector = 4;
   extvector = 0;
   create_attribute = 1;
+  dynamic_group_allow = 0;
 
   // optional args
 
diff --git a/src/info.cpp b/src/info.cpp
index 016a4c65261e1ed1819d470434d3906dd28e5881..459894ce487561ff6c328f3a93a07d2ee0d85660 100644
--- a/src/info.cpp
+++ b/src/info.cpp
@@ -327,9 +327,9 @@ void Info::command(int narg, char **arg)
               bstyles[domain->boundary[0][0]],bstyles[domain->boundary[0][1]],
               bstyles[domain->boundary[1][0]],bstyles[domain->boundary[1][1]],
               bstyles[domain->boundary[2][0]],bstyles[domain->boundary[2][1]]);
-      fprintf(out,"Xlo, zhi = %g, %g\n", domain->boxlo[0], domain->boxhi[0]);
-      fprintf(out,"Ylo, zhi = %g, %g\n", domain->boxlo[1], domain->boxhi[1]);
-      fprintf(out,"Zlo, zhi = %g, %g\n", domain->boxlo[2], domain->boxhi[2]);
+      fprintf(out,"xlo, xhi = %g, %g\n", domain->boxlo[0], domain->boxhi[0]);
+      fprintf(out,"ylo, yhi = %g, %g\n", domain->boxlo[1], domain->boxhi[1]);
+      fprintf(out,"zlo, zhi = %g, %g\n", domain->boxlo[2], domain->boxhi[2]);
       if (domain->triclinic)
           fprintf(out,"Xy, xz, yz = %g, %g, %g\n",
                   domain->xy, domain->xz, domain->yz);