# * **************************************************************************
# *
# * 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/>.
# *
# * **************************************************************************
# general imports
import sys
from math import sqrt
# from numpy
import numpy as np
from numpy import array, dot, zeros
from numpy import argmin, identity, outer
from numpy.linalg import solve, eigvals, eig
# from modules of crsurface package
from CDTK.Tools.NumericalDerivatives import gradient, hessian, _directionalDerivative
from CDTK.Tools.NumericalDerivatives import _secondDerivative, _crossedDerivative
from CDTK.Tools.Mathematics import gramSchmidt
from CDTK.Tools.Conversion import au2an, an2au
from CDTK.Tools import GaussianLog
############################################################################
[docs]
class Coordinate(object):
"""List of numbers representing a vector in Cartesian space
It can represent a geometry, 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.
- Atributes:
massWeighted -- Boolean, tells if the coordinates are in mass weighted cartesians
angstrom -- Boolean, tells if the coordinates are in angstroms
coordinates -- array, the coordinates
masses -- masses corresponding to each coordinate, in AMU
"""
def __init__(self,
massWeighted=False,
massWeightedAu=False,
angstrom=None,
masses=None,
coordinates=None,
symbols=None):
if type(coordinates) == type([]):
self.coordinates = array(coordinates)
else:
self.coordinates = coordinates
if type(masses) == type([]):
self.masses = array(masses)
else:
self.masses = masses
self.massWeighted = massWeighted
self.massWeightedAu = massWeightedAu
self.angstrom = angstrom
self.symbols = symbols
# Setters and getters of the object attributes
[docs]
def setMassWeighted(self, m):
"""Sets the massWeighted attribute
m -- False or True
"""
if type(m) != type(True):
sys.exit("Error in setMassWeighted, Boolean expected")
self.massWeighted = m
[docs]
def setMassWeightedAu(self, m):
"""Sets the massWeighted attribute
m -- False or True
"""
if type(m) != type(True):
sys.exit("Error in setMassWeightedAu, Boolean expected")
self.massWeightedAu = m
[docs]
def setAngstrom(self, m):
"""Sets the angstrom attribute
m -- False or True
"""
self.angstrom = m
[docs]
def setCoordinates(self, x):
"""Expects array o list. Stores it as the coordinates of mode
x -- array or list object
"""
if type(x) == type([]):
self.coordinates = array(x)
elif type(x) == type(array([])):
self.coordinates = x
else:
sys.exit("Error in setCoordinates, coordinates are not\
either of list or array type")
[docs]
def setMasses(self, x):
"""Expects array o list. Stores it as the masses of mode
x -- array or list object
"""
if type(x) == type([]):
self.masses = array(x)
elif type(x) == type(array([])):
self.masses = x
else:
sys.exit("Error in setMasses, masses are not\
either of list or array type")
[docs]
def setSymbols(self, s):
"""Set the list of atomic symbols
"""
if isinstance(s, list):
ss = []
for a in s:
ss.append(a)
self.symbols = ss
else:
print("in setSymbols, not a list as argument")
self.symbols = None
[docs]
def getCoordinates(self):
"""Return A COPY of the coordinates array"""
return array(self.coordinates)
[docs]
def getMasses(self):
"""Return A COPY of the masses array"""
return array(self.masses)
[docs]
def getSymbols(self):
"""Return A COPY of the list of atomic symbols
"""
if self.symbols is None:
return None
else:
s = []
for a in self.symbols:
s.append(a)
return s
# --------------------------------------------
[docs]
def massFromGaussian(self, f):
"""Get mass from a Gaussian log file. Store in self.mass
f -- full path to a gaussian log file
"""
self.masses = array(GaussianLog.getMasses(f))
[docs]
def symbolsFromGaussian(self, f):
"""Get the atomic symbols from a gaussian file, put them on list
"""
self.symbols = GaussianLog.getSymbols(f)
[docs]
def toMassWeighted(self, au=0):
"""
Transforms the coordinates to mass weighted
au = 0, use AMU as units of mass
au = 1, use au as units of mass, 1 AMU = 1822.89 au
"""
if not self.massWeighted:
if type(self.masses) == type(None):
sys.exit("Error in toMassWeighted, no masses present")
m = self.getMasses()
if au:
m = m * 1822.89
self.setMassWeightedAu(True)
x = self.getCoordinates() * m**0.5
self.setCoordinates(x)
self.setMassWeighted(True)
[docs]
def unweightCoordinates(self):
"""Transforms the coordinates from mass weighted to normal cartesian
"""
if self.massWeighted:
if type(self.masses) == type(None):
sys.exit("Error in unweightCoordinates, no masses present")
m = self.getMasses()
if self.massWeightedAu:
m = m * 1822.89
self.setMassWeightedAu(False)
x = self.getCoordinates() / m**0.5
self.setCoordinates(x)
self.massWeighted = False
[docs]
def toAngstroms(self):
"""Coordinates from au to angstroms
"""
if type(self.angstrom) == type(None):
sys.exit("Error in toAngstrom, don't know the type of coordinates")
if not self.angstrom:
self.coordinates = map(au2an, self.coordinates)
self.angstrom = True
[docs]
def toAu(self):
"""Coordinates from angstroms to au
"""
if type(self.angstrom) == type(None):
sys.exit("Error in toAu, type of coordinates unknown")
if self.angstrom:
self.coordinates = map(an2au, self.coordinates)
self.angstrom = False
[docs]
def module(self):
"""Return module of coordinates array"""
x = self.getCoordinates()
return dot(x, x)**0.5
# ---------------------------------------------------------------------------
############################################################################
[docs]
class NormalMode(Coordinate):
"""Stores a normal mode and provides methods to interact with it.
Uses numpy.array to store coordinates and other information.
"""
def __mul__(self, other):
"""Product of NormalMode with a number yelds an XYZGeometry
Product of two normal modes returns the projection
between them
"""
if isinstance(other, int):
other = float(other)
if isinstance(other, float):
if type(self.coordinates) != type(None):
c = self.coordinates*other
s = self.getSymbols()
x = XYZGeometry(coordinates=c,
angstrom=None,
massWeighted=self.massWeighted,
masses=self.masses,
symbols=s)
return x
if isinstance(other, NormalMode):
return dot(self.getCoordinates(), other.getCoordinates())
print("product could not be performed, check type of args")
return None
def __add__(self, other):
"""Add two normal modes
The two normal modes will keep normalized
"""
if (self.angstrom == other.angstrom) and \
(self.massWeighted == other.massWeighted):
x = self.getCoordinates() + other.getCoordinates()
m = NormalMode()
m.setCoordinates(x)
m.setAngstrom(self.angstrom)
m.setMasses(self.getMasses())
m.setMassWeighted(self.massWeighted)
m.setSymbols(self.getSymbols())
m.normalize()
return m
else:
print("These two modes cannot be added due to unit or weight diff")
return None
def __eq__(self, other):
if self is other:
return True
return False
[docs]
def getModeNumber(self):
"""Return mode number
"""
return self.modeNumber
[docs]
def setModeNumber(self, i):
"""Sets modenumber attribute"""
self.modeNumber = i
[docs]
def modeFromGaussian(self, file, i):
"""Get the coordinates of a mode from a gaussian log file
and stores them as an array them self.coordinates
file -- full path to the gaussian frequency file
i -- index of the normal mode to retrieve
"""
self.setCoordinates(array(GaussianLog.getNormalMode(file, i)))
self.massFromGaussian(file)
self.massWeighted = False
self.angstrom = None
self.modeNumber = i
self.symbolsFromGaussian(file)
self.frequency = GaussianLog.getNormalModeFreq(file, i)
[docs]
def modeReducedMass(self):
"""Return the mode reduced mass in AMU, the same as given
by Gaussian.
"""
c = self.getCoordinates()
m = self.getMasses()
c = c / (dot(c, c)**0.5)
if not self.massWeighted:
c = c * m**0.5
c = c / (dot(c, c)**0.5)
mvec = (c**2 / m)
mu = 1.0/mvec.sum()
return mu
[docs]
def toMassWeighted(self):
"""Transforms the coordinates to mass weighted.
Specializes by also normalizing the mode.
"""
Coordinate.toMassWeighted(self)
self.normalize()
[docs]
def unweightCoordinates(self):
"""Transforms the coordinates from mass weighted to normal cartesian
Specializes by also normalizing the mode.
"""
Coordinate.unweightCoordinates(self)
self.normalize()
[docs]
def normalize(self):
"""Normalizes the coordinates array to a length of 1"""
self.coordinates = self.coordinates/self.module()
# ---------------------------------------------------------------------------
############################################################################
[docs]
class XYZGeometry(Coordinate):
"""Stores a geometry.
Uses numpy.array
"""
def __add__(self, other):
if type(self) != type(other):
sys.exit("__add__ in XYZGeometry, operands are of different type")
a = []
for i in range(len(self.coordinates)):
a.append(self.coordinates[i] + other.coordinates[i])
a = array(a)
if type(self.angstrom) != type(None):
return XYZGeometry(coordinates=a,
massWeighted=self.massWeighted,
massWeightedAu=self.massWeightedAu,
angstrom=self.angstrom, masses=self.masses,
symbols=self.getSymbols())
else:
return XYZGeometry(coordinates=a,
massWeighted=self.massWeighted,
massWeightedAu=self.massWeightedAu,
angstrom=other.angstrom, masses=self.masses,
symbols=self.getSymbols())
[docs]
def geomFromGaussian(self, file):
"""Get the coordinates from a gaussian log file and stores
them in Coordinate.coordinates
file -- full path to the gaussian frequency file
"""
self.setCoordinates(array(GaussianLog.getCoordinates(file)))
self.massFromGaussian(file)
self.massWeighted = False
self.angstrom = True
self.symbolsFromGaussian(file)
[docs]
def energy(self, func):
"""Returns the energy corresponding to the current coordinates
Is resposibility of the caller to ensure that the XYZGeometry
object has the coordinates in the correct units before calling
the getEnergy method.
func -- function to be used in the evaluation of the energy
"""
if type(self.coordinates) == type(None):
sys.exit("Error in getEnergy, no coordinates available")
return func(self.coordinates)
[docs]
def gradient(self, func, modes, step=0.1):
"""Return a list with the numerical derivatives
in the directions of the modes stored in the list modes
func -- function returning the energy
modes -- list of NormalMode specifying the directions to follow
step -- step for the numerical derivative
"""
g = []
for m in modes:
g.append(self.gradientElement(func, m, step))
return g
[docs]
def hessian(self, func, modes, step=0.1):
"""Return a list of lists with the crossed derivatives
with respect to all the modes in the list modes
func -- function returning the energy
modes -- list of NormalMode specifying the directions to follow
step -- step for the numerical derivative
"""
h = [[0.0 for i in range(len(modes))] for j in range(len(modes))]
for i in range(len(modes)):
for j in range(i, len(modes)):
h[i][j] = self.hessianElement(
func, modes[i], modes[j], step)
h[j][i] = h[i][j]
return h
[docs]
def cartMassWeightHessian(self, func, step=0.1, au=0):
"""Return a matrix (array type) with the Hessian in mw coords
func -- returns the potential
step -- step to take in the derivatives (bohr)
au = 0, hessian is weighted in AMU (default)
au = 1, hessian is weighted in au, 1 AMU = 1822.89 au
"""
x = self.getCoordinates()
d = identity(len(x), float)
h = hessian(func, x, d, step)
m = self.getMasses()
if au:
m = m * 1822.89
for i in range(len(h)):
for j in range(len(h)):
h[i, j] = h[i, j]/(sqrt(m[i]*m[j]))
return h
[docs]
def cartHessian(self, func, step=0.01):
"""Return a matrix (array type) with the Hessian in mw coords
"""
x = self.getCoordinates()
d = identity(len(x), float)
h = hessian(func, x, d, step)
return h
[docs]
def gradientElement(self, func, mode, step=0.1):
"""Return the derivative in the direction of the given
normal mode
func -- function returning the energy
mode -- NormalMode specifying the direction to follow
step -- step for the numerical derivative
"""
x = self.getCoordinates()
xm = mode.getCoordinates()
return directionalDerivative(func, x, xm, step)
[docs]
def hessianElement(self, func, mode1, mode2, step=0.1):
"""Return the second crossed derivative of a function with
respect to two modes (Hessian matrix element).
func -- function returning the energy
mode1 -- NormalMode specifying the direction to follow
mode2 -- NormalMode specifying the direction to follow
step -- step for the numerical derivative
"""
x = self.getCoordinates()
m1 = mode1.getCoordinates()
if mode1 == mode2:
return secondDerivative(func, x, m1, step)
else:
m2 = mode2.getCoordinates()
return crossedDerivative(func, x, m1, m2, step)
[docs]
def frequencyFromMode(self, func, mode, step=0.1):
"""Frequency in cm-1 for a given mode at the current geometry
remember, this quantity is only menaniful in a stationary point
func -- function returning the energy
mode -- NormalMode specifying the direction to follow
step -- step for the numerical derivative
"""
m = mode.getCoordinates()
x = self.getCoordinates()
d = secondDerivative(func, x, m, step)
mu = mode.modeReducedMass()
f = aunits2wavenumber(d, mu)
return f
[docs]
def projectionToModeS(self, xref, modes):
"""Return an array with the projection of the current
geometry into the given normal modes
xref -- XYZGeometry object or array with the reference geometry
modes -- list of the normal modes to be projected
"""
if isinstance(xref, XYZGeometry):
xr = xref.getCoordinates()
else:
xr = array(xref)
if isinstance(modes, NormalMode):
mlist = []
mlist.append(modes)
else:
mlist = modes
d = self.getCoordinates() - xr
q = []
for mod in mlist:
qq = dot(mod.getCoordinates(), d)
q.append(qq)
return q
[docs]
def centerMolecule(self):
"""Center 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
"""
x = 0.0
y = 0.0
z = 0.0
mtot = 0.0
c = self.getCoordinates()
m = self.getMasses()
nat = len(m)/3
for i in range(nat):
x = x + c[3*i+0]*m[3*i+0]
y = y + c[3*i+1]*m[3*i+1]
z = z + c[3*i+2]*m[3*i+2]
mtot = mtot + m[3*i+2]
x = x/mtot
y = y/mtot
z = z/mtot
cm = array([x, y, z])
return cm
[docs]
def translate(self, vec):
"""Translate the system according to the given 3D vector
"""
x = self.getCoordinates()
for i in range(len(x)/3):
x[3*i+0] = x[3*i+0] + vec[0]
x[3*i+1] = x[3*i+1] + vec[1]
x[3*i+2] = x[3*i+2] + vec[2]
self.setCoordinates(x)
[docs]
def coordinatesToFile(self, f):
"""Write the coordinates to a file in XYZ format
The coordinates are written always in angstroms
"""
out = file(f, 'w')
coords = self.getCoordinates()
s = self.getSymbols()
if not self.angstrom:
coords = map(au2an, coords)
nat = len(coords)/3
for i in range(nat):
string = "%s %7.3f %7.3f %7.3f\n" % \
(s[i], coords[3*i], coords[3*i+1], coords[3*i+2])
out.write(string)
out.close()
[docs]
def projectorTransRot(self):
"""Return a projector to eliminate rotations and translations
"""
ptrans = self.projectorTrans()
prot = self.projectorRot()
ptot = ptrans + prot
return ptot
[docs]
def projectorTrans(self):
"""Return a projector to eliminate translations
"""
m = self.getMasses()
trx = zeros(len(m), float)
tri = zeros(len(m), float)
trz = zeros(len(m), float)
for i in range(len(m)/3):
trx[3*i+0] = sqrt(m[3*i+0])
tri[3*i+1] = sqrt(m[3*i+1])
trz[3*i+2] = sqrt(m[3*i+2])
ovecs = gramSchmidt([trx, tri, trz])
px = Projector(ovecs[0])
py = Projector(ovecs[1])
pz = Projector(ovecs[2])
ptrans = px + py + pz
return ptrans, ovecs
[docs]
def projectorRot(self):
"""Return a projector to eliminate rotations
"""
itv = self.momentsOfInertia()
x = self.getCoordinates()
cm = self.centerOfMass()
xc = zeros(len(x), float)
# set atom coordinates with respect to the center o mass
for i in range(len(x)/3):
xc[3*i+0] = x[3*i+0] - cm[0]
xc[3*i+1] = x[3*i+1] - cm[1]
xc[3*i+2] = x[3*i+2] - cm[2]
m = self.getMasses()
rot1 = zeros(len(m), float)
rot2 = zeros(len(m), float)
rot3 = zeros(len(m), float)
# Generate the vectors describing the (infinitesimal) rotations
for i in range(len(m)/3):
for j in range(3):
rot1[3*i+j] = ((dot(xc[3*i:3*i+3], itv[:, 1]) * itv[j, 2]
- dot(xc[3*i:3*i+3], itv[:, 2]) * itv[j, 1]) /
sqrt(m[i]))
rot2[3*i+j] = ((dot(xc[3*i:3*i+3], itv[:, 2]) * itv[j, 0]
- dot(xc[3*i:3*i+3], itv[:, 0]) * itv[j, 2]) /
sqrt(m[i]))
rot3[3*i+j] = ((dot(xc[3*i:3*i+3], itv[:, 0]) * itv[j, 1]
- dot(xc[3*i:3*i+3], itv[:, 1]) * itv[j, 0]) /
sqrt(m[i]))
ovecs = gramSchmidt([rot1, rot2, rot3])
prot = Projector(ovecs[0]) + Projector(ovecs[1]) + Projector(ovecs[2])
return prot, ovecs
[docs]
def momentsOfInertia(self):
"""Return a 3x3 matrix where the vectors v[:,i] are
the moments of inertia that diagonalize the inertia tensor
"""
x = self.getCoordinates()
m = self.getMasses()
t = inertiaTensor(x, m)
eval, evec = eig(t)
return eval, evec
############################################################################
[docs]
class Projector(object):
"""Define a projector, basically a matrix to be applied on the form
D=P.M.P where P is the projector and M a matrix
"""
def __init__(self, x=None):
if type(x) == type(None):
self.mat = None
elif type(x) == type([]):
self.fromArray(x)
elif type(x) == type(array([])):
if len(x.shape) == 1:
self.fromArray(x)
elif len(x.shape) == 2:
self.mat = x
else:
self.mat = None
else:
self.mat = None
def __add__(self, other):
"""Add two projectors by simply adding their matrices
"""
m = self.getMat() + other.getMat()
p = Projector(m)
return p
def __mul__(self, other):
"""Apply a projector to a matrix, so make M2 = P.M1.P
"""
if type(other) == type(array([])) and len(other) == len(self.getMat()):
mi = dot(other, self.getMat())
mf = dot(self.getMat(), mi)
return mf
else:
print("matrix and projector not compatible, check ranks")
[docs]
def setMat(self, m):
"""Set the atribute self.mat with the given matrix m"""
if type(m) == type(array([])) and len(m.shape) == 2:
self.mat = m
[docs]
def getMat(self):
"""Return a copy of self.mat"""
m = array(self.mat)
return m
[docs]
def fromArray(self, x):
"""Set the projector corresponding to an array as P=V.V+
where V is a column vector and V+ is a row vector
"""
xp = array(x)
mat = outer(xp, xp)
self.setMat(mat)
[docs]
def conjugateSpace(self):
"""Return a projector on the conjugate space of self
Return 1-P where P is the projector in self and 1 is the identity
"""
rank = len(self.getMat())
i = identity(rank, float)
m = i - self.getMat()
p = Projector(m)
return p
# ---------------------------------------------------------------------------
############################################################################
# MODULE METHODS
[docs]
def rangeOfModesFromGaussian(f, i, j):
"""Return a list of normalMode objects
f -- full path to the file containing the modes
i -- integer, first mode to retrieve
j -- integer, last mode to retrieve
Beware of the index conventions in Python, the last
index in a range(i,j) is not considered part of the sequence.
This is not so in this function.
"""
l = []
for k in range(i, j+1):
m = NormalMode()
m.modeFromGaussian(f, k)
l.append(m)
return l
[docs]
def rangeOfModesFromMatrix(h, m=None, freqs=0, omegas=0, mw=0, minnum=1):
"""From a matrix diagonalize and get the normal modes
Assume that h is a mass weighted hessian, diagonalize it
and return the corresponding normal modes as normalized
cartesian displacements (the same thing as gaussian does)
h -- hessian matrix
m -- masses of the system, necessary to transform from normal modes to cartesian displacements
freqs = 0, do not return the list of frequencies
freqs = 1, return as the second argument the list of frequencies in cm-1
omegas = 0, default
omegas = 1, return a list with the square roots of the eigenvalues
minnum = int, the mode from which ascendent numbering is assigned,
minmum = 7 means that 6 modes are not numbered
mw = 1, the normal modes are returned in mass-weighted coordinates
"""
if m is None:
m = np.ones(len(h), float)
h = array(h)
evals, evecs = eig(h)
# some extrange problem with projected hessians makes
# appear complex numbers in evals and evecs, without
# affecting the results, here we solve it "patatoidalment"
evals = evals.real
evecs = evecs.real
evecs = np.transpose(evecs)
evecs = np.take(evecs, np.argsort(abs(evals)))
evals = np.take(evals, np.argsort(abs(evals)))
nmodes = []
i = 1
ii = 1
olist = []
for v in evecs:
vv = NormalMode()
vv.setCoordinates(v)
vv.setMasses(m)
vv.setMassWeighted(True)
if not mw:
vv.unweightCoordinates()
vv.normalize()
if ii >= minnum:
vv.setModeNumber(i)
i = i+1
nmodes.append(vv)
olist.append(abs(evals[ii-1])**0.5)
ii = ii+1
if freqs:
flist = []
i = 0
for mod in nmodes:
f = aunits2wavenumber(evals[i])
flist.append(f)
i = i+1
return nmodes, flist
elif omegas:
return nmodes, olist
else:
return nmodes
[docs]
def modesOrthonormalization(modes):
"""Given a list of normal mode objects, orthonormalize them
"""
vectors = []
for m in modes:
vectors.append(m.getCoordinates())
onvects = gramSchmidt(vectors)
for i in range(len(modes)):
modes[i].setCoordinates(onvects[i])
[docs]
def inertiaTensor(x, m):
"""Return a 3x3 matrix with the inertia tensor corresponding
to the configuration x and masses m
"""
iT = zeros((3, 3), float)
for i in range(len(x)/3): # sum over all the atoms
ma = m[3*i]
xa = x[3*i+0]
ya = x[3*i+1]
za = x[3*i+2]
iT[0, 0] = iT[0, 0] + ma * (ya**2 + za**2)
iT[1, 1] = iT[1, 1] + ma * (xa**2 + za**2)
iT[2, 2] = iT[2, 2] + ma * (xa**2 + ya**2)
iT[0, 1] = iT[0, 1] - ma * (xa*ya)
iT[0, 2] = iT[0, 2] - ma * (xa*za)
iT[1, 2] = iT[1, 2] - ma * (ya*za)
iT[1, 0] = iT[0, 1]
iT[2, 0] = iT[0, 2]
iT[2, 1] = iT[1, 2]
return iT
[docs]
def aunits2wavenumber(d, mu=1.0):
"""Transform a second derivative in atomic units to the
corresponding vibrational frequency in cm-1
"""
img = 0
if d < 0.0:
d = d*(-1.0)
img = 1
f = (d/mu) * 9.37905*10**29 # hart/(bohr**2*amu) -> SI
f = f**0.5 # to Hz
f = f/(200.0 * 3.141592 * 299792458.0) # to cm-1
if img:
f = f*(-1.0)
return f
[docs]
def centerOfMass(xs, m):
"""Return 3D array with the center of mass
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
note: units are arbitrary
"""
xc = zeros(3, float)
mtot = 0.0
for i in range(len(xs)):
xc = xc + xs[i]*m[i]
mtot = mtot + m[i]
xc = xc/mtot
return xc
[docs]
def openOutputFiles(name, modes, access='w'):
"""
Return 3 lists with the energy, gradient, and hessian file objects
name -- basename to use to name the files
modes -- modes with respect to which gradients and hessians will be evaluated
"""
ener = []
grad = []
hess = []
enername = name + "_E"
gradname = name + "_G"
hessname = name + "_H"
ener.append(file(enername, access, 1))
for i in range(len(modes)):
try:
ki = modes[i].getModeNumber()
namei = gradname + "_" + str(ki)
except AttributeError:
ki = i+1
namei = gradname + "_n" + str(ki)
grad.append(file(namei, access, 1))
hf = []
for j in range(i+1):
try:
kj = modes[j].getModeNumber()
namej = hessname + "_" + str(ki) + "_" + str(kj)
except AttributeError:
kj = j+1
namej = hessname + "_n" + str(ki) + "_n" + str(kj)
hf.append(file(namej, access, 1))
hess.append(hf)
return ener, grad, hess
[docs]
def openOutputFilesProd(name, modes, access='w'):
"""Return 4 lists with the energy, hessian, hessian*q and
q*hessian*q file objects
name -- basename to use to name the files
modes -- modes with respect to which gradients and hessians will be evaluated
"""
ener = []
hess = []
hessq = []
hessqq = []
enername = name + "_E"
hessname = name + "_H"
hessnameq = name + "_Hq"
hessnameqq = name + "_Hqq"
ener.append(file(enername, access, 1))
hessqq.append(file(hessnameqq, access, 1))
for i in range(len(modes)):
try:
ki = modes[i].getModeNumber()
namei = hessnameq + "_" + str(ki)
except AttributeError:
ki = i+1
namei = hessnameq + "_n" + str(ki)
hessq.append(file(namei, access, 1))
hf = []
for j in range(i+1):
try:
kj = modes[j].getModeNumber()
namej = hessname + "_" + str(ki) + "_" + str(kj)
except AttributeError:
kj = j+1
namej = hessname + "_n" + str(ki) + "_n" + str(kj)
hf.append(file(namej, access, 1))
hess.append(hf)
return ener, hess, hessq, hessqq
[docs]
def openOutputFileList1(name, modes, access='w'):
"""Return list of files, one for each mode, ex. to store gradients
name -- basename to use to name the files
modes -- modes with respect to which gradients and hessians will be evaluated
"""
grad = []
gradname = name
for i in range(len(modes)):
try:
ki = modes[i].getModeNumber()
namei = gradname + "_" + str(ki)
except AttributeError:
ki = i+1
namei = gradname + "_n" + str(ki)
grad.append(file(namei, access, 1))
return grad
[docs]
def openOutputFileList0(name, modes, access='w'):
"""Return a file
name -- basename to use to name the files
modes -- modes with respect to which gradients and hessians will be evaluated
"""
ener = []
enername = name
ener.append(file(enername, access, 1))
return ener
[docs]
def outputEnerToFiles(coords, ener, file):
"""Prints energy to file"""
if isinstance(coords, float):
coords = [coords]
outstr = "%13.8f" % (ener)
for c in coords:
outstr = outstr + " %8.4f" % (c)
outstr = outstr + "\n"
if isinstance(file, list):
file[0].write(outstr)
else:
file.write(outstr)
[docs]
def outputGradToFiles(coords, grad, file):
"""Print derivative with respect to all modes to the corresponding files
"""
if isinstance(coords, float):
coords = [coords]
for i in range(len(grad)):
outstr = "%13.8f" % (grad[i])
for c in coords:
outstr = outstr + " %8.4f" % (c)
outstr = outstr + "\n"
file[i].write(outstr)
[docs]
def outputHessToFiles(coords, hess, file):
"""Print crossed derivative with respect to each pair of
mode to files
"""
if isinstance(coords, float):
coords = [coords]
for i in range(len(hess)):
for j in range(i+1):
outstr = "%13.8f" % (hess[i][j])
for c in coords:
outstr = outstr + " %8.4f" % (c)
outstr = outstr + "\n"
file[i][j].write(outstr)
[docs]
def stringToFiles(string, filelist):
"""Write string to all the files found in filelist
filelist can be a single file or a nested structure of lists
of lists of files down to any depth
"""
if not isinstance(filelist, list):
filelist.write(string)
else:
for l in filelist:
stringToFiles(string, l)
[docs]
def closeListOfFiles(filelist):
"""Closes all the files in a list of file objects
filelist can be a single file or a nested structure of lists
of lists of files down to any depth
"""
if not isinstance(filelist, list):
filelist.close()
else:
for l in filelist:
closeListOfFiles(l)
def _expectedDE(f, x, dir):
"""Debug routine
x -- cartesian coordinates
dir -- list of vectors (normal modes)
returns the expected energy change from position
x to some minimum energy point x0 assuming a quadratic
expansion in the directions contained in "dir"
"""
g0 = gradient(f, x, dir)
h0 = hessian(f, x, dir)
g0 = -1*g0
dq = solve(h0, g0)
# hi = inv(h0)
# dq = dot(hi,g0)
displ = sqrt(dot(dq, dq))
dE = -0.5 * dot(dq, dot(h0, dq))
fE = energy(f, x) + dE
ev = eigvals(h0)
minh = ev[argmin(ev)]
return dE, fE, displ, minh