Coverage for examples/structure_factor.py: 98%

1""" Example of computing the structure factor of a Lennard-Jones liquid. 

2 S(𝐪) = 1/N * |sum_{i=1}^{N} exp(-i𝐪•𝐫)|^2 

3 where 𝐪 is the wave vector, 𝐫 is the position of the particles, and N is the number of particles. 

4 The 𝐪 vectors are given by 𝐪 = 2π𝐧/L, where 𝐧 is a vector of integers and L is the box size. 

5""" 

6import matplotlib.pyplot as plt 

7import numpy as np 

8 

9import gamdpy as gp 

10 

11# Setup simulation of single-component Lennard-Jones liquid 

12temperature: float = 2.0 

13density: float = 0.973 

14configuration = gp.Configuration(D=3) 

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

16configuration['m'] = 1.0 

17configuration.randomize_velocities(temperature=1.44) 

18pair_func = gp.apply_shifted_force_cutoff(gp.LJ_12_6_sigma_epsilon) 

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

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

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

22runtime_actions = [gp.TrajectorySaver(), 

23 gp.ScalarSaver(), 

24 gp.MomentumReset(100)] 

25 

26 

27sim = gp.Simulation(configuration, pair_potential, integrator, runtime_actions, 

28 num_timeblocks=16, steps_per_timeblock=512, 

29 storage='memory') 

30 

31print("Equilibration run") 

32sim.run() 

33 

34print("Production run") 

35q_max = 18.0 

36calc_struct_fact = gp.CalculatorStructureFactor(configuration, backend='GPU') 

37calc_struct_fact.generate_q_vectors(q_max=q_max) 

38for block in sim.run_timeblocks(): 

39 print(sim.status(per_particle=True)) 

40 calc_struct_fact.update() 

41print(sim.summary()) 

42 

43struct_fact = calc_struct_fact.read(bins=100) 

44q = struct_fact['|q|'] 

45S = struct_fact['S(|q|)'] 

46 

47# Find maximum of S(q) and corresponding q 

48max_S_raw = np.max(S) 

49max_q_raw = q[np.argmax(S)] 

50# Fit polynomial to two nearest points to 

51n_range = 2 

52idx_max = np.argmax(S) 

53idx_fit = slice(idx_max-n_range, idx_max+n_range) 

54fit = np.polyfit(q[idx_fit], S[idx_fit], 2) 

55a, b, c = fit 

56max_q = -b / (2 * a) 

57max_S = a * max_q ** 2 + b * max_q + c 

58 

59plt.figure() 

60plt.title(f'Lennard-Jones liquid at $T={temperature}$ and $\\rho={density}$') 

61plt.plot(q, S, 'o') 

62x_fit = np.arange(q[idx_fit.start], q[idx_fit.stop], 0.01) 

63plt.plot(x_fit, a * x_fit ** 2 + b * x_fit + c, 'r--') 

64plt.text(max_q, max_S+1, f'Maximum: $S(q = {max_q:.3f}) = {max_S:.2f}$') 

65# Use S(q) in q→0 limit and estimate compressibility 

66plt.plot([0, 1], [S[0], S[0]], 'k--') 

67kappa_T = S[0]/(temperature*density) 

68plt.text(1, S[0]/2, r'$\kappa_T=\frac{S(q→0)}{\rho k_B T}=$' f'{kappa_T:.4f}') 

69plt.yscale('log') 

70plt.xlabel(r'Length of wave-vector, $|q|$ ($\sigma^{-1}$)') 

71plt.ylabel(r'Static structure factor, $S(|q|)$') 

72plt.ylim(1e-2, 10) 

73plt.xlim(0, max(q)) 

74if __name__ == "__main__": 

75 plt.show() 

76