.. -*- coding: utf-8 -*- .. _auxiliary: =============== Auxiliary files =============== Auxiliary readers allow you to read in timeseries data accompanying a trajectory that is not stored in the regular trajectory file. Supported formats ----------------- +-------------------------+--------+-----------+---------------------------------------+ | Reader | Format | Extension | Remarks | | | | (if file) | | +=========================+========+===========+=======================================+ | :class:`.XVGReader` | XVG | xvg | Produced by Gromacs during simulation | | | | | or analysis. | | | | (default) | Reads full file on initialisation. | +-------------------------+--------+-----------+---------------------------------------+ | :class:`.XVGFileReader` | XVG-F | xvg | Alternate xvg file reader, reading | | | | | each step from the file in turn for a | | | | | lower memory footprint. | | | | | :class:`.XVGReader` is the default | | | | | reader for .xvg files. | +-------------------------+--------+-----------+---------------------------------------+ | :class:`.EDRReader` | EDR | edr | Produced by Gromacs during simulation | | | | | or analysis. | | | | | Reads full file on initialisation. | +-------------------------+--------+-----------+---------------------------------------+ XVG Files --------- Reading data directly --------------------- .. ipython:: python import MDAnalysis as mda from MDAnalysis.tests.datafiles import XVG_BZ2 # cobrotoxin protein forces aux = mda.auxiliary.core.auxreader(XVG_BZ2) aux In stand-alone use, an auxiliary reader allows you to iterate over each step in a set of auxiliary data. .. ipython:: python for step in aux: print(step.data) Use slicing to skip steps. .. ipython:: python for step in aux[1:2]: print(step.time) The :func:`~MDAnalysis.auxiliary.core.auxreader` function uses the :func:`~MDAnalysis.auxiliary.core.get_auxreader_for` to return an appropriate class. This can guess the format either from a filename, ' .. ipython:: python mda.auxiliary.core.get_auxreader_for(XVG_BZ2) or return the reader for a specified format. .. ipython:: python mda.auxiliary.core.get_auxreader_for(format='XVG-F') Loading data into a Universe ---------------------------- Auxiliary data may be added to a trajectory Reader through the :meth:`~MDAnalysis.coordinates.base.ProtoReader.add_auxiliary` method. Auxiliary data may be passed in as a AuxReader instance, or directly as e.g. a filename, in which case :func:`~MDAnalysis.auxiliary.core.get_auxreader_for` is used to guess an appropriate reader. .. ipython:: python :okwarning: from MDAnalysis.tests.datafiles import PDB_xvf, TRR_xvf u = mda.Universe(PDB_xvf, TRR_xvf) u.trajectory.add_auxiliary('protein_force', XVG_BZ2) for ts in u.trajectory: print(ts.aux.protein_force) Passing arguments to auxiliary data ----------------------------------- For alignment with trajectory data, auxiliary readers provide methods to assign each auxiliary step to the nearest trajectory timestep, read all steps assigned to a trajectory timestep and calculate 'representative' value(s) of the auxiliary data for that timestep. To set a timestep or ?? 'Assignment' of auxiliary steps to trajectory timesteps is determined from the time of the auxiliary step, ``dt`` of the trajectory and time at the first frame of the trajectory. If there are no auxiliary steps assigned to a given timestep (or none within ``cutoff``, if set), the representative value(s) are set to ``np.nan``. Iterating over auxiliary data ----------------------------- Auxiliary data may not perfectly line up with the trajectory, or have missing data. .. ipython:: python :okwarning: from MDAnalysis.tests.datafiles import PDB, TRR u_long = mda.Universe(PDB, TRR) u_long.trajectory.add_auxiliary('protein_force', XVG_BZ2, dt=200) for ts in u_long.trajectory: print(ts.time, ts.aux.protein_force[:4]) The trajectory :class:`~MDAnalysis.coordinates.base.ProtoReader` methods :meth:`~MDAnalysis.coordinates.base.ProtoReader.next_as_aux` and :meth:`~MDAnalysis.coordinates.base.ProtoReader.iter_as_aux` allow for movement through only trajectory timesteps for which auxiliary data is available. .. ipython:: python for ts in u_long.trajectory.iter_as_aux('protein_force'): print(ts.time, ts.aux.protein_force[:4]) This may be used to avoid representative values set to ``np.nan``, particularly when auxiliary data is less frequent. Sometimes the auxiliary data is longer than the trajectory. .. ipython:: python :okwarning: u_short = mda.Universe(PDB) u_short.trajectory.add_auxiliary('protein_force', XVG_BZ2) for ts in u_short.trajectory: print(ts.time, ts.aux.protein_force) In order to acess auxiliary values at every individual step, including those outside the time range of the trajectory, :meth:`~MDAnalysis.coordinates.base.ProtoReader.iter_auxiliary` allows iteration over the auxiliary independent of the trajectory. .. ipython:: python :okwarning: for step in u_short.trajectory.iter_auxiliary('protein_force'): print(step.data) To iterate over only a certain section of the auxiliary: .. ipython:: python :okwarning: for step in u_short.trajectory.iter_auxiliary('protein_force', start=1, step=2): # every 2nd step from 1 print(step.time) The trajectory remains unchanged, and the auxiliary will be returned to the current timestep after iteration is complete. Accessing auxiliary attributes ------------------------------ To check the values of attributes of an added auxiliary, use :meth:`~MDAnalysis.coordinates.base.ProtoReader.get_aux_attribute`. .. ipython:: python u.trajectory.get_aux_attribute('protein_force', 'dt') If attributes are settable, they can be changed using :meth:`~MDAnalysis.coordinates.base.ProtoReader.set_aux_attribute`. .. ipython:: python u.trajectory.set_aux_attribute('protein_force', 'data_selector', [1]) The auxiliary may be renamed using ``set_aux_attribute`` with 'auxname', or more directly by using :meth:`~MDAnalysis.coordinates.base.ProtoReader.rename_aux`. .. ipython:: python u.trajectory.ts.aux.protein_force u.trajectory.rename_aux('protein_force', 'f') u.trajectory.ts.aux.f Recreating auxiliaries ---------------------- To recreate an auxiliary, the set of attributes necessary to replicate it can first be obtained with :meth:`~MDAnalysis.auxiliary.base.AuxReader.get_description`. The returned dictionary can then be passed to :func:`~MDAnalysis.auxiliary.core.auxreader` to load a new copy of the original auxiliary reader. .. ipython:: python description = aux.get_description() list(description.keys()) del aux mda.auxiliary.core.auxreader(**description) The 'description' of any or all the auxiliaries added to a trajectory can be obtained using :meth:`~MDAnalysis.coordinates.base.ProtoReader.get_aux_descriptions`. .. ipython:: python descriptions = u.trajectory.get_aux_descriptions(['f']) To reload, pass the dictionary into :meth:`~MDAnalysis.coordinates.base.ProtoReader.add_auxiliary`. .. ipython:: python :okwarning: u2 = mda.Universe(PDB, TRR) for desc in descriptions: u2.trajectory.add_auxiliary(**desc) EDR Files --------- EDR files are created by GROMACS during simulations and contain additional non-trajectory time-series data of the system, such as energies, temperature, or pressure. The EDRReader allows direct reading of these binary files and associating of the data with trajectory time steps just like the XVGReader does. As all functionality of the base AuxReaders (see XVGReader above) also works with the EDRReader, this section will highlight functionality unique to the EDRReader. Standalone Usage ---------------- The EDRReader is initialised by passing the path to an EDR file as an argument. .. ipython:: python :okwarning: import MDAnalysis as mda from MDAnalysisTests.datafiles import AUX_EDR aux = mda.auxiliary.EDR.EDRReader(AUX_EDR) # Or, for example # aux = mda.auxiliary.EDR.EDRReader("path/to/file/ener.edr") Dozens of terms can be defined in EDR files. A list of available data is conveniently available unter the `terms` attribute. .. ipython:: python aux.terms To extract data for plotting, the `get_data` method can be used. It can be used to extract a single, multiple, or all terms as follows: .. ipython:: python temp = aux.get_data("Temperature") print(temp.keys()) energies = aux.get_data(["Potential", "Kinetic En."]) print(energies.keys()) all_data = aux.get_data() print(all_data.keys()) The times of each data point are always part of the returned dictionary to facilicate plotting. .. ipython:: python import matplotlib.pyplot as plt plt.plot(temp["Time"], temp["Temperature"]) plt.ylabel("Temperature [K]") plt.xlabel("Time [ps]") plt.show() .. plot:: import matplotlib.pyplot as plt import MDAnalysis as mda from MDAnalysisTests.datafiles import AUX_EDR aux = mda.auxiliary.EDR.EDRReader(AUX_EDR) temp = aux.get_data("Temperature") plt.plot(temp["Time"], temp["Temperature"]) plt.ylabel("Temperature [K]") plt.xlabel("Time [ps]") plt.show() Unit Handling ------------- On object creation, the EDRReader creates a ``unit_dict`` attribute which contains information on the units of the data stored within. These units are read from the EDR file automatically and by default converted to MDAnalysis base units where such units are defined. The automatic unit conversion can be disabled by setting the ``convert_units`` kwarg to False. .. ipython:: python :okwarning: aux.unit_dict["Box-X"] aux_native = mda.auxiliary.EDR.EDRReader(AUX_EDR, convert_units=False) aux_native.unit_dict["Box-X"] Use with Trajectories --------------------- An arbitrary number of terms to be associated with a trajectory can be specified as a dictionary. The dictionary is chosen so that the name to be used in MDAnalysis to access the data is mapped to the name in `aux.terms`. .. ipython:: python :okwarning: from MDAnalysisTests.datafiles import AUX_EDR_TPR, AUX_EDR_XTC u = mda.Universe(AUX_EDR_TPR, AUX_EDR_XTC) term_dict = {"temp": "Temperature", "epot": "Potential"} u.trajectory.add_auxiliary(term_dict, aux) u.trajectory.ts.aux.epot Adding all data is possible by omitting the dictionary as follows. It is then necessary to specify that the EDRReader should be passed as the `auxdata` argument as such: .. ipython:: python :okwarning: u = mda.Universe(AUX_EDR_TPR, AUX_EDR_XTC) u.trajectory.add_auxiliary(auxdata=aux) u.trajectory.aux_list Selecting Trajectory Frames Based on Auxiliary Data --------------------------------------------------- One use case for the new auxiliary readers is the selection of frames based on auxiliary data. To select only those frames with a potential energy below a certain threshold, the following can be used: .. ipython:: python :okwarning: u = mda.Universe(AUX_EDR_TPR, AUX_EDR_XTC) term_dict = {"epot": "Potential"} u.trajectory.add_auxiliary(term_dict, aux) selected_frames = np.array([ts.frame for ts in u.trajectory if ts.aux.epot < -524600]) A slice of the trajectory can then be obtained from the list of frames as such: .. ipython:: python trajectory_slice = u.trajectory[selected_frames] print(len(u.trajectory)) print(len(trajectory_slice)) Memory Usage ------------ It is assumed that the EDR file is small enough to be read in full. However, since one EDRReader instance is created for each term added to a trajectory, memory usage monitoring was implemented. A warning will be issued if 1 GB of memory is used by auxiliary data. This warning limit can optionally be changed by defining the `memory_limit` (in bytes) when adding data to a trajectory. Below, the memory limit is set to 200 MB. .. ipython:: python :okwarning: u = mda.Universe(AUX_EDR_TPR, AUX_EDR_XTC) term_dict = {"temp": "Temperature", "epot": "Potential"} u.trajectory.add_auxiliary(term_dict, aux, memory_limit=2e+08)