Fixes and factory

main.py

@@ -1,53 +1,161 @@
 import sys
+import os
+import json
+import numpy as np
+
 from src.ljts.potential import LJTS
 from src.ljts.box import Box
 from src.ljts.orchestrator import Orchestrator, MetropolisMC
+from src.ljts.distortion import compute_distortion
 
 
 def run_with_orchestrator(config_file: str):
     """Run simulation using Orchestrator with JSON configuration."""
     try:
-        # Create orchestrator with configuration file
+        # Load JSON config directly for parameter logging
+        with open(config_file, 'r') as jf:
+            cfg = json.load(jf)
+
+        # Instantiate orchestrator & print its summary
         orchestrator = Orchestrator(config_file)
         orchestrator.print_config_summary()
-
-        # Create the inter-particle potential
+
+        # Extract setup parameters
+        setup = cfg["setup"]
+        Lx, Ly, Lz = setup["Lx"], setup["Ly"], setup["Lz"]
+        T = setup.get("temperature", 0.8)
+
+        # Build potential & box
         potential = LJTS(cutoff=2.5)
-
-        # Get box dimensions from config
-        setup = orchestrator.config["setup"]
         box = Box(
-            len_x=setup["Lx"], 
-            len_y=setup["Ly"], 
-            len_z=setup["Lz"], 
-            den_liq=setup["compartments"][1]["density"],
-            den_vap=setup["compartments"][0]["density"],
+            len_x   = Lx,
+            len_y   = Ly,
+            len_z   = Lz,
+            den_liq = setup["compartments"][1]["density"],
+            den_vap = setup["compartments"][0]["density"],
             potential=potential
         )
-
-        # Set the box in orchestrator
         orchestrator.box = box
-
-        # Print initial stats
+
+        # Extract step parameters
+        sim_cfg      = cfg["steps"]
+        n_pr         = sim_cfg["total"]
+        reset_points = set(sim_cfg.get("reset_sampling_at", []))
+
+        # Extract console + trajectory output parameters
+        console_freq = cfg.get("console_output", {}).get("frequency")
+        traj_cfg     = cfg.get("trajectory_output", {})
+        log_interval = traj_cfg.get("frequency")
+        traj_file    = traj_cfg.get("file")
+
+        # Extract configuration output parameters
+        conf_cfg   = cfg.get("configuration_output", {})
+        init_file  = conf_cfg.get("initial")
+        final_file = conf_cfg.get("final")
+
+        # Extract control parameters & distortions
+        ctrl        = cfg["control_parameters"]
+        max_disp    = ctrl["maximum_displacement"]
+        distortions = ctrl["test_area_distortion"]
+        d1, d2      = distortions[0], distortions[1]
+        sx1, sy1, sz1 = d1["sx"], d1["sy"], d1["sz"]
+        sx2, sy2, sz2 = d2["sx"], d2["sy"], d2["sz"]
+
+        # Initial stats & initial config dump
         print(f"Initial # molecules: {len(box._molecules)}")
         print(f"Initial E_pot:        {box.total_epot:.5f}")
-
-        T = 0.8
-
+        if init_file:
+            os.makedirs(os.path.dirname(init_file), exist_ok=True)
+            box.write_XYZ(init_file)
+
+        # Set up Metropolis MC
         orchestrator.setup_simulation(
             MetropolisMC,
             T=T,
             log_energy=True
         )
