Coverage for examples/analyze_dynamics.py: 90%

1""" Analyze and plot dynamics computed from a .h5 file  

2 Usage: 

3 

4 analyze_dynamics filename 

5""" 

6 

7import matplotlib.pyplot as plt 

8import gamdpy as gp 

9import numpy as np 

10import sys 

11import pickle 

12 

13argv = sys.argv.copy() 

14argv.pop(0) # remove scriptname 

15if __name__ == "__main__": 

16 if argv: 

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

18 else: 

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

20else: 

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

22 

23# Load existing data 

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

25 

26dynamics = gp.tools.calc_dynamics(output, 0, qvalues=7.5) 

27with open(filename+'_dynamics.pkl', 'wb') as f: 

28 pickle.dump(dynamics, f) 

29print(f"Wrote: {filename+'_dynamics.pkl'}") 

30 

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

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

33axs[0].set_ylabel('MSD') 

34axs[1].set_ylabel('Non Gaussian parameter') 

35axs[2].set_ylabel('Intermediate scattering function') 

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

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

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

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

40 

41num_types = dynamics['msd'].shape[1] 

42for i in range(num_types): 

43 axs[0].loglog(dynamics['times'], dynamics['msd'][:,i], 'o--', label=f'{i}') 

44 axs[1].semilogx(dynamics['times'], dynamics['alpha2'][:,i], 'o--') 

45 axs[2].semilogx(dynamics['times'], dynamics['Fs'][:,i], 'o--', label=f'{i}, q={dynamics["qvalues"][i]}') 

46factor = np.array([1, 30]) 

47axs[0].loglog(dynamics['times'][0]*factor, np.max(dynamics['msd'][0,:])*factor**2, 'k--', alpha=0.5, label='Slope 2') 

48axs[0].loglog(dynamics['times'][-1]/factor, np.min(dynamics['msd'][-1,:])/factor, 'k-.', alpha=0.5, label='Slope 1') 

49 

50axs[0].set_title(filename) 

51axs[0].legend() 

52axs[2].legend() 

53fig.savefig(filename+'_dynamics.pdf') 

54print(f"Wrote: {filename+'_dynamics.pdf'}") 

55 

56if __name__ == "__main__": 

57 plt.show(block=True) 

58