diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 5b7b502a..5f30e616 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -8,11 +8,23 @@ Change Log `_. 5.x -___ +--- Next release ^^^^^^^^^^^^^^^^^^ +*Added:* + +* ``gsd.hoomd`` schema now supports ``float64`` and ``float32`` values. + Choose the desired precision with the ``precision`` argument to ``open`` (#495). + +*Changed:* + +* ``Frame.validate`` no longer modifies its contents in place by default. + Pass ``inplace=True`` to recover the previous behavior (#495). +* New files created by ``gsd`` 5.0 are schema version 2.0 and may not + be readable by older software (#495). + *Removed*: * ``SIGTERM`` handler (#501) diff --git a/doc/credits.rst b/doc/credits.rst index 985cca37..49fcbedd 100644 --- a/doc/credits.rst +++ b/doc/credits.rst @@ -20,3 +20,4 @@ The following people contributed to GSD. * Charlotte Shiqi Zhao, University of Michigan * Tim Moore, University of Michigan * Joseph Burkhart, University of Michigan +* Nicholas Craven, Vanderbilt University diff --git a/doc/schema-hoomd.rst b/doc/schema-hoomd.rst index 3f2845a6..53d003d9 100644 --- a/doc/schema-hoomd.rst +++ b/doc/schema-hoomd.rst @@ -14,13 +14,21 @@ chunks. Any newer reader will initialize new data chunks with default values when they are not present in an older version file. :Schema name: ``hoomd`` -:Schema version: 1.4 +:Schema version: 2.0 .. seealso:: `hoomd.State` for a full description of how HOOMD interprets this data. +Changes +------- + +Version 2.0 +^^^^^^^^^^^ + +* Each ``float`` data chunk in 2.0 schema files may be either *float32* or *float64*. + Use-cases --------- @@ -138,7 +146,7 @@ Configuration .. chunk:: configuration/box - :Type: float + :Type: float32 *or* float64 :Size: 6x1 :Default: [1,1,1,0,0,0] :Units: *varies* @@ -150,6 +158,8 @@ Configuration * ``box[0:3]``: :math:`(l_x, l_y, l_z)` the box length in each direction, in length units * ``box[3:]``: :math:`(xy, xz, yz)` the tilt factors, dimensionless values + .. versionchanged:: 2.0 + Particle data ------------- @@ -208,31 +218,37 @@ Attributes .. chunk:: particles/mass - :Type: float (32-bit) + :Type: float32 *or* float64 :Size: Nx1 :Default: 1.0 :Units: mass Store the mass of each particle. + .. versionchanged:: 2.0 + .. chunk:: particles/charge - :Type: float (32-bit) + :Type: float32 *or* float64 :Size: Nx1 :Default: 0.0 :Units: charge Store the charge of each particle. + .. versionchanged:: 2.0 + .. chunk:: particles/diameter - :Type: float (32-bit) + :Type: float32 *or* float64 :Size: Nx1 :Default: 1.0 :Units: length Store the diameter of each particle. + .. versionchanged:: 2.0 + .. chunk:: particles/body :Type: int32 @@ -246,7 +262,7 @@ Attributes .. chunk:: particles/moment_inertia - :Type: float (32-bit) + :Type: float32 *or* float64 :Size: Nx3 :Default: 0,0,0 :Units: mass * length^2 @@ -255,12 +271,14 @@ Attributes This inertia tensor is diagonal in the body frame of the particle. The default value is for point particles. + .. versionchanged:: 2.0 + Properties ^^^^^^^^^^ .. chunk:: particles/position - :Type: float (32-bit) + :Type: float32 *or* float64 :Size: Nx3 :Default: 0,0,0 :Units: length @@ -281,9 +299,11 @@ Properties Where :math:`l_x`, :math:`l_y`, :math:`l_z`, :math:`xy`, :math:`xz`, and :math:`yz` are the simulation box parameters (:chunk:`configuration/box`). + .. versionchanged:: 2.0 + .. chunk:: particles/orientation - :Type: float (32-bit) + :Type: float32 *or* float64 :Size: Nx4 :Default: 1,0,0,0 :Units: unit quaternion @@ -293,21 +313,25 @@ Properties where the quaternion is :math:`q = r + a_xi + a_yj + a_zk`. A unit quaternion has the property: :math:`\sqrt{r^2 + a_x^2 + a_y^2 + a_z^2} = 1`. + .. versionchanged:: 2.0 + Momenta ^^^^^^^^ .. chunk:: particles/velocity - :Type: float (32-bit) + :Type: float32 *or* float64 :Size: Nx3 :Default: 0,0,0 :Units: length/time Store the velocity of each particle :math:`(v_x, v_y, v_z)`. + .. versionchanged:: 2.0 + .. chunk:: particles/angmom - :Type: float (32-bit) + :Type: float32 *or* float64 :Size: Nx4 :Default: 0,0,0,0 :Units: quaternion @@ -315,6 +339,8 @@ Momenta Store the angular momentum of each particle as a quaternion. See the HOOMD documentation for information on how to convert to a vector representation. + .. versionchanged:: 2.0 + .. chunk:: particles/image :Type: int32 @@ -509,7 +535,7 @@ Topology .. chunk:: constraints/value - :Type: float + :Type: float32 *or* float64 :Size: Nx1 :Default: 0 :Units: length @@ -517,6 +543,8 @@ Topology Store the distance of each constraint. Each constraint defines a fixed distance between two particles. + .. versionchanged:: 2.0 + .. chunk:: constraints/group :Type: uint32 diff --git a/gsd/hoomd.py b/gsd/hoomd.py index 073120ab..7022f130 100644 --- a/gsd/hoomd.py +++ b/gsd/hoomd.py @@ -89,7 +89,24 @@ def box(self, box): if self.dimensions is None: self.dimensions = 2 if Lz == 0 else 3 - def validate(self): + def cast(self): + """Lazy cast attributes to writeable dictionary. + + Convert every array attribute to a `numpy.ndarray` of the proper + type and check that all attributes have the correct dimensions. + + Ignore any attributes that are ``None``. + + Returns: + castedDict (dict): Stored attributes by attribute string name. + """ + castedDict = {} + if self.box is not None: + castedDict['box'] = numpy.ascontiguousarray(self.box, dtype=numpy.float32) + castedDict['box'] = castedDict['box'].reshape([6]) + return castedDict + + def validate(self, inplace=False): """Validate all attributes. Convert every array attribute to a `numpy.ndarray` of the proper @@ -97,13 +114,18 @@ def validate(self): Ignore any attributes that are ``None``. + Args: + inplace (bool): Whether or not the arrays are modified in place. + Warning: Array attributes that are not contiguous numpy arrays will be - replaced with contiguous numpy arrays of the appropriate type. + replaced with contiguous numpy arrays of the appropriate type + if inplace is ``True``. """ - if self.box is not None: - self.box = numpy.ascontiguousarray(self.box, dtype=numpy.float32) - self.box = self.box.reshape([6]) + if inplace: + dataDict = self.cast() + if self.box is not None: + self.box = dataDict['box'] class ParticleData: @@ -194,55 +216,126 @@ def __init__(self): self.image = None self.type_shapes = None - def validate(self): - """Validate all attributes. + def _validate_precision(self, data): + """Maintain floats in numpy arrays.""" + outFloat = numpy.float32 + if isinstance(data, list): + outFloat = numpy.float64 + elif data.dtype == numpy.float64: + outFloat = numpy.float64 + return outFloat + + def cast(self): + """Lazy cast attributes to writeable dictionary. Convert every array attribute to a `numpy.ndarray` of the proper type and check that all attributes have the correct dimensions. Ignore any attributes that are ``None``. - Warning: - Array attributes that are not contiguous numpy arrays will be - replaced with contiguous numpy arrays of the appropriate type. + Returns: + castedDict (dict): Stored attributes by attribute string name. """ + castedDict = {} if self.position is not None: - self.position = numpy.ascontiguousarray(self.position, dtype=numpy.float32) - self.position = self.position.reshape([self.N, 3]) + castedDict['position'] = numpy.ascontiguousarray( + self.position, dtype=self._validate_precision(self.position) + ) + castedDict['position'] = castedDict['position'].reshape([self.N, 3]) if self.orientation is not None: - self.orientation = numpy.ascontiguousarray( - self.orientation, dtype=numpy.float32 + castedDict['orientation'] = numpy.ascontiguousarray( + self.orientation, dtype=self._validate_precision(self.orientation) ) - self.orientation = self.orientation.reshape([self.N, 4]) + castedDict['orientation'] = castedDict['orientation'].reshape([self.N, 4]) if self.typeid is not None: - self.typeid = numpy.ascontiguousarray(self.typeid, dtype=numpy.uint32) - self.typeid = self.typeid.reshape([self.N]) + castedDict['typeid'] = numpy.ascontiguousarray( + self.typeid, dtype=numpy.uint32 + ) + castedDict['typeid'].reshape([self.N]) if self.mass is not None: - self.mass = numpy.ascontiguousarray(self.mass, dtype=numpy.float32) - self.mass = self.mass.reshape([self.N]) + castedDict['mass'] = numpy.ascontiguousarray( + self.mass, dtype=self._validate_precision(self.mass) + ) + castedDict['mass'] = castedDict['mass'].reshape([self.N]) if self.charge is not None: - self.charge = numpy.ascontiguousarray(self.charge, dtype=numpy.float32) - self.charge = self.charge.reshape([self.N]) + castedDict['charge'] = numpy.ascontiguousarray( + self.charge, dtype=self._validate_precision(self.charge) + ) + castedDict['charge'] = castedDict['charge'].reshape([self.N]) if self.diameter is not None: - self.diameter = numpy.ascontiguousarray(self.diameter, dtype=numpy.float32) - self.diameter = self.diameter.reshape([self.N]) + castedDict['diameter'] = numpy.ascontiguousarray( + self.diameter, dtype=self._validate_precision(self.diameter) + ) + castedDict['diameter'] = castedDict['diameter'].reshape([self.N]) if self.body is not None: - self.body = numpy.ascontiguousarray(self.body, dtype=numpy.int32) - self.body = self.body.reshape([self.N]) + castedDict['body'] = numpy.ascontiguousarray(self.body, dtype=numpy.int32) + castedDict['body'] = castedDict['body'].reshape([self.N]) if self.moment_inertia is not None: - self.moment_inertia = numpy.ascontiguousarray( - self.moment_inertia, dtype=numpy.float32 + castedDict['moment_inertia'] = numpy.ascontiguousarray( + self.moment_inertia, dtype=self._validate_precision(self.moment_inertia) + ) + castedDict['moment_inertia'] = castedDict['moment_inertia'].reshape( + [self.N, 3] ) - self.moment_inertia = self.moment_inertia.reshape([self.N, 3]) if self.velocity is not None: - self.velocity = numpy.ascontiguousarray(self.velocity, dtype=numpy.float32) - self.velocity = self.velocity.reshape([self.N, 3]) + castedDict['velocity'] = numpy.ascontiguousarray( + self.velocity, dtype=self._validate_precision(self.velocity) + ) + castedDict['velocity'] = castedDict['velocity'].reshape([self.N, 3]) if self.angmom is not None: - self.angmom = numpy.ascontiguousarray(self.angmom, dtype=numpy.float32) - self.angmom = self.angmom.reshape([self.N, 4]) + castedDict['angmom'] = numpy.ascontiguousarray( + self.angmom, dtype=self._validate_precision(self.angmom) + ) + castedDict['angmom'] = castedDict['angmom'].reshape([self.N, 4]) if self.image is not None: - self.image = numpy.ascontiguousarray(self.image, dtype=numpy.int32) - self.image = self.image.reshape([self.N, 3]) + castedDict['image'] = numpy.ascontiguousarray(self.image, dtype=numpy.int32) + castedDict['image'] = castedDict['image'].reshape([self.N, 3]) + + if self.types is not None and (not len(set(self.types)) == len(self.types)): + msg = 'Type names must be unique.' + raise ValueError(msg) + return castedDict + + def validate(self, inplace=False): + """Validate all attributes. + + Convert every array attribute to a `numpy.ndarray` of the proper + type and check that all attributes have the correct dimensions. + + Ignore any attributes that are ``None``. + + Args: + inplace (bool): Whether or not the arrays are modified in place. + + Warning: + Array attributes that are not contiguous numpy arrays will be + replaced with contiguous numpy arrays of the appropriate type + if inplace is ``True``. + """ + if inplace: + dataDict = self.cast() + if self.position is not None: + self.position = dataDict['position'] + if self.orientation is not None: + self.orientation = dataDict['orientation'] + if self.typeid is not None: + self.typeid = dataDict['typeid'] + if self.mass is not None: + self.mass = dataDict['mass'] + if self.charge is not None: + self.charge = dataDict['charge'] + if self.diameter is not None: + self.diameter = dataDict['diameter'] + if self.body is not None: + self.body = dataDict['body'] + if self.moment_inertia is not None: + self.moment_inertia = dataDict['moment_inertia'] + if self.velocity is not None: + self.velocity = dataDict['velocity'] + if self.angmom is not None: + self.angmom = dataDict['angmom'] + if self.image is not None: + self.image = dataDict['image'] if self.types is not None and (not len(set(self.types)) == len(self.types)): msg = 'Type names must be unique.' @@ -313,7 +406,33 @@ def __init__(self, M): self._default_value['typeid'] = numpy.uint32(0) self._default_value['group'] = numpy.array([0] * M, dtype=numpy.int32) - def validate(self): + def cast(self): + """Lazy cast attributes to writeable dictionary. + + Convert every array attribute to a `numpy.ndarray` of the proper + type and check that all attributes have the correct dimensions. + + Ignore any attributes that are ``None``. + + Returns: + castedDict (dict): Stored attributes by attribute string name. + """ + castedDict = {} + if self.typeid is not None: + castedDict['typeid'] = numpy.ascontiguousarray( + self.typeid, dtype=numpy.uint32 + ) + castedDict['typeid'] = castedDict['typeid'].reshape([self.N]) + if self.group is not None: + castedDict['group'] = numpy.ascontiguousarray(self.group, dtype=numpy.int32) + castedDict['group'] = castedDict['group'].reshape([self.N, self.M]) + + if self.types is not None and (not len(set(self.types)) == len(self.types)): + msg = 'Type names must be unique.' + raise ValueError(msg) + return castedDict + + def validate(self, inplace=False): """Validate all attributes. Convert every array attribute to a `numpy.ndarray` of the proper @@ -321,16 +440,20 @@ def validate(self): Ignore any attributes that are ``None``. + Args: + inplace (bool): Whether or not the arrays are modified in place. + Warning: Array attributes that are not contiguous numpy arrays will be - replaced with contiguous numpy arrays of the appropriate type. + replaced with contiguous numpy arrays of the appropriate type + if inplace is ``True``. """ - if self.typeid is not None: - self.typeid = numpy.ascontiguousarray(self.typeid, dtype=numpy.uint32) - self.typeid = self.typeid.reshape([self.N]) - if self.group is not None: - self.group = numpy.ascontiguousarray(self.group, dtype=numpy.int32) - self.group = self.group.reshape([self.N, self.M]) + if inplace: + dataDict = self.cast() + if self.typeid is not None: + self.typeid = dataDict['typeid'] + if self.group is not None: + self.group = dataDict['group'] if self.types is not None and (not len(set(self.types)) == len(self.types)): msg = 'Type names must be unique.' @@ -372,7 +495,38 @@ def __init__(self): self._default_value['value'] = numpy.float32(0) self._default_value['group'] = numpy.array([0] * self.M, dtype=numpy.int32) - def validate(self): + def _validate_precision(self, data): + """Maintain floats in numpy arrays.""" + outFloat = numpy.float32 + if isinstance(data, list): + outFloat = numpy.float64 + elif data.dtype == numpy.float64: + outFloat = numpy.float64 + return outFloat + + def cast(self): + """Lazy cast attributes to writeable dictionary. + + Convert every array attribute to a `numpy.ndarray` of the proper + type and check that all attributes have the correct dimensions. + + Ignore any attributes that are ``None``. + + Returns: + castedDict (dict): Stored attributes by attribute string name. + """ + castedDict = {} + if self.value is not None: + castedDict['value'] = numpy.ascontiguousarray( + self.value, dtype=self._validate_precision(self.value) + ) + castedDict['value'] = castedDict['value'].reshape([self.N]) + if self.group is not None: + castedDict['group'] = numpy.ascontiguousarray(self.group, dtype=numpy.int32) + castedDict['group'] = castedDict['group'].reshape([self.N, self.M]) + return castedDict + + def validate(self, inplace=False): """Validate all attributes. Convert every array attribute to a `numpy.ndarray` of the proper @@ -380,16 +534,21 @@ def validate(self): Ignore any attributes that are ``None``. + Args: + inplace (bool): Whether or not the arrays are modified in place. + Warning: Array attributes that are not contiguous numpy arrays will be - replaced with contiguous numpy arrays of the appropriate type. + replaced with contiguous numpy arrays of the appropriate type + if inplace is ``True``. """ - if self.value is not None: - self.value = numpy.ascontiguousarray(self.value, dtype=numpy.float32) - self.value = self.value.reshape([self.N]) - if self.group is not None: - self.group = numpy.ascontiguousarray(self.group, dtype=numpy.int32) - self.group = self.group.reshape([self.N, self.M]) + if inplace: + if self.value is not None: + self.value = numpy.ascontiguousarray(self.value, dtype=numpy.float32) + self.value = self.value.reshape([self.N]) + if self.group is not None: + self.group = numpy.ascontiguousarray(self.group, dtype=numpy.int32) + self.group = self.group.reshape([self.N, self.M]) class Frame: @@ -452,16 +611,25 @@ def __init__(self): 'hpmc/simple_polygon/vertices', ] - def validate(self): - """Validate all contained frame data.""" - self.configuration.validate() - self.particles.validate() - self.bonds.validate() - self.angles.validate() - self.dihedrals.validate() - self.impropers.validate() - self.constraints.validate() - self.pairs.validate() + def validate(self, inplace=False): + """Validate all contained data. + + Args: + inplace (bool): Whether or not the arrays are modified in place. + + Warning: + Array attributes that are not contiguous numpy arrays will be + replaced with contiguous numpy arrays of the appropriate type + if inplace is ``True``. + """ + self.configuration.validate(inplace) + self.particles.validate(inplace) + self.bonds.validate(inplace) + self.angles.validate(inplace) + self.dihedrals.validate(inplace) + self.impropers.validate(inplace) + self.constraints.validate(inplace) + self.pairs.validate(inplace) # validate HPMC state if self.particles.types is not None: @@ -674,17 +842,23 @@ class HOOMDTrajectory: Args: file (`gsd.fl.GSDFile`): File to access. + precision (str): Precision to use for floats on write. + Can be `'single'` or `'double'`. Open hoomd GSD files with `open`. """ - def __init__(self, file): + def __init__(self, file, precision='single'): if file.mode == 'ab': msg = 'Append mode not yet supported' raise ValueError(msg) self._file = file self._initial_frame = None + if precision not in ('single', 'double'): + message = "precision must be 'single' or 'double'" + raise ValueError(message) + self._precision = precision # Used to cache positive results when chunks exist in frame 0. self._chunk_exists_frame_0 = {} @@ -695,7 +869,7 @@ def __init__(self, file): raise RuntimeError('GSD file is not a hoomd schema file: ' + str(self.file)) valid = False version = self.file.schema_version - if version < (2, 0) and version >= (1, 0): + if version < (3, 0) and version >= (1, 0): valid = True if not valid: raise RuntimeError( @@ -705,6 +879,10 @@ def __init__(self, file): + str(self.file) ) + if version < (2, 0) and precision == 'double': + message = "schema 1.x files may are not compatible with precision='double'" + raise RuntimeError(message) + logger.info('found ' + str(len(self)) + ' frames') @property @@ -712,6 +890,11 @@ def file(self): """:class:`gsd.fl.GSDFile`: The file handle.""" return self._file + @property + def precision(self): + """:class:str: The object write precision.""" + return self._precision + def __len__(self): """The number of frames in the trajectory.""" return self.file.nframes @@ -737,6 +920,8 @@ def append(self, frame): if self._initial_frame is None and len(self) > 0: self._read_frame(0) + float_type = numpy.float64 if self.precision == 'double' else numpy.float32 + for path in [ 'configuration', 'particles', @@ -748,9 +933,14 @@ def append(self, frame): 'pairs', ]: container = getattr(frame, path) + dataDict = container.cast() for name in container._default_value: if self._should_write(path, name, frame): - data = getattr(container, name) + data = dataDict.get(name, getattr(container, name)) + if isinstance(data, numpy.ndarray) and numpy.issubdtype( + data.dtype, numpy.floating + ): + data = data.astype(float_type) if name == 'N': data = numpy.array([data], dtype=numpy.uint32) @@ -1062,7 +1252,7 @@ def flush(self): self._file.flush() -def open(name, mode='r'): # noqa: A001 - allow shadowing builtin open +def open(name, mode='r', precision='single'): # noqa: A001 - allow shadowing builtin open """Open a hoomd schema GSD file. The return value of `open` can be used as a context manager. @@ -1070,6 +1260,7 @@ def open(name, mode='r'): # noqa: A001 - allow shadowing builtin open Args: name (str): File name to open. mode (str): File open mode. + precision (str): Float precision to expect in File. Can be 'single' or 'double'. Returns: `HOOMDTrajectory` instance that accesses the file **name** with the @@ -1111,10 +1302,10 @@ def open(name, mode='r'): # noqa: A001 - allow shadowing builtin open mode=mode, application='gsd.hoomd ' + gsd.version.version, schema='hoomd', - schema_version=[1, 4], + schema_version=[2, 0], ) - return HOOMDTrajectory(gsdfileobj) + return HOOMDTrajectory(gsdfileobj, precision=precision) def read_log(name: str, scalar_only=False, glob_pattern='*'): @@ -1171,7 +1362,7 @@ def read_log(name: str, scalar_only=False, glob_pattern='*'): mode='r', application='gsd.hoomd ' + gsd.version.version, schema='hoomd', - schema_version=[1, 4], + schema_version=[2, 0], ) as gsdfileobj: logged_data_names = ( fnfilter(gsdfileobj.find_matching_chunk_names('log/'), glob_pattern) diff --git a/gsd/test/test_hoomd.py b/gsd/test/test_hoomd.py index a9ffb4af..4881f99d 100644 --- a/gsd/test/test_hoomd.py +++ b/gsd/test/test_hoomd.py @@ -214,7 +214,9 @@ def make_nondefault_frame(): frame0 = gsd.hoomd.Frame() frame0.configuration.step = 10000 frame0.configuration.dimensions = 2 - frame0.configuration.box = [4, 5, 6, 1.0, 0.5, 0.25] + frame0.configuration.box = numpy.array( + [4, 5, 6, 1.0, 0.5, 0.25], dtype=numpy.float64 + ) frame0.particles.N = 2 frame0.particles.types = ['A', 'B', 'C'] frame0.particles.type_shapes = [ @@ -222,53 +224,66 @@ def make_nondefault_frame(): {'type': 'Sphere', 'diameter': 3.0}, {'type': 'Sphere', 'diameter': 4.0}, ] - frame0.particles.typeid = [1, 2] - frame0.particles.mass = [2, 3] - frame0.particles.diameter = [3, 4] - frame0.particles.body = [10, 20] - frame0.particles.charge = [0.5, 0.25] - frame0.particles.moment_inertia = [[1, 2, 3], [3, 2, 1]] - frame0.particles.position = [[0.1, 0.2, 0.3], [-1.0, -2.0, -3.0]] - frame0.particles.orientation = [[1, 0.1, 0.2, 0.3], [0, -1.0, -2.0, -3.0]] - frame0.particles.velocity = [[1.1, 2.2, 3.3], [-3.3, -2.2, -1.1]] - frame0.particles.angmom = [[1, 1.1, 2.2, 3.3], [-1, -3.3, -2.2, -1.1]] - frame0.particles.image = [[10, 20, 30], [5, 6, 7]] + frame0.particles.typeid = numpy.array([1, 2]) + frame0.particles.mass = numpy.array([2, 3], dtype=numpy.float64) + frame0.particles.diameter = numpy.array([3, 4], dtype=numpy.float64) + frame0.particles.body = numpy.array([10, 20]) + frame0.particles.charge = numpy.array([0.5, 0.25], dtype=numpy.float64) + frame0.particles.moment_inertia = numpy.array( + [[1, 2, 3], [3, 2, 1]], dtype=numpy.float64 + ) + frame0.particles.position = numpy.array( + [[0.1, 0.2, 0.3], [-1.0, -2.0, -3.0]], dtype=numpy.float64 + ) + frame0.particles.orientation = numpy.array( + [[1, 0.1, 0.2, 0.3], [0, -1.0, -2.0, -3.0]], dtype=numpy.float64 + ) + frame0.particles.velocity = numpy.array( + [[1.1, 2.2, 3.3], [-3.3, -2.2, -1.1]], dtype=numpy.float64 + ) + frame0.particles.angmom = numpy.array( + [[1, 1.1, 2.2, 3.3], [-1, -3.3, -2.2, -1.1]], dtype=numpy.float64 + ) + frame0.particles.image = numpy.array([[10, 20, 30], [5, 6, 7]]) frame0.bonds.N = 1 frame0.bonds.types = ['bondA', 'bondB'] - frame0.bonds.typeid = [1] - frame0.bonds.group = [[0, 1]] + frame0.bonds.typeid = numpy.array([1]) + frame0.bonds.group = numpy.array([[0, 1]]) frame0.angles.N = 1 - frame0.angles.typeid = [2] + frame0.angles.typeid = numpy.array([2]) frame0.angles.types = ['angleA', 'angleB'] - frame0.angles.group = [[0, 1, 0]] + frame0.angles.group = numpy.array([[0, 1, 0]]) frame0.dihedrals.N = 1 - frame0.dihedrals.typeid = [3] + frame0.dihedrals.typeid = numpy.array([3]) frame0.dihedrals.types = ['dihedralA', 'dihedralB'] - frame0.dihedrals.group = [[0, 1, 1, 0]] + frame0.dihedrals.group = numpy.array([[0, 1, 1, 0]]) frame0.impropers.N = 1 - frame0.impropers.typeid = [4] + frame0.impropers.typeid = numpy.array([4]) frame0.impropers.types = ['improperA', 'improperB'] - frame0.impropers.group = [[1, 0, 0, 1]] + frame0.impropers.group = numpy.array([[1, 0, 0, 1]]) frame0.constraints.N = 1 - frame0.constraints.value = [1.1] - frame0.constraints.group = [[0, 1]] + frame0.constraints.value = numpy.array([1.1]) + frame0.constraints.group = numpy.array([[0, 1]]) frame0.pairs.N = 1 frame0.pairs.types = ['pairA', 'pairB'] - frame0.pairs.typeid = [1] - frame0.pairs.group = [[0, 3]] + frame0.pairs.typeid = numpy.array([1]) + frame0.pairs.group = numpy.array([[0, 3]]) - frame0.log['value'] = [1, 2, 4, 10, 12, 18, 22] + frame0.log['value'] = numpy.array([1, 2, 4, 10, 12, 18, 22]) return frame0 -def assert_frames_equal(s, frame0, check_position=True, check_step=True): +def assert_frames_equal( + s, frame0, check_position=True, check_step=True, precision='single' +): """Assert that two frames are equal.""" + dtype = numpy.float32 if precision == 'single' else numpy.float64 if check_step: assert s.configuration.step == frame0.configuration.step @@ -278,22 +293,40 @@ def assert_frames_equal(s, frame0, check_position=True, check_step=True): assert s.particles.types == frame0.particles.types assert s.particles.type_shapes == frame0.particles.type_shapes numpy.testing.assert_array_equal(s.particles.typeid, frame0.particles.typeid) - numpy.testing.assert_array_equal(s.particles.mass, frame0.particles.mass) - numpy.testing.assert_array_equal(s.particles.diameter, frame0.particles.diameter) + numpy.testing.assert_array_equal( + s.particles.mass.astype(dtype), + numpy.asarray(frame0.particles.mass).astype(dtype), + ) + numpy.testing.assert_array_equal( + s.particles.diameter.astype(dtype), + numpy.asarray(frame0.particles.diameter).astype(dtype), + ) numpy.testing.assert_array_equal(s.particles.body, frame0.particles.body) - numpy.testing.assert_array_equal(s.particles.charge, frame0.particles.charge) numpy.testing.assert_array_equal( - s.particles.moment_inertia, frame0.particles.moment_inertia + s.particles.charge.astype(dtype), + numpy.asarray(frame0.particles.charge).astype(dtype), + ) + numpy.testing.assert_array_equal( + s.particles.moment_inertia.astype(dtype), + numpy.asarray(frame0.particles.moment_inertia).astype(dtype), ) if check_position: numpy.testing.assert_array_equal( - s.particles.position, frame0.particles.position + s.particles.position.astype(dtype), + numpy.asarray(frame0.particles.position).astype(dtype), ) numpy.testing.assert_array_equal( - s.particles.orientation, frame0.particles.orientation + s.particles.orientation.astype(dtype), + numpy.asarray(frame0.particles.orientation).astype(dtype), + ) + numpy.testing.assert_array_equal( + s.particles.velocity.astype(dtype), + numpy.asarray(frame0.particles.velocity).astype(dtype), + ) + numpy.testing.assert_array_equal( + s.particles.angmom.astype(dtype), + numpy.asarray(frame0.particles.angmom).astype(dtype), ) - numpy.testing.assert_array_equal(s.particles.velocity, frame0.particles.velocity) - numpy.testing.assert_array_equal(s.particles.angmom, frame0.particles.angmom) numpy.testing.assert_array_equal(s.particles.image, frame0.particles.image) assert s.bonds.N == frame0.bonds.N @@ -317,7 +350,10 @@ def assert_frames_equal(s, frame0, check_position=True, check_step=True): numpy.testing.assert_array_equal(s.impropers.group, frame0.impropers.group) assert s.constraints.N == frame0.constraints.N - numpy.testing.assert_array_equal(s.constraints.value, frame0.constraints.value) + numpy.testing.assert_array_equal( + s.constraints.value, numpy.asarray(frame0.constraints.value).astype(dtype) + ) + # numpy.testing.assert_array_equal(s.constraints.value, frame0.constraints.value) numpy.testing.assert_array_equal(s.constraints.group, frame0.constraints.group) assert s.pairs.N == frame0.pairs.N @@ -326,13 +362,15 @@ def assert_frames_equal(s, frame0, check_position=True, check_step=True): numpy.testing.assert_array_equal(s.pairs.group, frame0.pairs.group) -def test_fallback(tmp_path, open_mode): +@pytest.mark.parametrize('precision', ['single', 'double']) +def test_fallback(tmp_path, open_mode, precision): """Test that properties fall back to defaults when the N changes.""" + floatType = numpy.float32 if precision == 'single' else numpy.float64 frame0 = make_nondefault_frame() frame1 = gsd.hoomd.Frame() frame1.particles.N = 2 - frame1.particles.position = [[-2, -1, 0], [1, 3.0, 0.5]] + frame1.particles.position = numpy.array([[-2, -1, 0], [1, 3.0, 0.5]]) frame1.bonds.N = None frame1.angles.N = None frame1.dihedrals.N = None @@ -355,7 +393,7 @@ def test_fallback(tmp_path, open_mode): frame2.pairs.N = 7 with gsd.hoomd.open( - name=tmp_path / 'test_fallback.gsd', mode=open_mode.write + name=tmp_path / 'test_fallback.gsd', mode=open_mode.write, precision=precision ) as hf: hf.extend([frame0, frame1, frame2]) @@ -363,14 +401,14 @@ def test_fallback(tmp_path, open_mode): assert len(hf) == 3 s = hf[0] - assert_frames_equal(s, frame0) + assert_frames_equal(s, frame0, precision=precision) assert 'value' in s.log numpy.testing.assert_array_equal(s.log['value'], frame0.log['value']) # test that everything but position remained the same in frame 1 s = hf[1] - assert_frames_equal(s, frame0, check_position=False) + assert_frames_equal(s, frame0, check_position=False, precision=precision) assert 'value' in s.log numpy.testing.assert_array_equal(s.log['value'], frame0.log['value']) @@ -385,40 +423,36 @@ def test_fallback(tmp_path, open_mode): s.particles.typeid, numpy.array([0, 0, 0], dtype=numpy.uint32) ) numpy.testing.assert_array_equal( - s.particles.mass, numpy.array([1, 1, 1], dtype=numpy.float32) + s.particles.mass, numpy.array([1, 1, 1], dtype=floatType) ) numpy.testing.assert_array_equal( - s.particles.diameter, numpy.array([1, 1, 1], dtype=numpy.float32) + s.particles.diameter, numpy.array([1, 1, 1], dtype=floatType) ) numpy.testing.assert_array_equal( - s.particles.body, numpy.array([-1, -1, -1], dtype=numpy.float32) + s.particles.body, numpy.array([-1, -1, -1], dtype=floatType) ) numpy.testing.assert_array_equal( - s.particles.charge, numpy.array([0, 0, 0], dtype=numpy.float32) + s.particles.charge, numpy.array([0, 0, 0], dtype=floatType) ) numpy.testing.assert_array_equal( s.particles.moment_inertia, - numpy.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]], dtype=numpy.float32), + numpy.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]], dtype=floatType), ) numpy.testing.assert_array_equal( s.particles.position, - numpy.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]], dtype=numpy.float32), + numpy.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]], dtype=floatType), ) numpy.testing.assert_array_equal( s.particles.orientation, - numpy.array( - [[1, 0, 0, 0], [1, 0, 0, 0], [1, 0, 0, 0]], dtype=numpy.float32 - ), + numpy.array([[1, 0, 0, 0], [1, 0, 0, 0], [1, 0, 0, 0]], dtype=floatType), ) numpy.testing.assert_array_equal( s.particles.velocity, - numpy.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]], dtype=numpy.float32), + numpy.array([[0, 0, 0], [0, 0, 0], [0, 0, 0]], dtype=floatType), ) numpy.testing.assert_array_equal( s.particles.angmom, - numpy.array( - [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=numpy.float32 - ), + numpy.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], dtype=floatType), ) numpy.testing.assert_array_equal( s.particles.image, @@ -481,7 +515,7 @@ def test_fallback(tmp_path, open_mode): assert s.constraints.N == 4 numpy.testing.assert_array_equal( - s.constraints.value, numpy.array([0, 0, 0, 0], dtype=numpy.float32) + s.constraints.value, numpy.array([0, 0, 0, 0], dtype=floatType) ) numpy.testing.assert_array_equal( s.constraints.group, @@ -516,7 +550,7 @@ def test_fallback_to_frame0(tmp_path, open_mode): frame1.pairs.N = None with gsd.hoomd.open( - name=tmp_path / 'test_fallback2.gsd', mode=open_mode.write + name=tmp_path / 'test_fallback2.gsd', mode=open_mode.write, precision='double' ) as hf: hf.extend([frame0, frame1]) @@ -527,7 +561,7 @@ def test_fallback_to_frame0(tmp_path, open_mode): s = hf[1] assert s.configuration.step == frame1.configuration.step - assert_frames_equal(s, frame0, check_step=False) + assert_frames_equal(s, frame0, check_step=False, precision='double') assert 'value' in s.log numpy.testing.assert_array_equal(s.log['value'], frame0.log['value']) @@ -987,7 +1021,9 @@ def test_read_log_warning(tmp_path): def test_initial_frame_copy(tmp_path, open_mode): """Ensure that the user does not unintentionally modify _initial_frame.""" with gsd.hoomd.open( - name=tmp_path / 'test_initial_frame_copy.gsd', mode=open_mode.write + name=tmp_path / 'test_initial_frame_copy.gsd', + mode=open_mode.write, + precision='double', ) as hf: frame = make_nondefault_frame() @@ -1075,3 +1111,131 @@ def test_initial_frame_copy(tmp_path, open_mode): for key in frame_1.log.keys(): assert frame_1.log[key] is initial.log[key] assert not frame_1.log[key].flags.writeable + + +def test_float64(tmp_path, open_mode): + """Test that charges can be optionally set to float32 or float64.""" + frame = gsd.hoomd.Frame() + frame.particles.N = 4 + charges32 = numpy.array([-0.1, 0.1 / 3, 0.1 / 3, 0.1 / 3], dtype=numpy.float32) + charges64 = numpy.array([-0.1, 0.1 / 3, 0.1 / 3, 0.1 / 3], dtype=numpy.float64) + frame.particles.charge = charges64 + numpy.testing.assert_array_equal(frame.particles.charge, charges64) + + with gsd.hoomd.open( + name=tmp_path / 'test_defaults.gsd', mode=open_mode.write, precision='double' + ) as hf: + hf.append(frame) + + with gsd.hoomd.open(name=tmp_path / 'test_defaults.gsd', mode=open_mode.read) as hf: + s = hf[0] + + numpy.testing.assert_array_equal(s.particles.charge, charges64) + numpy.testing.assert_raises( + AssertionError, + numpy.testing.assert_array_equal, + s.particles.charge, + charges32, + ) + + +def test_precision(tmp_path, open_mode): + """Test that single and double precision files can be stored.""" + frame = gsd.hoomd.Frame() + frame.particles.N = 4 + charges32 = numpy.array([-0.1, 0.1 / 3, 0.1 / 3, 0.1 / 3], dtype=numpy.float32) + charges64 = numpy.array([-0.1, 0.1 / 3, 0.1 / 3, 0.1 / 3], dtype=numpy.float64) + frame.particles.charge = charges64 + numpy.testing.assert_array_equal(frame.particles.charge, charges64) + + with gsd.hoomd.open( + name=tmp_path / 'double.gsd', mode=open_mode.write, precision='double' + ) as hf: + hf.append(frame) + + with gsd.hoomd.open( + name=tmp_path / 'single.gsd', mode=open_mode.write, precision='single' + ) as hf: + hf.append(frame) + + with gsd.hoomd.open(name=tmp_path / 'double.gsd', mode=open_mode.read) as hf: + s = hf[0] + + numpy.testing.assert_array_equal(s.particles.charge, charges64) + + with gsd.hoomd.open(name=tmp_path / 'single.gsd', mode=open_mode.read) as hf: + s = hf[0] + numpy.testing.assert_array_equal(s.particles.charge, charges32) + + +def test_all_precision(tmp_path, open_mode): + """Test both single and double precision on default Frame.""" + frame0 = make_nondefault_frame() + with gsd.hoomd.open( + name=tmp_path / 'double.gsd', mode=open_mode.write, precision='double' + ) as hf: + hf.append(frame0) + + with gsd.hoomd.open( + name=tmp_path / 'single.gsd', mode=open_mode.write, precision='single' + ) as hf: + hf.append(frame0) + + with gsd.hoomd.open( + name=tmp_path / 'double.gsd', mode=open_mode.read, precision='double' + ) as hf: + s = hf[0] + for name in ['position', 'orientation', 'velocity', 'angmom', 'charge']: + numpy.testing.assert_array_equal( + getattr(frame0.particles, name), getattr(s.particles, name) + ) + assert getattr(s.particles, name).dtype == numpy.float64 + + with gsd.hoomd.open(name=tmp_path / 'single.gsd', mode=open_mode.read) as hf: + s = hf[0] + for name in ['position', 'orientation', 'velocity', 'angmom', 'charge']: + numpy.testing.assert_array_equal( + numpy.ascontiguousarray( + getattr(frame0.particles, name), dtype=numpy.float32 + ), + getattr(s.particles, name), + ) + assert getattr(s.particles, name).dtype == numpy.float32 + + +def test_write_multiple_precision(tmp_path): + """Test single, then double precision writing on the particles position.""" + frame = gsd.hoomd.Frame() + frame.particles.N = 1 + frame.particles.position = numpy.array( + [[0.1, 0.2, 0.3]], dtype=numpy.float64 + ) # is a list of list of floats + with gsd.hoomd.open( # converts from float64 to float32 + name=tmp_path / 'single.gsd', mode='w', precision='single' + ) as hf: + hf.append(frame) + with gsd.hoomd.open( # converts from python float to dtype float64 + name=tmp_path / 'double.gsd', mode='w', precision='double' + ) as hf: + hf.append(frame) + + # frame is now fixed single precision, due to the conversion in hf.append(Frame) + with gsd.hoomd.open( + name=tmp_path / 'single.gsd', mode='r', precision='single' + ) as hf: + s = hf[0] # read as a single precision + numpy.testing.assert_array_equal( + numpy.ascontiguousarray(frame.particles.position, dtype=numpy.float32), + s.particles.position, + ) + assert s.particles.position.dtype == numpy.float32 + with gsd.hoomd.open( + name=tmp_path / 'double.gsd', mode='r', precision='double' + ) as hf: + s = hf[0] # read as double precision + numpy.testing.assert_array_equal( + # frame.particles.position, + numpy.array([[0.1, 0.2, 0.3]], dtype=numpy.float64), + s.particles.position, + ) + assert s.particles.position.dtype == numpy.float64