diff --git a/src/core/ActionWithVirtualAtom.cpp b/src/core/ActionWithVirtualAtom.cpp
index f1f71a8dde2496f6434e9f8c4139d5cf3f552e45..3e9ffd21a0456ec9d3918943cf2af258bda85b60 100644
--- a/src/core/ActionWithVirtualAtom.cpp
+++ b/src/core/ActionWithVirtualAtom.cpp
@@ -86,6 +86,21 @@ void ActionWithVirtualAtom::setBoxDerivatives(const std::vector<Tensor> &d){
   for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) boxDerivatives[j][i][j]+=pos[i];
 }
 
+void ActionWithVirtualAtom::setBoxDerivativesNoPbc(){
+  std::vector<Tensor> bd(3);
+  for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) for(unsigned k=0;k<3;k++){
+// Notice that this expression is very similar to the one used in Colvar::setBoxDerivativesNoPbc().
+// Indeed, we have the negative of a sum over dependent atoms (l) of the external product between positions
+// and derivatives. Notice that this only works only when Pbc have not been used to compute
+// derivatives.
+    for(unsigned l=0;l<getNumberOfAtoms();l++){
+      bd[k][i][j]-=getPosition(l)[i]*derivatives[l][j][k];
+    }
+  }
+  setBoxDerivatives(bd);
+}
+
+
 
 void ActionWithVirtualAtom::setGradientsIfNeeded(){
   if(isOptionOn("GRADIENTS")) { 
diff --git a/src/core/ActionWithVirtualAtom.h b/src/core/ActionWithVirtualAtom.h
index 69397e98601d2a88f2769e970984958e08317e64..10925335e5870fffd8dd570220160cfeee14cb14 100644
--- a/src/core/ActionWithVirtualAtom.h
+++ b/src/core/ActionWithVirtualAtom.h
@@ -66,6 +66,13 @@ protected:
 /// On the other hand if the vatom position is a non-linear function of atomic coordinates this
 /// should be called (see vatom::Ghost).
   void setBoxDerivatives(const std::vector<Tensor> &d);
+/// Set box derivatives automatically.
+/// It should be called after the settomsDerivatives has been used for all
+/// single atoms.
+/// \warning It only works for virtual atoms NOT using PBCs!
+///          This implies that all atoms used + the new virtual atom should be
+///          in the same periodic image.
+  void setBoxDerivativesNoPbc();
 public:
   void setGradients();
   const std::map<AtomNumber,Tensor> & getGradients()const;
diff --git a/src/vatom/Ghost.cpp b/src/vatom/Ghost.cpp
index 9d123eed31b017c12ece3892d6883179d6acb361..b434e591f23dcee8b601e5840108decae280dd90 100644
--- a/src/vatom/Ghost.cpp
+++ b/src/vatom/Ghost.cpp
@@ -90,11 +90,6 @@ void Ghost::calculate(){
   Vector pos;
   vector<Tensor> deriv(getNumberOfAtoms());
   vector<Vector> n;
-  Vector pp[3];
-
-  pp[0]=getPosition(0);
-  pp[1]=pp[0]+delta(getPosition(0), getPosition(1));
-  pp[2]=pp[0]+delta(getPosition(0), getPosition(2));
 
 // first versor 
   Vector n01 = delta(getPosition(0), getPosition(1));
@@ -184,17 +179,7 @@ void Ghost::calculate(){
   setAtomsDerivatives(deriv);
 
 // Virial contribution
-  std::vector<Tensor> bd(3);
-  for(unsigned i=0;i<3;i++) for(unsigned j=0;j<3;j++) for(unsigned k=0;k<3;k++){
-// Notice that this expression is very similar to the one used in Colvar::setBoxDerivativesNoPbc().
-// Indeed, we have the negative of a sum over dependent atoms (l) of the external product between positions
-// and derivatives. Notice that positions here should be computed in the same periodic image
-// as pos.
-    for(unsigned l=0;l<3;l++){
-      bd[k][i][j]-=pp[l][i]*deriv[l][j][k];
-    }
-  }
-  setBoxDerivatives(bd);
+  setBoxDerivativesNoPbc();
 }
 
 }