#* **************************************************************************
#*
#* 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/>.
#*
#* **************************************************************************
import numpy as np
import itertools as it
from .Mathematics import eulerMatrix
import CDTK.Tools.Mathematics as mt
[docs]
class Coordinate(object):
"""
List of numbers representing a coordinate vector
It can represent a normal mode, or in general anything
that can be represented by a set of coordinates. Can also contain
the masses related to those coordinates.
Uses numpy.array to store coordinates and other information.
coordinates -- array, the coordinates
masses -- masses corresponding to each coordinate, in AMU
labels -- list of labels representing each coordinate
"""
def __init__(self,**opts):
self._coordinates = None
self._masses = None
self._labels = None
pass
# Setters and getters of the object attributes
[docs]
def setCoordinates(self,x):
self._coordinates=np.asarray(x)
[docs]
def getCoordinates(self):
return self._coordinates
coordinates = property(getCoordinates,setCoordinates)
[docs]
def setMasses(self,x):
self._masses=np.asarray(x)
[docs]
def getMasses(self):
return self._masses
masses = property(getMasses,setMasses)
[docs]
def setLabels(self,s):
self._labels=s
[docs]
def getLabels(self):
return self._labels
labels = property(getLabels,setLabels)
[docs]
class XYZGeometry(object):
"""
Stores a cartesian geometry.
"""
def __init__(self,**opts):
self._atomPositions = None
self._masses = None
self._labels = None
self._pbcDims = None
[docs]
def setAtomPositions(self,x):
self._atomPositions=np.asarray(x)
[docs]
def getAtomPositions(self):
return self._atomPositions
atomPositions = property(getAtomPositions,setAtomPositions)
[docs]
def setMasses(self,x):
self._masses=np.asarray(x)
[docs]
def getMasses(self):
return self._masses
masses = property(getMasses,setMasses)
[docs]
def setLabels(self,s):
self._labels=s
[docs]
def getLabels(self):
return self._labels
labels = property(getLabels,setLabels)
[docs]
def setPbcDims(self,x):
self._pbcDims=np.asarray(x)
[docs]
def getPbcDims(self):
return self._pbcDims
pbcDims = property(getPbcDims,setPbcDims)
[docs]
def centerGeometry(self):
"""
Move the center of mass of the molecule to the origin
"""
cm = self.centerOfMass()
self.translate(-1.0*cm)
[docs]
def centerOfMass(self):
"""
Return a 3 components vector with the center of mass of the molecule
"""
return centerOfMass(self.atomPositions,self.masses)
[docs]
def translate(self,translation):
"""
Translate the system according translation
translation -- 3D vector
"""
self.atomPositions = self.atomPositions + translation
[docs]
def momentsOfInertia(self, sort=True, subgroup=None):
"""
Return a 3x3 matrix where the vectors v[:,i] are
the moments of inertia that diagonalize the inertia tensor
sort - if True (default), return sorted eigenvalues and eigenvectors
subgroup - index array, if it is given, the moments of inertia for a
rotation of the atoms specified by the indices is returned.
Default: None
"""
x = self.atomPositions
m = self.masses
if subgroup is not None:
t = inertiaTensor(x[subgroup], m[subgroup])
else:
t = inertiaTensor(x,m)
eval,evec = np.linalg.eig(t)
if sort:
s = eval.argsort()
evec = evec[:,s]
eval = eval[s]
return eval,evec
[docs]
def cartesianTranslationVector(self, ix):
"""
Returns the vector for displacement along axis ix in unweighted
cartesian coordinates
"""
self.centerGeometry()
x = self.atomPositions
m = self.masses
d = np.zeros_like(x)
lambda_t, S_t = self.momentsOfInertia(sort=False)
for i in range(d.shape[0]):
d[i,ix] = m[i]
return d
[docs]
def cartesianRotationVector(self, ix, subgroup=None):
"""
Returns the vector for minimal rotation around axis ix in unweighted
cartesian coordinates
ix - index of the vector
subgroup - index array of atoms forming a subgroup. If None (default),
the whole molecule is rotated. For subgroup rotations, regardless of
the provided index, the rotation with the smallest moment of inertia
is returned
"""
self.centerGeometry()
x = self.atomPositions
m = self.masses
d = np.zeros_like(x)
if subgroup is not None:
lambda_t, S_t = self.momentsOfInertia(sort=False, subgroup=subgroup)
ix = np.argmin(lambda_t)
for i in subgroup:
for j in range(3):
for k in range(3):
d[i,j] = d[i,j] + m[i] * mt.leviCivita(ix,k,j) * np.dot(S_t[k], x[i])
else:
lambda_t, S_t = self.momentsOfInertia(sort=False)
for i in range(d.shape[0]):
for j in range(3):
for k in range(3):
d[i,j] = d[i,j] + m[i] * mt.leviCivita(ix,k,j) * np.dot(S_t[k], x[i])
return d
[docs]
def cartesianVectors(self, subgroups=None):
"""
Returns a (3N,6+s) matrix where s is the number of subgroups.
The matrix contains, columnwise, orthonormal (3N) vectors. These vectors
form a basis of the space spanned by the vectors corresponding to
translation, overall rotation, and subgroup rotation around the axis
with the smallest moment of inertia of the respective group. They can
be used to project out normal modes corresponding to these movements
from a Hessian matrix. The vectors are given in unweighted, cartesian
coordinates.
subgroups - array with arrays of subgroup indices
"""
d = self.cartesianTranslationVector(0).flatten()
d = np.c_[d, self.cartesianTranslationVector(1).flatten()]
d = np.c_[d, self.cartesianTranslationVector(2).flatten()]
d = np.c_[d, self.cartesianRotationVector(0).flatten()]
d = np.c_[d, self.cartesianRotationVector(1).flatten()]
d = np.c_[d, self.cartesianRotationVector(2).flatten()]
if subgroups is not None:
for sg in subgroups:
d = np.c_[d, self.cartesianRotationVector(0,subgroup=sg).flatten()]
d = np.linalg.qr(d)[0] # Gram-Schmidt orthogonalization
return d
[docs]
def cartesianProjector(self, subgroups=None):
"""
Returns a (3N, 3N) matrix that projects out the subspace spanned by the
vectors corresponding to translation, overall rotation, and subgroup
rotation around the axis with the smallest moment of inertia. Can be
used to project Cartesian Hessian matrices prior to normal mode analysis.
subgroups -- array with arrays of subgroup indices
"""
D0 = self.cartesianVectors(subgroups=subgroups)
return np.identity(self.getAtomPositions().size) - np.dot(D0, D0.T)
[docs]
def radiusGyration(self):
"""
Returns the radius of gyration defined as
rg**2 = sum_i m_i * (r_i - r_COM)**2 / M
"""
com = self.centerOfMass()
mtot = np.sum(self.masses)
rg2 = np.dot(self.masses, np.linalg.norm( self.atomPositions - com, axis=1)**2) / mtot
return np.sqrt(rg2)
############################################################################
# MODULE METHODS
[docs]
def inertiaTensor(geom,m):
"""
Return a 3x3 matrix with the inertia tensor corresponding
to the configuration x and masses m
"""
iT = np.zeros((3,3),float)
for i,atom in enumerate(geom):
iT[0,0] = iT[0,0] + m[i]*(atom[1]**2 + atom[2]**2)
iT[1,1] = iT[1,1] + m[i]*(atom[0]**2 + atom[2]**2)
iT[2,2] = iT[2,2] + m[i]*(atom[0]**2 + atom[1]**2)
iT[0,1] = iT[0,1] - m[i]*(atom[0]*atom[1])
iT[0,2] = iT[0,2] - m[i]*(atom[0]*atom[2])
iT[1,2] = iT[1,2] - m[i]*(atom[1]*atom[2])
iT[1,0] = iT[0,1]
iT[2,0] = iT[0,2]
iT[2,1] = iT[1,2]
return iT
[docs]
def centerOfMass(geom,m):
"""
Return 3D array with the center of mass of geom
xs -- list or array of 3D arrays with the positions of the atoms
m -- list or array with the mass of each atom in xs
"""
xc=np.zeros(3,float)
mtot = m.sum()
for i,x in enumerate(geom):
xc = xc + x*m[i]
xc = xc/mtot
return xc
[docs]
def eckartFrame(x,xref,m):
"""
Calculate the euler angles that rotate x to the Eckart frame w.r.t xref
x,xref -- Nx3 array with atom coordinates
m -- array with masses, len(m) = N
"""
import scipy.optimize
# x and xref are copied and the COM of the copies is moved to the origin
xx = np.array(x)
mm = np.asarray(m)
cm = centerOfMass(xx,mm)
xx = xx-cm
xxref = np.array(xref)
cm = centerOfMass(xxref,mm)
xxref = xxref-cm
guess=np.zeros(3,float)
euler = scipy.optimize.fmin_bfgs(angularMomentumEckart,
guess,
args=(xx,xxref,mm),
gtol=1.0e-3,
maxiter=100,
epsilon=1.0e-6,
disp=0 # set to 1 for convergence info
)
return rotateGeometry(xx,euler)
[docs]
def angularMomentumEckart(euler,*args):
"""
Calculate the Eckart conditions given 3 euler angles
This is the functions to be zeroed (minimized) if x is to be
roted to the Eckart frame defined by xref
euler -- euler angles a,b,c
args[0] -- structure; (Nx3) array
args[1] -- reference structure; (Nx3) array
args[2] -- mass vector; N array
Returns the square of the norm of the vector:
amv = sum_i { mass[i] * (dx[i] *cross* xref[i]) }
"""
xx = np.array(args[0])
xxref= np.array(args[1])
mass=args[2]
xx = rotateGeometry(xx,euler)
dx = xx-xxref
amv = np.zeros(3,float)
for i,d in enumerate(dx):
amv = amv + np.cross(d,xxref[i])*mass[i]
am = np.dot(amv,amv)
return am
[docs]
def rotateGeometry(x,euler):
xx=np.array(x)
rotMat = eulerMatrix(euler[0],euler[1],euler[2])
for i in range(len(xx)):
xx[i] = np.dot(rotMat,xx[i])
return xx
[docs]
def pbc_translation_vectors(pbc_dims):
'''
Periodic boundary condition: calculate translation vectors from box dimensions
Parameters:
pbc_dims -- array containing the box dimensions along x/y/z
Returns:
[3,3] matrix containing the translation vectors along x/y/z
'''
return np.identity(3) * np.asarray(pbc_dims)
[docs]
def pbc_translation_matrix(pbc_dims):
'''
Periodic boundary condition: calculate translation matrix from box dimensions
Minimum image convention
Parameters:
pbc_dims -- array containing the box dimensions along x/y/z
Returns:
[3,3,3,3] matrix containing the translation vectors and their linear combinations
e.g. T_xyz_3D[0,2,1,:] is the vector that generates a translation in -x,+y,0z direction
'''
T_xyz = pbc_translation_vectors(pbc_dims)
T_xyz_3D = np.zeros([3,3,3,3])
factor = np.array([-1,0,1]) # subtract, ignore, add this translation vector
for ix in it.product(range(3),repeat=3):
T_xyz_3D[ix] = factor[ix[0]] * T_xyz[0] + factor[ix[1]] * T_xyz[1] \
+ factor[ix[2]] * T_xyz[2]
return T_xyz_3D