sef43 / openmm-mdanalysis-reporter

MIT License
9 stars 4 forks source link

MDAReporter cannot report to H5MD #5

Open orionarcher opened 5 months ago

orionarcher commented 5 months ago

I discovered this while investigating issue #4.

Expected behavior

MDAReporter successfully creates a H5MD file.

Actual behavior

MDAReporter is failing within the H5MD writer.

Traceback:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[70], line 13
     11 simulation = interchange.to_openmm_simulation(integrator=integrator)
     12 simulation.reporters.append(MDAReporter('traj.h5md',1))
---> 13 simulation.step(1)

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/openmm/app/simulation.py:147, in Simulation.step(self, steps)
    145 def step(self, steps):
    146     """Advance the simulation by integrating a specified number of time steps."""
--> 147     self._simulate(endStep=self.currentStep+steps)

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/openmm/app/simulation.py:247, in Simulation._simulate(self, endStep, endTime)
    244 # Generate the reports.
    246 if len(wrapped) > 0:
--> 247     self._generate_reports(wrapped, True)
    248 if len(unwrapped) > 0:
    249     self._generate_reports(unwrapped, False)

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/openmm/app/simulation.py:269, in Simulation._generate_reports(self, reports, periodic)
    265 state = self.context.getState(getPositions=getPositions, getVelocities=getVelocities, getForces=getForces,
    266                               getEnergy=getEnergy, getParameters=True, enforcePeriodicBox=periodic,
    267                               groups=self.context.getIntegrator().getIntegrationForceGroups())
    268 for reporter, next in reports:
--> 269     reporter.report(self, state)

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/mdareporter/mdareporter.py:126, in MDAReporter.report(self, simulation, state)
    123 self._mdaUniverse.dimensions[:3] = _sanitize_box_angles(self._mdaUniverse.dimensions[:3])
    125 # write to the trajectory file
--> 126 self._mdaWriter.write(self._atomGroup)
    128 self._nextModel += 1

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/MDAnalysis/coordinates/base.py:1609, in WriterBase.write(self, obj)
   1591 def write(self, obj):
   1592     """Write current timestep, using the supplied `obj`.
   1593 
   1594     Parameters
   (...)
   1607        removed. Use AtomGroup or Universe as an input instead.
   1608     """
-> 1609     return self._write_next_frame(obj)

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/MDAnalysis/coordinates/H5MD.py:1145, in H5MDWriter._write_next_frame(self, ag)
   1143     self._determine_units(ag)
   1144     self._open_file()
-> 1145     self._initialize_hdf5_datasets(ts)
   1147 return self._write_next_timestep(ts)

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/MDAnalysis/coordinates/H5MD.py:1284, in H5MDWriter._initialize_hdf5_datasets(self, ts)
   1282 if self.data_keys:
   1283     for key in self.data_keys:
-> 1284         self._create_observables_dataset(key, ts.data[key])

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/MDAnalysis/coordinates/H5MD.py:1350, in H5MDWriter._create_observables_dataset(self, group, data)
   1348 # guarantee ints and floats have a shape ()
   1349 data = np.asarray(data)
-> 1350 self._obsv.require_dataset(f'{group}/value',
   1351                            shape=(0,) + data.shape,
   1352                            maxshape=(None,) + data.shape,
   1353                            dtype=data.dtype)
   1354 if 'step' not in self._obsv[group]:
   1355     self._obsv[f'{group}/step'] = self._step

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/h5py/_hl/group.py:275, in Group.require_dataset(self, name, shape, dtype, exact, **kwds)
    273 with phil:
    274     if name not in self:
--> 275         return self.create_dataset(name, *(shape, dtype), **kwds)
    277     if isinstance(shape, int):
    278         shape = (shape,)

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/h5py/_hl/group.py:183, in Group.create_dataset(self, name, shape, dtype, data, **kwds)
    180         parent_path, name = name.rsplit(b'/', 1)
    181         group = self.require_group(parent_path)
--> 183 dsid = dataset.make_new_dset(group, shape, dtype, data, name, **kwds)
    184 dset = dataset.Dataset(dsid)
    185 return dset

File ~/miniconda3/envs/atomate_emmet/lib/python3.11/site-packages/h5py/_hl/dataset.py:86, in make_new_dset(parent, shape, dtype, data, name, chunks, compression, shuffle, fletcher32, maxshape, compression_opts, fillvalue, scaleoffset, track_times, external, track_order, dcpl, dapl, efile_prefix, virtual_prefix, allow_unknown_filter, rdcc_nslots, rdcc_nbytes, rdcc_w0)
     84     else:
     85         dtype = numpy.dtype(dtype)
---> 86     tid = h5t.py_create(dtype, logical=1)
     88 # Legacy
     89 if any((compression, shuffle, fletcher32, maxshape, scaleoffset)) and chunks is False:

File h5py/h5t.pyx:1658, in h5py.h5t.py_create()

File h5py/h5t.pyx:1682, in h5py.h5t.py_create()

File h5py/h5t.pyx:1742, in h5py.h5t.py_create()

TypeError: Object dtype dtype('O') has no native HDF5 equivalent

Code to reproduce the behavior

from openmm.app import *
from openmm import *
from openmm.unit import *
from mdareporter.data.files import VILLIN_PDB

pdb = PDBFile(VILLIN_PDB)
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(
    pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds
)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()

# relevant code that is failing
simulation.reporters.append(MDAReporter('traj.h5md',1))
simulation.step(1)

Current environment

orionarcher commented 5 months ago

Fixed by MDanalysis PR here

sef43 commented 5 months ago

Will leave this open untill: