以下は、tkcrystalbase.py ライブラリのSphinx対応Markdownドキュメントです。

このドキュメントは、Sphinxの myst_parser と sphinx.ext.autodoc 拡張機能が有効になっていることを前提としています。tkcrystalbase.py ファイルの各関数、クラス、変数には、Sphinxが解釈できるreStructuredText形式のdocstringが追加されている必要があります。

`tkcrystalbase` モジュールリファレンス

このドキュメントは、tkcrystalbase.py モジュールの機能について説明します。

.. module:: tkcrystalbase :synopsis: 結晶構造計算と描画に関連するユーティリティ関数を提供するモジュール。 :platform: Any :audience: 結晶構造解析を行う研究者、開発者 :stability: Beta

概要

このモジュールは、結晶の格子定数から実空間および逆空間の格子ベクトル、メトリックテンソル、単位格子の体積などを計算する機能を提供します。また、原子の分数座標とデカルト座標間の変換、原子間距離や角度の計算、そしてmatplotlibを使用した3D描画ヘルパー関数を含んでいます。物理定数や一般的な格子パラメータ、サイト情報もここで定義されています。

定数

.. currentmodule:: tkcrystalbase

汎用定数 ^^^^^^

.. autodata:: pi :annotation: = 3.141592653589793

.. autodata:: pi2 :annotation: = 6.283185307179586

.. autodata:: torad :annotation: = 0.01745329251944 rad/deg 度をラジアンに変換する係数 (π/180)。

.. autodata:: todeg :annotation: = 57.29577951472 deg/rad ラジアンを度に変換する係数 (180/π)。

.. autodata:: basee :annotation: = 2.71828183

物理定数 ^^^^^^^^

以下の物理定数は、CODATAなどの標準に基づいています。

.. autodata:: h :annotation: = 6.6260755e-34 Js

.. autodata:: h_bar :annotation: = 1.05459e-34 Js

.. autodata:: hbar :annotation: = 1.05459e-34 Js h_bar のエイリアス。

.. autodata:: c :annotation: = 2.99792458e8 m/s

.. autodata:: e :annotation: = 1.60218e-19 C

.. autodata:: me :annotation: = 9.1093897e-31 kg

.. autodata:: mp :annotation: = 1.6726231e-27 kg

.. autodata:: mn :annotation: = 1.67495e-27 kg

.. autodata:: u0 :annotation: = 4.0 * 3.14 * 1e-7 N/A^2 真空の透磁率。

.. autodata:: e0 :annotation: = 8.854418782e-12 F/m 真空の誘電率。

.. autodata:: e2_4pie0 :annotation: = 2.30711e-28 Nm^2 クーロン力定数 (e^2 / 4πε₀)。

.. autodata:: a0 :annotation: = 5.29177e-11 m

.. autodata:: kB :annotation: = 1.380658e-23 J/K

.. autodata:: NA :annotation: = 6.0221367e23 mol^-1

.. autodata:: R :annotation: = 8.31451 J/K/mol

.. autodata:: F :annotation: = 96485.3 C/mol

.. autodata:: g :annotation: = 9.81 m/s^2

格子パラメータとサイト情報 ^^^^^^^^^^^^^^^^^^^^^^^^^^

.. autodata:: lattice_parameters :annotation: = [ 5.62, 5.62, 5.62, 90.0, 90.0, 90.0]

.. autodata:: sites :annotation: = [...] 各要素は [原子名 (str), サイトラベル (str), 原子番号 (int), 原子質量 (float), 電荷 (float), 半径 (float), 色 (str), 分数座標 (numpy.ndarray)] の形式です。デフォルトはNaCl構造の例。

関数

.. autofunction:: reduce01

.. autofunction:: round01

.. autofunction:: round_parameter

.. autofunction:: cal_lattice_vectors

.. autofunction:: cal_metrics

.. autofunction:: cal_volume

.. autofunction:: cal_reciprocal_lattice_vectors

.. autofunction:: cal_reciprocal_lattice_parameters

.. autofunction:: fractional_to_cartesian

.. autofunction:: distance2

.. autofunction:: distance

.. autofunction:: angle