-
-        # Run the full simulation as configured
-        orchestrator.run_simulation()
-
-        # Final results
+        mc = orchestrator.simulation
+        setattr(mc, "b", max_disp)
+
+        # Prepare results directory & log filename with parameters
+        results_dir = "results"
+        os.makedirs(results_dir, exist_ok=True)
+        param_str = f"steps{n_pr}_d1({sx1},{sy1},{sz1})"
+        log_fname = os.path.join(results_dir, f"result_{param_str}.log")
+
+        # Open log and write parameters header
+        exp_s1, exp_s2 = [], []
+        with open(log_fname, "w") as f_log:
+            f_log.write("# Simulation Parameters\n")
+            f_log.write(f"# Box dimensions: {Lx} x {Ly} x {Lz}\n")
+            f_log.write(f"# Temperature: {T}\n")
+            f_log.write(f"# Total steps: {n_pr}\n")
+            f_log.write(f"# Log interval: {log_interval}\n")
+            f_log.write(f"# Max displacement: {max_disp}\n")
+            f_log.write(f"# Distortions: d1=({sx1},{sy1},{sz1}), d2=({sx2},{sy2},{sz2})\n")
+            f_log.write(f"# Reset sampling at: {sorted(reset_points)}\n")
+            f_log.write(f"# Console freq: {console_freq}\n")
+            f_log.write(f"# Trajectory freq/file: {log_interval}/{traj_file}\n")
+            f_log.write("\n# step   E_pot     acc     w1        avg1      gamma1    w2        avg2      gamma2\n")
+
+            # Monte Carlo loop with distortion logging
+            for step in range(1, n_pr + 1):
+                acc  = mc.step()
+                Epot = box.total_epot
+
+                # Reset accumulators if requested
+                if step in reset_points:
+                    exp_s1.clear()
+                    exp_s2.clear()
+
+                # Compute distortions each MC step
+                dU1, dA1 = compute_distortion(box, sx1, sy1, sz1)
+                dU2, dA2 = compute_distortion(box, sx2, sy2, sz2)
+
+                # Boltzmann weights
+                w1 = np.exp(-dU1 / T)
+                w2 = np.exp(-dU2 / T)
+                exp_s1.append(w1)
+                exp_s2.append(w2)
+
+                # Console output
+                if console_freq and (step % console_freq == 0):
+                    print(f"Step {step}: E_pot={Epot:.3f}, acc={acc:.3f}")
+
+                # Trajectory dump
+                if traj_file and (step % log_interval == 0):
+                    os.makedirs(os.path.dirname(traj_file), exist_ok=True)
+                    box.write_XYZ(traj_file)
+
+                # Logging to file at intervals
+                if step % log_interval == 0:
+                    avg1 = np.mean(exp_s1)
+                    avg2 = np.mean(exp_s2)
+                    gamma1 = -T * np.log(avg1) / dA1
+                    gamma2 = -T * np.log(avg2) / dA2
+
+                    f_log.write(
+                        f"{step:6d}  {Epot:8.3f}  {acc:5.3f}  "
+                        f"{w1:7.3e}  {avg1:7.3e}  {gamma1:8.3f}  "
+                        f"{w2:7.3e}  {avg2:7.3e}  {gamma2:8.3f}\n"
+                    )
+
+            # Final surface-tension summary
+            f_log.write("\n")
+            f_log.write(f"gamma(area increase) = {gamma1:.6f}\n")
+            f_log.write(f"gamma(area decrease) = {gamma2:.6f}\n")
+
+        # Final configuration dump
+        if final_file:
+            os.makedirs(os.path.dirname(final_file), exist_ok=True)
+            box.write_XYZ(final_file)
+
+        # Final console output
         print("\n=== Final Results ===")
         print(f"Final E_pot:        {box.total_epot:.5f}")
         print(f"Total # molecules:  {len(box._molecules)}")
-
+
+
     except FileNotFoundError as e:
         print(f"Error: {e}", file=sys.stderr)
         sys.exit(1)
@@ -57,12 +165,13 @@ def run_with_orchestrator(config_file: str):
 
 
 def main():
-    # Check if a JSON config file was provided as command line argument
     if len(sys.argv) == 2 and sys.argv[1].endswith('.json'):
-        config_file = sys.argv[1]
-        print(f"Running with JSON configuration: {config_file}")
-        run_with_orchestrator(config_file)
+        print(f"Running with JSON configuration: {sys.argv[1]}")
+        run_with_orchestrator(sys.argv[1])
+    else:
+        print("Usage: python main.py <config.json>", file=sys.stderr)
+        sys.exit(1)
 
 
 if __name__ == "__main__":
-    main()
+    main()

src/ljts/box.py

@@ -1,5 +1,5 @@
 import numpy as np
-from ljts.molecule import Molecule
+from src.ljts.molecule import Molecule
 from collections import defaultdict
 import itertools
 
@@ -14,7 +14,7 @@ def __init__(self, len_x, len_y, len_z, den_liq=None, den_vap=None, potential=No
         self.potential = potential
 
         if den_liq is not None and den_vap is not None:
-            self._populate_box(den_liq, den_vap)
+            self.populate_box(den_liq, den_vap)
 
         self._total_Epot = 0.0
         self.total_potential_energy()
@@ -126,3 +126,7 @@ def num_molecules(self) -> int:
     @property
     def box_size(self):
         return self._box_size
+
+    @property
+    def total_Epot(self):
+        return self._total_Epot

src/ljts/distortion.py

@@ -1,29 +1,33 @@
 import numpy as np
 
-def compute_distortion(molecules, box_size, potential, sx, sy, sz):
+def compute_distortion(box, sx, sy, sz):
     """
     Compute the change in potential energy (delta U) and interface area (delta A)
     for a small, volume-conserving distortion of the simulation box and coordinates.
     """
+    molecules = box.get_molecules
+    potential = box.potential
+    box_size = box.box_size
+
+    # Ensure volume-conserving distortions
+    if not (np.isclose(sx * sy * sz, 1.0)):
+        raise ValueError("Distortions must be volume-conserving (sx * sy * sz = 1.0)")
+
     # Undistorted energy
-    E0 = 0.0
-    N = len(molecules)
-    for i in range(N):
-        for j in range(i + 1, N):
-            p1 = molecules[i].position
-            p2 = molecules[j].position
-            E0 += potential.potential_energy(p1, p2, box_size)
+    E0 = box.total_Epot
 
     # Distorted box size and scale matrix
     new_box = box_size * np.array([sx, sy, sz])
     scale_matrix = np.diag([sx, sy, sz])
 
     # Distorted energy
     E1 = 0.0
+    N = len(molecules)
     for i in range(N):
         for j in range(i + 1, N):
-            p1 = scale_matrix @ molecules[i].position
-            p2 = scale_matrix @ molecules[j].position
+            center = box_size / 2.0
+            p1 = scale_matrix @ (molecules[i].position - center) + center
+            p2 = scale_matrix @ (molecules[j].position - center) + center
             E1 += potential.potential_energy(p1, p2, new_box)
 
     delta_U = E1 - E0