diff --git a/CMakeLists/myriad.CMakeLists.txt b/CMakeLists/myriad.CMakeLists.txt
index c940126..f8bd48d 100644
--- a/CMakeLists/myriad.CMakeLists.txt
+++ b/CMakeLists/myriad.CMakeLists.txt
@@ -1,7 +1,7 @@
 CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
 
 FIND_PACKAGE(deal.II 8.4.1 QUIET
-  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  HINTS /home/uccamva/opt/dealii-9.0.0/ ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
   )
 IF(NOT ${deal.II_FOUND})
   MESSAGE(FATAL_ERROR "\n"
@@ -26,7 +26,8 @@ DEAL_II_SETUP_TARGET(strain_md)
 
 ## Include LAMMPS sources repository
 INCLUDE_DIRECTORIES(
-  /shared/ucl/apps/lammps/16Mar18/USERINTEL/intel-2018/lammps-16Mar18/src/
+  /home/uccamva/sources/lammps-17Nov16/src/
+  /home/uccamva/.conda/envs/surrogate/include/python3.8/
   )
 
 ## Include BOOST sources repository
@@ -34,9 +35,10 @@ INCLUDE_DIRECTORIES(
 #INCLUDE_DIRECTORIES( ${Boost_INCLUDE_DIR} )
 
 ## Include LAMMPS library shared .so (or static .a)
-TARGET_LINK_LIBRARIES(dealammps ~/.local/lib64/liblammps.so)
-TARGET_LINK_LIBRARIES(init_material ~/.local/lib64/liblammps.so)
-TARGET_LINK_LIBRARIES(strain_md ~/.local/lib64/liblammps.so)
+TARGET_LINK_LIBRARIES(dealammps /home/uccamva/sources/lammps-17Nov16/src/liblammps.so)
+TARGET_LINK_LIBRARIES(dealammps /home/uccamva/.conda/envs/surrogate/lib/libpython3.8.so)
+TARGET_LINK_LIBRARIES(init_material /home/uccamva/sources/lammps-17Nov16/src/liblammps.so)
+TARGET_LINK_LIBRARIES(strain_md /home/uccamva/sources/lammps-17Nov16/src/liblammps.so)
 
 TARGET_LINK_LIBRARIES(dealammps LINK_PUBLIC ${Boost_LIBRARIES})
 TARGET_LINK_LIBRARIES(init_material LINK_PUBLIC ${Boost_LIBRARIES})
diff --git a/dealammps.cc b/dealammps.cc
index 3bb0dfc..bf731d8 100644
--- a/dealammps.cc
+++ b/dealammps.cc
@@ -159,7 +159,7 @@ namespace HMM
 		Tensor<1,dim> 			cg_dir;
 		boost::property_tree::ptree	input_config;
 
-		bool				activate_md_update;
+		int					stress_compute_method;
 		bool				approx_md_with_hookes_law;
 		bool				use_pjm_scheduler;
 
@@ -236,7 +236,7 @@ namespace HMM
 	    quadrature_formula 	= input_config.get<int>("continuum mesh.quadrature formula");
 
 	    // Scale-bridging parameters
-	    activate_md_update 	= input_config.get<bool>("scale-bridging.activate md update");
+	    stress_compute_method 	= input_config.get<int>("scale-bridging.stress computation method");
 	    approx_md_with_hookes_law	=input_config.get<bool>("scale-bridging.approximate md with hookes law");
 	    use_pjm_scheduler 	= input_config.get<bool>("scale-bridging.use pjm scheduler");
 
@@ -294,7 +294,7 @@ namespace HMM
 		
 		// Print a recap of all the parameters...
 		hcout << "Parameters listing:" << std::endl;
-		hcout << " - Activate MD updates (1 is true, 0 is false): "<< activate_md_update << std::endl;
+		hcout << " - Method to compute local stresses (0 - LAMMPS, 1 - Hooke's law, 2 - ML surrogate model): "<< stress_compute_method << std::endl;
 		hcout << " - Approximate MD sims with hookes law (1 is true, 0 is false): "<< approx_md_with_hookes_law << std::endl;
 		hcout << " - Use Pilot Job Manager to schedule MD jobs: "<< use_pjm_scheduler << std::endl;
 		hcout << " - FE timestep duration: "<< fe_timestep_length << std::endl;
@@ -523,7 +523,7 @@ namespace HMM
 										macrostatelocin, macrostatelocout,
 										macrostatelocres, macrologloc,
 										freq_checkpoint, freq_output_visu, freq_output_lhist, freq_output_lbcforce,
-										activate_md_update, mdtype, cg_dir,
+										stress_compute_method, mdtype, cg_dir,
 										twod_mesh_file, extrude_length, extrude_points, 
 										input_config, approx_md_with_hookes_law);
                                                                                 
diff --git a/examples/streched_polyhedron/inputs.json b/examples/streched_polyhedron/inputs.json
index c742d95..5801bf5 100644
--- a/examples/streched_polyhedron/inputs.json
+++ b/examples/streched_polyhedron/inputs.json
@@ -4,7 +4,7 @@
 		"strain rate": 0.002
  	},
   "scale-bridging":{
-    "activate md update": 1,
+    "stress computation method": 0,
     "approximate md with hookes law": 0,
     "use pjm scheduler": 0
   },