.. autofunction:: configure_axis_structure

.. autofunction:: draw_box

.. autofunction:: draw_unitcell

.. autofunction:: main :no-index:

.. note::

    この関数は通常、モジュールがスクリプトとして直接実行された場合に呼び出されます。
    ライブラリとしてインポートして利用する際には直接呼び出すことは推奨されません。
    また、内部で :py:func:`exit` を呼び出すため、プログラム全体が終了します。

`tkcrystalbase.py` (Sphinx対応docstringが追記されたコード)

"""
tkcrystalbase.py

概要: 結晶構造計算と描画に関連するユーティリティ関数を提供するモジュール。
詳細説明:
    このモジュールは、結晶の格子定数から実空間および逆空間の格子ベクトル、
    メトリックテンソル、単位格子の体積などを計算する機能を提供します。
    また、原子の分数座標とデカルト座標間の変換、原子間距離や角度の計算、
    そしてmatplotlibを使用した3D描画ヘルパー関数を含んでいます。
    物理定数や一般的な格子パラメータ、サイト情報もここで定義されています。

:関連リンク:
    :doc:`tkcrystalbase_usage`
"""

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

# -----------------------------------------------------------------------------
# 汎用定数
# -----------------------------------------------------------------------------

pi          = 3.14159265358979323846
"""float: 円周率."""
pi2         = pi + pi
"""float: 2π."""
torad       = 0.01745329251944 # rad/deg
"""float: 度をラジアンに変換する係数 (π/180)。"""
todeg       = 57.29577951472   # deg/rad
"""float: ラジアンを度に変換する係数 (180/π)。"""
basee       = 2.71828183
"""float: 自然対数の底 e."""

# -----------------------------------------------------------------------------
# 物理定数
# -----------------------------------------------------------------------------

h           = 6.6260755e-34    # Js
"""float: プランク定数 (Js)."""
h_bar       = 1.05459e-34      # Js
"""float: 換算プランク定数 (h / 2π) (Js)."""
hbar        = h_bar
"""float: 換算プランク定数 (h_bar のエイリアス) (Js)."""
c           = 2.99792458e8     # m/s
"""float: 真空中の光速 (m/s)."""
e           = 1.60218e-19      # C
"""float: 電気素量 (C)."""
me          = 9.1093897e-31    # kg
"""float: 電子質量 (kg)."""
mp          = 1.6726231e-27    # kg
"""float: 陽子質量 (kg)."""
mn          = 1.67495e-27      # kg
"""float: 中性子質量 (kg)."""
u0          = 4.0 * 3.14*1e-7  # Ns^2C^-2
"""float: 真空の透磁率 (N/A^2 または H/m)."""
e0          = 8.854418782e-12  # C^2N^-1m^-2
"""float: 真空の誘電率 (F/m)."""
e2_4pie0    = 2.30711e-28      # Nm^2
"""float: クーロン力定数 (e^2 / 4πε₀) (Nm^2)."""
a0          = 5.29177e-11      # m
"""float: ボーア半径 (m)."""
kB          = 1.380658e-23     # JK^-1
"""float: ボルツマン定数 (J/K)."""
NA          = 6.0221367e23     # mol^-1
"""float: アボガドロ定数 (mol^-1)."""
R           = 8.31451          # J/K/mol
"""float: 理想気体定数 (J/K/mol)."""
F           = 96485.3          # C/mol
"""float: ファラデー定数 (C/mol)."""
g           = 9.81             # m/s2
"""float: 標準重力加速度 (m/s^2)."""

# -----------------------------------------------------------------------------
# デフォルトの格子パラメータとサイト情報
# -----------------------------------------------------------------------------

# 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]
"""list of float: 単位格子の格子定数 `[a, b, c, alpha, beta, gamma]`。
長さはオングストローム、角度は度で与えられます。デフォルトは面心立方格子の辺長と角度。
"""

# 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])]
        ]
"""list of list: 単位格子内の原子サイト情報を定義するリスト。
各要素は `[原子名 (str), サイトラベル (str), 原子番号 (int), 原子質量 (float),
電荷 (float), 半径 (float), 色 (str), 分数座標 (numpy.ndarray)]` の形式です。
デフォルトはNaCl構造の例。
"""


