Source code for MDAnalysis.coordinates.FHIAIMS

# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the GNU Public Licence, v2 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#

"""
FHI-AIMS file format --- :mod:`MDAnalysis.coordinates.FHIAIMS`
==============================================================

Classes to read and write `FHI-AIMS`_ coordinate files.

The cell vectors are specified by the (optional) lines with the
``lattice_vector`` tag::

    lattice_vector x  y  z

where x, y, and z are expressed in ångström (Å).

.. Note::

   In the original `FHI-AIMS format`_, up to three lines with
   ``lattice_vector`` are allowed (order matters) where the absent line implies
   no periodicity in that direction.  In MDAnalysis, only the case of no
   ``lattice_vector`` or three ``lattice_vector`` lines are allowed.

Atomic positions and names are specified either by the ``atom`` or by the
``atom_frac`` tags::

    atom           x  y  z  name
    atom_frac      nx ny nz name

where x, y, and z are expressed in ångström, and nx, ny and nz are real numbers
in [0, 1] and are used to compute the atomic positions in units of the basic
cell.

Atomic velocities can be added on the line right after the corresponding
``atom`` in units of Å/ps using the ``velocity`` tag::

    velocity      vx vy vz


The field name is a string identifying the atomic species.  See also the
specifications in the official `FHI-AIMS format`_.

Classes
-------

.. autoclass:: Timestep
   :members:
   :inherited-members:

.. autoclass:: FHIAIMSReader
   :members:
   :inherited-members:

.. autoclass:: FHIAIMSWriter
   :members:
   :inherited-members:

Developer notes: ``FHIAIMSWriter`` format strings
-------------------------------------------------

The :class:`FHIAIMSWriter` class has a :attr:`FHIAIMSWriter.fmt`
attribute, which is a dictionary of different strings for writing
lines in ``.in`` files.  These are as follows:

``xyz``
  An atom line without velocities.  Requires that the `name` and
  `pos` keys be supplied.  E.g.::

     fmt['xyz'].format(pos=(0.0, 1.0, 2.0), name='O')

``vel``
  An line that specifies velocities::

     fmt['xyz'].format(vel=(0.1, 0.2, 0.3))

``box_triclinic``
  The (optional) initial lines of the file which gives box dimensions.
  Requires the `box` keyword, as a length 9 vector.  This is a flattened
  version of the (3, 3) triclinic vector representation of the unit
  cell.


.. Links

.. _FHI-AIMS: https://aimsclub.fhi-berlin.mpg.de/
.. _FHI-AIMS format: https://doi.org/10.6084/m9.figshare.12413477.v1

"""
from __future__ import absolute_import

import re
from six.moves import range, zip
from six import raise_from

import itertools
import warnings

import numpy as np

from . import base
from .core import triclinic_box, triclinic_vectors
from ..exceptions import NoDataError
from ..lib import util
from ..lib import mdamath


