Coverage for tests/test_structure_factor.py: 100%
92 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import numpy as np
2import matplotlib.pyplot as plt
4import gamdpy as gp
5from object_lib import pairpot_LJ as pair_potential
7# The following environment variable fixes a warning appearing during testing
8# The warning appears when the tbb library is installed with conda
9# https://github.com/numba/numba/issues/6350
10import os
11os.environ['NUMBA_THREADING_LAYER']='omp'
13def test_structure_factor():
14 # verbose = False
15 # plot = False
17 # Setup simulation of single-component Lennard-Jones liquid
18 temperature = 2.0
19 configuration = gp.Configuration(D=3)
20 configuration.make_lattice(gp.unit_cells.FCC, [8, 8, 8], rho=0.973)
21 configuration['m'] = 1.0
22 configuration.randomize_velocities(temperature=temperature * 2)
24 integrator = gp.integrators.NVT(temperature=temperature, tau=0.2, dt=0.005)
26 runtime_actions = [ gp.TrajectorySaver(),
27 gp.ScalarSaver(),
28 gp.MomentumReset(100)]
30 sim = gp.Simulation(configuration, pair_potential, integrator, runtime_actions,
31 steps_per_timeblock=1024, num_timeblocks=16,
32 storage='memory')
34 # if verbose:
35 # print('Equilibrating...')
36 # for _ in sim.run_timeblocks():
37 # if verbose:
38 # print(sim.status(per_particle=True))
39 #
40 # if verbose:
41 # print('Calculating structure factor in production run ...')
43 q_max: float = 16.0
44 calc_struct_fact = gp.CalculatorStructureFactor(configuration)
45 calc_struct_fact.generate_q_vectors(q_max=q_max)
46 for _ in sim.run_timeblocks():
47 calc_struct_fact.update()
48 # if verbose:
49 # print(sim.status(per_particle=True))
51 # Read the (binned) structure factor
52 # if verbose:
53 # print('Testing binned data output ...')
54 struc_fact = calc_struct_fact.read(bins=128)
56 # Assert that the structure factor is a dictionary
57 assert isinstance(struc_fact, dict)
58 # Assert that the dictionary contains the keys '|q|' and 'S(|q|)'
59 assert '|q|' in struc_fact, f"struc_fact.keys() = {struc_fact.keys()}"
60 assert 'S(|q|)' in struc_fact, f"struc_fact.keys() = {struc_fact.keys()}"
61 # Assert that the values of the dictionary are numpy arrays
62 assert isinstance(struc_fact['|q|'], type(np.array([])))
63 assert isinstance(struc_fact['S(|q|)'], type(np.array([])))
64 # Assert that the length of the arrays are the same
65 assert len(struc_fact['|q|']) == len(struc_fact['S(|q|)'])
66 # Assert that q values are between 0 and max_q
67 assert np.all(struc_fact['|q|'] >= 0)
68 assert np.all(struc_fact['|q|'] <= q_max)
69 # Assert that max S_q is between 2 and 3
70 assert np.max(struc_fact['S(|q|)']) < 3, f"max(S(q)) = {np.max(struc_fact['S(|q|)'])}"
71 assert np.max(struc_fact['S(|q|)']) > 2, f"max(S(q)) = {np.max(struc_fact['S(|q|)'])}"
72 # Assert that min S_q is between 0 and 1
73 assert np.min(struc_fact['S(|q|)']) < 1, f"min(S(q)) = {np.min(struc_fact['S(|q|)'])}"
74 assert np.min(struc_fact['S(|q|)']) >= 0, f"min(S(q)) = {np.min(struc_fact['S(|q|)'])}"
76 # if verbose:
77 # # Print shape of the arrays
78 # for key in struc_fact.keys():
79 # print(f"Shape of {key}: {struc_fact[key].shape}")
80 # print(f"Max of S(q): {np.max(struc_fact['S(|q|)'])}")
81 # print(f"Length of q at max of S(q): {struc_fact['|q|'][np.argmax(struc_fact['S(|q|)'])]}")
82 #
83 # if plot:
84 # plt.figure()
85 # plt.title('Structure factor')
86 # plt.plot(struc_fact['|q|'], struc_fact['S(|q|)'], 'o--')
87 # plt.xlabel('|q|')
88 # plt.ylabel('S(|q|)')
89 # plt.xlim(0, None)
90 # plt.ylim(0, None)
91 # plt.show()
92 #
93 # if verbose:
94 # print('Testing raw data output ...')
96 struc_fact_raw = calc_struct_fact.read(bins=None)
97 # Assert that the output is a dictionary
98 assert isinstance(struc_fact_raw, dict)
99 # Assert that the dictionary contains the keys 'q', '|q|', 'S(|q|)', 'rho_q', 'n_vectors'
100 assert 'q' in struc_fact_raw, f"struc_fact_raw.keys() = {struc_fact_raw.keys()}"
101 assert '|q|' in struc_fact_raw, f"struc_fact_raw.keys() = {struc_fact_raw.keys()}"
102 assert 'S(q)' in struc_fact_raw, f"struc_fact_raw.keys() = {struc_fact_raw.keys()}"
103 assert 'rho_q' in struc_fact_raw, f"struc_fact_raw.keys() = {struc_fact_raw.keys()}"
104 assert 'n_vectors' in struc_fact_raw, f"struc_fact_raw.keys() = {struc_fact_raw.keys()}"
105 # Assert that the values of the dictionary are numpy arrays
106 assert isinstance(struc_fact_raw['q'], type(np.array([])))
107 assert isinstance(struc_fact_raw['|q|'], type(np.array([])))
108 assert isinstance(struc_fact_raw['S(q)'], type(np.array([])))
109 assert isinstance(struc_fact_raw['rho_q'], type(np.array([])))
110 assert isinstance(struc_fact_raw['n_vectors'], type(np.array([])))
111 # Assert that the length of the arrays are the same
113 # if verbose:
114 # # Print dimensions of the arrays
115 # for key in struc_fact_raw.keys():
116 # print(f"Shape of {key}: {struc_fact_raw[key].shape}")
117 #
118 # if plot:
119 # # Plot image of 2D S(q) for q_z = 0
120 # q: np.ndarray = struc_fact_raw['q']
121 # S_q: np.ndarray = struc_fact_raw['S(q)']
122 # # Select q_vectors where q_z = 0
123 # mask = np.isclose(q[:, 2], 0)
124 # q = q[mask]
125 # q_x = q[:, 0]
126 # q_y = q[:, 1]
127 # S_q = S_q[mask]
128 # # Plot image of 2D S(q)
129 # plt.figure()
130 # plt.title('Structure factor')
131 # plt.scatter(q_x, q_y, c=S_q, s=50, alpha=0.8, cmap='viridis')
132 # plt.xlabel('q_x')
133 # plt.ylabel('q_y')
134 # plt.colorbar()
135 # plt.show()
136 #
137 # if verbose:
138 # print('All tests passed successfully')
140def test_structure_factor_backends():
141 # Test structure factor calculation with different backends
142 backends = ['CPU multi core', 'CPU single core']
143 sim = gp.get_default_sim()
144 for backend in backends:
145 print(f'Testing backend: {backend}')
147 # Test direct
148 n_vectors = [[5,6,7], [0, 0, 1], [2, -2, -5]]
149 calc_struct_fact = gp.CalculatorStructureFactor(sim.configuration, backend=backend, n_vectors=n_vectors)
150 calc_struct_fact.update()
151 struc_fact = calc_struct_fact.read(bins=8)
153 # Test generated q_vectors
154 calc_struct_fact = gp.CalculatorStructureFactor(sim.configuration, backend=backend)
155 calc_struct_fact.generate_q_vectors(q_max=10.0)
156 calc_struct_fact.update()
157 struc_fact = calc_struct_fact.read(bins=128)
159def test_atomic_form_factors():
160 D = 3
161 rho = 1.0
162 number_of_particles = 10_000
163 conf = gp.Configuration(D=D)
164 conf.make_positions(N=number_of_particles, rho=rho)
165 conf['m'] = 1.0
167 n_vectors = np.array([[5,6,7], [0, 0, 1], [2, -2, -5]])
168 atomic_form_factors = np.random.random(number_of_particles)
170 calc_single = gp.CalculatorStructureFactor(conf, n_vectors, atomic_form_factors, 'CPU single core')
171 number_of_updates = 4
172 for _ in range(number_of_updates):
173 conf['r'] = (np.random.rand(number_of_particles, D)-0.5) * conf.simbox.get_lengths() # Ideal gas configuration
174 conf.copy_to_device()
175 calc_single.update()
176 calc_multi = gp.CalculatorStructureFactor(conf, n_vectors, atomic_form_factors, 'CPU multi core')
177 for _ in range(number_of_updates):
178 conf['r'] = (np.random.rand(number_of_particles, D)-0.5) * conf.simbox.get_lengths() # Ideal gas configuration
179 conf.copy_to_device()
180 calc_multi.update()
182def test_structure_factor_gpu():
183 D = 3
184 rho = 1.0
185 number_of_particles = 10_000
186 conf = gp.Configuration(D=D)
187 conf.make_positions(N=number_of_particles, rho=rho)
188 conf['m'] = 1.0
190 # Ideal gas
191 conf['r'] = (np.random.rand(number_of_particles, D)-0.5) * conf.simbox.get_lengths()
193 # Test direct
194 calc_struct_fact = gp.CalculatorStructureFactor(conf, backend='GPU')
195 calc_struct_fact.generate_q_vectors(q_max=10.0)
196 print(f'{calc_struct_fact.q_vectors.shape = }')
197 calc_struct_fact.update()
198 struc_fact = calc_struct_fact.read(bins=32)
200 assert np.isclose(sum(struc_fact['S(|q|)'])/len(struc_fact['S(|q|)']), 1.0, atol=0.2)
203if __name__ == '__main__': # pragma: no cover
204 test_structure_factor()
205 test_structure_factor_backends()
206 test_atomic_form_factors()
207 test_structure_factor_gpu()