def reduce01(x: float) -> float:
    """数値を0以上1未満の範囲に変換する。

    入力された数値の小数部分を返します。
    例えば、3.75は0.75に、-2.3は0.7に変換されます。

    :param x: 変換する数値。
    :type x: float
    :returns: 0以上1未満の範囲に変換された値。
    :rtype: float
    """
    return x - int(x)


def round01(x: float) -> float:
    """数値を特定の閾値に基づいて0または1に丸めるか、0以上1未満の範囲に変換する。

    `x` が1.0に近い（絶対差が0.0002未満）場合は1.0を返します。
    `x` が0.0に近い（絶対差が0.0002未満）場合は0.0を返します。
    それ以外の場合、`x > 1.0` なら `x - int(x)` を返します。
    `x < 1.0` なら `x - int(x) + 1.0` を返します。
    この関数は、負の数や1を超える数値を0以上1未満の範囲に正規化するために使用されることがあります。

    :param x: 丸める、または変換する数値。
    :type x: float
    :returns: 丸められた、または0以上1未満の範囲に変換された値。
    :rtype: float
    """
    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: # This handles negative numbers and fractions less than 1.0. For positive fractions, it returns itself.
        return x - int(x) + 1.0
    return x # Should only be reached if x is exactly 1.0 or very close but not caught by abs(x-1.0) < 0.0002


def round_parameter(x: float, tol: float) -> float:
    """数値を指定された許容範囲 (tolerance) の整数倍に丸める。

    入力値 `x` を `tol` の整数倍に丸めます。
    丸め処理は、`(x + 0.1 * tol) / tol` を整数に変換し、`tol` を掛けることで行われます。
    これにより、例えば `tol=0.5` で `x=1.2` なら `1.0` に、`x=1.3` なら `1.5` に丸められます。

    :param x: 丸める数値。
    :type x: float
    :param tol: 丸めの基準となる許容範囲（単位）。
    :type tol: float
    :returns: `tol` の整数倍に丸められた値。
    :rtype: float
    """
    val = tol * int( (x+0.1*tol) / tol )
    return val 


def cal_lattice_vectors(latt: list[float]) -> np.ndarray:
    """単位格子の格子定数から実空間の格子ベクトルを計算する。

    格子定数 (a, b, c, α, β, γ) を受け取り、
    直交座標系における3つの格子ベクトル a_vec, b_vec, c_vec を含む3x3行列を計算します。
    格子ベクトルは、x軸に沿って a_vec が配置されるように構成されます。
    角度は度数で指定されますが、内部ではラジアンに変換して計算されます。

    :param latt: 格子定数のリスト `[a, b, c, alpha, beta, gamma]`。
                 長さはオングストローム、角度は度で与えられます。
    :type latt: list[float]
    :returns: 3x3 の格子ベクトル行列。
              `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :rtype: numpy.ndarray
    :raises ValueError: `sin(gamma)` がほぼゼロの場合、つまり格子が3次元でない場合に発生します。
    """
    cosa = cos(torad * latt[3])
    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
    
    # Robustly handle division by sin(gamma)
    if abs(sing) < 1.0e-8: # If gamma is 0 or 180 degrees, sin(gamma) is zero (degenerate lattice)
        if abs(latt[2] * (cosa - cosb * cosg)) < 1.0e-8: # If numerator is also near zero
            aij[2][1] = 0.0 # Assign 0 if gamma is ~0/180 and numerator is also ~0 (implies axes are aligned)
        else:
            raise ValueError("Lattice angle gamma is too close to 0 or 180 degrees, leading to a degenerate lattice.")
    else:
        aij[2][1] = latt[2] * (cosa - cosb * cosg) / sing

    # Calculate aij[2][2] and handle potential NaN from sqrt due to floating point errors or invalid lattice.
    val_for_sqrt = latt[2] * latt[2] - aij[2][0] * aij[2][0] - aij[2][1] * aij[2][1]
    if val_for_sqrt < -1.0e-8: # Allow slight negative for floating point errors
        raise ValueError("Invalid lattice parameters: c is too short or vectors are not physically possible.")
    aij[2][2] = sqrt(max(0.0, val_for_sqrt)) # Prevent sqrt of small negative number

    return aij

