import sys
import os
from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
import numpy as np
from numpy import linalg as la
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from tkcrystalbase import *
# 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])]
]
# Coefficient for atomic size to plot
kr = 100.0
# Coefficient to plot reciprocal unit cell w.r.t. real space unit cell
kRUC = 0.8
# Range of unit cells to draw crystal structure
nrange = [[-0.1, 1.1], [-0.1, 1.1], [-0.1, 1.1]]
# Figure configuration
figsize = (12, 12)
[ドキュメント]
def draw_box(ax, aij, nrange, color = 'black'):
# (0,0,0) -> ax
ax.plot([0.0, aij[0][0]],
[0.0, aij[0][1]],
[0.0, aij[0][2]], color = color)
# (0,0,0) -> ay
ax.plot([0.0, aij[1][0]],
[0.0, aij[1][1]],
[0.0, aij[1][2]], color = color)
# (0,0,0) -> az
ax.plot([0.0, aij[2][0]],
[0.0, aij[2][1]],
[0.0, aij[2][2]], color = color)
# ax -> ax + ay
ax.plot([aij[0][0], aij[0][0] + aij[1][0]],
[aij[0][1], aij[0][1] + aij[1][1]],
[aij[0][2], aij[0][2] + aij[1][2]], color = color)
# ax -> ax + az
ax.plot([aij[0][0], aij[0][0] + aij[2][0]],
[aij[0][1], aij[0][1] + aij[2][1]],
[aij[0][2], aij[0][2] + aij[2][2]], color = color)
# ay -> ay + ax
ax.plot([aij[1][0], aij[1][0] + aij[0][0]],
[aij[1][1], aij[1][1] + aij[0][1]],
[aij[1][2], aij[1][2] + aij[0][2]], color = color)
# ay -> ay + az
ax.plot([aij[1][0], aij[1][0] + aij[2][0]],
[aij[1][1], aij[1][1] + aij[2][1]],
[aij[1][2], aij[1][2] + aij[2][2]], color = color)
# az -> az + ax
ax.plot([aij[2][0], aij[2][0] + aij[0][0]],
[aij[2][1], aij[2][1] + aij[0][1]],
[aij[2][2], aij[2][2] + aij[0][2]], color = color)
# az -> ax + ay
ax.plot([aij[2][0], aij[2][0] + aij[1][0]],
[aij[2][1], aij[2][1] + aij[1][1]],
[aij[2][2], aij[2][2] + aij[1][2]], color = color)
# ax + ay -> ax + ay + az
ax.plot([aij[0][0] + aij[1][0], aij[0][0] + aij[1][0] + aij[2][0]],
[aij[0][1] + aij[1][1], aij[0][1] + aij[1][1] + aij[2][1]],
[aij[0][2] + aij[1][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
# ax + az -> ax + ay + az
ax.plot([aij[0][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]],
[aij[0][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]],
[aij[0][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
# ay + az -> ax + ay + az
ax.plot([aij[1][0] + aij[2][0], aij[0][0] + aij[1][0] + aij[2][0]],
[aij[1][1] + aij[2][1], aij[0][1] + aij[1][1] + aij[2][1]],
[aij[1][2] + aij[2][2], aij[0][2] + aij[1][2] + aij[2][2]], color = color)
[ドキュメント]
def draw_unitcell(ax, sites, aij, nrange, color = 'black'):
draw_box(ax, aij, nrange, color)
if sites is None:
return
for site in sites:
name, label, z, M, q, r, color, pos = site
pos01 = [reduce01(pos[0]), reduce01(pos[1]), reduce01(pos[2])]
for iz in range(int(nrange[2][0]) - 1, int(nrange[2][1]) + 1):
for iy in range(int(nrange[1][0]) - 1, int(nrange[1][1]) + 1):
for ix in range(int(nrange[0][0]) - 1, int(nrange[0][1]) + 1):
posn = [pos01[0] + ix, pos01[1] + iy, pos01[2] + iz]
if posn[0] < nrange[0][0] or nrange[0][1] < posn[0] \
or posn[1] < nrange[1][0] or nrange[1][1] < posn[1] \
or posn[2] < nrange[2][0] or nrange[2][1] < posn[2]:
continue
x, y, z = fractional_to_cartesian(posn, aij)
ax.scatter([x], [y], [z], marker = 'o', c = color, s = kr *r)
[ドキュメント]
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]))
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))
fig = plt.figure(figsize = figsize)
ax = fig.add_subplot(111, projection='3d')
# Real space unit cell
draw_unitcell(ax, sites, aij, nrange)
# Reciprocal space unit cell
k = max([*aij[0], *aij[1], *aij[2]]) / max([*Raij[0], *Raij[1], *Raij[2]]) * kRUC
kRaij = np.empty([3, 3])
for i in range(3):
for j in range(3):
kRaij[i][j] = k * Raij[i][j]
draw_unitcell(ax, None, kRaij, nrange, color = 'red')
# Note: set_aspect() is not implemented for 3D plots
# ax.set_aspect('equal','box')
xlim =ax.get_xlim()
ylim =ax.get_ylim()
zlim =ax.get_zlim()
lim = [min([xlim[0], ylim[0], zlim[0]]), max([xlim[1], ylim[1], zlim[1]])]
ax.set_xlim(lim)
ax.set_ylim(lim)
ax.set_zlim(lim)
# ax.set_xticks(np.linspace(*lim, 0))
# ax.set_yticks(np.linspace(*lim, 0))
# ax.set_zticks(np.linspace(*lim, 0))
plt.show()
print("")
exit()
if __name__ == '__main__':
main()