diff --git a/python/ester_wrap.cxx b/python/ester_wrap.cxx
index 1cdb2411..c0b9a35a 100644
--- a/python/ester_wrap.cxx
+++ b/python/ester_wrap.cxx
@@ -14844,37 +14844,6 @@ SWIGINTERN PyObject *_wrap_atm_onelayer(PyObject *SWIGUNUSEDPARM(self), PyObject
 }
 
 
-SWIGINTERN PyObject *_wrap_initial_composition(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
-  PyObject *resultobj = 0;
-  double arg1 ;
-  double arg2 ;
-  double val1 ;
-  int ecode1 = 0 ;
-  double val2 ;
-  int ecode2 = 0 ;
-  PyObject * obj0 = 0 ;
-  PyObject * obj1 = 0 ;
-  double_map result;
-  
-  if (!PyArg_ParseTuple(args,(char *)"OO:initial_composition",&obj0,&obj1)) SWIG_fail;
-  ecode1 = SWIG_AsVal_double(obj0, &val1);
-  if (!SWIG_IsOK(ecode1)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode1), "in method '" "initial_composition" "', argument " "1"" of type '" "double""'");
-  } 
-  arg1 = static_cast< double >(val1);
-  ecode2 = SWIG_AsVal_double(obj1, &val2);
-  if (!SWIG_IsOK(ecode2)) {
-    SWIG_exception_fail(SWIG_ArgError(ecode2), "in method '" "initial_composition" "', argument " "2"" of type '" "double""'");
-  } 
-  arg2 = static_cast< double >(val2);
-  result = initial_composition(arg1,arg2);
-  resultobj = SWIG_NewPointerObj((new double_map(static_cast< const double_map& >(result))), SWIGTYPE_p_double_map, SWIG_POINTER_OWN |  0 );
-  return resultobj;
-fail:
-  return NULL;
-}
-
-
 SWIGINTERN PyObject *_wrap_mapping_gl_set(PyObject *SWIGUNUSEDPARM(self), PyObject *args) {
   PyObject *resultobj = 0;
   mapping *arg1 = (mapping *) 0 ;
@@ -31325,7 +31294,6 @@ static PyMethodDef SwigMethods[] = {
 	 { (char *)"eos_idealrad", _wrap_eos_idealrad, METH_VARARGS, NULL},
 	 { (char *)"eos_opal", _wrap_eos_opal, METH_VARARGS, NULL},
 	 { (char *)"atm_onelayer", _wrap_atm_onelayer, METH_VARARGS, NULL},
-	 { (char *)"initial_composition", _wrap_initial_composition, METH_VARARGS, NULL},
 	 { (char *)"mapping_gl_set", _wrap_mapping_gl_set, METH_VARARGS, NULL},
 	 { (char *)"mapping_gl_get", _wrap_mapping_gl_get, METH_VARARGS, NULL},
 	 { (char *)"mapping_leg_set", _wrap_mapping_leg_set, METH_VARARGS, NULL},
diff --git a/python/ester_wrap.py b/python/ester_wrap.py
index 73436d59..862715a0 100644
--- a/python/ester_wrap.py
+++ b/python/ester_wrap.py
@@ -902,9 +902,6 @@ def atm_onelayer(X, Z, g, Teff, eos_name, opa_name, atm):
     return _ester_wrap.atm_onelayer(X, Z, g, Teff, eos_name, opa_name, atm)
 atm_onelayer = _ester_wrap.atm_onelayer
 
-def initial_composition(X, Z):
-    return _ester_wrap.initial_composition(X, Z)
-initial_composition = _ester_wrap.initial_composition
 MAP_BONAZZOLA = _ester_wrap.MAP_BONAZZOLA
 MAP_LINEAR = _ester_wrap.MAP_LINEAR
 class mapping(_object):
diff --git a/src/include/matrix.h b/src/include/matrix.h
index f2e1303f..5f61b753 100644
--- a/src/include/matrix.h
+++ b/src/include/matrix.h
@@ -289,7 +289,9 @@ class matrix_map : public std::map<std::string,matrix> {
 	inline matrix &operator[](const std::string& key) {return std::map<std::string,matrix>::operator[](key);}
 	inline const matrix &operator[](const std::string& key) const {return find(key)->second;};
 	matrix_map_elem operator()(int nfil, int ncol);
+	matrix_map_elem operator()(int ielem);
 	const double_map operator()(int nfil, int ncol) const;
+	const double_map operator()(int ielem) const;
 	const matrix_map row(int irow) const;
 	const matrix_map col(int icol) const;
 	const matrix_map block(int irow1,int irow2,int icol1,int icol2) const;
diff --git a/src/include/physics.h b/src/include/physics.h
index 223305cd..190f5caa 100644
--- a/src/include/physics.h
+++ b/src/include/physics.h
@@ -34,12 +34,12 @@ class composition_map : public matrix_map {
 
 int opa_calc(const matrix &X,double Z,const matrix &T,const matrix &rho,
 		opa_struct &opa);
-int eos_calc(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos);
-int nuc_calc(const matrix_map &X,const matrix &T,const matrix &rho,
+int eos_calc(const composition_map &, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos);
+int nuc_calc(const composition_map &, const matrix &T, const matrix &rho,
 		nuc_struct &nuc);
-int atm_calc(const matrix &X,double Z,const matrix &g,const matrix &Teff,
-		const char *eos_name,const char *opa_name,atm_struct &atm);
+int atm_calc(const composition_map &, const matrix &g, const matrix &Teff,
+		const char *eos_name, const char *opa_name, atm_struct &atm);
 
 int opa_opal(const matrix &X,double Z,const matrix &T,const matrix &rho,
 		opa_struct &opa);
@@ -58,19 +58,17 @@ int nuc_cesam(const composition_map &comp,const matrix &T,const matrix &rho,
 int nuc_cesam_dcomp(composition_map &comp,const matrix &T,const matrix &rho,
 		nuc_struct &nuc);
 
-int eos_ideal(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos);
-int eos_idealrad(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos);
-int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos);
-int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
+int eos_ideal(const composition_map &, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos);
+int eos_idealrad(const composition_map &, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos);
+int eos_opal(const composition_map &, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos);
+int eos_freeeos(const composition_map &, const matrix &T, const matrix &p,
 		matrix &rho, eos_struct &eos);
 
-int atm_onelayer(const matrix &X,double Z,const matrix &g,const matrix &Teff,
-		const char *eos_name,const char *opa_name,atm_struct &atm);
-
-double_map initial_composition(double X,double Z);
+int atm_onelayer(const composition_map &, const matrix &g, const matrix &Teff,
+		const char *eos_name, const char *opa_name, atm_struct &atm);
 
 #endif
 
diff --git a/src/include/star.h b/src/include/star.h
index 04e34b00..a0f7c313 100644
--- a/src/include/star.h
+++ b/src/include/star.h
@@ -37,6 +37,9 @@ class star2d {
         matrix rho, phi, p, T;
         matrix phiex;
         matrix vr, vt, G, w;
+        // the m_ prefix is used to explicit that it's a class' member https://stackoverflow.com/a/13018190/22158934 :
+        double m_He_isotopic_ratio; // He3/He4 hardcoded in star*d::init_comp // TODO: change it to be configurable
+        double_map m_metal_mix;
         composition_map comp;
         double X0, Y0, Z0;
         double R, M;
@@ -97,6 +100,7 @@ class star2d {
 
         virtual void dump_info();
 
+        virtual void init_metal_mix();
         virtual void init_comp();
 
         virtual solver *init_solver(int nvar_add=0);
diff --git a/src/matrix/matrix_map.cpp b/src/matrix/matrix_map.cpp
index 58336024..d756ac70 100644
--- a/src/matrix/matrix_map.cpp
+++ b/src/matrix/matrix_map.cpp
@@ -51,6 +51,17 @@ matrix_map_elem matrix_map::operator()(int nfil, int ncol) {
 	return elem;
 }
 
+matrix_map_elem matrix_map::operator()(int ielem) {
+
+	matrix_map_elem elem;
+	matrix_map::iterator it;
+	for(it=begin();it!=end();it++) {
+		elem[it->first]=&(it->second)(ielem);
+	}
+
+	return elem;
+}
+
 const double_map matrix_map::operator()(int nfil, int ncol) const {
 	
 	double_map elem;
@@ -62,6 +73,17 @@ const double_map matrix_map::operator()(int nfil, int ncol) const {
 	return elem;
 }
 
+const double_map matrix_map::operator()(int ielem) const {
+
+	double_map elem;
+	matrix_map::const_iterator it;
+	for(it=begin();it!=end();it++) {
+		elem[it->first]=(it->second)(ielem);
+	}
+
+	return elem;
+}
+
 const matrix_map matrix_map::row(int irow) const {
 
 	matrix_map res;
diff --git a/src/physics/EOS5_xtrin.F b/src/physics/EOS5_xtrin.F
index fe6d1043..d0af5778 100644
--- a/src/physics/EOS5_xtrin.F
+++ b/src/physics/EOS5_xtrin.F
@@ -121,7 +121,7 @@ subroutine esac (xh,ztab,t6,r,iorder,irad)
 c..... read the data files
         call readcoeos
         z=zz(1)
-        if (ztab .ne. z) then
+        if (abs(ztab-z) .ge. 1e-7) then
         write (*,'("requested z=",f10.6," EOS5_data is for z=",
      x             f10.6)')
      x    ztab,z      
diff --git a/src/physics/atm_onelayer.cpp b/src/physics/atm_onelayer.cpp
index db366430..9f6dbd18 100644
--- a/src/physics/atm_onelayer.cpp
+++ b/src/physics/atm_onelayer.cpp
@@ -7,8 +7,8 @@
 #include <string.h>
 #include <cmath>
 
-int atm_onelayer(const matrix &X,double Z,const matrix &g,const matrix &Teff,
-		const char *eos_name,const char *opa_name,atm_struct &atm) {
+int atm_onelayer(const composition_map &chemical_comp, const matrix &g, const matrix &Teff,
+		const char *eos_name, const char *opa_name, atm_struct &atm) {
 
 	int n=g.nrows();
 	int m=g.ncols();
@@ -43,8 +43,8 @@ int atm_onelayer(const matrix &X,double Z,const matrix &g,const matrix &Teff,
 				}
 				double F,dF,dlogps;
 				p=pow(10,logps)*ones(1,1);
-				if(eos_calc(X(i,j)*ones(1,1),Z,T,p,rho,eos)) return 1;
-				if(opa_calc(X(i,j)*ones(1,1),Z,T,rho,opa)) return 1;
+				if(eos_calc(chemical_comp(i,j)*ones(1,1), T, p, rho, eos)) return 1;
+				if(opa_calc(chemical_comp.X()(i,j)*ones(1,1), chemical_comp.Z()(i,j), T, rho, opa)) return 1;
 				F=logps+log10(opa.k)(0)-logg-log10(2./3.);
 				dF=1.-(1+opa.dlnxi_lnrho(0))/eos.chi_rho(0);
 				dlogps=-F/dF;
diff --git a/src/physics/composition.cpp b/src/physics/composition.cpp
index 1deacba7..4a6fa36c 100644
--- a/src/physics/composition.cpp
+++ b/src/physics/composition.cpp
@@ -1,30 +1,8 @@
 #ifndef WITH_CMAKE
 #include "ester-config.h"
 #endif
-#include "physics.h"
-
-double_map initial_composition(double X, double Z) {
-	double_map comp;
-	
-// The mixture below has an unidentified origin but is close to the one
-// of Grevesse & Sauval 1998.
-	comp["H"]=X;
-	comp["He3"]=3.15247417638132e-04*(1.-X-Z);
-	comp["He4"]=(1.-X-Z)-comp["He3"];
-	comp["C12"]=Z*1.71243418737847e-01;
-	comp["C13"]=Z*2.06388003380057e-03;
-	comp["N14"]=Z*5.29501630871695e-02;
-	comp["N15"]=Z*2.08372414940812e-04;
-	comp["O16"]=Z*4.82006487350336e-01;
-	comp["O17"]=Z*1.95126448826986e-04;
-	
-	double tot=comp.sum();
-	
-	comp["Ex"] =1-tot;
-	
-	return comp;
 
-}
+#include "physics.h"
 
 matrix composition_map::X() const {
 
diff --git a/src/physics/eos_freeeos.cpp b/src/physics/eos_freeeos.cpp
index 87da2cbc..253c693a 100644
--- a/src/physics/eos_freeeos.cpp
+++ b/src/physics/eos_freeeos.cpp
@@ -21,7 +21,33 @@ extern "C" {
             double *entropy, int *iteration_count);
 }
 
-int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
+void set_eps(double* eps, const double_map &local_chemical_mix) {
+        // TODO: maybe use atomic masses defined in src/global/global.cpp
+        // Only isotops that exists in non negligible quantities have been considered
+        // according to their wikipedia page on the 3rd of July 2023
+        eps[0] = local_chemical_mix["H"] / 1.008e0;                     // H
+        eps[1] = local_chemical_mix["He3"] / 3. + local_chemical_mix["He4"] / 4.;     // He3 + He4
+        eps[2] = local_chemical_mix["C12"] / 12. + local_chemical_mix["C13"] / 13.;   // C12 + C13
+        eps[3] = local_chemical_mix["N14"] / 14. + local_chemical_mix["N15"] / 15.;   // N14 + N15
+        eps[4] = local_chemical_mix["O16"] / 16. + local_chemical_mix["O17"] / 17.;   // O16 + O17
+        eps[5] = local_chemical_mix["Ne20"] / 20. + local_chemical_mix["Ne22"] / 22.; // Ne20 + Ne22
+        eps[6] = local_chemical_mix["Na23"] / 23.; // Na23
+        eps[7] = local_chemical_mix["Mg24"] / 24 + local_chemical_mix["Mg25"] / 25 + local_chemical_mix["Mg26"] / 26; // Mg24 + Mg25 + Mg26
+        eps[8] = local_chemical_mix["Al27"] / 27.0; // Al27
+        eps[9] = local_chemical_mix["Si28"] / 28.0; // Si28
+        eps[10] = local_chemical_mix["P31"] / 31.0; // P31
+        eps[11] = local_chemical_mix["S32"] / 32.0; // S32
+        eps[12] = local_chemical_mix["Cl35"] / 35.0 + local_chemical_mix["Cl37"] / 37.0; // Cl35 + Cl37
+        eps[13] = local_chemical_mix["A40"] / 40.0; // A40
+        eps[14] = local_chemical_mix["Ca40"] / 40.0; // Ca40
+        eps[15] = local_chemical_mix["Ti"] / 47.867; // Ti //there are too many isotops
+        eps[16] = local_chemical_mix["Cr"] / 52.0; // Cr //there are too many isotops
+        eps[17] = local_chemical_mix["Mn55"] / 55.0; // Mn
+        eps[18] = local_chemical_mix["Fe"] / 55.845; // Fe //there are too many isotops
+        eps[19] = local_chemical_mix["Ni"] / 58.693; // Ni //there are too many isotops
+}
+
+int eos_freeeos(const composition_map &chemical_comp, const matrix &T, const matrix &p,
         matrix &rho, eos_struct &eos) {
 
     double t;
@@ -84,29 +110,7 @@ int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
     eos.prad.dim(T.nrows(), T.ncols()); // added MR june 2023
 
     for (int i=0; i<N; i++) {
-
-        double_map comp = initial_composition(X(i), Z);
-
-        eps[0] = comp["H"] / 1.008e0;                       // H
-        eps[1] = (comp["He3"] + comp["He4"]) / 4.0026e0;   // He3 + He4
-        eps[2] = (comp["C12"] + comp["C13"]) / 12.0111e0; // C12 + C13
-        eps[3] = (comp["N14"] + comp["N15"]) / 14.0067e0; // N14 + N15
-        eps[4] = (comp["O16"] + comp["O17"]) / 15.9994e0; // O16 + O17
-        eps[5] = 0.0; // Ne
-        eps[6] = 0.0; // Na
-        eps[7] = 0.0; // Mg
-        eps[8] = 0.0; // AL
-        eps[9] = 0.0; // Si
-        eps[10] = .0; // P
-        eps[11] = .0; // S
-        eps[12] = .0; // Cl
-        eps[13] = .0; // A
-        eps[14] = .0; // Ca
-        eps[15] = .0; // Ti
-        eps[16] = .0; // Cr
-        eps[17] = .0; // Mn
-        eps[18] = .0; // Fe
-        eps[19] = .0; // Ni
+        set_eps(eps, chemical_comp(i));
 
         double pi = p(i);
         double match_variable = log(pi);
@@ -123,12 +127,13 @@ int eos_freeeos(const matrix &X, double Z, const matrix &T, const matrix &p,
                 &iteration_count);
         rho(i) = rhoi;
         if (iteration_count < 0) {
+            // TODO: this makes no sense freeEOS has no tables...
             ester_err(
                     "Values outside freeEOS eos table:\n"
                     "  X = %e\n"
                     "  Z = %e\n"
                     "  T = %e\n"
-                    "  p = %e", X(i), Z, t, p(i));
+                    "  p = %e", chemical_comp.X()(i), chemical_comp.Z()(i), t, p(i));
         }
         eos.s(i) = entropy[0];
         eos.G1(i) = gamma1;
diff --git a/src/physics/eos_ideal.cpp b/src/physics/eos_ideal.cpp
index 48169903..d2976d74 100644
--- a/src/physics/eos_ideal.cpp
+++ b/src/physics/eos_ideal.cpp
@@ -5,12 +5,13 @@
 #include "physics.h"
 #include "constants.h"
 
-int eos_ideal(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos) {
+int eos_ideal(const composition_map &chemical_comp, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos) {
 		
 	matrix mu,b;
-	
-	mu=4/(3+5*X-Z);
+
+	// TODO: this mix between X and Z implicitly assume a certain Z mixture
+	mu = 4/(3 + 5*chemical_comp.X() - chemical_comp.Z());
 	eos.prad=zeros(T.nrows(),T.ncols());
 	rho=p/T/K_BOL*mu*HYDROGEN_MASS;
 	b=ones(T.nrows(),T.ncols());
diff --git a/src/physics/eos_idealrad.cpp b/src/physics/eos_idealrad.cpp
index 439021c5..be163bcb 100644
--- a/src/physics/eos_idealrad.cpp
+++ b/src/physics/eos_idealrad.cpp
@@ -4,12 +4,13 @@
 #include "physics.h"
 #include "constants.h"
 
-int eos_idealrad(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos) {
+int eos_idealrad(const composition_map &chemical_comp, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos) {
 
 	matrix mu,b;
 	
-	mu=4/(3+5*X-Z);
+	// TODO: this mix between X and Z implicitly assume a certain Z mixture
+	mu = 4/(3 + 5*chemical_comp.X() - chemical_comp.Z());
 	eos.prad=A_RAD/3*pow(T,4);
 	rho=(p-eos.prad)/T/K_BOL*mu*HYDROGEN_MASS;
 	b=(p-eos.prad)/p;
diff --git a/src/physics/eos_opal.cpp b/src/physics/eos_opal.cpp
index c36b3978..1da56aaa 100644
--- a/src/physics/eos_opal.cpp
+++ b/src/physics/eos_opal.cpp
@@ -14,33 +14,30 @@ extern"C" {
 	extern struct{
 		double esact,eos[10];
 	} eeos_;
-	extern struct{
-		int itime;
-	} lreadco_;
 }
 
-int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p,
-		matrix &rho,eos_struct &eos) {
-    
+static double Z_OPAL = -1.;
+
+int eos_opal(const composition_map &chemical_comp, const matrix &T, const matrix &p,
+		matrix &rho, eos_struct &eos) {
+
     matrix t6,p_mb;
-    int i,N,error=0;
-    
-    static double Z_table=-99;
+    int N;
+
+    if(Z_OPAL == -1.){
+        Z_OPAL = chemical_comp.Z()(0,0);
+        ester_warn(
+            "OPAL EOS doesn't support non-Z-Uniform stars (interpolation table computing is too long).\n"
+            "For **EOS ONLY** Z will be uniform.\n"
+            "Taking Z = Z(0,0) = %.6f",
+            Z_OPAL
+        );
+        zfs_interp_eos5_(&Z_OPAL);
+    }
 
     t6=T*1e-6;
     p_mb=p*1e-12;
 
-    if(Z!=Z_table) {
-    	lreadco_.itime=0;
-    	zfs_interp_eos5_(&Z);
-	    Z_table=Z;
-    }
-    
-    if(error) {
-    	printf("Can't initialize OPAL EOS table\n");
-    	return error;
-    }
-
     rho.dim(T.nrows(),T.ncols());
     eos.s.dim(T.nrows(),T.ncols());
     eos.G1.dim(T.nrows(),T.ncols());
@@ -56,10 +53,15 @@ int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p,
     
     N=T.nrows()*T.ncols();
 
-	double Xi,Zi,t6i,p_mbi,rhoi;
-    for(i=0;i<N;i++) {
-    	Xi=X(i);Zi=Z;t6i=t6(i);p_mbi=p_mb(i);
-    	eos5_xtrin_(&Xi,&Zi,&t6i,&p_mbi,&rhoi);
+    double Xi, t6i, p_mbi, rhoi;
+    for(int i = 0; i < N; i++) {
+        Xi = chemical_comp.X()(i);
+
+        t6i = t6(i);
+        p_mbi = p_mb(i);
+
+        // NOTE: Z_OPAL is used, not Z(i) see ester_warn above
+        eos5_xtrin_(&Xi,&Z_OPAL,&t6i,&p_mbi,&rhoi);
     	rho(i)=rhoi;
     	eos.s(i)=1e6*(*(eeos_.eos+2));
 		eos.G1(i)=*(eeos_.eos+7);
@@ -76,7 +78,7 @@ int eos_opal(const matrix &X,double Z,const matrix &T,const matrix &p,
                     "  X = %e\n"
                     "  Z = %e\n"
                     "  T = %e\n"
-                    "  p = %e", Xi, Zi, t6i, p_mbi);
+                    "  p = %e", Xi, Z_OPAL, t6i, p_mbi);
         }
     }
     if(exist(rho==-9e99)) {
diff --git a/src/physics/physics.cpp b/src/physics/physics.cpp
index 457a4246..3f4bc5fa 100644
--- a/src/physics/physics.cpp
+++ b/src/physics/physics.cpp
@@ -26,19 +26,19 @@ int opa_calc(const matrix &X,double Z,const matrix &T,const matrix &rho,
 	return error;
 }
 
-int eos_calc(const matrix &X,double Z,const matrix &T,const matrix &p,
-        matrix &rho,eos_struct &eos) {
+int eos_calc(const composition_map &chemical_comp, const matrix &T, const matrix &p,
+        matrix &rho, eos_struct &eos) {
 
     int error=0;
 
     if(!strcmp(eos.name,"ideal"))
-        error=eos_ideal(X,Z,T,p,rho,eos);
+        error = eos_ideal(chemical_comp, T, p, rho, eos);
     else if(!strcmp(eos.name,"ideal+rad"))
-        error=eos_idealrad(X,Z,T,p,rho,eos);
+        error = eos_idealrad(chemical_comp, T, p, rho, eos);
     else if(!strcmp(eos.name,"opal"))
-        error=eos_opal(X,Z,T,p,rho,eos);
+        error = eos_opal(chemical_comp, T, p, rho, eos);
     else if(!strcmp(eos.name,"freeeos"))
-        error = eos_freeeos(X, Z, T, p, rho, eos);
+        error = eos_freeeos(chemical_comp, T, p, rho, eos);
     else {
         ester_err("Unknown equation of state: %s",eos.name);
         return 1;
@@ -48,15 +48,15 @@ int eos_calc(const matrix &X,double Z,const matrix &T,const matrix &p,
 
 }
 
-int nuc_calc(const matrix_map &X,const matrix &T,const matrix &rho,
+int nuc_calc(const composition_map &comp, const matrix &T, const matrix &rho,
 		nuc_struct &nuc) {
 
 	int error=0;
 
 	if(!strcmp(nuc.name,"simple")) {
-		error=nuc_simple(X,T,rho,nuc);
+		error = nuc_simple(comp, T, rho, nuc);
 	} else if(!strcmp(nuc.name,"cesam")) {
-		error=nuc_cesam(X,T,rho,nuc);
+		error = nuc_cesam(comp, T, rho, nuc);
     } else {
         ester_err("Unknown nuc. reac. type: %s",nuc.name);
     	return 1;
@@ -66,13 +66,13 @@ int nuc_calc(const matrix_map &X,const matrix &T,const matrix &rho,
 
 }
 
-int atm_calc(const matrix &X,double Z,const matrix &g,const matrix &Teff,
-		const char *eos_name,const char *opa_name,atm_struct &atm) {
+int atm_calc(const composition_map &chemical_comp, const matrix &g, const matrix &Teff,
+        const char *eos_name, const char *opa_name, atm_struct &atm) {
 
 	int error=0;
 
 	if(!strcmp(atm.name,"onelayer")) {
-		error=atm_onelayer(X,Z,g,Teff,eos_name,opa_name,atm);
+		error = atm_onelayer(chemical_comp, g, Teff, eos_name, opa_name, atm);
     } else {
         ester_err("Unknown atmosphere type: %s", atm.name);
     	return 1;
diff --git a/src/star/star1d_solvers.cpp b/src/star/star1d_solvers.cpp
index c83b8ff1..8ca02754 100644
--- a/src/star/star1d_solvers.cpp
+++ b/src/star/star1d_solvers.cpp
@@ -7,11 +7,14 @@
 #include "symbolic.h"
 
 #include "matplotlib.h"
-//---------------------------------------------------------------------
-void star1d::fill() {
-    Y0=1.-X0-Z0;
-    init_comp();
 
+// TODO: apply the following move (at the end of the work)
+// The code between [BEGIN MOVE] and [END MODE] should be placed in star2d_class.cpp
+// it hasn't anything to do with solving, but is required by reading, and ester info for exmaple
+// whereas ester info doesn't require any other methods from star2d_solvers.cpp
+// [BEGIN MOVE]
+
+void star1d::fill() {
     eq_state();
 
     opacity();
@@ -33,7 +36,9 @@ void star1d::fill() {
     Omega=0;Omega_bk=0;Ekman=0;Omegac=0;
 
 }
-//---------------------------------------------------------------------
+
+// [END MOVE]
+
 solver *star1d::init_solver(int nvar_add) {
     int nvar;
     solver *op;
@@ -49,7 +54,8 @@ solver *star1d::init_solver(int nvar_add) {
 
     return op;
 }
-//---------------------------------------------------------------------
+
+
 void star1d::register_variables(solver *op) {
     int i,var_nr[ndomains];
 
@@ -86,12 +92,14 @@ void star1d::register_variables(solver *op) {
     op->regvar_dep("nuc.eps");
 
 }
-//---------------------------------------------------------------------
+
+
 double star1d::solve(solver *op) {
     matrix_map error_map;
     return solve(op, error_map, 0);
 }
-//---------------------------------------------------------------------
+
+
 double star1d::solve(solver *op, matrix_map& error_map, int nit) {
     int info[5];
     matrix rho0;
@@ -184,7 +192,8 @@ double star1d::solve(solver *op, matrix_map& error_map, int nit) {
     return err;
 
 }
-//---------------------------------------------------------------------
+
+
 void star1d::update_map(matrix dR) {
     if(ndomains==1) return;
 
@@ -201,7 +210,8 @@ void star1d::update_map(matrix dR) {
         map.R=R0+h*dR;
     }
 }
-//---------------------------------------------------------------------
+
+
 void star1d::solve_definitions(solver *op) {
 // substitution of the dep variables by their expression
 // with the basic variables
@@ -255,7 +265,8 @@ void star1d::solve_definitions(solver *op) {
     op->add_d("nuc.eps","log_Tc",nuc.dlneps_lnT*nuc.eps);
 
 }
-//---------------------------------------------------------------------
+
+
 void star1d::solve_poisson(solver *op) {
     int n,j0;
     matrix rhs;
@@ -303,7 +314,8 @@ void star1d::solve_poisson(solver *op) {
     }
     op->set_rhs("Phi",rhs);
 }
-//---------------------------------------------------------------------
+
+
 void star1d::solve_pressure(solver *op) {
     int n,j0;
     matrix rhs_p,rhs_pi_c;
@@ -339,7 +351,8 @@ void star1d::solve_pressure(solver *op) {
     op->set_rhs(eqn,rhs_p);
     op->set_rhs("pi_c",rhs_pi_c);
 }
-//---------------------------------------------------------------------
+
+
 void star1d::solve_temp(solver *op) {
     int n,j0,j1;
     matrix q;
@@ -548,7 +561,8 @@ void star1d::solve_temp(solver *op) {
     op->set_rhs("Lambda",rhs_Lambda);
 
 }
-//---------------------------------------------------------------------
+
+
 void star1d::solve_dim(solver *op) {
 // write the equations defining the global variables:
 // Radius, R, Tc, Pc and non-dim mass.
@@ -613,7 +627,8 @@ void star1d::solve_dim(solver *op) {
     op->set_rhs("log_R",rhs);
 
 }
-//---------------------------------------------------------------------
+
+
 void star1d::solve_map(solver *op) {
     int n,j0,j1;
     matrix rhs=zeros(ndomains,1);
@@ -661,7 +676,8 @@ void star1d::solve_map(solver *op) {
     }
     op->set_rhs("Ri",rhs);
 }
-//---------------------------------------------------------------------
+
+
 void star1d::solve_gsup(solver *op) {
 // Comfort equation
 // Code the equation gsup=4*pi*G*rho_c*R*m
@@ -677,7 +693,8 @@ void star1d::solve_gsup(solver *op) {
 
     op->set_rhs("gsup",zeros(1,1));
 }
-//---------------------------------------------------------------------
+
+
 void star1d::solve_Teff(solver *op) {
     matrix q,Te,F;
     int n=ndomains-1;
@@ -703,7 +720,8 @@ void star1d::solve_Teff(solver *op) {
 */
     op->set_rhs("Teff",zeros(1,1));
 }
-//---------------------------------------------------------------------
+
+
 void star1d::check_jacobian(solver *op,const char *eqn) {
     star1d B;
     matrix rhs,drhs,drhs2,qq;
diff --git a/src/star/star2d_solvers.cpp b/src/star/star2d_solvers.cpp
index 832aeab5..e4efde0b 100644
--- a/src/star/star2d_solvers.cpp
+++ b/src/star/star2d_solvers.cpp
@@ -4,16 +4,21 @@
 #include "star.h"
 #include <stdlib.h>
 #include <sys/time.h>
+#include <set>
 #include <string.h>
 #include "symbolic.h"
+#include "utils.h"
 
-/// \brief Initialize star's chemical composition, equation of state, opacity,
+// TODO: apply the following move (at the end of the work)
+// The code between [BEGIN MOVE] and [END MODE] should be placed in star2d_class.cpp
+// it hasn't anything to do with solving, but is required by reading, and ester info for exmaple
+// whereas ester info doesn't require any other methods from star2d_solvers.cpp
+// [BEGIN MOVE]
+
+/// \brief Initialize star's equation of state, opacity,
 /// nuclear reaction, mass definition, pi_c, Lambda, velocity, units, atmosphere,
 /// flatness, scaled keplerian angular velocity
 void star2d::fill() {
-	Y0=1.-X0-Z0;
-	init_comp();
-
 	eq_state();
 	opacity();
 	nuclear();
@@ -36,56 +41,97 @@ void star2d::fill() {
 
 }
 
+
+void star2d::init_metal_mix() {
+	// a security to avoid initializing multiple time the metal mix
+	// it should be done only once, then m_metal_mix should be used
+	if(!m_metal_mix.empty()){
+		ester_err("init_metal_mix must be called only once");
+	}
+
+	file_parser fp;
+
+	// TODO: change this hardcoded path
+	char file[] = "metal-mix.cfg";
+
+	char* arg = NULL;
+	char* val = NULL;
+	std::set<std::string> metals = {"C12","C13","N14","N15","O16","O17","Ne20","Ne22","Na23",
+									"Mg24","Mg25","Mg26","Al27","Si28","P31","S32","Cl35","Cl37",
+									"A40","Ca40","Ti","Cr","Mn55","Fe","Ni"};
+	// Initialization of m_metal_mix
+	for(std::string metal: metals){
+		m_metal_mix[metal] = .0;
+	}
+	double metals_fraction = .0;
+
+	if(!fp.open(file)){
+		printf("Can't open configuration file %s\n", file);
+		perror("Error:");
+		exit(1);
+	} else {
+		int line;
+		while(line = fp.get(arg,val)) {
+			if(val == NULL){
+				printf("Syntax error in configuration file %s, line %d\n", file, line);
+				exit(1);
+			}
+			if(metals.find(arg) == metals.end()){
+				// the metal specified in the config file isn't supported
+				printf("%s is unknown, possible metals composition are: ", val);
+				for(std::string metal: metals){
+					printf("%s ", metal);
+				}
+				puts("\n");
+				exit(1);
+			}
+			metals_fraction += atof(val);
+			m_metal_mix[arg] = Z0 * atof(val);
+		}
+	}
+	fp.close();
+
+	// will be removed:
+	printf("A config file %s has been used to config metal composition\n", file);
+	printf("The sum of metals fractions (of Z sum of metal mass fraction) is %f, and should be as close to 1 as possible.\n", metals_fraction);
+}
+
+
 void star2d::init_comp() {
-// Update the object comp
+	Y0 = 1. - (X0 + Z0);
+	m_He_isotopic_ratio = 3.15247417638132e-04;
+
+	init_metal_mix();
 
-	comp=initial_composition(X0,Z0)*ones(nr,nth);
+	// make a copy of m_metal_mix
+	double_map temp_chemical_mix = double_map(m_metal_mix);
+	// Note: A global chemical_mix (analog to metal_mix) would make no sense
+	// because H and He quantities can be different in different star's layers
+	// that's not the case (at the time with a static version of ESTER) of the metal mixture
+
+	temp_chemical_mix["H"] = X0;
+	// TODO: should the following be hardcoded?
+	temp_chemical_mix["He3"] = m_He_isotopic_ratio * Y0;
+	temp_chemical_mix["He4"] = Y0 - temp_chemical_mix["He3"];
+	// Note that the previous way of setting He3 He4 enforce He3+He4 = Y
+
+	// Init the object comp
+	comp = temp_chemical_mix*ones(nr,nth);
 
 	if(!conv) return;
 
-// Count the number of point in the core:
-    int n = 0;
-    for (int i=0; i<conv; i++) {
-        n += map.gl.npts[i];
-    }
-// and put some fraction Xc of X_envelope in the core 
-    comp.setblock(0,n-1,0,-1,initial_composition(Xc*X0,Z0)*ones(n,nth));
-
-// Put a special composition in the last domain: Daniel's request
-/*	printf("WARNING : Daniel's superstar computed!\n");
-
-    n = 0; // Count the number of points to the before before last domain
-    for (int i=0; i<ndomains-1; i++) {
-        n += map.gl.npts[i];
-    }
-	printf(" in init_comp n=%d, nr=%d\n",n,nr);
-//	double x=0.70,z=0.02;
-	double x=X0*0.10,z=0.062;
-        comp.setblock(0,n-1,0,-1,initial_composition(x,z)*ones(n,nth));
-	comp.setblock(n,nr-1,0,-1,initial_composition(X0,Z0)*ones(nr-n,nth));
-	printf(" in init_comp comp.X=%e, \n",comp.X()(0,0));
-*/
-
-/* Delphine's attempt
-    if(stratified_comp == 0) {
-        comp.setblock(0,n-1,0,-1,initial_composition(Xc*X0,Z0)*ones(n,nth));
-    }
-    else {
-        comp.setblock(0, n-1, 0, -1,
-                initial_composition(Xc*X0, Z0)*ones(n, nth));
-        int m = 0;
-        int l = n;
-        double a = (1.-exp((1.-(nr/n))))/(X0*(1.-Xc));
-        for(int i=conv+1; i<ndomains; i++) {
-            m = map.gl.npts[i];
-            comp.setblock(l, l+m-1, 0, -1,
-                    initial_composition(((Xc*X0)+((1./a)*(1.-exp((1.-(l/n)))))),Z0)*ones(m,nth));
-            l += m;
-        }
-    }
-*/
+	// Count the number of point in the core:
+	int n = 0;
+	for (int i=0; i<conv; i++) {
+		n += map.gl.npts[i];
+	}
+	// and put some fraction Xc of X_envelope in the core
+	temp_chemical_mix["H"] = Xc * X0;
+	comp.setblock(0, n-1, 0, -1, (temp_chemical_mix*ones(n,nth)));
 }
 
+// [END MOVE]
+
 void star2d::calc_veloc() {
 // vr=rz*V^zeta+rt*V^theta,  vt=r*V^theta
 	//vr=(G,map.leg.D_11)/r+(map.rt/r+cos(th)/sin(th))/r*G;
@@ -98,6 +144,7 @@ void star2d::calc_veloc() {
 	vt/=rho;
 }
 
+
 solver *star2d::init_solver(int nvar_add) {
 	int nvar;
 	solver *op;
@@ -113,6 +160,7 @@ solver *star2d::init_solver(int nvar_add) {
 	return op;
 }
 
+
 void star2d::register_variables(solver *op) {
 	int i,var_nr[ndomains+1];
 
@@ -158,11 +206,13 @@ void star2d::register_variables(solver *op) {
 
 }
 
+
 double star2d::solve(solver *op) {
     matrix_map error_map;
     return solve(op, error_map, 0);
 }
 
+
 /// \brief Performs one step of the Newton algorithm to compute the star's
 /// internal structure.
 double star2d::solve(solver *op, matrix_map& error_map, int nit) {
@@ -281,9 +331,11 @@ double star2d::solve(solver *op, matrix_map& error_map, int nit) {
 
 }
 
+
 // Special treatment for updating R_i because we need
 // a different relaxation parameter "h".
 
+
 void star2d::update_map(matrix dR) {
 	double h=1,dmax=config.newton_dmax;
 
@@ -299,6 +351,7 @@ void star2d::update_map(matrix dR) {
 
 }
 
+
 /// \brief insert the definitions depending on opacity and eos tables into the solver,
 /// and the definitions used by the mapping (eta,deta,Ri,dRi,...), and the entropy
 void star2d::solve_definitions(solver *op) {
@@ -919,7 +972,6 @@ void star2d::solve_temp(solver *op) {
 
 
 /// \brief Writes the equations for the dimensional quantities (T_c, rho_c, R, etc.)
-
 void star2d::solve_dim(solver *op) {
 	int n,j0,j1;
 	matrix q,rhs;
@@ -1109,8 +1161,6 @@ void star2d::solve_map(solver *op) {
 }
 
 
-
-
 /// \brief Equation setting the equatorial angular velocity
 void star2d::solve_Omega(solver *op) {
 	int n;
@@ -1135,7 +1185,6 @@ void star2d::solve_Omega(solver *op) {
 }
 
 
-
 /// \brief Equation giving the effective surface gravity gsup
 /// gsup=(-\vn\cdot\grad P)/rho
 void star2d::solve_gsup(solver *op) {
@@ -1172,6 +1221,7 @@ void star2d::solve_gsup(solver *op) {
 
 }
 
+
 /// \brief Equation setting the surface effective temperature
 /// Derived from sigma T_e^4 = -xi\vn\cdot\gradT
 void star2d::solve_Teff(solver *op) {
@@ -1209,6 +1259,7 @@ void star2d::solve_Teff(solver *op) {
 
 }
 
+
 /// \brief Equation setting the 'simple' atmosphere model equations
 /// 'simple' == the polytropic model for the atmosphere
 /// checked but approximate on ... "p",-1/p.row(-1)/(n_atm+1) ...
diff --git a/src/star/star_phys.cpp b/src/star/star_phys.cpp
index 0a9153fc..7d1435c8 100644
--- a/src/star/star_phys.cpp
+++ b/src/star/star_phys.cpp
@@ -35,10 +35,12 @@ void star2d::eq_state() {
 	int error;
 
 	matrix rhoc_m(1,1);
-	eos_calc(comp.X()(0,0)*ones(1,1),Z0,ones(1,1)*Tc,ones(1,1)*pc,rhoc_m,eos);
-	rhoc=rhoc_m(0);
+	// TODO: is the following really usefull???
+	// couldn't we just get rhoc from rho(0) ??? instead of this weird call only computing central value
+	eos_calc(comp(0,0)*ones(1,1), Tc*ones(1,1), pc*ones(1,1), rhoc_m, eos);
+	rhoc = rhoc_m(0);
 
-	error=eos_calc(comp.X(),Z0,T*Tc,p*pc,rho,eos);
+	error = eos_calc(comp, T*Tc, p*pc, rho, eos);
 
 
 //	rhoc=rho(0);
@@ -61,7 +63,7 @@ void star2d::atmosphere() {
 	}
 
 
-	error=atm_calc(comp.X(),Z0,gsup(),Teff(),eos.name,opa.name,atm);
+	error = atm_calc(comp, gsup(), Teff(), eos.name, opa.name, atm);
 
 	if(error) exit(1);
 
diff --git a/utils/starR/star1dR_class.cpp b/utils/starR/star1dR_class.cpp
index 91962f8c..9b1243ed 100644
--- a/utils/starR/star1dR_class.cpp
+++ b/utils/starR/star1dR_class.cpp
@@ -58,7 +58,6 @@ double star1dR::solve(solver *op) {
 
 void star1dR::fill() {
 
-	init_comp();
 	eq_state();
 	m=4*PI*(map.gl.I,rho*r*r)(0);
 	M=m*rhoc*R*R*R;
@@ -112,7 +111,6 @@ void star1dR::solve_definitions(solver *op) {
 	double dXc=1e-8;
 	
 	Xc+=dXc;
-	init_comp();
 	nuclear();
 	opacity();
 	eq_state();
diff --git a/utils/starR/star2dR_class.cpp b/utils/starR/star2dR_class.cpp
index d3c551ac..db5fa3b7 100644
--- a/utils/starR/star2dR_class.cpp
+++ b/utils/starR/star2dR_class.cpp
@@ -65,7 +65,6 @@ double star2dR::solve(solver *op) {
 
 void star2dR::fill() {
 
-	init_comp();
 	eq_state();
 	m=2*PI*(map.gl.I,rho*r*r*map.rz,map.leg.I_00)(0);
 	M=m*rhoc*R*R*R;
@@ -124,7 +123,6 @@ void star2dR::solve_definitions(solver *op) {
 	double dXc=1e-8;
 
 	Xc+=dXc;
-	init_comp();
 	nuclear();
 	opacity();
 	eq_state();