Convert unit cell to different base vectors and draw those unit cells
Requirement: tkcrystalbase.py

Download script from .\crystal_convert_cell.py

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
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()