@@ -96,53 +96,69 @@ chosen to make the calculation faster.
9696
9797 import numpy as np
9898 from MonteCarlo import MonteCarlo
99-
10099 from scipy import constants as cst
101100 from pint import UnitRegistry
101+ from reader import import_data
102+
103+ # Initialize the unit registry
102104 ureg = UnitRegistry()
103105
104106 # Constants
105- kB = cst.Boltzmann* ureg.J/ ureg.kelvin # boltzman constant
106- Na = cst.Avogadro/ ureg.mole # avogadro
107- R = kB* Na # gas constant
107+ kB = cst.Boltzmann * ureg.J / ureg.kelvin # Boltzmann constant
108+ Na = cst.Avogadro / ureg.mole # Avogadro's number
109+ R = kB * Na # Gas constant
108110
109111 # Parameters taken from Wood1957
110- tau = 2 # ratio between volume / reduced volume
111- epsilon = (119.76 * ureg.kelvin* kB * Na).to(ureg.kcal/ ureg.mol) # kcal/mol
112- r_star = 3.822 * ureg.angstrom # angstrom
113- sigma = r_star / 2 ** (1 / 6 ) # angstrom
114- N_atom = 50 # no units
115- m_argon = 39.948 * ureg.gram/ ureg.mol # g/mol
116- T = 328.15 * ureg.degK # 328 K or 55°C
112+ tau = 2 # Ratio between volume / reduced volume
113+ epsilon = (119.76 * ureg.kelvin * kB * Na).to(ureg.kcal / ureg.mol) # kcal/mol
114+ r_star = 3.822 * ureg.angstrom # Angstrom
115+ sigma = r_star / 2 ** (1 / 6 ) # Angstrom
116+ N_atom = 50 # Number of atoms
117+ m_argon = 39.948 * ureg.gram / ureg.mol # g/mol
118+ T = 328.15 * ureg.degK # 328 K or 55°C
117119 volume_star = r_star** 3 * Na * 2 ** (- 0.5 )
118- cut_off = sigma* 2.5 # angstrom
119- displace_mc = sigma/ 5 # angstrom
120- volume = N_atom* volume_star* tau/ Na # angstrom**3
121- box_size = volume** (1 / 3 ) # angstrom
122-
123- mc = MonteCarlo(maximum_steps = 15000 ,
120+ cut_off = sigma * 2.5 # Angstrom
121+ displace_mc = sigma / 5 # Angstrom
122+ volume = N_atom * volume_star * tau / Na # Angstrom^3
123+ box_size = volume** (1 / 3 ) # Angstrom
124+
125+ # Initialize and run the Monte Carlo simulation
126+ mc = MonteCarlo(
127+ maximum_steps = 15000 ,
124128 dumping_period = 1000 ,
125129 thermo_period = 1000 ,
126- thermo_outputs = " Epot-press" ,
130+ thermo_outputs = " Epot-press" ,
127131 neighbor = 50 ,
128132 number_atoms = [N_atom],
129133 epsilon = [epsilon.magnitude],
130134 sigma = [sigma.magnitude],
131135 atom_mass = [m_argon.magnitude],
132136 box_dimensions = [box_size.magnitude, box_size.magnitude, box_size.magnitude],
133- displace_mc = displace_mc.magnitude,
134- desired_temperature = T.magnitude,
135- cut_off = cut_off.magnitude,
136- data_folder = " Outputs/" ,
137- )
137+ displace_mc = displace_mc.magnitude,
138+ desired_temperature = T.magnitude,
139+ cut_off = cut_off.magnitude,
140+ data_folder = " Outputs/" ,
141+ )
138142 mc.run()
139143
140- # Import the data and calculate p V / R T
141- # output = np.mean(np.loadtxt("Outputs/pressure.dat")[:,1][10:])
142- # pressure = (output*ureg.atm).to(ureg.pascal)
143- # volume = (volume_star * tau / Na).to(ureg.meter**3)
144- # pV_over_RT = np.round((pressure * volume / (R * T) * Na).magnitude,2)
145- # print("p v / R T =", pV_over_RT, " --- (The expected value from Wood1957 is 1.5)")
144+ # Test function using pytest
145+ def test_pV_over_RT ():
146+ # Import the data and calculate pV / RT
147+ pressure_data = np.array(import_data(" Outputs/simulation.log" )[1 :])[0 ][:,2 ][10 :] # Skip initial values
148+ pressure_atm = np.mean(pressure_data) # atm
149+ pressure_Pa = (pressure_atm * ureg.atm).to(ureg.pascal) # Pa
150+ volume = (volume_star * tau / Na).to(ureg.meter** 3 ) # m3
151+ pV_over_RT = (pressure_Pa * volume / (R * T) * Na).magnitude
152+
153+ # Assert that pV_over_RT is close to 1.5
154+ assert np.isclose(pV_over_RT, 1.5 , atol = 1.0 ), f " Test failed: pV/Rt = { pV_over_RT} , expected close to 1.5 "
155+ print (f " pV/RT = { pV_over_RT:.2f } --- (The expected value from Wood1957 is 1.5) " )
156+
157+ # If the script is run directly, execute the tests
158+ if __name__ == " __main__" :
159+ import pytest
160+ # Run pytest programmatically
161+ pytest.main([" -s" , __file__ ])
146162
147163 .. label :: end_test_7a_class
148164
@@ -154,4 +170,4 @@ of 1.5 by Wood and Parker for :math:`\tau = V/V^* = 2` (see Fig. 4 in Ref. :cite
154170 (...)
155171 p v / R T = 1.56 --- (The expected value from Wood1957 is 1.5)
156172
157- The exact value will varie from one simulation to the other due to noise.
173+ The exact value will vary from one simulation to the other due to noise.
0 commit comments