Coverage for examples/analyze_structure.py: 89%

1""" Minimal example for calculing rdf from existing data  

2 

3 Usage: 

4 

5 analyze_structure filename 

6""" 

7 

8import gamdpy as gp 

9import numpy as np 

10import numba 

11import matplotlib.pyplot as plt 

12import sys 

13import pickle 

14 

15gp.select_gpu() 

16 

17argv = sys.argv.copy() 

18argv.pop(0) # remove scriptname 

19if __name__ == "__main__": 

20 if argv: 

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

22 else: 

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

24else: 

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

26 

27# Load existing data 

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

29# Read number of particles N and dimensions from data 

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

31 

32# Create configuration object 

33configuration = gp.Configuration(D=D, N=N) 

34configuration.simbox = gp.Orthorhombic(D, output['initial_configuration'].attrs['simbox_data']) 

35configuration.ptype = output['initial_configuration/ptype'] 

36configuration.copy_to_device() 

37# Call the rdf calculator 

38calc_rdf = gp.CalculatorRadialDistribution(configuration, bins=300) 

39 

40# NOTE: the structure of the block is (outer_block, inner_steps, pos&img, npart, dimensions) 

41# the zero is to select the position array and discard images 

42positions = output['trajectory_saver/positions'][:,:,:,:] 

43positions = positions.reshape(nblocks*nconfs,N,D) 

44# Loop over saved configurations 

45for pos in positions[nconfs-1::nconfs]: 

46 configuration['r'] = pos 

47 configuration.copy_to_device() 

48 calc_rdf.update() 

49 

50rdf_data = calc_rdf.read() 

51with open(filename+'_rdf.pkl', 'wb') as f: 

52 pickle.dump(rdf_data, f) 

53print(f"Wrote: {filename+'_rdf.pkl'}") 

54 

55num_types = rdf_data['rdf_ptype'].shape[1] 

56plt.figure(figsize=(8, 4)) 

57for i in range(num_types): 

58 for j in range(i, num_types): 

59 rdf_ij = np.mean(rdf_data['rdf_ptype'][:,i,j,:], axis=0) 

60 plt.plot(rdf_data['distances'], rdf_ij, label=f'{i}-{j}') 

61if num_types > 1: 

62 plt.legend() 

63plt.title(filename) 

64plt.xlabel('Distance') 

65plt.ylabel('Radial Distribution Function') 

66plt.savefig(filename+'_rdf.pdf') 

67print(f"Wrote: {filename+'_rdf.pdf'}") 

68if __name__ == "__main__": 

69 plt.show() 

70