diff --git a/.circleci/config.yml b/.circleci/config.yml
index 249067ac..53621ef0 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -56,9 +56,10 @@ jobs:
             pip3 install pytest
             pip3 install pytest-mpi
             pip3 install pytest-cov
-            COVERAGE_FILE=.coverage_1 pytest -vv --disable-pytest-warnings --cov=ensemble_md --cov-report=xml --color=yes ensemble_md/tests/
-            COVERAGE_FILE=.coverage_2 mpirun -np 4 pytest -vv --disable-pytest-warnings --cov=ensemble_md --cov-report=xml --color=yes ensemble_md/tests/test_mpi_func.py --with-mpi
-            coverage combine .coverage_*
+            pytest -vv --disable-pytest-warnings --cov=ensemble_md --cov-report=xml --color=yes ensemble_md/tests/
+            # COVERAGE_FILE=.coverage_1 pytest -vv --disable-pytest-warnings --cov=ensemble_md --cov-report=xml --color=yes ensemble_md/tests/
+            # COVERAGE_FILE=.coverage_2 mpirun -np 4 pytest -vv --disable-pytest-warnings --cov=ensemble_md --cov-report=xml --color=yes ensemble_md/tests/test_mpi_func.py --with-mpi
+            # coverage combine .coverage_*
 
       - run:
           name: CodeCov