@@ -30,7 +30,7 @@
     "md":{
       "min quadrature strain norm": 1.0e-10
     },
-    "clustering":{
+  "clustering":{
       "spline points": 10,
       "min steps": 5000,
       "diff threshold": 0.000001,
@@ -56,14 +56,14 @@
   },
   "computational resources":{
     "machine cores per node": 24,
-    "maximum number of cores for FEM simulation": 1,
+		"maximum number of cores for FEM simulation": 1,
     "minimum number of cores for MD simulation": 1
   },
   "output data":{
     "checkpoint frequency": 1,
     "visualisation output frequency": 1,
     "analytics output frequency": 1,
-    "loaded boundary force output frequency": 1,
+		"loaded boundary force output frequency": 1,
     "homogenization output frequency": 1000
   },
   "directory structure":{
diff --git a/headers/FE.h b/headers/FE.h
index 9805dec..eada133 100644
--- a/headers/FE.h
+++ b/headers/FE.h
@@ -96,7 +96,7 @@ namespace HMM
 		SymmetricTensor<2,dim> upd_strain;
 		SymmetricTensor<2,dim> newton_strain;
 		MatHistPredict::Strain6D hist_strain;
-		bool to_be_updated;
+		bool to_be_updated_with_md;
 
 		// Characteristics
 		unsigned int qpid;
@@ -233,7 +233,7 @@ namespace HMM
 
 							void init (int sstp, double tlength, std::string mslocin, std::string mslocout,
 											std::string mslocres, std::string mlogloc, int fchpt, int fovis, int folhis, int folbcf,
-											bool actmdup, std::vector<std::string> mdt, Tensor<1,dim> cgd, 
+											int stress_method, std::vector<std::string> mdt, Tensor<1,dim> cgd,
 											std::string twodmfile, double extrudel, int extrudep, 
 											boost::property_tree::ptree inconfig, bool hookeslaw);
 							void beginstep (int tstp, double ptime);
@@ -272,6 +272,9 @@ namespace HMM
 							void gather_qp_update_list(ScaleBridgingData &scale_bridging_data);
 							template <typename T>
 							std::vector<T> gather_vector(std::vector<T> local_vector);
+							SymmetricTensor<2,dim> compute_stress_with_surrogate(SymmetricTensor<2,dim> old_strain,
+																				 SymmetricTensor<2,dim> new_strain,
+																				 SymmetricTensor<2,dim> old_stress);
 							void update_stress_quadrature_point_history
 									(const Vector<double>& displacement_update, ScaleBridgingData scale_bridging_data);
 							void clean_transfer();
@@ -360,7 +363,7 @@ namespace HMM
 							int									freq_output_lhist;
 							int									freq_output_lbcforce;
 
-							bool 								activate_md_update;
+							int 								stress_compute_method;
 
 							std::string		twod_mesh_file;
 							double                  extrude_length;
diff --git a/headers/FE_problem.h b/headers/FE_problem.h
index ff1be0b..756b1b7 100644
--- a/headers/FE_problem.h
+++ b/headers/FE_problem.h
@@ -83,6 +83,9 @@
 #include "compact_tension.h"
 #include "FE.h"
 
+// For Dongwei's surrogate function
+#include "Python.h"
+
 namespace HMM
 {
 	using namespace dealii;
@@ -477,7 +480,7 @@ namespace HMM
 											{
 													local_quadrature_points_history[q].new_strain = 0;
 													local_quadrature_points_history[q].upd_strain = 0;
-													local_quadrature_points_history[q].to_be_updated = false;
+													local_quadrature_points_history[q].to_be_updated_with_md = false;
 													local_quadrature_points_history[q].new_stress = 0;
 													local_quadrature_points_history[q].qpid = cell->active_cell_index()*quadrature_formula.size() + q;
 
@@ -1138,20 +1141,20 @@ namespace HMM
 					//   (i) all cells,
 					//  (ii) cells in given location,
 					// (iii) cells based on their id
-					if (activate_md_update
+					if (stress_compute_method == 0
 						// otherwise MD simulation unecessary, because no significant volume change and MD will fail
 										&& (local_quadrature_points_history[q].upd_strain.norm() >= min_qp_strain//> 1.0e-10
-											|| local_quadrature_points_history[q].to_be_updated == true)
+											|| local_quadrature_points_history[q].to_be_updated_with_md == true)
 						)
-					//if (activate_md_update && cell->barycenter()(1) <  3.0*tt && cell->barycenter()(0) <  1.10*(ww - aa) && cell->barycenter()(0) > 0.0*(ww - aa))
-					/*if (activate_md_update && (cell->active_cell_index() == 2922 || cell->active_cell_index() == 2923
+					//if (stress_compute_method && cell->barycenter()(1) <  3.0*tt && cell->barycenter()(0) <  1.10*(ww - aa) && cell->barycenter()(0) > 0.0*(ww - aa))
+					/*if (stress_compute_method && (cell->active_cell_index() == 2922 || cell->active_cell_index() == 2923
 						|| cell->active_cell_index() == 2924 || cell->active_cell_index() == 2487
 						|| cell->active_cell_index() == 2488 || cell->active_cell_index() == 2489))*/ // For debug...
 					{
-						local_quadrature_points_history[q].to_be_updated = true;
+						local_quadrature_points_history[q].to_be_updated_with_md = true;
 					}
 					else{
-						local_quadrature_points_history[q].to_be_updated = false;
+						local_quadrature_points_history[q].to_be_updated_with_md = false;
 					}
 				}
 			}
@@ -1212,7 +1215,7 @@ namespace HMM
 						ExcInternalError());
 				for (unsigned int q=0; q<quadrature_formula.size(); ++q)
 				{
-					if(local_quadrature_points_history[q].to_be_updated)
+					if(local_quadrature_points_history[q].to_be_updated_with_md)
 					{
 						histories.push_back(&local_quadrature_points_history[q].hist_strain);
 						//local_quadrature_points_history[q].hist_strain.print();
@@ -1312,7 +1315,7 @@ namespace HMM
 						&quadrature_point_history.back(),
 						ExcInternalError());
 				for (unsigned int q=0; q<quadrature_formula.size(); ++q)
-					if(local_quadrature_points_history[q].to_be_updated
+					if(local_quadrature_points_history[q].to_be_updated_with_md
 							&& local_quadrature_points_history[q].hist_strain.run_new_md())
 					{
 						// The cell will get its stress from MD, but should it run an MD simulation?
@@ -1483,9 +1486,151 @@ namespace HMM
 		return qp;
 	}
 
+	template <int dim>
+	SymmetricTensor<2,dim> FEProblem<dim>::compute_stress_with_surrogate(SymmetricTensor<2,dim> old_strain,
+																		 SymmetricTensor<2,dim> new_strain,
+																		 SymmetricTensor<2,dim> old_stress)
+	{
+		SymmetricTensor<2,dim> new_stress;
+
+		// Create sudo input data
+		std::vector<float> strain_pre; strain_pre.assign(2*dim, 1.0e-20);
+		std::vector<float> strain_cur; strain_cur.assign(2*dim, 1.0e-20);
+		std::vector<float> stress_pre; stress_pre.assign(2*dim, 1.0e-20);
+		//std::vector<float> stress_pre = {1.0e-20, 1.0e-20, 1.0e-20, 1.0e-20, 1.0e-20, 1.0e-20};
+		std::vector<float> inputs;
+
+		for(unsigned int k=0;k<2*dim;k++)
+		{
+			if (k<dim){
+				strain_cur[k] = new_strain[0][k];
+				strain_pre[k] = old_strain[0][k];
+				stress_pre[k] = old_stress[0][k];
+			}
+			else if(k==3){
+				strain_cur[k] = new_strain[1][1];
+				strain_pre[k] = old_strain[1][1];
+				stress_pre[k] = old_stress[1][1];
+			}
+			else if(k==4){
+				strain_cur[k] = new_strain[1][2];
+				strain_pre[k] = old_strain[1][2];
+				stress_pre[k] = old_stress[1][2];
+			}
+			else if(k==5){
+				strain_cur[k] = new_strain[2][2];
+				strain_pre[k] = old_strain[2][2];
+				stress_pre[k] = old_stress[2][2];
+			}
+		}
+
+		// Concatenate the data together, easy for passing to python
+		// I do it a bit urgly here... but it works
+		// You know c++ better than me, you can improve it
+		inputs.insert(inputs.end(), strain_cur.begin(), strain_cur.end());
+		inputs.insert(inputs.end(), strain_pre.begin(), strain_pre.end());
+		inputs.insert(inputs.end(), stress_pre.begin(), stress_pre.end());
+
+		//Initialize the py function if it is not activated
+		if ( !Py_IsInitialized() )
+		{
+			Py_Initialize();
+		}
+		else
+		{
+			//std::cout<<"Pyfunction activated already"<<std::endl;
+		}
+
+		/////////// From C++ to python ////////////////
+		// Start to execute the python file
+		Py_Initialize();
+		if ( !Py_IsInitialized() )
+		{
+			std::cerr << "Python kernel could not be initialised to execute the surrogate "
+					"model to compute the stress of a given QP." << std::endl;
+			exit(1);
+		}
+
+		// import some python lib
+		PyRun_SimpleString("import sys");
+		PyRun_SimpleString("import os");
+		// Check current directory
+		//PyRun_SimpleString("print(os.getcwd())");
+		// Add current directory to the directory
+		PyRun_SimpleString("sys.path.append('./')");
+		PyRun_SimpleString("sys.path.append('./surrogate_model/')");
+
+		// Set up Null variable for data passing
+		PyObject* pModule =NULL;
+		PyObject* pList = NULL;
+		PyObject* pFunc = NULL;
+		PyObject* pArgs = NULL;
+
+		// Open the .py under current direcotory
+		pModule = PyImport_ImportModule("surrogate");
+		// Load function 'surrogate_model'
+		pFunc = PyObject_GetAttrString(pModule, "surrogate_model");
+		// Open up an argument with Tuple for python to take in (size of tuple)
+		pArgs = PyTuple_New(1);
+		// Open up a list (size of len, start from size 0 for append)
+		pList = PyList_New(0);
+
+		// Add value to the list
+		for (int i = 0; i < inputs.size(); i++)
+		{
+			// Append the python list with values, 'f' means datatype float
+			PyList_Append(pList, Py_BuildValue("f", inputs[i]));
+		}
+		// Pass the list to Args
+		PyTuple_SetItem(pArgs, 0, pList);
+		// Call the surrogate_model function with args: inputs
+		PyObject* pRet = PyEval_CallObject(pFunc, pArgs);
+
+		////////////////// From python to c++ /////////////////
+		// Return value has to be a 1D list
+		// Check size of list
+		int SizeOfList = PyList_Size(pRet);
+
+		// Create a empty c++ array with size of return list:
+		std::vector<double> result;
+
+		for(int i = 0; i < SizeOfList; i++)
+		{
+			// Get every element in the list to ListItem using PyList_GetItem
+			PyObject *ListItem = PyList_GetItem(pRet, i);
+			// Take data from python float to C++ double
+			result.push_back(PyFloat_AS_DOUBLE(ListItem));
+		}
+
+		Py_DECREF(pModule);
+		Py_DECREF(pList);
+		Py_DECREF(pFunc);
+		Py_DECREF(pArgs);
+		Py_DECREF(pRet);
+		//Py_Finalize();
+
+		for(unsigned int k=0;k<2*dim;k++)
+		{
+			if (k<dim){
+				new_stress[0][k] = result[k];
+			}
+			else if(k==3){
+				new_stress[1][1] = result[k];
+			}
+			else if(k==4){
+				new_stress[1][2] = result[k];
+			}
+			else if(k==5){
+				new_stress[2][2] = result[k];
+			}
+		}
+
+		return new_stress;
+	}
+
 	template <int dim>
 	void FEProblem<dim>::update_stress_quadrature_point_history(const Vector<double>& displacement_update,
-																															ScaleBridgingData scale_bridging_data)
+			ScaleBridgingData scale_bridging_data)
 	{
 		FEValues<dim> fe_values (fe, quadrature_formula,
 				update_values | update_gradients);
@@ -1520,97 +1665,90 @@ namespace HMM
 				for (unsigned int q=0; q<quadrature_formula.size(); ++q)
 				{
 					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.;
-
-					if (local_quadrature_points_history[q].to_be_updated){
 
-						// Updating stiffness tensor
-						/*SymmetricTensor<4,dim> stmp_stiff;
-						sprintf(filename, "%s/last.%s.stiff", macrostatelocout.c_str(), cell_id);
-						read_tensor<dim>(filename, stmp_stiff);
-
-						// Rotate the output stiffness wrt the flake angles
-						local_quadrature_points_history[q].new_stiff =
-								rotate_tensor(stmp_stiff, transpose(local_quadrature_points_history[q].rotam));
-						 */
+					if (newtonstep == 0) local_quadrature_points_history[q].inc_stress = 0.;
 
-						// Updating stress tensor
-						//bool load_stress = true;
+					if (stress_compute_method==0){
+						if (local_quadrature_points_history[q].to_be_updated_with_md){
 
-						/*SymmetricTensor<4,dim> loc_stiffness;
-						sprintf(filename, "%s/last.%s.stiff", macrostatelocout.c_str(), cell_id);
-						read_tensor<dim>(filename, loc_stiffness);*/
+							QP qp;
+							qp = get_qp_with_id(qp_id, scale_bridging_data);
+							//sprintf(filename, "%s/last.%s.stress", macrostatelocout.c_str(), cell_id);
+							//load_stress = read_tensor<dim>(filename, loc_stress);
+							//std::cout << "Putting stress into quadrature_points_history" << std::endl;
+							//for (int i=0; i<6; i++){
+							//  std::cout << " " << qp.update_stress[i];
+							//} std::cout << std::endl;
+							SymmetricTensor<2,dim> loc_stress(qp.update_stress);
 
-						QP qp;
-						qp = get_qp_with_id(qp_id, scale_bridging_data);
-						//sprintf(filename, "%s/last.%s.stress", macrostatelocout.c_str(), cell_id);
-						//load_stress = read_tensor<dim>(filename, loc_stress);
-			//std::cout << "Putting stress into quadrature_points_history" << std::endl;
-      //for (int i=0; i<6; i++){
-      //  std::cout << " " << qp.update_stress[i];
-      //} std::cout << std::endl;
-						SymmetricTensor<2,dim> loc_stress(qp.update_stress);
+							// Rotate the output stress wrt the flake angles
+							loc_stress = rotate_tensor(loc_stress, transpose(local_quadrature_points_history[q].rotam));
 
-						// Rotate the output stress wrt the flake angles
-						loc_stress = rotate_tensor(loc_stress, transpose(local_quadrature_points_history[q].rotam));
+							if (approx_md_with_hookes_law == false){
+								local_quadrature_points_history[q].new_stress = loc_stress;
+							}
+							else {
+								local_quadrature_points_history[q].new_stress = loc_stress + local_quadrature_points_history[q].old_stress;
+							}
 
-						if (approx_md_with_hookes_law == false){
-							local_quadrature_points_history[q].new_stress = loc_stress;
+							// Resetting the update strain tensor
+							local_quadrature_points_history[q].upd_strain = 0;
 						}
 						else {
-							local_quadrature_points_history[q].new_stress = loc_stress + local_quadrature_points_history[q].old_stress;
-						}
-
-						// Resetting the update strain tensor
-						local_quadrature_points_history[q].upd_strain = 0;
+							local_quadrature_points_history[q].new_stress +=                                                                                                    local_quadrature_points_history[q].new_stiff*local_quadrature_points_history[q].newton_strain;                                                                                                                                            						}
 					}
-					else{
-					// Tangent stiffness computation of the new stress tensor
+					else if (stress_compute_method==1
+							//|| (local_quadrature_points_history[q].qpid != 0)
+							){
+						// Tangent stiffness computation of the new stress tensor
 						local_quadrature_points_history[q].new_stress +=
-							local_quadrature_points_history[q].new_stiff*local_quadrature_points_history[q].newton_strain;
+								local_quadrature_points_history[q].new_stiff*local_quadrature_points_history[q].newton_strain;
+					}
+					else if (stress_compute_method==2){
+						local_quadrature_points_history[q].new_stress = compute_stress_with_surrogate(local_quadrature_points_history[q].old_strain,
+								local_quadrature_points_history[q].new_strain,
+								local_quadrature_points_history[q].old_stress);
+					}
+					else {
+						std::cerr << "Local stress computation method not implemented." << std::endl;
+						exit(1);
 					}
 				}
-			
-					// Secant stiffness computation of the new stress tensor
-					//local_quadrature_points_history[q].new_stress =
-					//		local_quadrature_points_history[q].new_stiff*local_quadrature_points_history[q].new_strain;
-
-					// Write stress tensor for each gauss point
-					/*sprintf(filename, "%s/last.%s-%d.stress", macrostatelocout.c_str(), cell_id,q);
-					write_tensor<dim>(filename, local_quadrature_points_history[q].new_stress);*/
-
-					// Apply rotation of the sample to the new state tensors.
-					// Only needed if the mesh is modified...
-					/*const Tensor<2,dim> rotation
-					= get_rotation_matrix (displacement_update_grads[q]);
-
-					const SymmetricTensor<2,dim> rotated_new_stress
-					= symmetrize(transpose(rotation) *
-							static_cast<Tensor<2,dim> >
-					(local_quadrature_points_history[q].new_stress) *
-					rotation);
-
-					const SymmetricTensor<2,dim> rotated_new_strain
-					= symmetrize(transpose(rotation) *
-							static_cast<Tensor<2,dim> >
-					(local_quadrature_points_history[q].new_strain) *
-					rotation);
-
-					const SymmetricTensor<2,dim> rotated_upd_strain
-					= symmetrize(transpose(rotation) *
-							static_cast<Tensor<2,dim> >
-					(local_quadrature_points_history[q].upd_strain) *
-					rotation);
-
-					local_quadrature_points_history[q].new_stress
-					= rotated_new_stress;
-					local_quadrature_points_history[q].new_strain
-					= rotated_new_strain;
-					local_quadrature_points_history[q].upd_strain
-					= rotated_upd_strain;*/
+			}
+		}
+		/*MPI_Barrier(FE_communicator);
+		// Retrieving all quadrature points computation and storing them in the
+		// quadrature_points_history structure
+		for (typename DoFHandler<dim>::active_cell_iterator
+				cell = dof_handler.begin_active();
+				cell != dof_handler.end(); ++cell){
+			if (cell->is_locally_owned())
+			{
+				SymmetricTensor<2,dim> avg_upd_strain_tensor;
+				//SymmetricTensor<2,dim> avg_stress_tensor;
+
+				PointHistory<dim> *local_quadrature_points_history
+				= reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
+				Assert (local_quadrature_points_history >=
+						&quadrature_point_history.front(),
+						ExcInternalError());
+				Assert (local_quadrature_points_history <
+						&quadrature_point_history.back(),
+						ExcInternalError());
+				fe_values.reinit (cell);
+				fe_values.get_function_gradients (displacement_update,
+						displacement_update_grads);
+
+				for (unsigned int q=0; q<quadrature_formula.size(); ++q)
+				{
+					if (local_quadrature_points_history[q].qpid % 100 != 0){
+						std::cout << local_quadrature_points_history[q].qpid << " " << std::flush;
+						std::cout << local_quadrature_points_history[q].new_strain << " " << std::flush;
+						std::cout << local_quadrature_points_history[q].new_stress << " " << std::endl;
+					}
 				}
 			}
+		}*/
 	}
 
 
@@ -1639,7 +1777,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].qpid);
 
