Coverage for examples/thermodynamics.py: 97%

1""" Investigation of thermodynamic properties 

2 

3This example show how thermodynamic data can be extracted 

4using the `extract_scalars` function from the `gamdpy` package. 

5 

6The script runs a NVT simulation of a Lennard-Jones crystal. 

7The potential energy, virial, and kinetic energy are extracted 

8from the simulation output and the mean values are printed. 

9Derived quantities such as kinetic temperature and pressure are computed. 

10The energy trajectory is plotted and the error estimate of potential energy 

11is estimated using the blocking method. 

12 

13""" 

14 

15import numpy as np 

16import matplotlib.pyplot as plt 

17from scipy import stats 

18 

19import gamdpy as gp 

20 

21 

22############################### 

23# Setup and run simulation # 

24############################### 

25 

26# Setup configuration: FCC Lattice 

27configuration = gp.Configuration(D=3) 

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

29configuration['m'] = 1.0 

30configuration.randomize_velocities(temperature=0.7) 

31 

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

33pair_func = gp.apply_shifted_potential_cutoff(gp.LJ_12_6_sigma_epsilon) 

34sig, eps, cut = 1.0, 1.0, 2.5 

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

36 

37# Setup integrator: NVT 

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

39 

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

41runtime_actions = [gp.TrajectorySaver(), 

42 gp.ScalarSaver(4), 

43 gp.MomentumReset(100)] 

44 

45# Setup Simulation. 

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

47 num_timeblocks=64, steps_per_timeblock=2048, 

48 storage='memory') 

49 

50# Run simulation 

51sim.run() 

52 

53# Basic information of NVT simulation 

54N = sim.configuration.N # Number of particles 

55T = sim.integrator.temperature # Temperature 

56V = configuration.get_volume() 

57rho = N / V # Density 

58 

59 

60######################## 

61# Extracting scalars # 

62######################## 

63 

64# Extract potential energy (U), virial (W), and kinetic energy (K) 

65# We use first_block=1 to skip the initial "equilibration" block. 

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

67 

68# Print mean values 

69print(f"Mean potential energy per particle: {np.mean(U) / N}") 

70print(f"Mean virial per particle: {np.mean(W) / N}") 

71print(f"Mean kinetic energy per particle: {np.mean(K) / N}") 

72 

73 

74######################## 

75# Derived quantities # 

76######################## 

77 

78# Time 

79dt = sim.integrator.dt 

80time = np.arange(len(U)) * dt * sim.output['scalar_saver'].attrs["steps_between_output"] 

81print(f"Total time of analysed trajectory: {time[-1]}") 

82 

83# Compute kinetic temperature 

84D = sim.configuration.D # dimension of space 

85dof = D * N - D # degrees of freedom 

86T_kin = 2 * K / dof 

87print(f"Mean kinetic temperature: {np.mean(T_kin)}") 

88 

89# Compute instantaneous pressure 

90P = rho * T_kin + W / V 

91print(f"Mean pressure: {np.mean(P)}") 

92 

93 

94################################ 

95# Plotting energy trajectory # 

96################################ 

97 

98# Plot potential energy per particle 

99plt.figure() 

100plt.plot(time, U / N, label='Potential energy') 

101plt.axhline(np.mean(U) / N, color='k', linestyle='--', label='Mean') 

102plt.xlabel('Time') 

103plt.ylabel('Potential energy') 

104plt.legend() 

105if __name__ == "__main__": 

106 plt.show() 

107 

108 

109######################################## 

110# Statistical analysis using blocking # 

111######################################## 

112 

113def error_estimate_by_blocking(data, blocks, confidence=0.95): 

114 """ Estimate error on mean using block averaging. """ 

115 block_length = len(data) // blocks 

116 data = data[:blocks * block_length].reshape(blocks, block_length) 

117 block_means = np.mean(data, axis=1) # Mean of each block 

118 sem = np.std(block_means) / np.sqrt(blocks) # Standard error of the mean 

119 t_value = stats.t.ppf((1 + confidence) / 2, blocks - 1) # t-value for confidence interval 

120 return t_value * sem # Error estimate 

121 

122 

123# Find optimal block size 

124test_num_blocks = list(set(np.logspace(0, np.log10(len(U) / 4), num=32, dtype=int))) 

125test_num_blocks.sort() 

126test_num_blocks = test_num_blocks[2:] 

127errors = [error_estimate_by_blocking(U/N, blocks) for blocks in test_num_blocks] 

128 

129# Good choice for number of blocks (see figure below) 

130# If num_blocks is too small, there is large uncertainty on the error estimate. 

131# If num_blocks is too large, blocks are not statistical independent. 

132num_blocks = 64 

133error = error_estimate_by_blocking(U/N, num_blocks) 

134 

135plt.figure() 

136plt.title('Error estimate of potential energy per particle') 

137plt.plot(test_num_blocks, errors, 'o') 

138plt.plot(num_blocks, error, 'ro', label='Estimated error') 

139plt.xscale('log') 

140plt.xlim(1, None) 

141plt.ylim(0, None) 

142plt.xlabel('Number of blocks') 

143plt.ylabel('Estimated error (95% confidence interval)') 

144if __name__ == "__main__": 

145 plt.show() 

146 

147print(f'Potential energy per particle {np.mean(U) / N:.4f} ± {error:.4f} (95% confidence interval)')