#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
#* Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
#*
#* Copyright (C) 2017, 2018, 2019
#* Ralph Welsch, DESY, <ralph.welsch@desy.de>
#*
#* Copyright (C) 2020, 2021, 2022, 2023
#* Ludger Inhester, DESY, ludger.inhester@cfel.de
#*
#* This file is part of CDTK.
#*
#* CDTK is free software: you can redistribute it and/or modify
#* it under the terms of the GNU General Public License as published by
#* the Free Software Foundation, either version 3 of the License, or
#* (at your option) any later version.
#*
#* This program is distributed in the hope that it will be useful,
#* but WITHOUT ANY WARRANTY; without even the implied warranty of
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#* GNU General Public License for more details.
#*
#* You should have received a copy of the GNU General Public License
#* along with this program. If not, see <http://www.gnu.org/licenses/>.
#*
#* **************************************************************************
"""
Collection of tools used in the Newton package
Note on the implementation:
positions, velocities, etc., are expected as N x 3 arrays
[ [x1, y1, z1],
[x2, y2, z2],
...
]
masses are expected as an N array
[ m1, m2, .... ]
"""
import math
import numpy as np
from functools import reduce
[docs]
def calcCM(positions,masses,**opts):
ndim = opts.get('ndim',3) # spatial dimension
cm = np.zeros(ndim,float)
for i in range(len(masses)):
cm = cm + (positions[i]*masses[i])/masses.sum()
return cm
[docs]
def etrans(vel,masses,**opts):
ndim = opts.get('ndim',3) # spatial dimension
ek = 0.0
vcm = np.zeros(ndim,float)
mtot = masses.sum()
for i in range(len(masses)):
vcm = vcm + (vel[i]*masses[i])
vcm = vcm/mtot
ek = 0.5*mtot*np.dot(vcm,vcm)
return ek
[docs]
def erot(pos,vel,masses):
cm = calcCM(pos,masses)
vcm = calcCM(vel,masses)
poscm = pos - cm
velcm = vel - vcm
J = 0.0
for i in range(len(masses)):
J = J + masses[i]*np.cross(poscm[i],velcm[i])
ii,ix = principalAxes(poscm,masses)
ii = ii + 1.e-8 # regularization
ii = np.diag(ii**-1)
iii = np.dot(ix,np.dot(ii,np.transpose(ix)))
return 0.5*np.dot(J,np.dot(iii,J))
[docs]
def evib(pos,vel,masses,**opts):
ndim = opts.get('ndim',3) # spatial dimension
ekintot = 0.0
for i in range(len(masses)):
ekintot = ekintot + 0.5*masses[i]*np.dot(vel[i],vel[i])
return ekintot - erot(pos,vel,masses) - etrans(vel,masses,ndim=ndim)
[docs]
def pvib(pos,vel,masses,**opts):
ndim = opts.get('ndim',3) # spatial dimension
# [vel_vib] = [vel] - [vel_trans] - [vel_rot]
# compute translational velocity [vel_trans] for the total system
vel_trans = np.zeros(ndim,float)
mtot = masses.sum()
for i in range(len(masses)):
vel_trans = vel_trans + (vel[i]*masses[i])
vel_trans = vel_trans / mtot
if ndim == 3:
# compute rotational velocity [vel_rot] for each particle
cm = calcCM(pos,masses,ndim=ndim)
vcm = calcCM(vel,masses,ndim=ndim)
poscm = pos - cm
velcm = vel - vcm
J = 0.0
for i in range(len(masses)):
J = J + masses[i]*np.cross(poscm[i],velcm[i])
ii,ix = principalAxes(poscm,masses,ndim=ndim)
ii = ii + 1.e-8 # regularization
ii = np.diag(ii**-1)
iii = np.dot(ix,np.dot(ii,np.transpose(ix)))
omega = np.dot(iii,J) # angular velocity of the molecular frame
vel_rot = np.zeros((len(masses),ndim),float)
for i in range(len(masses)):
vel_rot[i] = np.cross(omega,poscm[i])
elif ndim == 1:
# one dimensional particles carry zero rotational energy
vel_rot = np.zeros((len(masses),ndim),float)
# compute the vibrational momentum for each particle
vel_vib = np.zeros((len(masses),ndim),float)
pv = []
for i in range(len(masses)):
vel_vib[i] = vel[i] - vel_trans - vel_rot[i]
pv.append(masses[i] * vel_vib[i])
return np.array(pv)
[docs]
def energy_partition(pos,vel,masses):
etr = etrans(vel,masses)
ero = erot(pos,vel,masses)
ekintot = 0.0
for i in range(len(masses)):
ekintot = ekintot + 0.5*masses[i]*np.dot(vel[i],vel[i])
evi = ekintot - etr - ero
return etr,ero,evi
[docs]
def jtot(pos,vel,masses):
cm = calcCM(pos,masses)
vcm = calcCM(vel,masses)
poscm = pos - cm
velcm = vel - vcm
J = 0.0
for i in range(len(masses)):
J = J + masses[i]*np.cross(poscm[i],velcm[i])
# v = waterAxes(poscm,masses)
# JR = np.dot(np.transpose(v),J)
# return J,JR,np.dot(JR,JR)**0.5
return J,np.dot(JR,JR)**0.5
[docs]
def principalAxes(q,masses,**opts):
ndim = opts.get('ndim',3) # spatial dimension
it = inertiaTensor(q,masses)
eval,evec = np.linalg.eig(it)
for i in range(ndim):
evec[:,i] = evec[:,i]/np.dot(evec[:,i],evec[:,i])**0.5
evec[:] = evec.take(eval.argsort(),axis=1) # axis=1: eigenvectors are expected in columns
eval[:] = eval.take(eval.argsort())
return eval,evec
[docs]
def inertiaTensor(q,m):
"""
Return a 3x3 matrix with the inertia tensor corresponding
to the configuration x and masses m
"""
qcm = q - calcCM(q,m)
iT = np.zeros((3,3),float)
for i in range(len(m)): # sum over all the atoms
iT[0,0] = iT[0,0] + m[i]*(qcm[i,1]**2 + qcm[i,2]**2)
iT[1,1] = iT[1,1] + m[i]*(qcm[i,0]**2 + qcm[i,2]**2)
iT[2,2] = iT[2,2] + m[i]*(qcm[i,0]**2 + qcm[i,1]**2)
iT[0,1] = iT[0,1] - m[i]*(qcm[i,0]*qcm[i,1])
iT[0,2] = iT[0,2] - m[i]*(qcm[i,0]*qcm[i,2])
iT[1,2] = iT[1,2] - m[i]*(qcm[i,1]*qcm[i,2])
iT[1,0] = iT[0,1]
iT[2,0] = iT[0,2]
iT[2,1] = iT[1,2]
return iT
[docs]
def veltoprincipal(v,x,m):
mtot = m.sum()
# Position CM of fragment
xcm = 0.0
for i in range(len(m)):
xcm = xcm + (x[i]*m[i])
xcm = xcm/mtot
xx = x - xcm
# Velocity CM of fragment
vcm = 0.0
for i in range(len(m)):
vcm = vcm + (v[i]*m[i])
vcm = vcm/mtot
vv = v - vcm
# Principal axes in the CM system
eva,eve = principalAxes(xx,m)
# Transform velocity to principal axes
vv = np.dot(np.transpose(eve),vv)
return vv
[docs]
def veltoCM(v,m):
vcm = np.zeros(len(m),float)
for i in range(len(m)):
vcm = vcm + v[i]*m[i]
vcm = vcm/m.sum()
return v - vcm
[docs]
def velautocor(v0,v):
aut = np.dot(v0.flatten(),v.flatten())
aut = aut/np.dot(v0.flatten(),v0.flatten())
return aut
[docs]
def internalcoords(x):
v1 = x[0] - x[1]
v2 = x[0] - x[2]
d1 = np.dot(v1,v1)**0.5
d2 = np.dot(v2,v2)**0.5
a = math.acos(np.dot(v1,v2)/d1/d2)
return 0.5*(d1+d2),0.5*(d1-d2),a*180.0/math.pi
[docs]
def flatten_list(list):
"""
Return flattened list
"""
flattened_list = reduce(lambda x,y: x+y, list)
return flattened_list