diff --git a/src/MOLECULE/angle_charmm.cpp b/src/MOLECULE/angle_charmm.cpp
index 84578bae6e..e69846ad97 100644
--- a/src/MOLECULE/angle_charmm.cpp
+++ b/src/MOLECULE/angle_charmm.cpp
@@ -87,7 +87,7 @@ void AngleCharmm::compute(int eflag, int vflag)
     delx1 = x[i1][0] - x[i2][0];
     dely1 = x[i1][1] - x[i2][1];
     delz1 = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx1,&dely1,&delz1);
+    domain->minimum_image(delx1,dely1,delz1);
 
     rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
     r1 = sqrt(rsq1);
@@ -97,7 +97,7 @@ void AngleCharmm::compute(int eflag, int vflag)
     delx2 = x[i3][0] - x[i2][0];
     dely2 = x[i3][1] - x[i2][1];
     delz2 = x[i3][2] - x[i2][2];
-    domain->minimum_image(&delx2,&dely2,&delz2);
+    domain->minimum_image(delx2,dely2,delz2);
 
     rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
     r2 = sqrt(rsq2);
@@ -107,7 +107,7 @@ void AngleCharmm::compute(int eflag, int vflag)
     delxUB = x[i3][0] - x[i1][0];
     delyUB = x[i3][1] - x[i1][1];
     delzUB = x[i3][2] - x[i1][2];
-    domain->minimum_image(&delxUB,&delyUB,&delzUB);
+    domain->minimum_image(delxUB,delyUB,delzUB);
 
     rsqUB = delxUB*delxUB + delyUB*delyUB + delzUB*delzUB;
     rUB = sqrt(rsqUB);
@@ -285,19 +285,19 @@ double AngleCharmm::single(int type, int i1, int i2, int i3, double rfactor)
   double delx1 = x[i1][0] - x[i2][0];
   double dely1 = x[i1][1] - x[i2][1];
   double delz1 = x[i1][2] - x[i2][2];
-  domain->minimum_image(&delx1,&dely1,&delz1);
+  domain->minimum_image(delx1,dely1,delz1);
   double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
   
   double delx2 = x[i3][0] - x[i2][0];
   double dely2 = x[i3][1] - x[i2][1];
   double delz2 = x[i3][2] - x[i2][2];
-  domain->minimum_image(&delx2,&dely2,&delz2);
+  domain->minimum_image(delx2,dely2,delz2);
   double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
 
   double delxUB = x[i3][0] - x[i1][0];
   double delyUB = x[i3][1] - x[i1][1];
   double delzUB = x[i3][2] - x[i1][2];
-  domain->minimum_image(&delxUB,&delyUB,&delzUB);
+  domain->minimum_image(delxUB,delyUB,delzUB);
   double rUB = sqrt(delxUB*delxUB + delyUB*delyUB + delzUB*delzUB);
 
   double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
diff --git a/src/MOLECULE/angle_cosine.cpp b/src/MOLECULE/angle_cosine.cpp
index 2226fc07c8..0fd2899b31 100644
--- a/src/MOLECULE/angle_cosine.cpp
+++ b/src/MOLECULE/angle_cosine.cpp
@@ -79,7 +79,7 @@ void AngleCosine::compute(int eflag, int vflag)
     delx1 = x[i1][0] - x[i2][0];
     dely1 = x[i1][1] - x[i2][1];
     delz1 = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx1,&dely1,&delz1);
+    domain->minimum_image(delx1,dely1,delz1);
 
     rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
     r1 = sqrt(rsq1);
@@ -89,7 +89,7 @@ void AngleCosine::compute(int eflag, int vflag)
     delx2 = x[i3][0] - x[i2][0];
     dely2 = x[i3][1] - x[i2][1];
     delz2 = x[i3][2] - x[i2][2];
-    domain->minimum_image(&delx2,&dely2,&delz2);
+    domain->minimum_image(delx2,dely2,delz2);
 
     rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
     r2 = sqrt(rsq2);
@@ -227,13 +227,13 @@ double AngleCosine::single(int type, int i1, int i2, int i3, double rfactor)
   double delx1 = x[i1][0] - x[i2][0];
   double dely1 = x[i1][1] - x[i2][1];
   double delz1 = x[i1][2] - x[i2][2];
-  domain->minimum_image(&delx1,&dely1,&delz1);
+  domain->minimum_image(delx1,dely1,delz1);
   double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
   
   double delx2 = x[i3][0] - x[i2][0];
   double dely2 = x[i3][1] - x[i2][1];
   double delz2 = x[i3][2] - x[i2][2];
-  domain->minimum_image(&delx2,&dely2,&delz2);
+  domain->minimum_image(delx2,dely2,delz2);
   double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
 
   double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
diff --git a/src/MOLECULE/angle_cosine_squared.cpp b/src/MOLECULE/angle_cosine_squared.cpp
index 0b353d12cf..9b7fa3bf9d 100644
--- a/src/MOLECULE/angle_cosine_squared.cpp
+++ b/src/MOLECULE/angle_cosine_squared.cpp
@@ -84,7 +84,7 @@ void AngleCosineSquared::compute(int eflag, int vflag)
     delx1 = x[i1][0] - x[i2][0];
     dely1 = x[i1][1] - x[i2][1];
     delz1 = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx1,&dely1,&delz1);
+    domain->minimum_image(delx1,dely1,delz1);
 
     rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
     r1 = sqrt(rsq1);
@@ -94,7 +94,7 @@ void AngleCosineSquared::compute(int eflag, int vflag)
     delx2 = x[i3][0] - x[i2][0];
     dely2 = x[i3][1] - x[i2][1];
     delz2 = x[i3][2] - x[i2][2];
