diff --git a/integration_tests/test_monte_carlo_move.py b/integration_tests/test_monte_carlo_move.py index 0548b37..0dd00f4 100644 --- a/integration_tests/test_monte_carlo_move.py +++ b/integration_tests/test_monte_carlo_move.py @@ -1,5 +1,8 @@ import sys import os +import time +import pstats +import cProfile import unittest import numpy as np from pint import UnitRegistry @@ -30,7 +33,7 @@ def setUp(self): # Initialize the MonteCarlo object self.mc = MonteCarlo( ureg = ureg, - maximum_steps = 10000, + maximum_steps = 100000, thermo_period = 1000, dumping_period = 1000, number_atoms = [nmb_1], @@ -41,20 +44,58 @@ def setUp(self): cut_off = rc, thermo_outputs = "Epot-press", desired_temperature = T, - neighbor = 20, + neighbor = 50, displace_mc = displace_mc, ) + """ def test_monte_carlo_run(self): - """Test if the Monte Carlo simulation runs without errors.""" + # Test if the Monte Carlo simulation runs without errors. try: # Run the Monte Carlo simulation (this should not raise an exception) + ti = time.time() + + # Profile the run() method self.mc.run() + # self.mc.run() # If it runs successfully, assert True + tf = time.time() + + print("Duration:", np.round(tf-ti, 2), "s") self.assertTrue(True) except Exception as e: # If any exception occurs, fail the test and print the error self.fail(f"Monte Carlo simulation failed with error: {e}") + """ + + def test_monte_carlo_run(self): + """Test if the Monte Carlo simulation runs without errors.""" + profiler = cProfile.Profile() + try: + # Start profiling before running the Monte Carlo simulation + profiler.enable() + + # Run the Monte Carlo simulation + ti = time.time() + self.mc.run() # Assuming self.mc is your Monte Carlo simulation object + + # Stop profiling after the run + profiler.disable() + + tf = time.time() + print("Duration:", np.round(tf - ti, 2), "s") + + # Convert the profiler stats into a readable format + stats = pstats.Stats(profiler) + stats.strip_dirs() # Remove extraneous directory information + stats.sort_stats('time') # Sort by time spent in function + stats.print_stats(10) # Print top 10 slowest functions + + self.assertTrue(True) + except Exception as e: + # If any exception occurs, fail the test and print the error + self.fail(f"Monte Carlo simulation failed with error: {e}") + if __name__ == "__main__": unittest.main() diff --git a/molecular_simulation_code/measurements_utilities.py b/molecular_simulation_code/measurements_utilities.py index 56f30de..7b730b2 100644 --- a/molecular_simulation_code/measurements_utilities.py +++ b/molecular_simulation_code/measurements_utilities.py @@ -1,6 +1,6 @@ from numba import njit import numpy as np -from typing import List +from numba.typed import List from distances_utilities import compute_vector_matrix from forces_utilities import compute_force_matrix diff --git a/molecular_simulation_code/monte_carlo.py b/molecular_simulation_code/monte_carlo.py index 31f9fb5..d6caa55 100644 --- a/molecular_simulation_code/monte_carlo.py +++ b/molecular_simulation_code/monte_carlo.py @@ -6,6 +6,7 @@ from initialize_simulation import InitializeSimulation from measurements_utilities import compute_epot from monte_carlo_utilities import calculate_Lambda +# from numba_test import numba_copy import warnings warnings.filterwarnings('ignore') @@ -52,7 +53,9 @@ def monte_carlo_move(self): self.Epot = compute_epot(self.neighbor_lists, self.atoms_positions, self.box_mda, self.cross_coefficients) # Make a copy of the initial atom positions and initial energy initial_Epot = self.Epot - initial_positions = copy.deepcopy(self.atoms_positions) + # initial_positions = copy.deepcopy(self.atoms_positions) + # initial_positions = np.copy(self.atoms_positions) + # initial_positions = numba_copy(self.atoms_positions) # Pick an atom id randomly atom_id = np.random.randint(np.sum(self.number_atoms)) # Move the chosen atom in a random direction @@ -62,9 +65,12 @@ def monte_carlo_move(self): move = np.append(move, 0.0) # Pad with zero for the z-component else: # 3D case move = (np.random.random(3) - 0.5) * self.displace_mc + initial_position_atom = np.copy(self.atoms_positions[atom_id]) self.atoms_positions[atom_id] += move + # Measure the potential energy of the new configuration trial_Epot = compute_epot(self.neighbor_lists, self.atoms_positions, self.box_mda, self.cross_coefficients) + # Evaluate whether the new configuration should be kept or not beta = 1/self.desired_temperature delta_E = trial_Epot-initial_Epot @@ -81,7 +87,8 @@ def monte_carlo_move(self): self.Epot = trial_Epot self.successful_move += 1 else: # Reject new position - self.atoms_positions = initial_positions # Revert to initial positions + self.atoms_positions[atom_id] = initial_position_atom + # self.atoms_positions = initial_positions # Revert to initial positions self.failed_move += 1 def monte_carlo_swap(self): diff --git a/molecular_simulation_code/numba_test.py b/molecular_simulation_code/numba_test.py new file mode 100644 index 0000000..236c9be --- /dev/null +++ b/molecular_simulation_code/numba_test.py @@ -0,0 +1,12 @@ +import numpy as np +from numba import njit + +@njit +def numba_copy(array): + # Create an empty array of the same shape and type + new_array = np.empty_like(array) + # Copy elements manually + for i in range(array.shape[0]): + for j in range(array.shape[1]): + new_array[i, j] = array[i, j] + return new_array \ No newline at end of file diff --git a/molecular_simulation_code/simulation_handler.py b/molecular_simulation_code/simulation_handler.py index ccd30c0..b4752a4 100644 --- a/molecular_simulation_code/simulation_handler.py +++ b/molecular_simulation_code/simulation_handler.py @@ -1,4 +1,5 @@ import numpy as np +from numba.typed import List from contacts_utilities import contact_matrix, compute_neighbor_lists @@ -18,7 +19,12 @@ def update_neighbor_lists(self, force_update=False): box=self.box_mda) # Compute the neighbor lists from the contact matrix - self.neighbor_lists = compute_neighbor_lists(matrix) + python_neighbor_lists = compute_neighbor_lists(matrix) + + # Convert Python list to numba.typed.List + self.neighbor_lists = List() + for neighbors in python_neighbor_lists: + self.neighbor_lists.append(neighbors) def update_cross_coefficients(self, force_update=False): """Update the Lennard-Jones cross-coefficients for all atom pairs.""" diff --git a/molecular_simulation_code/simulation_logger.py b/molecular_simulation_code/simulation_logger.py index 1f2e838..f12f4c6 100644 --- a/molecular_simulation_code/simulation_logger.py +++ b/molecular_simulation_code/simulation_logger.py @@ -9,8 +9,10 @@ def setup_logger(folder_name, overwrite=False): logger.setLevel(logging.INFO) logger.propagate = False # Disable propagation to prevent double logging - # Clear any existing handlers if this function is called again + # Close and clear existing handlers, if any if logger.hasHandlers(): + for handler in logger.handlers: + handler.close() # Ensure handlers are properly closed logger.handlers.clear() # Create handlers for console and file