diff --git a/CHANGES/v2.3.txt b/CHANGES/v2.3.txt
index 73dac9a448361a2246c8943172460e34be810a48..84294d28952b0881d02e06728936c898b5f4765a 100644
--- a/CHANGES/v2.3.txt
+++ b/CHANGES/v2.3.txt
@@ -137,6 +137,7 @@ Version 2.3.1 (Mar 31, 2017)
 - fixed `plumed-config` that was not working.
 - log file points to the `config.txt` files to allow users to check which features were available in that compiled version.
 - `make clean` in root dir now also cleans `vim` subdirectory.
+- Resolved problem with nan in SMAC with SPECIESA and SPECIESB involving molecules that are the same
 - Updated gromacs patch to version 2016.3 
 
 For developers:
diff --git a/regtest/crystallization/rt-smac/colvar2.reference b/regtest/crystallization/rt-smac/colvar2.reference
new file mode 100644
index 0000000000000000000000000000000000000000..e62860edc84192c3c2251dd72be3ca1e0ba07ab5
--- /dev/null
+++ b/regtest/crystallization/rt-smac/colvar2.reference
@@ -0,0 +1,12 @@
+#! FIELDS time s3.morethan
+ 0.000000   0.1436
+ 1.000000   0.1815
+ 2.000000   0.1786
+ 3.000000   0.3008
+ 4.000000   0.1060
+ 5.000000   0.2212
+ 6.000000   0.2735
+ 7.000000   0.2153
+ 8.000000   0.2188
+ 9.000000   0.1001
+ 10.000000   0.2119
diff --git a/regtest/crystallization/rt-smac/plumed.dat b/regtest/crystallization/rt-smac/plumed.dat
index bfab738c22eb932cc73e12ab6677e3159eacdcd6..3b4f348f227cf152760c3334a7fb14daba3e9313 100644
--- a/regtest/crystallization/rt-smac/plumed.dat
+++ b/regtest/crystallization/rt-smac/plumed.dat
@@ -340,3 +340,7 @@ s2: SMAC SPECIES=m3 LOWMEM KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} KERNEL2={GAUS
 # s2num: SMAC DATA=m3 KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} KERNEL2={GAUSSIAN CENTER=pi SIGMA=0.480} SWITCH={RATIONAL R_0=0.6} MORE_THAN={RATIONAL R_0=0.7} SWITCH_COORD={EXP R_0=4} NUMERICAL_DERIVATIVES
 
 DUMPDERIVATIVES ARG=s2.* FILE=deriv FMT=%8.4f
+
+s3: SMAC LOWMEM SPECIESA=m3 SPECIESB=m1 KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} KERNEL2={GAUSSIAN CENTER=pi SIGMA=0.480} SWITCH={RATIONAL R_0=0.6} MORE_THAN={RATIONAL R_0=0.7} SWITCH_COORD={EXP R_0=4}
+
+PRINT ARG=s3.* FILE=colvar2 FMT=%8.4f
diff --git a/src/crystallization/MoleculeOrientation.cpp b/src/crystallization/MoleculeOrientation.cpp
index ab0a04c6dfbeb7c186b510abcd6e3ab4f057f1ba..3b74b69b2e16e6a088170cb366ea4b21d3c1b5e6 100644
--- a/src/crystallization/MoleculeOrientation.cpp
+++ b/src/crystallization/MoleculeOrientation.cpp
@@ -55,6 +55,7 @@ private:
 public:
   static void registerKeywords( Keywords& keys );
   explicit MoleculeOrientation( const ActionOptions& ao );
+  AtomNumber getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const ;
   void calculateVector( multicolvar::AtomValuePack& myatoms ) const;
   void normalizeVector( std::vector<double>& vals ) const ;
   void normalizeVectorDerivatives( MultiValue& myvals ) const ;
@@ -89,6 +90,12 @@ VectorMultiColvar(ao)
   } 
 }
 
+AtomNumber MoleculeOrientation::getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const {
+  plumed_dbg_assert( iatom<atom_lab.size() );
+  plumed_assert( atom_lab[iatom].first==0 );
+  return ActionAtomistic::getAbsoluteIndex( ablocks[0][atom_lab[iatom].second] );
+}
+
 void MoleculeOrientation::calculateVector( multicolvar::AtomValuePack& myatoms ) const {
   for(unsigned i=0;i<nvectors;++i){
       Vector distance; distance=getSeparation( myatoms.getPosition(2*i+0), myatoms.getPosition(2*i+1) );
diff --git a/src/crystallization/MoleculePlane.cpp b/src/crystallization/MoleculePlane.cpp
index 618b4edb43c966cd5e2a9a9645a7a86751fb9bc0..063d8e380ba33bf558a85dfa74af2c30eca669f0 100644
--- a/src/crystallization/MoleculePlane.cpp
+++ b/src/crystallization/MoleculePlane.cpp
@@ -41,6 +41,7 @@ private:
 public:
   static void registerKeywords( Keywords& keys );
   explicit MoleculePlane( const ActionOptions& ao );
+  AtomNumber getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const ;
   void calculateVector( multicolvar::AtomValuePack& myatoms ) const ;
 };
 
@@ -69,6 +70,12 @@ VectorMultiColvar(ao)
   setVectorDimensionality( 3 ); setupMultiColvarBase( all_atoms );
 }
 
+AtomNumber MoleculePlane::getAbsoluteIndexOfCentralAtom( const unsigned& iatom ) const {
+  plumed_dbg_assert( iatom<atom_lab.size() );
+  plumed_assert( atom_lab[iatom].first==0 );
+  return ActionAtomistic::getAbsoluteIndex( ablocks[0][atom_lab[iatom].second] );
+}
+
 void MoleculePlane::calculateVector( multicolvar::AtomValuePack& myatoms ) const { 
   Vector d1, d2, cp; 
   if( myatoms.getNumberOfAtoms()==3 ){
diff --git a/src/multicolvar/MultiColvarBase.cpp b/src/multicolvar/MultiColvarBase.cpp
index a4e6bd7f106fe35a5ee78fdfbb85168b018edb32..19bec3bfe1c94ed55f1b42edbd23d83ed03ed596 100644
--- a/src/multicolvar/MultiColvarBase.cpp
+++ b/src/multicolvar/MultiColvarBase.cpp
@@ -376,8 +376,8 @@ void MultiColvarBase::setupMultiColvarBase( const std::vector<AtomNumber>& atoms
             bool found=false; unsigned inum;
             for(unsigned j=0;j<nat1;++j){
                 if( atom_lab[nat1+i].first>0 && atom_lab[j].first>0 ){
-                    if( mybasemulticolvars[atom_lab[nat1+i].first-1]->getLabel()==mybasemulticolvars[atom_lab[j].first-1]->getLabel() &&
-                        atom_lab[nat1+i].second==atom_lab[j].second ){ found=true; inum=j; break; }
+                    if( mybasemulticolvars[atom_lab[nat1+i].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[nat1+i].second)==
+                        mybasemulticolvars[atom_lab[j].first-1]->getAbsoluteIndexOfCentralAtom(atom_lab[j].second) ){ found=true; inum=j; break; }
                 } else if( all_atoms[atom_lab[nat1+i].second]==all_atoms[atom_lab[j].second] ){ found=true; inum=j; break; }
             }
             // This prevents mistakes being made in colvar setup