-    domain->minimum_image(&delx2,&dely2,&delz2);
+    domain->minimum_image(delx2,dely2,delz2);
 
     rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
     r2 = sqrt(rsq2);
@@ -248,13 +248,13 @@ double AngleCosineSquared::single(int type, int i1, int i2, int i3,
   double delx1 = x[i1][0] - x[i2][0];
   double dely1 = x[i1][1] - x[i2][1];
   double delz1 = x[i1][2] - x[i2][2];
-  domain->minimum_image(&delx1,&dely1,&delz1);
+  domain->minimum_image(delx1,dely1,delz1);
   double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
   
   double delx2 = x[i3][0] - x[i2][0];
   double dely2 = x[i3][1] - x[i2][1];
   double delz2 = x[i3][2] - x[i2][2];
-  domain->minimum_image(&delx2,&dely2,&delz2);
+  domain->minimum_image(delx2,dely2,delz2);
   double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
 
   double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
diff --git a/src/MOLECULE/angle_harmonic.cpp b/src/MOLECULE/angle_harmonic.cpp
index e209a2de49..e7a4869032 100644
--- a/src/MOLECULE/angle_harmonic.cpp
+++ b/src/MOLECULE/angle_harmonic.cpp
@@ -80,7 +80,7 @@ void AngleHarmonic::compute(int eflag, int vflag)
     delx1 = x[i1][0] - x[i2][0];
     dely1 = x[i1][1] - x[i2][1];
     delz1 = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx1,&dely1,&delz1);
+    domain->minimum_image(delx1,dely1,delz1);
 
     rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
     r1 = sqrt(rsq1);
@@ -90,7 +90,7 @@ void AngleHarmonic::compute(int eflag, int vflag)
     delx2 = x[i3][0] - x[i2][0];
     dely2 = x[i3][1] - x[i2][1];
     delz2 = x[i3][2] - x[i2][2];
-    domain->minimum_image(&delx2,&dely2,&delz2);
+    domain->minimum_image(delx2,dely2,delz2);
 
     rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
     r2 = sqrt(rsq2);
@@ -247,13 +247,13 @@ double AngleHarmonic::single(int type, int i1, int i2, int i3, double rfactor)
   double delx1 = x[i1][0] - x[i2][0];
   double dely1 = x[i1][1] - x[i2][1];
   double delz1 = x[i1][2] - x[i2][2];
-  domain->minimum_image(&delx1,&dely1,&delz1);
+  domain->minimum_image(delx1,dely1,delz1);
   double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
   
   double delx2 = x[i3][0] - x[i2][0];
   double dely2 = x[i3][1] - x[i2][1];
   double delz2 = x[i3][2] - x[i2][2];
-  domain->minimum_image(&delx2,&dely2,&delz2);
+  domain->minimum_image(delx2,dely2,delz2);
   double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
 
   double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
diff --git a/src/MOLECULE/atom_vec_angle.cpp b/src/MOLECULE/atom_vec_angle.cpp
index ec86151237..2d5b2a1377 100644
--- a/src/MOLECULE/atom_vec_angle.cpp
+++ b/src/MOLECULE/atom_vec_angle.cpp
@@ -14,7 +14,6 @@
 #include "stdlib.h"
 #include "atom_vec_angle.h"
 #include "atom.h"
-#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -208,12 +207,13 @@ void AtomVecAngle::copy(int i, int j)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecAngle::pack_comm(int n, int *list, double *buf, int *pbc_flags)
+int AtomVecAngle::pack_comm(int n, int *list, double *buf,
+			    int pbc_flag, double *pbc_dist)
 {
   int i,j,m;
 
   m = 0;
-  if (pbc_flags[0] == 0) {
+  if (pbc_flag == 0) {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = x[j][0];
@@ -221,14 +221,11 @@ int AtomVecAngle::pack_comm(int n, int *list, double *buf, int *pbc_flags)
       buf[m++] = x[j][2];
     }
   } else {
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_flags[1]*xprd;
-      buf[m++] = x[j][1] + pbc_flags[2]*yprd;
-      buf[m++] = x[j][2] + pbc_flags[3]*zprd;
+      buf[m++] = x[j][0] + pbc_dist[0];
+      buf[m++] = x[j][1] + pbc_dist[1];
+      buf[m++] = x[j][2] + pbc_dist[2];
     }
   }
   return m;
