diff --git a/binarycpython/utils/distribution_functions.py b/binarycpython/utils/distribution_functions.py
index e7a7891165269e7e163b3a4666016d05c81b850a..2f0ad116667454609ccbadc9129ebea849849532 100644
--- a/binarycpython/utils/distribution_functions.py
+++ b/binarycpython/utils/distribution_functions.py
@@ -1125,7 +1125,7 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
 
     if options["multiplicity_model"] == "Poisson":
         multiplicity = Moecache["rinterpolator_multiplicity"].interpolate(
-            [np.log10(options["M1"])]
+            [np.log10(options["M_1"])]
         )[0]
 
         # Fill the multiplicity array
@@ -1151,7 +1151,7 @@ def Moe_de_Stefano_2017_multiplicity_fractions(options, verbosity=0):
         for n in range(3):
             full_fractions_array[n] = (
                 Moecache["rinterpolator_multiplicity"].interpolate(
-                    [np.log10(options["M1"])]
+                    [np.log10(options["M_1"])]
                 )[n + 1])
 
         # Set last value
@@ -1272,7 +1272,7 @@ def build_q_table(options, m, p, verbosity=0):
         #     0
         # ]  # TODO: this lower range must not be lower than Mmin.
         # print("build_q_table qmin: {}".format(qmin))
-        qmin = options['Mmin']/options['M1']
+        qmin = options['Mmin']/options['M_1']
         print("build_q_table qmin: {}".format(qmin))
 
         # qmax = maximum_mass_ratio_for_RLOF(options[m], options[p])
