diff --git a/clustering/coarsegrain_dependency_network.py b/clustering/coarsegrain_dependency_network.py
index a8d4ca3..9d118c4 100644
--- a/clustering/coarsegrain_dependency_network.py
+++ b/clustering/coarsegrain_dependency_network.py
@@ -89,9 +89,9 @@ def coarsegrain_dependency_network(input_folder, out_mapping_fname, num_gps):
         for i, mapp in enumerate(mapping):
             outfile_map.write(str(i) + " " + str(mapp) + "\n");
 
-    print("Converged in", iterations, "iterations")
-    print("Number of gauss points to be udpated: ", num_gp_tbu)
-    print("Number of simulations required: ", num_gp_tbu-neighbour_removed)
+    print("              Converged in", iterations, "iterations")
+    print("              Number of gauss points to be udpated: ", num_gp_tbu)
+    print("              Number of simulations required: ", num_gp_tbu-neighbour_removed)
 
 
 if __name__ == "__main__":
diff --git a/headers/FE_problem.h b/headers/FE_problem.h
index 3b07edb..5fcb097 100644
--- a/headers/FE_problem.h
+++ b/headers/FE_problem.h
@@ -1093,9 +1093,9 @@ namespace HMM
 									local_quadrature_points_history[q].new_strain[0][1],
 									local_quadrature_points_history[q].new_strain[0][2],
 									local_quadrature_points_history[q].new_strain[1][2]);
-            // remember the last ID used to get results from
-						local_quadrature_points_history[q].hist_strain.set_most_recent_ID_to_get_results_from(local_quadrature_points_history[q].hist_strain.get_ID_to_update_from());
-            // Default to get results from self
+						// remember the last ID used to get results from (WARNING: what about first timestep?)
+						local_quadrature_points_history[q].hist_strain.set_most_recent_ID_to_get_results_from(local_quadrature_points_history[q].hist_strain.get_ID_to_get_results_from());
+						// Default to get results from self
 						local_quadrature_points_history[q].hist_strain.set_ID_to_get_results_from(local_quadrature_points_history[q].qpid);
 					}
 				}
@@ -1273,7 +1273,7 @@ namespace HMM
 	
                 // might be interesting to set them global and not reload them at every iteration?		
 		num_spline_points = input_config.get<int>("model precision.clustering.spline points");
-    min_num_steps_before_spline = input_config.get<int>("model precision.clustering.min steps");
+		min_num_steps_before_spline = input_config.get<int>("model precision.clustering.min steps");
 		acceptable_diff_threshold = input_config.get<double>("model precision.clustering.diff threshold");
 		clusteringscriptsloc = input_config.get<std::string>("model precision.clustering.scripts directory");
 
@@ -1339,7 +1339,7 @@ namespace HMM
 								qp.update_strain[i] = rot_avg_upd_strain_tensor.access_raw_entry(i); 
 							}
 							qp.id = local_quadrature_points_history[q].qpid;
-              qp.most_recent_id = local_quadrature_points_history[q].hist_strain.get_most_recent_ID_to_update_from();
+							qp.most_recent_id = local_quadrature_points_history[q].hist_strain.get_most_recent_ID_to_get_results_from();
 							qp.material = celldata.get_composition(cell->active_cell_index());
 							scale_bridging_data.update_list.push_back(qp);
 							//sprintf(filename, "%s/last.%s.upstrain", macrostatelocout.c_str(), cell_id);
@@ -1516,8 +1516,7 @@ namespace HMM
 
 				for (unsigned int q=0; q<quadrature_formula.size(); ++q)
 				{
-					char cell_id[1024]; sprintf(cell_id, "%d", local_quadrature_points_history[q].hist_strain.get_ID_to_update_from());
-					int qp_id = local_quadrature_points_history[q].hist_strain.get_ID_to_update_from();
+					int qp_id = local_quadrature_points_history[q].hist_strain.get_ID_to_get_results_from();
 				
 					if (newtonstep == 0) local_quadrature_points_history[q].inc_stress = 0.;
 
@@ -1817,7 +1816,7 @@ namespace HMM
            	 	        }
 
                     std::ofstream ofile;
