#* **************************************************************************
#*
#* 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/>.
#*
#* **************************************************************************
"""
SECTIONS
- LINEAR ALGEBRA
- TENSOR ALGEBRA
- VECTORS in 3D SPACE
- INTERPOLATION
- NUMERICAL
- MONTE-CARLO
- OTHER
- internal functions
"""
import os
import sys
import random
import math
from math import sqrt,pi,sin,cos,asin,acos,exp
import scipy.integrate as si
import scipy.interpolate as sip
from scipy.special import wofz
import numpy as np
import CDTK
from CDTK.Tools import Conversion as cv
from CDTK.Tools.Iterators import MultiCounter
from CDTK.Tools.Utils import ReverseArray, ShiftArray
##############################################################################
# --- LINEAR ALGEBRA: Matrix and vector related utilities ---
##############################################################################
[docs]
def krylov(func,v,**opts):
"""Construct a Krylov space
func -- takes one vector as argument and returns the action of
the Hamiltonian matrix on that vector
v -- Start vector for the construction of the Krylov space
argdict -- dictionary:
order: order (dimension) of the Krylov space to be constructed
kvecs : np.array, space to store the krylov vectors
kmat(order x order) : space to store the krylov matrix
work : scratch array of length len(v)
Return the tridiagonal matrix, kmat, and a list of the vectors
spanning the Krylov space kvecs
Test: The condition K^* A K = t has been tested for this routine
using an Hermitian matrix 40x40 and a Krylov space 5x5
Orthonormality of the Kvecs was also tested
19.01.2007
"""
n = len(v)
order = opts.get("order",10)
kvecs = opts.get("kvecs",None)
kmat = opts.get("kmat",None)
t = opts.get("work",None)
if not kvecs:
kvecs = np.zeros((order,n),type(v[0]))
if not kmat:
kmat = np.zeros((order,order),type(v[0]))
if not t:
t = np.zeros(n,type(v[0]))
# Normalize starting vector
kvecs[0] = v/(np.dot(v.conjugate(),v)**0.5)
a = 0.0
b = 0.0
for j in range(order):
# Calculate H(v) on the current vector
t[:] = func(kvecs[j])
if j > 0:
# Schmidt orthogonalize with respect to prev. prev. vector
t[:] = t - b*kvecs[j-1]
# Calculate diagonal term
a = np.dot(kvecs[j].conjugate(),t)
# Schmidt orthogonalize with respect to prev. vector
t[:] = t - a*kvecs[j]
# Calculate non-diagonal term
b = np.dot(t.conjugate(),t)**0.5
kmat[j,j] = a
if j < order-1:
kvecs[j+1] = (1./b)*t
kmat[j+1,j] = b
kmat[j,j+1] = b
return kmat,kvecs
[docs]
def ovrMatrix(vecs1,vecs2,nvecs,veclen,matrix):
"""
Compute the overlap matrix between two lists of vectors
Parameters
----------
vecs1 : array_like
must be conformable to shape=(nvecs,veclen)
vecs2 : array_like
must be conformable to shape=(nvecs,veclen)
nvecs : integer, number of vectors
veclen: integer, length of each vector
matrix: 2D array, output, overlap matrix
"""
saveshape = vecs1.shape
vecs1.shape = (nvecs,veclen)
vecs2.shape = (nvecs,veclen)
for i,v2 in enumerate(vecs2):
for j,v1 in enumerate(vecs1):
matrix[j,i] = np.vdot(v1,v2)
vecs1.shape = saveshape
vecs2.shape = saveshape
return matrix
[docs]
def ovrMatrix2(vecs1,vecs2,nvecs1,nvecs2,veclen,matrix):
"""
Compute the overlap matrix between two lists of vectors
Parameters
----------
vecs1 : array_like
must be conformable to shape=(nvecs1,veclen1)
vecs2 : array_like
must be conformable to shape=(nvecs2,veclen2)
nvecs1 : integer, number of vectors in vecs1
nvecs2 : integer, number of vectors in vecs2
veclen: integer, length of each vector
matrix: 2D array, output, overlap matrix
"""
saveshape1 = vecs1.shape
saveshape2 = vecs2.shape
vecs1.shape = (nvecs1,veclen)
vecs2.shape = (nvecs2,veclen)
for j,v1 in enumerate(vecs1):
for i,v2 in enumerate(vecs2):
matrix[j,i] = np.vdot(v1,v2)
vecs1.shape = saveshape1
vecs2.shape = saveshape2
return matrix
[docs]
def sortEigValsVecs(evals,evecs):
"""
Return the pair (eigvals,eigvecs) sorted from low to high values
"""
evecs[:] = evecs.take(evals.argsort(),axis=1) # axis=1: eigenvectors are expected in columns
evals[:] = evals.take(evals.argsort())
[docs]
def distributeMatrix(mat1,mat2,oper):
"""
Given two square matrices that operate on spaces S1 and S2,
return a single matrix that operates on the direct product
space S = S1 x S2
mat1 -- square matrix n*n
mat2 -- square matrix m*m
oper -- 's' for summ or 'p' for product
whether the resulting matrices are summed or
matrix multiplied
Example:
Given the matrices:
| 1 2 3 4| |a b c|
mat1 = | 5 6 7 8| mat2 = |d e f|
| 9 10 11 12| |g h i|
|13 14 15 16|
The matrix of dimension 12 that operates on the product space is
generated from matrices A and B where
| 1 2 3 4 | |a b c |
| 1 2 3 4 | |d e f |
| 1 2 3 4 | |g h i |
| 5 6 7 8 | | a b c |
| 5 6 7 8 | | d e f |
A = | 5 6 7 8 | B = | g h i |
| 9 10 11 12 | | a b c |
| 9 10 11 12 | | d e f |
| 9 10 11 12 | | g h i |
|13 14 15 16 | | a b c|
| 13 14 15 16 | | d e f|
| 13 14 15 16 | | g h i|
"""
n = len(mat1)
m = len(mat2)
nm = n*m
nmat1 = np.zeros((nm,nm),type(mat1[0,0]))
nmat2 = np.zeros((nm,nm),type(mat2[0,0]))
for k in range(nm):
for l in range(nm):
i = k%m
j = l%m
ii = k/m
jj = l/m
if i == j:
nmat1[k,l] = mat1[ii,jj]
if ii == jj:
nmat2[k,l] = mat2[i,j]
if oper == 's': return nmat1 + nmat2
if oper == 'p': return np.dot(nmat1,nmat2)
[docs]
def distributeMatrixList(ml,oper):
m2 = ml.pop()
while ml:
m1 = ml.pop()
m3 = distributeMatrix(m1,m2,oper)
m2 = m3
return m2
[docs]
def isPositiveDefinite(m):
"""Return 1 if m is positive definite, 0 if it's not"""
if type(m) != type(np.array([])):
return 0
elif np.rank(m) != 2:
return 0
try:
ch=np.linalg.cholesky(m)
except np.linalg.LinAlgError:
return 0
return 1
[docs]
def toPositiveDefinite(m,min=0.0,fac=1.0):
"""
Convert matrix m to positive definite
changes only negative or slightly positive eigenvalues.
1. Diagonalize M to get M=U.L.U+
2. Supress the negative eigenvalues of L, L->Lt
3. Get the new potivive definite matrix, Mt=U.Lt.U+
m -- matrix to treat
min -- at the end all eigenvalues will have a value at least of min
fac -- scale factor for the eigenvalues found negative
The default behaviour leaves any positive definite
matrix unaltered.
"""
fac = abs(fac)
eval,evecm = eig(m)
for i in range(len(eval)):
if eval[i] < 0.0:
eval[i] = -1.0*fac*eval[i]
if eval[i] < min:
eval[i] = min
evalm = np.identity(len(eval))*eval
mt = np.dot(np.dot(evecm,evalm),np.transpose(evecm))
return mt
[docs]
def invert_reg(mat,tol=1.0e-8):
"""
Invert matrix mat, regularizing it in case it has eigenvalues close to 0
"""
(evals,evecs) = np.linalg.eig(mat)
evalsreg = np.empty(evals.shape,evals.dtype)
evalsreg[:] = tol
evalsreg = np.fmax(evals,evalsreg)
ievalsreg = evalsreg**-1
ievalsmat = np.diag(ievalsreg)
inv = np.dot(evecs,np.dot(ievalsmat,evecs.transpose()))
return (inv,evals,evalsreg,evecs)
[docs]
def schroedinger_to_interaction(VS,H0,t):
"""
Return interaction picture matrix VI using diagonal matrix H0
VI = exp(i*H0*t) VS exp(-i*H0*t) Sakurai 5.5.7
"""
exp_p = np.diag( np.exp( 1j*np.diag(H0)*t ) )
exp_m = np.diag( np.exp(-1j*np.diag(H0)*t ) )
return np.dot( exp_p, np.dot(VS,exp_m) )
##############################################################################
#--- TENSOR ALGEBRA: functions related to operations with multiindex arrays --
##############################################################################
[docs]
def hprod_phi(hlist,nphi,phidim,sdim,phi,hole=None):
"""
Apply matrices in hlist to phi
|phi[i]> <= h1.h2.h2....hn |phi[i]>
hlist - list (or array) of 2D square matrices
nphi - number of tensors in phi
phidim- total length of each tensor phi
sdim - shape of each tensor phi[i]
phi - 1D array, a list of tensors phi[i] (input/output)
hole - integer, hlist[hole] is skipped
The operations are accumulated in phi, i.e. phi is modified in place.
- note: phi must be conformable to (nphi,phidim) or (nphi,prod(sdim))
"""
saveshape = phi.shape
phi.shape = (nphi,phidim)
for i in range(nphi):
for l,hmat in enumerate(hlist):
if l is hole: continue
matrix_tensor(l,sdim,hmat,phi[i])
phi.shape = saveshape
[docs]
def hprod_phi2(hlist,nphi,sdim1,sdim2,phi,hole=None):
"""
Apply matrices in hlist to phi
|traphi[i]> = h1.h2.h2....hn |phi[i]>
hlist - list (or array) of 2D matrices. NEED NOT BE SQUARE.
nphi - number of tensors in phi
phidim- total length of each tensor phi
sdim1 - shape of each tensor on output (first dimension of each matrix)
sdim2 - shape of each tensor on input (second dimension of each matrix)
phi - 1D array, a list of tensors phi[i] (input)
hole - integer, hlist[hole] is skipped
traphi - transformed phi (output)
"""
phic = phi.copy()
philen1 = 1
for i in sdim1: philen1 = philen1*i
philen2 = 1
for i in sdim2: philen2 = philen2*i
phic.shape = (nphi,philen2)
traphi = []
for i in range(nphi):
phiold = phic[i]
sdimold = list(sdim2)
for l,hmat in enumerate(hlist):
if l is hole: continue
phinew = matrix_tensor2(l,sdimold,hmat,phiold)
phiold = phinew
sdimold[l] = sdim1[l]
traphi.append(phiold)
traphi = np.array(traphi)
return traphi.flatten()
[docs]
def matrix_tensor(m,sdim,mat,phi):
"""
Perform a matrix tensor multiplication on index m of tensor phi
m - input, index summed over
sdim - input, dimensions (shape) of phi when in tesorized form
mat - input, 2D square array
phi - input/output, 1D array conformable to the sdim shape
phi(I,i,K) <= sum_j mat(i,j)*phi(I,j,K)
"""
ijk = ijk_shape(m,sdim)
saveshape = phi.shape
phi.shape = ijk
for i in range(ijk[0]):
for k in range(ijk[2]):
phi[i,:,k] = np.dot(mat,phi[i,:,k])
phi.shape = saveshape
[docs]
def matrix_tensor2(m,sdim,mat,phi):
"""
Perform a matrix tensor multiplication on index m of tensor phi
m - input, index summed over
sdim - input, dimensions (shape) of phi when in tesorized form
mat - input, 2D array, need not be square
phi - input/output, 1D array conformable to the sdim shape
The length of the multiplied index changes if the first dimension
of mat is not equal to the second dimension.
phiout(I,i,K) <= sum_j mat(i,j)*phi(I,j,K)
"""
ijk = ijk_shape(m,sdim)
ilk = list(ijk)
ilk[1] = mat.shape[0]
ilk = tuple(ilk)
phiout = np.empty(ilk,complex)
saveshape = phi.shape
phi.shape = ijk
for i in range(ijk[0]):
for k in range(ijk[2]):
phiout[i,:,k] = np.dot(mat,phi[i,:,k])
phi.shape = saveshape
return phiout.flatten()
[docs]
def vector_tensor(m,sdim,vec,phi):
"""
Perform a vector tensor multiplication on index m of tensor phi
m - input, index summed over
sdim - input, dimensions (shape) of phi when in tesorized form
mat - input, 1D
phi - input/output, 1D array conformable to the sdim shape
hphi(I,j,K) <= vec(j)*phi(I,j,K)
"""
ijk = ijk_shape(m,sdim)
saveshape = phi.shape
phi.shape = ijk
for i in range(ijk[0]):
for k in range(ijk[2]):
phi[i,:,k] = vec*phi[i,:,k]
phi.shape = saveshape
return phi
[docs]
def compute_density(dm,f,sdim,ten1,ten2=None):
"""
Calculate the density matrix for axis f of tensors ten1 and ten2
dm: density matrix, output
f: integer, axis to be used for the density
sdim: dimensions of ten1 and ten2
ten1: array conformable to ten1.shape=sdim
ten2: if not given ten2=ten1
"""
if ten2 is None: ten2=ten1
oldshape = ten1.shape
l = len(sdim)
ax = range(f)+range(f+1,l)
ten1.shape = sdim
ten2.shape = sdim
dm[:] = np.tensordot(ten1.conjugate(),ten2,axes=(ax,ax))
ten1.shape = oldshape
ten2.shape = oldshape
[docs]
def overlap(avecs,bvecs):
"""
Compute the overlap matrix from two two-dimensional np.array
avecs and bvecs need only to conform in the second index (index 1)
"""
return np.tensordot(avecs.conjugate(),bvecs,axes=(1,1))
[docs]
def ijk_shape(f,oshape):
"""
change tuple oshape representing a shape to a 3 index tuple collecting
indices before and after index f
"""
nmodes = len(oshape)
a = np.array(oshape)
nshape=[1,1,1]
nshape[1] = oshape[f]
if f > 0:
nshape[0] = a[:f].prod()
if f < nmodes - 1:
nshape[2] = a[f+1:].prod()
return tuple(nshape)
##############################################################################
# --- VECTORS in 3D SPACE: vectors in cartesian space ---
##############################################################################
[docs]
def eulerMatrix(a,b,c,inv=False):
"""Return a 3x3 rotation matrix
xyz is the body fixed frame
XYZ is the space fixed frame
The normal convention (default) is the rotation XYZ -> xyz
First rotate around Z by an ammount of xhi,
then around Y by an ammount of theta
and finally around Z by an ammount of phi.
The theta and phi angles can also be seen as the two spherical angles
of a vector with respect to the Z vector.
phi azimutal spherical angle (0,2pi)
theta polar spherical angle (0,pi)
phi -- rotation around Z axis
theta -- rotation around Y axis
xhi -- rotation around Z axis
The inverse convention corresponds to xyz -> XYZ
The vector is rotated by the specified Euler angles that align it
with the space fixed frame.
Normal and Inverse conventions are inverse operations with respect to
each other
"""
A=np.zeros((3,3),float)
B=np.zeros((3,3),float)
C=np.zeros((3,3),float)
A[0,0] = cos(a)
A[1,1] = cos(a)
A[2,2] = 1.0
B[0,0] = cos(b)
B[1,1] = 1.0
B[2,2] = cos(b)
C[0,0] = cos(c)
C[1,1] = cos(c)
C[2,2] = 1.0
if not inv:
A[0,1] = -sin(a)
A[1,0] = sin(a)
B[0,2] = sin(b)
B[2,0] = -sin(b)
C[0,1] = -sin(c)
C[1,0] = sin(c)
r=np.dot(np.dot(A,B),C)
else:
A[0,1] = sin(a)
A[1,0] = -sin(a)
B[0,2] = -sin(b)
B[2,0] = sin(b)
C[0,1] = sin(c)
C[1,0] = -sin(c)
r=np.dot(np.dot(C,B),A)
return r
[docs]
def alignMatrix(v,axis='z'):
"""Return a 3x3 rotation matrix
Given a 3D vector, the routine returns a matrix that,
when multiplied by v, rotates it to the specified axis.
The routine computes the rotation axis around which the
rotation has to be performed, the angle that has to be
rotated, and then use the Rodrigues formula obtain the
desired rotation matrix.
"""
if axis == 'x':
v1=np.array([1.,0.,0.])
elif axis == 'y':
v1=np.array([0.,1.,0.])
elif axis == 'z':
v1=np.array([0.,0.,1.])
else:
return np.zeros((3,3),float)
v2 = np.cross(v,v1)
if sqrt(np.dot(v2,v2)) < 1.e-8:
return np.identity(3,float)
cosphi = np.dot(v,v1)/sqrt(np.dot(v,v))
phi = acos(cosphi)
return rodriguesMatrix(v2,phi)
[docs]
def rodriguesMatrix(v,phi):
"""Return a 3x3 rotation matrix
v -- 3D vector, the rotation axis (a copy will be
normalized in the routine)
phi -- ammount of rotation
Given the vector v=[v0,v1,v2] and an angle phi, it returns the
rotation matrix around that axis by the ammount phi.
The routine uses the Rodrigues Formula:
M = I + sin(phi)*A + (1-cos(phi))*A^2
where I is the unit matrix and A is given by:
| 0 -v2 v1|
A = | v2 0 -v0|
|-v1 v0 0|
"""
vn = v/sqrt(np.dot(v,v))
id=np.identity(3,float)
a=np.zeros((3,3),float)
a[0,1]=-vn[2]
a[0,2]= vn[1]
a[1,0]= vn[2]
a[1,2]=-vn[0]
a[2,0]=-vn[1]
a[2,1]= vn[0]
aa=np.dot(a,a)
m = id + sin(phi)*a + (1-cos(phi))*aa
return m
[docs]
def pointsDistance(p,q):
"""
Return the distance (Euclidean norm of connecting vector) between two points
p - vector or list
q - vector or list
"""
d=0.0
for i in range(len(p)):
d = d + (p[i]-q[i])**2
return d**0.5
[docs]
def pointsVector(p,q):
"""
Return the vector q ---> p (p-q) between points p and q
p - vector or list
q - vector or list
on return v is a list
"""
v = []
for i in range(len(p)):
v.append(p[i]-q[i])
return v
[docs]
def sphericalAngles(v):
"""
return the spherical angles theta(polar) and phi(azimutal) of v
v - 3D vector (x,y,z)
"""
r=pointsDistance(v,[0.0,0.0,0.0])
theta = acos(v[2]/r)
rxy = r*sin(theta)
phi = acos(v[0]/rxy)
if v[1] < 0.0:
phi=2*pi-phi
return theta,phi
[docs]
def leviCivita(i,j,k):
"""
Return the value of the antisymmetric permutation tensor for indices
i,j,k
i,j,k - ints
"""
idty = np.identity(3)
return np.dot(idty[i], np.cross(idty[j], idty[k]))
##############################################################################
# --- INTERPOLATION on a regular grid ---
##############################################################################
[docs]
def ndimIntGrid(xgrid1,ygrid1,xgrid2,intfunctions):
"""Interpolate xgrid2 into the grid defined by xgrid1 and ygrid1
xgrid1 -- 2D list containing arrays or lists of the original coordinates
ygrid1 -- ND list or array containing the values of the function at
each point.
ygrid1[i1,i2,i3,...,in] = y( xgrid1[1,i1],xgrid1[2,i2],...,xgrid1[n,in] )
xgrid2 -- 2D list or array containing the arrays of the coordinates
to be interpolated.
intfunctions -- list containing the pointers to the functions that will handle
the interpolation of each degree of freedom
Each element of the list is also a list, the first
element being the pointer and the rest the arguments
expected by the interpolation routine, i.e.
((f1,op1a,op1b),(f2,op2a,op2b,op2c),...)
"""
# declare the ygrid2 output array
# l contains the number of points in each dimension
l = []
for xvec in xgrid2:
l.append(len(xvec))
ygrid2 = np.zeros(l,float)
# loop over all the points in xgrid2 and interpolate
ygrid1a = np.array(ygrid1)
p = MultiCounter()
p.setPointersN(0,len(l))
p.setMaxIter(l)
for ic in p:
i=0
xpoint = []
for xvec in xgrid2:
xpoint.append(xvec[ic[i]])
i = i+1
ygrid2[ic] = ndimIntPoint(xgrid1,ygrid1a,xpoint,intfunctions)
return ygrid2
[docs]
def ndimIntPoint(xgrid1,ygrid1,xpoint,intfunctions):
"""Iterpolate xpoint into the grid defined by xgrid1 and ygrid1
xgrid1 -- 2D list containing arrays or lists of the original coordinates
ygrid1 -- ND np.array with the function values
xpoint -- list or array of the coordinates of the point to be
interpolated
intfunctions -- list with pointers to the functions that will handle
the interpolation of each degree of freedom
Each element of the list is also a list, the first
element being the pointer and the rest the arguments
expected by the interpolation routine
Example:
We have a 3D grid, each coordinate (q1,q2,q3) is defined between
0.0 and 5.0. Each coordinate is sampled at 11 equispaced points.
We have a mapping f(q1,q2,q3) from the grid points to the real numbers.
The specific form of the mapping is of course irrelevant here, we only
are interested on the value corresponding to each grid point.
Thus xgrid1 contains:
xgrid1 = [ [0.0, 0.5, 1.0, 1.5, ..., 5.0],
[0.0, 0.5, 1.0, 1.5, ..., 5.0],
[0.0, 0.5, 1.0, 1.5, ..., 5.0]
]
And ygrid1 contains:
ygrid1[i][j][k] = f(xgrid1[0][i],xgrid1[1][j],xgrid1[2][k])
And we provide also an interpolation point from which we want to
know the interpolation f_intp(xpoint):
xpoint = [1.3,2.2,3.7]
"""
ygridtmpA = np.array(ygrid1)
ndim = len(xpoint)
l = ygrid1.shape
# at each step we get a new grid of 1 dimension less until we reach
# the final interpolated number
for d in range(ndim):
xg = np.array(xgrid1[d])
xp = xpoint[d]
if l[d+1:]:
ygridtmpB = zeros(l[d+1:],float)
p = MultiCounter()
p.setPointersN(0,len(l[d+1:]))
p.setMaxIter(l[d+1:])
# loop over the (ndim-d+1) dimensional grid
for ic in p:
yg = zeros(l[d],float)
# assign y values for the 1D interpolation
for j in range(l[d]):
k = (j,)
yg[j] = ygridtmpA[k+ic]
ygridtmpB[ic] = onedimInt(xg,yg,xp,intfunctions[d])
ygridtmpA = ygridtmpB
else: # last iteration
ygridtmpB = 0.0
yg = ygridtmpA
ygridtmpB = onedimInt(xg,yg,xp,intfunctions[d])
return ygridtmpB
return None
[docs]
def scattered2grid2D(points,xgrid,ygrid,kernel,**args):
"""
Interpolate scattered 2D data into a regular grid
Arguments:
kernel -- w(d)
z = Sum_i w(d_i)
"""
griddata = np.zeros((len(xgrid),len(ygrid)),float)
for i,x in enumerate(xgrid):
for j,y in enumerate(ygrid):
z = 0.0
for p in points:
d = (p[0]-x)**2 + (p[1]-y)**2
d = d**0.5
z = z + kernel(d)
griddata[i,j] = z
return griddata
##############################################################################
# 1D interpolators
##############################################################################
[docs]
def onedimInt(xg,yg,xp,intfunc):
"""Interpolate point xp using the yg[xg] function
this routine is a thin wrapper over 1D interpolation routines
xg -- array with tabulated x points
yg -- array with tabulated y points
xp -- point where the interpolation is desired
func -- list of one or more elements, the first is the pointer
to the interpolation function, the rest are possible
arguments expected by such function
If the last argument is "periodic", a periodic interpolation
will be atempted enlarging the xg segment, usually up to 4
more points per side to be safe with respect to the cubic
spline interpolation
"""
f = intfunc[0]
ops = intfunc[1:]
if ops:
if ops[len(ops)-1] == 'periodic':
xg,yg = _toPeriodic(xg,yg)
return f(xg,yg,xp,ops[:-1])
return f(xg,yg,xp,ops)
else:
return f(xg,yg,xp)
[docs]
def linealInterpolation(xg,yg,xp,args=None):
"""1D lineal interpolation of x
xg -- array with tabulated x points
yg -- array with tabulated y points
xp -- point where the interpolation is desired
args -- arguments, it is ignored and does not need to be set,
only keeps the interface consistency of 1D interpolators
"""
imin = _segment(xg,xp)
y = yg[imin] + ((yg[imin+1]-yg[imin])/(xg[imin+1]-xg[imin]))*(xp-xg[imin])
return y
[docs]
def cubicSplineInterpolation(xg,yg,xp,args=None):
"""Cubic Spline interpolation
xg -- array with tabulated x points
yg -- array with tabulated y points
xp -- point where the interpolation is desired
args -- ignored by now
"""
d11 = 1e+40
d1n = 1e+40
y2 = np.zeros(len(xg),float)
u = np.zeros(len(xg),float)
splinenm(xg,yg,d11,d1n,y2,u)
return spline1d(xg,yg,y2,xp)
[docs]
def ENOInterpolation(xg,yg,xp,args=None):
"""Essentially Non-oscilatory Interpolation
xg -- array with tabulated x points
yg -- array with tabulated y points
xp -- point where the interpolation is desired
args -- list of arguments, can be a list of 0,1,2 or three velues
args = [ [[[min. order], max. order], factor ]
min. order -- minimum order of interpolation, def=1
max. order -- maximum order of interpolation, def=3
factor -- balance factor, the larger it is, the more likely
a symmetric interpolation is made, loosing the gains
of ENO. If too small, small changes in slope are always
avoided, thus giving interpolants with discontinous
derivatives.
def: 50.0
"""
if args:
try: min = args[0]
except IndexError: min = 1
try: max = args[1]
except IndexError: max = 3
try: fact = float(args[2])
except IndexError: fact = 50.0
else:
min = 1
max = 3
fact = 50.0
diff=np.zeros((len(xg),len(xg)),float)
return enoi(xg,yg,xp,diff,max,min,fact)
# A dictionary is defined to translate from the input file
# code of an interpolation method (usually 3 characters)
# and the pointer to the real function, add entries always
# when a new 1D interpolator is implemented
keyToFunction = {}
keyToFunction["default"] = linealInterpolation
keyToFunction["lin"] = linealInterpolation
keyToFunction["eno"] = ENOInterpolation
keyToFunction["spl"] = cubicSplineInterpolation
##############################################################################
# --- NUMERICAL functions ---
##############################################################################
[docs]
def dkron(i,j):
"""Return the delta of Kroneker operation of integers i and j
i,j - integers, in case they are not integers the function
silently returns 0
"""
if isinstance(i,int) and isinstance(j,int):
if i==j:
return 1
else:
return 0
else:
return 0
[docs]
def richardson(f,hi,con=1.4,niter=10,safe=2.0,fast=False):
"""
Implements Richardson extrapolation to calculate numerical estimates
of values obtained by finite difference formulas.
f -- function depending on 1 variable, the finite difference, what f is,
a derivative, a numerical integral... does not play any role here
hi -- initial estimate for the finite difference
con -- constant factor by which hi is divided at each step
niter -- maximum number of iterations
safe -- when error grows by a factor more than safe, terminate and return
"""
if fast:
if hi > 0.05:
return f(0.05),0.0,hi
else:
return f(hi),0.0,hi
con2 = con*con
big = 1.e30
a=np.zeros((niter,niter),float)
hh=hi
a[0,0] = f(hh)
err = big
for i in range(1,niter):
hh = hh/con
a[0,i] = f(hh)
fac = con2
for j in range(1,i+1):
a[j,i] = (a[j-1,i]*fac - a[j-1,i-1])/(fac-1.0)
fac = con2*fac
errt = max(abs(a[j,i]-a[j-1,i]),abs(a[j,i]-a[j-1,i-1]))
if errt <= err:
err = errt
res = a[j,i]
hf = hh
if abs(a[i,i]-a[i-1,i-1]) >= safe*err:
return res,err,hf
return res,err,hf
##############################################################################
# --- Sampling & Monte-Carlo functions ---
##############################################################################
[docs]
def WignerStep(M_nmode,w_nmode,qnum_nmode,xref,normalModes,min=None,max=None):
"""
Make a random movement for generating Wigner distribution for sampling normal modes.
M_nmode - N masses of normal modes
w_nmode - N energies of normal modes
qnum_nmode - N quantum numbers of normal modes
return
xnext - coordinates and velocities
"""
ncoords = len(xref)
nmodes = len(w_nmode)
xnext = np.zeros(2*ncoords)
qq = np.zeros(nmodes)
vv = np.zeros(nmodes)
for i in range(nmodes):
q,v = WignerStep1D(M_nmode[i], w_nmode[i], qnum_nmode[i])
qq[i] = q
vv[i] = v
xnext[0:ncoords] = np.asarray(xref) + np.dot(normalModes.transpose(), qq)
xnext[ncoords:] = np.dot(normalModes.transpose(), vv)
return xnext
[docs]
def WignerStep1D(m, w, qnum):
"""
Make a random movement to generate a 1D Wigner distribution of a harmonic oscillator.
Input:
m - mass of the normal mode
w - frequency of the normal mode
qnum - quantum number (not used!)
Output:
QQ - new position (still in normal mode representation)
VV - new velocity (still in normal mode representation)
"""
while True:
P = np.random.uniform(-3.0,3.0,1)
Q = np.random.uniform(-3.0,3.0,1)
W = np.exp(-(P**2)) * np.exp(-(Q**2))
r = np.random.uniform(0.0,1.0,1)
if W > r:
QQ = Q / np.sqrt( w * m )
VV = P * np.sqrt( w / m)
return QQ, VV
[docs]
def randomWalkStep(xvec,min=None,max=None,stepnorm=1.0):
"""
Make a random movement on a vector of real (float) coordinates
xvec - current state of the system in coordinate space
min - minimum value for each coordinate
max - maximum value for each coordinate
stepnorm - length of random steps to be performed
return
- next random position
"""
lx = len(xvec)
dual_step = False
try:
if len(stepnorm) == 2:
dual_step = True
except:
dual_step = False
rstep = np.random.random(lx) - 0.5
norm = np.linalg.norm(np.split(rstep,2), axis=1)
if dual_step == False:
if lx < 5:
rstep[0:lx//2] = (stepnorm)*rstep[0:lx//2]
rstep[lx//2:] = (stepnorm)*rstep[lx//2:]
else:
rstep[0:lx//2] = (stepnorm/norm[0])*rstep[0:lx//2]
rstep[lx//2:] = (stepnorm/norm[1])*rstep[lx//2:]
else:
if lx < 5:
rstep[0:lx//2] = (stepnorm[0])*rstep[0:lx//2]
rstep[lx//2:] = (stepnorm[1])*rstep[lx//2:]
else:
rstep[0:lx/2] = (stepnorm[0]/norm[0])*rstep[0:lx//2]
rstep[lx//2:] = (stepnorm[1]/norm[1])*rstep[lx//2:]
xvecnew = xvec + rstep
if not max is None:
delta = max - xvecnew
hit = (delta < 0)
if hit.any():
for i in xrange(lx):
if hit[i]: xvecnew[i] = max[i] + delta[i]
if not min is None:
delta = min - xvecnew
hit = (delta > 0)
if hit.any():
for i in xrange(lx):
if hit[i]: xvecnew[i] = min[i] + delta[i]
return xvecnew
[docs]
def metropolisStep(xvec,pval,pfunc,pfuncargs=None,
min=None,max=None,stepnorm=1.0):
"""
Make a random movement, accept or reject accordint to metropolis alg.
xvec - xvec state of the vector of coordinates (list, tuple or array)
pval - probability value associated to the xvec position
pfunc - maps the xvec position in N-Dim space to a probability
pfuncargs - passed to pfunc, put pfunc arguments there
return
1. True if accepted, False if rejected
2. accepted or old xvec position
3. probability accepted point, pval retuned if not acc.
"""
xnext = randomWalkStep(xvec,min=min,max=max,stepnorm=stepnorm)
if pfuncargs is None:
pnext = pfunc(xnext)
else:
pnext = pfunc(xnext,pfuncargs)
if pnext > pval:
return True,xnext,pnext
else:
p = pnext/pval
r = random.random()
if r < p:
return True,xnext,pnext
else:
return False,xvec,pval
[docs]
def randomWalkIntStep(xvec,xvecsize):
"""
Make a random movement on a vector of integer coordinates
xvec - current state of the vector of integers (list, tuple or array)
xvecsize - vector with the maximum allowed value per each coordinate
return
1. list with the next position
2. index of the changed coordinate
3. +1 if changed in positive direction, -1 if negative
"""
next = list(xvec)
n = len(xvec) - 1
i = random.randint(0,n)
s = 2*random.randint(1,2) - 3
next[i] = next[i] + s
if next[i] > xvecsize[i] or next[i] < 0:
next,i,s = randomWalkStep(xvec,xvecsize)
return next,i,s
[docs]
def metropolisIntStep(xvec,xvecsize,xvecVal,pfunc,pfuncargs=None):
"""
Make a random movement, accept or reject accordint to metropolis alg.
xvec - xvec state of the vector of integers (list, tuple or array)
xvecsize - vector with the maximum allowed value per each coordinate
xvecVal - probability value associated to the xvec position
pfunc - maps the xvec position in N-Dim space to a probability
pfunc must accept 2 arguments, xvec and pfuncargs
pfuncargs - passed to pfunc, put pfunc arguments there
return
1. True if accepted, False if rejected
2. list with the accepted or old position
3. index of the changed coordinate, (-1 if no change)
4. +1 if changed in positive direction, -1 if negative, 0 if no change
5. function value of the accepted point, xvecVal retuned if not acc.
"""
xnext,i,s = randomWalkStep(xvec,xvecsize)
xnextVal = pfunc(xnext,pfuncargs)
if xnextVal > xvecVal:
return True,xnext,i,s,xnextVal
else:
p = xnextVal/xvecVal
r = random.random()
if r < p:
return True,xnext,i,s,xnextVal
else:
return False,xvec,-1,0,xvecVal
[docs]
def WignerWalkQuantumStep(qnum,M_amu,w,radius2,min=None,max=None):
"""
Make a random movement on a vector of real (float) coordinates
bound on a 1D sphere of radius [radius]
radius2 - radius squared of the sphere
min - minimum value for each coordinate
max - maximum value for each coordinate
return
- next random position
"""
M = M_amu * au2am
radius = sqrt(2.0 *radius2)
accepted = False
if qnum == 0:
qref = 0
else:
qref = radius # [T]=[V] for HO
pref = qref
while not accepted:
pos = random.random() - 0.5
q1D = radius * sin(pi*pos)
p1D = 1.0/radius * cos(pi*pos) * np.sign(pos)
# von Neumann scheme
#print (q1D**2 + p1D**2)
if qnum == 0:
ratio = (1/np.pi) * exp( - q1D**2 - p1D**2 )
else:
ratio = (1/np.pi) * exp( - ( q1D - qref )**2 - (p1D - pref)**2 )
if ratio > random.random():
accepted = True
xvecnext1D = np.array([q1D, p1D], dtype=float)
return xvecnext1D
[docs]
def WignerQuantumStep(M_nmode,w_nmode,qnum_nmode,min=None,max=None):
"""
Make a random movement for generating Wigner distribution
M_nmode - N masses of normal modes
w_nmode - N energies of normal modes
qnum_nmode - N quantum numbers of normal modes
return
xnext - coordinates and momenta
"""
hbar = 1.0
ncoords = len(w_nmode)
xnext = np.zeros(2*ncoords, dtype=float) # [q,p] 2N vector
for i in range(ncoords):
# energy of the normal modes
ET_nmode = hbar * (float(qnum_nmode[i]) + 0.5) * float(w_nmode[i])
#print ET_nmode
# generate a 2-unit vector on sphere of radius squared [ET_nmode]
# [q,p]
xnext1D = WignerWalkQuantumStep(qnum_nmode[i],M_nmode[i],w_nmode[i],ET_nmode,min=min,max=max)
# scale unit of q and p backwards
xnext[i] = xnext1D[0] #* sqrt(2.0) # q
xnext[i+ncoords] = xnext1D[1] #* sqrt(2.0) # p
#
return xnext
##############################################################################
# --- Statistical Functions ---
##############################################################################
[docs]
def rmsError(list1,list2,energy=30.0):
global table
pardict = parametersS0.paramsFromList(parlist)
err = 0.0
n = 0
for q in table:
if(q[7]<energy):
err = err + (potevb(q[0:7],pardict)-q[7])**2
n = n+1
err = err/float(n)
err = err**0.5
return err,n
##############################################################################
# --- Orthonormal functions ---
##############################################################################
[docs]
def generateOrthonormalSpace(vs,n):
"""
Given an initial vector generate n-1 additional orthonormal vectors
Discrete orthonormality is preserved:
dot(v_n,v_m) = delta(n,m)
v0 - initial vector
n - size of orthonormal space generated
"""
m = len(vs)
if n > m:
raise ValueError("n > m, maximum m orthonormal vectors can be generated")
v0 = np.asarray(vs)
vspace = np.zeros((n,m),type(v0[0]))
vspace[0] = v0
x = np.linspace(0.0,1.0,m+2)[1:-1]
for i in range(1,n):
vspace[i] = np.sin((i+1)*np.pi*x) # particle in a box excited states
gramSchmidt(vspace)
return vspace
[docs]
def gramSchmidt(vecs):
"""
Perform a Gram-Schmidt orthonormalization on a given list of arrays
vecs -- list of arrays to orthonormalize
"""
vecs[0] = vecs[0] / sqrt(np.vdot(vecs[0],vecs[0]))
if len(vecs) > 1:
for i in range(1,len(vecs)):
for j in range(0,i):
vecs[i] = vecs[i] - np.vdot(vecs[j],vecs[i])*vecs[j]
vecs[i] = vecs[i] / np.sqrt(np.vdot(vecs[i],vecs[i]))
##############################################################################
# --- Function analysis ---
##############################################################################
[docs]
def convolutionFilter(input,fwhm,step=None,lbound=None,hbound=None,
filter='gaussian',externalFilter=None):
"""
Filter the x,y input array with a gaussian filter
input -- x,y input array
fwhm -- full width at half maximum of the filter
step -- discretization interval in the abcisas axis
lbound -- lower bound of for the output data
hbound -- higher bound for the output data
filter -- string: type of function used for filtering, gaussian
and lorentzian shapes are implemented
externalFilter -- externally provided function for filtering,
must take two parameters:
x - value of the abcisa
b - center of the filter
return: an x,y array with the convolved function
"""
if externalFilter:
f=externalFilter
else:
if filter == 'gaussian':
def f(x,b):
k = 4.0*math.log(2.0)
N = 2.0*math.log(2.0)/np.sqrt(math.pi)
N = N/fwhm
return N*np.exp(-k*((x-b)/fwhm)**2)
elif filter == 'ngaussian':
def f(x,b):
sigma = fwhm / (2.*np.sqrt(2.*math.log(2.)))
N = 1./sigma/np.sqrt(2*np.pi)
return N*np.exp(- (x-b)**2 / 2./ sigma**2)
elif filter == 'lorentzian':
def f(x,b):
k=0.5*fwhm
k=k*k
return k/((x-b)**2 + k)
elif filter == 'nlorentzian':
def f(x,b):
k=0.5*fwhm
return k / ((x-b)**2 + k**2) / np.pi
if lbound is None: lbound=input[0][0]
if hbound is None: hbound=input[0][-1]
if step is None: step=fwhm*0.05
xvec = np.arange(lbound,hbound+step,step)
nstep = len(xvec)
xrange = xvec[-1] - xvec[0]
yinp = np.zeros(nstep,float)
# Map input data to the inp array
for i,x in enumerate(input[0]):
xi = (x-lbound)/step
ii = int(np.round(xi))
yinp[ii] = yinp[ii] + input[1][i]
# I(w) = Sum_t Func(t)*Filt(w-t)
yout = np.array([np.dot(f(t,xvec),yinp) for t in xvec])
return np.array([xvec,yout])
[docs]
def voigtConvolution(input, fwhm_gauss, fwhm_lor, step=None, lbound=None, hbound=None):
"""
Compute the convolution with the Voigt profile
input -- x,y input array
fwhm_gauss -- full width at half maximum of the Gaussian function
fwhm_lor -- full width at half maximum of the Lorentz function
step -- discretization interval in the abcisas axis
lbound -- lower bound of for the output data
hbound -- higher bound for the output data
return: an x,y array with the convolved function
"""
def f(x,b):
# Voigt profile
sigma = fwhm_gauss / (2.*np.sqrt(2.*math.log(2.)))
z = ((x-b) + 1j*fwhm_lor) / sigma / np.sqrt(2.)
# complex error function from scipy.special.wofz
return np.real(wofz(z)) / sigma / np.sqrt(2*np.pi)
if lbound is None: lbound=input[0][0]
if hbound is None: hbound=input[0][-1]
if step is None: step=min(fwhm_gauss, fwhm_lor)*0.05
xvec = np.arange(lbound,hbound+step,step)
nstep = len(xvec)
xrange = xvec[-1] - xvec[0]
yinp = np.zeros(nstep,float)
# Map input data to the inp array
for i,x in enumerate(input[0]):
xi = (x-lbound)/step
ii = int(np.round(xi))
yinp[ii] = yinp[ii] + input[1][i]
# I(w) = Sum_t Func(t)*Filt(w-t)
yout = np.array([np.dot(f(t,xvec),yinp) for t in xvec])
return np.array([xvec,yout])
[docs]
def Convolution(f, g, dt, **option):
"""
Compute convolution integral. It can be used to convolve two functions of
different length (number of points), but the initial and final domain points
(x or t coordinate) must coincide. The resulting function has the same
initial and final domain points (x or t coordinate) as the input functions.
If the input functions (f and g) have different spacing than output function
(specified by dt), spline interpolation of order kspline will be used to
approximate the necessary points.
Parameters:
-----------
f, g : 1D array_like
function to be convolved.
dt : float
interval of the co-domain (convolved) function.
method : string, default: 'simp', optional
integration method to be used.
'trapz' for trapezoidal rule,
'simp' for Simpson rule
npoints : int, default: max(len(f), len(g))
length of output function.
kspline : int {1-5}, default: 3, optional
spline order of interpolation of f and/or g (if the interpolation is to
be done). Interpolation is to be done on f/g if npoints is different
compare to len(f)/len(g).
Returns
-------
fog : 1D numpy.ndarray
convolution integral of f and g function.
"""
# evaluating options
method = option.get('method', 'simps')
npoints = option.get('npoints',max(len(f), len(g)))
kspline = option.get('kspline', 3)
# making domain of f and g which has the dt spacing
tspan = dt * (npoints - 1)
tline = np.linspace(0.0, tspan, npoints)
# ensure that the functions are in numpy array
func_f = np.asarray(f)
func_g = np.asarray(g)
# do interpolation if necessary
if len(func_g) != npoints:
gline = np.linspace(0.0, tspan, len(func_g))
fspline = sip.InterpolatedUnivariateSpline(gline, func_g, k=kspline)
func_g = fspline(tline)
if len(func_f) != npoints:
fline = np.linspace(0.0, tspan, len(func_f))
fspline = sip.InterpolatedUnivariateSpline(fline, func_f, k=kspline)
func_f = fspline(tline)
rfunc_g = ReverseArray(func_g)
# integration method
int_method = {'simp':si.simps, 'trapz':np.trapz}
if method != 'simp' and method != 'trapz':
method = 'simp'
convolution = np.zeros(npoints)
# if number of points is odd
if npoints % 2:
shifting_const = (npoints - 1) / 2
else:
# if number of points is even
shifting_const = npoints / 2
# do integration for each shifting
for i in xrange(npoints):
shifting = i - shifting_const
integrand = func_f * ShiftArray(rfunc_g, shifting, periodic=False)
convolution[i] = int_method[method](integrand, dx = dt)
return convolution
##############################################################################
# --- OTHER functions ---
##############################################################################
[docs]
def wignerHarmonic1D(q,p,m,w,hbar=1.0,q0=0.0,p0=0.0):
"""
Compute a probability density in phase space for the harmonic case.
Given a point in phase space compute the probability density using
the wigner transform of the vibrational ground state in harmonic
approximation. The wigner transform of a gaussian is another gaussian.
If hbar is not specified atomic units are assumed.
q -- position in phase space
p -- momentum in phase space
m -- mass
w -- frequency of the harmonic oscillator
"""
qexp = math.exp(-(m*w/hbar)*(q-q0)**2)
pexp = math.exp(-(1.0/(m*w*hbar))*(p-p0)**2)
return (1.0/(hbar*math.pi))*qexp*pexp
[docs]
def wignerE(q,p,m,w,hbar=1.0):
e = 0.0
for i in range(len(q)):
e += (m*w/hbar)*(q[i])**2
e += (1.0/(m*w*hbar))*(p[i])**2
return e
[docs]
def wignerEFT(q,p,m,w,temp,hbar=1.0):
e = 0.0
for i in range(len(q)):
f = math.tanh((hbar*w)/(2.0*KBAU*TEMP))
e += (f*m*w/hbar)*(q[i])**2
e += (f/(m*w*hbar))*(p[i])**2
return e
[docs]
def wignerHarmonicND(q,p,m,w,hbar=1.0):
"""
Compute a probability density in phase space for the N-Dim harmonic case.
Given a point in phase space compute the probability density using
the wigner transform of the vibrational ground state in harmonic
approximation. The wigner transform of a gaussian is another gaussian.
If hbar is not specified atomic units are assumed.
q -- list of positions in phase space
p -- list of momenta in phase space
m -- list of masses
w -- list of harmonic frequencies
"""
wprob = 1.0
for i in range(len(q)):
wprob = wprob*wignerHarmonic1D(q[i],p[i],m[i],w[i],hbar)
return wprob
[docs]
def wignerHarmonicFT1D(q,p,m,w,beta,hbar=1.0,q0=0.0,p0=0.0):
"""
Compute a probability density in phase space at finite temperature
If hbar is not specified atomic units are assumed.
q -- position in phase space
p -- momentum in phase space
m -- mass
w -- frequency of the harmonic oscillator
beta -- 1/kbT
"""
f = math.tanh((beta*hbar*w)/(2.0))
qexp = math.exp(-(f*m*w/hbar)*(q-q0)**2)
pexp = math.exp(-(1.0*f/(m*w*hbar))*(p-p0)**2)
return (1.0*f/(hbar*math.pi))*qexp*pexp
[docs]
def wignerHarmonicFTND(q,p,m,w,beta,hbar=1.0):
"""
Compute a probability density in phase space for the N-Dim harmonic case
at finite temperature
If hbar is not specified atomic units are assumed.
q -- list of positions in phase space
p -- list of momenta in phase space
m -- list of masses
w -- list of harmonic frequencies
beta -- 1/kbT
"""
wprob = 1.0
for i in range(len(q)):
wprob = wprob*wignerHarmonicFT1D(q[i],p[i],m[i],w[i],beta,hbar)
return wprob
[docs]
def corrspec(corr,dt,**opts):
"""
Compute spectrum of a correlation function
Parameters
----------
corr : array_like
1D complex array with correlation function.
dt : float
time interval between correlation values; [au]
Options
-------
eref : float; 0.0
energy shift, e.g. zero point energy; [au]
emin : float; 0.0
start of spectrum counting from eref; [au]
emax : float; 0.0 (calculated from dt if not given)
end of spectrum counting from eref; [au]
npoints : int; 500
number of discrete points where spectrum is evaluated
between emin and emax
EP : bool; False
if True, multiplication of spectrum with energy factor
ncos : int; 1
exponent of the cosine damping function
"""
eref = opts.get('eref',0.0)
emin = opts.get('emin',0.0)
emax = opts.get('emax',0.0)
npoints = opts.get('npoints',500)
EP = opts.get('EP',False)
ncos = opts.get('ncos',1)
a = np.asarray(corr)
maxwidth = 2.0*math.pi/dt
eminspec = eref + emin
emaxspec = eref + emax
if emax <= emin:
emaxspec = eminspec + maxwidth
T = (len(a) - 1)*dt
tlist = np.linspace(0.0,T,len(a))
elist = np.linspace(eminspec,emaxspec,npoints)
spect = np.zeros(len(elist),float)
cos_damp = np.cos(0.5*math.pi*tlist/T)**ncos
adamp = a*cos_damp
for i,energy in enumerate(elist):
tmp = energy*tlist*1.0j
tmp = np.exp(tmp)
tmp = tmp*adamp
spect[i] = tmp.sum()*dt
spect[i] = spect[i].real
if EP:
spect = spect*(elist-eref)
return elist-eref,spect
[docs]
def makeMesh(f,grid,ctype=float):
"""
make a mesh of f at the positions of a product grid
f -- function depending on N variables
grid -- tuple or list of N 1D arrays which define the primitive grid
on return: a N-D array with the values of f at the grid points.
"""
from .Iterators import MultiCounter
dim = len(grid)
indexes = []
for i in range(dim):
indexes.append(len(grid[i]))
z=np.zeros(indexes,ctype)
x=np.zeros(dim,float)
mc = MultiCounter()
mc.setPointersN(0,dim)
mc.setMaxIter(indexes)
for p in mc:
for i in range(dim):
x[i]=grid[i][p[i]]
z[p]=f(x)
return z
[docs]
def bssrSumOfStates(freqs,emax,m,E0=0.0):
"""
Beyer-Swinehart-Stein-Rabinovitch (BSSR) sum of states for harmonic systems
freqs: list of vibrational frequencies (hbar * omega)
emax: maximum energy related to the ZPE
m: discretization for the sum
E0: energy offset (can be a reaction barrier if dealing with a TS)
On return four lists containing:
E : binned energy units
N(E) : number of states at energy E (+- dE)
rho(E) : density of states at energy E
S(E) : sum of states up to energy E
"""
t = np.zeros(m+1,float)
ta = np.zeros(m+1,float)
e = np.zeros(m+1,float)
f = np.array(freqs)
fzpe = f/2.0 # zero point energies of the modes
zpe = fzpe.sum() # zero point energy of the system
EMAX = emax # Absolute maximum energy
dE = emax/m
E0 = E0 + zpe
R0 = int(divmod(E0,dE)[0])
t[R0] = 1.0
ta[R0] = 1.0
for i in range(m+1):
e[i] = i*dE
for l,fmode in enumerate(f):
k = 1
Ek = fmode*k
Rk = int(divmod(Ek,dE)[0])
while Ek < emax:
for j in range(m-Rk+1):
ta[Rk+j] = ta[Rk+j] + t[j]
k += 1
Ek = fmode*k
Rk = int(divmod(Ek,dE)[0])
t[:] = ta
for i in range(1,m+1):
ta[i] = ta[i] + ta[i-1]
return e,t,t/dE,ta
[docs]
def classicalSumOfStates(freqs,emax,m,E0=0.0):
"""
Classical approximation to the sum of states for harmonic systems
See: W.H. Miller, JACS 101, 6810 (1979)
freqs: list of vibrational frequencies (hbar * omega)
emax: maximum energy related to the ZPE
m: discretization for the sum
E0: energy offset (can be a reaction barrier if dealing with a TS)
On return four lists containing:
E : binned energy units
N(E) : number of states at energy E (+- dE)
rho(E) : density of states at energy E
S(E) : sum of states up to energy E
"""
from scipy.misc import factorial
E = np.zeros(m+1,float)
f = np.array(freqs)
fzpe = f/2.0 # zero point energies of the modes
zpe = fzpe.sum() # zero point energy of the system
EMAX = emax # Absolute maximum energy
dE = emax/m # discretization
s = len(f)
for i in range(m+1):
E[i] = i*dE
E[0] = 1e-8
sfac = factorial(s,exact=0)
den = sfac*f.prod()
N = (E-E0)**s/den
for i in range(m*1):
if E[i] >= E0: continue
if E[i] < E0: N[i] = 0.0
RhoE = s*N/E
nE = RhoE*dE
return E,nE,RhoE,N
[docs]
def GaussianFunction(x, **option):
"""
Gaussian function
Parameters:
-----------
x : 1D array_like
variable domain.
xbar : float, default: 0.0, optional
mean value
std : float, default: 1.0, optional
standard deviation
fwhm : float, default: None, optional
full width at half maximum. If fwhm is None (or non-positive number),
the std (standard deviation) will be used. If fwhm is present, then it
is used in to calculate std. fwhm = 2.0 * std * sqrt(2.0 * ln 2.0)
norm : bool, default: True, optional
if norm == True, the gaussian function is normalized
else, its pre-factor is 1.0
Returns:
--------
Gaussian function
"""
xbar = option.get('xbar', 0.0)
std = option.get('std', 1.0)
norm = option.get('norm', True)
fwhm = option.get('fwhm', None)
if fwhm != None and fwhm > 0.0:
std = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0)))
A = 1.0
if norm:
A = 1.0 / (std * np.sqrt(2.0 * np.pi))
gaussian = A * np.exp(-0.5 * ((x - xbar)/std)**2)
return gaussian
[docs]
def gaussian_moments(x, y, mode='area'):
"""
First and second moment (mean value and width) and height for Gaussian peak fitting
The return values of this function can be used as a starting point in fit optimization.
Parameters:
x : 1D array_like (x values)
y : 1D array_like (y / data values)
mode : string, default: None, optional
'area' - height is determined such that the area below the peak is preserved
'max' - height is determined to maximum data value
'mean' - height is set to mean data value
Returns:
mean, width, height
"""
mean = np.sum( x * y ) / np.sum ( y )
width = np.sqrt( np.abs( np.sum( ( x - mean )**2 * y ) / np.sum( y )))
if mode=='area':
if x.size==1:
height = np.max(y)
else:
dx = x[1]-x[0]
height = np.sum( dx * y ) / width / np.sqrt( 2.0 * np.pi )
elif mode=='max':
height = np.max( y )
elif mode=='mean':
height = np.mean( y )
return mean, width, height
[docs]
def boltzmann_weights(e_tot, T_K=300):
"""
Calculates the Boltzmann weights of an ensemble
Parameters:
e_tot - array of energies of the ensemble
T_K - ensemble temperature [Kelvin, default: 300]
Returns:
array of Boltzmann weights for the ensemble
"""
e_tot = e_tot - np.min(e_tot)
beta = 1. / ( cv.FPC_K_B_EV * T_K ) # Boltzmann constant [eV/K]
zsum = np.sum( np.exp( - e_tot * beta ))
return np.exp( - e_tot * beta) / zsum
[docs]
def calc_rdf(r_xyz, ix_A, ix_B, dr=0.1, a=15.96455):
'''
Calculate radial distribution function (pair correlation function) for Type A - Type B
from a trajectory snapshot obeying the minimum image convention
Cf. Levine et al, J Comp Phys 230, 3556 (2011)
Normalized to box volume and particle number(s)
Bin centered (like VMD output)
Parameters:
r_xyz - array with xyz atom coordinates (single frame)
ix_A - row indices of atom type A in r_xyz
ix_B - row indices of atom type B in r_xyz
dr - bin size [default: 0.1 an]
a - box length [default: 15.96455 an from CP2K trajectories]
Returns:
r - r values (ordinate)
g - radial distribution function (abscissa)
'''
rmax = a / 2. # minimum image convention PBC
dist = _calc_distances_levine(r_xyz, ix_A, ix_B, a=a)
r = np.arange(0, rmax, dr)
p = np.zeros_like(r)
# To ensure correct normalization, calculate histogram by hand
for i in range(r.size):
for j in range(r.size + 1):
r_kappa = j * dr
if r_kappa <= r[i] < r_kappa + dr:
for k in range(dist.shape[0]):
for l in range(dist.shape[1]):
if k==l:
continue
if r_kappa <= dist[k,l] < r_kappa + dr:
p[i] += 1
# Normalize RDF to volume
V = a**3
g = p / 4.0 / np.pi / dr / (r + dr/2.)**2 * V
# Normalize RDF to particle number
if np.all(ix_A == ix_B):
g = g / ix_A.size / (ix_A.size - 1)
else:
g = g / ix_A.size / ix_B.size
return r + dr/2., g
[docs]
def angle(v1, v2, deg=False):
'''
Return the angle between vectors v1, v2
deg - if True, use degree, else rad
'''
a = np.arccos( np.dot(v1, v2) / np.linalg.norm(v1) / np.linalg.norm(v2) )
if deg:
return a*180./np.pi
return a
##############################################################################
# --- INTERNAL functions ---
##############################################################################
def _calc_distances_levine(r_xyz, ix_A, ix_B, a=15.96455):
'''
Calculate pair-wise distances between atom types A and B
Parameters:
r_xyz - array with xyz atom coordinates (single frame)
ix_A - row indices of atom type A in r_xyz
ix_B - row indices of atom type B in r_xyz
a - box length [default: 15.96455 an from CP2K trajectories]
Returns:
pairwise distances, array-shaped ( len(ix_A), len(ix_B) )
'''
r = np.zeros( [ len(ix_A), len(ix_B), 3 ] )
for j in range(ix_A.size):
for k in range(ix_B.size):
r_jk = np.abs( r_xyz[ix_A[j]] - r_xyz[ix_B[k]] )
for i in range(3):
if r_jk[i] <= a/2.:
r[j,k,i] = r_jk[i]
else:
r[j,k,i] = a - r_jk[i]
return np.sqrt( r[:,:,0]**2 + r[:,:,1]**2 + r[:,:,2]**2 )
def _segment(vect,point):
"""Given a list of ordered values and a value, what is the index
of the point at the left of the given point?
vect -- list of ordered values
point -- value
on return, an index (integer). The returned index is always
at least 0 and at maximum len(vect)-2
"""
if point < vect[0]: return 0
if point > vect[len(vect)-1]: return len(vect)-2
for i in range(len(vect)-1):
if point > vect[i] and point < vect[i+1]:
return i
return 0
def _toPeriodic(xg,yg,order=4):
"""Modify xg and yg vectors for a periodic interpolation of order "order"
xg -- list of x values
yg -- list of y(x) values
order -- maximum order of the interpolation
"""
if len(xg) < order:
error = "error in periodic interpolation, too few points"
raise ValueError(error)
xgnew = []
ygnew = []
xlast = xg[len(xg)-1]
dx = xg[1]-xg[0]
for i in range(len(xg)):
xgnew.append(xg[i])
ygnew.append(yg[i])
for i in range(order):
xgnew.append(xlast+dx*(i+1))
xgnew.insert(0,xg[0]-dx*(i+1))
ygnew.append(yg[i])
ygnew.insert(0,yg[len(yg)-1-i])
return xgnew,ygnew
[docs]
def lag(x, n):
"""
Compute the hypergeometric Laguerre polynomial Ln(x)
"""
# Lm(x) [0]-[n] ord.
poly = np.zeros(n+1, dtype=float)
poly[0] = 1
poly[1] = 1 - x
m = 2
while m < (n + 1):
mf = float(m)
poly[m] = ((2*mf-1-x) * poly[m-1] - (mf-1) * poly[m-2]) / mf
m += 1
return poly[n]
[docs]
def distributionFromList(l,fwhm,nstep=100,lbound=None,hbound=None):
"""
Generate a distribution from a 1D list of values using a Gaussian kernel
l - 1D array with data
fwhm - full width at half maximum of the Gaussian kernel
nstep - number of discrete output points
lbound - lower bound of output;
hbound - higher bound of output
"""
lrange = l.max() - l.min()
if lbound is None:
lbound = l.min() - 0.1*lrange
if hbound is None:
hbound = l.max() + 0.1*lrange
x = np.linspace(lbound,hbound,nstep)
y = np.zeros(len(x),float)
for i,xp in enumerate(x):
for xl in l:
y[i] = y[i] + _gaussian(xl-xp,fwhm)
return x,y
__k = 4.0*math.log(2.0)
__N = 4.0*math.log(2.0)/math.pi
def _gaussian(x,fwhm=1.0,norm=True):
N = __N/fwhm
if norm: return N*np.exp(-__k*((x)/fwhm)**2)
else: return np.exp(-__k*((x)/fwhm)**2)