Coverage for examples/consistency_NPT.py: 97%

1""" Consistency check for NPT integrators. 

2 

3Simulation of a Lennard-Jones liquid in the NPT ensemble. 

4It's possible to switch between Langevin and Atomic NPT integrators. 

5This script verifies that several thermodynamics quantity are properly produced. 

6 

7""" 

8 

9import numpy as np 

10 

11import gamdpy as gp 

12 

13# Here you can decide to use "NPT_Atomic" or "NPT_Langevin" 

14flag = "Atomic" 

15my_T, my_rho, my_p = 2.0, 0.754289412611, 4.7 # Pressure should be P=4.7 for T=2.0 at this density 

16 

17# Choose integrator 

18if flag=="Atomic": 

19 integrator = gp.integrators.NPT_Atomic (temperature=my_T, tau=0.4, pressure=my_p, tau_p=20, dt=0.001) 

20elif flag=="Langevin": 

21 gp.integrators.NPT_Langevin(temperature=my_T, pressure=my_p, alpha=0.1, alpha_baro=0.0001, mass_baro=0.0001, 

22 volume_velocity=0.0, barostatModeISO=True, boxFlucCoord=2, dt=0.001, seed=2023) 

23 

24# Setup configuration: FCC Lattice 

25configuration = gp.Configuration(D=3, compute_flags={'Vol':True}) 

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

27configuration['m'] = 1.0 

28configuration.randomize_velocities(temperature=my_T, seed=0) 

29 

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

31pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

32sig, eps, cut = 1.0, 1.0, 2.5 

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

34 

35# NVT equilibration for calculation of c_V and Thermal pressure coefficient 

36integratorNVT = gp.integrators.NVT(temperature=my_T, tau=0.2, dt=0.001) 

37 

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

39runtime_actions = [gp.TrajectorySaver(), 

40 gp.ScalarSaver(32), 

41 gp.MomentumReset(100)] 

42 

43sim = gp.Simulation(configuration, pair_pot, integratorNVT, runtime_actions, 

44 num_timeblocks=8, steps_per_timeblock=16384, 

45 storage='memory') 

46# Equilibration run 

47print("Running NVT simulation at (\\rho, T) = ({my_rho},{my_T})") 

48print("Equilibration") 

49sim.run() 

50# Data run 

51print("Data run") 

52sim.run() 

53U, W, K, Vol = gp.extract_scalars(sim.output, ['U', 'W', 'K', 'Vol'], first_block=1) 

54# Full c_V (not excess) 

55c_V = np.std(U+K)**2/my_T**2/configuration.N 

56dU = U - np.mean(U) 

57dW = W - np.mean(W) 

58# Thermal pressure coefficient https://en.wikipedia.org/wiki/Thermal_pressure 

59# from Eq. 89 in http://glass.ruc.dk/pdf/articles/2008_J_Chem_Phys_129_184508.pdf 

60P_th = my_rho + np.mean(dW*dU)/my_T**2/np.mean(Vol) 

61print(f"NVT simulation at T={my_T} and \\rho={my_rho}") 

62print(f"Mean values of U, W and T_kin: {np.mean(U)/configuration.N} {np.mean(W)/configuration.N} {2*np.mean(K)/3/configuration.N}") 

63print(f"Standard dev of U, W and T_kin: {np.std(U)/configuration.N} {np.std(W)/configuration.N} {2*np.std(K)/3/configuration.N}") 

64print(f"Pressure (mean): {my_rho*(2*np.mean(K)/3+np.mean(W))/configuration.N}") 

65print(f"Thermal pressure coefficient: {P_th}") 

66print(f"Specific heat at constant volume: {c_V}") 

67print() 

68 

69# NPT Simulation  

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

71 num_timeblocks=8, steps_per_timeblock=16384, 

72 storage='memory') 

73sim.run() 

74sim.run() 

75U, W, K, Vol = gp.extract_scalars(sim.output, ['U', 'W', 'K', 'Vol'], first_block=1) 

76H = K + U + my_p * Vol # enthalpy H = Etot + PV  

77dH = H - np.mean(H) 

78dV = Vol - np.mean(Vol) 

79# Full c_P (not excess) 

80c_P = np.std(H)**2/my_T**2/configuration.N 

81print(f"NPT simulation at T={my_T} and p={my_p}") 

82print(f"Mean values of U, W and T_kin: {np.mean(U)/configuration.N} {np.mean(W)/configuration.N} {2*np.mean(K)/3/configuration.N}") 

83print(f"Standard dev of U, W and T_kin: {np.std(U)/configuration.N} {np.std(W)/configuration.N} {2*np.std(K)/3/configuration.N}") 

84print(f"Mean and std of Volume: {np.mean(Vol)} {np.std(Vol)}") 

85print(f"Pressure (mean): {(2*np.mean(K)/3+np.mean(W))/np.mean(Vol)}") 

86print(f"Density should be {my_rho} and is {configuration.N/np.mean(Vol)}") 

87K_T = my_T*np.mean(Vol)/np.std(Vol)**2 

88print(f"Isothermal bulk modulus K_T: {K_T}") 

89print() 

90print(f"Check 1: Thermal expansion coefficient") 

91beta_P = P_th/K_T # from https://en.wikipedia.org/wiki/Thermal_pressure (it's named \alpha on wiki) 

92print(f"Thermal expansion coefficient from NVT fluctuations \\beta_P: {beta_P}") 

93print(f"Thermal expansion coefficient form NPT fluctuations \\beta_P: {np.mean(dH*dV)/my_T**2/np.mean(Vol)}") # 1/V dV/dT|_P Allen-Tildesley Eq. 2.87 pg 53 

94print() 

95print(f"Check 2: Consistency relations in Appendix B of http://dx.doi.org/10.1063/1.467468") 

96print(f"Equation B3 <P_int> = P_ext + T<V^-1> : {np.mean((2*K/3+W)/Vol)} ?= {my_p + my_T*np.mean(1/Vol)}") 

97print(f"Equation B4 <P_int V> = P_exp <V> : {np.mean( 2*K/3+W )/configuration.N} ?= {my_p * np.mean(Vol)/configuration.N}") 

98print() 

99print(f"Check 3: Specific heat at constant pressure (this needs a longer run, use num_timeblocks=128)") 

100print(f"c_P from enthalpy fluctuations in the NPT ensamble : {c_P}") 

101# Verify that C_P = C_V + T V \beta_P**2 * K_T from problem 5.14 part c of Daniel V. Schroeder - An Introduction to Thermal Physics (1999) 

102print(f"c_P obtained from c_P = c_V + T V \\beta_P**2 * K_T: {c_V + my_T * np.mean(Vol)/configuration.N * beta_P**2 * K_T}" ) 

103print() 

104