Merge pull request #10 from mdcourse/some-improvements

Some improvements

Author: simongravelle

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()

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

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):

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

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."""

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