#!/usr/bin/env python
"""Contains functions for import and generation of initial conditions for MD simulations."""
from __future__ import print_function
import os
import numpy as np
import matplotlib.pyplot as plt
import json
from . import atoms as at
from . import xyz
from . import constants as c
[docs]def get_initial_random(symbols, coords, temp):
"""
Sample random initial conditions at a given temperature.
Parameters
----------
symbols : list of str
list of atomic symbols of the molecule
coords : np.2darray
atomic coordinates of the molecule (shape: (N, 3))
temp : float
average temperature to be sampled
Returns
-------
symbols : list of str
list of atomic symbols of the molecule
coords : np.2darray
atomic coordinates of the molecule (shape: (N, 3))
velos : np.2darray
atomic velocities of the molecule (shape: (N, 3))
"""
nat = len(symbols)
if nat == 1:
print("WARNING: temperature is not defined for single atoms")
f = 0
elif nat == 2:
f = 1
elif nat > 2:
f = (3*nat)-6
energy = temp*f*c.kb/2
masses = np.array([at.atomicmasses[s.capitalize()] for s in symbols])
xdelta = np.sqrt(2*energy/(3*nat*masses))
p = np.random.normal(size=(nat, 3))
velos = np.vstack((xdelta, xdelta, xdelta)).transpose()*p
return symbols, coords, velos
[docs]def get_temperature(symbols, velos):
"""
Get the temperature corresponding to a set of atomic velocities.
Parameters
----------
symbols : list of str
list of atomic symbols of the molecule
velos : np.2darray
atomic velocities of the molecule (shape: (N, 3))
Returns
-------
temp : float
temperature
"""
nat = len(symbols)
masses = np.array([at.atomicmasses[s.capitalize()] for s in symbols])
ekin = 0.0
for i in range(nat):
ekin += masses[i]*0.5*np.sum(velos[i]**2)
if nat == 1:
print("WARNING: temperature is not defined for single atoms")
elif nat == 2:
f = 1
elif nat > 2:
f = (3*nat)-6
return 2*ekin/(f*c.kb)
[docs]def get_temperature_distribution(symbols, coords, temp, n=10000):
"""
Get the temperature distribution of a set of randomly sampled initial conditions.
Parameters
----------
symbols : list of str
list of atomic symbols of the molecule
coords : np.2darray
atomic coordinates of the molecule (shape: (N, 3))
temp : float
average temperature to be sampled
n : int
number of initial conditions to be sampled to determine the temperature distribution
Returns
-------
temps : list of float
list of temperatures that have been determined for a set of initial conditions
"""
temps = []
for i in range(n):
symb, coord, velos = get_initial_random(symbols, coords, temp)
temps.append(get_temperature(symbols, velos))
return temps
[docs]def plot_temperature_distribution(temps):
"""
Plot a histogram of the sampled temperatures.
Parameters
----------
temps : list of float
list of temperatures that have been determined for a set of initial conditions
"""
plt.hist(temps, bins=50)
plt.show()
[docs]def get_initial_from_log(step, logdir="."):
"""
Get initial coordinates and velocities of a given time step from a previous calculation.
Parameters
----------
step : int
time step from which from which the initial conditions should be taken
logdir : str
path to the results directory that contains the respective logfile (default: '.')
Returns
-------
symbols : list of str
list of atomic symbols of the molecule
coords : np.2darray
atomic coordinates of the molecule (shape: (N, 3))
velos : np.2darray
atomic velocities of the molecule (shape: (N, 3))
"""
filendx = 1
while os.path.exists(logdir+"/dynamics-%d.log" % (filendx)):
filendx += 1
for ndx in range(filendx):
with open(logdir+"/dynamics-%d.log" % ndx) as f:
for line in f:
if "step: %d" % step in line:
while not "coordinates" in line:
line = f.next()
nat = int(line.split()[-1])
symbcoord = np.array([f.next().split() for i in range(nat)])
symbols = list(symbcoord[:, 0])
coord = symbcoord[:, 1:].astype(float)
while not "velocities" in line:
line = f.next()
velos = np.array([f.next().split() for i in range(nat)], dtype=float)
return symbols, coord, velos
[docs]def get_laststep_from_log(logdir="."):
"""
Determine the latest step number that has been performed in an existing dynamics simulation.
Parameters
----------
logdir : str
path to the results directory that contains the respective logfile (default: '.')
Returns
-------
step : int
latest step number
"""
if not os.path.exists(logdir+"/dynamics-0.log"):
return 0
filendx = 1
while os.path.exists(logdir+"/dynamics-%d.log" % (filendx)):
filendx += 1
for line in reversed(open(logdir+"/dynamics-%d.log" % (filendx-1)).readlines()):
if "step: " in line:
return int(line.split()[-1])
# if log-file contains no steps, delete it and try again
os.remove("dynamics-%d.log" % (filendx-1))
return get_laststep_from_log()
[docs]def read_initial():
"""
Read initial coordinates and velocities from files in xyz-format.
Returns
-------
symbols : list of str
list of atomic symbols of the molecule
coords : np.2darray
atomic coordinates of the molecule (shape: (N, 3))
velos : np.2darray
atomic velocities of the molecule (shape: (N, 3))
"""
try:
symbols, coords = xyz.read_xyz("initial.xyz")
symbols, velos = xyz.read_xyz("initial.vel")
return symbols, coords, velos
except IOError:
print("'initial.xyz and 'initial.vel' not found; trying to create new initial conditions.")
with open("meta-config.json") as f:
config = json.load(f)
generate_initial(config)
symbols, coords = xyz.read_xyz("initial.xyz")
symbols, velos = xyz.read_xyz("initial.vel")
return symbols, coords, velos
[docs]def write_initial(symbols, coords, velos):
"""
Write initial coordinates and velocities in xyz-format for use in the dynamics program.
Parameters
----------
symbols : list of str
list of atomic symbols of the molecule
coords : np.2darray
atomic coordinates of the molecule (shape: (N, 3))
velos : np.2darray
atomic velocities of the molecule (shape: (N, 3))
"""
xyz.write_xyz("initial.xyz", symbols, coords, comment="initial coordinates", fmt="%18.9e")
xyz.write_xyz("initial.vel", symbols, velos, comment="initial velocities", fmt="%18.9e")
[docs]def generate_initial(config):
"""
Generate initial coordinates and velocities from the information given in the config dict.
Parameters
----------
config : dict
molecular dynamics configuration dictionary as read from 'meta-config.json'
"""
if "molecular dynamics" in config.keys() and "initial" in config["molecular dynamics"].keys():
config = config["molecular dynamics"]["initial"]
else:
return
if config["source"] == "new":
if "xyz" in config.keys() and "temperature" in config.keys():
symbols, coords = xyz.read_xyz(config["xyz"])
write_initial(*get_initial_random(symbols, coords, config["temperature"]))
elif config["source"] == "log":
if "logdir" in config.keys():
logdir = config["logdir"]
else:
logdir = "."
if "step" in config.keys():
step = config["step"]
else:
step = get_laststep_from_log(logdir)
write_initial(*get_initial_from_log(step, logdir))