#!/usr/bin/env python
"""Contains a class for molecular dynamics simulation runs."""
from __future__ import print_function
import sys
import os
import datetime
import numpy as np
import json
from . import initial
from . import meta
from . import constants as c
from . import atoms as at
from .interfaces.types import types as ifaces
[docs]class moleculardynamics:
"""
Class to carry out molecular dynamics simulations using the velocity-Verlet algorithm.
Attributes
----------
symbols : list of str
atomic symbols of the molecule
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
velos : np.2darray
cartesian velocities of the molecule (shape: (N, 3))
accel : np.2darray
cartesian accelerations determined dy some quantum chemistry code (shape: (N, 3))
masses : np.array
masses of the atoms in the molecule in atomic mass units (shape: (N))
energy : float
electronic energy in the current time step
filendx : int
index of the current dynamics-\*.log file
restart : bool
if `True`, the previous results are read in and extended (default: False)
start : int
first dynamics step to start with when doing a restart
dometa : bool
modify the forces in each the dynamics step by the metadynamics scheme (default: False)
updatemeta : bool
if `False`, the existing metadynamics potential is not updated (default: True)
"""
def __init__(self):
"""Construct object of class `moleculardynamics`."""
self.symbols = []
self.coords = []
self.velos = []
self.accel = []
self.masses = []
self.energy = 0.0
self.filendx = 0
# self.restart = False
self.start = 0
self.dometa = False
self.updatemeta = True
[docs] def get_cfg(self, config, key, default=None):
try:
return config[key]
except KeyError:
if default is not None:
return default
else:
print("Error: mandatory keyword %s missing in the input file!" % key)
[docs] def set_coordinates(self, coords):
"""
Set coordinates attribute.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
"""
self.coords = coords
[docs] def set_velocities(self, velos):
"""
Set velocities attribute.
Parameters
----------
velos : np.2darray
cartesian velocities of the molecule (shape: (N, 3))
"""
self.velos = velos
[docs] def get_centerofmass(self, coords):
"""
Return center of mass of the given molecule.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
Returns
-------
centerofmass : np.array
x-, y-, z-coordinates of the center of mass (shape: (3))
"""
return np.sum((coords.T*self.masses).T/self.totalmass, axis=0)
[docs] def shift_to_centerofmass(self, coords):
"""
Translate molecule into its center of mass.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
Returns
-------
newcoords : np.2darray
shifted cartesian coordinates (shape: (N, 3))
"""
return coords - self.get_centerofmass(coords)
[docs] def momentum(self, velos):
"""
Return the total momentum.
Parameters
----------
velos : np.2darray
cartesian velocities of the molecule (shape: (N, 3))
Returns
-------
mom : np.array
x-, y-, z-components of the total momentum of the molecule (shape: (3))
"""
return np.sum((velos.T*self.masses).T, axis=0)
[docs] def eliminate_trans(self, velos, mom):
"""
Eliminate translation from velocities.
Parameters
----------
velos : np.2darray
cartesian velocities of the molecule (shape: (N, 3))
mom : np.array
x-, y-, z-components of the total momentum of the molecule (shape: (3))
Returns
-------
newvelos : np.2darray
new velocities with eliminated translation
"""
return velos - mom.T/self.totalmass
[docs] def angular_momentum(self, coords, velos):
"""
Return the total angular momentum.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
velos : np.2darray
cartesian velocities of the molecule (shape: (N, 3))
Returns
-------
angmom : np.array
x-, y-, z-components of the total angular momentum of the molecule (shape: (3))
"""
return np.sum((np.cross(coords, velos, axis=1).T*self.masses).T, axis=0)
[docs] def moment_of_inertia(self, coords):
"""
Return the moment of inertia.
Parameters
----------
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
Returns
-------
inertia : np.array
x-, y-, z-components of the moment of inertia of the molecule (shape: (3))
"""
inertia = np.zeros((3, 3), dtype=float)
Emat = np.identity(3, dtype=float)
for i, mass in enumerate(self.masses):
inertia += mass*(np.inner(coords[i], coords[i])*Emat-np.outer(coords[i], coords[i]))
return inertia
[docs] def get_angular_velocity(self, angmom, inertia):
"""
Return the angular velocity.
Parameters
----------
angmom : np.array
x-, y-, z-components of the total angular momentum of the molecule (shape: (3))
inertia : np.array
x-, y-, z-components of the moment of inertia of the molecule (shape: (3))
Returns
-------
angvel : np.array
x-, y-, z-components of the angular velocity of the molecule (shape: (3))
"""
return np.dot(np.linalg.inv(inertia), angmom)
[docs] def eliminate_rot(self, velos, angvel, coords):
"""
Eliminate rotation from velocities.
Parameters
----------
velos : np.2darray
cartesian velocities of the molecule (shape: (N, 3))
angvel : np.array
x-, y-, z-components of the angular velocity of the molecule (shape: (3))
coords : np.2darray
cartesian coordinates of the molecule (shape: (N, 3))
Returns
-------
newvelos : np.2darray
new velocities with eliminated rotation (shape: (N, 3))
"""
for i in range(self.nat):
velos[i] = velos[i]-np.cross(angvel, coords[i])
return velos
[docs] def get_ekin(self, velos):
"""
Return the total kinetic energy.
Parameters
----------
velos : np.2darray
cartesian velocities of the molecule (shape: (N, 3))
Returns
-------
ekin : float
kinetic energy
"""
return 0.5*np.sum(((velos**2).T*self.masses).T)
[docs] def get_temperature(self, ekin):
"""
Return the temperature corresponding to the kinetic energy.
Parameters
----------
ekin : float
current kinetic energy
Returns
-------
temp : float
current temperature
"""
if self.nat > 2:
f = 3*self.nat-6
else:
f = 3*self.nat-5
return (2*ekin)/(f*c.kb)
[docs] def set_energy_and_accel(self):
"""Use the QC-calculator to determine the current energy and accelerations."""
self.calculator.set_coords(self.coords)
self.calculator.calculate()
self.energy = self.calculator.get_energy()
forces = -self.calculator.get_gradient()
self.accel = (forces.transpose()/np.array(self.masses)).transpose()
[docs] def initiate_trajectory(self):
"""Initialize trajectory for a molecular dynamics simulation."""
self.step = 0
self.write_initial_step()
self.coords = self.shift_to_centerofmass(self.coords)
self.velos = self.eliminate_trans(self.velos, self.momentum(self.velos))
angmom = self.angular_momentum(self.coords, self.velos)
inertia = self.moment_of_inertia(self.coords)
angvel = self.get_angular_velocity(angmom, inertia)
self.velos = self.eliminate_rot(self. velos, angvel, self.coords)
self.calculator = ifaces[self.iface](self.symbols, self.coords, **self.iface_cfg)
self.set_energy_and_accel()
ekin = self.get_ekin(self.velos)
self.write_step(self.start, ekin, self.energy, self.energy+ekin, 0.0)
### debugging hook
try:
sys.path.insert(0, os.getcwd())
self.debug = __import__("debug_hook")
except:
self.debug = False
### end debugging hook
[docs] def write_initial_step(self):
"""Initialize the dynamics-\*.log file."""
if self.restart:
while os.path.exists("dynamics-%s.log" % self.filendx):
self.filendx += 1
else:
os.system("rm -f dynamics-*.log")
with open("dynamics-%s.log" % self.filendx, "w") as f:
f.write("=================================\n")
f.write("=== MOLECULAR DYNAMICS OUTPUT ===\n")
f.write("=================================\n\n")
f.write("=== INITIAL SYSTEM ===\n")
[docs] def write_step(self, step, ekin, epot, etot, conservation):
"""Write results of current dynamics step to the dynamics-\*.log file."""
time = step*self.tstep/c.fs2au
with open("dynamics-%s.log" % self.filendx, "a") as f:
f.write("\nstep: %i\n" % step)
f.write("time: %.2f fs\n" % time)
f.write("date: %s\n" % datetime.datetime.now().strftime("%Y-%m-%d %H:%M"))
f.write("coordinates: %i\n" % self.nat)
for symb, coord in zip(self.symbols, self.coords):
fmt = "%2s %16.9e %16.9e %16.9e\n"
f.write(fmt % (symb, coord[0], coord[1], coord[2]))
f.write("\nvelocities: %i\n" % self.nat)
for velo in self.velos:
fmt = " %16.9e %16.9e %16.9e\n"
f.write(fmt % (velo[0], velo[1], velo[2]))
fmt = "%s %20.12f\n"
f.write(fmt % ("\nTemperature : ", self.get_temperature(ekin)))
f.write(fmt % ("Kinetic energy : ", ekin))
f.write(fmt % ("Electronic energy : ", epot))
f.write(fmt % ("Total energy : ", etot))
f.write(fmt % ("Conservation : ", conservation))
# ### checkpoint files
# np.save("coords.npy", self.coords)
# np.save("velos.npy", self.velos)
[docs] def verlet(self):
"""Perform molecular dynamics simulation using the velocity-Verlet algorithm."""
self.initiate_trajectory()
### metadynamics ###
if self.dometa:
metaob = meta.metadynamics(symbols=self.symbols, restart=self.restart,
update=self.updatemeta)
metaob.read_input()
metaob.initialize_vg()
if self.iface == "dftbaby":
metaob.set_dftb_object(self.calculator.dftbpes)
metaforces = metaob.get_forces(self.coords, self.start)
self.accel += (metaforces.transpose()/self.masses).transpose()
### end metadynamics ###
with open("dynamics-%s.log" % self.filendx, "a") as f:
f.write("\n=== DYNAMICS BEGIN ===\n")
self.coords = self.coords+self.tstep*self.velos+0.5*((self.tstep)**2)*self.accel
self.coords = self.shift_to_centerofmass(self.coords)
ekin = self.get_ekin(self.velos)
for i in range(self.start+1, self.nsteps+1):
self.step = i
old_accel = self.accel
old_etot = self.energy+ekin
self.set_energy_and_accel()
### metadynamics ###
if self.dometa:
print("original accelerations step %i: \n" % i, self.accel)
oldaccel = np.copy(self.accel)
metaforces = metaob.get_forces(self.coords, i)
self.accel += (metaforces.transpose()/self.masses).transpose()
print("metadynamics accelerations step %i: \n" % i,
self.accel-oldaccel)
print("total accelerations step %i: \n" % i, self.accel)
### end metadynamics ###
### end dynamics if calculation has not converged
if not self.calculator.converged:
print("QM calculation not converged, end dynamics")
break
### end end dynamics
self.velos = self.velos+self.tstep*0.5*(old_accel+self.accel)
self.velos = self.eliminate_trans(self.velos, self.momentum(self.velos))
angmom = self.angular_momentum(self.coords, self.velos)
inertia = self.moment_of_inertia(self.coords)
angvel = self.get_angular_velocity(angmom, inertia)
self.velos = self.eliminate_rot(self.velos, angvel, self.coords)
ekin = self.get_ekin(self.velos)
### thermostat
if self.thermostat != "none":
t = self.get_temperature(ekin)
if self.thermostat == "berendsen":
self.velos *= np.sqrt(1+(self.tstep/self.tau)*(self.tref/t-1))
elif self.thermostat == "bussi":
if self.tau > 0.01:
efunc = np.exp(-self.tstep/self.tau)
else:
efunc = 0.
f = 3*self.nat-6
r = np.random.normal(size=f)
self.velos *= np.sqrt(efunc+self.tref/(f*t)*(1-efunc)*np.sum(r**2)+
2*efunc*np.sqrt(self.tref/(f*t)*(1-efunc))*r[0])
### end thermostat
ekin = self.get_ekin(self.velos)
etot = ekin+self.energy
self.write_step(i, ekin, self.energy, etot, etot-old_etot)
self.coords = self.coords+self.tstep*self.velos+((self.tstep)**2)*self.accel*0.5
self.coords = self.shift_to_centerofmass(self.coords)
### debugging hook
if self.debug:
self.debug.run(i)
### end debugging hook
with open("dynamics-%s.log" % self.filendx,"a") as f:
f.write("\n=== DYNAMICS END ===")
# def modbeeman(self):
# self.initiate_trajectory()
# with open("dynamics-%s.log" % self.filendx,"a") as output:
# output.write("\n=== DYNAMICS BEGIN ===\n")
# #---- r0,v0
# ekin=self.get_ekin(self.velos)
# self.coords=self.coords+self.tstep*self.velos+float(0.5)*((self.tstep)**2)*self.accel #F0,r1
# self.coords=self.shift_to_centerofmass(self.coords)
# oldforces=self.accel
# oldtotal_energy=self.energy+ekin
# self.energy, self.accel=self.set_energy_and_accel() #F1
# self.velos=self.velos+self.tstep*float(0.5)*(oldforces+self.accel) #v1
# self.velos=self.eliminate_trans(self.velos,self.momentum(self.velos))
# angmom=self.angular_momentum(self.coords,self.velos)
# inertia=self.moment_of_inertia(self.coords,self.masses)
# angvel=self.get_angular_velocity(angmom,inertia)
# self.velos=self.eliminate_rot(self.velos,angvel,self.coords)
# ekin=self.get_ekin(self.velos)
# total_energy=ekin+self.energy
# self.write_step(self.start+1,self.symbols,self.coords,self.velos,ekin,self.energy,total_energy,total_energy-oldtotal_energy)
#
# for i in range(self.start+2,self.nsteps+1): # Beeman requires r1,v1,F0,F1
# olderforces=oldforces # Fn-1
# oldforces=self.accel # Fn
# oldtotal_energy=self.energy+ekin
# self.coords=self.coords+self.tstep*self.velos+((self.tstep)**2)*(oldforces*float(2.0/3.0)-float(1.0/6.0)*olderforces) # r2 with Beeman using Fn and Fn-1
# self.coords=self.shift_to_centerofmass(self.coords)
# self.energy, self.accel=self.set_energy_and_accel()
# # Beeman: vn+1=vn + 1/3*dt*F(n+1) + 5/6*dt*Fn - 1/6*dt*Fn-1
# self.velos=self.velos+float(1.0/3.0)*self.tstep*self.accel+float(5.0/6.0)*self.tstep*oldforces-float(1.0/6.0)*self.tstep*olderforces # vn+1
# self.velos=self.eliminate_trans(self.velos,self.momentum(self.velos))
# angmom=self.angular_momentum(self.coords,self.velos)
# inertia=self.moment_of_inertia(self.coords,self.masses)
# angvel=self.get_angular_velocity(angmom,inertia)
# self.velos=self.eliminate_rot(self.velos,angvel,self.coords)
# angmom=self.angular_momentum(self.coords,self.velos)
# ekin=self.get_ekin(self.velos)
# total_energy=ekin+self.energy
# self.write_step(i,self.symbols,self.coords,self.velos,ekin,self.energy,total_energy,total_energy-oldtotal_energy)
if __name__ == "__main__":
import argparse
parser = argparse.ArgumentParser(description="Start molecular dynamics program. No options \
available.")
args = parser.parse_args()
dyn = moleculardynamics()
dyn.read_input()
dyn.verlet()