added benchmark

Author: simongravelle

integration_tests/benchmark/configurationA/README.md

@@ -0,0 +1,11 @@
+# Configuration A
+
+Particle in the NVT ensemble.
+
+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

integration_tests/benchmark/configurationA/nvt.lmp

@@ -0,0 +1,49 @@
+variable dump equal 5000
+variable thermo equal 5000
+variable steps equal 100000
+
+variable nmb_1 equal 50  # Define atom number
+variable sig_1 equal 3  # Define LJ parameters (sigma)
+variable eps_1 equal 0.1  # Define LJ parameters (epsilon)
+variable mss_1 equal 10  # Define atom mass        
+variable L equal 20  # Define box size
+variable rc equal 2.5*${sig_1}  # Define cut_off (angstrom)
+variable T equal 300  # Pick the desired temperature (kelvin)
+
+units real
+dimension 3
+atom_style atomic
+pair_style lj/cut ${rc}
+boundary p p p
+
+#read_data twoparticle.data
+
+variable L2 equal ${L}/2
+region myreg block -${L2} ${L2} -${L2} ${L2} -${L2} ${L2}
+create_box 1 myreg
+create_atoms 1 random ${nmb_1} 32141 myreg
+
+mass 1 ${mss_1}
+pair_coeff 1 1 ${eps_1} ${sig_1}
+
+velocity all create ${T} 4928459
+fix mynve all nve
+fix myber all temp/berendsen 300 300 100
+timestep 0.25
+
+thermo ${thermo}
+dump mydmp all custom ${dump} dump.lammpstrj id type x y z vx vy vz
+
+run ${steps}  # 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
+run ${steps}

integration_tests/test_monte_carlo_move.py

@@ -17,30 +17,31 @@ def setUp(self):
 
         ureg = UnitRegistry()
 
+        # Configuration A
         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
-        displace_mc = sig_1/4  # choose the displace_mc
+        displace_mc = sig_1/2  # choose the displace_mc
 
         # Initialize the MonteCarlo object
         self.mc = MonteCarlo(
             ureg = ureg,
-            maximum_steps=100,
-            thermo_period=10,
-            dumping_period=10,
-            number_atoms=[nmb_1],
-            epsilon=[eps_1], # kcal/mol
-            sigma=[sig_1], # A
-            atom_mass=[mss_1], # g/mol
-            box_dimensions=[L, L, L], # A
-            cut_off=rc,
-            thermo_outputs="Epot-press",
-            desired_temperature=T, # K
-            neighbor=20,
+            maximum_steps = 10000,
+            thermo_period = 1000,
+            dumping_period = 1000,
+            number_atoms = [nmb_1],
+            epsilon = [eps_1],
+            sigma = [sig_1],
+            atom_mass = [mss_1],
+            box_dimensions = [L, L, L],
+            cut_off = rc,
+            thermo_outputs = "Epot-press",
+            desired_temperature = T,
+            neighbor = 20,
             displace_mc = displace_mc,
         )
 

molecular_simulation_code/forces_utilities.py

@@ -61,6 +61,6 @@ def compute_force_matrix(neighbor_lists: List[np.ndarray], atoms_positions: np.n
 
             # Update force matrix
             force_matrix[atom_i, atom_j] += fij_xyz * rij_xyz / rij
-            force_matrix[atom_j, atom_i] += fij_xyz * rij_xyz / rij
+            force_matrix[atom_j, atom_i] -= fij_xyz * rij_xyz / rij
     return force_matrix
 

molecular_simulation_code/simulation_logger.py

@@ -60,6 +60,7 @@ def log_simulation_data(code):
                 pressure = compute_pressure(code.atoms_positions, code.box_mda,
                                               code.neighbor_lists, code.cross_coefficients,
                                               code.potential_type, code.desired_temperature)
+
                 press = pressure * code.ref_pressure  # Atm
                 logger.info(f"{code.step} {Epot.magnitude:.2f} {press.magnitude:.2f}")