-					if(local_quadrature_points_history[q].to_be_updated
+					if(local_quadrature_points_history[q].to_be_updated_with_md
 							&& local_quadrature_points_history[q].hist_strain.run_new_md()){
 						// Removing stiffness passing file
 						//sprintf(filename, "%s/last.%s.stiff", macrostatelocout.c_str(), cell_id);
@@ -2202,7 +2340,7 @@ namespace HMM
 	void FEProblem<dim>::init (int sstp, double tlength,
 							   std::string mslocin, std::string mslocout,
 							   std::string mslocres, std::string mlogloc,
-							   int fchpt, int fovis, int folhis, int folbcf, bool actmdup,
+							   int fchpt, int fovis, int folhis, int folbcf, int stress_method,
 							   std::vector<std::string> mdt, Tensor<1,dim> cgd,
 							   std::string twodmfile, double extrudel, int extrudep,
 						     boost::property_tree::ptree inconfig, 
@@ -2218,7 +2356,7 @@ namespace HMM
 		freq_output_lbcforce = folbcf;
 
 		// Setting up usage of MD to update constitutive behaviour
-		activate_md_update = actmdup;
+		stress_compute_method = stress_method;
 
 		// Setting up starting timestep number and timestep length
 		start_timestep = sstp;
diff --git a/input_configurations/inputs_compact.json b/input_configurations/inputs_compact.json
index 946b80a..a842145 100644
--- a/input_configurations/inputs_compact.json
+++ b/input_configurations/inputs_compact.json
@@ -4,7 +4,7 @@
     "velocity": 0.00010
  	},
   "scale-bridging":{
-    "activate md update": 1,
+    "stress computation method": 0,
     "approximate md with hookes law": 0,
     "use pjm scheduler": 0
   },
diff --git a/input_configurations/inputs_dogbone_cuboid.json b/input_configurations/inputs_dogbone_cuboid.json
index d4f1d03..e6b108e 100644
--- a/input_configurations/inputs_dogbone_cuboid.json
+++ b/input_configurations/inputs_dogbone_cuboid.json
@@ -4,7 +4,7 @@
 		"strain rate": 0.002
  	},
   "scale-bridging":{
-    "activate md update": 1,
+    "stress computation method": 0,
 		"approximate md with hookes law": 0,
     "use pjm scheduler": 0
   },
diff --git a/input_configurations/inputs_dogbone_file3D.json b/input_configurations/inputs_dogbone_file3D.json
index 2421ce5..b36b96b 100644
--- a/input_configurations/inputs_dogbone_file3D.json
+++ b/input_configurations/inputs_dogbone_file3D.json
@@ -4,7 +4,7 @@
     "strain rate": 0.010
   },
   "scale-bridging":{
-    "activate md update": 1,
+    "stress computation method": 0,
     "approximate md with hookes law": 0,
     "use pjm scheduler": 0
   },
