Refactoring

MCMC/MCMC-flow/init.py

@@ -38,9 +38,13 @@ def get_parameters():
     parameters["trajectory_write_freq"] = [10000]
 
     # run parameters
-    parameters["n_steps"] = [[1e7, 1e8]]
-    parameters["kT"] = [[10, 1.5]]
-    parameters["max_trans"] = [[3.0, 0.5]]
+    parameters["mixing_steps"] = [2000]
+    parameters["mixing_kT"] = [10]
+    parameters["mixing_max_trans"] = [0.5]
+
+    parameters["n_steps"] = [[2000]]
+    parameters["kT"] = [[1.5]]
+    parameters["max_trans"] = [[0.5]]
     parameters["seed"] = [20]
 
     return list(parameters.keys()), list(product(*parameters.values()))
@@ -58,24 +62,14 @@ def main():
         parent_job = project.open_job(parent_statepoint)
         parent_job.init()
         parent_job.doc.setdefault("done", False)
+        parent_job.doc.setdefault("current_run", 0)
+        parent_job.doc.setdefault("mixed", False)
         parent_job.doc.setdefault("timestep", [])
         parent_job.doc.setdefault("accepted_moves", [])
         parent_job.doc.setdefault("rejected_moves", [])
         parent_job.doc.setdefault("acceptance_ratio", [])
         parent_job.doc.setdefault("tps", [])
         parent_job.doc.setdefault("energy", [])
-        # each pair of (`n_steps`, `kT`) defines a phase of simulation run. The `phase_{i}` key in job doc determines
-        # whether that phase is already done or not. False means phase is not done.
-        n_steps_size = len(parent_statepoint['n_steps'])
-        kT_size = len(parent_statepoint['kT'])
-        max_trans_size = len(parent_statepoint['max_trans'])
-        if n_steps_size == kT_size == max_trans_size:
-            for i in range(1, n_steps_size + 1):
-                parent_job.doc.setdefault("phase_{}".format(i), False)
-        else:
-            raise ValueError(
-                "These lists must have the same length: `n_steps`, `kT` and `max_trans` for each job! \n "
-                "`n_step` size: {}, `kT` size: {}, `max_trans` size: {}".format(n_steps_size, kT_size, max_trans_size))
     if custom_job_doc:
         for key in custom_job_doc:
             parent_job.doc.setdefault(key, custom_job_doc[key])

MCMC/MCMC-flow/project.py

@@ -96,7 +96,7 @@ def sample(job):
         print("----------------------")
 
         # Setting up the system