-                    char fname[1024]; sprintf(fname, "%s/resultforce.csv", macrologloc.c_str());
+                    char fname[1024]; sprintf(fname, "%s/loadedbc_force.csv", macrologloc.c_str());
 
                     // writing the header of the file
                     if (timestep == start_timestep){
diff --git a/headers/scale_bridging_data.h b/headers/scale_bridging_data.h
index 756ca0a..a673bba 100644
--- a/headers/scale_bridging_data.h
+++ b/headers/scale_bridging_data.h
@@ -12,7 +12,7 @@ namespace HMM {
 	struct QP // quadradture point
 	{
 		int 	id;
-    int   most_recent_id;
+		int   	most_recent_id;
 		int 	material;
 		double 	update_strain[6];
 		double	update_stress[6];
diff --git a/headers/stmd_problem.h b/headers/stmd_problem.h
index f2247f5..11d592e 100644
--- a/headers/stmd_problem.h
+++ b/headers/stmd_problem.h
@@ -30,346 +30,362 @@
 
 namespace HMM
 {
-	using namespace dealii;
-	using namespace LAMMPS_NS;
+using namespace dealii;
+using namespace LAMMPS_NS;
 
 
-	template <int dim>
-	class STMDProblem
-	{
-	public:
-		STMDProblem (MPI_Comm mdcomm, int pcolorr);
-		~STMDProblem ();
-		void strain (MDSim<dim>& md_sim, bool approx_md_with_hookes_law);
+template <int dim>
+class STMDProblem
+{
+public:
+	STMDProblem (MPI_Comm mdcomm, int pcolorr);
+	~STMDProblem ();
+	void strain (MDSim<dim>& md_sim, bool approx_md_with_hookes_law);
 
-	private:
+private:
 
-		SymmetricTensor<2,dim> lammps_straining(MDSim<dim> md_sim);
-  		SymmetricTensor<2,dim> stress_from_hookes_law (SymmetricTensor<2,dim> strain, SymmetricTensor<4,dim> stiffness);
+	SymmetricTensor<2,dim> lammps_straining(MDSim<dim> md_sim);
+	SymmetricTensor<2,dim> stress_from_hookes_law (SymmetricTensor<2,dim> strain, SymmetricTensor<4,dim> stiffness);
 
-		MPI_Comm 							md_batch_communicator;
-		const int 							md_batch_n_processes;
-		const int 							this_md_batch_process;
-		int 								md_batch_pcolor;
+	MPI_Comm 							md_batch_communicator;
+	const int 							md_batch_n_processes;
+	const int 							this_md_batch_process;
+	int 								md_batch_pcolor;
 
-		ConditionalOStream 					mdcout;
+	ConditionalOStream 					mdcout;
 
-	};
+};
 
 
-	template <int dim>
-	STMDProblem<dim>::STMDProblem (MPI_Comm mdcomm, int pcolor)
-	:
-		md_batch_communicator (mdcomm),
-		md_batch_n_processes (Utilities::MPI::n_mpi_processes(md_batch_communicator)),
-		this_md_batch_process (Utilities::MPI::this_mpi_process(md_batch_communicator)),
-		md_batch_pcolor (pcolor),
-		mdcout (std::cout,(this_md_batch_process == 0))
-	{}
+template <int dim>
+STMDProblem<dim>::STMDProblem (MPI_Comm mdcomm, int pcolor)
+:
+md_batch_communicator (mdcomm),
+md_batch_n_processes (Utilities::MPI::n_mpi_processes(md_batch_communicator)),
+this_md_batch_process (Utilities::MPI::this_mpi_process(md_batch_communicator)),
+md_batch_pcolor (pcolor),
+mdcout (std::cout,(this_md_batch_process == 0))
+{}
 
 
-	template <int dim>
-	STMDProblem<dim>::~STMDProblem ()
-	{}
+template <int dim>
+STMDProblem<dim>::~STMDProblem ()
+{}
 
 
-	// The straining function is ran on every quadrature point which
-	// requires a stress_update. Since a quandrature point is only reached*
-	// by a subset of processes N, we should automatically see lammps be
-	// parallelized on the N processes.
-	template <int dim>
-	SymmetricTensor<2,dim> STMDProblem<dim>::lammps_straining (MDSim<dim> md_sim)
-	{
-		char locff[1024]; /*reaxff*/
-		if (md_sim.force_field == "reax"){
-			sprintf(locff, "%s/ffield.reax.2", md_sim.scripts_folder.c_str()); /*reaxff*/
+// The straining function is ran on every quadrature point which
+// requires a stress_update. Since a quandrature point is only reached*
+// by a subset of processes N, we should automatically see lammps be
+// parallelized on the N processes.
+template <int dim>
+SymmetricTensor<2,dim> STMDProblem<dim>::lammps_straining (MDSim<dim> md_sim)
+{
+	char locff[1024]; /*reaxff*/
+	if (md_sim.force_field == "reax"){
+		sprintf(locff, "%s/ffield.reax.2", md_sim.scripts_folder.c_str()); /*reaxff*/
+	}
+	// Name of nanostate binary files
+	char mdstate[1024];
+	sprintf(mdstate, "%s_%d", md_sim.matid.c_str(), md_sim.replica);
+
+	char initdata[1024];
+	sprintf(initdata, "%s/init.%s.bin", md_sim.output_folder.c_str(), mdstate);
+
+	char homogdata_time[1024];
+	sprintf(homogdata_time, "%s/%s.%d.%s.lammpstrj", md_sim.log_file.c_str(),
+			md_sim.time_id.c_str(), md_sim.qp_id, mdstate);
+
+	char straindata_lcts[1024];
+	sprintf(straindata_lcts, "%s/lcts.%d.%s.dump", md_sim.restart_folder.c_str(),
+			md_sim.qp_id, mdstate);
+
+	// Find relevant dump file,
+	// Check if this qp made a dump file of or it is branching from a previous qp
+	char straindata_last_load[1024];
+	char straindata_last_write[1024];
+	if (md_sim.qp_id != md_sim.most_recent_qp_id){
+		sprintf(straindata_last_load, "%s/last.%d.%s.dump", md_sim.output_folder.c_str(),
+				md_sim.most_recent_qp_id, mdstate);
+		sprintf(straindata_last_write, "%s/last.%d.%s.dump", md_sim.output_folder.c_str(),
+						md_sim.qp_id, mdstate);
+		// Check if the 'most_recent_qp_id' is still set as its default value,
+		// if it is then there should not be any dumpfile to be found which will cause
+		// the MD simulation should be performed using the initial position of the atoms,
+		// else the dump file should exist
+		std::ifstream ifile_most_recent(straindata_last_load);
+		if(md_sim.most_recent_qp_id==std::numeric_limits<uint32_t>::max()){
+			assert (ifile_most_recent.good() == false);
 		}
-		// Name of nanostate binary files
-		char mdstate[1024];
-		sprintf(mdstate, "%s_%d", md_sim.matid.c_str(), md_sim.replica);
+		else{
+			assert (ifile_most_recent.good() == true);
+			ifile_most_recent.close();
+		}
+	}
+	else {
+		sprintf(straindata_last_load, "%s/last.%d.%s.dump", md_sim.output_folder.c_str(),
+				md_sim.qp_id, mdstate);
+		sprintf(straindata_last_write, "%s", straindata_last_load);
+	}
 
-		char initdata[1024];
-		sprintf(initdata, "%s/init.%s.bin", md_sim.output_folder.c_str(), mdstate);
+	char cline[1024];
+	char cfile[1024];
 
-		char homogdata_time[1024];
-		sprintf(homogdata_time, "%s/%s.%d.%s.lammpstrj", md_sim.log_file.c_str(),
-				md_sim.time_id.c_str(), md_sim.qp_id, mdstate);
+	// Specifying the command line options for screen and log output file
+	int nargs = 5;
+	char **lmparg = new char*[nargs];
+	lmparg[0] = NULL;
+	lmparg[1] = (char *) "-screen";
+	lmparg[2] = (char *) "none";
+	lmparg[3] = (char *) "-log";
+	lmparg[4] = new char[1024];
+	sprintf(lmparg[4], "%s/log.stress_strain", md_sim.log_file.c_str());
 
-		char straindata_lcts[1024];
-		sprintf(straindata_lcts, "%s/lcts.%d.%s.dump", md_sim.restart_folder.c_str(),
-				md_sim.qp_id, mdstate);
+	// Creating LAMMPS instance
+	LAMMPS *lmp = NULL;
+	lmp = new LAMMPS(nargs,lmparg,md_batch_communicator);
+
+	// Passing location for output as variable
+	sprintf(cline, "variable mdt string %s", md_sim.matid.c_str()); lammps_command(lmp,cline);
+	sprintf(cline, "variable loco string %s", md_sim.log_file.c_str()); lammps_command(lmp,cline);
+	sprintf(cline, "variable locs string %s", md_sim.scripts_folder.c_str()); lammps_command(lmp,cline);
+
+	// Setting testing temperature
+	sprintf(cline, "variable tempt equal %f", md_sim.temperature); lammps_command(lmp,cline);
+
+	// Setting general parameters for LAMMPS independentely of what will be
+	// tested on the sample next.
+
+	sprintf(cfile, "%s/%s", md_sim.scripts_folder.c_str(), "in.set.lammps");
+	lammps_file(lmp,cfile);
 
-    // Find relevant dump file,
-    // Check if this qp made a dump file of or it is branching from a previous qp
-		char straindata_last[1024];
-    if (md_sim.qp_id != md_sim.most_recent_qp_id){
-	  	sprintf(straindata_last, "%s/last.%d.%s.dump", md_sim.output_folder.c_str(),
-		  		md_sim.most_recent_qp_id, mdstate);
-    }
-    else {
-		  sprintf(straindata_last, "%s/last.%d.%s.dump", md_sim.output_folder.c_str(),
-			  	md_sim.qp_id, mdstate);
-    }
-
-
-		char cline[1024];
-		char cfile[1024];
-
-		// Specifying the command line options for screen and log output file
-		int nargs = 5;
-		char **lmparg = new char*[nargs];
-		lmparg[0] = NULL;
-		lmparg[1] = (char *) "-screen";
-		lmparg[2] = (char *) "none";
-		lmparg[3] = (char *) "-log";
-		lmparg[4] = new char[1024];
-		sprintf(lmparg[4], "%s/log.stress_strain", md_sim.log_file.c_str());
-
-		// Creating LAMMPS instance
-		LAMMPS *lmp = NULL;
-		lmp = new LAMMPS(nargs,lmparg,md_batch_communicator);
-
-		// Passing location for output as variable
-		sprintf(cline, "variable mdt string %s", md_sim.matid.c_str()); lammps_command(lmp,cline);
-		sprintf(cline, "variable loco string %s", md_sim.log_file.c_str()); lammps_command(lmp,cline);
-		sprintf(cline, "variable locs string %s", md_sim.scripts_folder.c_str()); lammps_command(lmp,cline);
-		
-		// Setting testing temperature
-		sprintf(cline, "variable tempt equal %f", md_sim.temperature); lammps_command(lmp,cline);
-
-		// Setting general parameters for LAMMPS independentely of what will be
-		// tested on the sample next.
-
-		sprintf(cfile, "%s/%s", md_sim.scripts_folder.c_str(), "in.set.lammps");
-		lammps_file(lmp,cfile);
-
-		/*mdcout << "               "
+	/*mdcout << "               "
 				<< "(MD - " << timeid <<"."<< cellid << " - repl " << repl << ") "
 				<< "Compute current state data...       " << std::endl;*/
 
-		// Check if a previous state has already been computed specifically for
-		// this quadrature point, otherwise use the initial state (which is the
-		// last state of this quadrature point)
-		/*mdcout << "               "
+	// Check if a previous state has already been computed specifically for
+	// this quadrature point, otherwise use the initial state (which is the
+	// last state of this quadrature point)
+	/*mdcout << "               "
 				<< "(MD - " << timeid <<"."<< cellid << " - repl " << repl << ") "
 				<< "   ... from previous state data...   " << std::flush;*/
 
-		// Check the presence of a dump file to restart from
-		std::ifstream ifile(straindata_last);
-		if (ifile.good()){
-			/*mdcout << "  specifically computed." << std::endl;*/
-			ifile.close();
-
-			if (md_sim.force_field == "reax") {
-				sprintf(cline, "read_restart %s", initdata); lammps_command(lmp,cline); /*reaxff*/
-				sprintf(cline, "rerun %s dump x y z vx vy vz ix iy iz box yes scaled yes wrapped yes format native", straindata_last); /*reaxff*/
-				lammps_command(lmp,cline); /*reaxff*/
-			}
-			else if (md_sim.force_field == "opls") {
-				sprintf(cline, "read_restart %s", straindata_last); /*opls*/
-				lammps_command(lmp,cline); /*opls*/
-			}
-
-			sprintf(cline, "print 'specifically computed'"); lammps_command(lmp,cline);
-		}
-		else{
-			/*mdcout << "  initially computed." << std::endl;*/
-			sprintf(cline, "read_restart %s", initdata); 
-			lammps_command(lmp,cline);
-			sprintf(cline, "print 'initially computed'"); lammps_command(lmp,cline);
+	// Check the presence of a dump file to restart from
+	// (either from current or previous ID_to_get_results_from)
+	std::ifstream ifile(straindata_last_load);
+	if (ifile.good()){
+		/*mdcout << "  specifically computed." << std::endl;*/
+		ifile.close();
+
+		if (md_sim.force_field == "reax") {
+			sprintf(cline, "read_restart %s", initdata); lammps_command(lmp,cline); /*reaxff*/
+			sprintf(cline, "rerun %s dump x y z vx vy vz ix iy iz box yes scaled yes wrapped yes format native", straindata_last_load); /*reaxff*/
+			lammps_command(lmp,cline); /*reaxff*/
 		}
-		
-
-		// Query box dimensions
-		char vdir[1024];
-		std::vector<double> lbdim (dim);
-		sprintf(cline, "variable ll1 equal lx"); lammps_command(lmp,cline);
-		sprintf(cline, "variable ll2 equal ly"); lammps_command(lmp,cline);
-		sprintf(cline, "variable ll3 equal lz"); lammps_command(lmp,cline);
-		for(unsigned int i=0;i<dim;i++){
-			sprintf(vdir, "ll%d",i+1);
-			lbdim[i] = *((double *) lammps_extract_variable(lmp,vdir,NULL));
+		else if (md_sim.force_field == "opls") {
+			sprintf(cline, "read_restart %s", straindata_last_load); /*opls*/
+			lammps_command(lmp,cline); /*opls*/
 		}
 
-		// Correction of strain tensor with actual box dimensions
-		for (unsigned int i=0; i<dim; i++){
-			md_sim.strain[i][i] /= lbdim[i];
-			md_sim.strain[i][(i+1)%dim] /= lbdim[(i+2)%dim];
-		}
+		sprintf(cline, "print 'specifically computed'"); lammps_command(lmp,cline);
+	}
+	else{
+		/*mdcout << "  initially computed." << std::endl;*/
+		sprintf(cline, "read_restart %s", initdata);
+		lammps_command(lmp,cline);
+		sprintf(cline, "print 'initially computed'"); lammps_command(lmp,cline);
+	}
+
+
+	// Query box dimensions
+	char vdir[1024];
+	std::vector<double> lbdim (dim);
+	sprintf(cline, "variable ll1 equal lx"); lammps_command(lmp,cline);
+	sprintf(cline, "variable ll2 equal ly"); lammps_command(lmp,cline);
+	sprintf(cline, "variable ll3 equal lz"); lammps_command(lmp,cline);
+	for(unsigned int i=0;i<dim;i++){
+		sprintf(vdir, "ll%d",i+1);
+		lbdim[i] = *((double *) lammps_extract_variable(lmp,vdir,NULL));
+	}
+
+	// Correction of strain tensor with actual box dimensions
+	for (unsigned int i=0; i<dim; i++){
+		md_sim.strain[i][i] /= lbdim[i];
+		md_sim.strain[i][(i+1)%dim] /= lbdim[(i+2)%dim];
+	}
 
-		// Number of timesteps in the MD simulation, rounding to nearest 10, enforcing at least 10
- //int nts = std::max(int(std::ceil(md_sim.strain.norm()/(md_sim.timestep_length*md_sim.strain_rate)/10)*10),1);
-		int nts;
-		double strain_time = md_sim.strain.norm() / md_sim.strain_rate;
-    nts = std::ceil( (strain_time/md_sim.timestep_length) /10.0) * 10;// rounded to the nearest 10
-		nts = std::max(nts,10); // check that its not 0
+	// Number of timesteps in the MD simulation, rounding to nearest 10, enforcing at least 10
+	//int nts = std::max(int(std::ceil(md_sim.strain.norm()/(md_sim.timestep_length*md_sim.strain_rate)/10)*10),1);
+	int nts;
+	double strain_time = md_sim.strain.norm() / md_sim.strain_rate;
+	nts = std::ceil( (strain_time/md_sim.timestep_length) /10.0) * 10;// rounded to the nearest 10
+	nts = std::max(nts,10); // check that its not 0
 
 
-		sprintf(cline, "variable dts equal %f", md_sim.timestep_length); lammps_command(lmp,cline);
-		sprintf(cline, "variable nts equal %d", nts); lammps_command(lmp,cline);
+	sprintf(cline, "variable dts equal %f", md_sim.timestep_length); lammps_command(lmp,cline);
+	sprintf(cline, "variable nts equal %d", nts); lammps_command(lmp,cline);
 
-		for(unsigned int k=0;k<dim;k++)
-			for(unsigned int l=k;l<dim;l++)
-			{
-				sprintf(cline, "variable ceeps_%d%d equal %.6e", k, l, 
-												md_sim.strain[k][l]/(nts*md_sim.timestep_length));
-				lammps_command(lmp,cline);
-			}
+	for(unsigned int k=0;k<dim;k++)
+		for(unsigned int l=k;l<dim;l++)
+		{
+			sprintf(cline, "variable ceeps_%d%d equal %.6e", k, l,
+					md_sim.strain[k][l]/(nts*md_sim.timestep_length));
+			lammps_command(lmp,cline);
+		}
 
-		// Run the NEMD simulations of the strained box
-		/*mdcout << "               "
+	// Run the NEMD simulations of the strained box
+	/*mdcout << "               "
 				<< "(MD - " << timeid <<"."<< cellid << " - repl " << repl << ") "
 				<< "   ... reading and executing in.strain.lammps.       " << std::endl;*/
-		sprintf(cfile, "%s/%s", md_sim.scripts_folder.c_str(), "in.strain.lammps");
-		lammps_file(lmp,cfile);
+	sprintf(cfile, "%s/%s", md_sim.scripts_folder.c_str(), "in.strain.lammps");
+	lammps_file(lmp,cfile);
 
-		/*mdcout << "               "
+	/*mdcout << "               "
 				<< "(MD - " << timeid <<"."<< cellid << " - repl " << repl << ") "
 				<< "Saving state data...       " << std::endl;*/
-		// Save data to specific file for this quadrature point
-		if (md_sim.force_field == "opls"){
-			sprintf(cline, "write_restart %s", straindata_last); /*opls*/
-			lammps_command(lmp,cline); /*opls*/
-		}
-		else if (md_sim.force_field == "reax"){
-			sprintf(cline, "write_dump all custom %s id type xs ys zs vx vy vz ix iy iz", straindata_last); /*reaxff*/
-			lammps_command(lmp,cline); /*reaxff*/
-		}
+	// Save data to specific file for this quadrature point
+	if (md_sim.force_field == "opls"){
+		sprintf(cline, "write_restart %s", straindata_last_write); /*opls*/
+		lammps_command(lmp,cline); /*opls*/
+	}
+	else if (md_sim.force_field == "reax"){
+		sprintf(cline, "write_dump all custom %s id type xs ys zs vx vy vz ix iy iz", straindata_last_write); /*reaxff*/
+		lammps_command(lmp,cline); /*reaxff*/
+	}
 
-		if(md_sim.checkpoint){
-			if (md_sim.force_field == "opls") {
-				sprintf(cline, "write_restart %s", straindata_lcts); lammps_command(lmp,cline); /*opls*/
-			}
-			else if (md_sim.force_field == "reax") {
-				sprintf(cline, "write_dump all custom %s id type xs ys zs vx vy vz ix iy iz", straindata_lcts); lammps_command(lmp,cline); /*reaxff*/
-			}
+	if(md_sim.checkpoint){
+		if (md_sim.force_field == "opls") {
+			sprintf(cline, "write_restart %s", straindata_lcts); lammps_command(lmp,cline); /*opls*/
+		}
+		else if (md_sim.force_field == "reax") {
+			sprintf(cline, "write_dump all custom %s id type xs ys zs vx vy vz ix iy iz", straindata_lcts); lammps_command(lmp,cline); /*reaxff*/
 		}
-		// close down LAMMPS
-		delete lmp;
+	}
+	// close down LAMMPS
+	delete lmp;
 
-		/*mdcout << "               "
+	/*mdcout << "               "
 				<< "(MD - " << timeid <<"."<< cellid << " - repl " << repl << ") "
 				<< "Homogenization of stiffness and stress using in.elastic.lammps...       " << std::endl;*/
 
-		// Creating LAMMPS instance
-		sprintf(lmparg[4], "%s/log.homogenization", md_sim.log_file.c_str());
-		lmp = new LAMMPS(nargs,lmparg,md_batch_communicator);
+	// Creating LAMMPS instance
+	sprintf(lmparg[4], "%s/log.homogenization", md_sim.log_file.c_str());
+	lmp = new LAMMPS(nargs,lmparg,md_batch_communicator);
 
-        	sprintf(cline, "variable loco string %s", md_sim.log_file.c_str()); lammps_command(lmp,cline);
-		sprintf(cline, "variable locs string %s", md_sim.scripts_folder.c_str()); lammps_command(lmp,cline);
-		
-		// Setting testing temperature
-		sprintf(cline, "variable tempt equal %f", md_sim.temperature); lammps_command(lmp,cline);
+	sprintf(cline, "variable loco string %s", md_sim.log_file.c_str()); lammps_command(lmp,cline);
+	sprintf(cline, "variable locs string %s", md_sim.scripts_folder.c_str()); lammps_command(lmp,cline);
 
-		// Setting general parameters for LAMMPS independentely of what will be
-		// tested on the sample next.
-		sprintf(cfile, "%s/%s", md_sim.scripts_folder.c_str(), "in.set.lammps");
-		lammps_file(lmp,cfile);
+	// Setting testing temperature
+	sprintf(cline, "variable tempt equal %f", md_sim.temperature); lammps_command(lmp,cline);
 
-		if (md_sim.force_field == "reax"){
-			sprintf(cline, "read_restart %s", initdata); /*reaxff*/
-			lammps_command(lmp,cline); /*reaxff*/
-			sprintf(cline, "rerun %s dump x y z vx vy vz ix iy iz box yes scaled yes wrapped yes format native", straindata_last); /*reaxff*/
-			lammps_command(lmp,cline); /*reaxff*/
+	// Setting general parameters for LAMMPS independentely of what will be
+	// tested on the sample next.
+	sprintf(cfile, "%s/%s", md_sim.scripts_folder.c_str(), "in.set.lammps");
+	lammps_file(lmp,cfile);
 
-		}
-		else if (md_sim.force_field == "opls"){
-			sprintf(cline, "read_restart %s", straindata_last); /*opls*/
-			lammps_command(lmp,cline); /*opls*/
-		}
-    
-		sprintf(cline, "reset_timestep 0"); lammps_command(lmp,cline);
+	if (md_sim.force_field == "reax"){
+		sprintf(cline, "read_restart %s", initdata); /*reaxff*/
+		lammps_command(lmp,cline); /*reaxff*/
+		sprintf(cline, "rerun %s dump x y z vx vy vz ix iy iz box yes scaled yes wrapped yes format native", straindata_last_write); /*reaxff*/
+		lammps_command(lmp,cline); /*reaxff*/
 
-		sprintf(cline, "variable dts equal %f", md_sim.timestep_length); lammps_command(lmp,cline);
+	}
+	else if (md_sim.force_field == "opls"){
+		sprintf(cline, "read_restart %s", straindata_last_write); /*opls*/
+		lammps_command(lmp,cline); /*opls*/
+	}
 
-		if(md_sim.output_homog){
-			// Setting dumping of atom positions for post analysis of the MD simulation
-			// DO NOT USE CUSTOM DUMP: WRONG ATOM POSITIONS...
-			sprintf(cline, "dump atom_dump all atom %d %s", 1, homogdata_time); lammps_command(lmp,cline);
-		}
+	sprintf(cline, "reset_timestep 0"); lammps_command(lmp,cline);
 
-		// Compute the secant stiffness tensor at the given stress/strain state
-		sprintf(cline, "variable locbe string %s/%s", md_sim.scripts_folder.c_str(), "ELASTIC");
-		lammps_command(lmp,cline);
+	sprintf(cline, "variable dts equal %f", md_sim.timestep_length); lammps_command(lmp,cline);
 
-		// Set sampling and straining time-lengths
-		sprintf(cline, "variable nssample0 equal %d", md_sim.nsteps_sample); lammps_command(lmp,cline);
-		sprintf(cline, "variable nssample  equal %d", md_sim.nsteps_sample); lammps_command(lmp,cline);
-
-		// Using a routine based on the example ELASTIC/ to compute the stress tensor
-		sprintf(cfile, "%s/%s", md_sim.scripts_folder.c_str(), "ELASTIC/in.homogenization.lammps");
-		lammps_file(lmp,cfile);
-
-		// Filling 3x3 stress tensor and conversion from ATM to Pa
-		// Useless at the moment, since it cannot be used in the Newton-Raphson algorithm.
-		// The MD evaluated stress is flucutating too much (few MPa), therefore prevents
-		// the iterative algorithm to converge...
-		for(unsigned int k=0;k<dim;k++)
-			for(unsigned int l=k;l<dim;l++)
-			{
-				char vcoef[1024];
-				sprintf(vcoef, "pp%d%d", k+1, l+1);
-				md_sim.stress[k][l] = *((double *) lammps_extract_variable(lmp,vcoef,NULL))*(-1.0)*1.01325e+05;
-			}
-
-		if(md_sim.output_homog){
-			// Unetting dumping of atom positions
-			sprintf(cline, "undump atom_dump"); lammps_command(lmp,cline);
-		}
+	if(md_sim.output_homog){
+		// Setting dumping of atom positions for post analysis of the MD simulation
+		// DO NOT USE CUSTOM DUMP: WRONG ATOM POSITIONS...
+		sprintf(cline, "dump atom_dump all atom %d %s", 1, homogdata_time); lammps_command(lmp,cline);
+	}
 
-		// close down LAMMPS
-		delete lmp;
+	// Compute the secant stiffness tensor at the given stress/strain state
+	sprintf(cline, "variable locbe string %s/%s", md_sim.scripts_folder.c_str(), "ELASTIC");
+	lammps_command(lmp,cline);
 
-		return md_sim.stress;
-	}
+	// Set sampling and straining time-lengths
+	sprintf(cline, "variable nssample0 equal %d", md_sim.nsteps_sample); lammps_command(lmp,cline);
+	sprintf(cline, "variable nssample  equal %d", md_sim.nsteps_sample); lammps_command(lmp,cline);
 
+	// Using a routine based on the example ELASTIC/ to compute the stress tensor
+	sprintf(cfile, "%s/%s", md_sim.scripts_folder.c_str(), "ELASTIC/in.homogenization.lammps");
+	lammps_file(lmp,cfile);
 
-	template <int dim>
-  SymmetricTensor<2,dim> STMDProblem<dim>::stress_from_hookes_law (SymmetricTensor<2,dim> strain, SymmetricTensor<4,dim> stiffness)
-	{
-		SymmetricTensor<2,dim> stress;
-		stress = stiffness * strain;
-		return stress;
+	// Filling 3x3 stress tensor and conversion from ATM to Pa
+	// Useless at the moment, since it cannot be used in the Newton-Raphson algorithm.
+	// The MD evaluated stress is flucutating too much (few MPa), therefore prevents
+	// the iterative algorithm to converge...
+	for(unsigned int k=0;k<dim;k++)
+		for(unsigned int l=k;l<dim;l++)
+		{
+			char vcoef[1024];
+			sprintf(vcoef, "pp%d%d", k+1, l+1);
+			md_sim.stress[k][l] = *((double *) lammps_extract_variable(lmp,vcoef,NULL))*(-1.0)*1.01325e+05;
+		}
+
+	if(md_sim.output_homog){
+		// Unetting dumping of atom positions
+		sprintf(cline, "undump atom_dump"); lammps_command(lmp,cline);
 	}
 
+	// close down LAMMPS
+	delete lmp;
 
-	template <int dim>
-  void STMDProblem<dim>::strain (MDSim<dim>& md_sim, bool approx_md_with_hookes_law)
-	{
+	return md_sim.stress;
+}
 
-		if (md_sim.force_field != "opls" && md_sim.force_field != "reax"){
-			std::cerr << "Error: Force field is " << md_sim.force_field
-					  << " but only 'opls' and 'reax' are implemented... "
-					  << std::endl;
-			exit(1);
-		}
 
-		// Then the lammps function instanciates lammps, starting from an initial
-		// microstructure and applying the complete new_strain or starting from
-		// the microstructure at the old_strain and applying the difference between
-		// the new_ and _old_strains, returns the new_stress state.
-		if(this_md_batch_process == 0)
-		{
-			std::cout << " \t" << md_sim.qp_id <<"-"<< md_sim.replica<<"-start" << std::endl << std::flush;
-		}
-		MPI_Barrier(md_batch_communicator);
+template <int dim>
+SymmetricTensor<2,dim> STMDProblem<dim>::stress_from_hookes_law (SymmetricTensor<2,dim> strain, SymmetricTensor<4,dim> stiffness)
+{
+	SymmetricTensor<2,dim> stress;
+	stress = stiffness * strain;
+	return stress;
+}
 
-		if (approx_md_with_hookes_law == true){
-			// this option is meant for testing
-			md_sim.stress = stress_from_hookes_law(md_sim.strain, md_sim.stiffness);
-			md_sim.stress_updated = true;
-		}
-		else {
-			md_sim.stress = lammps_straining(md_sim);
-			md_sim.stress_updated = true;
-		}
 
-		MPI_Barrier(md_batch_communicator);
-		if(this_md_batch_process == 0)
-		{
-			std::cout << " \t" << md_sim.qp_id <<"-"<< md_sim.replica << std::endl << std::flush;
-		}
+template <int dim>
+void STMDProblem<dim>::strain (MDSim<dim>& md_sim, bool approx_md_with_hookes_law)
+{
+
+	if (md_sim.force_field != "opls" && md_sim.force_field != "reax"){
+		std::cerr << "Error: Force field is " << md_sim.force_field
+				<< " but only 'opls' and 'reax' are implemented... "
+				<< std::endl;
+		exit(1);
 	}
+
+	// Then the lammps function instanciates lammps, starting from an initial
+	// microstructure and applying the complete new_strain or starting from
+	// the microstructure at the old_strain and applying the difference between
+	// the new_ and _old_strains, returns the new_stress state.
+	if(this_md_batch_process == 0)
+	{
+		std::cout << " \t" << md_sim.qp_id <<"-"<< md_sim.replica<<"-start" << std::endl << std::flush;
+	}
+	MPI_Barrier(md_batch_communicator);
+
+	if (approx_md_with_hookes_law == true){
+		// this option is meant for testing
+		md_sim.stress = stress_from_hookes_law(md_sim.strain, md_sim.stiffness);
+		md_sim.stress_updated = true;
+	}
+	else {
+		md_sim.stress = lammps_straining(md_sim);
+		md_sim.stress_updated = true;
+	}
+
+	MPI_Barrier(md_batch_communicator);
+	if(this_md_batch_process == 0)
+	{
+		std::cout << " \t" << md_sim.qp_id <<"-"<< md_sim.replica << std::endl << std::flush;
+	}
+}
 }
 
 #endif
diff --git a/headers/stmd_sync.h b/headers/stmd_sync.h
index b8e6718..e786629 100644
--- a/headers/stmd_sync.h
+++ b/headers/stmd_sync.h
@@ -9,6 +9,7 @@
 #include <string>
 #include <sys/stat.h>
 #include <math.h>
+#include <assert.h>
 
 #include "boost/archive/text_oarchive.hpp"
 #include "boost/archive/text_iarchive.hpp"
@@ -490,7 +491,7 @@ std::vector< MDSim<dim> > STMDSync<dim>::prepare_md_simulations(ScaleBridgingDat
 
 	// Setting up batch of processes
 	set_md_procs(nmdruns);
-	mcout << "set md procs"<<std::endl;
+
 	for (uint32_t qp=0; qp<n_qp; ++qp)
 	{
 		for(uint32_t repl=0; repl<nrepl; repl++)
@@ -500,7 +501,7 @@ std::vector< MDSim<dim> > STMDSync<dim>::prepare_md_simulations(ScaleBridgingDat
 
 			MDSim<dim> md_sim;
 			md_sim.qp_id = update_list[qp].id;
-      md_sim.most_recent_qp_id = update_list[qp].most_recent_id;
+			md_sim.most_recent_qp_id = update_list[qp].most_recent_id;
 
 			md_sim.replica = repl + 1; // +1 to match input file lables... fix
 			md_sim.material = update_list[qp].material;
@@ -509,7 +510,7 @@ std::vector< MDSim<dim> > STMDSync<dim>::prepare_md_simulations(ScaleBridgingDat
 			md_sim.matid = replica_data[replica_data_index].mat;
 			md_sim.time_id = time_id;
 
-			md_sim.force_field 			= md_force_field;
+			md_sim.force_field 		= md_force_field;
 			md_sim.timestep_length  = md_timestep_length;
 			md_sim.temperature      = md_temperature;
 			md_sim.nsteps_sample    = md_nsteps_sample;
@@ -543,9 +544,9 @@ std::vector< MDSim<dim> > STMDSync<dim>::prepare_md_simulations(ScaleBridgingDat
 			md_sim.stiffness = replica_data[replica_data_index].init_stiff;
 
 			request_simulations.push_back(md_sim);
-			//}
 		}
 	}
+
 	return request_simulations;
 }
 
diff --git a/headers/strain2spline.h b/headers/strain2spline.h
index 1f7b6e0..5f67172 100644
--- a/headers/strain2spline.h
+++ b/headers/strain2spline.h
@@ -384,12 +384,12 @@ namespace MatHistPredict {
             /* Returns the ID of the gauss point whose MD simulation(s) results should be used
              * also by this gauss point.
              */
-            uint32_t get_ID_to_update_from()
+            uint32_t get_ID_to_get_results_from()
             {
                 return ID_to_get_results_from;
             }
 
-            uint32_t get_most_recent_ID_to_update_from()
+            uint32_t get_most_recent_ID_to_get_results_from()
             {
                 return most_recent_ID_to_get_results_from;
             }