def cal_metrics(latt: list[float]) -> dict:
    """格子定数から格子ベクトルとメトリックテンソルを計算する。

    単位格子の格子定数から実空間の格子ベクトル `aij` を計算し、
    それらの内積からメトリックテンソル `gij` を構築します。
    結果は辞書形式で返されます。

    :param latt: 格子定数のリスト `[a, b, c, alpha, beta, gamma]`。
                 長さはオングストローム、角度は度で与えられます。
    :type latt: list[float]
    :returns: 計算された格子ベクトル行列 (`'aij'`) とメトリックテンソル (`'gij'`)
              をキーに持つ辞書。
              `{'aij': numpy.ndarray, 'gij': numpy.ndarray}` の形式。
    :rtype: dict
    """
    inf = {}

    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: np.ndarray) -> float:
    """格子ベクトルから単位格子の体積を計算する。

    3つの格子ベクトル (a_vec, b_vec, c_vec) のスカラー三重積 (a_vec・(b_vec × c_vec)) を用いて、
    単位格子の体積を計算します。

    :param aij: 3x3 の格子ベクトル行列。
                `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :type aij: numpy.ndarray
    :returns: 単位格子の体積。
    :rtype: float
    """
    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: np.ndarray) -> list[np.ndarray]:
    """実空間格子ベクトルから逆空間格子ベクトルを計算する。

    実空間の格子ベクトル `aij` (a_vec, b_vec, c_vec) を用いて、
    逆空間の格子ベクトル (a*_vec, b*_vec, c*_vec) を計算します。
    逆空間ベクトルは、単位格子の体積 V と実空間ベクトルの外積によって定義されます。

    :param aij: 3x3 の実空間格子ベクトル行列。
                `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :type aij: numpy.ndarray
    :returns: 逆空間の3つの格子ベクトル `[Ra_vec, Rb_vec, Rc_vec]` のリスト。
    :rtype: list[numpy.ndarray]
    :raises ValueError: 単位格子の体積がゼロの場合、つまり実空間格子ベクトルが線形従属である場合に発生します。
    """
    V = cal_volume(aij)
    if V == 0 or abs(V) < 1.0e-10: # Check for near-zero volume
        raise ValueError("Lattice volume is zero or too small, reciprocal lattice vectors cannot be calculated.")
    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: list[np.ndarray]) -> list[float]:
    """逆空間格子ベクトルから逆空間の格子定数を計算する。

    逆空間の格子ベクトル `Raij` (a*_vec, b*_vec, c*_vec) を用いて、
    逆空間の格子定数 (a*, b*, c*, α*, β*, γ*) を計算します。
    長さはベクトルのノルムから、角度はベクトルの内積から導出されます。

    :param Raij: 逆空間の3つの格子ベクトル `[Ra_vec, Rb_vec, Rc_vec]` のリスト。
    :type Raij: list[numpy.ndarray]
    :returns: 逆空間の格子定数 `[Ra, Rb, Rc, Ralpha, Rbeta, Rgamma]` のリスト。
    :rtype: list[float]
    """
    Ra = la.norm(Raij[0])
    Rb = la.norm(Raij[1])
    Rc = la.norm(Raij[2])

    Ralpha, Rbeta, Rgamma = 0.0, 0.0, 0.0 # Default values in case of degenerate vectors

    # Only calculate angles if norms are non-zero to avoid division by zero or NaN
    if Ra > 1.0e-10 and Rb > 1.0e-10 and Rc > 1.0e-10:
        # np.clip is used to prevent arccos from receiving values outside [-1, 1] due to floating point inaccuracies.
        Ralpha = todeg * arccos(np.clip(np.dot(Raij[1], Raij[2]) / Rb / Rc, -1.0, 1.0))
        Rbeta  = todeg * arccos(np.clip(np.dot(Raij[2], Raij[0]) / Rc / Ra, -1.0, 1.0))
        Rgamma = todeg * arccos(np.clip(np.dot(Raij[0], Raij[1]) / Ra / Rb, -1.0, 1.0))

    return [Ra, Rb, Rc, Ralpha, Rbeta, Rgamma]