@@ -282,12 +279,13 @@ void AtomVecAngle::unpack_reverse(int n, int *list, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecAngle::pack_border(int n, int *list, double *buf, int *pbc_flags)
+int AtomVecAngle::pack_border(int n, int *list, double *buf,
+			      int pbc_flag, double *pbc_dist)
 {
   int i,j,m;
 
   m = 0;
-  if (pbc_flags[0] == 0) {
+  if (pbc_flag == 0) {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = x[j][0];
@@ -299,14 +297,11 @@ int AtomVecAngle::pack_border(int n, int *list, double *buf, int *pbc_flags)
       buf[m++] = molecule[j];
     }
   } else {
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_flags[1]*xprd;
-      buf[m++] = x[j][1] + pbc_flags[2]*yprd;
-      buf[m++] = x[j][2] + pbc_flags[3]*zprd;
+      buf[m++] = x[j][0] + pbc_dist[0];
+      buf[m++] = x[j][1] + pbc_dist[1];
+      buf[m++] = x[j][2] + pbc_dist[2];
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
@@ -589,21 +584,20 @@ int AtomVecAngle::unpack_restart(double *buf)
 }
 
 /* ----------------------------------------------------------------------
-   create one atom of type itype at x0,y0,z0
+   create one atom of itype at coord
    set other values to defaults
 ------------------------------------------------------------------------- */
 
-void AtomVecAngle::create_atom(int itype, double x0, double y0, double z0,
-			       int ihybrid)
+void AtomVecAngle::create_atom(int itype, double *coord, int ihybrid)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = 0;
   type[nlocal] = itype;
-  x[nlocal][0] = x0;
-  x[nlocal][1] = y0;
-  x[nlocal][2] = z0;
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
   mask[nlocal] = 1;
   image[nlocal] = (512 << 20) | (512 << 10) | 512;
   v[nlocal][0] = 0.0;
@@ -622,8 +616,8 @@ void AtomVecAngle::create_atom(int itype, double x0, double y0, double z0,
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-void AtomVecAngle::data_atom(double xtmp, double ytmp, double ztmp,
-			     int imagetmp, char **values, int ihybrid)
+void AtomVecAngle::data_atom(double *coord, int imagetmp, char **values,
+			     int ihybrid)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -638,9 +632,9 @@ void AtomVecAngle::data_atom(double xtmp, double ytmp, double ztmp,
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one("Invalid atom type in Atoms section of data file");
 
-  x[nlocal][0] = xtmp;
-  x[nlocal][1] = ytmp;
-  x[nlocal][2] = ztmp;
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
 
   image[nlocal] = imagetmp;
 
diff --git a/src/MOLECULE/atom_vec_angle.h b/src/MOLECULE/atom_vec_angle.h
index d81a05a998..0ca4f087f3 100644
--- a/src/MOLECULE/atom_vec_angle.h
+++ b/src/MOLECULE/atom_vec_angle.h
@@ -26,11 +26,11 @@ class AtomVecAngle : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int *);
+  int pack_comm(int, int *, double *, int, double *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int *);
+  int pack_border(int, int *, double *, int, double *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
@@ -40,8 +40,8 @@ class AtomVecAngle : public AtomVec {
   int size_restart_one(int);
   int pack_restart(int, double *);
   int unpack_restart(double *);
-  void create_atom(int, double, double, double, int);
-  void data_atom(double, double, double, int, char **, int);
+  void create_atom(int, double *, int);
+  void data_atom(double *, int, char **, int);
   int memory_usage();
 
  private:
diff --git a/src/MOLECULE/atom_vec_bond.cpp b/src/MOLECULE/atom_vec_bond.cpp
index 5fe1144095..22fbd9f1e2 100644
--- a/src/MOLECULE/atom_vec_bond.cpp
+++ b/src/MOLECULE/atom_vec_bond.cpp
@@ -14,7 +14,6 @@
 #include "stdlib.h"
 #include "atom_vec_bond.h"
 #include "atom.h"
-#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -177,12 +176,13 @@ void AtomVecBond::copy(int i, int j)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecBond::pack_comm(int n, int *list, double *buf, int *pbc_flags)
+int AtomVecBond::pack_comm(int n, int *list, double *buf,
+			   int pbc_flag, double *pbc_dist)
 {
   int i,j,m;
 
   m = 0;
-  if (pbc_flags[0] == 0) {
+  if (pbc_flag == 0) {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = x[j][0];
@@ -190,14 +190,11 @@ int AtomVecBond::pack_comm(int n, int *list, double *buf, int *pbc_flags)
       buf[m++] = x[j][2];
     }
   } else {
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_flags[1]*xprd;
-      buf[m++] = x[j][1] + pbc_flags[2]*yprd;
-      buf[m++] = x[j][2] + pbc_flags[3]*zprd;
+      buf[m++] = x[j][0] + pbc_dist[0];
+      buf[m++] = x[j][1] + pbc_dist[1];
+      buf[m++] = x[j][2] + pbc_dist[2];
     }
   }
   return m;
@@ -251,12 +248,13 @@ void AtomVecBond::unpack_reverse(int n, int *list, double *buf)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecBond::pack_border(int n, int *list, double *buf, int *pbc_flags)
+int AtomVecBond::pack_border(int n, int *list, double *buf,
+			     int pbc_flag, double *pbc_dist)
 {
   int i,j,m;
 
   m = 0;
-  if (pbc_flags[0] == 0) {
+  if (pbc_flag == 0) {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = x[j][0];
@@ -268,14 +266,11 @@ int AtomVecBond::pack_border(int n, int *list, double *buf, int *pbc_flags)
       buf[m++] = molecule[j];
     }
   } else {
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_flags[1]*xprd;
-      buf[m++] = x[j][1] + pbc_flags[2]*yprd;
-      buf[m++] = x[j][2] + pbc_flags[3]*zprd;
+      buf[m++] = x[j][0] + pbc_dist[0];
+      buf[m++] = x[j][1] + pbc_dist[1];
+      buf[m++] = x[j][2] + pbc_dist[2];
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
@@ -526,21 +521,20 @@ int AtomVecBond::unpack_restart(double *buf)
 }
 
 /* ----------------------------------------------------------------------
-   create one atom of type itype at x0,y0,z0
+   create one atom of itype at coord
    set other values to defaults
 ------------------------------------------------------------------------- */
 
-void AtomVecBond::create_atom(int itype, double x0, double y0, double z0,
-			      int ihybrid)
+void AtomVecBond::create_atom(int itype, double *coord, int ihybrid)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = 0;
   type[nlocal] = itype;
-  x[nlocal][0] = x0;
-  x[nlocal][1] = y0;
-  x[nlocal][2] = z0;
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
   mask[nlocal] = 1;
   image[nlocal] = (512 << 20) | (512 << 10) | 512;
   v[nlocal][0] = 0.0;
@@ -558,8 +552,8 @@ void AtomVecBond::create_atom(int itype, double x0, double y0, double z0,
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-void AtomVecBond::data_atom(double xtmp, double ytmp, double ztmp,
-			    int imagetmp, char **values, int ihybrid)
+void AtomVecBond::data_atom(double *coord, int imagetmp, char **values,
+			      int ihybrid)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -574,9 +568,9 @@ void AtomVecBond::data_atom(double xtmp, double ytmp, double ztmp,
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one("Invalid atom type in Atoms section of data file");
 
-  x[nlocal][0] = xtmp;
-  x[nlocal][1] = ytmp;
-  x[nlocal][2] = ztmp;
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
 
   image[nlocal] = imagetmp;
 
diff --git a/src/MOLECULE/atom_vec_bond.h b/src/MOLECULE/atom_vec_bond.h
index 96f69b98ce..9da8f73fd7 100644
--- a/src/MOLECULE/atom_vec_bond.h
+++ b/src/MOLECULE/atom_vec_bond.h
@@ -26,11 +26,11 @@ class AtomVecBond : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int *);
+  int pack_comm(int, int *, double *, int, double *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int *);
+  int pack_border(int, int *, double *, int, double *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
@@ -40,8 +40,8 @@ class AtomVecBond : public AtomVec {
   int size_restart_one(int);
   int pack_restart(int, double *);
   int unpack_restart(double *);
-  void create_atom(int, double, double, double, int);
-  void data_atom(double, double, double, int, char **, int);
+  void create_atom(int, double *, int);
+  void data_atom(double *, int, char **, int);
   int memory_usage();
 
  private:
diff --git a/src/MOLECULE/atom_vec_full.cpp b/src/MOLECULE/atom_vec_full.cpp
index 7de7dfb748..2fa18dbd5a 100644
--- a/src/MOLECULE/atom_vec_full.cpp
+++ b/src/MOLECULE/atom_vec_full.cpp
@@ -14,7 +14,6 @@
 #include "stdlib.h"
 #include "atom_vec_full.h"
 #include "atom.h"
-#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -286,12 +285,13 @@ void AtomVecFull::copy(int i, int j)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecFull::pack_comm(int n, int *list, double *buf, int *pbc_flags)
+int AtomVecFull::pack_comm(int n, int *list, double *buf,
+			   int pbc_flag, double *pbc_dist)
 {
   int i,j,m;
 
   m = 0;
-  if (pbc_flags[0] == 0) {
+  if (pbc_flag == 0) {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = x[j][0];
@@ -299,14 +299,11 @@ int AtomVecFull::pack_comm(int n, int *list, double *buf, int *pbc_flags)
       buf[m++] = x[j][2];
     }
   } else {
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_flags[1]*xprd;
-      buf[m++] = x[j][1] + pbc_flags[2]*yprd;
-      buf[m++] = x[j][2] + pbc_flags[3]*zprd;
+      buf[m++] = x[j][0] + pbc_dist[0];
+      buf[m++] = x[j][1] + pbc_dist[1];
+      buf[m++] = x[j][2] + pbc_dist[2];
     }
   }
   return m;
@@ -361,12 +358,12 @@ void AtomVecFull::unpack_reverse(int n, int *list, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecFull::pack_border(int n, int *list, double *buf,
-				  int *pbc_flags)
+			     int pbc_flag, double *pbc_dist)
 {
   int i,j,m;
 
   m = 0;
-  if (pbc_flags[0] == 0) {
+  if (pbc_flag == 0) {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = x[j][0];
@@ -379,14 +376,11 @@ int AtomVecFull::pack_border(int n, int *list, double *buf,
       buf[m++] = molecule[j];
     }
   } else {
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_flags[1]*xprd;
-      buf[m++] = x[j][1] + pbc_flags[2]*yprd;
-      buf[m++] = x[j][2] + pbc_flags[3]*zprd;
+      buf[m++] = x[j][0] + pbc_dist[0];
+      buf[m++] = x[j][1] + pbc_dist[1];
+      buf[m++] = x[j][2] + pbc_dist[2];
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
@@ -751,21 +745,20 @@ int AtomVecFull::unpack_restart(double *buf)
 }
 
 /* ----------------------------------------------------------------------
-   create one atom of type itype at x0,y0,z0
+   create one atom of itype at coord
    set other values to defaults
 ------------------------------------------------------------------------- */
 
-void AtomVecFull::create_atom(int itype, double x0, double y0, double z0,
-				   int ihybrid)
+void AtomVecFull::create_atom(int itype, double *coord, int ihybrid)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = 0;
   type[nlocal] = itype;
-  x[nlocal][0] = x0;
-  x[nlocal][1] = y0;
-  x[nlocal][2] = z0;
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
   mask[nlocal] = 1;
   image[nlocal] = (512 << 20) | (512 << 10) | 512;
   v[nlocal][0] = 0.0;
@@ -787,8 +780,8 @@ void AtomVecFull::create_atom(int itype, double x0, double y0, double z0,
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-void AtomVecFull::data_atom(double xtmp, double ytmp, double ztmp,
-				 int imagetmp, char **values, int ihybrid)
+void AtomVecFull::data_atom(double *coord, int imagetmp, char **values,
+			    int ihybrid)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -805,9 +798,9 @@ void AtomVecFull::data_atom(double xtmp, double ytmp, double ztmp,
 
   q[nlocal] = atof(values[3]);
 
-  x[nlocal][0] = xtmp;
-  x[nlocal][1] = ytmp;
-  x[nlocal][2] = ztmp;
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
 
   image[nlocal] = imagetmp;
 
diff --git a/src/MOLECULE/atom_vec_full.h b/src/MOLECULE/atom_vec_full.h
index 81f2a6fcca..b4ed270763 100644
--- a/src/MOLECULE/atom_vec_full.h
+++ b/src/MOLECULE/atom_vec_full.h
@@ -26,11 +26,11 @@ class AtomVecFull : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int *);
+  int pack_comm(int, int *, double *, int, double *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int *);
+  int pack_border(int, int *, double *, int, double *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
@@ -40,8 +40,8 @@ class AtomVecFull : public AtomVec {
   int size_restart_one(int);
   int pack_restart(int, double *);
   int unpack_restart(double *);
-  void create_atom(int, double, double, double, int);
-  void data_atom(double, double, double, int, char **, int);
+  void create_atom(int, double *, int);
+  void data_atom(double *, int, char **, int);
   int memory_usage();
 
  private:
diff --git a/src/MOLECULE/atom_vec_molecular.cpp b/src/MOLECULE/atom_vec_molecular.cpp
index d86eede108..1f6488fa73 100644
--- a/src/MOLECULE/atom_vec_molecular.cpp
+++ b/src/MOLECULE/atom_vec_molecular.cpp
@@ -14,7 +14,6 @@
 #include "stdlib.h"
 #include "atom_vec_molecular.h"
 #include "atom.h"
-#include "domain.h"
 #include "modify.h"
 #include "fix.h"
 #include "memory.h"
@@ -280,12 +279,13 @@ void AtomVecMolecular::copy(int i, int j)
 
 /* ---------------------------------------------------------------------- */
 
-int AtomVecMolecular::pack_comm(int n, int *list, double *buf, int *pbc_flags)
+int AtomVecMolecular::pack_comm(int n, int *list, double *buf, 
+				int pbc_flag, double *pbc_dist)
 {
   int i,j,m;
 
   m = 0;
-  if (pbc_flags[0] == 0) {
+  if (pbc_flag == 0) {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = x[j][0];
@@ -293,14 +293,11 @@ int AtomVecMolecular::pack_comm(int n, int *list, double *buf, int *pbc_flags)
       buf[m++] = x[j][2];
     }
   } else {
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_flags[1]*xprd;
-      buf[m++] = x[j][1] + pbc_flags[2]*yprd;
-      buf[m++] = x[j][2] + pbc_flags[3]*zprd;
+      buf[m++] = x[j][0] + pbc_dist[0];
+      buf[m++] = x[j][1] + pbc_dist[1];
+      buf[m++] = x[j][2] + pbc_dist[2];
     }
   }
   return m;
@@ -355,12 +352,12 @@ void AtomVecMolecular::unpack_reverse(int n, int *list, double *buf)
 /* ---------------------------------------------------------------------- */
 
 int AtomVecMolecular::pack_border(int n, int *list, double *buf,
-				  int *pbc_flags)
+				  int pbc_flag, double *pbc_dist)
 {
   int i,j,m;
 
   m = 0;
-  if (pbc_flags[0] == 0) {
+  if (pbc_flag == 0) {
     for (i = 0; i < n; i++) {
       j = list[i];
       buf[m++] = x[j][0];
@@ -372,14 +369,11 @@ int AtomVecMolecular::pack_border(int n, int *list, double *buf,
       buf[m++] = molecule[j];
     }
   } else {
-    double xprd = domain->xprd;
-    double yprd = domain->yprd;
-    double zprd = domain->zprd;
     for (i = 0; i < n; i++) {
       j = list[i];
-      buf[m++] = x[j][0] + pbc_flags[1]*xprd;
-      buf[m++] = x[j][1] + pbc_flags[2]*yprd;
-      buf[m++] = x[j][2] + pbc_flags[3]*zprd;
+      buf[m++] = x[j][0] + pbc_dist[0];
+      buf[m++] = x[j][1] + pbc_dist[1];
+      buf[m++] = x[j][2] + pbc_dist[2];
       buf[m++] = tag[j];
       buf[m++] = type[j];
       buf[m++] = mask[j];
@@ -736,21 +730,20 @@ int AtomVecMolecular::unpack_restart(double *buf)
 }
 
 /* ----------------------------------------------------------------------
-   create one atom of type itype at x0,y0,z0
+   create one atom of itype at coord
    set other values to defaults
 ------------------------------------------------------------------------- */
 
-void AtomVecMolecular::create_atom(int itype, double x0, double y0, double z0,
-				   int ihybrid)
+void AtomVecMolecular::create_atom(int itype, double *coord, int ihybrid)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
 
   tag[nlocal] = 0;
   type[nlocal] = itype;
-  x[nlocal][0] = x0;
-  x[nlocal][1] = y0;
-  x[nlocal][2] = z0;
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
   mask[nlocal] = 1;
   image[nlocal] = (512 << 20) | (512 << 10) | 512;
   v[nlocal][0] = 0.0;
@@ -771,8 +764,8 @@ void AtomVecMolecular::create_atom(int itype, double x0, double y0, double z0,
    initialize other atom quantities
 ------------------------------------------------------------------------- */
 
-void AtomVecMolecular::data_atom(double xtmp, double ytmp, double ztmp,
-				 int imagetmp, char **values, int ihybrid)
+void AtomVecMolecular::data_atom(double *coord, int imagetmp, char **values,
+				 int ihybrid)
 {
   int nlocal = atom->nlocal;
   if (nlocal == nmax) grow(0);
@@ -787,9 +780,9 @@ void AtomVecMolecular::data_atom(double xtmp, double ytmp, double ztmp,
   if (type[nlocal] <= 0 || type[nlocal] > atom->ntypes)
     error->one("Invalid atom type in Atoms section of data file");
 
-  x[nlocal][0] = xtmp;
-  x[nlocal][1] = ytmp;
-  x[nlocal][2] = ztmp;
+  x[nlocal][0] = coord[0];
+  x[nlocal][1] = coord[1];
+  x[nlocal][2] = coord[2];
 
   image[nlocal] = imagetmp;
 
diff --git a/src/MOLECULE/atom_vec_molecular.h b/src/MOLECULE/atom_vec_molecular.h
index 795e62418d..2b2e3a49fe 100644
--- a/src/MOLECULE/atom_vec_molecular.h
+++ b/src/MOLECULE/atom_vec_molecular.h
@@ -26,11 +26,11 @@ class AtomVecMolecular : public AtomVec {
   void zero_owned(int);
   void zero_ghost(int, int);
   void copy(int, int);
-  int pack_comm(int, int *, double *, int *);
+  int pack_comm(int, int *, double *, int, double *);
   void unpack_comm(int, int, double *);
   int pack_reverse(int, int, double *);
   void unpack_reverse(int, int *, double *);
-  int pack_border(int, int *, double *, int *);
+  int pack_border(int, int *, double *, int, double *);
   int pack_border_one(int, double *);
   void unpack_border(int, int, double *);
   int unpack_border_one(int, double *);
@@ -40,8 +40,8 @@ class AtomVecMolecular : public AtomVec {
   int size_restart_one(int);
   int pack_restart(int, double *);
   int unpack_restart(double *);
-  void create_atom(int, double, double, double, int);
-  void data_atom(double, double, double, int, char **, int);
+  void create_atom(int, double *, int);
+  void data_atom(double *, int, char **, int);
   int memory_usage();
 
  private:
diff --git a/src/MOLECULE/bond_fene.cpp b/src/MOLECULE/bond_fene.cpp
index 3369dfef25..534ea9d27e 100644
--- a/src/MOLECULE/bond_fene.cpp
+++ b/src/MOLECULE/bond_fene.cpp
@@ -79,7 +79,7 @@ void BondFENE::compute(int eflag, int vflag)
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
     delz = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx,&dely,&delz);
+    domain->minimum_image(delx,dely,delz);
 
     // force from log term
 
diff --git a/src/MOLECULE/bond_fene_expand.cpp b/src/MOLECULE/bond_fene_expand.cpp
index c9a4f738e0..e8900a57ca 100644
--- a/src/MOLECULE/bond_fene_expand.cpp
+++ b/src/MOLECULE/bond_fene_expand.cpp
@@ -81,7 +81,7 @@ void BondFENEExpand::compute(int eflag, int vflag)
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
     delz = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx,&dely,&delz);
+    domain->minimum_image(delx,dely,delz);
 
     // force from log term
 
diff --git a/src/MOLECULE/bond_harmonic.cpp b/src/MOLECULE/bond_harmonic.cpp
index dc56435432..c839a5a942 100644
--- a/src/MOLECULE/bond_harmonic.cpp
+++ b/src/MOLECULE/bond_harmonic.cpp
@@ -73,7 +73,7 @@ void BondHarmonic::compute(int eflag, int vflag)
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
     delz = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx,&dely,&delz);
+    domain->minimum_image(delx,dely,delz);
 
     rsq = delx*delx + dely*dely + delz*delz;
     r = sqrt(rsq);
diff --git a/src/MOLECULE/bond_morse.cpp b/src/MOLECULE/bond_morse.cpp
index 355206cae5..2c3dbee48c 100644
--- a/src/MOLECULE/bond_morse.cpp
+++ b/src/MOLECULE/bond_morse.cpp
@@ -78,7 +78,7 @@ void BondMorse::compute(int eflag, int vflag)
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
     delz = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx,&dely,&delz);
+    domain->minimum_image(delx,dely,delz);
 
     rsq = delx*delx + dely*dely + delz*delz;
     r = sqrt(rsq);
diff --git a/src/MOLECULE/bond_nonlinear.cpp b/src/MOLECULE/bond_nonlinear.cpp
index ee7f47252d..9a801d32e8 100644
--- a/src/MOLECULE/bond_nonlinear.cpp
+++ b/src/MOLECULE/bond_nonlinear.cpp
@@ -74,7 +74,7 @@ void BondNonlinear::compute(int eflag, int vflag)
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
     delz = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx,&dely,&delz);
+    domain->minimum_image(delx,dely,delz);
 
     rsq = delx*delx + dely*dely + delz*delz;
     r = sqrt(rsq);
diff --git a/src/MOLECULE/bond_quartic.cpp b/src/MOLECULE/bond_quartic.cpp
index 5d958f7da1..b1b52549fb 100755
--- a/src/MOLECULE/bond_quartic.cpp
+++ b/src/MOLECULE/bond_quartic.cpp
@@ -92,7 +92,7 @@ void BondQuartic::compute(int eflag, int vflag)
     delx = x[i1][0] - x[i2][0];
     dely = x[i1][1] - x[i2][1];
     delz = x[i1][2] - x[i2][2];
-    domain->minimum_image(&delx,&dely,&delz);
+    domain->minimum_image(delx,dely,delz);
 
     rsq = delx*delx + dely*dely + delz*delz;
 
diff --git a/src/MOLECULE/dihedral_charmm.cpp b/src/MOLECULE/dihedral_charmm.cpp
index 7587fa705e..4741bd29f5 100644
--- a/src/MOLECULE/dihedral_charmm.cpp
+++ b/src/MOLECULE/dihedral_charmm.cpp
@@ -106,26 +106,26 @@ void DihedralCharmm::compute(int eflag, int vflag)
     vb1x = x[i1][0] - x[i2][0];
     vb1y = x[i1][1] - x[i2][1];
     vb1z = x[i1][2] - x[i2][2];
-    domain->minimum_image(&vb1x,&vb1y,&vb1z);
+    domain->minimum_image(vb1x,vb1y,vb1z);
 
     // 2nd bond
 
     vb2x = x[i3][0] - x[i2][0];
     vb2y = x[i3][1] - x[i2][1];
     vb2z = x[i3][2] - x[i2][2];
-    domain->minimum_image(&vb2x,&vb2y,&vb2z);
+    domain->minimum_image(vb2x,vb2y,vb2z);
 
     vb2xm = -vb2x;
     vb2ym = -vb2y;
     vb2zm = -vb2z;
-    domain->minimum_image(&vb2xm,&vb2ym,&vb2zm);
+    domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
     // 3rd bond
 
     vb3x = x[i4][0] - x[i3][0];
     vb3y = x[i4][1] - x[i3][1];
     vb3z = x[i4][2] - x[i3][2];
-    domain->minimum_image(&vb3x,&vb3y,&vb3z);
+    domain->minimum_image(vb3x,vb3y,vb3z);
     
     ax = vb1y*vb2zm - vb1z*vb2ym;
     ay = vb1z*vb2xm - vb1x*vb2zm;
@@ -270,7 +270,7 @@ void DihedralCharmm::compute(int eflag, int vflag)
       delx = x[i1][0] - x[i4][0];
       dely = x[i1][1] - x[i4][1];
       delz = x[i1][2] - x[i4][2];
-      domain->minimum_image(&delx,&dely,&delz);
+      domain->minimum_image(delx,dely,delz);
       rsq = delx*delx + dely*dely + delz*delz;
       r2inv = 1.0/rsq;
       r6inv = r2inv*r2inv*r2inv;
diff --git a/src/MOLECULE/dihedral_harmonic.cpp b/src/MOLECULE/dihedral_harmonic.cpp
index 6034c3d60c..98a8a94fd0 100644
--- a/src/MOLECULE/dihedral_harmonic.cpp
+++ b/src/MOLECULE/dihedral_harmonic.cpp
@@ -97,26 +97,26 @@ void DihedralHarmonic::compute(int eflag, int vflag)
     vb1x = x[i1][0] - x[i2][0];
     vb1y = x[i1][1] - x[i2][1];
     vb1z = x[i1][2] - x[i2][2];
-    domain->minimum_image(&vb1x,&vb1y,&vb1z);
+    domain->minimum_image(vb1x,vb1y,vb1z);
 
     // 2nd bond
 
     vb2x = x[i3][0] - x[i2][0];
     vb2y = x[i3][1] - x[i2][1];
     vb2z = x[i3][2] - x[i2][2];
-    domain->minimum_image(&vb2x,&vb2y,&vb2z);
+    domain->minimum_image(vb2x,vb2y,vb2z);
 
     vb2xm = -vb2x;
     vb2ym = -vb2y;
     vb2zm = -vb2z;
-    domain->minimum_image(&vb2xm,&vb2ym,&vb2zm);
+    domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
     // 3rd bond
 
     vb3x = x[i4][0] - x[i3][0];
     vb3y = x[i4][1] - x[i3][1];
     vb3z = x[i4][2] - x[i3][2];
-    domain->minimum_image(&vb3x,&vb3y,&vb3z);
+    domain->minimum_image(vb3x,vb3y,vb3z);
     
     // c,s calculation
 
diff --git a/src/MOLECULE/dihedral_helix.cpp b/src/MOLECULE/dihedral_helix.cpp
index bf1b601a0b..0ca292eeab 100644
--- a/src/MOLECULE/dihedral_helix.cpp
+++ b/src/MOLECULE/dihedral_helix.cpp
@@ -101,26 +101,26 @@ void DihedralHelix::compute(int eflag, int vflag)
     vb1x = x[i1][0] - x[i2][0];
     vb1y = x[i1][1] - x[i2][1];
     vb1z = x[i1][2] - x[i2][2];
-    domain->minimum_image(&vb1x,&vb1y,&vb1z);
+    domain->minimum_image(vb1x,vb1y,vb1z);
 
     // 2nd bond
 
     vb2x = x[i3][0] - x[i2][0];
     vb2y = x[i3][1] - x[i2][1];
     vb2z = x[i3][2] - x[i2][2];
-    domain->minimum_image(&vb2x,&vb2y,&vb2z);
+    domain->minimum_image(vb2x,vb2y,vb2z);
 
     vb2xm = -vb2x;
     vb2ym = -vb2y;
     vb2zm = -vb2z;
-    domain->minimum_image(&vb2xm,&vb2ym,&vb2zm);
+    domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
     // 3rd bond
 
     vb3x = x[i4][0] - x[i3][0];
     vb3y = x[i4][1] - x[i3][1];
     vb3z = x[i4][2] - x[i3][2];
-    domain->minimum_image(&vb3x,&vb3y,&vb3z);
+    domain->minimum_image(vb3x,vb3y,vb3z);
 
     // c0 calculation
         
diff --git a/src/MOLECULE/dihedral_multi_harmonic.cpp b/src/MOLECULE/dihedral_multi_harmonic.cpp
index 864625fc38..1d047ffde0 100644
--- a/src/MOLECULE/dihedral_multi_harmonic.cpp
+++ b/src/MOLECULE/dihedral_multi_harmonic.cpp
@@ -100,26 +100,26 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag)
     vb1x = x[i1][0] - x[i2][0];
     vb1y = x[i1][1] - x[i2][1];
     vb1z = x[i1][2] - x[i2][2];
-    domain->minimum_image(&vb1x,&vb1y,&vb1z);
+    domain->minimum_image(vb1x,vb1y,vb1z);
 
     // 2nd bond
 
     vb2x = x[i3][0] - x[i2][0];
     vb2y = x[i3][1] - x[i2][1];
     vb2z = x[i3][2] - x[i2][2];
-    domain->minimum_image(&vb2x,&vb2y,&vb2z);
+    domain->minimum_image(vb2x,vb2y,vb2z);
 
     vb2xm = -vb2x;
     vb2ym = -vb2y;
     vb2zm = -vb2z;
-    domain->minimum_image(&vb2xm,&vb2ym,&vb2zm);
+    domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
     // 3rd bond
 
     vb3x = x[i4][0] - x[i3][0];
     vb3y = x[i4][1] - x[i3][1];
     vb3z = x[i4][2] - x[i3][2];
-    domain->minimum_image(&vb3x,&vb3y,&vb3z);
+    domain->minimum_image(vb3x,vb3y,vb3z);
 
     // c0 calculation
         
diff --git a/src/MOLECULE/dihedral_opls.cpp b/src/MOLECULE/dihedral_opls.cpp
index efb84c41d6..d50602f1c7 100644
--- a/src/MOLECULE/dihedral_opls.cpp
+++ b/src/MOLECULE/dihedral_opls.cpp
@@ -100,26 +100,26 @@ void DihedralOPLS::compute(int eflag, int vflag)
     vb1x = x[i1][0] - x[i2][0];
     vb1y = x[i1][1] - x[i2][1];
     vb1z = x[i1][2] - x[i2][2];
-    domain->minimum_image(&vb1x,&vb1y,&vb1z);
+    domain->minimum_image(vb1x,vb1y,vb1z);
 
     // 2nd bond
 
     vb2x = x[i3][0] - x[i2][0];
     vb2y = x[i3][1] - x[i2][1];
     vb2z = x[i3][2] - x[i2][2];
-    domain->minimum_image(&vb2x,&vb2y,&vb2z);
+    domain->minimum_image(vb2x,vb2y,vb2z);
 
     vb2xm = -vb2x;
     vb2ym = -vb2y;
     vb2zm = -vb2z;
-    domain->minimum_image(&vb2xm,&vb2ym,&vb2zm);
+    domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
     // 3rd bond
 
     vb3x = x[i4][0] - x[i3][0];
     vb3y = x[i4][1] - x[i3][1];
     vb3z = x[i4][2] - x[i3][2];
-    domain->minimum_image(&vb3x,&vb3y,&vb3z);
+    domain->minimum_image(vb3x,vb3y,vb3z);
 
     // c0 calculation
 
diff --git a/src/MOLECULE/improper_cvff.cpp b/src/MOLECULE/improper_cvff.cpp
index ad344d5901..081892e105 100644
--- a/src/MOLECULE/improper_cvff.cpp
+++ b/src/MOLECULE/improper_cvff.cpp
@@ -92,26 +92,26 @@ void ImproperCvff::compute(int eflag, int vflag)
     vb1x = x[i1][0] - x[i2][0];
     vb1y = x[i1][1] - x[i2][1];
     vb1z = x[i1][2] - x[i2][2];
-    domain->minimum_image(&vb1x,&vb1y,&vb1z);
+    domain->minimum_image(vb1x,vb1y,vb1z);
 
     // 2nd bond
 
     vb2x = x[i3][0] - x[i2][0];
     vb2y = x[i3][1] - x[i2][1];
     vb2z = x[i3][2] - x[i2][2];
-    domain->minimum_image(&vb2x,&vb2y,&vb2z);
+    domain->minimum_image(vb2x,vb2y,vb2z);
 
     vb2xm = -vb2x;
     vb2ym = -vb2y;
     vb2zm = -vb2z;
-    domain->minimum_image(&vb2xm,&vb2ym,&vb2zm);
+    domain->minimum_image(vb2xm,vb2ym,vb2zm);
 
     // 3rd bond
 
     vb3x = x[i4][0] - x[i3][0];
     vb3y = x[i4][1] - x[i3][1];
     vb3z = x[i4][2] - x[i3][2];
-    domain->minimum_image(&vb3x,&vb3y,&vb3z);
+    domain->minimum_image(vb3x,vb3y,vb3z);
 
     // c0 calculation
         
diff --git a/src/MOLECULE/improper_harmonic.cpp b/src/MOLECULE/improper_harmonic.cpp
index b8af76dc8c..e4a7fca1b4 100644
--- a/src/MOLECULE/improper_harmonic.cpp
+++ b/src/MOLECULE/improper_harmonic.cpp
@@ -88,17 +88,17 @@ void ImproperHarmonic::compute(int eflag, int vflag)
     v1x = x[i2][0] - x[i1][0];
     v1y = x[i2][1] - x[i1][1];
     v1z = x[i2][2] - x[i1][2];
-    domain->minimum_image(&v1x,&v1y,&v1z);
+    domain->minimum_image(v1x,v1y,v1z);
 
     v2x = x[i3][0] - x[i2][0];
     v2y = x[i3][1] - x[i2][1];
     v2z = x[i3][2] - x[i2][2];
-    domain->minimum_image(&v2x,&v2y,&v2z);
+    domain->minimum_image(v2x,v2y,v2z);
 
     v3x = x[i4][0] - x[i3][0];
     v3y = x[i4][1] - x[i3][1];
     v3z = x[i4][2] - x[i3][2];
-    domain->minimum_image(&v3x,&v3y,&v3z);
+    domain->minimum_image(v3x,v3y,v3z);
 
     ss1 = 1.0 / (v1x*v1x + v1y*v1y + v1z*v1z);
     ss2 = 1.0 / (v2x*v2x + v2y*v2y + v2z*v2z);