diff --git a/input_configurations/inputs_dropweight_cuboid.json b/input_configurations/inputs_dropweight_cuboid.json
index a3253e0..1a607b8 100644
--- a/input_configurations/inputs_dropweight_cuboid.json
+++ b/input_configurations/inputs_dropweight_cuboid.json
@@ -6,7 +6,7 @@
 		"steps to accelerate": 5
  	},
   "scale-bridging":{
-    "activate md update": 1,
+    "stress computation method": 0,
 		"approximate md with hookes law": 0,
     "use pjm scheduler": 0
   },
diff --git a/surrogate_model/README b/surrogate_model/README
new file mode 100644
index 0000000..1a6db22
--- /dev/null
+++ b/surrogate_model/README
@@ -0,0 +1,7 @@
+Load standard modules required to build and run deal.ii and lammps
+
+Set the python environment in order to be able to import Keras and Tensorflow:
+source env_surrogate # which should include keras and tensorflow > 2.2
+
+Set a symlink in the execution directory to surrogate.py:
+ln -s /path/to/SCEMa/surrogate_model/surrogate.py
diff --git a/surrogate_model/model_small_uniaxial.bin b/surrogate_model/model_small_uniaxial.bin
new file mode 100644
index 0000000..140f0bb
Binary files /dev/null and b/surrogate_model/model_small_uniaxial.bin differ
diff --git a/surrogate_model/scaler.pkl b/surrogate_model/scaler.pkl
new file mode 100644
index 0000000..8b5d804
Binary files /dev/null and b/surrogate_model/scaler.pkl differ
diff --git a/surrogate_model/surrogate.py b/surrogate_model/surrogate.py
new file mode 100644
index 0000000..0e32017
--- /dev/null
+++ b/surrogate_model/surrogate.py
@@ -0,0 +1,29 @@
+# The data passed to py is a list
+
+import numpy as np
+from keras import models
+from pickle import load
+
+#import tensorflow as tf
+def surrogate_model(nlist):
+  # Convert inputs list to numpy array
+  #print(type(nlist[0]))
+  #print(nlist)
+  inputs = np.array(nlist)
+  inputs_reshape = inputs.reshape(1, -1)
+
+  # Load scaler and transform the X_test
+  scaler = load(open('./surrogate_model/scaler.pkl', 'rb'))
+  inputs_scaled = scaler.transform(inputs_reshape)
+  #print(inputs_scaled)
+
+  # Load model and predict
+  model = models.load_model("./surrogate_model/model_small_uniaxial.bin")
+  output_pred = model.predict(inputs_scaled)
+  #print(output_pred)
+
+  # convert the data back to a list
+  output_list = output_pred.tolist()[0]
+  #print(output_list)
+
+  return output_list