def fractional_to_cartesian(pos: np.ndarray | list[float], aij: np.ndarray) -> tuple[float, float, float]:
    """分数座標からデカルト座標に変換する。

    単位格子内の原子の分数座標 (x_frac, y_frac, z_frac) と格子ベクトル `aij` を用いて、
    原子のデカルト座標 (x_cart, y_cart, z_cart) を計算します。

    :param pos: 変換する分数座標 `[x_frac, y_frac, z_frac]`。
    :type pos: numpy.ndarray or list[float]
    :param aij: 3x3 の格子ベクトル行列。
                `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :type aij: numpy.ndarray
    :returns: 計算されたデカルト座標 `(x_cart, y_cart, z_cart)`。
    :rtype: tuple[float, float, float]
    """
    pos_array = np.asarray(pos) # Ensure pos is treated as a numpy array for vectorized operations
    x = pos_array[0] * aij[0][0] + pos_array[1] * aij[1][0] + pos_array[2] * aij[2][0]
    y = pos_array[0] * aij[0][1] + pos_array[1] * aij[1][1] + pos_array[2] * aij[2][1]
    z = pos_array[0] * aij[0][2] + pos_array[1] * aij[1][2] + pos_array[2] * aij[2][2]

    return x, y, z

def distance2(pos0: np.ndarray | list[float], pos1: np.ndarray | list[float], gij: np.ndarray) -> float:
    """2つの分数座標間の距離の2乗を計算する。

    2つの原子の分数座標 `pos0` と `pos1` の間の距離の2乗を、
    メトリックテンソル `gij` を使用して計算します。
    この関数は周期境界条件を考慮しません。

    :param pos0: 1点目の分数座標 `[x0, y0, z0]`。
    :type pos0: numpy.ndarray or list[float]
    :param pos1: 2点目の分数座標 `[x1, y1, z1]`。
    :type pos1: numpy.ndarray or list[float]
    :param gij: 3x3 のメトリックテンソル。
    :type gij: numpy.ndarray
    :returns: 2点間の距離の2乗。
    :rtype: float
    """
    dx = np.asarray(pos1) - np.asarray(pos0)
    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: np.ndarray | list[float], pos1: np.ndarray | list[float], gij: np.ndarray) -> float:
    """2つの分数座標間の距離を計算する。

    2つの原子の分数座標 `pos0` と `pos1` の間の距離を、
    メトリックテンソル `gij` を使用して計算します。
    この関数は周期境界条件を考慮しません。

    :param pos0: 1点目の分数座標 `[x0, y0, z0]`。
    :type pos0: numpy.ndarray or list[float]
    :param pos1: 2点目の分数座標 `[x1, y1, z1]`。
    :type pos1: numpy.ndarray or list[float]
    :param gij: 3x3 のメトリックテンソル。
    :type gij: numpy.ndarray
    :returns: 2点間の距離。
    :rtype: float
    """
    r2 = distance2(pos0, pos1, gij)
    r = sqrt(r2)

    return r

def angle(pos0: np.ndarray | list[float], pos1: np.ndarray | list[float], pos2: np.ndarray | list[float], gij: np.ndarray) -> float:
    """3つの分数座標で定義される角度を計算する。

    `pos0` を頂点とし、`pos0-pos1` と `pos0-pos2` の2つのベクトルがなす角度を計算します。
    角度の計算にはメトリックテンソル `gij` が使用されます。
    結果は度数で返され、180度を超える場合は `360.0 - angle` として調整されます。

    :param pos0: 角度の頂点となる分数座標 `[x0, y0, z0]`。
    :type pos0: numpy.ndarray or list[float]
    :param pos1: 1つ目のベクトルを定義する分数座標 `[x1, y1, z1]`。
    :type pos1: numpy.ndarray or list[float]
    :param pos2: 2つ目のベクトルを定義する分数座標 `[x2, y2, z2]`。
    :type pos2: numpy.ndarray or list[float]
    :param gij: 3x3 のメトリックテンソル。
    :type gij: numpy.ndarray
    :returns: 3点間で形成される角度（度）。
              いずれかの距離が0の場合は0.0を返します。
    :rtype: float
    """
    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 = np.asarray(pos1) - np.asarray(pos0)
    dx02 = np.asarray(pos2) - np.asarray(pos0)
    
    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 = np.clip(ip / dis01 / dis02, -1.0, 1.0) # Clip to prevent arccos from receiving values outside [-1, 1]
    angle = todeg * arccos(cosa)
    if angle > 180.0: # Adjust angle to be between 0 and 180 (or 0 and 360 if explicitly wanted, but 180 is common for bond angles)
        angle = 360.0 - angle 

    return angle;

