CDTK.Tools.Mathematics module

SECTIONS

• LINEAR ALGEBRA

• TENSOR ALGEBRA

• VECTORS in 3D SPACE

• INTERPOLATION

• NUMERICAL

• MONTE-CARLO

• OTHER

• internal functions

CDTK.Tools.Mathematics.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, g1D array_like

function to be convolved.

dtfloat

interval of the co-domain (convolved) function.

methodstring, default: ‘simp’, optional

integration method to be used. ‘trapz’ for trapezoidal rule, ‘simp’ for Simpson rule

npointsint, default: max(len(f), len(g))

length of output function.

ksplineint {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 – convolution integral of f and g function.

rtype:

1D numpy.ndarray

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.FourierTransform(y, dt, **option)

Compute Fourier transform

Parameters:

y1D array_like

function to be Fourier transformed.

dtfloat

interval of the domain function.

methodstring, default: ‘simp’, optional

integration method to be used. ‘trapz’ for trapezoidal rule, ‘simp’ for Simpson rule

t_initfloat, default: 0.0, optional

initial point of the domain function.

windowstring, default: None, optional

Window function or apodization or tapering function to minimize the leakage of Fourier transformation. see https://en.wikipedia.org/wiki/Window_function for a quick reference. Currently, the window functions implemented here are: Blackman function: ‘blackman’ see numpy.blackman Bartlett function: ‘bartlett’ see numpy.bartlett Hamming function: ‘hamming’ see numpy.hamming Hann function: ‘hanning’ or ‘hann’ see numpy.hanning

usage{‘phys’, ‘math’}, default: ‘phys’, optional

usage of the Fourier transform: ‘phys’: mostly for physics application, it means using angular frequency, (2pi*f),

thus the inverse is asymmetric.

‘math’: mostly for mathematics application, it means using normal frequency,

thus the inverse is symmetric.

f_rangetuple of float, default: if usage=’phys’ (-pi/dt, pi/dt), if usage=’math’ (-0.5/dt, 0.5/dt), optional

minimum frequency and maximum frequency.

npointsint, default: len(y), optional

number of points of transformed domain function.

symmetry{‘gen’, ‘even’, ‘odd’}, default: ‘gen’, optional

symmetry of the function y to be transformed. ‘gen’: general, the exponential function will be used. ‘even’: even function, cosine transform will be used. ‘odd’: odd function, sine transform will be used. if symmetry is ‘even’ or ‘odd’, the t_init will be set to 0.0

returns:

• fline, FT_func ((list))

• fline (1D numpy.ndarray) – Domain of transformed function (frequency domain)

• FT_func (1D numpy.ndarray) – Fourier transform of y function.

CDTK.Tools.Mathematics.GaussianFunction(x, **option)

Gaussian function

Parameters:

x1D array_like

variable domain.

xbarfloat, default: 0.0, optional

mean value

stdfloat, default: 1.0, optional

standard deviation

fwhmfloat, 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)

normbool, default: True, optional

if norm == True, the gaussian function is normalized else, its pre-factor is 1.0

Returns:

Gaussian function

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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)

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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.

CDTK.Tools.Mathematics.angle(v1, v2, deg=False)

Return the angle between vectors v1, v2

deg - if True, use degree, else rad

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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)

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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|

CDTK.Tools.Mathematics.distributeMatrixList(ml, oper)

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.gramSchmidt(vecs)

Perform a Gram-Schmidt orthonormalization on a given list of arrays

vecs – list of arrays to orthonormalize

CDTK.Tools.Mathematics.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))

CDTK.Tools.Mathematics.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)

CDTK.Tools.Mathematics.ijk_shape(f, oshape)

change tuple oshape representing a shape to a 3 index tuple collecting indices before and after index f

CDTK.Tools.Mathematics.invert_reg(mat, tol=1e-08)

Invert matrix mat, regularizing it in case it has eigenvalues close to 0

CDTK.Tools.Mathematics.isPositiveDefinite(m)

Return 1 if m is positive definite, 0 if it’s not

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.lag(x, n)

Compute the hypergeometric Laguerre polynomial Ln(x)

CDTK.Tools.Mathematics.leviCivita(i, j, k)

Return the value of the antisymmetric permutation tensor for indices i,j,k

i,j,k - ints

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.makeMesh(f, grid, ctype=<class '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.

CDTK.Tools.Mathematics.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)

CDTK.Tools.Mathematics.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)

CDTK.Tools.Mathematics.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.

CDTK.Tools.Mathematics.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.

CDTK.Tools.Mathematics.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),…)

CDTK.Tools.Mathematics.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]

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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)

CDTK.Tools.Mathematics.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)

CDTK.Tools.Mathematics.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)

CDTK.Tools.Mathematics.pointsDistance(p, q)

Return the distance (Euclidean norm of connecting vector) between two points p - vector or list q - vector or list

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.rmsError(list1, list2, energy=30.0)

CDTK.Tools.Mathematics.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|

CDTK.Tools.Mathematics.scattered2grid2D(points, xgrid, ygrid, kernel, **args)

Interpolate scattered 2D data into a regular grid

Parameters:

w (kernel --) – z = Sum_i w(d_i)

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.sortEigValsVecs(evals, evecs)

Return the pair (eigvals,eigvecs) sorted from low to high values

CDTK.Tools.Mathematics.sphericalAngles(v)

return the spherical angles theta(polar) and phi(azimutal) of v

v - 3D vector (x,y,z)

CDTK.Tools.Mathematics.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.

CDTK.Tools.Mathematics.uniformPointsWithinSphere(center, n=100, r=1.0)

Draw a uniform distribution of points within a sphere around center

Draw uniform vector within cube containing sphere, and discard if not within sphere

Parameters: center - center of the sphere n - number of points [default: 100] r - sphere radius [default: 1.0]

Return: Coordinates of n points

CDTK.Tools.Mathematics.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)

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.wignerE(q, p, m, w, hbar=1.0)

CDTK.Tools.Mathematics.wignerEFT(q, p, m, w, temp, hbar=1.0)

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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

CDTK.Tools.Mathematics.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