diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0e248919..96fb4e9a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -42,6 +42,7 @@ set(SUB_DIRS
     src/NFtest/complex
     src/NFtest/templateMolecule
     src/NFtest/mappingSet
+    src/NFtest/reactantTree
     src/NFscheduler
     src/NFreactions/transformations
     src/NFreactions/reactions
diff --git a/CMakeLists.x86.txt b/CMakeLists.x86.txt
index 2c6d79a1..7231a76e 100644
--- a/CMakeLists.x86.txt
+++ b/CMakeLists.x86.txt
@@ -33,6 +33,10 @@ set(SUB_DIRS
     src/NFtest/reactionClass
     src/NFtest/observable
     src/NFtest/molecule
+    src/NFtest/mappingSet
+    src/NFtest/templateMolecule
+    src/NFtest/complex
+    src/NFtest/reactantTree
     src/NFscheduler 
     src/NFreactions/transformations 
     src/NFreactions/reactions 
diff --git a/src/NFcore/moleculeType.cpp b/src/NFcore/moleculeType.cpp
index ab063b88..9c83766d 100644
--- a/src/NFcore/moleculeType.cpp
+++ b/src/NFcore/moleculeType.cpp
@@ -238,10 +238,10 @@ void MoleculeType::addEquivalentComponents(vector <vector <string> > &identicalC
 
 
 bool MoleculeType::isIntegerComponent(const string& cName) const {
-	for(int c=0; c<numOfComponents; c++)
-			if(compName[c].compare(cName)==0) {
-				return this->isIntegerCompState[c];
-			}
+	auto it = compNameMap.find(cName);
+	if (it != compNameMap.end()) {
+		return this->isIntegerCompState[it->second];
+	}
 	cerr<<"!!! error !!! cannot find site name "<< cName << " in MoleculeType: "<<name;
 	cerr<<"in function isIntegerComponent(string cName).  "<<endl;
 	this->printDetails();
diff --git a/src/NFcore/reactionClass.cpp b/src/NFcore/reactionClass.cpp
index 6e11c9f8..d5315c75 100755
--- a/src/NFcore/reactionClass.cpp
+++ b/src/NFcore/reactionClass.cpp
@@ -406,19 +406,6 @@ string ReactionClass::fire(double random_A_number, bool track) {
 	}
 
 
-	// // output something if the reaction was tagged
-	// if(tagged) {
-	// 	for(unsigned int k=0; k<n_reactants; k++) {
-	// 		cout<<" [";
-	// 		for(unsigned int p=0; p<mappingSet[k]->getNumOfMappings();p++) {
-	// 			Molecule *mForTag = mappingSet[k]->get(p)->getMolecule();
-	// 			cout<<" "<<mForTag->getMoleculeTypeName()<<mForTag->getUniqueID();
-	// 		}
-	// 		cout<<" ]";
-	// 	}
-	// 	cout<<endl;
-	// }
-
 	// Generate the set of possible products that we need to update
 	// (excluding new molecules, we'll get those later --Justin)
 	this->transformationSet->getListOfProducts(mappingSet,products,traversalLimit);
diff --git a/src/NFcore/system.cpp b/src/NFcore/system.cpp
index 0393d936..943ffcec 100644
--- a/src/NFcore/system.cpp
+++ b/src/NFcore/system.cpp
@@ -660,12 +660,6 @@ void System::prepareForSimulation()
 	//Note!!  : the order of preparing the system matters!  You have to prepare
 	//some things before others, because certain things require other
 
-	//First, set the observables up correctly, so when functions evaluate, they get the
-	//correct values
-	//for(molTypeIter = allMoleculeTypes.begin(); molTypeIter != allMoleculeTypes.end(); molTypeIter++ ) {
-	//	(*molTypeIter)->addAllToObservables();
-	//}
-
   	//First, we have to prep all the functions...
   	for( functionIter = globalFunctions.begin(); functionIter != globalFunctions.end(); functionIter++ )
   		(*functionIter)->prepareForSimulation(this);
@@ -754,14 +748,6 @@ void System::prepareForSimulation()
 
   	//cout<<"here 7..."<<endl;
 
-  	//add all the molecules to the appropriate observables
-  	//NOT NECESSARY - molecules are added to observables when they are prepared
-  	//for(obsIter=obsToOutput.begin(); obsIter != obsToOutput.end(); obsIter++)
-  	//	(*obsIter)->clear();
-  	//for(molTypeIter = allMoleculeTypes.begin(); molTypeIter != allMoleculeTypes.end(); molTypeIter++ ) {
-  	//	(*molTypeIter)->addAllToObservables();
-  	//}
-
   	//Add the complexes to Species observables
   	int match = 0;
   	for(obsIter = speciesObservables.begin(); obsIter != speciesObservables.end(); obsIter++)
@@ -906,9 +892,6 @@ void System::update_A_tot(ReactionClass *r, double old_a, double new_a)
 {
 	a_tot = selector->update(r,old_a,new_a);
 
-	//BUILT IN DIRECT SEARCH
-	//a_tot-=old_a;
-	//a_tot+=new_a;
 }
 
 
@@ -917,16 +900,6 @@ double System::recompute_A_tot()
 	a_tot = selector->refactorPropensities();
 	return a_tot;
 
-
-//  BUILT IN DIRECT SEARCH
-//	//Loop through the reactions and add up the rates
-//	a_tot = 0;
-//	for(rxnIter = allReactions.begin(); rxnIter != allReactions.end(); rxnIter++ )
-//	{
-//		a_tot += (*rxnIter)->update_a();
-//		if(DEBUG) (*rxnIter)->printDetails();
-//	}
-//	return a_tot;
 }
 
 
diff --git a/src/NFfunction/muParser/muParserFixes.h b/src/NFfunction/muParser/muParserFixes.h
index 0e79d37f..d024df1a 100644
--- a/src/NFfunction/muParser/muParserFixes.h
+++ b/src/NFfunction/muParser/muParserFixes.h
@@ -58,117 +58,5 @@
 #endif
 
 
-//---------------------------------------------------------------------------
-//
-// MSVC6
-//
-//---------------------------------------------------------------------------
-
-
-#if defined(_MSC_VER) && _MSC_VER==1200
-
-/** \brief Macro to replace the MSVC6 auto_ptr with the _my_auto_ptr class.
-
-    Hijack auto_ptr and replace it with a version that actually does
-    what an auto_ptr normally does. If you use std::auto_ptr in your other code
-    might either explode or work much better. The original crap created
-    by Microsoft, called auto_ptr and bundled with MSVC6 is not standard compliant.
-*/
-#define auto_ptr _my_auto_ptr
-
-// This is another stupidity that needs to be undone in order to de-pollute
-// the global namespace!
-#undef min
-#undef max
-
-
-namespace std
-{
-  typedef ::size_t size_t;
-
-  //---------------------------------------------------------------------------
-  /** \brief MSVC6 fix: Put rand into namespace std. */
-  using ::rand;
-
-  //---------------------------------------------------------------------------
-  /** \brief MSVC6 fix: Put strlen into namespace std. */
-  using ::strlen;
-
-  //---------------------------------------------------------------------------
-  /** \brief MSVC6 fix: Put strncmp into namespace std. */
-  using ::strncmp;
-
-  //---------------------------------------------------------------------------
-  template<typename T>
-  T max(T a, T b)
-  {
-    return (a>b) ? a : b;
-  }
-
-  //---------------------------------------------------------------------------
-  template<typename T>
-  T min(T a, T b)
-  {
-    return (a<b) ? a : b;
-  }
-
-  //---------------------------------------------------------------------------
-  /** Standard compliant auto_ptr redefinition for MSVC6.
-
-    The code is taken from VS.NET 2003, slightly modified to reduce
-    it's dependencies from other classes.
-  */
-  template<class _Ty>
-  class _my_auto_ptr
-  {
-  public:
-    typedef _Ty element_type;
-
-	  explicit _my_auto_ptr(_Ty *_Ptr = 0)
-	    :_Myptr(_Ptr)
-	  {}
-
-	  _my_auto_ptr(_my_auto_ptr<_Ty>& _Right)
-	    :_Myptr(_Right.release())
-	  {}
-
-	  template<class _Other>
-    operator _my_auto_ptr<_Other>()
-	  {
-      return (_my_auto_ptr<_Other>(*this));
-	  }
-
-	  template<class _Other>
-	  _my_auto_ptr<_Ty>& operator=(_my_auto_ptr<_Other>& _Right)
-	  {
-      reset(_Right.release());
-	    return (*this);
-	  }
-
-	 ~auto_ptr()              { delete _Myptr;    }
-	  _Ty& operator*() const  { return (*_Myptr); }
-	  _Ty *operator->() const	{ return (&**this);	}
-	  _Ty *get() const        { return (_Myptr);	}
-
-    _Ty *release()
-    {
-		  _Ty *_Tmp = _Myptr;
-		  _Myptr = 0;
-		  return (_Tmp);
-		}
-
-	  void reset(_Ty* _Ptr = 0)
-		{
-		  if (_Ptr != _Myptr)
-			  delete _Myptr;
-		  _Myptr = _Ptr;
-    }
-
-  private:
-	    _Ty *_Myptr;
-	}; // class _my_auto_ptr
-} // namespace std
-
-#endif // Microsoft Visual Studio Version 6.0
 
 #endif // include guard
diff --git a/src/NFinput/NFinput.cpp b/src/NFinput/NFinput.cpp
index 943fdbf3..363f7d6b 100644
--- a/src/NFinput/NFinput.cpp
+++ b/src/NFinput/NFinput.cpp
@@ -667,12 +667,6 @@ bool NFinput::initMoleculeTypes(
 			firstSymSiteToAppend.clear();
 		}
 
-		// prints out allowed state map
-		//for ( std::map< string, int, std::less< int > >::const_iterator iter = allowedStates.begin();
-		//      iter != allowedStates.end(); ++iter )
-		//      cout << iter->first << '\t' << iter->second << '\n';
-
-
 		//Getting here means we read everything we could successfully
 		return true;
 	} catch (...) {
@@ -996,12 +990,6 @@ static bool processSingleSpecies(
 				usedComponentNames.clear();
 
 
-				//We dont' have to do this anymore, because we handled it earlier!
-				//int eqClassCount = mt->getNumOfEquivalencyClasses();
-				//int *currentCount = new int[eqClassCount];
-				//for(int i=0; i<eqClassCount; i++) { currentCount[i]=1; }
-				//string *eqCompNames = mt->getEquivalencyClassCompNames();
-
 				//loop to create the actual molecules of this type
 				vector <Molecule *> currentM;
 				molecules.push_back(currentM);
@@ -1016,21 +1004,11 @@ static bool processSingleSpecies(
 						mids.push_back(mol->getMoleculeType()->getTypeID());
 						mgids.push_back(mol->getUniqueID());
 
-						//for(int i=0; i<eqClassCount; i++) { currentCount[i]=1; }
-
 						//Loop through the states and set the ones we need to set
 						int k=0;
 						for(snIter = stateName.begin(); snIter != stateName.end(); k++, snIter++ )
 						{
-							//if(mt->isEquivalentComponent((*snIter))) {
-							//	int eqNum = mt->getEquivalencyClassNumber((*snIter));
-							//	std::stringstream numStream; numStream << currentCount[eqNum];
-							//	string postFix = numStream.str();
-							//	m->setComponentState((*snIter)+postFix, (int)stateValue.at(k));
-							//	currentCount[eqNum]++;
-							//} else {
-								mol->setComponentState((*snIter), (int)stateValue.at(k));
-							//}
+							mol->setComponentState((*snIter), (int)stateValue.at(k));
 							
 							// AS2023 - this is here to reduce the number of operations written
 							// note that the default molecule starts the component state at the 
@@ -1064,8 +1042,6 @@ static bool processSingleSpecies(
 
 				}
 
-				//delete [] currentCount;
-
 				//Reset the states for the next wave...
 				stateName.clear();
 				stateValue.clear();
diff --git a/src/NFinput/parseSymRxns.cpp b/src/NFinput/parseSymRxns.cpp
index 3ddbc51a..e58f9cdc 100644
--- a/src/NFinput/parseSymRxns.cpp
+++ b/src/NFinput/parseSymRxns.cpp
@@ -372,7 +372,7 @@ void createFullSymMaps(
 //saves all possible names for all possible components
 void assembleFullSymmetryList(
 		vector <vector <vector <component> > > &symmetries, //for output
-		vector <string> &moleculeIds,    //also for output
+		map <string, int> &moleculeIds,    //also for output
 		map<string,component> &symComps, //the input of symmetric components
 		bool isRxnCenter  //set to true if you are looking at reaction centers
 		)
@@ -392,13 +392,13 @@ void assembleFullSymmetryList(
 			string thisMoleculeId = id.substr(0,length);
 			int moleculeIndex = -1;
 
-			for(unsigned int i=0; i<moleculeIds.size(); i++) {
-				if(moleculeIds.at(i).compare(thisMoleculeId)==0)
-					moleculeIndex = i; break;
+			map<string, int>::iterator mIt = moleculeIds.find(thisMoleculeId);
+			if (mIt != moleculeIds.end()) {
+				moleculeIndex = mIt->second;
 			}
 			if(moleculeIndex==-1) {
 				moleculeIndex = moleculeIds.size();
-				moleculeIds.push_back(thisMoleculeId);
+				moleculeIds[thisMoleculeId] = moleculeIndex;
 
 				//Create the vector to store all of our potential permutations
 				vector <vector <component> > v;
@@ -430,7 +430,7 @@ void assembleFullSymmetryList(
 //saves all possible names for all possible components
 void assembleFullSymmetryListOnRxnCenter(
 		vector <vector <vector <component> > > &symmetries, //for output
-		vector <string> &moleculeIds,    //also for output
+		map <string, int> &moleculeIds,    //also for output
 		map<string,component> &symComps //the input of symmetric components
 		)
 {
@@ -449,14 +449,13 @@ void assembleFullSymmetryListOnRxnCenter(
 			string thisMoleculeId = id.substr(0,length);
 			int moleculeIndex = -1;
 
-			for(unsigned int i=0; i<moleculeIds.size(); i++) {
-				if(moleculeIds.at(i).compare(thisMoleculeId)==0) {
-					moleculeIndex = i; break;
-				}
+			map<string, int>::iterator mIt = moleculeIds.find(thisMoleculeId);
+			if (mIt != moleculeIds.end()) {
+				moleculeIndex = mIt->second;
 			}
 			if(moleculeIndex==-1) {
 				moleculeIndex = moleculeIds.size();
-				moleculeIds.push_back(thisMoleculeId);
+				moleculeIds[thisMoleculeId] = moleculeIndex;
 
 				//Create the vector to store all of our potential permutations
 				vector <vector <component> > v;
@@ -528,7 +527,7 @@ bool isMoleculePermuationValid(
 //
 void assembleOffRxnCenterSymClasses(
 		vector <vector <vector <string> > > &offRxnCenterSymClasses, //the output
-		vector <string> &moleculeIds, //input list of molecule names
+		map <string, int> &moleculeIds, //input list of molecule names
 		map<string,component> &symComps)  //input list of symmetric components off the rxn center
 {
 	offRxnCenterSymClasses.clear();
@@ -550,12 +549,9 @@ void assembleOffRxnCenterSymClasses(
 		string thisMoleculeId = id.substr(0,length);
 
 		int mIndex = -1;
-		for(unsigned int i=0; i<moleculeIds.size(); i++) {
-			if(thisMoleculeId.compare(moleculeIds.at(i))==0) {
-				//cout<<" found at index: "<<i<<endl;
-				mIndex=i;
-				break;
-			}
+		map<string, int>::iterator mIt = moleculeIds.find(thisMoleculeId);
+		if (mIt != moleculeIds.end()) {
+			mIndex = mIt->second;
 		}
 		if(mIndex==-1) { cout<<"ERROR in parseSymRxns.cpp - in assebmly of off rxn center sym classes"<<endl; exit(1); }
 
@@ -662,7 +658,7 @@ bool NFinput::generateRxnPermutations(vector<map<string,component> > &permutatio
 	if(verbose) cout<<"\t\t\tGenerating symmetric permutations..."<<endl;
 
 	vector <vector <vector <component> > > symmetries;
-	vector <string> moleculeIds;
+	map <string, int> moleculeIds;
 
 	//Assemble the list of possible components for each symmetric class on a reaction center
 	assembleFullSymmetryListOnRxnCenter(symmetries,moleculeIds,symRxnCenter);
diff --git a/src/NFreactions/reactions/DORreaction.cpp b/src/NFreactions/reactions/DORreaction.cpp
index 415d89e8..936e0bef 100644
--- a/src/NFreactions/reactions/DORreaction.cpp
+++ b/src/NFreactions/reactions/DORreaction.cpp
@@ -64,9 +64,6 @@ DORRxnClass::DORRxnClass(
 	//Step 2: Some bookkeeping so that we can quickly get the function values from a mapping set
 	// Now that we have found the DOR reactant, which can potentially have multiple functions, lets
 	// figure out which functions apply to which
-	// vector <int> indexIntoMappingSet;  //list of the index into the transformations for each of the local functions
-	//vector <double> localFunctionValue;  //list of the value of each of the local functions needed to evaluate
-	                                    //the rate law
 	//Array to double check that we have used all pointer references we have created
 	bool *hasMatched = new bool [transformationSet->getNumOfTransformations(DORreactantIndex)];
 	for(int i=0; i<transformationSet->getNumOfTransformations(DORreactantIndex); i++) hasMatched[i]=false;
@@ -103,10 +100,6 @@ DORRxnClass::DORRxnClass(
 					argMappedMolecule[i] = 0;
 					argScope[i] = lfr->getFunctionScope();
 
-
-					//this->lfList.push_back(lfList.at(i));
-					//localFunctionValue.push_back(0);
-					//indexIntoMappingSet.push_back(k);
 					hasMatched[k]=true;
 					match=true;
 				}
@@ -520,15 +513,7 @@ double DORRxnClass::evaluateLocalFunctions(MappingSet *ms)
 
 	return value;
 
-	/*Molecule
 
-	for(int i=0; i<(signed)lfList.size(); i++) {
-		Molecule *molObject = ms->get(this->indexIntoMappingSet.at(i))->getMolecule();
-		int index = lfList.at(i)->getIndexOfTypeIFunctionValue(molObject);
-		this->localFunctionValue.at(i)=molObject->getLocalFunctionValue(index);
-	}
-	return this->localFunctionValue.at(0);
-	*/
 }
 
 
@@ -846,8 +831,6 @@ DOR2RxnClass::DOR2RxnClass(
 	//Step 2: Some bookkeeping so that we can quickly get the function values from a mapping set
 	// Now that we have found the DOR reactant, which can potentially have multiple functions, lets
 	// figure out which functions apply to which
-	// vector <int> indexIntoMappingSet;    //list of the index into the transformations for each of the local functions
-	// vector <double> localFunctionValue;  //list of the value of each of the local functions needed to evaluate the rate law
 
 	// DOR reactant1
 	//Array to double check that we have used all pointer references we have created
diff --git a/src/NFreactions/reactions/reaction.hh b/src/NFreactions/reactions/reaction.hh
index 13664e71..d9962bc4 100644
--- a/src/NFreactions/reactions/reaction.hh
+++ b/src/NFreactions/reactions/reaction.hh
@@ -153,15 +153,6 @@ namespace NFcore
 			Molecule ** argMappedMolecule;
 			int * argScope;
 
-
-			//vector <int> argIndexIntoMappingSet;
-
-
-
-			//vector <LocalFunction *> lfList;
-			//vector <int> indexIntoMappingSet;
-			//vector <double> localFunctionValue;
-
 	};
 
 	/* A reaction class with DOR calculations on two reactants.
diff --git a/src/NFreactions/transformations/transformationSet.cpp b/src/NFreactions/transformations/transformationSet.cpp
index 22df5684..27356cba 100644
--- a/src/NFreactions/transformations/transformationSet.cpp
+++ b/src/NFreactions/transformations/transformationSet.cpp
@@ -855,7 +855,6 @@ Molecule * TransformationSet::getPopulationPointer( unsigned int r ) const
 }
 
 bool TransformationSet::getListOfAddedMolecules(MappingSet **mappingSets, list <Molecule *> &products, int traversalLimit)
-// bool TransformationSet::getListOfAddedMolecules(MappingSet **mappingSets, vector <Molecule *> &products, int traversalLimit)
 {
 	std::unordered_set<Molecule*> product_set(products.begin(), products.end());
 
diff --git a/src/NFscheduler/Scheduler.cpp b/src/NFscheduler/Scheduler.cpp
index 243cbee4..0f2533fd 100644
--- a/src/NFscheduler/Scheduler.cpp
+++ b/src/NFscheduler/Scheduler.cpp
@@ -408,16 +408,6 @@ void slave_work(int rank, job& jnow) {
 
 	NFstream& strm = s->getOutputFileStream();
 	push_stream(rank, strm);
-	
-// 	s->prepareForSimulation();
-// 	s->updateSystemWithNewparameters();	
-
-// 	double eqTime = 0;
-// 	double sTime  = 10;
-// 	int oSteps = 10;
-	
-// 	s->equilibrate(eqTime);
-// 	s->sim(sTime, oSteps);
 }
 
 void master_init(int size) {
diff --git a/src/NFsim.cpp b/src/NFsim.cpp
index 24e7c842..546e1fa9 100644
--- a/src/NFsim.cpp
+++ b/src/NFsim.cpp
@@ -175,6 +175,7 @@
 #include "NFtest/compartment/test_compartment.hh"
 #include "NFtest/input/test_input.hh"
 #include "NFtest/mappingSet/mappingSet_test.hh"
+#include "NFtest/reactantTree/reactantTree_test.hh"
 
 #include <iostream>
 #include <string>
@@ -394,6 +395,10 @@ int runNFsimMain(int argc, char *argv[])
 					NFtest_compartment::run();
 					foundATest=true;
 				}
+				if(test=="reactantTree") {
+					NFtest_reactantTree::run();
+					foundATest=true;
+				}
 				if(test=="mappingSet") {
 					NFtest_mappingSet::run();
 					foundATest=true;
diff --git a/src/NFtest/compartment/test_compartment.cpp b/src/NFtest/compartment/test_compartment.cpp
index 5d79229a..f4842a99 100644
--- a/src/NFtest/compartment/test_compartment.cpp
+++ b/src/NFtest/compartment/test_compartment.cpp
@@ -21,22 +21,32 @@ void NFtest_compartment::run()
 
 	cout << "  Testing Compartment::isInside..." << endl;
 
+    // Test false return paths
     if (root->isInside(nullptr) != false) {
         throw std::runtime_error("isInside(nullptr) did not return false");
     }
 
+    // Test early return for identity
     if (!root->isInside(root)) {
         throw std::runtime_error("isInside(this) did not return true");
     }
 
+    // Test early return for identity on a non-root compartment
+    if (!grandchild1->isInside(grandchild1)) {
+        throw std::runtime_error("isInside(this) did not return true");
+    }
+
+    // Test pointer traversal logic (is inside parent)
     if (!grandchild1->isInside(child1)) {
         throw std::runtime_error("isInside(parent) did not return true");
     }
 
+    // Test pointer traversal logic (is inside grandparent)
     if (!grandchild1->isInside(root)) {
         throw std::runtime_error("isInside(grandparent) did not return true");
     }
 
+    // Test false return paths: passing a child to check if parent is inside
     if (child1->isInside(grandchild1)) {
         throw std::runtime_error("parent isInside(child) returned true, expected false");
     }
@@ -45,6 +55,7 @@ void NFtest_compartment::run()
         throw std::runtime_error("grandparent isInside(grandchild) returned true, expected false");
     }
 
+    // Test false return paths: checking siblings
     if (child1->isInside(child2)) {
         throw std::runtime_error("sibling isInside(sibling) returned true, expected false");
     }
diff --git a/src/NFtest/mappingSet/mappingSet_test.cpp b/src/NFtest/mappingSet/mappingSet_test.cpp
index c6d7f6cc..115631e1 100644
--- a/src/NFtest/mappingSet/mappingSet_test.cpp
+++ b/src/NFtest/mappingSet/mappingSet_test.cpp
@@ -45,17 +45,42 @@ void NFtest_mappingSet::run()
     MoleculeType* mt = new MoleculeType("TestMol", compNames, defaultStates, possibleStates, sys);
 
     // We need to pass mt, listId=0, and compartment=NULL (or appropriate compartment)
-    Molecule* mol = new Molecule(mt, 0, NULL);
+    Molecule* mol1 = new Molecule(mt, 0, NULL);
+    Molecule* mol2 = new Molecule(mt, 1, NULL);
+    Molecule* mol3 = new Molecule(mt, 2, NULL);
 
     // set molecule
-    ms->set(0, mol);
-    ms->set(1, mol);
+    ms->set(0, mol1);
+    ms->set(1, mol2);
 
-    if (ms->get(0)->getMolecule() != mol) {
+    if (ms->get(0)->getMolecule() != mol1) {
         cerr << "Failed mapping set to molecule" << endl;
         failCount++;
     }
 
+    // Test MappingSet::checkForCollisions
+    MappingSet *ms2 = new MappingSet(3, transformations);
+
+    // Setup for NO collision test
+    ms->set(0, mol1);
+    ms->set(1, mol2);
+    ms2->set(0, mol3);
+    ms2->set(1, mol3);
+
+    if (MappingSet::checkForCollisions(ms, ms2)) {
+        cerr << "Failed MappingSet::checkForCollisions: incorrectly detected a collision when there was none" << endl;
+        failCount++;
+    }
+
+    // Setup for collision test (both have mol2)
+    ms2->set(0, mol3);
+    ms2->set(1, mol2);
+
+    if (!MappingSet::checkForCollisions(ms, ms2)) {
+        cerr << "Failed MappingSet::checkForCollisions: failed to detect a collision" << endl;
+        failCount++;
+    }
+
     // Call clear
     ms->clear();
 
@@ -73,9 +98,12 @@ void NFtest_mappingSet::run()
         failCount++;
     }
 
-    delete mol;
+    delete mol1;
+    delete mol2;
+    delete mol3;
     delete sys; // Deletes molType too
     delete ms;
+    delete ms2;
     delete msClone;
 
     // Clean up transformations
diff --git a/src/NFtest/moleculeType/test_moleculeType.cpp b/src/NFtest/moleculeType/test_moleculeType.cpp
index e14995fe..8eea9283 100644
--- a/src/NFtest/moleculeType/test_moleculeType.cpp
+++ b/src/NFtest/moleculeType/test_moleculeType.cpp
@@ -115,6 +115,75 @@ void NFtest_moleculeType::run()
 
 	cout << "  MoleculeType::printDetails tests passed!" << endl;
 
+	cout << "  Testing MoleculeType::addEquivalentComponents..." << endl;
+
+	vector<string> compNames3;
+	compNames3.push_back("site1");
+	compNames3.push_back("site2");
+	compNames3.push_back("site3");
+	compNames3.push_back("otherSite");
+
+	vector<string> defaultStates3;
+	defaultStates3.push_back("u");
+	defaultStates3.push_back("u");
+	defaultStates3.push_back("u");
+	defaultStates3.push_back("u");
+
+	vector<vector<string>> allowedStates3;
+	vector<string> compAllowedStates3;
+	compAllowedStates3.push_back("u");
+	compAllowedStates3.push_back("p");
+	allowedStates3.push_back(compAllowedStates3);
+	allowedStates3.push_back(compAllowedStates3);
+	allowedStates3.push_back(compAllowedStates3);
+	allowedStates3.push_back(compAllowedStates3);
+
+	MoleculeType* mt3 = new MoleculeType("testMT_eq", compNames3, defaultStates3, allowedStates3, s);
+
+	vector<vector<string>> identicalComponents;
+	vector<string> eqGroup;
+	eqGroup.push_back("site1");
+	eqGroup.push_back("site2");
+	eqGroup.push_back("site3");
+	identicalComponents.push_back(eqGroup);
+
+	mt3->addEquivalentComponents(identicalComponents);
+
+	if (mt3->getNumOfEquivalencyClasses() != 1) {
+		throw runtime_error("addEquivalentComponents did not set the correct number of equivalency classes. Expected 1, got " + to_string(mt3->getNumOfEquivalencyClasses()));
+	}
+
+	if (mt3->getEquivalencyClassCompNames()[0] != "site") {
+		throw runtime_error("addEquivalentComponents did not set the correct generic component name. Expected 'site', got '" + mt3->getEquivalencyClassCompNames()[0] + "'");
+	}
+
+	if (mt3->getEquivalencyClassNumber("site") != 0) {
+		throw runtime_error("getEquivalencyClassNumber('site') returned " + to_string(mt3->getEquivalencyClassNumber("site")) + " instead of 0");
+	}
+
+	if (mt3->getEquivalenceClassNumber(0) != 0 ||
+	    mt3->getEquivalenceClassNumber(1) != 0 ||
+	    mt3->getEquivalenceClassNumber(2) != 0) {
+		throw runtime_error("getEquivalenceClassNumber did not map site1, site2, site3 to class 0 properly.");
+	}
+
+	if (mt3->getEquivalenceClassNumber(3) != -1) {
+		throw runtime_error("getEquivalenceClassNumber did not map otherSite to class -1 properly, got " + to_string(mt3->getEquivalenceClassNumber(3)));
+	}
+
+	int* components;
+	int n_components;
+	mt3->getEquivalencyClass(components, n_components, "site");
+
+	if (n_components != 3) {
+		throw runtime_error("getEquivalencyClass returned " + to_string(n_components) + " components for 'site' instead of 3");
+	}
+
+	if (components[0] != 0 || components[1] != 1 || components[2] != 2) {
+		throw runtime_error("getEquivalencyClass did not return the correct component indices for 'site'");
+	}
+
+	cout << "  MoleculeType::addEquivalentComponents tests passed!" << endl;
 
 	cout << "NFcore::MoleculeType tests completed successfully." << endl;
 
diff --git a/src/NFtest/reactantTree/CMakeLists.txt b/src/NFtest/reactantTree/CMakeLists.txt
new file mode 100644
index 00000000..dad500a5
--- /dev/null
+++ b/src/NFtest/reactantTree/CMakeLists.txt
@@ -0,0 +1 @@
+include_directories( ${CMAKE_SOURCE_DIR}/src )
diff --git a/src/NFtest/reactantTree/reactantTree_test.cpp b/src/NFtest/reactantTree/reactantTree_test.cpp
new file mode 100644
index 00000000..fc7c33e3
--- /dev/null
+++ b/src/NFtest/reactantTree/reactantTree_test.cpp
@@ -0,0 +1,78 @@
+#include "reactantTree_test.hh"
+#include "../../NFreactions/reactantLists/reactantTree.hh"
+#include "../../NFreactions/transformations/transformationSet.hh"
+#include <iostream>
+#include <vector>
+
+#ifndef _WIN32
+#include <unistd.h>
+#include <sys/wait.h>
+#endif
+
+using namespace std;
+using namespace NFcore;
+
+void NFtest_reactantTree::run()
+{
+    cout << "Running ReactantTree tests..." << endl;
+
+    int failCount = 0;
+
+    // Test for exit(1) on removeMappingSet from empty ReactantTree
+    cout << "  Testing removeMappingSet on empty tree (expecting exit(1))..." << endl;
+
+#ifndef _WIN32
+    // We will fork a process to test that it calls exit(1)
+    pid_t pid = fork();
+    if (pid == 0) {
+        // In the child process
+        // Redirect cerr so we don't spam the console if not needed, but here it's expected
+        if (freopen("/dev/null", "w", stderr) == nullptr) {
+            // Ignore if freopen fails
+        }
+
+        // Create an empty transformation set
+        vector<TemplateMolecule*> tempMols;
+        TransformationSet ts(tempMols);
+        ts.finalize(); // Finalize to prevent "TransformationSet cannot generate blank mapping if it is not finalized!"
+
+        // Create ReactantTree
+        ReactantTree* tree = new ReactantTree(0, &ts, 10);
+
+        // Call the method that should exit
+        tree->removeMappingSet(123);
+
+        // If we get here, the test failed (exit was not called)
+        exit(0); // Return 0 to indicate failure of the test
+    } else if (pid > 0) {
+        // In the parent process
+        int status;
+        waitpid(pid, &status, 0);
+
+        if (WIFEXITED(status)) {
+            int exit_status = WEXITSTATUS(status);
+            if (exit_status == 1) {
+                cout << "    Success: Empty tree removeMappingSet exited with code 1." << endl;
+            } else {
+                cout << "    Failure: Child process exited with code " << exit_status << " instead of 1." << endl;
+                failCount++;
+            }
+        } else {
+            cout << "    Failure: Child process did not exit normally." << endl;
+            failCount++;
+        }
+    } else {
+        cerr << "Fork failed!" << endl;
+        failCount++;
+    }
+#else
+    cout << "    Skipping exit(1) test on Windows as fork() is not available." << endl;
+#endif
+
+    if (failCount == 0) {
+        cout << "All ReactantTree tests passed successfully!" << endl;
+    } else {
+        cout << "ReactantTree tests failed with " << failCount << " errors." << endl;
+        exit(1);
+    }
+}
diff --git a/src/NFtest/reactantTree/reactantTree_test.hh b/src/NFtest/reactantTree/reactantTree_test.hh
new file mode 100644
index 00000000..10c77770
--- /dev/null
+++ b/src/NFtest/reactantTree/reactantTree_test.hh
@@ -0,0 +1,11 @@
+#ifndef REACTANTTREE_TEST_HH_
+#define REACTANTTREE_TEST_HH_
+
+#include "../../NFcore/NFcore.hh"
+
+namespace NFtest_reactantTree
+{
+    void run();
+}
+
+#endif
diff --git a/src/NFtest/transformations/test_transformations.cpp b/src/NFtest/transformations/test_transformations.cpp
index 2d7b7fa7..9851d275 100644
--- a/src/NFtest/transformations/test_transformations.cpp
+++ b/src/NFtest/transformations/test_transformations.cpp
@@ -8,10 +8,20 @@
 #include <stdexcept>
 #include <vector>
 #include <string>
+#include <sstream>
 
 using namespace std;
 using namespace NFcore;
 
+class TestTransformationSet : public TransformationSet {
+public:
+	TestTransformationSet(vector<TemplateMolecule*> &reactants) : TransformationSet(reactants) {}
+	bool testCanReach(Molecule *m1, Molecule *m2, int excludeComp) {
+		return canReachExcludingBond(m1, m2, excludeComp);
+	}
+};
+
+
 void NFtest_transformations::run()
 {
 	cout << "Running transformations tests..." << endl;
@@ -105,6 +115,22 @@ void NFtest_transformations::run()
 		throw runtime_error("TransformationSet getNmappingSets failed");
 	}
 
+	{
+		TemplateMolecule *tMissing = new TemplateMolecule(molX);
+		std::ostringstream localCerr;
+		std::streambuf* oldCerr = std::cerr.rdbuf(localCerr.rdbuf());
+		bool result = ts->addStateChangeTransform(tMissing, "p", "P");
+		std::cerr.rdbuf(oldCerr);
+
+		if (result) {
+			throw runtime_error("TransformationSet addStateChangeTransform should have failed for missing template");
+		}
+		if (localCerr.str().find("Couldn't find the template you gave me") == string::npos) {
+			throw runtime_error("TransformationSet addStateChangeTransform did not output expected error message");
+		}
+		delete tMissing;
+	}
+
 	ts->addStateChangeTransform(tx, "p", "P");
 	if (ts->getNumOfTransformations(0) != 1) {
 	    throw runtime_error("TransformationSet getNumOfTransformations failed for state change");
@@ -169,6 +195,76 @@ void NFtest_transformations::run()
 	ts2->finalize();
 	delete ts2;
 
+
+
+
+
+
+	// --- Testing canReachExcludingBond ---
+	cout << "  Testing canReachExcludingBond..." << endl;
+	vector<string> ringComps;
+	ringComps.push_back("s1");
+	ringComps.push_back("s2");
+	ringComps.push_back("s3");
+	vector<string> ringStates;
+	ringStates.push_back("No State");
+	ringStates.push_back("No State");
+	ringStates.push_back("No State");
+	vector<vector<string> > ringAllowedStates(3); // 3 components, empty lists means no states
+	vector<string> noStates;
+	ringAllowedStates[0] = noStates;
+	ringAllowedStates[1] = noStates;
+	ringAllowedStates[2] = noStates;
+	MoleculeType *molRing = new MoleculeType("Ring", ringComps, ringStates, ringAllowedStates, s);
+	s->addMoleculeType(molRing);
+
+	Molecule *m1 = molRing->genDefaultMolecule();
+	Molecule *m2 = molRing->genDefaultMolecule();
+	Molecule *m3 = molRing->genDefaultMolecule();
+	Molecule *m4 = molRing->genDefaultMolecule();
+
+	// Topology:
+	// m1(s1) - m2(s1)
+	// m2(s2) - m3(s1)
+	// m3(s2) - m4(s1)
+
+	Molecule::bind(m1, 0, m2, 0);
+	Molecule::bind(m2, 1, m3, 0);
+	Molecule::bind(m3, 1, m4, 0);
+
+	TemplateMolecule *tm1 = new TemplateMolecule(molRing);
+	vector<TemplateMolecule*> ringReactants;
+	ringReactants.push_back(tm1);
+	TestTransformationSet *testTS = new TestTransformationSet(ringReactants);
+
+	// Test on the line m1 - m2. Exclude bond at m1's s1 (index 0).
+	if (testTS->testCanReach(m1, m2, 0) != false) {
+		throw runtime_error("canReachExcludingBond failed on line topology (should be false)");
+	}
+
+	// Close the ring: m4(s2) - m1(s2)
+	Molecule::bind(m4, 1, m1, 1);
+
+	// Now m1, m2, m3, m4 are in a ring. Test excluding bond at m1's s1 (index 0).
+	if (testTS->testCanReach(m1, m2, 0) != true) {
+		throw runtime_error("canReachExcludingBond failed on ring topology (should be true)");
+	}
+
+	// Add another branch to test BFS robustness
+	Molecule *m5 = molRing->genDefaultMolecule();
+	Molecule::bind(m3, 2, m5, 0); // m3(s3) - m5(s1)
+
+	if (testTS->testCanReach(m1, m2, 0) != true) {
+		throw runtime_error("canReachExcludingBond failed on ring topology with branch (should be true)");
+	}
+
+	delete testTS;
+	delete tm1;
+
+	cout << "  canReachExcludingBond tests passed!" << endl;
+
+
+
 	cout << "  TransformationSet basic tests passed!" << endl;
 
 	cout << "  Testing SpeciesCreator..." << endl;
@@ -221,5 +317,33 @@ void NFtest_transformations::run()
 
 	cout << "  SpeciesCreator tests passed!" << endl;
 
+    // Test addExcludeReactant
+    TemplateMolecule *filterPattern = new TemplateMolecule(molX);
+    map<string, TemplateMolecule*> parsedTemplates;
+
+    // We can reuse ts3 or create a new one. Let's create a new one.
+    TemplateMolecule *tx4 = new TemplateMolecule(molX);
+    vector<TemplateMolecule*> reactants4;
+    reactants4.push_back(tx4);
+    TransformationSet *ts4 = new TransformationSet(reactants4);
+
+    // The filter expects it to match to return false.
+    ts4->addExcludeReactant(0, filterPattern, parsedTemplates);
+
+    Molecule *molX_test = molX->genDefaultMolecule();
+    bool checkFilter = ts4->checkReactantFilters(0, molX_test);
+    if (checkFilter) {
+        throw runtime_error("checkReactantFilters failed to exclude molecule matching pattern");
+    }
+
+    // Test that a filter on a different reactant index doesn't exclude it
+    bool checkFilterDiffIndex = ts4->checkReactantFilters(1, molX_test);
+    if (!checkFilterDiffIndex) {
+        throw runtime_error("checkReactantFilters excluded molecule when index didn't match");
+    }
+
+    delete ts4;
+    cout << "  TransformationSet::addExcludeReactant tests passed!" << endl;
+
 	cout << "Transformations tests completed successfully." << endl;
 }