import sys
import os
import copy
from numpy import sin, cos, tan, arcsin, arccos, arctan, exp, log, sqrt
import numpy as np
from numpy import linalg as la
from tkcrystalbase import *
"""
Convert unit cell to different base vectors and draw those unit cells
Requirement: tkcrystalbase.py
"""
# Crystal name
crystal_name = 'FCC'
# Conversion mode
conversion_mode = 'FCCPrim'
# 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)
crystal_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])]
]
# Coefficient for atomic size to plot
kr = 1000.0
# Distance to judge identical atom site, in angstrom
rmin = 0.1
# 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)
# Atom color for converted cell
converted_atom_color = 'gray'
[ドキュメント]
def get_crystal(name):
if name == 'FCC':
return lattice_parameters, crystal_sites
if name == 'BCC':
return lattice_parameters, \
[ ['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.5, 0.5, 0.5])] ]
if name == 'Hex':
return [ 5.62, 5.62, 4.00, 90.0, 90.0, 120.0], \
[ ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])] ]
if name == 'Rhomb':
return [ 5.62, 5.62, 5.62, 70.0, 70.0, 70.0], \
[ ['Na', 'Na1', 11, 22.98997, +1.0, 0.7, 'red', np.array([0.0, 0.0, 0.0])] ]
[ドキュメント]
def get_conversion_matrix(key):
if key == 'FCCPrim':
return np.array([ [ 0.5, 0.5, 0],
[ 0, 0.5, 0.5],
[ 0.5, 0, 0.5] ])
elif key == 'BCCPrim':
return np.array([ [-0.5, 0.5, 0.5],
[ 0.5,-0.5, 0.5],
[ 0.5, 0.5,-0.5] ])
elif key == 'ACenterPrim':
return np.array([ [ 1.0, 0.0, 0.0],
[ 0.0, 0.5, 0.5],
[ 0.0,-0.5, 0.5] ])
elif key == 'BCenterPrim':
return np.array([ [ 0.5, 0.0, 0.5],
[ 0.0, 1.0, 0.0],
[-0.5, 0.0, 0.5] ])
elif key == 'CCenterPrim':
return np.array([ [ 0.5, 0.5, 0.0],
[-0.5, 0.5, 0.0],
[ 0.0, 0.0, 1.0] ])
elif key == 'RhombHex':
return np.array([ [ 1.0, -1.0, 0.0],
[ 0.0, 1.0, -1.0],
[ 1.0, 1.0, 1.0] ])
elif key == 'HexRhomb':
return np.array([ [ 2.0/3, 1.0/3, 1.0/3],
[-1.0/3, 1.0/3, 1.0/3],
[-1.0/3,-2.0/3, 1.0/3] ])
elif key == 'HexOrtho':
return np.array([ [ 1.0, 0.0, 0.0],
[ 1.0, 2.0, 0.0],
[ 0.0, 0.0, 1.0] ])
# Treat arguments
if __name__ == "__main__":
argv = sys.argv
narg = len(argv)
if narg >= 2:
crystal_name = argv[1]
if narg >= 3:
conversion_mode = argv[2]
if narg >= 4:
kr = float(argv[3])
[ドキュメント]
def usage():
print("")
print("Usage: python {} crystal_name conversion_mode kRatom".format(argv[0]))
print(" crystal_name : 'FCC', 'BCC', 'Hex', 'Rhomb'")
print(" conversion_mode: 'FCCPrim', 'BCCPrim', 'HexOrtho', 'HexRhomb', 'RhombHex'")
print(" kRatom : Coefficient of atomic radius to draw")
print("ex: python {} {} {} {}".format(argv[0], crystal_name, conversion_mode, kr))
print("")
[ドキュメント]
def terminate():
usage()
exit()
[ドキュメント]
def main():
latt, sites = get_crystal(crystal_name)
tij = get_conversion_matrix(conversion_mode)
print("")
print("Crystal name:", crystal_name)
print(" Lattice parameters: {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g}".format(*latt))
aij = cal_lattice_vectors(latt)
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(latt)
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("Sites:")
for s in sites:
print(" {:4}: {:4}: ({:8.3g}, {:8.3g}, {:8.3g}) Z={:6.3g}".format(s[0], s[1], s[7][0], s[7][1], s[7][2], s[4]))
print("")
print("Conversion mode:", conversion_mode)
print("(tij) = |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij[0][0], tij[0][1], tij[0][2]))
print(" |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij[1][0], tij[1][1], tij[1][2]))
print(" |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij[2][0], tij[2][1], tij[2][2]))
acij = convert_lattice_vectors(aij, tij)
print(" Converted lattice vectors:")
print(" acx: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(acij[0][0], acij[0][1], acij[0][2]))
print(" acy: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(acij[1][0], acij[1][1], acij[1][2]))
print(" acz: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(acij[2][0], acij[2][1], acij[2][2]))
lattc = cal_lattice_parameters_from_lattice_vectors(acij)
print(" Lattice parameters: {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g} {:8.4g}".format(*lattc))
infc = cal_metrics(lattc)
gcij = infc['gij']
print(" Metric tensor:")
print(" gcij: ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gcij[0]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gcij[1]))
print(" ({:10.4g}, {:10.4g}, {:10.4g}) A".format(*gcij[2]))
volumec = cal_volume(acij)
print(" Volumec: {:12.4g} A^3".format(volumec))
print(" Sites:")
tij_x2xc = get_conversion_matrix_from_tij(tij, 'x2xc')
tij_x2xc = np.linalg.inv(tij).transpose()
print(" (tij:x to xc) = |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij_x2xc[0][0], tij_x2xc[0][1], tij_x2xc[0][2]))
print(" |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij_x2xc[1][0], tij_x2xc[1][1], tij_x2xc[1][2]))
print(" |{:8.3g}, {:8.3g}, {:8.3g}|".format(tij_x2xc[2][0], tij_x2xc[2][1], tij_x2xc[2][2]))
print(" Converted sites:")
csites = []
for site in sites:
name, label, z, M, q, r, color, pos = copy.deepcopy(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):
posc = convert_fractional_coordinates_by_tij([pos01[0] + ix, pos01[1] + iy, pos01[2] + iz], tij)
posn = [reduce01(posc[0]), reduce01(posc[1]), reduce01(posc[2])]
if add_site(csites, [name, label, z, M, q, r * 0.3, color, posn], gcij, rmin):
print(" {:4}: {:4}: ({:8.3g}, {:8.3g}, {:8.3g}) Z={:6.3g}".format(name, label, *posn, q))
allsites = []
for site in sites:
name, label, z, M, q, r, color, pos = copy.deepcopy(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 -0.1 <= posn[0] <= 1.1 and -0.1 <= posn[1] <= 1.1 and -0.1 <= posn[2] <= 1.1:
add_site(allsites, [name, label, z, M, q, r, color, posn], gij, rmin)
allcsites = []
for site in csites:
name, label, z, M, q, r, color, pos = copy.deepcopy(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 -0.1 <= posn[0] <= 1.1 and -0.1 <= posn[1] <= 1.1 and -0.1 <= posn[2] <= 1.1:
add_site(allcsites, [name, label, z, M, q, r, color, posn], gcij, rmin)
print("")
print("All sites to draw:")
for s in allsites:
print(" {:4}: {:4}: ({:8.3g}, {:8.3g}, {:8.3g}) Z={:6.3g}".format(s[0], s[1], s[7][0], s[7][1], s[7][2], s[4]))
print("")
print("All converted sites to draw:")
for s in allcsites:
print(" {:4}: {:4}: ({:8.3g}, {:8.3g}, {:8.3g}) Z={:6.3g}".format(s[0], s[1], s[7][0], s[7][1], s[7][2], s[4]))
fig = plt.figure(figsize = figsize)
ax = fig.add_subplot(111, projection='3d')
# Real space unit cell
draw_unitcell(ax, allsites, aij, nrange, kr, linecolor = 'black')
draw_unitcell(ax, allcsites, acij, nrange, kr, linecolor = 'blue', atomcolor = converted_atom_color)
# 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("")
terminate()
if __name__ == '__main__':
main()