#!/usr/bin/env python
"""Contains useful functions for the work with xyz-files and cartesian coordinates."""
import numpy as np
from . import constants as c
[docs]def read_xyz(filename):
"""
Read single structure from xyz-file and return symbols and coordinates
Parameters
----------
filename : string
path to xyz-file to be read
Returns
-------
symbols : list of strings
element symbols
coords : numpy array (shape: (N, 3))
coordinates in atomic units
"""
with open(filename) as f:
nat = int(f.readline())
symcoords = [line.split() for line in f.readlines()[1:nat+1]]
symbols = [sc[0] for sc in symcoords]
coords = np.array([sc[1:] for sc in symcoords], dtype=float)/c.a0
return symbols, coords
[docs]def read_trajectory(filename, first=None, last=None):
"""
Read trajectory from xyz-file and return symbols and coordinates
Parameters
----------
filename : string
path to xyz-file to be read
first : int
first frame to be read (default: 0)
last : int
last frame to be read (default: -1)
Returns
-------
symbols : list of lists of strings
element symbols
coords : list of numpy arrays (shape: (N, 3))
coordinates in atomic units
"""
with open(filename) as f:
lines = f.readlines()
nat = int(lines[0])
nsteps = len(lines)/(nat+2)
symbols = []
coords = []
for i in range(nsteps):
symcoords = [line.split() for line in lines[i*(nat+2)+2:(i+1)*(nat+2)]]
symbols.append([sc[0] for sc in symcoords])
coords.append(np.array([sc[1:4] for sc in symcoords], dtype=float)/c.a0)
return symbols[first:last], coords[first:last]
[docs]def write_trajectory(filename, symbols, coords, comments=[], fmt="%12.8f", mode="w"):
"""
Write trajectory to xyz-file
Parameters
----------
filename : string
path to xyz-file to be saved
symbols : list of strings
element symbols
coords : list of numpy arrays (shape: (N, 3))
coordinates in atomic units
comments : list of str
comments to be inserted into the second line of every structure
fmt : str
formatstring to be used for float numbers in the file
mode : str
file opening mode (w for new file, a for append to existing; default: w
"""
nat = len(symbols)
with open(filename, mode) as f:
for i, coord in enumerate(coords):
coord = np.copy(coord)*c.a0
f.write("%d\n" % nat)
if len(comments)>i and comments[i]:
f.write("%s\n" % comments[i])
else:
f.write("\n")
for s, crd in zip(symbols, coord):
line = "%2s "+fmt+" "+fmt+" "+fmt+" "+"\n"
f.write(line % (s, crd[0], crd[1], crd[2]))
[docs]def write_xyz(filename, symbols, coords, comment="", fmt="%12.8f", mode="w"):
"""
Write single structure to xyz-file
Parameters
----------
filename : string
path to xyz-file to be saved
symbols : list of strings
element symbols
coords : numpy array (shape: (N, 3))
coordinates in atomic units
comment : str
comment to be inserted into the second line of the xyz-file
fmt : str
formatstring to be used for float numbers in the file
mode : str
file opening mode (w for new file, a for append to existing; default: w)
"""
coords = np.copy(coords)*c.a0
with open(filename, mode) as f:
f.write("%d\n" % len(symbols))
f.write("%s\n" % comment)
for i, s in enumerate(symbols):
line = "%2s "+fmt+" "+fmt+" "+fmt+" "+"\n"
f.write(line % (s, coords[i, 0], coords[i, 1], coords[i, 2]))
[docs]def rot_x(coords, angle):
"""
Rotate molecule around x-axis
Parameters
----------
coords : numpy array (shape: (N, 3))
coordinates of the molecule
angle : float
angle in radians
Returns
-------
coords : numpy array (shape: (N, 3))
rotated coordinates of the molecule
"""
rotmat = np.identity(3)
rotmat[1, 1] = np.cos(angle)
rotmat[1, 2] = -np.sin(angle)
rotmat[2, 1] = np.sin(angle)
rotmat[2, 2] = np.cos(angle)
return np.dot(coords, rotmat)
[docs]def rot_y(coords, angle):
"""
Rotate molecule around y-axis
Parameters
----------
coords : numpy array (shape: (N, 3))
coordinates of the molecule
angle : float
angle in radians
Returns
-------
coords : numpy array (shape: (N, 3))
rotated coordinates of the molecule
"""
rotmat = np.identity(3)
rotmat[0, 0] = np.cos(angle)
rotmat[0, 2] = np.sin(angle)
rotmat[2, 0] = -np.sin(angle)
rotmat[2, 2] = np.cos(angle)
return np.dot(coords, rotmat)
[docs]def rot_z(coords, angle):
"""
Rotate molecule around z-axis
Parameters
----------
coords : numpy array (shape: (N, 3))
coordinates of the molecule
angle : float
angle in radians
Returns
-------
coords : numpy array (shape: (N, 3))
rotated coordinates of the molecule
"""
rotmat = np.identity(3)
rotmat[0, 0] = np.cos(angle)
rotmat[0, 1] = -np.sin(angle)
rotmat[1, 0] = np.sin(angle)
rotmat[1, 1] = np.cos(angle)
return np.dot(coords, rotmat)
[docs]def rot_rodrigues(coords, angle, n, center=np.zeros(3)):
"""
Rotate molecule around arbitrary axis n
Parameters
----------
coords : numpy array (shape: (N, 3))
coordinates of the molecule
angle : float
angle in radians
n : numpy array (shape: 3)
normed rotation vector
center : numpy array (shape: 3)
rotation center (default: origin)
Returns
-------
coords : numpy array (shape: (N, 3))
rotated coordinates of the molecule
"""
coords -= center
rotmat = np.zeros((3, 3))
cos1 = 1-np.cos(angle)
cos = np.cos(angle)
sin = np.sin(angle)
rotmat[0, 0] = n[0]**2*cos1+cos
rotmat[0, 1] = n[0]*n[1]*cos1-n[2]*sin
rotmat[0, 2] = n[0]*n[2]*cos1+n[1]*sin
rotmat[1, 0] = n[1]*n[0]*cos1+n[2]*sin
rotmat[1, 1] = n[1]**2*cos1+cos
rotmat[1, 2] = n[1]*n[2]*cos1-n[0]*sin
rotmat[2, 0] = n[2]*n[0]*cos1-n[1]*sin
rotmat[2, 1] = n[2]*n[1]*cos1+n[0]*sin
rotmat[2, 2] = n[2]**2*cos1+cos
coords = np.dot(coords, rotmat)
coords += center
return coords
[docs]def trans_x(coords, t):
"""
Translate molecule in x-direction
Parameters
----------
coords : numpy array (shape: (N, 3))
coordinates of the molecule
t : float
translation
Returns
-------
coords : numpy array (shape: (N, 3))
translated coordinates of the molecule
"""
coords[:, 0] += t
return coords
[docs]def trans_y(coords, t):
"""
Translate molecule in y-direction
Parameters
----------
coords : numpy array (shape: (N, 3))
coordinates of the molecule
t : float
translation
Returns
-------
coords : numpy array (shape: (N, 3))
translated coordinates of the molecule
"""
coords[:, 1] += t
return coords
[docs]def trans_z(coords, t):
"""
Translate molecule in z-direction
Parameters
----------
coords : numpy array (shape: (N, 3))
coordinates of the molecule
t : float
translation
Returns
-------
coords : numpy array (shape: (N, 3))
translated coordinates of the molecule
"""
coords[:, 2] += t
return coords