Merge pull request #9 from mdcourse/improve-benchmark-with-LAMMPS

improved benchmark

integration_tests/benchmark/configurationA/README.md

@@ -1,11 +1,16 @@
-# Configuration A
+# Configuration A - Particles in the NVT ensemble
 
-Particle in the NVT ensemble.
+## Parameters
 
-nmb_1= 50  # Define atom number
-sig_1 = 3 * ureg.angstrom  # Define LJ parameters (sigma)
-eps_1 = 0.1 * ureg.kcal/ureg.mol  # Define LJ parameters (epsilon)
-mss_1 = 10 * ureg.gram/ureg.mol  # Define atom mass        
-L = 20 * ureg.angstrom  # Define box size
-rc = 2.5 * sig_1  # Define cut_off
-T = 300 * ureg.kelvin  # Pick the desired temperature
+nmb_1= 50  # Define atom number  
+sig_1 = 3 * ureg.angstrom  # Define LJ parameters (sigma)  
+eps_1 = 0.1 * ureg.kcal/ureg.mol  # Define LJ parameters (epsilon)  
+mss_1 = 10 * ureg.gram/ureg.mol  # Define atom mass          
+L = 20 * ureg.angstrom  # Define box size  
+rc = 2.5 * sig_1  # Define cut_off  
+T = 300 * ureg.kelvin  # Pick the desired temperature  
+
+## Measurements
+
+Epot = -3.532 ± 0.021 kcal/mol
+press = 298.069 ± 0.832 atm

integration_tests/benchmark/configurationA/nvt.lmp

@@ -1,6 +1,8 @@
-variable dump equal 5000
-variable thermo equal 5000
-variable steps equal 100000
+variable dump equal 50000
+variable thermo equal 50000
+variable steps equal 1000000
+variable eqs equal 100000
+variable av equal ${dump}/10
 
 variable nmb_1 equal 50  # Define atom number
 variable sig_1 equal 3  # Define LJ parameters (sigma)
@@ -34,16 +36,17 @@ timestep 0.25
 thermo ${thermo}
 dump mydmp all custom ${dump} dump.lammpstrj id type x y z vx vy vz
 
-run ${steps}  # equilibration
+run ${eqs}  # equilibration
 
 variable Epot equal pe
 variable Ekin equal ke
 variable Etot equal v_Epot+v_Ekin
 variable pressure equal press
 variable temperature equal temp
-fix myat1 all ave/time ${dump} 1 ${dump} v_Epot file Epot.dat
-fix myat2 all ave/time ${dump} 1 ${dump} v_Ekin file Ekin.dat
-fix myat3 all ave/time ${dump} 1 ${dump} v_Etot file Etot.dat
-fix myat4 all ave/time ${dump} 1 ${dump} v_pressure file pressure.dat
-fix myat5 all ave/time ${dump} 1 ${dump} v_temperature file temperature.dat
+fix myat1 all ave/time 10 ${av} ${dump} v_Epot file Epot.dat
+fix myat2 all ave/time 10 ${av} ${dump} v_Ekin file Ekin.dat
+fix myat3 all ave/time 10 ${av} ${dump} v_Etot file Etot.dat
+fix myat4 all ave/time 10 ${av} ${dump} v_pressure file pressure.dat
+fix myat5 all ave/time 10 ${av} ${dump} v_temperature file temperature.dat
+
 run ${steps}

integration_tests/benchmark/configurationA/update_README.py

@@ -0,0 +1,44 @@
+import numpy as np
+import re
+
+# Function to calculate the average and standard error of the second column in a data file
+def calculate_statistics(filename):
+    try:
+        data = np.loadtxt(filename, usecols=1)  # Load the second column
+        mean = np.mean(data)
+        std_error = np.std(data, ddof=1) / np.sqrt(len(data))  # Standard error of the mean
+        return mean, std_error
+    except Exception as e:
+        print(f"Error reading {filename}: {e}")
+        return None, None
+
+# Filenames
+config_file = "README.md"  # Your configuration file
+energy_file = "Epot.dat"
+pressure_file = "pressure.dat"
+
+# Calculate statistics
+average_epot, error_epot = calculate_statistics(energy_file)
+average_press, error_press = calculate_statistics(pressure_file)
+
+if average_epot is None or average_press is None:
+    print("Failed to calculate averages or errors. Exiting.")
+    exit(1)
+
+# Update the configuration file
+try:
+    with open(config_file, 'r') as file:
+        config_lines = file.readlines()
+
+    with open(config_file, 'w') as file:
+        for line in config_lines:
+            if re.match(r"^Epot\s*=", line):
+                file.write(f"Epot = {average_epot:.3f} ± {error_epot:.3f} kcal/mol\n")
+            elif re.match(r"^press\s*=", line):
+                file.write(f"press = {average_press:.3f} ± {error_press:.3f} atm\n")
+            else:
+                file.write(line)
+
+    print(f"Updated {config_file} with Epot = {average_epot:.3f} ± {error_epot:.3f} kcal/mol and press = {average_press:.3f} ± {error_press:.3f} atm")
+except Exception as e:
+    print(f"Error updating {config_file}: {e}")

integration_tests/test_monte_carlo_move.py

@@ -28,12 +28,12 @@ def setUp(self):
         L = 20 * ureg.angstrom  # Define box size
         rc = 2.5 * sig_1  # Define cut_off
         T = 300 * ureg.kelvin  # Pick the desired temperature
-        displace_mc = sig_1/2  # choose the displace_mc
+        displace_mc = sig_1/4  # choose the displace_mc
 
         # Initialize the MonteCarlo object
         self.mc = MonteCarlo(
             ureg = ureg,
-            maximum_steps = 10000,
+            maximum_steps = 50000,
             thermo_period = 1000,
             dumping_period = 1000,
             number_atoms = [nmb_1],