diff --git a/ensemble_md/tests/data/mdp/compare_1.mdp b/ensemble_md/tests/data/mdp/compare_1.mdp
new file mode 100644
index 00000000..bddd8923
--- /dev/null
+++ b/ensemble_md/tests/data/mdp/compare_1.mdp
@@ -0,0 +1,122 @@
+; Run control
+integrator = md-vv
+tinit = 0
+dt = 0.002
+nsteps = 100000000
+nstcomm = 10
+
+; Output control
+nstlog = 1000
+nstcalcenergy = 10
+nstenergy = 1000
+nstxout_compressed = 1000
+
+; Neighborsearching and short-range nonbonded interactions
+nstlist = 10
+ns_type = grid
+pbc = xyz
+rlist = 1.3
+
+; Electrostatics
+coulombtype = PME-switch
+rcoulomb_switch = 0.88
+rcoulomb = 0.9
+
+; van der Waals
+vdw_type = switch
+rvdw_switch = 0.85
+rvdw = 0.9
+
+; Apply long range dispersion corrections for Energy and Pressure
+DispCorr = AllEnerPres
+
+; Spacing for the PME/PPPM FFT grid
+fourierspacing = 0.12
+
+; EWALD/PME/PPPM parameters
+pme_order = 4
+ewald_rtol = 1e-05
+ewald_geometry = 3d
+epsilon_surface = 0
+optimize_fft = yes
+
+; Temperature coupling
+tcoupl = nose-hoover
+nsttcouple = 10
+tc_grps = System
+tau_t = 1.0
+ref_t = 300
+
+; Pressure coupling is on for NPT
+pcoupl = no
+
+gen_vel = yes
+gen_temp = 300
+gen_seed = -1
+
+; options for bonds
+constraints = h-bonds
+
+; Type of constraint algorithm
+constraint_algorithm = shake
+shake_tol = 1e-05
+
+; Free energy calculation
+free_energy = expanded
+calc_lambda_neighbors = -1
+sc_alpha = 0.5
+sc_power = 1
+sc_sigma = 0.5
+couple_moltype = MOL
+couple_lambda0 = vdw-q
+couple_lambda1 = none
+couple_intramol = no
+init_lambda_state = 0
+nstdhdl = 100
+dhdl_print_energy = total
+
+; Seed for Monte Carlo in lambda space
+symmetrized_transition_matrix = no
+nst_transition_matrix = 100000
+; wl-scale                 = 0.8
+; wl-ratio                 = 0.7
+; init-wl-delta            = 10
+
+; expanded ensemble variables
+nstexpanded = 100
+lmc_stats = no
+lmc_move = metropolized-gibbs
+; lmc-weights-equil       = wl-delta
+; weight-equil-wl-delta   = 0.001
+; wl-oneovert             = yes
+
+; lambda-states          = 1      2      3      4      5      6      7      8      9      10     11     12     13     14     15     16     17     18     19     20     21     22     23     24     25     26     27     28     29     30     31     32     33     34     35     36     37     38     39     40
+
+coul_lambdas = 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+vdw_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.45 0.5 0.55 0.6 0.63 0.66 0.69 0.72 0.75 0.78 0.81 0.84 0.88 0.92 1.0
+
+; PULL CODE
+pull = yes
+pull_ngroups = 2
+pull_ncoords = 1
+pull_group1_name = HOS
+pull_group2_name = MOL
+pull_pbc_ref_prev_step_com = yes
+
+pull_coord1_groups = 1 2
+pull_coord1_type = umbrella
+pull_coord1_geometry = distance
+pull_coord1_dim = Y Y Y
+pull_coord1_origin = 0.0 0.0 0.0
+pull_coord1_vec = 0.0 0.0 0.0
+pull_coord1_start = yes
+pull_coord1_init = 0
+pull_coord1_rate = 0
+pull_coord1_k = 0
+pull_coord1_kB = 1000
+pull_nstfout = 400000
+pull_nstxout = 1000
+pull-print-ref-value = yes
+
+restraint_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+init_lambda_weights = 0.0 57.88597 112.71883 163.84425 210.48097 253.80261 294.79849 333.90408 370.82669 406.02515 438.53116 468.53751 496.24649 521.58417 544.57404 565.26697 583.7337 599.60651 613.43958 624.70471 633.95947 638.29785 642.44977 646.33551 649.91626 651.54779 652.93359 654.13263 654.94073 655.13086 655.07239 654.66443 653.68683 652.32123 650.72308 649.2381 647.94586 646.599 645.52063 643.99133
diff --git a/ensemble_md/tests/data/mdp/compare_2.mdp b/ensemble_md/tests/data/mdp/compare_2.mdp
new file mode 100644
index 00000000..baa4e137
--- /dev/null
+++ b/ensemble_md/tests/data/mdp/compare_2.mdp
@@ -0,0 +1,122 @@
+; Run control
+integrator = md-vv
+tinit = 0
+dt = 0.002
+nsteps = 100000000
+nstcomm = 10
+
+; Output control
+nstlog = 1000
+nstcalcenergy = 10
+nstenergy = 1000
+nstxout_compressed = 1000
+
+; Neighborsearching and short-range nonbonded interactions
+nstlist = 10
+ns_type = grid
+pbc = xyz
+rlist = 1.3
+
+; Electrostatics
+coulombtype = PME-switch
+rcoulomb_switch = 0.88
+rcoulomb = 0.9
+
+; van der Waals
+vdw_type = switch
+rvdw_switch = 0.85
+rvdw = 0.9
+
+; Apply long range dispersion corrections for Energy and Pressure
+DispCorr = AllEnerPres
+
+; Spacing for the PME/PPPM FFT grid
+fourierspacing = 0.12
+
+; EWALD/PME/PPPM parameters
+pme_order = 4
+ewald_rtol = 1e-05
+ewald_geometry = 3d
+epsilon_surface = 0
+optimize_fft = yes
+
+; Temperature coupling
+tcoupl = nose-hoover
+nsttcouple = 10
+tc_grps = System
+tau_t = 1.0
+ref_t = 300
+
+; Pressure coupling is on for NPT
+pcoupl = no
+
+gen_vel = yes
+gen_temp = 300
+gen_seed = -1
+
+; options for bonds
+constraints = h-bonds
+
+; Type of constraint algorithm
+constraint_algorithm = shake
+shake_tol = 1e-05
+
+; Free energy calculation
+free_energy = expanded
+calc_lambda_neighbors = -1
+sc_alpha = 0.5
+sc_power = 1
+sc_sigma = 0.5
+couple_moltype = MOL
+couple_lambda0 = vdw-q
+couple_lambda1 = none
+couple_intramol = no
+init_lambda_state = 0
+nstdhdl = 100
+dhdl_print_energy = total
+
+; Seed for Monte Carlo in lambda space
+symmetrized_transition_matrix = no
+nst_transition_matrix = 100000
+; wl-scale                 = 0.8
+; wl-ratio                 = 0.7
+; init-wl-delta            = 10
+
+; expanded ensemble variables
+nstexpanded = 100
+lmc_stats = no
+lmc_move = metropolized-gibbs
+; lmc-weights-equil       = wl-delta
+; weight-equil-wl-delta   = 0.001
+; wl-oneovert             = yes
+
+; lambda-states          = 1      2      3      4      5      6      7      8      9      10     11     12     13     14     15     16     17     18     19     20     21     22     23     24     25     26     27     28     29     30     31     32     33     34     35     36     37     38     39     40
+
+coul_lambdas = 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+vdw_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.45 0.5 0.55 0.6 0.63 0.66 0.69 0.72 0.75 0.78 0.81 0.84 0.88 0.92 1.0
+
+; PULL CODE
+pull = yes
+pull_ngroups = 2
+pull_ncoords = 1
+pull_group1_name = HOS
+pull_group2_name = MOL
+pull_pbc_ref_prev_step_com = yes
+
+pull_coord1_groups = 1 2
+pull_coord1_type = umbrella
+pull_coord1_geometry = distance
+pull_coord1_dim = Y Y Y
+pull_coord1_origin = 0.0 0.0 0.0
+pull_coord1_vec = 0.0 0.0 0.0
+pull_coord1_start = yes
+pull_coord1_init = 0
+pull_coord1_rate = 0
+pull_coord1_k = 0
+pull_coord1_kB = 1000
+pull_nstfout = 400000
+pull_nstxout = 1000
+pull_print_ref_value = yes
+
+restraint_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+init_lambda_weights = 0.0 57.88597 112.71883 163.84425 210.48097 253.80261 294.79849 333.90408 370.82669 406.02515 438.53116 468.53751 496.24649 521.58417 544.57404 565.26697 583.7337 599.60651 613.43958 624.70471 633.95947 638.29785 642.44977 646.33551 649.91626 651.54779 652.93359 654.13263 654.94073 655.13086 655.07239 654.66443 653.68683 652.32123 650.72308 649.2381 647.94586 646.599 645.52063 643.99133
diff --git a/ensemble_md/tests/data/mdp/compare_3.mdp b/ensemble_md/tests/data/mdp/compare_3.mdp
new file mode 100644
index 00000000..4b36dabd
--- /dev/null
+++ b/ensemble_md/tests/data/mdp/compare_3.mdp
@@ -0,0 +1,123 @@
+; Run control
+integrator = md-vv
+tinit = 0
+dt = 0.002
+nsteps = 100000000
+nstcomm = 10
+
+; Output control
+nstlog = 1000
+nstcalcenergy = 10
+nstenergy = 1000
+nstxout_compressed = 1000
+
+; Neighborsearching and short-range nonbonded interactions
+nstlist = 10
+ns_type = grid
+pbc = xyz
+rlist = 1.3
+
+; Electrostatics
+coulombtype = PME-switch
+rcoulomb_switch = 0.88
+rcoulomb = 0.9
+
+; van der Waals
+vdw_type = switch
+rvdw_switch = 0.85
+rvdw = 0.9
+
+; Apply long range dispersion corrections for Energy and Pressure
+DispCorr = AllEnerPres
+
+; Spacing for the PME/PPPM FFT grid
+fourierspacing = 0.12
+
+; EWALD/PME/PPPM parameters
+pme_order = 4
+ewald_rtol = 1e-05
+ewald_geometry = 3d
+epsilon_surface = 0
+optimize_fft = yes
+
+; Temperature coupling
+tcoupl = nose-hoover
+nsttcouple = 10
+tc_grps = System
+tau_t = 1.0
+ref_t = 300
+
+; Pressure coupling is on for NPT
+pcoupl = no
+
+gen_vel = yes
+gen_temp = 300
+gen_seed = -1
+
+; options for bonds
+constraints = h-bonds
+
+; Type of constraint algorithm
+constraint_algorithm = shake
+shake_tol = 1e-05
+
+; Free energy calculation
+free_energy = expanded
+calc_lambda_neighbors = -1
+sc_alpha = 0.5
+sc_power = 1
+sc_sigma = 0.5
+couple_moltype = MOL
+couple_lambda0 = vdw-q
+couple_lambda1 = none
+couple_intramol = no
+init_lambda_state = 0
+nstdhdl = 10
+dhdl_print_energy = total
+
+; Seed for Monte Carlo in lambda space
+lmc_seed = -1
+symmetrized_transition_matrix = no
+nst_transition_matrix = 100000
+wl_scale = 0.8
+wl_ratio = 0.7
+init_wl_delta = 10
+
+; expanded ensemble variables
+nstexpanded = 10
+lmc_stats = wang-landau
+lmc_move = metropolized-gibbs
+lmc_weights_equil = wl-delta
+weight_equil_wl_delta = 0.001
+wl_oneovert = yes
+
+; lambda-states          = 1      2      3      4      5      6      7      8      9      10     11     12     13     14     15     16     17     18     19     20     21     22     23     24     25     26     27     28     29     30     31     32     33     34     35     36     37     38     39     40
+
+coul_lambdas = 0.0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+vdw_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.1 0.2 0.3 0.4 0.45 0.5 0.55 0.6 0.63 0.66 0.69 0.72 0.75 0.78 0.81 0.84 0.88 0.92 1.0
+
+; PULL CODE
+pull = yes
+pull_ngroups = 2
+pull_ncoords = 1
+pull_group1_name = HOS
+pull_group2_name = MOL
+pull_pbc_ref_prev_step_com = yes
+
+pull_coord1_groups = 1 2
+pull_coord1_type = umbrella
+pull_coord1_geometry = distance
+pull_coord1_dim = Y Y Y
+pull_coord1_origin = 0.0 0.0 0.0
+pull_coord1_vec = 0.0 0.0 0.0
+pull_coord1_start = yes
+pull_coord1_init = 0
+pull_coord1_rate = 0
+pull_coord1_k = 0
+pull_coord1_kB = 1000
+pull_nstfout = 400000
+pull_nstxout = 1000
+pull_print_ref_value = yes
+
+restraint_lambdas = 0.0 0.0 0.0 0.0 0.0 0.0 0.01 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+
diff --git a/ensemble_md/tests/test_gmx_parser.py b/ensemble_md/tests/test_gmx_parser.py
index dcb5f35c..0caef7ee 100644
--- a/ensemble_md/tests/test_gmx_parser.py
+++ b/ensemble_md/tests/test_gmx_parser.py
@@ -98,3 +98,25 @@ def test_write(self):
         mdp.write('test_2.mdp', skipempty=True)
         os.remove('test_1.mdp')
         os.remove('test_2.mdp')