def configure_axis_structure(ax, xrange: tuple[float, float], yrange: tuple[float, float], zrange: tuple[float, float], fontsize: int = 12, legend_fontsize: int = 12) -> None:
    """Matplotlibの3D軸オブジェクトの表示設定を行う。

    MatplotlibのAxes3Dオブジェクト (`ax`) に対して、軸ラベルのフォントサイズ、
    軸の線、目盛り、グリッド、ペインの表示設定をカスタマイズします。
    デフォルトでは、軸の線、目盛り、ラベルは非表示に設定され、
    箱型のアスペクト比で視覚的に整えられます。

    :param ax: 設定を行うMatplotlibの3D軸オブジェクト。
    :type ax: matplotlib.axes.Axes
    :param xrange: x軸の表示範囲 `(xmin, xmax)`。
    :type xrange: tuple[float, float]
    :param yrange: y軸の表示範囲 `(ymin, ymax)`。
    :type yrange: tuple[float, float]
    :param zrange: z軸の表示範囲 `(zmin, zmax)`。
    :type zrange: tuple[float, float]
    :param fontsize: 軸ラベルおよび目盛りのフォントサイズ。
                     0以下の場合は軸ラベルと目盛りを非表示にします。デフォルトは12。
    :type fontsize: int, optional
    :param legend_fontsize: 凡例のフォントサイズ（この関数では直接使用されません）。デフォルトは12。
    :type legend_fontsize: int, optional
    :returns: None
    :rtype: None
    """
    ax.tick_params(axis = 'both', which = 'major', labelsize = fontsize)
    if fontsize > 0:
        ax.set_xlabel(f'$x$', fontsize = fontsize, labelpad = -5)
        ax.set_ylabel(f'$y$', fontsize = fontsize, labelpad = -5)
        ax.set_zlabel(f'$z$', fontsize = fontsize, labelpad = 0)

# x軸、y軸、z軸の線を非表示にする
    ax.xaxis.line.set_visible(False)
    ax.yaxis.line.set_visible(False)
    ax.zaxis.line.set_visible(False)
    
    ax.grid(False)
    ax.xaxis.pane.fill = False
    ax.yaxis.pane.fill = False
    ax.zaxis.pane.fill = False
    ax.xaxis.pane.set_edgecolor('w')
    ax.yaxis.pane.set_edgecolor('w')
    ax.zaxis.pane.set_edgecolor('w')
    ax.xaxis.pane.set_alpha(0)
    ax.yaxis.pane.set_alpha(0)
    ax.zaxis.pane.set_alpha(0)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_zticks([])
    ax.set_xticklabels([]) 
    ax.set_yticklabels([]) 
    ax.set_zticklabels([]) 

    ax.set_xlim(xrange)
    ax.set_ylim(yrange)
    ax.set_zlim(zrange)
    ax.set_aspect('equal','box')

def draw_box(ax, aij: np.ndarray, nrange: list[list[int]], color: str = 'black') -> None:
    """3D軸上に単位格子のボックスを描画する。

    格子ベクトル `aij` で定義される単位格子の12本の辺を、Matplotlibの3D軸 `ax` に描画します。
    この関数は、単一の単位格子を描画するために設計されており、`nrange` は現在の実装では使用されません。

    :param ax: 描画対象のMatplotlib 3D軸オブジェクト。
    :type ax: matplotlib.axes.Axes
    :param aij: 3x3 の格子ベクトル行列。
                `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :type aij: numpy.ndarray
    :param nrange: 描画するセルの範囲 `[[xmin, xmax], [ymin, ymax], [zmin, zmax]]`。
                   この関数では使用されません。
    :type nrange: list[list[int]]
    :param color: ボックスの線の色。デフォルトは 'black'。
    :type color: str, optional
    :returns: None
    :rtype: None
    """
