tools/i-pi/ipi/utils/io/io_xyz.py

   1 """Contains the functions used to print the trajectories and read input
   2 configurations with xyz formatting.
   3
   4 Copyright (C) 2013, Joshua More and Michele Ceriotti
   5
   6 This program is free software: you can redistribute it and/or modify
   7 it under the terms of the GNU General Public License as published by
   8 the Free Software Foundation, either version 3 of the License, or
   9 (at your option) any later version.
  10
  11 This program is distributed in the hope that it will be useful,
  12 but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  14 GNU General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with this program. If not, see <http.//www.gnu.org/licenses/>.
  18
  19
  20 Functions:
  21    print_xyz_path: Prints all the bead configurations.
  22    print_xyz: Prints the centroid configurations.
  23    read_xyz: Reads the cell parameters and atom configurations from a xyz file.
  24 """
  25
  26 __all__ = ['print_xyz_path', 'print_xyz', 'read_xyz', 'iter_xyz']
  27
  28 import numpy as np
  29 import math, sys
  30 import ipi.utils.mathtools as mt
  31 from ipi.utils.depend import depstrip
  32 from ipi.engine.atoms import Atoms
  33 from ipi.utils.units import *
  34
  35 def print_xyz_path(beads, cell, filedesc = sys.stdout):
  36    """Prints all the bead configurations, into a xyz formatted file.
  37
  38    Prints all the replicas for each time step separately, rather than all at
  39    once.
  40
  41    Args:
  42       beads: A beads object giving the bead positions.
  43       cell: A cell object giving the system box.
  44       filedesc: An open writable file object. Defaults to standard output.
  45    """
  46
  47    a, b, c, alpha, beta, gamma = mt.h2abc_deg(cell.h)
  48
  49    natoms = beads.natoms
  50    nbeads = beads.nbeads
  51    for j in range(nbeads):
  52       filedesc.write("%d\n# bead: %d CELL(abcABC): %10.5f  %10.5f  %10.5f  %10.5f  %10.5f  %10.5f \n" % (natoms, j, a, b, c, alpha, beta, gamma))
  53       for i in range(natoms):
  54          qs = depstrip(beads.q)
  55          lab = depstrip(beads.names)
  56          filedesc.write("%8s %12.5e %12.5e %12.5e\n" % (lab[i], qs[j][3*i], qs[j][3*i+1], qs[j][3*i+2]))
  57
  58 def print_xyz(atoms, cell, filedesc = sys.stdout, title=""):
  59    """Prints the centroid configurations, into a xyz formatted file.
  60
  61    Args:
  62       atoms: An atoms object giving the centroid positions.
  63       cell: A cell object giving the system box.
  64       filedesc: An open writable file object. Defaults to standard output.
  65       title: This gives a string to be appended to the comment line.
  66    """
  67
  68    a, b, c, alpha, beta, gamma = mt.h2abc_deg(cell.h)
  69
  70    natoms = atoms.natoms
  71    filedesc.write("%d\n# CELL(abcABC): %10.5f  %10.5f  %10.5f  %10.5f  %10.5f  %10.5f  %s\n" % ( natoms, a, b, c, alpha, beta, gamma, title))
  72    # direct access to avoid unnecessary slow-down
  73    qs = depstrip(atoms.q)
  74    lab = depstrip(atoms.names)
  75    for i in range(natoms):
  76       filedesc.write("%8s %12.5e %12.5e %12.5e\n" % (lab[i], qs[3*i], qs[3*i+1], qs[3*i+2]))
  77
  78 def read_xyz(filedesc):
  79    """Takes a xyz-style file and creates an Atoms object.
  80
  81    Args:
  82       filedesc: An open readable file object from a xyz formatted file.
  83
  84    Returns:
  85       An Atoms object with the appropriate atom labels, masses and positions.
  86    """
  87
  88    natoms = filedesc.readline()
  89    if natoms == "":
  90       raise EOFError("The file descriptor hit EOF.")
  91    natoms = int(natoms)
  92    comment = filedesc.readline()
  93
  94    qatoms = []
  95    names = []
  96    masses = []
  97    iat = 0
  98    while (iat < natoms):
  99       body = filedesc.readline()
 100       if body.strip() == "":
 101          break
 102       body = body.split()
 103       name = body[0]
 104       names.append(name)
 105       masses.append(Elements.mass(name))
 106       x = float(body[1])
 107       y = float(body[2])
 108       z = float(body[3])
 109       qatoms.append(x)
 110       qatoms.append(y)
 111       qatoms.append(z)
 112       iat += 1
 113
 114    if natoms != len(names):
 115       raise ValueError("The number of atom records does not match the header of the xyz file.")
 116
 117    atoms = Atoms(natoms)
 118 #   for i in range(natoms):
 119 #      nat = atoms[i]
 120 #      nat.q = qatoms[i]
 121 #      nat.name = names[i]
 122 #      nat.m = Elements.mass(names[i])
 123    atoms.q = np.asarray(qatoms)
 124    atoms.names = np.asarray(names, dtype='|S4')
 125    atoms.m = np.asarray(masses)
 126
 127    return atoms
 128
 129 def iter_xyz(filedesc):
 130    """Takes a xyz-style file and yields one Atoms object after another.
 131
 132    Args:
 133       filedesc: An open readable file object from a xyz formatted file.
 134
 135    Returns:
 136       Generator over the xyz trajectory, that yields
 137       Atoms objects with the appropriate atom labels, masses and positions.
 138    """
 139
 140    try:
 141       while 1:
 142          atoms = read_xyz(filedesc)
 143          yield atoms
 144    except EOFError:
 145       pass