Coverage for gamdpy/configuration/Configuration.py: 80%
360 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import numpy as np
2import numba
3import math
4from numba import cuda
5from .colarray import colarray
6from ..simulation_boxes import Orthorhombic, LeesEdwards
7from .topology import Topology, duplicate_topology, replicate_topologies
8from ..simulation.get_default_compute_flags import get_default_compute_flags
10# IO
11import h5py
12import gzip
14# TODO: add possibility of "with ... as conf:" TypeError: 'Configuration' object does not support the context manager protocol
16class Configuration:
17 """ The configuration class
19 Store particle vectors (positions, velocities, forces) and scalars (energy, virial, mass ...).
20 Also store particle type, image coordinates, and the simulation box.
22 Parameters
23 ----------
24 D : int
25 Spatial dimension for the configuration.
27 N : int [Optional]
28 Number of particles.
29 If not set, this will be determined the first time particle data is written to the configuration.
31 Examples
32 -------
34 >>> import gamdpy as gp
35 >>> conf = gp.Configuration(D=3, N=1000)
36 >>> print(conf.vector_columns) # Print names of vector columns
37 ['r', 'v', 'f']
38 >>> print(conf.scalar_columns) # Print names of scalar columns
39 ['U', 'W', 'K', 'm']
40 >>> print(conf['r'].shape) # Vectors are stored as (N, D) numpy arrays
41 (1000, 3)
42 >>> print(conf['m'].shape) # Scalars are stored as (N,) numpy arrays
43 (1000,)
46 Data can be accessed via string keys (similar to dataframes in pandas):
48 >>> conf['r'] = np.ones((1000, 3))
49 >>> conf['v'] = 2 # Broadcast by numpy to correct shape
50 >>> print(conf['r'] + 0.01*conf['v'])
51 [[1.02 1.02 1.02]
52 [1.02 1.02 1.02]
53 [1.02 1.02 1.02]
54 ...
55 [1.02 1.02 1.02]
56 [1.02 1.02 1.02]
57 [1.02 1.02 1.02]]
60 A configuration can be specified without setting the number particles, N.
61 In that case N is determined the first time the particle data is written to the configuration:
63 >>> import numpy as np
64 >>> conf = gp.Configuration(D=3)
65 >>> conf['r'] = np.zeros((400, 3))
66 >>> print(conf['r'].shape)
67 (400, 3)
69 """
71 scalar_parameters = ['m']
72 scalar_computables_interactions = ['U', 'W', 'lapU']
73 scalar_computables_integrator = ['K', 'Fsq']
74 scalar_decriptions = {'m': 'Particle mass.',
75 'U': 'Potential energy.',
76 'W': 'Virial.',
77 'lapU': 'Laplace(U).',
78 'K': 'Kinetic energy.',
79 'Fsq': 'Squared length of force vector.',
80 }
83 def __init__(self, D: int, N: int = None, type_names=None, compute_flags=None, ftype=np.float32, itype=np.int32) -> None:
84 self.D = D
85 self.N = N
87 self.type_names = type_names
88 self.index_from_type_name = {}
89 if type_names:
90 for index, type_name in enumerate(type_names):
91 self.index_from_type_name[type_name] = index
93 self.compute_flags = get_default_compute_flags()
94 if compute_flags != None:
95 # only keys present in the default are processed
96 for k in compute_flags:
97 if k in self.compute_flags:
98 self.compute_flags[k] = compute_flags[k]
99 else:
100 raise ValueError('Unknown key in compute_flags:%s' %k)
102 self.vector_columns = ['r', 'v', 'f'] # Should be user modifiable
103 if self.compute_flags['stresses']:
104 if self.D > 4:
105 raise ValueError("compute_flags['stresses'] should not be set for D>4")
106 self.vector_columns += ['sx', 'sy', 'sz','sw'][:self.D]
109 self.num_cscalars = 0
110 self.sid = {}
111 self.scalar_columns = []
112 sid_index = 0
114 for label in self.scalar_computables_interactions:
115 if self.compute_flags[label]:
116 self.sid[label] = sid_index
117 self.scalar_columns.append(label)
118 sid_index += 1
119 self.num_cscalars += 1
121 for label in self.scalar_computables_integrator:
122 if self.compute_flags[label]:
123 self.sid[label] = sid_index
124 self.scalar_columns.append(label)
125 sid_index += 1
128 for label in self.scalar_parameters:
129 self.sid[label] = sid_index
130 self.scalar_columns.append(label)
131 sid_index += 1
133 self.simbox = None
134 self.topology = Topology()
135 self.ptype_function = self.make_ptype_function()
136 self.ftype = ftype
137 self.itype = itype
138 if self.N != None:
139 self.__allocate_arrays()
141 def __allocate_arrays(self):
142 self.vectors = colarray(self.vector_columns, size=(self.N, self.D), dtype=self.ftype)
143 self.scalars = np.zeros((self.N, len(self.scalar_columns)), dtype=self.ftype)
144 self.r_im = np.zeros((self.N, self.D), dtype=self.itype) # Move to vectors
145 self.ptype = np.zeros(self.N, dtype=self.itype) # Move to scalars
146 return
148 def __repr__(self):
149 return f'Configuration(D={self.D}, N={self.N}, compute_flags={self.compute_flags})'
151 def __code__(self):
152 import sys
153 np.set_printoptions(threshold=sys.maxsize)
154 code_str = "# Define configuration class\n"
155 code_str += f"from gamdpy import Configuration\n"
156 code_str += f"configuration = Configuration(D={self.D}, N={self.N}, compute_flags={self.compute_flags})\n"
157 # Following part needs to be done with a read function from the .h5
158 for key in self.vector_columns:
159 code_str += f"configuration['{key}'] = {self[key]}\n"
160 for key in self.scalar_columns:
161 code_str += f"configuration['{key}'] = {self[key]}\n"
162 return code_str
164 def __str__(self):
165 if self.N == None:
166 return f'{self.D} dimensional configuration. Particles not yet assigned.'
167 str = f'{self.N} particles in {self.D} dimensions. Number density (atomic): {self.N/self.get_volume():.3f}'
168 num_types = np.max(self.ptype)+1
169 if num_types==1:
170 str += '. Single component. '
171 else:
172 str += f'. {num_types} components with fractions '
173 for ptype in range(num_types):
174 str += f'{np.mean(self.ptype==ptype):.3f}, '
175 str += '\nCurrent scalar data per particle:'
176 for key in self.sid:
177 str += f'\n{key+",":5} {np.mean(self.scalars[:,self.sid[key]]):.3f}'
178 if key in self.scalar_decriptions:
179 str += '\t' + self.scalar_decriptions[key]
180 return str
182 def __setitem__(self, key, data):
183 if self.N is None: # First time setting particle data, so allocate arrays
184 if type(data) != np.ndarray:
185 raise (TypeError)(
186 'Number of particles, N, not determined yet, so assignment needs to be with a numpy array')
187 self.N = data.shape[0]
188 self.__allocate_arrays()
190 if key in self.vector_columns:
191 self.__set_vector(key, data)
192 return
193 if key in self.scalar_columns:
194 self.__set_scalar(key, data)
195 return
196 raise ValueError(f'Unknown key {key}. Vectors: {self.vector_columns}, Scalars: {self.scalar_columns}')
198 def __set_vector(self, key: str, data: np.ndarray) -> None:
199 """ Set new vector data """
201 if type(data) == np.ndarray: # Allow for possibility of using scalar float, which is then broadcast by numpy
202 N, D = data.shape
203 if N != self.N:
204 raise ValueError(f'Inconsistent number of particles, {N} <> {self.N}')
205 if D != self.D:
206 raise ValueError(f'Inconsistent number of dimensions, {D} <> {self.D}')
207 self.vectors[key] = data
208 return
210 def __set_scalar(self, key: str, data) -> None:
211 """ Set new scalar data """
213 if type(data) == np.ndarray: # Allow for possibility of using scalar float, which is then broadcast by numpy
214 N, = data.shape
215 if N != self.N:
216 raise ValueError(f'Inconsistent number of particles, {N} <> {self.N}')
217 self.scalars[:, self.sid[key]] = data
218 return
220 def __getitem__(self, key):
221 if key in self.vector_columns:
222 return self.vectors[key]
223 if key in self.scalar_columns:
224 return self.scalars[:, self.sid[key]]
225 raise ValueError(f'Unknown key {key}. Vectors: {self.vector_columns}, Scalars: {self.scalar_columns}')
227 def copy_to_device(self):
228 """ Copy all data to device memory """
229 self.d_vectors = cuda.to_device(self.vectors.array)
230 self.d_scalars = cuda.to_device(self.scalars)
231 self.d_r_im = cuda.to_device(self.r_im)
232 self.d_ptype = cuda.to_device(self.ptype)
233 self.simbox.copy_to_device()
234 return
236 def copy_to_host(self):
237 """ Copy all data to host memory """
238 self.vectors.array = self.d_vectors.copy_to_host()
239 self.scalars = self.d_scalars.copy_to_host()
240 self.r_im = self.d_r_im.copy_to_host()
241 self.ptype = self.d_ptype.copy_to_host()
242 self.simbox.copy_to_host()
243 return
245 def make_ptype_function(self) -> callable:
246 def ptype_function(pid, ptype_array):
247 ptype = ptype_array[pid] # Default: read from ptype_array
248 return ptype
250 return ptype_function
252 def get_potential_energy(self) -> float:
253 """ Get total potential energy of the configuration """
254 return float(np.sum(self['U']))
256 def get_volume(self):
257 """ Get volume of simulation box associated with configuration """
258 return self.simbox.get_volume()
260 def set_kinetic_temperature(self, temperature, ndofs=None):
261 if ndofs is None:
262 ndofs = self.D * (self.N - 1)
264 T_ = np.sum(np.dot(self['m'], np.sum(self['v'] ** 2, axis=1))) / ndofs
265 if T_ == 0:
266 raise ValueError('Cannot rescale velocities when all equal to zero')
267 self['v'] *= (temperature / T_) ** 0.5
269 def randomize_velocities(self, temperature, seed=None, ndofs=None):
270 """ Randomize velocities according to a given temperature. If T <= 0, set all velocities to zero. """
271 if self.D is None:
272 raise ValueError('Cannot randomize velocities. Start by assigning positions.')
273 masses = self['m']
274 if np.any(masses == 0):
275 raise ValueError('Cannot randomize velocities when any mass is zero')
276 if temperature > 0.0:
277 self['v'] = generate_random_velocities(self.N, self.D, T=temperature, seed=seed, m=self['m'])
278 # rescale to get the kinetic temperature exactly right
279 self.set_kinetic_temperature(temperature=temperature, ndofs=ndofs)
280 else:
281 self['v'] = np.zeros((self.N, self.D), np.float32)
283 def make_lattice(self, unit_cell: dict, cells: list, rho: float = None) -> None:
284 """ Generate a lattice configuration
286 The lattice is constructed by replicating the unit cell in all directions.
287 Unit cell is a dictonary with `fractional_coordinates` for particles, and
288 the `lattice_constants` as a list of unit cell lengths in all directions.
290 Unit cells are available in gamdpy.unit_cells
292 Example
293 -------
295 >>> import gamdpy as gp
296 >>> conf = gp.Configuration(D=3)
297 >>> conf.make_lattice(gp.unit_cells.FCC, cells=[8, 8, 8], rho=1.0)
298 >>> print(gp.unit_cells.FCC) # Example of a unit cell dict
299 {'fractional_coordinates': [[0.0, 0.0, 0.0], [0.5, 0.5, 0.0], [0.5, 0.0, 0.5], [0.0, 0.5, 0.5]], 'lattice_constants': [1.0, 1.0, 1.0]}
301 """
302 from .make_lattice import make_lattice
303 positions, box_vector = make_lattice(unit_cell=unit_cell, cells=cells, rho=rho)
304 self['r'] = positions
305 self.simbox = Orthorhombic(self.D, box_vector)
306 return
308 def make_positions(self, N, rho):
309 """
310 Generate particle positions in D dimensions.
312 Positions are generated in a simple cubic configuration in D dimensions.
313 Takes the number of particles N and the density rho as inputs
315 Example
316 -------
318 >>> import gamdpy as gp
319 >>> conf = gp.Configuration(D=3)
320 >>> conf.make_positions(N=27, rho=0.2)
321 """
323 D = self.D
324 part_per_line = np.ceil(pow(N, 1./D))
326 box_length = pow(N/rho, 1./D)
327 box_vector = np.array(D*[box_length])
329 index = 0
330 x = [] # empty list
332 # This loop places particles in a simple cubic configuration
333 # The first particle is in D*[0]
334 while index < N:
335 dcurrent = D - 1
336 i_d = D*[float(0)]
337 i_d[dcurrent] = (index / pow(part_per_line, dcurrent))
338 rest = index % (pow(part_per_line, dcurrent))
339 while dcurrent != 0:
340 dcurrent = dcurrent - 1
341 i_d[dcurrent] = (rest/pow(part_per_line,dcurrent))
342 rest = index % (pow(part_per_line, dcurrent))
343 x.append(i_d)
344 index = index + 1
345 pos = np.array(x)
346 # Centering the array
347 dcurrent = 1
348 remove = 0
349 while dcurrent < D:
350 remove += D**(D-dcurrent)
351 pos[:, dcurrent] -= remove/N
352 dcurrent = dcurrent + 1
353 pos -= np.array(D*[int(0.5*part_per_line)]) # center cube at 0
354 # Scaling for density
355 pos *= box_length/part_per_line
356 # Saving to Configuration object
357 self['r'] = pos
358 self.simbox = Orthorhombic(self.D, box_vector)
359 # Check all particles are in the box (-L/2, L/2)
360 assert np.any(np.abs(pos))<0.5*box_length
362 return
364 def atomic_scale(self, density):
365 actual_rho = self.N / self.get_volume()
366 scale_factor = (actual_rho / density)**(1/3)
367 self.vectors['r'] *= scale_factor
368 self.simbox.scale(scale_factor)
371 def save(self, output: h5py.File, group_name: str, mode: str="w", include_topology: bool=False) -> None:
372 """ Write a configuration to a HDF5 file
374 Parameters
375 ----------
377 configuration : gamdpy.Configuration
378 a gamdpy configuration object
380 output : h5py.File
381 h5 file
383 group_name : str
384 name of the group which will be created in the h5 and in which
385 the configuration will be saved
387 mode: str
388 default value is "w" and corresponds to replacing existing dataset
390 Example
391 -------
393 >>> import os
394 >>> import h5py
395 >>> import gamdpy as gp
396 >>> conf = gp.Configuration(D=3)
397 >>> conf.make_positions(N=10, rho=1.0)
398 >>> conf.save(output=h5py.File("final.h5", "w"), group_name="configuration", mode="w")
399 >>> os.remove("final.h5") # Removes file (for doctests)
400 >>> with h5py.File("manyconfs.h5", "a") as fout:
401 ... conf.save(output=fout, group_name="restarts/restart0000", mode="w")
402 ... conf.save(output=fout, group_name="restarts/restart0001", mode="w")
403 ... conf.save(output=fout, group_name="restarts/restart0002", mode="w")
404 >>> os.remove("manyconfs.h5") # Removes file (for doctests)
406 """
408 # Sanity:
409 #print(f"output {isinstance(output, h5py.File)} {output}")
410 #print(f"group_name {isinstance(group_name, str)} {group_name}")
411 # Creating group group_name in h5 root
412 if group_name in output.keys() and mode=="w":
413 print(f"{group_name} already present in h5 root, replacing it")
414 del output[f'{group_name}']
415 output.create_group(group_name)
416 # Checks if group group_name exists in case mode="append"
417 elif group_name not in output.keys() and mode=="a":
418 output.create_group(group_name)
419 elif group_name not in output.keys() and mode=="w":
420 output.create_group(group_name)
421 elif group_name in output.keys() and mode=="a":
422 print(f"append data to {group_name} in h5 root")
423 else:
424 raise ValueError("Unexpected combination of input in save method of Configuration")
426 # Save attributes of group group_name
427# output[group_name].attrs['simbox'] = self.simbox.get_lengths()
429 # Saving vectors separately
430 #output[group_name].create_dataset('r', data=self['r'], dtype=np.float32)
431 #output[f"{group_name}/r"].attrs['simbox'] = self.simbox.get_lengths()
432 #output[group_name].create_dataset('v', data=self['v'], dtype=np.float32)
433 #output[group_name].create_dataset('f', data=self['f'], dtype=np.float32)
435 # Saving vectors all together
436 #output[group_name].create_dataset('vectors', data=np.hstack([self['r'], self['v'], self['f']]), dtype=np.float32)
437 output[group_name].create_dataset('vectors', data=self.vectors.array, dtype=np.float32)
438 output[f"{group_name}/vectors"].attrs['vector_columns'] = self.vector_columns
440 # Saving other things
441 output[group_name].create_dataset('ptype', data=self.ptype, dtype=np.int32)
442 #output[group_name].create_dataset('m', data=self['m'], dtype=np.float32) # included in scalars
443 output[group_name].create_dataset('r_im', data=self.r_im, dtype=np.int32)
444 output[group_name].create_dataset('scalars', data=self.scalars, dtype=np.float32)
445 output[f"{group_name}/scalars"].attrs['scalar_columns'] = self.scalar_columns
447 # save simulation box
448 output[group_name].attrs['simbox_name'] = self.simbox.get_name()
449 #output[group_name].attrs['simbox_data'] = self.simbox.get_lengths()
450 output[group_name].attrs['simbox_data'] = self.simbox.data_array
452 # save topology, depending on flag
453 if include_topology:
454 output[group_name].create_group('topology')
455 self.topology.save(output[f'{group_name}/topology'])
457 # The following is equivalent to overloading in c++ : https://stackoverflow.com/questions/12179271/meaning-of-classmethod-and-staticmethod-for-beginner
458 # cls stands for class, in this case the Configuration class
459 @classmethod
460 def from_h5(cls, h5file: h5py.File, group_name: str, reset_images: bool=False, compute_flags: bool=None):
461 """ Read a configuration from an open HDF5 file identified by group-name
463 Parameters
464 ----------
465 h5file : HDF5 File
466 open HDF5 open, as returned by h5py.File()
468 group_name : str
469 string corresponding to a key in the h5py.File containing a saved gamdpy configuration
471 reset_images : bool
472 if True set the images to zero (default False)
474 compute_flags : bool
475 NOTE: still to be developed, should be possible to define compute flags from dictionary
476 compute_flags defining what will be stored in the configuration (default None)
478 Returns
479 -------
481 configuration : gamdpy.Configuration
482 a gamdpy configuration object
485 Example:
486 --------
488 >>> import gamdpy as gp
489 >>> output_file = h5py.File('examples/Data/LJ_r0.973_T0.70_toread.h5')
490 >>> conf = Configuration.from_h5(output_file, 'restarts/restart0000')
491 >>> print(conf.D, conf.N, conf['r'][0]) # Print number of dimensions D, number of particles N and position of first particle
492 3 2048 [-6.384221 -6.3622074 -6.3125153]
494 """
497 vectors_array = h5file[group_name]['vectors'][:]
498 _, N, D = vectors_array.shape
499 configuration = cls(D=D, N=N, compute_flags=compute_flags)
501 configuration.vector_columns = h5file[group_name]['vectors'].attrs['vector_columns']
502 configuration.scalar_columns = h5file[group_name]['scalars'].attrs['scalar_columns']
504 configuration.ptype = h5file[group_name]['ptype'][:]
506 scalars_array = h5file[group_name]['scalars'][:]
507 configuration.scalars = scalars_array
508 configuration.vectors.array = vectors_array
510 simbox_name = h5file[group_name].attrs['simbox_name']
511 simbox_data = h5file[group_name].attrs['simbox_data']
513 if simbox_name == 'Orthorhombic':
514 configuration.simbox = Orthorhombic(D, simbox_data)
515 elif simbox_name == 'LeesEdwards':
516 box_shift_image = {True:0, False: int(simbox_data[D+1])} [reset_images]
517 configuration.simbox = LeesEdwards(D, simbox_data[:D], simbox_data[D], box_shift_image)
518 else:
519 raise ValueError('simbox_name %s not recognized in group %s' % (simbox_name, group_name))
521 if reset_images:
522 configuration.r_im = np.zeros((N, D), dtype=np.int32)
523 else:
524 configuration.r_im = h5file[group_name]['r_im'][:]
525 return configuration
528# Helper functions
530def generate_random_velocities(N, D, T, seed, m=1, dtype=np.float32):
531 """ Generate random velocities according to a given temperature. """
532 v = np.zeros((N, D), dtype=dtype)
533 # default value of seed is None and random.seed(None) has no effect
534 np.random.seed(seed)
535 for k in range(D):
536 # to cover the case that m is a 1D array of length N, need to
537 # generate one column at a time, passing the initial zeros as the
538 # mean to avoid problems with inferring the correct shape
539 v[:, k] = np.random.normal(v[:, k], (T / m) ** 0.5)
540 CM_drift = np.mean(m * v[:, k]) / np.mean(m)
541 v[:, k] -= CM_drift
542 return dtype(v)
544@numba.njit
545def generate_fcc_positions(nx, ny, nz, rho, dtype=np.float32):
546 # This function is not recommended to use, and should be considered deprecated
547 # raise DeprecationWarning('Use Configuration.make_lattice() instead')
549 D = 3
550 conf = np.zeros((nx * ny * nz * 4, D), dtype=dtype)
551 count = 0
552 for ix in range(nx):
553 for iy in range(ny):
554 for iz in range(nz):
555 conf[count + 0, :] = [ix + 0.25, iy + 0.25, iz + 0.25]
556 conf[count + 1, :] = [ix + 0.75, iy + 0.75, iz + 0.25]
557 conf[count + 2, :] = [ix + 0.75, iy + 0.25, iz + 0.75]
558 conf[count + 3, :] = [ix + 0.25, iy + 0.75, iz + 0.75]
559 count += 4
560 for k in range(D):
561 conf[:, k] -= np.mean(conf[:, k]) # put sample in the middle of the box
562 sim_box = np.array((nx, ny, nz), dtype=dtype)
563 rho_initial = 4.0
564 scale_factor = dtype((rho_initial / rho) ** (1 / D))
566 return conf * scale_factor, sim_box * scale_factor
569def make_configuration_fcc(nx, ny, nz, rho, N=None):
570 """
571 Generate Configuration for particle positions and simbox of a FCC lattice with a given density
572 (nx x ny x nz unit cells),
573 and default types ('0') and masses ('1.')
574 If N is given, only N particles will be in the configuration
575 (needs to be equal to or smaller than number of particle in generated crystal)
576 """
578 # This function is not recommended to use, and should be considered deprecated
579 # raise DeprecationWarning('Use Configuration.make_lattice() instead')
581 positions, simbox_data = generate_fcc_positions(nx, ny, nz, rho)
582 N_, D = positions.shape
583 if N == None:
584 N = N_
585 else:
586 if N > N_:
587 raise ValueError(
588 f'N ({N}) needs to be equal to or smaller than number of particle in generated crystal ({N_})')
589 scale_factor = (N / N_) ** (1 / 3)
590 positions *= scale_factor
591 simbox_data *= scale_factor
593 configuration = Configuration(D=3)
594 configuration['r'] = positions[:N, :]
595 configuration.simbox = Orthorhombic(D, simbox_data)
596 configuration['m'] = np.ones(N, dtype=np.float32) # Set masses
597 configuration.ptype = np.zeros(N, dtype=np.int32) # Set types
599 return configuration
602def replicate_molecules(molecule_dicts, num_molecules_each_type_list, safety_distance, random_rotations=True, compute_flags=None):
603 """ Construct a configuration containing different molecules, with the numbers of each type specified
605 Parameters:
606 moleculde_dicts (list): A list of dictionaries, each of which contains keys "positions", "particle_types", "masses" and "topology", whose values are corresponding lists of data for that molecule
607 num_molecules_each_type_list (list): A list of integers, specifying how many molecules of each type are to be included
608 safety_distance (float): A length to be added in all directions to the size of the bounding box to be used for each molecule when placing them initially on a lattice
609 random_rotation (Bool): Whether the x,y,z coordinates in each molecule should be randomly permutated to give a simple randomization of orientations.
610 Returns:
611 configuration (Configuration): the resulting configuration with all molecules replicated and including the corresponding replicated topology
612 """
613 D = 3
614 num_molecule_types = len(molecule_dicts)
615 total_num_particles = 0
616 total_num_molecules = 0
617 mol_types = []
618 positions_array_list = []
619 cell_length_list = []
620 size_molecule_type_list = []
622 # unpack the list of molecule dictionaries and make lists of positions, particle_types, masses and topologies
623 mol_topology_list = []
624 mol_positions_list = []
625 particle_type_list = []
626 masses_list = []
627 for idx in range(num_molecule_types):
628 mol_positions_list.append(molecule_dicts[idx]["positions"])
629 particle_type_list.append(molecule_dicts[idx]["particle_types"])
630 masses_list.append(molecule_dicts[idx]["masses"])
631 mol_topology_list.append(molecule_dicts[idx]["topology"])
633 # tally the total numbers of particles and molecules, make shifted arrays of possitions for each molecule type
634 for idx in range(num_molecule_types):
635 num_mol_this_type = num_molecules_each_type_list[idx]
636 total_num_particles += len(mol_positions_list[idx]) * num_mol_this_type
637 total_num_molecules += num_mol_this_type
638 mol_types += [idx] * num_mol_this_type
639 positions_array = np.array(mol_positions_list[idx])
640 positions_array -= np.min(positions_array, axis=0)
641 positions_array_list.append(positions_array)
642 size_molecule_type_list.append( len(mol_positions_list[idx]) )
643 cell_length = np.max(positions_array) + safety_distance
644 cell_length_list.append(cell_length)
646 # shuffle molecule types randomly
647 np.random.shuffle(mol_types)
650 configuration = Configuration(D=D, N=total_num_particles, compute_flags=compute_flags)
651 configuration.topology = replicate_topologies(mol_topology_list, num_molecules_each_type_list, mol_types, size_molecule_type_list)
653 max_cell_length = max(cell_length_list)
654 # make a cubic box big enough to hold the total number of molecules
655 num_cells_axis = math.ceil(total_num_molecules**(1/3))
656 simbox_data = np.ones(D) * (num_cells_axis * max_cell_length)
657 configuration.simbox = Orthorhombic(D, simbox_data)
659 mol_count = 0
660 particle_count = 0
661 for ix in range(num_cells_axis):
662 for iy in range(num_cells_axis):
663 for iz in range(num_cells_axis):
664 if mol_count < total_num_molecules:
665 # add a molecule
666 this_mol_type = mol_types[mol_count]
667 particles_this_molecule = size_molecule_type_list[this_mol_type]
668 arr = np.arange(D)
669 if random_rotations:
670 np.random.shuffle(arr)
672 configuration['r'][particle_count:particle_count+particles_this_molecule,0] = positions_array_list[this_mol_type][:,arr[0]] + ix*max_cell_length
673 configuration['r'][particle_count:particle_count+particles_this_molecule,1] = positions_array_list[this_mol_type][:,arr[1]] + iy*max_cell_length
674 configuration['r'][particle_count:particle_count+particles_this_molecule,2] = positions_array_list[this_mol_type][:,arr[2]] + iz*max_cell_length
675 configuration['m'][particle_count:particle_count+particles_this_molecule] = masses_list[this_mol_type]
676 configuration.ptype[particle_count:particle_count+particles_this_molecule] = particle_type_list[this_mol_type]
677 particle_count += particles_this_molecule
679 mol_count += 1
681 assert mol_count == total_num_molecules
682 assert particle_count == total_num_particles
684 for i in range(configuration.D):
685 configuration['r'][:,i] -= configuration.simbox.get_lengths()[i]/2
687 return configuration