[docs]class Timestep(base.Timestep): def _init_unitcell(self): return np.zeros((3, 3), dtype=np.float32) @property def dimensions(self): """unitcell dimensions (A, B, C, alpha, beta, gamma)""" return triclinic_box(self._unitcell[0], self._unitcell[1], self._unitcell[2]) @dimensions.setter def dimensions(self, new): self._unitcell[:] = triclinic_vectors(new)
[docs]class FHIAIMSReader(base.SingleFrameReaderBase): """Reader for the FHIAIMS geometry format. Single frame reader for the `FHI-AIMS`_ input file format. Reads geometry (3D and molecules only), positions (absolut or fractional), velocities if given, all according to the `FHI-AIMS format`_ specifications """ format = ['IN', 'FHIAIMS'] units = {'time': 'ps', 'length': 'Angstrom', 'velocity': 'Angstrom/ps'} _Timestep = Timestep def _read_first_frame(self): with util.openany(self.filename, 'rt') as fhiaimsfile: relative, positions, velocities, lattice_vectors = [], [], [], [] skip_tags = ["#", "initial_moment"] oldline = '' for line in fhiaimsfile: line = line.strip() if line.startswith("atom"): positions.append(line.split()[1:-1]) relative.append('atom_frac' in line) oldline = line continue if line.startswith("velocity"): if not 'atom' in oldline: raise ValueError( 'Non-conforming line (velocity must follow atom): ({0})in FHI-AIMS input file {0}'.format(line, self.filename)) velocities.append(line.split()[1:]) oldline = line continue if line.startswith("lattice"): lattice_vectors.append(line.split()[1:]) oldline = line continue if any([line.startswith(tag) for tag in skip_tags]): oldline = line continue raise ValueError( 'Non-conforming line: ({0})in FHI-AIMS input file {0}'.format(line, self.filename)) # positions and velocities are lists of lists of strings; they will be # cast to np.arrays(..., dtype=float32) during assignment to ts.positions/ts.velocities lattice_vectors = np.asarray(lattice_vectors, dtype=np.float32) if len(velocities) not in (0, len(positions)): raise ValueError( 'Found incorrect number of velocity tags ({0}) in the FHI-AIMS file, should be {1}.'.format( len(velocities), len(positions))) if len(lattice_vectors) not in (0, 3): raise ValueError( 'Found partial periodicity in FHI-AIMS file. This cannot be handled by MDAnalysis.') if len(lattice_vectors) == 0 and any(relative): raise ValueError( 'Found relative coordinates in FHI-AIMS file without lattice info.') # create Timestep self.n_atoms = n_atoms = len(positions) self.ts = ts = self._Timestep(n_atoms, **self._ts_kwargs) ts.positions = positions if len(lattice_vectors) > 0: ts._unitcell[:] = lattice_vectors ts.positions[relative] = np.matmul( ts.positions[relative], lattice_vectors) if len(velocities) > 0: ts.velocities = velocities self.ts.frame = 0 # 0-based frame number
[docs] def Writer(self, filename, n_atoms=None, **kwargs): """Returns a FHIAIMSWriter for *filename*. Parameters ---------- filename: str filename of the output FHI-AIMS file Returns ------- :class:`FHIAIMSWriter` """ if n_atoms is None: n_atoms = self.n_atoms return FHIAIMSWriter(filename, n_atoms=n_atoms, **kwargs)
[docs]class FHIAIMSWriter(base.WriterBase): """FHI-AIMS Writer. Single frame writer for the `FHI-AIMS`_ format. Writes geometry (3D and molecules only), positions (absolut only), velocities if given, all according to the `FHI-AIMS format`_ specifications. If no atom names are given, it will set each atom name to "X". """ format = ['IN', 'FHIAIMS'] units = {'time': None, 'length': 'Angstrom', 'velocity': 'Angstrom/ps'} #: format strings for the FHI-AIMS file (all include newline) fmt = { # coordinates output format, see https://doi.org/10.6084/m9.figshare.12413477.v1 'xyz': "atom {pos[0]:12.8f} {pos[1]:12.8f} {pos[2]:12.8f} {name:<3s}\n", 'vel': "velocity {vel[0]:12.8f} {vel[1]:12.8f} {vel[2]:12.8f}\n", # unitcell 'box_triclinic': "lattice_vector {box[0]:12.8f} {box[1]:12.8f} {box[2]:12.8f}\nlattice_vector {box[3]:12.8f} {box[4]:12.8f} {box[5]:12.8f}\nlattice_vector {box[6]:12.8f} {box[7]:12.8f} {box[8]:12.8f}\n" } def __init__(self, filename, convert_units=True, n_atoms=None, **kwargs): """Set up the FHI-AIMS Writer Parameters ----------- filename : str output filename n_atoms : int (optional) number of atoms """ self.filename = util.filename(filename, ext='.in', keep=True) self.n_atoms = n_atoms def _write_next_frame(self, obj): """Write selection at current trajectory frame to file. Parameters ----------- obj : AtomGroup or Universe """ # write() method that complies with the Trajectory API # TODO 2.0: Remove timestep logic try: # make sure to use atoms (Issue 46) ag_or_ts = obj.atoms # can write from selection == Universe (Issue 49) except AttributeError: if isinstance(obj, base.Timestep): ag_or_ts = obj.copy() else: raise_from( TypeError("No Timestep found in obj argument"), None) # Check for topology information missing_topology = [] try: names = ag_or_ts.names except (AttributeError, NoDataError): names = itertools.cycle(('X',)) missing_topology.append('names') try: atom_indices = ag_or_ts.ids except (AttributeError, NoDataError): atom_indices = range(1, ag_or_ts.n_atoms+1) missing_topology.append('ids') if missing_topology: warnings.warn( "Supplied AtomGroup was missing the following attributes: " "{miss}. These will be written with default values. " "".format(miss=', '.join(missing_topology))) positions = ag_or_ts.positions try: velocities = ag_or_ts.velocities has_velocities = True except (AttributeError, NoDataError): has_velocities = False with util.openany(self.filename, 'wt') as output_fhiaims: # Lattice try: # for AtomGroup/Universe tri_dims = obj.universe.trajectory.ts.triclinic_dimensions except AttributeError: # for Timestep tri_dims = obj.triclinic_dimensions # full output if np.any(tri_dims != 0): output_fhiaims.write( self.fmt['box_triclinic'].format(box=tri_dims.flatten())) # Atom descriptions and coords # Dont use enumerate here, # all attributes could be infinite cycles! for atom_index, name in zip( range(ag_or_ts.n_atoms), names): output_fhiaims.write(self.fmt['xyz'].format( pos=positions[atom_index], name=name)) if has_velocities: output_fhiaims.write(self.fmt['vel'].format( vel=velocities[atom_index]))