diff --git a/CHANGES/v2.1.txt b/CHANGES/v2.1.txt
index e88da95dde0289c9b52d16f91ff61de7767f6c7e..58e35a00931259cc194d2c69656fbb9a49cb3995 100644
--- a/CHANGES/v2.1.txt
+++ b/CHANGES/v2.1.txt
@@ -189,6 +189,7 @@ For users:
   See notes <a href="http://github.com/plumed/plumed2/issues/132"> here </a>
   This fix requires gromacs to be repatched and could be very important if you run biased simulations in the NPT ensemble.
 - Fixed a bug in the virial computed with \ref FIT_TO_TEMPLATE when the reference pdb had center non located at the origin.
+- Fixed a bug in the the forces computed with \ref FIT_TO_TEMPLATE when used in combination with \ref COM, \ref CENTER, or \ref GHOST
 - Updated internal molfile plugins to VMD 1.9.2.
 - Included crd and crdbox formats to internal molfile.
 - Added --natoms to \ref driver . This is required to read coordinate
diff --git a/regtest/basic/rt63/plumed.dat b/regtest/basic/rt63/plumed.dat
index 3c2bd175e213d789d71559757a5a9edf5e7d80d9..5abd34da8322fc9fb4fcc6731e7d7aa5b8e399d0 100644
--- a/regtest/basic/rt63/plumed.dat
+++ b/regtest/basic/rt63/plumed.dat
@@ -7,7 +7,8 @@ MOLINFO STRUCTURE=helix.pdb
 # ps: DISTANCE ATOMS=c,3 SCALED_COMPONENTS
 # new way (resetting reference frame)
 FIT_TO_TEMPLATE REFERENCE=align.pdb TYPE=SIMPLE
-p: POSITION ATOM=3
+s: COM ATOMS=3
+p: POSITION ATOM=s
 ps: POSITION SCALED_COMPONENTS ATOM=3
 # using one or the other should give the same result and same forces
 
diff --git a/src/core/ActionWithVirtualAtom.cpp b/src/core/ActionWithVirtualAtom.cpp
index ae05fb724dc4f3bc117813b6f277b38d93971377..f1f71a8dde2496f6434e9f8c4139d5cf3f552e45 100644
--- a/src/core/ActionWithVirtualAtom.cpp
+++ b/src/core/ActionWithVirtualAtom.cpp
@@ -46,10 +46,13 @@ ActionWithVirtualAtom::~ActionWithVirtualAtom(){
 }
 
 void ActionWithVirtualAtom::apply(){
-  const Vector & f(atoms.forces[index.index()]);
+  Vector & f(atoms.forces[index.index()]);
   for(unsigned i=0;i<getNumberOfAtoms();i++) modifyForces()[i]=matmul(derivatives[i],f);
   Tensor & v(modifyVirial());
   for(unsigned i=0;i<3;i++) v+=boxDerivatives[i]*f[i];
+  f.zero(); // after propagating the force to the atoms used to compute the vatom, we reset this to zero
+            // this is necessary to avoid double counting if then one tries to compute the total force on the c.o.m. of the system.
+            // notice that this is currently done in FIT_TO_TEMPLATE
 }
 
 void ActionWithVirtualAtom::requestAtoms(const std::vector<AtomNumber> & a){