Coverage for examples/analyze_thermodynamics.py: 90%

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 

6 Usage: 

7 

8 analyze_thermodynamics filename 

9""" 

10 

11import matplotlib.pyplot as plt 

12import gamdpy as gp 

13import numpy as np 

14import sys 

15 

16 

17max_plot_points = 100_000 

18 

19argv = sys.argv.copy() 

20argv.pop(0) # remove scriptname 

21if __name__ == "__main__": 

22 if argv: 

23 filename = argv.pop(0) # get filename (.h5 added by script) 

24 else: 

25 filename = 'Data/LJ_r0.973_T0.70_toread' # Used in testing 

26else: 

27 filename = 'Data/LJ_r0.973_T0.70_toread' 

28 

29output = gp.tools.TrajectoryIO(filename+'.h5').get_h5() 

30 

31nblocks, nconfs, N, D = output['trajectory_saver/positions'].shape 

32simbox = output['initial_configuration'].attrs['simbox_data'] 

33volume = np.prod(simbox) 

34rho = N/volume 

35 

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

37# first_block can be used to skip the initial "equilibration". 

38U, W, K = gp.extract_scalars(output, ['U', 'W', 'K'], first_block=0) 

39 

40mU = np.mean(U) 

41mW = np.mean(W) 

42mK = np.mean(K) 

43 

44# Hack to find parts of data not valid 

45#print(np.mean(K>0)) 

46 

47# Time 

48dt = output.attrs['dt'] 

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

50 

51# Compute mean kinetic temperature 

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

53T_kin = 2 * mK / dof 

54 

55# Compute mean pressure 

56P = rho * T_kin + mW / volume 

57 

58# Compute W-U correlations 

59dU = U - mU 

60dW = W - mW 

61gamma = np.dot(dW,dU)/np.dot(dU,dU) 

62R = np.dot(dW,dU)/(np.dot(dW,dW)*np.dot(dU,dU))**0.5 

63 

64# Plot  

65plotindex = range(len(U)) 

66#print(len(plotindex)) 

67if len(U)>max_plot_points: 

68 step = int(len(U)/max_plot_points+1) 

69 plotindex = plotindex[::step] 

70#print(len(plotindex)) 

71 

72title = f'N={N}, rho={rho:.3f}, Tkin={np.mean(T_kin):.3f}, P={np.mean(P):.3f}, R={R:.3f}, gamma={gamma:.3f}' 

73 

74fig, axs = plt.subplots(3, 1, figsize=(8, 9), sharex=True) 

75fig.subplots_adjust(hspace=0.00) # Remove vertical space between axes 

76axs[0].set_title(title) 

77axs[0].set_ylabel('U/N') 

78axs[1].set_ylabel('W/N') 

79axs[2].set_ylabel('K/N') 

80axs[2].set_xlabel('Time') 

81axs[0].grid(linestyle='--', alpha=0.5) 

82axs[1].grid(linestyle='--', alpha=0.5) 

83axs[2].grid(linestyle='--', alpha=0.5) 

84 

85label = f'mean: {mU/N:.3f} std: {np.std(U/N):.3f}' 

86axs[0].plot(time[plotindex], U[plotindex] / N, label=label) 

87axs[0].axhline(mU / N, color='k', linestyle='--') 

88axs[0].legend(loc= 'upper right') 

89 

90label = f'mean: {mW/N:.3f} std: {np.std(W/N):.3f}' 

91axs[1].plot(time[plotindex], W[plotindex] / N, label=label) 

92axs[1].axhline(mW / N, color='k', linestyle='--') 

93axs[1].legend(loc= 'upper right') 

94 

95label = f'mean: {mK/N:.3f} std: {np.std(K/N):.3f}' 

96axs[2].plot(time[plotindex], K[plotindex] / N, label=label) 

97axs[2].axhline(mK / N, color='k', linestyle='--') 

98axs[2].legend(loc= 'upper right') 

99 

100fig.savefig(filename+'_thermodynamics.pdf') 

101if __name__ == "__main__": 

102 plt.show(block=True)