Base library for calculating crystal properties
Download script from .\tkcrystalbase.py
from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
import numpy as np
from numpy import linalg as la
from pprint import pprint
"""
Base library for calculating crystal properties
"""
pi = 3.14159265358979323846
pi2 = pi + pi
torad = 0.01745329251944 # rad/deg";
todeg = 57.29577951472 # deg/rad";
basee = 2.71828183
h = 6.6260755e-34 # Js";
h_bar = 1.05459e-34 # "Js";
hbar = h_bar
c = 2.99792458e8 # m/s";
e = 1.60218e-19 # C";
me = 9.1093897e-31 # kg";
mp = 1.6726231e-27 # kg";
mn = 1.67495e-27 # kg";
u0 = 4.0 * 3.14*1e-7; # . "Ns2C-2";
e0 = 8.854418782e-12; # C2N-1m-2";
e2_4pie0 = 2.30711e-28 # Nm2";
a0 = 5.29177e-11 # m";
kB = 1.380658e-23 # JK-1";
NA = 6.0221367e23 # mol-1";
R = 8.31451 # J/K/mol";
F = 96485.3 # C/mol";
g = 9.81 # m/s2";
# Lattice parameters (angstrom and degree)
#lattice_parameters = [ 5.62, 5.62, 5.62, 60.0, 60.0, 60.0]
lattice_parameters = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]
# Site information (atom name, site label, atomic number, atomic mass, charge, radius, color, position)
sites = [
['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])]
,['Na', 'Na2', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.5, 0.5])]
,['Na', 'Na3', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.0, 0.5])]
,['Na', 'Na4', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.5, 0.5, 0.0])]
,['Cl', 'Cl1', 17, 35.4527, +1.0, 1.4, 'blue', np.array([0.5, 0.0, 0.0])]
,['Cl', 'Cl2', 17, 35.4527, +1.0, 1.4, 'blue', np.array([0.5, 0.5, 0.5])]
,['Cl', 'Cl3', 17, 35.4527, +1.0, 1.4, 'blue', np.array([0.0, 0.0, 0.5])]
,['Cl', 'Cl4', 17, 35.4527, +1.0, 1.4, 'blue', np.array([0.0, 0.5, 0.0])]
]
def reduce01(x):
return x - int(x)
return x
def round01(x):
if abs(x - 1.0) < 0.0002:
return 1.0
if abs(x) < 0.0002:
return 0.0
if x > 1.0:
return x - int(x)
if x < 1.0:
return x - int(x) + 1.0
return x
def round(x, xmin, xmax):
if x > xmax:
return x - int(x - xmin)
if x < xmin:
return x - int(x - xmin) + 1.0
return x
def round_parameter(x, tol):
val = tol * int( (x+0.1*tol) / tol )
return val
def cal_lattice_vectors(latt):
cosb = cos(torad * latt[4])
cosg = cos(torad * latt[5])
sing = sin(torad * latt[5])
aij = np.empty([3, 3], dtype = float)
aij[0][0] = latt[0]
aij[0][1] = 0.0;
aij[0][2] = 0.0;
aij[1][0] = latt[1] * cosg
aij[1][1] = latt[1] * sing
aij[1][2] = 0.0;
aij[2][0] = latt[2] * cosb
aij[2][1] = (latt[2] * cosb - aij[2][0] * cosg)
if abs(aij[2][1]) < 1.0e-8:
aij[2][1] = 0.0
else:
aij[2][1] = aij[2][1] / sing
aij[2][2] = sqrt(latt[2] * latt[2] - aij[2][0] * aij[2][0] - aij[2][1] * aij[2][1])
return aij
def cal_metrics(latt):
inf = {}
cosa = cos(torad * latt[3])
cosb = cos(torad * latt[4])
cosg = cos(torad * latt[5])
aij = cal_lattice_vectors(latt)
inf['aij'] = aij
gij = np.empty([3, 3], dtype = float)
for i in range(3):
for j in range(i, 3):
gij[i][j] = np.dot(aij[i], aij[j])
gij[j][i] = gij[i][j]
inf['gij'] = gij
return inf
def cal_volume(aij):
axb = np.cross(aij[0], aij[1]) # Outner product
vol = np.dot(axb, aij[2]) # Inner product
return vol
def cal_reciprocal_lattice_vectors(aij):
V = cal_volume(aij)
Ra = np.cross(aij[1], aij[2]) / V
Rb = np.cross(aij[2], aij[0]) / V
Rc = np.cross(aij[0], aij[1]) / V
return [Ra, Rb, Rc]
def cal_reciprocal_lattice_parameters(Raij):
Ra = la.norm(Raij[0])
Rb = la.norm(Raij[1])
Rc = la.norm(Raij[2])
Ralpha = todeg * arccos(np.dot(Raij[1], Raij[2]) / Rb / Rc)
Rbeta = todeg * arccos(np.dot(Raij[2], Raij[0]) / Rc / Ra)
Rgamma = todeg * arccos(np.dot(Raij[0], Raij[1]) / Ra / Rb)
return [Ra, Rb, Rc, Ralpha, Rbeta, Rgamma]
def get_conversion_matrix_from_tij(tij, mode = 'x2xc'):
if mode == 'a2ac':
return tij
if mode == 'x2xc':
return np.linalg.inv(tij).transpose()
if mode == 'xc2x':
return tij.transpose()
if mode == 'h2hc':
return tij
if mode == 'a*2a*c':
return np.linalg.inv(tij).transpose()
if mode == 'ac*2a*':
return tij.transpose()
print("")
print("Error in tkcrystalbase.get_conversion_matrix: Invalid mode [{}]".format(mode))
print("")
exit()
def convert_lattice_vectors(aij, tij):
return tij @ aij
def cal_lattice_parameters_from_lattice_vectors(aij):
a = np.linalg.norm(aij[0])
b = np.linalg.norm(aij[1])
c = np.linalg.norm(aij[2])
alpha = np.dot(aij[1], aij[2]) / b / c * todeg
beta = np.dot(aij[2], aij[0]) / c / a * todeg
gamma = np.dot(aij[0], aij[1]) / a / b * todeg
return a, b, c, alpha, beta, gamma
def convert_fractional_coordinates_by_tij(pos, tij):
tijc = get_conversion_matrix_from_tij(tij, 'x2xc')
return tijc @ pos
def fractional_to_cartesian(pos, aij):
x = pos[0] * aij[0][0] + pos[1] * aij[1][0] + pos[2] * aij[2][0]
y = pos[0] * aij[0][1] + pos[1] * aij[1][1] + pos[2] * aij[2][1]
z = pos[0] * aij[0][2] + pos[1] * aij[1][2] + pos[2] * aij[2][2]
return x, y, z
def distance2(pos0, pos1, gij):
# dx = pos1 - pos0
dx = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
r2 = gij[0][0] * dx[0]*dx[0] + gij[1][1] * dx[1]*dx[1] + gij[2][2] * dx[2]*dx[2] \
+ 2.0 * (gij[0][1] * dx[0]*dx[1] + gij[0][2] * dx[0]*dx[2] + gij[1][2] * dx[1]*dx[2])
return r2
def distance(pos0, pos1, gij):
dx = pos1 - pos0
# dx = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
r2 = gij[0][0] * dx[0]*dx[0] + gij[1][1] * dx[1]*dx[1] + gij[2][2] * dx[2]*dx[2] \
+ 2.0 * (gij[0][1] * dx[0]*dx[1] + gij[0][2] * dx[0]*dx[2] + gij[1][2] * dx[1]*dx[2])
r = sqrt(r2)
return r
def add_site(sites, newsite, gij, rmin = 0.1):
rmin2 = rmin * rmin
name, label, z, M, q, r, color, pos = newsite
for s in sites:
name2, label2, z2, M2, q2, r2, color2, pos2 = s
dis2 = distance2(pos2, pos, gij)
if dis2 < rmin2:
return False
sites.append(newsite)
return True
def angle(pos0, pos1, pos2, gij):
dis01 = distance(pos0, pos1, gij)
if dis01 == 0.0:
return 0.0
dis02 = distance(pos0, pos2, gij)
if dis02 == 0.0:
return 0.0
dx01 = pos1 - pos0
dx02 = pos2 - pos0
# dx01 = [pos1[0] - pos0[0], pos1[1] - pos0[1], pos1[2] - pos0[2]]
# dx02 = [pos2[0] - pos0[0], pos2[1] - pos0[1], pos2[2] - pos0[2]]
ip = gij[0][0] * dx01[0]*dx02[0] + gij[1][1] * dx01[1]*dx02[1] + gij[2][2] * dx01[2]*dx02[2] \
+ 2.0 * (gij[0][1] * dx01[0]*dx02[1] + gij[0][2] * dx01[0]*dx02[2] + gij[1][2] * dx01[1]*dx02[2])
cosa = ip / dis01 / dis02
angle = todeg * arccos(cosa)
if angle > 180.0:
angle = 360.0 - angle
return angle;
def main():
print("")
print("Lattice parameters:", lattice_parameters)
aij = cal_lattice_vectors(lattice_parameters)
print("Lattice vectors:")
print(" ax: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[0][0], aij[0][1], aij[0][2]))
print(" ay: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[1][0], aij[1][1], aij[1][2]))
print(" az: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(aij[2][0], aij[2][1], aij[2][2]))
inf = cal_metrics(lattice_parameters)
gij = inf['gij']
print("Metric tensor:")
print(" gij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[0]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[1]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gij[2]))
print("")
volume = cal_volume(aij)
print("Volume: {:12.4g} A^3".format(volume))
print("")
print("Unit cell volume: {:12.4g} A^3".format(volume))
Raij = cal_reciprocal_lattice_vectors(aij)
Rlatt = cal_reciprocal_lattice_parameters(Raij)
Rinf = cal_metrics(Rlatt)
Rgij = Rinf['gij']
print("Reciprocal lattice parameters:", Rlatt)
print("Reciprocal lattice vectors:")
print(" Rax: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[0]))
print(" Ray: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[1]))
print(" Raz: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Raij[2]))
print("Reciprocal lattice metric tensor:")
print(" Rgij: ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[0]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[1]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A^-1".format(*Rgij[2]))
Rvolume = cal_volume(Raij)
print("Reciprocal unit cell volume: {:12.4g} A^-3".format(Rvolume))
print("")
print("dis=", distance(np.array([0,0,0]), np.array([1,1,1]), gij))
print("angle=", angle (np.array([0,0,0]), np.array([1,1,1]), np.array([1,0,0]), gij))
print("angle=", angle (np.array([0,0,0]), np.array([1,0,0]), np.array([0,1,0]), gij))
print("")
exit()
if __name__ == '__main__':
main()