-        restart = job.isfile("system.txt")
+        restart = job.isfile("restart.gsd")
         sim = Simulation(n_particles=job.sp.n_particles, n_density=job.sp.n_density, r=job.sp.r, r_cut=job.sp.r_cut,
                          energy_write_freq=job.sp.energy_write_freq, trajectory_write_freq=job.sp.trajectory_write_freq,
                          energy_func=ENERGY_FUNCS[job.sp.energy_func], hard_sphere=job.sp.hard_sphere, restart=restart,
@@ -106,27 +106,44 @@ def sample(job):
         print("----------------------")
         print("System generated...")
         print("----------------------")
+
+        # running mixing step if not done before
+        if not job.doc["mixed"]:
+            print("----------------------")
+            print("Starting mixing simulation...")
+            print("----------------------")
+            sim.run(n_steps=job.sp.mixing_steps, kT=job.sp.mixing_kT, max_trans=job.sp.mixing_max_trans)
+            sim.save_trajectory(fname="trajectory_mixing.gsd")
+            job.doc["timestep"].append(sim.timestep)
+            job.doc["accepted_moves"].append(sim.accepted_moves)
+            job.doc["rejected_moves"].append(sim.rejected_moves)
+            job.doc["acceptance_ratio"].append(sim.acceptance_ratio)
+            job.doc["tps"].append(sim.tps)
+            job.doc["energy"].append(sim.energy)
+            job.doc["mixed"] = True
+            sim.reset_system()
+
         print("----------------------")
         print("Starting simulation...")
         print("----------------------")
         # Running the simulation
-        for i, (n_steps, kT, max_trans) in enumerate(zip(job.sp.n_steps, job.sp.kT, job.sp.max_trans)):
-            current_run = i+1
-            if not job.doc["phase_{}".format(current_run)]:
-                sim.run(n_steps=n_steps, kT=kT, max_trans=max_trans)
-                sim.save_trajectory(fname="trajectory_{}.gsd".format(i))
-                sim.save_system()
-                job.doc["timestep"].append(sim.timestep)
-                job.doc["accepted_moves"].append(sim.accepted_moves)
-                job.doc["rejected_moves"].append(sim.rejected_moves)
-                job.doc["acceptance_ratio"].append(sim.acceptance_ratio)
-                job.doc["tps"].append(sim.tps)
-                job.doc["energy"].append(sim.energy)
-                job.doc["phase_{}".format(current_run)] = True
-                sim.reset_system()
-
-        if sum(job.doc["timestep"]) == sum(job.sp.n_steps):
-            job.doc["done"] = True
+        for (n_steps, kT, max_trans) in zip(job.sp.n_steps, job.sp.kT, job.sp.max_trans):
+            sim.run(n_steps=n_steps, kT=kT, max_trans=max_trans)
+            sim.save_trajectory(fname="trajectory_{}.gsd".format(job.doc["current_run"]))
+            sim.save_snapshot('restart.gsd')
+            job.doc["timestep"].append(sim.timestep)
+            job.doc["accepted_moves"].append(sim.accepted_moves)
+            job.doc["rejected_moves"].append(sim.rejected_moves)
+            job.doc["acceptance_ratio"].append(sim.acceptance_ratio)
+            job.doc["tps"].append(sim.tps)
+            job.doc["energy"].append(sim.energy)
+            job.doc["phase_{}".format(job.doc["current_run"])] = True
+            sim.reset_system()
+            job.doc["current_run"] += 1
+
+
+        # TODO: Find a way to check the timestep based on all the run attempts from PT.
+        job.doc["done"] = True
         print("-----------------------------")
         print("Simulation finished completed")
         print("-----------------------------")

MCMC/simulation.py

@@ -76,7 +76,7 @@ def _init_system(self, restart):
         x and y coordinates are between -L/2 and L/2.
         :return: A 2D numpy array of shape (self.n_particles, 2).
         """
-        if restart and os.path.isfile("system.txt"):
+        if restart and os.path.isfile("restart.gsd"):
             system = self._load_system()
         else:
             disks_per_row = math.floor((self.L - self.r) / (2 * self.r))
@@ -211,13 +211,18 @@ def visualize(self, frame_number=-1, save_path=None):
             plt.savefig(os.path.join(save_path, fig_name))
         return plt
 
-    def save_system(self, frame=-1):
+    def save_snapshot(self,fname="restart.gsd"):
         """
-        Save a single snapshot of system_history to a .txt file.
-        :param frame: Index number of system_history
+        Save a snapshot of system to a .gsd file.
+        :param fname: name of the file.
         """
-        system_array = self.system_history[frame]
-        np.savetxt(fname="system.txt", X=system_array)
+        with gsd.hoomd.open(fname, 'wb') as traj:
+            snap = gsd.hoomd.Snapshot()
+            snap.particles.N = self.n_particles
+            snap.configuration.box = [self.L, self.L, self.L, 0, 0, 0]
+            snap.particles.type = ['A']
+            snap.particles.position = np.append(self.system, np.zeros((self.system.shape[0], 1)), axis=1)
+            traj.append(snap)
 
     def save_trajectory(self, fname="traj.gsd"):
         with gsd.hoomd.open(fname, 'wb') as traj:
@@ -248,10 +253,18 @@ def _update_log_file(self):
                         file.write(f"{e},{t}" + "\n")
 
     def _load_system(self):
-        system = np.loadtxt(fname="system.txt")
+        snapshot = gsd.hoomd.open("restart.gsd")[0]
         try:
-            assert system.shape[0] == self.n_particles
+            assert snapshot.particles.N == self.n_particles
         except AssertionError:
             raise AssertionError(
                 "Number of particles in the saved system is not equal to the specified number of particles!")
-        return system
+
+        try:
+            assert snapshot.configuration.box[0] == self.L
+        except AssertionError:
+            raise AssertionError(
+                "Box size in the saved system is not equal to the box size!")
+
+        sys = snapshot.particles.position[:, :2]
+        return sys