Coverage for tests/test_widom.py: 100%

1""" Test Widom's particle insertion method """ 

2 

3import numpy as np 

4import matplotlib.pyplot as plt 

5 

6import gamdpy as gp 

7 

8def test_widom_insertion(): 

9 np.random.seed(0) 

10 

11 # Setup configuration: FCC Lattice 

12 configuration = gp.Configuration(D=3) 

13 configuration.make_lattice(gp.unit_cells.FCC, cells=[8, 8, 8], rho=0.4) 

14 configuration['m'] = 1.0 

15 configuration.randomize_velocities(temperature=2.0, seed=0) 

16 

17 # Setup pair potential: Single component 12-6 Lennard-Jones 

18 pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

19 sig, eps, cut = 1.0, 1.0, 2.5 

20 pair_pot = gp.PairPotential(pair_func, params=[sig, eps, cut], max_num_nbs=1000) 

21 

22 # Setup integrator: NVT 

23 temperature = 1.0 

24 integrator = gp.integrators.NVT(temperature=1.0, tau=0.2, dt=0.0) # dummy dt 

25 

26 # Setup runtime actions, i.e. actions performed during simulation of timeblocks 

27 runtime_actions = [gp.TrajectorySaver(), 

28 gp.ScalarSaver(), 

29 gp.MomentumReset(16)] 

30 

31 # Setup Simulation 

32 sim = gp.Simulation(configuration, pair_pot, integrator, runtime_actions, 

33 num_timeblocks=2, steps_per_timeblock=32, 

34 storage='memory') 

35 

36 # Setup the Widom's particle insertion calculator 

37 num_ghost_particles = 500_000 

38 ghost_positions = np.random.rand(num_ghost_particles, configuration.D) * configuration.simbox.get_lengths() 

39 calc_widom = gp.CalculatorWidomInsertion(sim.configuration, pair_pot, temperature, ghost_positions) 

40 print('Production run') 

41 for block in sim.run_timeblocks(): 

42 calc_widom.update() 

43 

44 calc_widom_data = calc_widom.read() 

45 

46 # Test if the excess chemical potential in in expected range 

47 mu = calc_widom_data['chemical_potential'] 

48 print(f"Excess chemical potential: {mu}") 

49 mu_expected = 0.4469328390273799 

50 print(f"Expected excess chemical potential: {mu_expected}") 

51 print(f"Error: {mu - mu_expected}") 

52 my_tol = 1e-4 

53 assert np.isclose(mu, mu_expected, rtol=my_tol), f"mu = {mu}, mu_expected = {mu_expected}" 

54 

55if __name__ == '__main__': # pragma: no cover 

56 test_widom_insertion()