diff --git a/MCMC/MCMC-flow/init.py b/MCMC/MCMC-flow/init.py index 15a4142..d99869a 100644 --- a/MCMC/MCMC-flow/init.py +++ b/MCMC/MCMC-flow/init.py @@ -25,7 +25,7 @@ def get_parameters(): parameters["r"] = [0.5] parameters["r_cut"] = [2.5] parameters["energy_func"] = ["lj"] - parameters["hard_sphere"] = [True] + parameters["hard_sphere"] = [False] # LJ energy parameters parameters["epsilon"] = [1.0] @@ -35,12 +35,16 @@ def get_parameters(): # logging parameters parameters["energy_write_freq"] = [1000] - parameters["trajectory_write_freq"] = [10000] + parameters["trajectory_write_freq"] = [1000] # run parameters - parameters["n_steps"] = [[1e7, 1e8]] - parameters["kT"] = [[10, 1.5]] - parameters["max_trans"] = [[3.0, 0.5]] + parameters["mixing_steps"] = [10e5] + parameters["mixing_kT"] = [10] + parameters["mixing_max_trans"] = [0.5] + + parameters["n_steps"] = [[10e5]] + 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]) diff --git a/MCMC/MCMC-flow/project.py b/MCMC/MCMC-flow/project.py index 22b0018..92c1322 100644 --- a/MCMC/MCMC-flow/project.py +++ b/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("-----------------------------") diff --git a/MCMC/simulation.py b/MCMC/simulation.py index 24c7e22..3221ab3 100644 --- a/MCMC/simulation.py +++ b/MCMC/simulation.py @@ -22,7 +22,7 @@ def __init__( trajectory_write_freq=10000, seed=20, energy_func=None, - hard_sphere=True, + hard_sphere=False, restart=False, **kwargs): """ @@ -53,7 +53,7 @@ def __init__( self.energies = [copy.deepcopy(self.energy)] self.temperatures = [] self._tps = [] - + random.seed(seed) @property @@ -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: @@ -232,6 +237,8 @@ def save_trajectory(self, fname="traj.gsd"): def reset_system(self): """Clear the system history and reset the system.""" self.system_history.clear() + self.energies.clear() + self.energies.clear() self.timestep = 0 self.accepted_moves = 0 self.rejected_moves = 0 @@ -248,10 +255,17 @@ 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 np.isclose(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 diff --git a/MCMC/utils.py b/MCMC/utils.py index 6da56ce..8f71fe4 100644 --- a/MCMC/utils.py +++ b/MCMC/utils.py @@ -74,7 +74,7 @@ def inverse_distance_repulsive(distances): return (1. / np.power(distances, 2)).sum() -def lj_energy(distances, epsilon=1.0, sigma=1.0, n=12, m=6): +def lj_energy(distances, epsilon=1.0, sigma=1.0, n=12, m=6, e_factor=1): """ Calculates Lennard-Jones pair potential. :param distances: 1D array of distances. @@ -82,9 +82,10 @@ def lj_energy(distances, epsilon=1.0, sigma=1.0, n=12, m=6): :param sigma: Particle size sigma. :param n: Repulsive power factor. :param m: Attractive power factor. + :param e_factor: epsilon coefficient. :return: Total energy. """ - return 4 * epsilon * (np.power(sigma / distances, n) - np.power(sigma / distances, m)).sum() + return 4 * (e_factor* epsilon) * (np.power(sigma / distances, n) - np.power(sigma / distances, m)).sum() @jit(nopython=True)