#
#@BEGIN LICENSE
#
# PSI4: an ab initio quantum chemistry software package
#
# This program 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 2 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, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
#@END LICENSE
#
r"""File for accessory procedures in the chem module.
Credit for the libmints vector3 class to Justin M. Turney and
incremental improvements by other psi4 developers.
"""
from __future__ import absolute_import
from __future__ import print_function
import math
import copy
from .exceptions import *
ZERO = 1.0E-14
[docs]def norm(v):
"""Compute the magnitude of vector *v*."""
return math.sqrt(sum(v[i] * v[i] for i in range(len(v))))
[docs]def add(v, u):
"""Compute sum of vectors *v* and *u*."""
return [u[i] + v[i] for i in range(len(v))]
[docs]def sub(v, u):
"""Compute difference of vectors *v* - *u*."""
return [v[i] - u[i] for i in range(len(v))]
[docs]def dot(v, u):
"""Compute dot product of vectors *v* and *u*."""
return sum(u[i] * v[i] for i in range(len(v)))
[docs]def scale(v, d):
"""Compute by-element scale by *d* of vector *v*."""
return [d * v[i] for i in range(len(v))]
[docs]def naivemult(v, u):
"""Compute by-element multiplication of vectors *v* and *u*."""
if len(u) != len(v):
raise ValidationError('naivemult() only defined for vectors of same length \n')
return [u[i] * v[i] for i in range(len(v))]
[docs]def normalize(v):
"""Compute normalized vector *v*."""
vmag = norm(v)
return [v[i] / vmag for i in range(len(v))]
[docs]def distance(v, u):
"""Compute the distance between points defined by vectors *v* and *u*."""
return norm(sub(v, u))
[docs]def cross(v, u):
"""Compute cross product of length 3 vectors *v* and *u*."""
if len(u) != 3 or len(v) != 3:
raise ValidationError('cross() only defined for vectors of length 3\n')
return [v[1] * u[2] - v[2] * u[1],
v[2] * u[0] - v[0] * u[2],
v[0] * u[1] - v[1] * u[0]]
[docs]def rotate(v, theta, axis):
"""Rotate length 3 vector *v* about *axis* by *theta* radians."""
if len(v) != 3 or len(axis) != 3:
raise ValidationError('rotate() only defined for vectors of length 3\n')
unitaxis = normalize(copy.deepcopy(axis))
# split into parallel and perpendicular components along axis
parallel = scale(axis, dot(v, axis) / dot(axis, axis))
perpendicular = sub(v, parallel)
# form unit vector perpendicular to parallel and perpendicular
third_axis = perp_unit(axis, perpendicular)
third_axis = scale(third_axis, norm(perpendicular))
result = add(parallel, add(scale(perpendicular, math.cos(theta)), scale(third_axis, math.sin(theta))))
for item in range(len(result)):
if math.fabs(result[item]) < ZERO:
result[item] = 0.0
return result
[docs]def perp_unit(u, v):
"""Compute unit vector perpendicular to length 3 vectors *u* and *v*."""
if len(u) != 3 or len(v) != 3:
raise ValidationError('perp_unit() only defined for vectors of length 3\n')
# try cross product
result = cross(u, v)
resultdotresult = dot(result, result)
if resultdotresult < 1.E-16:
# cross product is too small to normalize
# find the largest of this and v
dotprodt = dot(u, u)
dotprodv = dot(v, v)
if dotprodt < dotprodv:
d = copy.deepcopy(v)
dotprodd = dotprodv
else:
d = copy.deepcopy(u)
dotprodd = dotprodt
# see if d is big enough
if dotprodd < 1.e-16:
# choose an arbitrary vector, since the biggest vector is small
result = [1.0, 0.0, 0.0]
return result
else:
# choose a vector perpendicular to d
# choose it in one of the planes xy, xz, yz
# choose the plane to be that which contains the two largest components of d
absd = [math.fabs(d[0]), math.fabs(d[1]), math.fabs(d[2])]
if (absd[1] - absd[0]) > 1.0e-12:
#if absd[0] < absd[1]:
axis0 = 1
if (absd[2] - absd[0]) > 1.0e-12:
#if absd[0] < absd[2]:
axis1 = 2
else:
axis1 = 0
else:
axis0 = 0
if (absd[2] - absd[1]) > 1.0e-12:
#if absd[1] < absd[2]:
axis1 = 2
else:
axis1 = 1
result = [0.0, 0.0, 0.0]
# do the pi/2 rotation in the plane
result[axis0] = d[axis1]
result[axis1] = -1.0 * d[axis0]
result = normalize(result)
return result
else:
# normalize the cross product and return the result
result = scale(result, 1.0 / math.sqrt(resultdotresult))
return result
[docs]def determinant(mat):
"""Given 3x3 matrix *mat*, compute the determinat
"""
if len(mat) != 3 or len(mat[0]) != 3 or len(mat[1]) != 3 or len(mat[2]) != 3:
raise ValidationError('determinant() only defined for arrays of dimension 3x3\n')
det = mat[0][0] * mat[1][1] * mat[2][2] - mat[0][2] * mat[1][1] * mat[2][0] + \
mat[0][1] * mat[1][2] * mat[2][0] - mat[0][1] * mat[1][0] * mat[2][2] + \
mat[0][2] * mat[1][0] * mat[2][1] - mat[0][0] * mat[1][2] * mat[2][1]
return det
[docs]def diagonalize3x3symmat(M):
"""Given an real symmetric 3x3 matrix *M*, compute the eigenvalues
"""
if len(M) != 3 or len(M[0]) != 3 or len(M[1]) != 3 or len(M[2]) != 3:
raise ValidationError('diagonalize3x3symmat() only defined for arrays of dimension 3x3\n')
A = copy.deepcopy(M) # Symmetric input matrix
Q = [[1, 0, 0], [0, 1, 0], [0, 0, 1]] # Storage buffer for eigenvectors
w = [A[0][0], A[1][1], A[2][2]] # Storage buffer for eigenvalues
# sd, so # Sums of diagonal resp. off-diagonal elements
# s, c, t # sin(phi), cos(phi), tan(phi) and temporary storage
# g, h, z, theta # More temporary storage
# Calculate SQR(tr(A))
sd = 0.0
for i in range(3):
sd += math.fabs(w[i])
sd = sd * sd
# Main iteration loop
for nIter in range(50):
# Test for convergence
so = 0.0
for p in range(3):
for q in range(p + 1, 3):
so += math.fabs(A[p][q])
if so == 0.0:
return w, Q # return eval, evec
if nIter < 4:
thresh = 0.2 * so / (3 * 3)
else:
thresh = 0.0
# Do sweep
for p in range(3):
for q in range(p + 1, 3):
g = 100.0 * math.fabs(A[p][q])
if nIter > 4 and (math.fabs(w[p]) + g == math.fabs(w[p])) and \
(math.fabs(w[q]) + g == math.fabs(w[q])):
A[p][q] = 0.0
elif math.fabs(A[p][q]) > thresh:
# Calculate Jacobi transformation
h = w[q] - w[p]
if math.fabs(h) + g == math.fabs(h):
t = A[p][q] / h
else:
theta = 0.5 * h / A[p][q]
if theta < 0.0:
t = -1.0 / (math.sqrt(1.0 + theta * theta) - theta)
else:
t = 1.0 / (math.sqrt(1.0 + theta * theta) + theta)
c = 1.0 / math.sqrt(1.0 + t * t)
s = t * c
z = t * A[p][q]
# Apply Jacobi transformation
A[p][q] = 0.0
w[p] -= z
w[q] += z
for r in range(p):
t = A[r][p]
A[r][p] = c * t - s * A[r][q]
A[r][q] = s * t + c * A[r][q]
for r in range(p + 1, q):
t = A[p][r]
A[p][r] = c * t - s * A[r][q]
A[r][q] = s * t + c * A[r][q]
for r in range(q + 1, 3):
t = A[p][r]
A[p][r] = c * t - s * A[q][r]
A[q][r] = s * t + c * A[q][r]
# Update eigenvectors
for r in range(3):
t = Q[r][p]
Q[r][p] = c * t - s * Q[r][q]
Q[r][q] = s * t + c * Q[r][q]
return None
[docs]def zero(m, n):
""" Create zero matrix"""
new_matrix = [[0 for row in range(n)] for col in range(m)]
return new_matrix
[docs]def identity(m):
"""Create identity matrix"""
new_matrix = zero(m, m)
for i in range(m):
new_matrix[i][i] = 1.0
return new_matrix
[docs]def show(matrix):
""" Print out matrix"""
for col in matrix:
print(col)
[docs]def mscale(matrix, d):
"""Return *matrix* scaled by scalar *d*"""
for i in range(len(matrix)):
for j in range(len(matrix[0])):
matrix[i][j] *= d
return matrix
[docs]def mult(matrix1, matrix2):
""" Matrix multiplication"""
if len(matrix1[0]) != len(matrix2):
# Check matrix dimensions
raise ValidationError('Matrices must be m*n and n*p to multiply!')
else:
# Multiply if correct dimensions
try:
new_matrix = zero(len(matrix1), len(matrix2[0]))
for i in range(len(matrix1)):
for j in range(len(matrix2[0])):
for k in range(len(matrix2)):
new_matrix[i][j] += matrix1[i][k] * matrix2[k][j]
except TypeError:
new_matrix = zero(len(matrix1), 1)
for i in range(len(matrix1)):
for k in range(len(matrix2)):
new_matrix[i][0] += matrix1[i][k] * matrix2[k]
return new_matrix
[docs]def transpose(matrix):
"""Return matrix transpose"""
if len(matrix[0]) != len(matrix):
# Check matrix dimensions
raise ValidationError('Matrices must be square.')
tmat = [list(i) for i in zip(*matrix)]
return tmat
[docs]def matadd(matrix1, matrix2, fac1=1.0, fac2=1.0):
"""Matrix addition"""
if (len(matrix1[0]) != len(matrix2[0])) or (len(matrix1) != len(matrix2)):
raise ValidationError('Matrices must be same dimension to add.')
new_matrix = zero(len(matrix1), len(matrix1[0]))
for i in range(len(matrix1)):
for j in range(len(matrix1[0])):
new_matrix[i][j] = fac1 * matrix1[i][j] + fac2 * matrix2[i][j]
return new_matrix