Coverage for examples/widoms_particle_insertion.py: 100%

1""" Widom's particle insertion method for calculating the chemical potential 

2 

3In this example we use the Widom's particle insertion method 

4to calculate the chemical potential of a Lennard-Jones fluid. 

5 

6The excess chemical potential, μᵉˣ, of a not to dense fluid can be computed as: 

7 

8 μᵉˣ = -kT ln 〈exp(-ΔU/kT)〉 

9 

10where ΔU is the energy difference between the system with and without a ghost particle. 

11And 〈...〉 is an average for all possible positions of the ghost particle  

12(sometimes writte as an integral over space). 

13""" 

14 

15import numpy as np 

16 

17import gamdpy as gp 

18 

19# Setup configuration: FCC Lattice 

20rho = 0.4 # Number density 

21configuration = gp.Configuration(D=3) 

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

23configuration['m'] = 1.0 

24configuration.randomize_velocities(temperature=2.0) 

25 

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

27pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

28sig, eps, cut = 1.0, 1.0, 2.5 

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

30 

31# Setup integrator: NVT 

32temperature = 1.0 

33integrator = gp.integrators.NVT(temperature=temperature, tau=0.2, dt=0.005) 

34 

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

36runtime_actions = [gp.TrajectorySaver(), 

37 gp.ScalarSaver(), 

38 gp.MomentumReset(100)] 

39 

40# Setup Simulation 

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

42 num_timeblocks=20, steps_per_timeblock=1024, 

43 storage='memory') 

44 

45# Equilibrate the system 

46sim.run() 

47 

48# Setup the Widom's particle insertion calculator 

49num_ghost_particles = 500_000 

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

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

52print('Production run:') 

53for block in sim.run_timeblocks(): 

54 calc_widom.update() 

55 print('.', end='', flush=True) 

56 

57calc_widom_data = calc_widom.read() 

58print(f"\nExcess chemical potential: {calc_widom_data['chemical_potential']}") 

59 

60# Error estimation assuming that timeblocks are statistically independent 

61mu = calc_widom_data['chemical_potential'] 

62sigma = np.std(calc_widom_data['chemical_potentials']) / np.sqrt(len(calc_widom_data['chemical_potentials'])) 

63print(f"95 % confidence interval: [{mu - 1.96*sigma}, {mu + 1.96*sigma}]") 

64