interatomic_distance.py ダウンロード/コピー

interatomic_distance.py をダウンロード

interatomic_distance.py

interatomic_distance.py

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 *


"""
概要: 原子間距離を計算するためのスクリプト。
詳細説明:
    結晶構造の格子パラメータとサイト情報に基づいて、指定された距離範囲内の原子間距離を計算し、表示します。
    TKCrystalBaseモジュールに依存しており、結晶学的な計算（格子ベクトル、計量テンソルなど）はそのモジュールが提供する関数を使用します。
    計算された原子間距離は、サイトの組み合わせと周期境界条件を考慮した単位胞の並進ベクトルと共にソートして出力されます。
関連リンク: :doc:`interatomic_distance_usage`
"""


# 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])]
        ]

# Distance range
rmin = 0.1  # angstrom. used to judge identical site
rmax = 4.5  # angstrom.


def main():
    """
    概要: プログラムの主要な処理を実行し、原子間距離を計算して出力する。
    詳細説明:
        この関数は以下のステップを実行します。
        1.  与えられた格子パラメータから格子ベクトル、計量テンソル、単位胞体積、およびそれらの逆格子版を計算し、標準出力に表示します。
        2.  原子間距離の計算範囲 `rmax` に基づいて、考慮すべき単位胞の最大並進数 (`nxmax`, `nymax`, `nzmax`) を決定します。
        3.  `sites` リストに定義された全ての原子サイトのペアについて、計算された単位胞の並進範囲内で原子間距離を計算します。
        4.  計算された距離が `rmin` と `rmax` の間にある場合、その距離と関連するサイト情報、並進ベクトルを `rlist` に追加します。
        5.  `rlist` を距離でソートした後、結果を整形して標準出力に表示します。
    """
    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))

# Calculate the range of unit cells
    nxmax = int(rmax / lattice_parameters[0]) + 1
    nymax = int(rmax / lattice_parameters[1]) + 1
    nzmax = int(rmax / lattice_parameters[2]) + 1
    print("")
    print("nmax:", nxmax, nymax, nzmax)

# Calculate interatomic distances and store them in rlist list variable
    rlist = []
    for isite0 in range(len(sites)):
        site0 = sites[isite0]
        name0, label0, z0, M0, q0, r0, color0, pos0 = site0
        for isite1 in range(len(sites)):
            site1 = sites[isite1]
            name1, label1, z1, M1, q1, r1, color1, pos1 = site1
            for iz in range(-nzmax, nzmax+1):
             for iy in range(-nymax, nymax+1):
              for ix in range(-nxmax, nxmax+1):
                dis = distance(pos0, pos1 + np.array([ix, iy, iz]), gij)
                if dis < rmin or rmax < dis:
                    continue

                rlist.append([dis, site0, site1, ix, iy, iz])

#                print("  {:4} ({:8.4g}, {:8.4g}, {:8.4g}) - {:4} ({:8.4g}, {:8.4g}, {:8.4g}) + ({:2d}, {:2d}, {:2d}): dis = {:10.4g} A"
#                    .format(label0, pos0[0], pos0[1], pos0[2], label1, pos1[0], pos1[1], pos1[2], ix, iy, iz, dis))

# Sort rlist by distance (x[0] priority)
    rlist.sort(key = lambda x: (x[0], x[1], x[3], x[4], x[5], x[2]))

    print("")
    print("Interatomic distances:")
    for rinf in rlist:
        dis, site0, site1, ix, iy, iz = rinf
        name0, label0, z0, M0, q0, r0, color0, pos0 = site0
        name1, label1, z1, M1, q1, r1, color1, pos1 = site1
        print("  {:4} ({:8.4g}, {:8.4g}, {:8.4g}) - {:4} ({:8.4g}, {:8.4g}, {:8.4g}) + ({:2d}, {:2d}, {:2d}): dis = {:10.4g} A"
            .format(label0, pos0[0], pos0[1], pos0[2], label1, pos1[0], pos1[1], pos1[2], ix, iy, iz, dis))

    print("")
    exit()


if __name__ == '__main__':
    main()