@@ -1817,7 +1817,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
     prob_dict['multiplicity'] = multiplicity_probability
     verbose_print(
         "\tM&S: Moe_de_Stefano_2017_pdf: Appended multiplicity (mass1 = {}) probability ({}) to the prob dict ({})".format(
-            options["M1"], prob_dict['multiplicity'], prob_dict
+            options["M_1"], prob_dict['multiplicity'], prob_dict
         ),
         verbosity,
         _MS_VERBOSITY_LEVEL,
@@ -1831,11 +1831,11 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
 
     # TODO: Create an n-part-powerlaw method that can have breakpoints and slopes. I'm using a three-part powerlaw now.
     # TODO: is this actually the correct way? putting the M1 in there? Do we sample in logspace?
-    M1_probability = Kroupa2001(options["M1"]) * options["M1"]
-    prob_dict['M1'] = M1_probability
+    M1_probability = Kroupa2001(options["M_1"]) * options["M_1"]
+    prob_dict['M_1'] = M1_probability
     verbose_print(
         "\tM&S: Moe_de_Stefano_2017_pdf: Appended Mass (m={}) probability ({}) to the prob dict ({})".format(
-            options["M1"], prob_dict['M1'], prob_dict
+            options["M_1"], prob_dict['M_1'], prob_dict
         ),
         verbosity,
         _MS_VERBOSITY_LEVEL,
@@ -1870,7 +1870,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
             )
 
         # Make a table storing Moe's data for q distributions
-        if options.get("M2", None) or options.get("M3", None) or options.get("M4", None):
+        if options.get("M_2", None) or options.get("M_3", None) or options.get("M_4", None):
             if not Moecache.get("rinterpolator_q", None):
                 Moecache["rinterpolator_q"] = py_rinterpolate.Rinterpolate(
                     table=Moecache["q_distributions"],  # Contains the table of data
@@ -1907,30 +1907,30 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
         ###############
         # Calculation for period of the binary
 
-        if options.get("M2", None):
+        if options.get("M_2", None):
             # Separation of the inner binary
             options["sep"] = calc_sep_from_period(
-                options["M1"], options["M2"], options["P"]
+                options["M_1"], options["M_2"], options["P"]
             )
             # TODO: add check for min_P with instant RLOF?
             # TODO: Actually use the value above.
             # Total mass inner binary:
-            options["M1+M2"] = options["M1"] + options["M2"]
+            options["M_1+M_2"] = options["M_1"] + options["M_2"]
 
         # Calculate P integral or use cached value
         calc_P_integral(
-            options, 1, "P_integrals", "rinterpolator_log10P", "M1", verbosity
+            options, 1, "P_integrals", "rinterpolator_log10P", "M_1", verbosity
         )
 
         # Set probabilty for P1
         p_val = Moecache["rinterpolator_log10P"].interpolate(
-            [np.log10(options["M1"]), np.log10(options["P"])]
+            [np.log10(options["M_1"]), np.log10(options["P"])]
         )[0]
-        p_val = p_val / Moecache["P_integrals"][options["M1"]]
+        p_val = p_val / Moecache["P_integrals"][options["M_1"]]
         prob_dict['P'] = p_val
         verbose_print(
             "\tM&S: Moe_de_Stefano_2017_pdf: Appended period (m={}, P={}) probability ({}) to the prob list ({})".format(
-                options["M1"], options["P"], prob_dict['P'], prob_dict
+                options["M_1"], options["P"], prob_dict['P'], prob_dict
             ),
             verbosity,
             _MS_VERBOSITY_LEVEL,
@@ -1945,11 +1945,11 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
         # we need to construct the q table for the given M1
         # subject to qmin = Mmin/M1
 
-        if options.get("M2", None):
+        if options.get("M_2", None):
             # Build the table for q
-            primary_mass = options["M1"]
-            secondary_mass = options["M2"]
-            m_label = "M1"
+            primary_mass = options["M_1"]
+            secondary_mass = options["M_2"]
+            m_label = "M_1"
             p_label = "P"
 
             # Construct the q table
@@ -1969,9 +1969,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
             prob_dict['q'] = q_prob
             verbose_print(
                 "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M={} P={} q={}) probability ({}) to the prob list ({}) ".format(
-                    options["M1"],
+                    options["M_1"],
                     options["P"],
-                    options["M2"] / options["M1"],
+                    options["M_2"] / options["M_1"],
                     prob_dict['q'],
                     prob_dict,
                 ),
@@ -1991,19 +1991,19 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
             if not options['ecc']==None:
                 # Calculate ecc integral or use cached value
                 calc_e_integral(
-                    options, "ecc_integrals", "rinterpolator_e", "M1", "P", verbosity
+                    options, "ecc_integrals", "rinterpolator_e", "M_1", "P", verbosity
                 )
-                mass_period_string = "{}_{}".format(options["M1"], options["P"])
+                mass_period_string = "{}_{}".format(options["M_1"], options["P"])
 
                 # Set probabilty for ecc
                 ecc_val = Moecache["rinterpolator_e"].interpolate(
-                    [np.log10(options["M1"]), np.log10(options["P"]), options["ecc"]]
+                    [np.log10(options["M_1"]), np.log10(options["P"]), options["ecc"]]
                 )[0]
                 ecc_val = ecc_val / Moecache["ecc_integrals"][mass_period_string]
                 prob_dict['ecc'] = ecc_val
                 verbose_print(
                     "\tM&S: Moe_de_Stefano_2017_pdf: Appended eccentricity (m={}, P={}, ecc={}) probability ({}) to the prob list ({})".format(
-                        options["M1"], options["P"], options["ecc"], prob_dict['ecc'], prob_dict
+                        options["M_1"], options["P"], options["ecc"], prob_dict['ecc'], prob_dict
                     ),
                     verbosity,
                     _MS_VERBOSITY_LEVEL,
@@ -2029,7 +2029,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
 
             # TODO: Is this a correct assumption?
             max_sep = 10.0 * options["sep"] * (1.0 + options["ecc"])
-            min_P2 = calc_period_from_sep(options["M1+M2"], options["mmin"], max_sep)
+            min_P2 = calc_period_from_sep(options["M_1+M_2"], options["mmin"], max_sep)
 
             if options["P2"] < min_P2:
                 # period is too short : system is not hierarchical
@@ -2050,7 +2050,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
                 # hence the separation between the outer star
                 # and inner binary
                 options["sep2"] = calc_sep_from_period(
-                    options["M3"], options["M1+M2"], options["P2"]
+                    options["M_3"], options["M_1+M_2"], options["P2"]
                 )
 
                 # Check for cached value of P integral or calculate
@@ -2059,19 +2059,19 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
                     min_P2,
                     "P2_integrals",
                     "rinterpolator_log10P",
-                    "M1+M2",
+                    "M_1+M_2",
                     verbosity,
                 )
 
                 # Add the probability
                 p_val = Moecache["rinterpolator_log10P"].interpolate(
-                    [np.log10(options["M1+M2"]), np.log10(options["P2"])]
+                    [np.log10(options["M_1+M_2"]), np.log10(options["P2"])]
                 )[0]
-                p_val = p_val / Moecache["P2_integrals"][options["M1+M2"]]
+                p_val = p_val / Moecache["P2_integrals"][options["M_1+M_2"]]
                 prob_dict['P2'] = p_val
                 verbose_print(
                     "\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (m1={} m2={}, P2={}) probability ({}) to the prob list ({})".format(
-                        options["M1"], options["M2"], options["P2"], prob_dict['P2'], prob_dict
+                        options["M_1"], options["M_2"], options["P2"], prob_dict['P2'], prob_dict
                     ),
                     verbosity,
                     _MS_VERBOSITY_LEVEL,
@@ -2087,9 +2087,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
                 # subject to qmin = Mmin/(M1+M2)
 
                 # Set the variables for the masses and their names
-                primary_mass = options["M1+M2"]
-                secondary_mass = options["M3"]
-                m_label = "M1+M2"
+                primary_mass = options["M_1+M_2"]
+                secondary_mass = options["M_3"]
+                m_label = "M_1+M_2"
                 p_label = "P2"
 
                 # Build q table
@@ -2106,9 +2106,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
                 ].interpolate([secondary_mass / primary_mass])[0]
                 prob_dict['q2'] = q2_val
                 verbose_print(
-                    "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M1+M2={} M3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
-                        options["M1+M2"],
-                        options["M3"],
+                    "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M_1+M_2={} M_3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
+                        options["M_1+M_2"],
+                        options["M_3"],
                         options["P"],
                         secondary_mass / primary_mass,
                         prob_dict['q2'],
@@ -2136,7 +2136,7 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
                     # TODO: fix this here
                     max_sep3 = 0.2 * options["sep2"] * (1.0 + options["ecc2"])
                     max_per3 = calc_period_from_sep(
-                        options["M1+M2"], options["mmin"], max_sep3
+                        options["M_1+M_2"], options["mmin"], max_sep3
                     )
 
                     # Calculate P integral or use the cached value
@@ -2146,19 +2146,19 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
                         min_P2,
                         "P2_integrals",
                         "rinterpolator_log10P",
-                        "M1+M2",
+                        "M_1+M_2",
                         verbosity,
                     )
 
                     # Set probability
                     p_val = Moecache["rinterpolator_log10P"].interpolate(
-                        [np.log10(options["M1+M2"]), np.log10(options["P2"])]
+                        [np.log10(options["M_1+M_2"]), np.log10(options["P2"])]
                     )[0]
-                    p_val = p_val / Moecache["P2_integrals"][options["M1+M2"]]
+                    p_val = p_val / Moecache["P2_integrals"][options["M_1+M_2"]]
                     prob_dict['P3'] = p_val
                     verbose_print(
-                        "\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (M=4) (m1={} m2={}, P2={}) probability ({}) to the prob list ({})".format(
-                            options["M1"], options["M2"], options["P2"], prob_dict['P3'], prob_dict
+                        "\tM&S: Moe_de_Stefano_2017_pdf: Appended period2 (M=4) (M_1={} M_2={}, P2={}) probability ({}) to the prob list ({})".format(
+                            options["M_1"], options["M_2"], options["P2"], prob_dict['P3'], prob_dict
                         ),
                         verbosity,
                         _MS_VERBOSITY_LEVEL,
@@ -2175,9 +2175,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
                     # Make a table storing Moe's data for q distributions
 
                     # Build the table for q2
-                    primary_mass = options["M1+M2"]
-                    secondary_mass = options["M3"]
-                    m_label = "M1+M2"
+                    primary_mass = options["M_1+M_2"]
+                    secondary_mass = options["M_3"]
+                    m_label = "M_1+M_2"
                     p_label = "P2"
 
                     # Calculate new q table
@@ -2194,9 +2194,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
                     ].interpolate([secondary_mass / primary_mass])[0]
                     prob_dict['q3'] = q3_prob
                     verbose_print(
-                        "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M1+M2={} M3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
-                            options["M1+M2"],
-                            options["M3"],
+                        "\tM&S: Moe_de_Stefano_2017_pdf: appended mass ratio (M_1+M_2={} M_3={} P={} q={}) probability ({}) to the prob list ({}) ".format(
+                            options["M_1+M_2"],
+                            options["M_3"],
                             options["P"],
                             secondary_mass / primary_mass,
                             prob_dict['q3'],
@@ -2229,8 +2229,8 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
     # Some info
     if multiplicity == 1:
         verbose_print(
-            "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q=N/A log10P=N/A ({}): {} -> {}\n".format(
-                options["M1"],
+            "\tM&S: Moe_de_Stefano_2017_pdf: M_1={} q=N/A log10P=N/A ({}): {} -> {}\n".format(
+                options["M_1"],
                 len(prob_dict),
                 str(prob_dict),
                 prob_dict['total_probdens']
@@ -2240,9 +2240,9 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
         )
     elif multiplicity == 2:
         verbose_print(
-            "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ecc={} ({}): {} -> {}\n".format(
-                options["M1"],
-                options["M2"] / options["M1"] if options.get("M2", None) else "N/A",
+            "\tM&S: Moe_de_Stefano_2017_pdf: M_1={} q={} log10P={} ecc={} ({}): {} -> {}\n".format(
+                options["M_1"],
+                options["M_2"] / options["M_1"] if options.get("M_2", None) else "N/A",
                 np.log10(options["P"]),
                 options['ecc'] if options.get("ecc", None) else "N/A",
                 len(prob_dict),
@@ -2254,12 +2254,12 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
         )
     elif multiplicity == 3:
         verbose_print(
-            "\tM&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ecc={} M3={} log10P2={} ecc2={} ({}): {} -> {}".format(
-                options["M1"],
-                options["M2"] / options["M1"] if options.get("M2", None) else "N/A",
+            "\tM&S: Moe_de_Stefano_2017_pdf: M_1={} q={} log10P={} ecc={} M_3={} log10P2={} ecc2={} ({}): {} -> {}".format(
+                options["M_1"],
+                options["M_2"] / options["M_1"] if options.get("M_2", None) else "N/A",
                 np.log10(options["P"]),
                 options["ecc"] if options.get("ecc", None) else "N/A",
-                options["M3"],
+                options["M_3"],
                 np.log10(options["P2"]),
                 options["ecc2"] if options.get("ecc2", None) else "N/A",
                 len(prob_dict),
@@ -2271,15 +2271,15 @@ def Moe_de_Stefano_2017_pdf(options, verbosity=0):
         )
     elif multiplicity == 4:
         verbose_print(
-            "M&S: Moe_de_Stefano_2017_pdf: M1={} q={} log10P={} ecc={} M3={} log10P2={} ecc2={} M4={} log10P3={} ecc3={} ({}) : {} -> {}".format(
-                options["M1"],
-                options["M2"] / options["M1"] if options.get("M2", None) else "N/A",
+            "M&S: Moe_de_Stefano_2017_pdf: M_1={} q={} log10P={} ecc={} M_3={} log10P2={} ecc2={} M_4={} log10P3={} ecc3={} ({}) : {} -> {}".format(
+                options["M_1"],
+                options["M_2"] / options["M_1"] if options.get("M_2", None) else "N/A",
                 np.log10(options["P"]),
                 options['ecc'] if options.get("ecc", None) else "N/A",
-                options["M3"],
+                options["M_3"],
                 np.log10(options["P2"]),
                 options["ecc2"] if options.get("ecc2", None) else "N/A",
-                options["M4"],
+                options["M_4"],
                 np.log10(options["P3"]),
                 options["ecc3"] if options.get("ecc3", None) else "N/A",
                 len(prob_dict), str(prob_dict), prob_dict['total_probdens']
diff --git a/binarycpython/utils/functions.py b/binarycpython/utils/functions.py
index 310ff4bb64767fd62780a0fa8ca7914ac5e24d0a..8c3d9d2866fe7b845da2465c90320e10ffce4b8f 100644
--- a/binarycpython/utils/functions.py
+++ b/binarycpython/utils/functions.py
@@ -15,6 +15,7 @@ import copy
 import inspect
 import sys
 import time
+import types
 
 from io import StringIO
 from typing import Union, Any
@@ -170,17 +171,17 @@ def subtract_dicts(dict_1: dict, dict_2: dict) -> dict:
                     del new_dict[key]
 
             else:
-                print(
-                    "Error key: {} value: {} type: {} and key: {} value: {} type: {} are not of the same type and cannot be merged".format(
-                        key,
-                        dict_1[key],
-                        type(dict_1[key]),
-                        key,
-                        dict_2[key],
-                        type(dict_2[key]),
-                    )
+                msg = "Error key: {} value: {} type: {} and key: {} value: {} type: {} are not of the same type and cannot be merged".format(
+                    key,
+                    dict_1[key],
+                    type(dict_1[key]),
+                    key,
+                    dict_2[key],
+                    type(dict_2[key]),
                 )
-                raise ValueError
+
+                print(msg)
+                raise ValueError(msg)
 
         # This is where the keys are the same
         else:
@@ -211,26 +212,37 @@ def subtract_dicts(dict_1: dict, dict_2: dict) -> dict:
     return new_dict
 
 
-def get_moe_distefano_dataset(options):
+def get_moe_distefano_dataset(options, verbosity=0):
     """
     Function to get the default moe and Distefano dataset or accept a userinput.
     """
 
     if not options.get("file", None):
-        print("Using the default Moe and de Stefano 2017 datafile")
+        verbose_print(
+            "Using the default Moe and de Stefano 2017 datafile",
+            verbosity,
+            1,
+        )
 
         json_data = copy.deepcopy(moe_distefano_data.moe_distefano_2017_data)
 
     else:
         if not os.path.isfile(options["file"]):
-            print(
+            verbose_print(
                 "The provided 'file' Moe and de Stefano JSON file does not seem to exist at {}".format(
                     options["file"]
-                )
+                ),
+                verbosity,
+                1,
             )
+
             raise ValueError
         if not options["file"].endswith(".json"):
-            print("Provided filename does not end with .json")
+            verbose_print(
+                "Provided filename is not a json file",
+                verbosity,
+                1,
+            )
 
         else:
             # Read input data and Clean up the data if there are whitespaces around the keys
diff --git a/binarycpython/utils/grid.py b/binarycpython/utils/grid.py
index 454f7e6e03eec6d8915ed6889f61114d9872137b..2317fbddcaa5abcf38e9631f0326341e9cfa4a9e 100644
--- a/binarycpython/utils/grid.py
+++ b/binarycpython/utils/grid.py
@@ -3264,9 +3264,9 @@ class Population:
         """
 
         # Only if the grid is loaded and moecache contains information
-        if (not self.grid_options['_loaded_ms_data']) and Moecache:
+        if (not self.grid_options['_loaded_ms_data']) and not Moecache:
             # Load the json data
-            json_data = get_moe_distefano_dataset(self.grid_options['m&s_options'])
+            json_data = get_moe_distefano_dataset(self.grid_options['m&s_options'], verbosity=self.grid_options['verbosity'])
 
             # entry of log10M1 is a list containing 1 dict. We can take the dict out of the list
             json_data["log10M1"] = json_data["log10M1"][0]
@@ -3571,9 +3571,9 @@ class Population:
             ],
             gridtype="centred",
             dphasevol="dlnm1",
-            precode='M_1 = np.exp(lnm1); options["M1"]=M_1',
+            precode='M_1 = np.exp(lnm1); options["M_1"]=M_1',
             probdist="Moe_de_Stefano_2017_pdf({{{}, {}, {}}}, verbosity=self.grid_options['verbosity'])['total_probdens'] if multiplicity == 1 else 1".format(
-                str(dict(self.grid_options['m&s_options']))[1:-1], "'multiplicity': multiplicity", "'M1': M_1"
+                str(dict(self.grid_options['m&s_options']))[1:-1], "'multiplicity': multiplicity", "'M_1': M_1"
             ),
         )
 
@@ -3690,7 +3690,6 @@ q2max=maximum_mass_ratio_for_RLOF(M_1+M_2, orbital_period_triple)
 
                 # Triples: mass ratio
                 # Note, the mass ratio is M_outer/M_inner
-                print("M2 resolution {}".format(self.grid_options['m&s_options']["resolutions"]["M"][2]))
                 self.add_grid_variable(
                     name="q2",
                     parameter_name="M_3",
@@ -3775,7 +3774,6 @@ q3max=maximum_mass_ratio_for_RLOF(M_3, orbital_period_quadruple)
                     )
 
                     # Quadruple: mass ratio : M_outer / M_inner
-                    print("M3 resolution {}".format(self.grid_options['m&s_options']["resolutions"]["M"][3]))
                     self.add_grid_variable(
                         name="q3",
                         parameter_name="M_4",
@@ -3783,7 +3781,7 @@ q3max=maximum_mass_ratio_for_RLOF(M_3, orbital_period_quadruple)
                         valuerange=[
                             self.grid_options['m&s_options']["ranges"]["q"][0]
                             if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
-                            else "options.get('Mmin', 0.07)/(m3)",
+                            else "options.get('Mmin', 0.07)/(M_3)",
                             self.grid_options['m&s_options']["ranges"]["q"][1]
                             if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
                             else "q3max",
@@ -3796,11 +3794,11 @@ q3max=maximum_mass_ratio_for_RLOF(M_3, orbital_period_quadruple)
 M_4 = q3 * M_3
 sep3 = calc_sep_from_period((M_3), M_4, orbital_period_quadruple)
 eccentricity3=0
-        """,
+""",
                         spacingfunc="const({}, {}, {})".format(
                             self.grid_options['m&s_options']["ranges"]["q"][0]
                             if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
-                            else "options.get('Mmin', 0.07)/(m3)",
+                            else "options.get('Mmin', 0.07)/(M_3)",
                             self.grid_options['m&s_options']["ranges"]["q"][1]
                             if self.grid_options['m&s_options'].get("ranges", {}).get("q", None)
                             else "q3max",
@@ -3844,10 +3842,10 @@ eccentricity3=0
         updated_options = "{{{}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}, {}}}".format(
             str(dict(self.grid_options['m&s_options']))[1:-1],
             '"multiplicity": multiplicity',
-            '"M1": M_1',
-            '"M2": M_2',
-            '"M3": M_3',
-            '"M4": M_4',
+            '"M_1": M_1',
+            '"M_2": M_2',
+            '"M_3": M_3',
+            '"M_4": M_4',
             '"P": orbital_period',
             '"P2": orbital_period_triple',
             '"P3": orbital_period_quadruple',
@@ -3928,7 +3926,6 @@ eccentricity3=0
         grid_options['multiplicity_fraction_function'] will be checked for the choice
         """
 
-
         # Just return 1 if no option has been chosen
         if self.grid_options['multiplicity_fraction_function'] == 0:
             print("Chosen not to use any multiplicity fraction.")
@@ -3942,15 +3939,23 @@ eccentricity3=0
         # Go over the chosen options
         if self.grid_options['multiplicity_fraction_function'] == 1:
             # Arenou 2010 will be used
-            # TODO: add verbose logging
-            print("Using Arenou 2010 to calculate multiplicity fractions")
+            verbose_print(
+                "_calculate_multiplicity_fraction: Using Arenou 2010 to calculate multiplicity fractions",
+                self.grid_options["verbosity"],
+                1,
+            )
+
             binary_fraction = Arenou2010_binary_fraction(system_dict['M_1'])
             multiplicity_fraction_dict = {1: 1-binary_fraction, 2: binary_fraction, 3: 0, 4: 0}
 
         elif self.grid_options['multiplicity_fraction_function'] == 2:
             # Raghavan 2010 will be used
-            # TODO: add verbose logging
-            print("Using Rhagavan 2010 to calculate multiplicity fractions")
+            verbose_print(
+                "_calculate_multiplicity_fraction: Using Rhagavan 2010 to calculate multiplicity fractions",
+                self.grid_options["verbosity"],
+                1,
+            )
+
             binary_fraction = raghavan2010_binary_fraction(system_dict['M_1'])
             multiplicity_fraction_dict = {1: 1-binary_fraction, 2: binary_fraction, 3: 0, 4: 0}
 
@@ -3964,7 +3969,11 @@ eccentricity3=0
 
             # Second: is the moecache filled.
             if not Moecache:
-                print("Moecache is empty. It needs to be filled with the data for the interpolators. Loading the data now")
+                verbose_print(
+                    "_calculate_multiplicity_fraction: Moecache is empty. It needs to be filled with the data for the interpolators. Loading the data now",
+                    self.grid_options["verbosity"],
+                    1,
+                )
 
                 # Load the data
                 self._load_moe_distefano_data()
@@ -3976,18 +3985,23 @@ eccentricity3=0
             self.grid_options['m&s_options']['M_1'] = system_dict['M_1']
 
             # Calculate the multiplicity fraction
-            ding = Moe_de_Stefano_2017_multiplicity_fractions(self.grid_options['m&s_options'], self.grid_options["verbosity"])
-            print(ding)
+            multiplicity_fraction_list = Moe_de_Stefano_2017_multiplicity_fractions(self.grid_options['m&s_options'], self.grid_options["verbosity"])
+
+            # Turn into dict
+            multiplicity_fraction_dict = {el+1:multiplicity_fraction_list[el] for el in range(len(multiplicity_fraction_list))}
 
             # Set the prev value back
-            self.grid_options['m&s_options'][prev_M1_value_ms]
+            self.grid_options['m&s_options']['M_1'] = prev_M1_value_ms
 
         # we don't know what to do next
         else:
             msg = "Chosen value for the multiplicity fraction function is not known."
             raise ValueError(msg)
 
-        # TODO: add extra verbose logging
-        print("Multiplicity: {} multiplicity_fraction: {}".format(system_dict['multiplicity'], multiplicity_fraction_dict[system_dict['multiplicity']]))
+        verbose_print(
+            "Multiplicity: {} multiplicity_fraction: {}".format(system_dict['multiplicity'], multiplicity_fraction_dict[system_dict['multiplicity']]),
+            self.grid_options["verbosity"],
+            1,
+        )
 
-        return multiplicity_fraction_dict[system_dict['multiplicity']]
\ No newline at end of file
+        return multiplicity_fraction_dict[system_dict['multiplicity']]