#* **************************************************************************
#*
#* 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 CDTK.Tools.Mathematics as mt
import CDTK.Tools.Geometry as geom
[docs]
def normal_modes_from_HessianCart(Hc,m,x0=None,proj=False,subgroups=None,**opts):
"""
Normal modes analysis from the Cartesian Hessian matrix
Hc -- Hessian matrix in Cartesian coordinates (in general non mass-weighted)
np.array, shape=(N,N) where N is the number of coordinates
m -- Mass of each coordinate, 3N array for molecules in 3D space
np.array, shape=(N,)
x0 -- Reference geometry on which the Hessian was calculated [optional]
proj -- Requires x0. If True, the overall translation and rotation of the molecule, as
well as the rotation of any given >subgroups are projected out
before normal mode analysis [optional]
subgroups -- Requires x0 and proj=True. May contain arrays of atomic indices that
each specify a rotating subgroup with respect to the atom index
positions given in m. The rotation of this subgroup with the smallest
moment of inertia is projected out before normal mode analysis
[optional]
output:
omega -- array with angular frequencies of each mode
L -- mass weighted normal modes (in columns)
"""
if proj:
# Get minimal displacement vectors from Geometry
geomX = geom.XYZGeometry()
geomX.setAtomPositions(x0.reshape(x0.size / 3,3))
geomX.setMasses(m[::3]) # Geometry requires one mass per xyz triple
geomX.centerGeometry()
# Project out the translations and rotations
P0 = geomX.cartesianProjector(subgroups=subgroups)
Hc = np.dot(P0, np.dot(Hc, P0))
M = np.diag( 1.0/np.sqrt(m) ) # Diagonal matrix with entries 1/sqrt(m_i)
# Construct mass-weighted Hessian
Hw = np.dot(M,np.dot(Hc,M))
# Diagonalize it and sort eigenvalues and eigenvectors by increasing eigenvalue
l,L = np.linalg.eig(Hw)
mt.sortEigValsVecs(l,L)
# Find negative eigenvalues
b = l < 0.0 # Bool array with True for negative entries
l[b] = np.abs(l[b])
omega = np.sqrt(l)
return omega,L
[docs]
def q_mw_to_x(m,L):
"""
Normal mode Mass weighted displacements to unweighted Cartesian
m -- array with mass of each coordinate
L -- mass weighted normal modes
Return:
2D array with normalized Cartesian normal modes in _rows_
This array can now be used to perform displacements in the unweighted
Cartesian space.
"""
# Conversion to Cartesian normalized modes of unprojected Hessian.
# This is a somewhat dirty procedure and ideally translations and rotations
# should be projected out first
modes = L.transpose()
mode_masses = np.zeros(len(modes))
# Unweight mass-weighted modes to Cartesian
modesX = np.dot(modes, np.diag(1.0/np.sqrt(m)))
# Normalize Cartesian modes; mode masses = square of normalization factor
for i,mode in enumerate(modesX):
nfactor = 1.0/np.linalg.norm(mode)
modesX[i] = mode*nfactor
mode_masses[i] = nfactor*nfactor
return modesX,mode_masses