diff --git a/model.cpp b/model.cpp
index 2fcce0e861b2eec7be78e3d6b35d85fe6e1f50c3..f84a692079329a8fe5047bdf2128b0b135d1fae5 100644
--- a/model.cpp
+++ b/model.cpp
@@ -182,20 +182,21 @@ int bioem_model::readModel(const char* filemodel)
 	cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
 	exampleReadCoor.close();
 
-	//Moving to Model to its center of mass:
+	//Moving to Model to its center of density mass:
 	myfloat3_t r_cm;
 
 	for(int n = 0; n < 3; n++)r_cm.pos[n] = 0.0;
 
 	for(int n = 0; n < nPointsModel; n++)
 	{
-		r_cm.pos[0] += points[n].point.pos[0];
-		r_cm.pos[1] += points[n].point.pos[1];
-		r_cm.pos[2] += points[n].point.pos[2];
+		r_cm.pos[0] += points[n].point.pos[0]*points[n].density;
+		r_cm.pos[1] += points[n].point.pos[1]*points[n].density;
+		r_cm.pos[2] += points[n].point.pos[2]*points[n].density;
 	}
-	r_cm.pos[0] = r_cm.pos[0] / (myfloat_t) nPointsModel;
-	r_cm.pos[1] = r_cm.pos[1] / (myfloat_t) nPointsModel;
-	r_cm.pos[2] = r_cm.pos[2] / (myfloat_t) nPointsModel;
+	r_cm.pos[0] = r_cm.pos[0] / NormDen;
+	r_cm.pos[1] = r_cm.pos[1] / NormDen;
+	r_cm.pos[2] = r_cm.pos[2] / NormDen;
+
 	for(int n = 0; n < nPointsModel; n++)
 	{
 		points[n].point.pos[0] -= r_cm.pos[0];