# (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: list | None, aij: np.ndarray, nrange: list[list[int]], color: str = 'black', facecolor: str = 'black', edgecolor: str = 'white', alpha: float = 0.7, kr: float = 1.0) -> None:
    """3D軸上に単位格子とその内部のサイト原子を描画する。

    まず :py:func:`draw_box` 関数を呼び出して単位格子の境界線を描画します。
    次に、`sites` リストに含まれる各原子について、その分数座標をデカルト座標に変換し、
    指定された `nrange` の範囲内の周期的な位置に原子を散布図として描画します。
    原子の描画スタイル（色、透明度、サイズ）は引数でカスタマイズ可能です。

    :param ax: 描画対象のMatplotlib 3D軸オブジェクト。
    :type ax: matplotlib.axes.Axes
    :param sites: サイト情報（原子名、ラベル、原子番号、原子質量、電荷、半径、色、位置）のリスト。
                  各要素は `[name, label, z, M, q, r, color, pos]` の形式。
                  `None` の場合、原子は描画されません。
    :type sites: list of list or None
    :param aij: 3x3 の格子ベクトル行列。
                `[[ax_x, ax_y, ax_z], [ay_x, ay_y, ay_z], [az_x, az_y, az_z]]` の形式。
    :type aij: numpy.ndarray
    :param nrange: 描画するセルの範囲 `[[xmin, xmax], [ymin, ymax], [zmin, zmax]]`。
                   原子の繰り返し構造を描画するために使用されます。
    :type nrange: list[list[int]]
    :param color: 単位格子ボックスの線の色。デフォルトは 'black'。
    :type color: str, optional
    :param facecolor: サイト原子の面の描画色。デフォルトは 'black'。
    :type facecolor: str, optional
    :param edgecolor: サイト原子のエッジの描画色。デフォルトは 'white'。
    :type edgecolor: str, optional
    :param alpha: サイト原子の透明度（0.0から1.0）。デフォルトは 0.7。
    :type alpha: float, optional
    :param kr: サイト原子の描画サイズに対するスケール因子。デフォルトは 1.0。
    :type kr: float, optional
    :returns: None
    :rtype: None
    """
    draw_box(ax, aij, nrange, color)

    if sites is None: return

    for isite, site in enumerate(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]
            
            # Check if the periodic image is within the specified nrange's boundaries
            # The condition `nrange[i][1] < posn[i]` ensures that points *strictly* beyond the upper boundary are excluded.
            # Points *on* the upper boundary (e.g., 1.0 if nrange is [0,1]) are included in the default behavior
            # of `posn[i] < nrange[i][0]` or `nrange[i][1] < posn[i]`.
            # This is because `posn[i] < nrange[i][0]` is false for 1.0 < 0.
            # And `nrange[i][1] < posn[i]` is false for 1 < 1.0.
            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 = facecolor, edgecolors = edgecolor, alpha = alpha, s = kr *r)


def main() -> None:
    """モジュールの主要な計算と結果の表示を実行する。

    定義された :py:data:`lattice_parameters` を使用して、実空間および逆空間の
    格子ベクトル、メトリックテンソル、単位格子の体積を計算します。
    さらに、サンプルとして原子間距離と角度の計算も行い、
    その結果を標準出力に表示します。
    この関数は通常、スクリプトが直接実行された場合に呼び出されます。

    :returns: None
    :rtype: None
    """
    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() # This will terminate the Python interpreter. Generally not recommended in library functions.


if __name__ == '__main__':
    main()

tkcrystalbase モジュールリファレンス

概要

関連リンク

定数

関数

tkcrystalbase.py (Sphinx対応docstringが追記されたコード)

`tkcrystalbase` モジュールリファレンス

`tkcrystalbase.py` (Sphinx対応docstringが追記されたコード)