diff --git a/hoomd-blue/research/diblock-domain-fission/equi.py b/hoomd-blue/research/diblock-domain-fission/equi.py
new file mode 100644
index 0000000..35fc701
--- /dev/null
+++ b/hoomd-blue/research/diblock-domain-fission/equi.py
@@ -0,0 +1,55 @@
+
+import sys
+import time
+import h5py
+import shutil
+import numpy as np
+
+import hoomd
+hoomd.context.initialize()
+import hoomd.hdf5
+import hoomd.md
+
+
+def main(argv):
+    with hoomd.context.SimulationContext():
+        system = hoomd.init.read_gsd("start.gsd", time_step=0)
+
+        nl = hoomd.md.nlist.tree()
+        nl.reset_exclusions([])
+
+        qr = []
+        qr += ['temperature','potential_energy','kinetic_energy']
+        qr += ["pressure_xx","pressure_yy","pressure_zz","pressure_xy","pressure_xz","pressure_yz"]
+        qr += ["bond_harmonic_energy"]
+        qr += ["pair_dpd_energy",]
+
+        harmonic = hoomd.md.bond.harmonic()
+        harmonic.bond_coeff.set("backbone", k=16/3, r0=0)
+
+        seed = int(time.time())
+        dpd = hoomd.md.pair.dpd(r_cut=1.0, nlist=nl, seed=seed, kT=1.)
+        dpd.pair_coeff.set("A", "A", A=5., gamma=1.0)
+        dpd.pair_coeff.set("A", "B", A=7., gamma=1.0)
+        dpd.pair_coeff.set("B", "B", A=5., gamma=1.0)
+
+        hoomd.md.integrate.nve(group=hoomd.group.all())
+        hoomd.md.integrate.mode_standard(dt=1e-4)
+
+        # writes out trajectory every 1000 MD steps
+        gsd = hoomd.dump.gsd("equi_trajectory.gsd", group=hoomd.group.all(), period=1000, overwrite=True)
+
+        with hoomd.hdf5.File("equi_data.h5", "w") as h5file:
+            hoomd.hdf5.log(h5file, quantities=qr, period=1000)
+            hoomd.run(1e6, limit_multiple=10000)
+        final = hoomd.dump.gsd("equi_final.gsd", group=hoomd.group.all(), period=None, truncate=True, overwrite=True)
+
+        # Copy result of equi to preprocess for run.py
+        Nreplica = 14
+        for replica_num in range(Nreplica):
+            shutil.copy("equi_final.gsd",
+                        f"final_{replica_num}.gsd")
+
+
+if __name__ == "__main__":
+    main(sys.argv)
diff --git a/hoomd-blue/research/diblock-domain-fission/execute.sh b/hoomd-blue/research/diblock-domain-fission/execute.sh
new file mode 100644
index 0000000..45bda79
--- /dev/null
+++ b/hoomd-blue/research/diblock-domain-fission/execute.sh
@@ -0,0 +1,6 @@
+#!/usr/bin/bash
+
+python gen_gsd.py
+python equi.py
+python run.py
+python plot_result.py
diff --git a/hoomd-blue/research/diblock-domain-fission/gen_gsd.py b/hoomd-blue/research/diblock-domain-fission/gen_gsd.py
new file mode 100644
index 0000000..616d3dd
--- /dev/null
+++ b/hoomd-blue/research/diblock-domain-fission/gen_gsd.py
@@ -0,0 +1,59 @@
+import sys
+import numpy as np
+import gsd
+import gsd.hoomd
+
+class System:
+
+    def __init__(self,n,N=256,b0=0.75):
+        self.L = np.asarray([10, 10, 10])
+        self.N = N
+        self.n = n
+        self.b0 = b0
+
+def main(argv):
+    if len(argv) != 1:
+        print("Usage: ./gen_gsd.py")
+        return
+
+    n = 200
+    system = System(n)
+
+    snapshot = gsd.hoomd.Snapshot()
+
+    snapshot.configuration.box = np.concatenate((system.L,[0,0,0]))
+
+    snapshot.particles.N = system.n * system.N
+
+    snapshot.particles.types = ["A", "B"]
+    snapshot.particles.position = np.zeros( (snapshot.particles.N,3) )
+    snapshot.particles.velocity = np.random.standard_normal( (snapshot.particles.N,3) )
+    snapshot.particles.image = np.zeros( (snapshot.particles.N,3),dtype=np.int )
+    snapshot.particles.typeid = np.zeros( snapshot.particles.N, dtype = np.int )
+    snapshot.bonds.N = system.n * (system.N-1)
+    snapshot.bonds.types = ["backbone"]
+    snapshot.bonds.typeid = np.zeros( snapshot.bonds.N )
+    snapshot.bonds.group = []
+
+    for poly in range(system.n):
+        mono = 0
+        snapshot.particles.position[ poly * system.N + mono ] = [0,-system.L[1]/2,0] + np.random.standard_normal(3)*1e-3
+
+        for mono in range(1, system.N):
+            snapshot.particles.position[ poly * system.N +mono] = snapshot.particles.position[poly*system.N + mono-1] + np.random.standard_normal( 3 ) * system.b0/3
+
+            snapshot.bonds.group.append( [ poly*system.N+mono-1 , poly*system.N+mono])
+            if mono > 16:
+                snapshot.particles.typeid[poly * system.N +mono] = 1
+
+    snapshot.particles.image += np.rint(snapshot.particles.position/system.L).astype(np.int64)
+    snapshot.particles.position -= np.rint(snapshot.particles.position/system.L) * system.L
+
+    snapshot.particles.validate()
+
+    with gsd.hoomd.open( "start.gsd","wb") as f:
+        f.append(snapshot)
+
+
+if __name__ == "__main__":
+    main(sys.argv)
diff --git a/hoomd-blue/research/diblock-domain-fission/plot_result.py b/hoomd-blue/research/diblock-domain-fission/plot_result.py
new file mode 100644
index 0000000..396195a
--- /dev/null
+++ b/hoomd-blue/research/diblock-domain-fission/plot_result.py
@@ -0,0 +1,91 @@
+#!/usr/bin/env python
+import sys
+sys.path.append("..")
+import numpy as np
+import importlib
+
+from droplet import *
+import hoomd
+import hoomd.md
+
+import pysages
+import pickle
+import matplotlib.pyplot as plt
+
+from gyana.plot_util import savefig, plot_data_to_str
+
+def plot_energy(energies, centers, epsilons):
+    fig, ax = plt.subplots()
+
+    ax.set_xlabel(r"distance R $[\sigma]$")
+    ax.set_ylabel(r"free energy $[\epsilon]$")
+
+    # ax.set_ylim((0,850))
+    ax.set_xlim((0,6))
+
+    attachement_string = ""
+    for i in range(len(energies)):
+        center = centers[i]
+        free_energy = energies[i]
+        free_energy = free_energy
+        free_energy -= free_energy[0]
+        epsilon = epsilons[i]
+        ax.plot(center, free_energy, '.-', label=rf"$\Delta A = {epsilon}$")
+        attachement_string += plot_data_to_str((center, free_energy), r"$\Delta\epsilon = {epsilon-1}k_BT$")
+    ax.legend(loc="best")
+    savefig(fig, attachement_string, "energy.pdf")
+    plt.close(fig)
+
+# def plot_hist(result, bins=50):
+#     fig, ax = plt.subplots(2, 2)
+
+#     # ax.set_xlabel("CV")
+#     # ax.set_ylabel("p(CV)")
+
+#     counter = 0
+#     hist_per = len(result["centers"]) // 4 + 1
+#     for x in range(2):
+#         for y in range(2):
+#             for i in range(hist_per):
+#                 if counter + i < len(result["centers"]):
+#                     center = np.asarray(result["centers"][counter + i])
+#                     histo, edges = result["histograms"][counter + i].get_histograms(bins=bins)
+#                     edges = np.asarray(edges)[0]
+#                     edges = (edges[1:] + edges[:-1]) / 2
+#                     ax[x, y].plot(edges, histo, label=f"center {center}")
+#                     ax[x, y].legend(loc="best", fontsize="xx-small")
+#                     ax[x, y].set_yscale("log")
+#             counter += hist_per
+#     while counter < len(result["centers"]):
+#         center = np.asarray(result["centers"][counter])
+#         histo, edges = result["histograms"][counter].get_histograms(bins=bins)
+#         edges = np.asarray(edges)[0]
+#         edges = (edges[1:] + edges[:-1]) / 2
+#         ax[1, 1].plot(edges, histo, label=f"center {center}")
+#         counter += 1
+
+#     attachement_string = str(result)
+#     savefig(fig, attachement_string, "hist.pdf")
+#     plt.close(fig)
+
+
+def main(argv):
+    if len(argv) != 0:
+        print("Usage: ./plot_result.py")
+        return
+
+
+    deltaAs = [0.1, 0.2, 0.3, 0.4]
+    centers = []
+    energies = []
+    for delta in deltaAs:
+        with open(f"result_{delta}.pkl", "rb") as pkl_file:
+            result = pickle.load(pkl_file)
+        centers.append(np.asarray(result["centers"]))
+        energies.append(np.asarray(result["free_energy"]))
+
+    plot_energy(energies, centers, deltaAs)
+
+
+if __name__ == "__main__":
+    main(sys.argv[1:])
diff --git a/hoomd-blue/research/diblock-domain-fission/run.py b/hoomd-blue/research/diblock-domain-fission/run.py
new file mode 100644
index 0000000..0ed980e
--- /dev/null
+++ b/hoomd-blue/research/diblock-domain-fission/run.py
@@ -0,0 +1,116 @@
+import sys
+import time
+import datetime
+import h5py
+import numpy as np
+import importlib
+
+from droplet import *
+
+import hoomd
+import hoomd.md
+
+import pysages
+from pysages.colvars import Distance
+from pysages.methods import UmbrellaIntegration, SerialExecutor
+
+import gsd
+import gsd.hoomd
+
+import pickle
+
+def generate_context(**kwargs):
+    # if kwargs.get("mpi", False):
+    #     MPI = importlib.import_module("mpi4py.MPI")
+    #     init_kwargs = {"mpi_comm": MPI.COMM_SELF}
+    # else:
+    #     init_kwargs = {}
+    # hoomd.context.initialize("--single-mpi", **init_kwargs)
+    hoomd.context.initialize("--single-mpi")
+
+    context = hoomd.context.SimulationContext()
+    with context:
+        deltaA = kwargs.get("deltaA")
+        print(f"Operating replica {kwargs.get('replica_num')}: {deltaA}")
+        system = hoomd.init.read_gsd(f"final_{kwargs.get('replica_num')}.gsd", time_step=0)
+        box = system.box
+        nl = hoomd.md.nlist.tree()
+        nl.reset_exclusions([])
+
+        # Note that pysages energies from the bias potential can not be logged from hoomd
+        qr = []
+        qr += ['temperature','potential_energy','kinetic_energy']
+        qr += ["bond_harmonic_energy"]
+        qr += ["pair_dpd_energy",]
+
+        harmonic = hoomd.md.bond.harmonic()
+        harmonic.bond_coeff.set("backbone", k=16/3, r0=0)
+
+        seed = int(time.time()) ^ kwargs.get("replica_num")
+        dpd = hoomd.md.pair.dpd(r_cut=1.0, nlist=nl, seed=seed, kT=1.)
+        dpd.pair_coeff.set("A", "A", A=5., gamma=1.0)
+        dpd.pair_coeff.set("A", "B", A=5.+deltaA, gamma=1.0)
+        dpd.pair_coeff.set("B", "B", A=5., gamma=1.0)
+
+        hoomd.md.integrate.nve(group=hoomd.group.all())
+        hoomd.md.integrate.mode_standard(dt=1e-3)
+
+        hoomd.analyze.log(f"data_{kwargs.get('replica_num')}_{deltaA}.dat", qr, period=10000, overwrite=True)
+    return context
+
+def post_run_action(**kwargs):
+    hoomd.util.quiet_status()
+    deltaA = kwargs.get("deltaA")
+    hoomd.dump.gsd(
+        filename=f"final_{kwargs.get('replica_num'):03d}_{deltaA}.gsd",
+        overwrite=True,
+        period=None,
+        group=hoomd.group.all(),
+    )
+    hoomd.util.unquiet_status()
+
+def get_executor(mpi):
+    if mpi:
+        futures = importlib.import_module("mpi4py.futures")
+        return futures.MPIPoolExecutor()
+    return SerialExecutor()
+
+
+def main(argv):
+
+    # Find indeces of the droplet forming polymer part
+    with gsd.hoomd.open("../start.gsd", "rb") as start_file:
+        frame = start_file[0]
+        idx = np.where(frame.particles.typeid == frame.particles.types.index("A"))[0].astype("int")
+        box = frame.configuration.box[:3]
+
+    Nidx = len(idx)
+    idxA = list(idx)[:Nidx//2]
+    idxB = list(idx)[Nidx//2:]
+
+    start = datetime.datetime.now()
+
+    for deltaA in [0.0, 0.1, 0.2, 0.3, 0.4]:
+        cvs = [Distance((idxA, idxB))]
+        centers = list(np.linspace(0, 6, 14))
+        method = UmbrellaIntegration(cvs, 1000., centers, 100, int(1e6))
+
+        context_args = {"mpi": True}
+        context_args["deltaA"] = deltaA
+
+        executor = get_executor(context_args["mpi"])
+
+        raw_result = pysages.run(method, generate_context, 2e6,
+                                 context_args=context_args,
+                                 post_run_action=post_run_action,
+                                 executor=executor)
+        result = pysages.analyze(raw_result)
+        print("end", time.asctime())
+        end = datetime.datetime.now()
+        diff = end - start
+        with open(f"result_{deltaA}.pkl", "wb") as file_handle:
+            pickle.dump(result, file_handle)
+
+
+if __name__ == "__main__":
+    main(sys.argv)