+
+def test_compare_MDPs():
+    mdp_list = ['ensemble_md/tests/data/mdp/compare_1.mdp', 'ensemble_md/tests/data/mdp/compare_2.mdp', 'ensemble_md/tests/data/mdp/compare_3.mdp']
+    result_1 = gmx_parser.compare_MDPs(mdp_list[:2], print_diff=True)
+    result_2 = gmx_parser.compare_MDPs(mdp_list[1:], print_diff=True)
+    dict_1 = {}  # the first two are the same but just in different formats
+    dict_2 = {
+        'nstdhdl': [100, 10],
+        'wl_oneovert': [None, 'yes'],
+        'weight_equil_wl_delta': [None, 0.001],
+        'init_lambda_weights': [[0.0, 57.88597, 112.71883, 163.84425, 210.48097, 253.80261, 294.79849, 333.90408, 370.82669, 406.02515, 438.53116, 468.53751, 496.24649, 521.58417, 544.57404, 565.26697, 583.7337, 599.60651, 613.43958, 624.70471, 633.95947, 638.29785, 642.44977, 646.33551, 649.91626, 651.54779, 652.93359, 654.13263, 654.94073, 655.13086, 655.07239, 654.66443, 653.68683, 652.32123, 650.72308, 649.2381, 647.94586, 646.599, 645.52063, 643.99133], None],
+        'wl_ratio': [None, 0.7],
+        'lmc_weights_equil': [None, 'wl_delta'],
+        'lmc_stats': ['no', 'wang_landau'],
+        'wl_scale': [None, 0.8],
+        'init_wl_delta': [None, 10],
+        'lmc_seed': [None, -1],
+        'nstexpanded': [100, 10]
+    }
+
+    assert result_1 == dict_1
+    assert result_2 == dict_2
diff --git a/ensemble_md/tests/test_mpi_func.py b/ensemble_md/tests/test_mpi_func.py
index 889dd22e..07557847 100644
--- a/ensemble_md/tests/test_mpi_func.py
+++ b/ensemble_md/tests/test_mpi_func.py
@@ -14,12 +14,8 @@
 import yaml
 import shutil
 import pytest
