Coverage for gamdpy/tools/calc_dynamics.py: 54%
96 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import math
2import numpy as np
3import matplotlib.pyplot as plt
6def calc_dynamics_(positions, images, ptype, simbox, block0, conf_index0, block1, conf_index1, time_index, msd, m4d, qvalues=None, Fs=None):
7 """
8 Calculate contribution to dynamical properties from conf_index1 in block1 using conf_index0 in block0
9 as initial time, and add it at time_index in appropiate arrays.
11 TODO: Allow more than one q-value per type
12 """
13 dR = positions[block1, conf_index1, :, :] - positions[block0, conf_index0, :, :]
14 dR += (images[block1, conf_index1, :, :] - images[block0, conf_index0, :, :]) * simbox
15 for i in range(np.max(ptype) + 1):
16 dR_type = dR[ptype == i, :]
17 dR_i_sq = np.sum( dR_type**2, axis=1)
18 msd[time_index, i] += np.mean(dR_i_sq)
19 m4d[time_index, i] += np.mean(dR_i_sq ** 2)
20 Fs[time_index, i] += np.mean(np.cos(dR_type*qvalues[i]))
22 return msd, m4d, Fs
25def calc_dynamics(trajectory, first_block, qvalues=None):
26 """Compute dynamical properties from a saved trajectory HDF5 file.
28 This function processes blocks of saved configurations to evaluate time‐dependent
29 dynamical observables, including the mean square displacement (MSD), the non‐Gaussian
30 parameter (alpha2), and the self‐intermediate scattering function (Fs), for one or more
31 particle types.
33 Parameters
34 ----------
35 trajectory : h5py.File object in the gamdpy style
37 first_block : int
38 Index of the first block to use as the reference origin.
40 qvalues : float or array‐like of shape (num_types,), optional
41 Wavevector magnitudes at which to compute the self‐intermediate scattering
42 function Fs. If a single float is provided, it is broadcast to all particle types.
43 If None, Fs is not computed (remains zero).
45 Returns
46 -------
47 results : dict
48 Dictionary containing dynamcal data.
50 Examples
51 --------
52 For command‐line usage, see:
53 $ python -m gamdpy.tools.calc_dynamics -h
55 Usage within a Python script:
57 >>> import gamdpy as gp
58 >>> sim = gp.get_default_sim() # Replace with your simulation object
59 >>> for block in sim.run_timeblocks(): pass
60 >>> dynamics = gp.calc_dynamics(sim.output, first_block=0, qvalues=7.5)
61 >>> dynamics.keys()
62 dict_keys(['times', 'msd', 'alpha2', 'qvalues', 'Fs', 'count'])
64 """
65 ptype = trajectory['initial_configuration/ptype'][:].copy()
66 attributes = trajectory.attrs
68 simbox = trajectory['initial_configuration'].attrs['simbox_data'].copy()
69 num_types = np.max(ptype) + 1
70 if isinstance(qvalues, float):
71 qvalues = np.ones(num_types)*qvalues
72 num_blocks, conf_per_block, N, D = trajectory['trajectory_saver/positions'].shape
73 blocks = trajectory['trajectory_saver/positions'] # If picking out dataset in inner loop: Very slow!
74 images = trajectory['trajectory_saver/images']
76 #print(num_types, first_block, num_blocks, conf_per_block, _, N, D, qvalues)
77 if first_block > num_blocks - 1:
78 print("Warning [calc_dynamics] first_block greater than number of blocks. Remainder will be taken")
79 first_block = first_block % num_blocks # necessary to allow the pythonic idiom of negative indexes
81 extra_times = int(math.log2(num_blocks - first_block)) - 1
82 total_times = conf_per_block - 1 + extra_times
83 count = np.zeros((total_times, 1), dtype=np.int32)
84 msd = np.zeros((total_times, num_types))
85 m4d = np.zeros((total_times, num_types))
86 Fs = np.zeros((total_times, num_types))
88 times = attributes['dt'] * 2 ** np.arange(total_times)
90 for block in range(first_block, num_blocks):
91 for i in range(conf_per_block - 1):
92 count[i] += 1
93 calc_dynamics_(blocks, images, ptype, simbox, block, i + 1, block, 0, i, msd, m4d, qvalues, Fs)
95 # Compute times longer than blocks
96 for block in range(first_block, num_blocks):
97 for i in range(extra_times):
98 index = conf_per_block - 1 + i
99 other_block = block + 2 ** (i + 1)
100 # print(other_block, end=' ')
101 if other_block < num_blocks:
102 count[index] += 1
103 calc_dynamics_(blocks, images, ptype, simbox, other_block, 0, block, 0, index, msd, m4d, qvalues, Fs)
105 msd /= count
106 m4d /= count
107 Fs /= count
108 alpha2 = 3 * m4d / (5 * msd ** 2) - 1
109 return {'times': times, 'msd': msd, 'alpha2': alpha2, 'qvalues':qvalues, 'Fs':Fs, 'count': count}
112def create_msd_plot(dynamics, figsize=(8, 6)):
113 fig, axs = plt.subplots(1, 1, figsize=figsize)
114 for dyn in dynamics:
115 axs.loglog(dyn['times'], dyn['msd'], '.-', label=dyn['name'])
116 axs.set_xlabel('Time')
117 axs.set_ylabel('MSD')
118 axs.legend()
119 return fig, axs
121def create_alpha2_plot(dynamics, figsize=(8, 6)):
122 fig, axs = plt.subplots(1, 1, figsize=figsize)
123 for dyn in dynamics:
124 axs.semilogx(dyn['times'], dyn['alpha2'], '.-', label=dyn['name'])
125 axs.set_xlabel('Time')
126 axs.set_ylabel('alpha2')
127 axs.legend()
128 return fig, axs
131def main(argv: list = None) -> None:
132 """ Command line interface for calc_dynamics """
133 import sys
134 import h5py
136 help_message = """gamdpy: calc_dynamics
138Calculate and show the mean square displacement (MSD)
140Usage: python -m gamdpy.tools.calc_dynamics [options] <input filename> [<input filename> ...]
142Options:
143 -h, --help Print this help message and exit.
144 -f <int> First block to use. Default is 0.
145 -o <filename> Output filename. Default is no output.
147Example: python -m gamdpy.tools.calc_dynamics -f 4 -o msd.pdf LJ*.h5
148 """
150 if argv is None:
151 argv = sys.argv
153 # Print help if '-h' or '--help' is in argv or if no arguments are given
154 if '-h' in argv or '--help' in argv or len(argv) == 1:
155 print(help_message)
156 return
158 argv.pop(0) # remove name
160 first_block = 0
161 output_filename = ''
162 while argv[0][0] == '-':
163 if argv[0] == '-f':
164 argv.pop(0) # remove '-f'
165 first_block = int(argv.pop(0)) # read and remove parameter
166 if argv[0] == '-o':
167 argv.pop(0) # remove '-o'
168 output_filename = argv.pop(0) # read and remove parameter
170 # The rest should be filenames...
171 dynamics = []
172 for filename in argv:
173 print(filename, ':', end=' ')
174 with h5py.File(filename, "r") as f:
175 dynamics.append(calc_dynamics(f, first_block))
176 dynamics[-1]['name'] = filename[:-3]
178 fig, axs = create_msd_plot(dynamics)
179 if not output_filename == '':
180 plt.savefig(output_filename)
181 plt.show()
184if __name__ == '__main__':
185 main()