#!/usr/bin/env python
"""Contains a class for the analysis of molecular dynamics trajectories."""
import numpy as np
import os
import json
from .. import md
from ..cvs import cvfunctions as cvf
from .. import constants as c
from .. import xyz
[docs]class analyze_trajectory(md.moleculardynamics):
def __init__(self):
md.moleculardynamics.__init__(self)
self.read_energies = False
self.filendx = 0
while os.path.exists("dynamics-%s.log" % str(self.filendx+1)):
self.filendx += 1
self.energytypes = ["pot", "kin", "temp", "cons"]
self.structuretypes = ["bond", "angle", "torsion", "cn"]
self.read_input()
def _parse_input(self):
"""Prepare analysis run from the information in self.config."""
analyze_cfg = self.get_cfg(self.config, "analyze", default={})
self.stride = self.get_cfg(analyze_cfg, "stride", default=1)
self.start = self.get_cfg(analyze_cfg, "start", default=0)
self.end = self.get_cfg(analyze_cfg, "end", default=-1)
self.nlog = self.get_cfg(analyze_cfg, "nlog", default=-1)
self.energies = self.get_cfg(analyze_cfg, "properties", default=[])
self.structure = self.get_cfg(analyze_cfg, "structure", default=[])
if len(self.energies):
self.read_energies = True
properties_dicts = self.energies + self.structure
self.properties = self.energies + [[self.get_cfg(prop, "type")] for prop in self.structure]
for i, prop in enumerate(self.properties):
if isinstance(prop, list) and prop[0] in ["bond", "angle", "torsion"]:
prop += properties_dicts[i]["atoms"]
if isinstance(prop, list) and prop[0] == "cn":
atom = properties_dicts[i]["atom"]
ref = properties_dicts[i]["reference"]
m = properties_dicts[i]["m"]
n = properties_dicts[i]["n"]
d = properties_dicts[i]["d"]
prop += [atom, n, m, d, ref]
[docs] def get_properties(self):
"""
Return all properties that have been specified to be analyzed.
See also
--------
add_property, remove_property, reset_properties
"""
props = []
if "analyze" in self.config.keys():
if "properties" in self.config["analyze"].keys():
props += self.config["analyze"]["properties"]
if "structure" in self.config["analyze"].keys():
props += self.config["analyze"]["structure"]
return props
[docs] def add_property(self, ptype, **kwargs):
"""
Add a property to be analyzed along the calculated trajectory.
Parameters
----------
ptype : str
property type to be analyzed (see below for all possibilities)
kwargs
see below
Keyword Arguments
-----------------
atoms : list
atom indices that specify the atoms to be used for bond, angle or torsion analysis
atom : int
atom index for the calculation of the coordination number
reference : list of int
indices of reference atoms to be included (default: all atoms other than ``atom``)
n : int
first exponent
m : int
second exponent
d : 'auto' or float or list of float
single distance or list of distances for all atoms in ``reference`` (default: 'auto')
Notes
-----
The following properties are available for analysis:
============== ============================================================================
Property types Explanation
============== ============================================================================
pot potential energy
kin kinetic energy
temp temperature
cons energy conservation compared to last step (only NVE simulations)
bond bond length (requires ``atoms`` kwarg)
angle bond angle (requires ``atoms`` kwarg)
torsion torsion angle (requires ``atoms`` kwarg)
cn coordination number (requires ``atom``, ``reference``, ``m``, ``n``, and
``d`` kwargs)
============== ============================================================================
Example
-------
In order to analyze the bond length between the atoms 0 and 1, do:
>>> from metafalcon.analysis.analyze import analyze_trajectory
>>> anl = analyze_trajectory()
>>> anl.read_input()
>>> anl.add_property("bond", atoms=[0, 1])
>>> anl.analyze()
See Also
--------
remove_property, reset_properties, get_properties
"""
if ptype in self.energytypes:
self._append_enproperty(ptype)
elif ptype in self.structuretypes:
self._append_strproperty(ptype, **kwargs)
self.save_input()
def _append_enproperty(self, ptype):
if "analyze" not in self.config.keys():
self.config["analyze"] = {}
if "properties" not in self.config["analyze"].keys():
self.config["analyze"]["properties"] = []
if ptype not in self.config["analyze"]["properties"]:
self.config["analyze"]["properties"].append(ptype)
def _append_strproperty(self, ptype, **kwargs):
if "analyze" not in self.config.keys():
self.config["analyze"] = {}
if "structure" not in self.config["analyze"].keys():
self.config["analyze"]["structure"] = []
strprop = kwargs
strprop["type"] = ptype
if strprop not in self.config["analyze"]["structure"]:
self.config["analyze"]["structure"].append(strprop)
[docs] def remove_property(self, ptype, **kwargs):
"""
Remove a property from the properties to be analyzed along the calculated trajectory. If
the given parameters are not unique, all matching properties will be removed.
Parameters
----------
ptype : str
property type to be removed (see below for all possibilities)
kwargs
see below
Keyword Arguments
-----------------
atoms : list
atom indices that specify the atoms to be used for bond, angle or torsion analysis
atom : int
atom index for the calculation of the coordination number
reference : list of int
indices of reference atoms to be included (default: all atoms other than ``atom``)
n : int
first exponent
m : int
second exponent
d : 'auto' or float or list of float
single distance or list of distances for all atoms in ``reference`` (default: 'auto')
Notes
-----
The following properties are available for analysis:
============== ============================================================================
Property types Explanation
============== ============================================================================
pot potential energy
kin kinetic energy
temp temperature
cons energy conservation compared to last step (only NVE simulations)
bond bond length (requires ``atoms`` kwarg)
angle bond angle (requires ``atoms`` kwarg)
torsion torsion angle (requires ``atoms`` kwarg)
cn coordination number (requires ``atom``, ``reference``, ``m``, ``n``, and
``d`` kwargs)
============== ============================================================================
Examples
--------
Suppose you have added two bond lengths and a bond angle for analysis:
>>> from metafalcon.analysis.analyze import analyze_trajectory
>>> anl = analyze_trajectory()
>>> anl.read_input()
>>> anl.add_property("bond", atoms=[0, 1])
>>> anl.add_property("bond", atoms=[0, 2])
>>> anl.add_property("angle", atoms=[0, 1, 2])
You can remove a single property specifically
>>> anl.remove_property("bond", atoms=[0, 2])
or remove all bond lengths:
>>> anl.remove_property("bond")
See Also
--------
add_property, reset_properties, get_properties
"""
if ptype in self.energytypes:
props = [prop for prop in self.get_properties() if prop == ptype]
# for i, oldprop in enumerate(self.config["analyze"]["properties"]):
# if self.config["analyze"]["properties"] == delprop:
# self.config["analyze"]["properties"].pop(i)
for prop in props:
self.config["analyze"]["properties"].remove(prop)
elif ptype in self.structuretypes:
props = [prop for prop in self.get_properties() if isinstance(prop, dict) and
prop["type"] == ptype]
for key, value in kwargs.items():
props = [prop for prop in props if key in prop.keys() and prop[key] == value]
for prop in props:
self.config["analyze"]["structure"].remove(prop)
self.save_input()
[docs] def reset_properties(self):
"""
Remove all properties that have been configured before.
See also
--------
add_property, remove_property, get_properties
"""
if "analyze" in self.config.keys():
self.config["analyze"]["properties"] = []
self.config["analyze"]["structure"] = []
self.save_input()
def _set_config(self, key, value):
if "analyze" not in self.config.keys():
self.config["analyze"] = {}
self.config["analyze"][key] = value
self.save_input()
[docs] def set_stride(self, stride):
"""
Set the stride attribute to speed up the analysis.
Parameters
----------
stride : int
number of frames to be skipped after each analysis step (default: 1)
"""
self.stride = stride
self._set_config("stride", stride)
[docs] def set_start(self, start):
"""
Set the first dynamics step to be analyzed.
Parameters
----------
start : int
first dynamics step to be analyzed (default: 0)
"""
self.start = start
self._set_config("start", start)
[docs] def set_end(self, end):
"""
Set the last dynamics step to be analyzed.
Parameters
----------
end : int
last dynamics step to be analyzed (default: last)
"""
self.end = end
self._set_config("end", end)
[docs] def set_nlog(self, nlog):
"""
Set the log-file to be analyzed.
Parameters
----------
nlog : int
log-file index to be analyzed
"""
self.nlog = nlog
[docs] def set_coordinates(self, coords):
"""
Set the coordinates attribute.
Parameters
----------
coords : np.2darray
molecular cartesian coordinates (shape: (N, 3))
"""
self.coords = coords
def _get_coordlist(self):
"""Read a list of coordinates from the log-file."""
self.coordlist = []
self.enlist = []
self.steps = []
if self.nlog < 0:
indices = range(self.filendx+1)
else:
indices = [self.nlog]
for i in indices:
with open("dynamics-%i.log" % i) as f:
for line in f:
if "step:" in line:
step = int(line.split()[1])
if step % self.stride or step in self.steps or step < self.start or (step > self.end and self.end > 0):
continue
self.steps.append(step)
[f.next() for a in range(3)]
tmp = np.array(
[f.next().split()[1:] for at in range(self.nat)],
dtype=float)
self.coordlist.append(tmp)
if self.read_energies:
[f.next() for a in range(self.nat+3)]
self.enlist.append([f.next() for en in range(5)])
def _get_nat(self):
"""Read number of atoms from the dynamics-\*.log file."""
with open("dynamics-0.log") as f:
for line in f:
if "coordinates" in line:
return int(line.split()[-1])
def _get_bond(self, ndx):
"""Return bond length for the given atom indices."""
p = [cvf.bond(coords[ndx]) for coords in self.coordlist]
return np.array(p)*c.a0
def _get_bond_single(self, ndx):
"""Return bond length for the given atom indices."""
return cvf.bond(self.coords[ndx])
def _get_angle(self, ndx):
"""Return bond angle for the given atom indices."""
p = [cvf.angle(coords[ndx]) for coords in self.coordlist]
return np.array(p)
def _get_angle_single(self, ndx):
"""Return bond angle for the given atom indices."""
return cvf.angle(self.coords[ndx])
def _get_torsion(self, ndx):
"""Return torsion angle for the given atom indices."""
p = [cvf.torsion(coords[ndx]) for coords in self.coordlist]
return np.array(p)
def _get_torsion_single(self, ndx):
"""Return torsion angle for the given atom indices."""
return cvf.torsion(self.coords[ndx])
def _get_cn(self, prm):
"""Return coordination number evaluated with parameters `prm`."""
p = [cvf.cn(self.symbols, coords, *prm)[0] for coords in self.coordlist]
return np.array(p)
def _get_cn_single(self, prm):
"""Return coordination number evaluated with parameters `prm`."""
return cvf.cn(self.symbols, self.coords, *prm)[0]
def _get_custom(self, prm):
"""Return custom variable evaluated with parameters `prm`."""
p = [md.meta.cvc.customcv(self.symbols, coords, prm) for coords in self.coordlist]
return np.array(p)
def _get_custom_single(self, prm):
"""Return custom variable evaluated with parameters `prm`."""
return md.meta.cvc.customcv(self.symbols, self.coords, prm)
def _get_grep(self, en):
"""Return properties grepped from the log file."""
ndx = {"temp": 0, "kin": 1, "pot": 2, "total": 3, "cons": 4}
p = [float(e[ndx[en]].split()[-1]) for e in self.enlist]
return np.array(p)
def _write_result(self, prop, result):
"""Write analyze_\*.dat file that contains the result."""
if prop[0] == "cn":
prop = prop[:2]
if isinstance(prop, list):
prop = "-".join([str(p) for p in prop])
with open("analyze_"+prop+".dat", "w") as f:
f.write("### ANALYZE TRAJECTORY OUTPUT ###\n")
f.write("# units are fs, angstrom, degrees, kelvin, hartree\n")
f.write("# time, "+prop+"\n")
fmt = "%8.2f %16.8f\n"
for s, r in zip(self.steps, result):
f.write(fmt % (s*self.tstep/c.fs2au, r))
def _write_init(self):
"""Initialize analyze_\*.dat file that will contain the result."""
for i, prop in enumerate(self.properties):
if isinstance(prop, list):
prop = "-".join(prop)
if not isinstance(prop, basestring):
with open("analyze_"+prop+".dat", "w") as f:
f.write("### ANALYZE TRAJECTORY OUTPUT ###\n")
f.write("# units are fs, angstrom, degrees, kelvin, hartree\n")
f.write("# time, "+prop+"\n")
def _write_result_single(self, prop, result, step):
"""Write results to existing analyze_\*.dat file."""
if isinstance(prop, list):
prop = "-".join(prop)
with open("analyze_"+prop+".dat", "a") as f:
fmt = "%8.2f %16.8f\n"
f.write(fmt % (step*self.tstep/c.fs2au, result))
def _write_xyz(self):
"""Write dynamics.xyz file."""
if self.nlog < 0:
filename = "dynamics.xyz"
else:
filename = "dynamics-%d.xyz" % self.nlog
comments = ["step: %i, t: %.2f fs" % (s, s*self.tstep/c.fs2au) for s in self.steps]
xyz.write_trajectory(filename, self.symbols, self.coordlist, comments=comments)
[docs] def log2xyz(self):
"""Wrapper function for the coordinate conversion from the log-file to `dynamics.xyz`."""
self._parse_input()
self._get_coordlist()
self._write_xyz()
[docs] def analyze(self):
"""Wrapper function for the trajectory analysis of structure and energy properties."""
self._parse_input()
self._get_coordlist()
for i, prop in enumerate(self.properties):
if prop[0] == "bond":
p = self._get_bond(prop[1:])
elif prop[0] == "angle":
p = self._get_angle(prop[1:])
elif prop[0] == "torsion":
p = self._get_torsion(prop[1:])
elif prop[0] == "cn":
p = self._get_cn(prop[1:])
elif prop[0] == "custom":
prm = json.loads(" ".join(prop[1:]))
p = self._get_custom(prm)
elif isinstance(prop, basestring):
p = self._get_grep(prop)
self._write_result(prop, p)
def _analyze_single(self, step):
"""Wrapper function for the trajectory analysis of structure and energy properties."""
for i, prop in enumerate(self.properties):
if prop[0] == "bond":
p = self._get_bond_single([int(p) for p in prop[1:]])
elif prop[0] == "angle":
p = self._get_angle_single([int(p) for p in prop[1:]])
elif prop[0] == "torsion":
p = self._get_torsion_single([int(p) for p in prop[1:]])
elif prop[0] == "cn":
p = self._get_cn_single([int(p) for p in prop[1:]]+self.cnprm[i])
elif prop[0] == "custom":
prm = json.loads(" ".join(prop[1:]))
p = self._get_custom_single(prm)
if not isinstance(prop, basestring):
self._write_result_single(prop, p, step)