-from mpi4py import MPI
 from ensemble_md.replica_exchange_EE import ReplicaExchangeEE
 
-comm = MPI.COMM_WORLD
-rank = comm.Get_rank()
-
 
 @pytest.fixture
 def params_dict():
@@ -30,7 +26,7 @@ def params_dict():
         'gmx_executable': 'gmx',
         'gro': 'ensemble_md/tests/data/sys.gro',
         'top': 'ensemble_md/tests/data/sys.top',
-        'mdp': 'ensemble_md/tests/data/expanded.mdp',
+        'mdp': 'ensemble_md/tests/data/mdp/expanded.mdp',
         'n_sim': 4,
         'n_iter': 10,
         's': 1,
@@ -50,6 +46,7 @@ def get_REXEE_instance(input_dict, yml_file='params.yaml'):
     with open(yml_file, 'w') as f:
         yaml.dump(input_dict, f)
     REXEE = ReplicaExchangeEE(yml_file)
+    print(os.path.abspath(yml_file))
     return REXEE
 
 
@@ -97,6 +94,10 @@ def get_gmx_cmd_from_output(output):
 
 @pytest.mark.mpi
 def test_run_grompp(params_dict):
+    from mpi4py import MPI
+    comm = MPI.COMM_WORLD
+    rank = comm.Get_rank()
+
     params_dict['grompp_args'] = {'-maxwarn': '1'}
 
     # Case 1: The first iteration, i.e., n = 0