1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331
|
# -*- 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 Lesser GNU Public Licence, v2.1 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:: 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
"""
import re
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
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"}
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.dimensions = triclinic_box(*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
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)
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
.. versionchanged:: 2.0.0
Support for deprecated Timestep argument has now been removed.
Use AtomGroup or Universe as an input instead.
"""
# write() method that complies with the Trajectory API
try:
# make sure to use atoms (Issue 46)
ag = obj.atoms
# can write from selection == Universe (Issue 49)
except AttributeError:
errmsg = "Input obj is neither an AtomGroup or Universe"
raise TypeError(errmsg) from None
# Check for topology information
missing_topology = []
try:
names = ag.names
except (AttributeError, NoDataError):
names = itertools.cycle(("X",))
missing_topology.append("names")
try:
atom_indices = ag.ids
except (AttributeError, NoDataError):
atom_indices = range(1, ag.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.positions
try:
velocities = ag.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 tri_dims is not None:
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.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])
)
|