"""
結晶構造を管理するためのルートパッケージ。
このモモジュールは、`tkCrystalObject` クラスを提供し、結晶の格子パラメータ、空間群、
原子タイプ、原子サイトなどの情報を管理するための基本的な機能を提供します。
結晶の幾何学的特性や化学組成の分析、回折計算、座標変換など、様々な結晶学的な操作をサポートします。
:doc:`tkcrystalobject_usage`
"""
from numpy import exp, log, sin, cos, tan, arcsin, sqrt
import numpy as np
from numpy import linalg as la
from pprint import pprint
from tklib.tkobject import tkObject
import tklib.tkre as tkre
from tklib.tkre import DelQuote
from tklib.tkutils import SplitFilePath
from tklib.tksci.tksci import Round
from tklib.tksci.tksci import a0, NA, pi, pi2, torad, todeg, Round, Factorize, IsFactor, Factors
from tklib.tksci.tksci import degcos, degsin, degtan, acos, asin, atan, degacos, degasin, degatan
from tklib.tkcrystal.tkspacegroup import tkSpaceGroup
from tklib.tkcrystal.tkatomtypeobject import SplitAtomName, Charge2f
from tklib.tkcrystal.tkatomsiteobject import RoundParameter, Round01
from tklib.tkcrystal.tkatomtype import tkAtomType
from tklib.tkcrystal.tkatomsite import tkAtomSite
[ドキュメント]
class tkCrystalObject(tkObject):
"""
結晶構造データを管理するクラス。
結晶の格子パラメータ、空間群、原子タイプ、原子サイトなどの情報を一元的に管理します。
`tklib.tkobject.tkObject` を継承し、オブジェクトのパス管理などの基本機能も利用できます。
このクラスは、結晶の物理的・化学的特性を計算し、様々な結晶学的な分析を行うための基盤を提供します。
"""
def __init__(self, **args):
"""
tkCrystalObjectのコンストラクタ。
結晶構造オブジェクトを初期化し、内部変数を設定します。
"""
self.InitializeRoot()
def __del__(self):
"""
tkCrystalObjectのデストラクタ。
オブジェクトが破棄される際に特別な処理は行いません。
"""
pass
def __str__(self):
"""
オブジェクトの文字列表現を返します。
:returns: str - このオブジェクトのクラスパス。
"""
return self.ClassPath()
[ドキュメント]
def InitializeRoot(self):
"""
オブジェクトの内部変数を初期化します。
結晶名、サンプル名、化学組成、格子ベクトル、計量テンソル、空間群、
原子タイプリスト、原子サイトリスト、温度などをリセットします。
"""
self.__CrystalName = None
self.__SampleName = None
self.__chemicalcomposition = ''
self.__sumchemicalComposition = ''
self.__aij = None
self.__gij = None
self.__Latt = None
self.__Raij = None
self.__Rgij = None
self.__RLatt = None
self.__SpaceGroup = None
self.ClearSpaceGroup()
self.ClearAtomTypeList()
self.ClearAtomSiteList()
self.SetTemperature(25.0)
[ドキュメント]
def ClearAtomSiteList(self):
"""
原子サイトのリストをクリアします。
:returns: list - 空になった原子サイトのリスト。
"""
self.__AtomSiteList = []
return self.__AtomSiteList
[ドキュメント]
def ClearAtomTypeList(self):
"""
原子タイプのリストをクリアします。
:returns: list - 空になった原子タイプのリスト。
"""
self.__AtomTypeList = []
return self.__AtomTypeList
[ドキュメント]
def ClearSpaceGroup(self):
"""
空間群情報をクリアし、P1空間群に初期化します。
:returns: tklib.tkcrystal.tkspacegroup.tkSpaceGroup - 新しく初期化された空間群オブジェクト。
"""
self.__SpaceGroup = tkSpaceGroup()
self.ConventionalSPGName = ''
self.__SpaceGroup.SetP1()
return self.__SpaceGroup
[ドキュメント]
def RoundLatticeParameters(self, val, eps = 1.0e-8):
"""
格子定数の角度を特定の標準値（60, 90, 120度）に丸めます。
指定された値が90度、60度、120度に非常に近い場合、その標準値に丸めます。
それ以外の場合は、`RoundParameter` 関数を使用して丸めます。
:param val: float - 丸めたい格子定数の角度（度）。
:param eps: float, optional - 許容誤差。デフォルトは1.0e-8。
:returns: float - 丸められた格子定数の角度。
"""
if abs(val - 90.0) < eps:
return 90.0
elif abs(val - 60.0) < eps:
return 60.0
elif abs(val - 120.0) < eps:
return 120.0
else:
RoundParameter(val, eps)
return val
[ドキュメント]
def SetSampleName(self, name = None, MakeFromPath = 0):
"""
サンプル名を設定します。
`name` が指定された場合、それをサンプル名として設定します。
`MakeFromPath` が1の場合、オブジェクトのパスからファイル名を抽出し、それをサンプル名として設定します。
:param name: str, optional - 設定するサンプル名。Noneの場合、`MakeFromPath` の値に応じて処理されます。
:param MakeFromPath: int, optional - オブジェクトのパスからサンプル名を生成するかどうかを示すフラグ。
1の場合、パスから生成します。デフォルトは0。
:returns: None
"""
if name is not None:
self.__SampleName = name
return
if MakeFromPath:
dir, filename, filebody, ext = SplitFilePath(self.path)
self.__SampleName = filebody
[ドキュメント]
def SampleName(self, Create = 0):
"""
現在のサンプル名を取得します。
サンプル名が設定されていない場合、`Create` が1であればオブジェクトのパスからファイル名を抽出し、
それをサンプル名として返します。
:param Create: int, optional - サンプル名が未設定の場合にパスから生成するかどうかを示すフラグ。
1の場合、生成を試みます。デフォルトは0。
:returns: str or None - サンプル名、または生成に失敗した場合はNone。
"""
if self.__SampleName is not None:
return self.__SampleName
if Create:
try:
dir, filename, filebody, ext = SplitFilePath(self.path)
return filebody
except:
return None
[ドキュメント]
def SetCrystalName(self, name = None, MakeFromPath = 0):
"""
結晶名を設定します。
`name` が指定された場合、それを結晶名として設定します。
`MakeFromPath` が1の場合、オブジェクトのパスからファイル名を抽出し、それを結晶名として設定します。
:param name: str, optional - 設定する結晶名。Noneの場合、`MakeFromPath` の値に応じて処理されます。
:param MakeFromPath: int, optional - オブジェクトのパスから結晶名を生成するかどうかを示すフラグ。
1の場合、パスから生成します。デフォルトは0。
:returns: None
"""
if name is not None:
self.__CrystalName = name
return
if MakeFromPath:
dir, filename, filebody, ext = SplitFilePath(self.path)
self.__CrystalName = filebody
[ドキュメント]
def CrystalName(self, Create = 0):
"""
現在の結晶名を取得します。
結晶名が設定されていない場合、`Create` が1であれば化学組成またはオブジェクトのパスから
ファイル名を抽出し、それを結晶名として返します。
:param Create: int, optional - 結晶名が未設定の場合に生成するかどうかを示すフラグ。
1の場合、生成を試みます。デフォルトは0。
:returns: str or None - 結晶名、または生成に失敗した場合はNone。
"""
if self.__CrystalName is not None:
return self.__CrystalName
if Create:
comp = self.ChemicalComposition()
if comp != '':
return comp
dir, filename, filebody, ext = SplitFilePath(self.path)
return filebody
[ドキュメント]
def SetTemperature(self, T):
"""
温度を設定します。
:param T: float - 設定する温度。
:returns: None
"""
self.__temperature = T
[ドキュメント]
def Temperature(self,):
"""
現在の温度を取得します。
:returns: float - 現在設定されている温度。
"""
return self.__temperature
[ドキュメント]
def nAtomType(self):
"""
登録されている原子タイプの数を取得します。
:returns: int - 原子タイプの数。
"""
return len(self.__AtomTypeList)
[ドキュメント]
def AtomTypeList(self):
"""
登録されている原子タイプのリストを取得します。
:returns: list of tklib.tkcrystal.tkatomtype.tkAtomType - 原子タイプオブジェクトのリスト。
"""
return self.__AtomTypeList
[ドキュメント]
def natom_site(self):
"""
非対称単位内の原子サイトの数を取得します。
:returns: int - 原子サイトの数。
"""
return len(self.__AtomSiteList)
[ドキュメント]
def nAtomSite(self):
"""
非対称単位内の原子サイトの数を取得します。(エイリアス)
`natom_site` メソッドと同じです。
:returns: int - 原子サイトの数。
"""
return self.natom_site()
[ドキュメント]
def atom_site_list(self, mode = 0):
"""
原子サイトのリストを取得します。
`mode` が0の場合、非対称単位内の原子サイトのリストを返します。
それ以外の場合、展開された原子サイトのリストを返します。
:param mode: int, optional - 取得するリストの種類を指定するモード。
0: 非対称単位内の原子サイトリスト。
それ以外: 展開された原子サイトリスト。デフォルトは0。
:returns: list of tklib.tkcrystal.tkatomsite.tkAtomSite - 原子サイトオブジェクトのリスト。
"""
if mode == 0:
return self.__AtomSiteList
else:
return self.__ExpandedAtomSiteList
[ドキュメント]
def AtomSiteList(self, mode = 0):
"""
原子サイトのリストを取得します。(エイリアス)
`atom_site_list` メソッドと同じです。
:param mode: int, optional - 取得するリストの種類を指定するモード。
0: 非対称単位内の原子サイトリスト。
それ以外: 展開された原子サイトリスト。デフォルトは0。
:returns: list of tklib.tkcrystal.tkatomsite.tkAtomSite - 原子サイトオブジェクトのリスト。
"""
return self.atom_site_list(mode = mode)
[ドキュメント]
def nexpanded_atom_site(self):
"""
展開された原子サイトの総数を取得します。
:returns: int - 展開された原子サイトの数。
"""
return len(self.__ExpandedAtomSiteList)
[ドキュメント]
def nExpandedAtomSiteList(self):
"""
展開された原子サイトの総数を取得します。(エイリアス)
`nexpanded_atom_site` メソッドと同じです。
:returns: int - 展開された原子サイトの数。
"""
return self.nexpanded_atom_site()
[ドキュメント]
def expanded_atom_site_list(self):
"""
展開された（空間群操作によって生成された）原子サイトのリストを取得します。
:returns: list of tklib.tkcrystal.tkatomsite.tkAtomSite - 展開された原子サイトオブジェクトのリスト。
"""
return self.__ExpandedAtomSiteList
[ドキュメント]
def ExpandedAtomSiteList(self):
"""
展開された（空間群操作によって生成された）原子サイトのリストを取得します。(エイリアス)
`expanded_atom_site_list` メソッドと同じです。
:returns: list of tklib.tkcrystal.tkatomsite.tkAtomSite - 展開された原子サイトオブジェクトのリスト。
"""
return self.expanded_atom_site_list()
[ドキュメント]
def ChemicalComposition(self):
"""
計算された化学組成を取得します。
例: 'H2O'
:returns: str - 化学組成文字列。
"""
return self.__chemicalcomposition
[ドキュメント]
def SetChemicalComposition(self, composition):
"""
化学組成を設定します。
:param composition: str - 設定する化学組成文字列。
:returns: None
"""
self.__chemicalcomposition = composition
[ドキュメント]
def SumChemicalComposition(self):
"""
単位胞内の全原子に基づく合計化学組成を取得します。
例: 単位胞内にHが2個、Oが1個の場合 'H2O'。
この値は、`AnalyzeChemicalComposition` メソッドによって設定されます。
:returns: str - 合計化学組成文字列。
"""
try:
return self.__sumchemicalComposition
except:
return ''
#mode: all|asym|type
[ドキュメント]
def get_atom_name_list(self, mode = 'all', NameOnly = True):
"""
指定されたモードに基づいて原子名のリストを取得します。
:param mode: str, optional - 取得する原子サイトのモード。
'aym': 非対称単位内の原子サイト。 (typo: 'asym' の意図)
'all': 展開された全原子サイト。
'type': 原子タイプリスト。
デフォルトは 'all'。
:param NameOnly: bool, optional - 原子名のみを返すかどうか。Trueの場合、原子名（例: 'C'）を返します。
Falseの場合、完全な原子名（例: 'C1'）を返すことがあります。デフォルトはTrue。
:returns: list of str - 原子名のリスト。
"""
if mode == 'aym':
sites = self.AtomSiteList()
elif mode == 'all':
sites = self.ExpandedAtomSiteList()
elif mode == 'type':
sites = self.AtomTypeList()
else:
print("\nError in tkcrystal.get_atom_name_list: Invalide mode [{}]\n".format(mode))
exit()
namelist = []
for i in range(len(sites)):
site = sites[i]
if NameOnly:
namelist.append(sites[i].AtomName())
else:
namelist.append(sites[i].AtomNameOnly())
return namelist
#mode: all|asym
[ドキュメント]
def get_position_list(self, mode = 'all', IsReduce01 = False):
"""
指定されたモードに基づいて原子位置のリストを取得します。
:param mode: str, optional - 取得する原子サイトのモード。
'aym': 非対称単位内の原子サイト。 (typo: 'asym' の意図)
'all': 展開された全原子サイト。
デフォルトは 'all'。
:param IsReduce01: bool, optional - 位置座標を0から1の範囲に正規化するかどうか。
Trueの場合、正規化します。デフォルトはFalse。
:returns: list of list of float - 各原子サイトの位置座標 (例: [[x1, y1, z1], [x2, y2, z2], ...])。
"""
if mode == 'aym':
sites = self.AtomSiteList()
elif mode == 'all':
sites = self.ExpandedAtomSiteList()
else:
print("\nError in tkcrystal.get_atom_name_list: Invalide mode [{}]\n".format(mode))
exit()
poslist = []
for i in range(len(sites)):
site = sites[i]
poslist.append(sites[i].Position(IsReduce01))
return poslist
[ドキュメント]
def SetSumChemicalComposition(self, composition):
"""
単位胞内の全原子に基づく合計化学組成を設定します。
この値は、通常 `AnalyzeChemicalComposition` メソッドによって設定されます。
:param composition: str - 設定する合計化学組成文字列。
:returns: None
"""
self.__sumchemicalcomposition = composition
[ドキュメント]
def IsAnionic(self, name):
"""
指定された原子名が陰イオン性かどうかを判定します。
酸素(O), 硫黄(S), セレン(Se), テルル(Te), フッ素(F), 塩素(Cl), ヨウ素(I), 窒素(N),
ヒ素(As), アンチモン(Sb) を陰イオン性と見なします。
:param name: str - 判定する原子名（例: 'O', 'Cl'）。
:returns: int - 陰イオン性であれば1、そうでなければ0。
"""
if tkre.Match(r'(O|S|Se|Te|F|Cl|I|N|As|Sb)$', name, 'i'):
return 1
return 0
[ドキュメント]
def IsCationic(self, name):
"""
指定された原子名が陽イオン性かどうかを判定します。
`IsAnionic` メソッドの逆の結果を返します。
:param name: str - 判定する原子名。
:returns: bool - 陽イオン性であればTrue、そうでなければFalse。
"""
return not self.IsAnionic(name)
[ドキュメント]
def AnalyzeChemicalComposition(self, sort = ''):
"""
結晶の化学組成を分析し、設定します。
展開された原子サイトのリストから各原子の数をカウントし、化学組成と
単位胞あたりの化学式単位 (Z値) を計算して設定します。
ソートオプションにより、出力される組成の順序を制御できます。
:param sort: str, optional - 原子をソートする基準。
'atom_name': 原子名でソート。
'atomic_number': 原子番号でソート。
'atomic_mass': 原子量でソート。
'' (空文字列): 陽イオン性/陰イオン性に基づいてソート（デフォルト）。
:returns: str - 計算された化学組成文字列。
"""
atoms_dict = {}
AtomTypes = self.AtomTypeList()
for t in AtomTypes:
typea = t.AtomType()
typeo = t.AtomTypeOnly()
charge = t.Charge()
AN = t.AtomicNumber()
M = t.AtomicMass()
if AN is None: AN = 9999
if M == 0.0: M = 9999
atoms_dict[typeo] = {"atom_name": typeo, "charge": charge, "atomic_number": AN, "atomic_mass": M}
if sort == 'atom_name':
atoms_dict = dict(sorted(atoms_dict.items()))
elif sort == 'atomic_number':
atoms_dict = dict(sorted(atoms_dict.items(), key = lambda item: item[1]["atomic_number"]))
elif sort == 'atomic_mass':
atoms_dict = dict(sorted(atoms_dict.items(), key = lambda item: item[1]["atomic_mass"]))
ExpandedAtomSites = self.ExpandedAtomSiteList()
nExpandedAtomSite = self.nExpandedAtomSiteList()
AtomCount = {name: 0.0 for name, val in atoms_dict.items()}
for i in range(nExpandedAtomSite):
site = ExpandedAtomSites[i]
name = site.AtomNameOnly(1);
if '?' in name or '.' in name or '[' in name: continue
# if name not in AtomCount.keys():
# AtomCount[name] = 0.0
AtomCount[name] += site.Occupancy()
Names = list(AtomCount.keys())
for key in Names:
AtomCount[key] = RoundParameter(AtomCount[key], 1.0e-4)
# print("AtomCount=", AtomCount)
nk = len(Names)
MinN = 1000000
C1 = ''
C2 = ''
for name, d in atoms_dict.items():
n = AtomCount[name]
MinN = min([n, MinN])
if n == 1:
n = ''
else:
n = str(n)
n = tkre.Sub('0+$', '', n)
n = tkre.Sub(r'\.$', '', n)
if sort == '':
if self.IsCationic(name):
C1 += name + n
else:
C2 += name + n
else:
C1 += name + n
'''
for i in range(nk):
name = Names[i]
n = AtomCount[name]
MinN = min([n, MinN])
if n == 1:
n = ''
else:
n = str(n)
n = tkre.Sub('0+$', '', n)
n = tkre.Sub(r'\\.$', '', n)
if self.IsCationic(name):
C1 += name + n
else:
C2 += name + n
'''
Composition = C1 + C2
self.SetSumChemicalComposition(Composition)
# print("Composition=", Composition)
elements = Factors(MinN);
# elements = Factorize(MinN);
fac = 1
ne = len(elements)
for i in range(ne):
num = elements[i]
IsCommonFactor = 1
for j in range(nk):
name = Names[j];
n = AtomCount[name]
if IsFactor(n, num * fac) == 0:
IsCommonFactor = 0
if IsCommonFactor:
fac *= num
self.SetFormulaUnit(fac)
C1 = ''
C2 = ''
for name, d in atoms_dict.items():
n = AtomCount[name] / fac
if n == 1:
n = ''
else:
n = str(n)
n = tkre.Sub('0+$', '', n)
n = tkre.Sub(r'\.$', '', n)
if sort == '':
if self.IsAnionic(name) == 0:
C1 += name + n
else:
C2 += name + n
else:
C1 += name + n
'''
for i in range(nk):
name = Names[i]
n = AtomCount[name] / fac
if n == 1:
n = ''
else:
n = str(n)
n = tkre.Sub('0+$', '', n)
n = tkre.Sub(r'\\.$', '', n)
if self.IsAnionic(name) == 0:
C1 += name + n
else:
C2 += name + n
'''
Composition = C1 + C2
self.SetChemicalComposition(Composition)
return Composition
[ドキュメント]
def UnitcellWeight(self):
"""
単位胞の分子量を取得します。
この値は、通常 `CalculateDensity` メソッドによって設定されます。
:returns: float - 単位胞の分子量。未設定の場合は0.0。
"""
try:
return self.__UnitcellWeight
except:
return 0.0
[ドキュメント]
def SetUnitcellWeight(self, MW):
"""
単位胞の分子量を設定します。
:param MW: float - 設定する分子量。
:returns: None
"""
self.__UnitcellWeight = MW
[ドキュメント]
def Density(self):
"""
結晶の密度を取得します。
この値は、通常 `CalculateDensity` メソッドによって設定されます。
:returns: float - 結晶の密度。未設定の場合は0.0。
"""
try:
return self.__density
except:
return 0.0
[ドキュメント]
def AtomDensity(self):
"""
結晶の原子密度を計算して取得します。
原子密度は、展開された原子サイトの総数を単位胞体積（pm^3）で割った値です。
:returns: float - 結晶の原子密度 (A^-3)。
"""
AtomSites = self.ExpandedAtomSiteList()
v = self.Volume()
self.__AtomDensity = len(AtomSites) / (v * 1.0e-24)
return self.__AtomDensity
[ドキュメント]
def SetDensity(self, density):
"""
結晶の密度を設定します。
:param density: float - 設定する密度。
:returns: None
"""
self.__density = density
[ドキュメント]
def SetAtomDensity(self, density):
"""
結晶の原子密度を設定します。
:param density: float - 設定する原子密度。
:returns: None
"""
self.__AtomDensity = density
[ドキュメント]
def CalculateDensity(self):
"""
単位胞の分子量と結晶密度を計算し、設定します。
原子タイプリストと非対称単位内の原子サイトリストに基づいて、
単位胞あたりの分子量 (MW) と結晶密度を計算します。
計算された密度は4桁に丸められます。
:returns: float - 計算された結晶密度。
"""
AtomTypes = self.AtomTypeList()
AtomSites = self.AtomSiteList()
nAtomSite = self.nAtomSite()
MW = 0.0;
self.dprint(1, "tkcrystal.tkcrystal.CalculateDensity: nS=", nAtomSite)
for i in range(nAtomSite):
AtomSite = AtomSites[i]
AtomName = AtomSite.AtomNameOnly(0)
mult = AtomSite.Multiplicity()
iAtomType = self.FindiAtomType(AtomName)
if iAtomType < 0:
print("Warning in tkcrystalobject.CalculateDensity: Invalid iAtomType={iAtomType}. Skip")
# exit()
continue
AtomType = AtomTypes[iAtomType]
occ = AtomSite.Occupancy()
AtomicMass = AtomType.AtomicMass()
if occ is None:
occ = 1.0
if mult is None:
mult = 1
MW += AtomicMass * mult * occ
self.dprint(1, "tkcrystal.tkcrystal.CalculateDensity: "
+ "iAtomSite={}: M={:6.4f} m={:6.4f} occ={:6.4f} MW={:6.4f}"
.format(i, AtomicMass, mult, occ, MW))
self.SetUnitcellWeight(MW)
vol = self.Volume()
self.dprint(1, "tkcrystal.tkcrystal.CalculateDensity:MW=", MW)
self.dprint(1, "tkcrystal.tkcrystal.CalculateDensity:vol=", vol)
if vol == 0.0:
self.__density = 0.0
else:
d = MW / NA / vol * 1.0e24
self.__density = Round(d, 4)
return self.__density
#iAtomTypeは0から始まる
[ドキュメント]
def FindAtomType(self, atomtype):
"""
指定された原子タイプ名に対応する原子タイプリストのインデックスを検索します。
:param atomtype: str - 検索する原子タイプ名。
:returns: int - 見つかった場合はそのインデックス (0から始まる)。見つからない場合は-1。
"""
for i in range(len(self.__AtomTypeList)):
if atomtype == self.__AtomTypeList[i].AtomType():
return i
else:
return -1
#iAtomSiteは0から始まる
[ドキュメント]
def FindAtomSite(self, atomname, mode = 0):
"""
指定された原子名に対応する原子サイトオブジェクトを検索します。
`mode` が0の場合、非対称単位内の原子サイトリストから検索します。
:param atomname: str - 検索する原子名。
:param mode: int, optional - 検索するリストのモード。
0: 非対称単位内の原子サイトリスト。
デフォルトは0。
:returns: tklib.tkcrystal.tkatomsite.tkAtomSite or None - 見つかった場合はその原子サイトオブジェクト。見つからない場合はNone。
"""
if mode == 0:
AtomSites = self.AtomSiteList(mode)
else:
pass
for i in range(len(AtomSites)):
atom = AtomSites[i]
name = atom.AtomNameOnly()
if atomname.upper() == name.upper():
return atom
return None
#iAtomTypeは0から始まる
[ドキュメント]
def FindiAtomType(self, atomname):
"""
指定された原子名に対応する原子タイプリストのインデックスを検索します。(エイリアス)
`FindAtomType` メソッドと似ていますが、`SplitAtomName` を内部で使用して原子名を解析します。
:param atomname: str - 検索する原子名。
:returns: int - 見つかった場合はそのインデックス (0から始まる)。見つからない場合は-1。
"""
at = tkAtomType()
name, nametype, type, charge = SplitAtomName(atomname)
# name, nametype, type, charge = at.SplitAtomName(atomname)
name = name.upper()
AtomTypeList = self.AtomTypeList()
nAtomType = self.nAtomType()
for i in range(nAtomType):
name1 = AtomTypeList[i].AtomTypeOnly().upper()
if name == name1:
return i
return -1
[ドキュメント]
def find_nearest_site(self, pos0, irange = [1, 1, 1]):
"""
指定された位置 `pos0` に最も近い原子サイトを、展開された全原子サイトから検索し、そのインデックスを返します。
:param pos0: numpy.ndarray or list of float - 検索の基準となる分数座標 [x, y, z]。
:param irange: list of int, optional - 探索する隣接単位胞の範囲 [ix, iy, iz]。
各要素は中心単位胞からの最大オフセットを示します。デフォルトは[1, 1, 1]。
:returns: int - 最も近い原子サイトのリスト内インデックス。見つからなかった場合（リストが空の場合など）は-1。
"""
AtomSites = self.AtomSiteList(1)
minidx = -1
mindis = 1.0e10;
for i in range(len(AtomSites)):
atom = AtomSites[i]
pos1 = atom.Position(1)
dis = self.get_nearest_interatomic_distance(pos0, pos1, irange = irange)
if dis < mindis:
mindis = dis
minidx = i
return minidx
[ドキュメント]
def FindNearestSite(self, pos0, irange = [1, 1, 1]):
"""
指定された位置 `pos0` に最も近い原子サイトを、展開された全原子サイトから検索し、そのインデックスを返します。(エイリアス)
`find_nearest_site` メソッドと同じです。
:param pos0: numpy.ndarray or list of float - 検索の基準となる分数座標 [x, y, z]。
:param irange: list of int, optional - 探索する隣接単位胞の範囲 [ix, iy, iz]。
各要素は中心単位胞からの最大オフセットを示します。デフォルトは[1, 1, 1]。
:returns: int - 最も近い原子サイトのリスト内インデックス。見つからなかった場合（リストが空の場合など）は-1。
"""
return self.find_nearest_site(pos0, irange = [1, 1, 1])
[ドキュメント]
def BurstToP1(self):
"""
現在の結晶構造をP1空間群に変換し、すべての原子を非対称単位に展開します。
この操作により、すべての空間群対称性が除去され、すべての原子が独立したサイトとして扱われます。
現在の空間群がP1に設定され、展開された原子サイトリストが非対称単位の原子サイトリストとして再設定され、
その後座標が再度展開されます。
"""
# my @ExpandedAtomSiteList = $this->GetCExpandedAtomSiteList();
# SPG = $this->GetCSpaceGroup();
self.SpaceGroup().SetP1()
# self.SetSpaceGroup("P 1", 1, 1)
# self.SetLatticeSystem('triclinic', 0)
# 警告: self.SetAtomSiteList メソッドは定義されていません。__AtomSiteListを直接操作する意図かもしれません。
self.__AtomSiteList = self.ExpandedAtomSiteList()
self.ExpandCoordinates();
[ドキュメント]
def GetBravaisLattice(self):
"""
現在の空間群のブラベー格子名を取得します。
:returns: str - ブラベー格子名。
"""
SPG = self.SpaceGroup()
return SPG.BravaisLattice()
[ドキュメント]
def expand_coordinates(self, DoTranslation = 1):
"""
非対称単位内の原子サイトを空間群操作と並進操作に基づいて展開し、全原子サイトのリストを生成します。
このメソッドは、`__ExpandedAtomSiteList` を生成し、各非対称原子サイトの多重度を設定し、
さらに単位胞の分子量、密度、化学組成を再計算します。
:param DoTranslation: int, optional - 並進対称操作を適用するかどうかを示すフラグ。
1の場合、並進操作も考慮します。デフォルトは1。
:returns: None
"""
#初期化
self.__nTotalExpandedAtomSite = 0;
self.__Density = 0.0;
#同一原子の判断半径を決める
latt = self.LatticeParameters()
OverlapCriteriaRadius = 0.001 * min([latt[0], latt[1], latt[2]])
#非対称構造中のパラメータ
AtomTypes = self.AtomTypeList()
AtomSites = self.AtomSiteList()
nAtomType = len(AtomTypes)
nAtomSite = len(AtomSites)
SPG = self.SpaceGroup()
if DoTranslation == 1:
nTranslation = 1
else:
nTranslation = SPG.nTranslation()
nSymmetryOperation = SPG.nSymmetryOperation()
nTotalSymmetry = nTranslation * nSymmetryOperation
#座標展開を開始
ExpandedAtomSiteList = []
nExpandedAtomSiteList = []
count = 0;
for ia in range(nAtomSite):
nExpandedAtomSiteList.append(0)
atom = AtomSites[ia]
label = atom.Label()
AtomName = atom.AtomName()
Charge = atom.Charge()
iAtomType = self.FindiAtomType(AtomName)
AtomType = AtomTypes[iAtomType]
pos = atom.Position(1)
pos[0] = self.RoundSymmetricPosition(pos[0])
pos[1] = self.RoundSymmetricPosition(pos[1])
pos[2] = self.RoundSymmetricPosition(pos[2])
Occupancy = atom.Occupancy()
force = atom.Force()
v = atom.Velocity()
atom.SetId(ia)
# print("ia={}: {} ({}, {}, {})".format(ia, label, pos[0], pos[1], pos[2]))
ExpandedAtomSiteList.append(tkAtomSite())
ExpandedAtomSiteList[count].SetIdAsymmetricAtomSite(ia)
ExpandedAtomSiteList[count].SetiAtomType(iAtomType)
ExpandedAtomSiteList[count].SetAtomType(AtomType)
ExpandedAtomSiteList[count].SetLabel(label)
ExpandedAtomSiteList[count].SetAtomName(AtomName)
ExpandedAtomSiteList[count].SetCharge(Charge)
ExpandedAtomSiteList[count].SetPosition(pos)
if force is not None:
ExpandedAtomSiteList[count].SetForce(force)
ExpandedAtomSiteList[count].SetOccupancy(Occupancy)
if v is not None:
ExpandedAtomSiteList[count].SetVelocity(v)
nExpandedAtomSiteList[ia] += 1
count += 1
#print "($ia) x,y,z: $x,$y,$z<br>\n";
for iss in range(nTotalSymmetry):
symop = SPG.SymmetryOperation(iss)
poss = SPG.DoSymmetryOperation(iss, pos, 1)
forces = SPG.DoVelocitySymmetryOperation(iss, force, 1)
vs = SPG.DoVelocitySymmetryOperation(iss, v, 1)
# print(" iss={}: [{}] to ({}, {}, {})".format(iss, symop, poss[0], poss[1], poss[2]))
IsPassed = 1
for it in range(count):
id1 = ExpandedAtomSiteList[it].Id()
# print(" id1={}".format(id1))
if ia != id1:
continue
# iat = ExpandedAtomSiteList[it].iAtomType()
# if iat != iAtomType:
# continue
pos0 = ExpandedAtomSiteList[it].Position(1)
# print(" id1={}: ({}, {}, {})".format(id1, pos0[0], pos0[1], pos0[2]))
dis = self.GetNearestInterAtomicDistance(pos0, poss)
# print(" dis=", dis, " < {}?".format(OverlapCriteriaRadius))
if dis < OverlapCriteriaRadius:
IsPassed = 0
break
if IsPassed == 0:
continue
ExpandedAtomSiteList.append(tkAtomSite())
ExpandedAtomSiteList[count].SetIdAsymmetricAtomSite(ia)
ExpandedAtomSiteList[count].SetiAtomType(iAtomType)
ExpandedAtomSiteList[count].SetAtomType(AtomType)
ExpandedAtomSiteList[count].SetLabel(label)
ExpandedAtomSiteList[count].SetAtomName(AtomName)
ExpandedAtomSiteList[count].SetCharge(Charge)
ExpandedAtomSiteList[count].SetPosition(poss)
ExpandedAtomSiteList[count].SetForce(forces)
ExpandedAtomSiteList[count].SetOccupancy(Occupancy)
ExpandedAtomSiteList[count].SetVelocity(vs)
#Idは1からはじまる
ExpandedAtomSiteList[count].SetId(ia)
#print "v: ($vx,$vy,$vz) ($vxs, $vys, $vzs)<br>\n";
nExpandedAtomSiteList[ia] += 1
count += 1
self.__ExpandedAtomSiteList = ExpandedAtomSiteList
self.__nTotalExpandedAtomSite = count
self.__nExpandedAtomSite = nExpandedAtomSiteList
#多重度を、非対称原子に設定
for ia in range(nAtomSite):
atom = AtomSites[ia]
atom.SetMultiplicity(nExpandedAtomSiteList[ia])
#多重度を、全原子に設定
self.__MultiplitictyList = []
for i in range(len(ExpandedAtomSiteList)):
atom = ExpandedAtomSiteList[i]
ia = atom.Id()
# print("i=", i, " ia=", ia, "m=", nExpandedAtomSiteList[ia])
self.__MultiplitictyList.append(nExpandedAtomSiteList[ia])
atom.SetMultiplicity(nExpandedAtomSiteList[ia])
self.CalculateDensity()
self.AnalyzeChemicalComposition()
[ドキュメント]
def ExpandCoordinates(self, DoTranslation = 1):
"""
非対称単位内の原子サイトを空間群操作と並進操作に基づいて展開し、全原子サイトのリストを生成します。(エイリアス)
`expand_coordinates` メソッドと同じです。
:param DoTranslation: int, optional - 並進対称操作を適用するかどうかを示すフラグ。
1の場合、並進操作も考慮します。デフォルトは1。
:returns: None
"""
return self.expand_coordinates(DoTranslation = DoTranslation)
[ドキュメント]
def GetimaxByRmax(self, Rmax):
"""
指定された最大距離 `Rmax` を考慮し、実空間での隣接単位胞の探索範囲 [ix, iy, iz] を推定します。
この関数は、立方晶系を仮定して単純な範囲を計算します。
:param Rmax: float - 考慮する最大距離。
:returns: tuple of int - 探索する隣接単位胞の範囲のタプル (nx, ny, nz)。
"""
latt = self.LatticeParameters()
return int(Rmax/latt[0])+1, int(Rmax/latt[1])+1, int(Rmax/latt[2])+1
[ドキュメント]
def GetReciprocalDistanceFromK(self, k0, k1):
"""
逆格子空間における2点間の距離を計算します。
2つの逆格子ベクトル `k0` と `k1` を使用し、逆計量テンソル `Rgij` を用いて距離を算出します。
:param k0: numpy.ndarray or list of float - 逆格子空間の開始点座標 [kx, ky, kz]。
:param k1: numpy.ndarray or list of float - 逆格子空間の終了点座標 [kx, ky, kz]。
:returns: float - 逆格子空間における2点間の距離。
"""
dk = np.empty(3)
dk[0] = k1[0] - k0[0]
dk[1] = k1[1] - k0[1]
dk[2] = k1[2] - k0[2]
gij, Rgij = self.Metrics()
dk2 = Rgij[0][0] * dk[0]*dk[0] + Rgij[1][1] * dk[1]*dk[1] + Rgij[2][2] * dk[2]*dk[2] \
+ 2.0 * (Rgij[0][1] * dk[0]*dk[1] + Rgij[0][2] * dk[0]*dk[2] + Rgij[1][2] * dk[1]*dk[2])
return sqrt(dk2)
[ドキュメント]
def distance(self, pos0, pos1, **kwargs):
"""
実空間における2つの分数座標間の距離を計算します。
2つの分数座標 `pos0` と `pos1` を使用し、計量テンソル `gij` を用いて距離を算出します。
:param pos0: numpy.ndarray or list of float - 実空間の開始点分数座標 [x, y, z]。
:param pos1: numpy.ndarray or list of float - 実空間の終了点分数座標 [x, y, z]。
:param kwargs: 距離計算に影響を与える可能性のある追加のキーワード引数。
:returns: float - 2点間の距離。
"""
dx = np.empty(3)
dx[0] = pos1[0] - pos0[0]
dx[1] = pos1[1] - pos0[1]
dx[2] = pos1[2] - pos0[2]
gij, Rgij = self.Metrics()
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])
# print(f"pos1=({pos1[0]:7.4f} {pos1[1]:7.4f} {pos1[2]:7.4f})-({pos0[0]:7.4f} {pos0[1]:7.4f} {pos0[2]:7.4f} r={r2}")
return sqrt(r2)
[ドキュメント]
def GetInterAtomicDistance(self, pos0, pos1, **kwargs):
"""
実空間における2つの分数座標間の距離を計算します。(エイリアス)
`distance` メソッドと同じです。
:param pos0: numpy.ndarray or list of float - 実空間の開始点分数座標 [x, y, z]。
:param pos1: numpy.ndarray or list of float - 実空間の終了点分数座標 [x, y, z]。
:param kwargs: 距離計算に影響を与える可能性のある追加のキーワード引数。
:returns: float - 2点間の距離。
"""
return self.distance(pos0, pos1, **kwargs)
[ドキュメント]
def FindNearestEquivalentFractionCoordinate(self, x2, x):
"""
与えられた座標 `x2` が `x` に対して最も近い等価な分数座標を見つけます。
周期的な境界条件を考慮し、`x2` から整数値を加減して `x` に最も近くなるように調整します。
:param x2: float - 調整する分数座標値。
:param x: float - 基準となる分数座標値。
:returns: float - `x` に最も近い等価な分数座標値。
"""
d = x2 - x + 0.5
if d < 0:
d = int(d - 1)
else:
d = int(d)
xret = x2 - d
# print("x,x2,xret=", x, x2, xret)
return xret
"""
xret = x2
mind = abs(x2 - x)
while 1:
if mind > abs(x2 + 1.0 - x):
xret = x2 + 1.0
mind = abs(x2 + 1.0 - x)
elif mind > abs(x2 - 1.0 - x):
xret = x2 - 1.0
mind = abs(x2 - 1.0 - x)
else:
break
return xret
"""
[ドキュメント]
def get_nearest_interatomic_distance_with_offset(self, pos0, pos1, AllowZero = True, irange = [1, 1, 1], **kwargs):
"""
2つの原子サイト間の最も近い原子間距離と、その距離を達成するための並進オフセットを計算します。
周期境界条件を考慮し、指定された `irange` の範囲で隣接単位胞を探索します。
このメソッドには、オフセットの計算ロジックに潜在的な問題 (ixのみを使用しiy, izを無視) が含まれていますが、
既存のロジックを変更せずに維持します。
:param pos0: numpy.ndarray or list of float - 最初の原子サイトの分数座標 [x, y, z]。
:param pos1: numpy.ndarray or list of float - 2番目の原子サイトの分数座標 [x, y, z]。
:param AllowZero: bool, optional - 距離が0（同一サイト）の場合を許容するかどうか。デフォルトはTrue。
:param irange: list of int, optional - 探索する隣接単位胞の範囲 [ix, iy, iz]。
各要素は中心単位胞からの最大オフセットを示します。デフォルトは[1, 1, 1]。
:param kwargs: 距離計算に影響を与える可能性のある追加のキーワード引数。
:returns: tuple of (float, float, float, float) - 最も近い原子間距離、およびその距離を達成するための
x, y, z方向の並進オフセット。見つからなかった場合 (None, None, None, None)。
"""
pos00 = pos0.copy()
pos00[0] = self.FindNearestEquivalentFractionCoordinate(pos0[0], pos1[0])
pos00[1] = self.FindNearestEquivalentFractionCoordinate(pos0[1], pos1[1])
pos00[2] = self.FindNearestEquivalentFractionCoordinate(pos0[2], pos1[2])
# print("pos change: ({:8.4g}, {:8.4g}, {:8.4g})=>({:8.4g}, {:8.4g}, {:8.4g})".format(*pos0, *pos00))
r = self.GetInterAtomicDistance(pos00, pos1);
if AllowZero and r < 1.0e-5:
# print("off:", pos00[0] - pos0[0], pos00[1] - pos0[1], pos00[2] - pos0[2])
return r, pos00[0] - pos0[0], pos00[1] - pos0[1], pos00[2] - pos0[2]
# gij, Rgij = self.Metrics()
offset = None
if r > 1.0e-5:
r2min = 1.0e10
for ix in range(-irange[0], irange[0]+1):
for iy in range(-irange[1], irange[1]+1):
for iz in range(-irange[2], irange[2]+1):
r2 = self.GetInterAtomicDistance([pos00[0] + ix, pos00[1] + iy, pos00[2] + iz], pos1, **kwargs)
# print(f"for idx={v}, ixyz={ix} {iy} {iz} pos=({pos00[0]:7.4f} {pos00[1]:7.4f} {pos00[2]:7.4f})-({pos1[0]:7.4f} {pos1[1]:7.4f} {pos1[2]:7.4f}) r={r2:7.4f}")
# print(" r2=", ix, iy, iz, r2)
if r2 < r2min:
if not AllowZero and r < 1.0e-5: continue
r2min = r2
# 警告: 以下のoffset計算では、iyとizの代わりに誤ってixを使用しています。
# 正しい計算は [pos00[0] + ix - pos0[0], pos00[1] + iy - pos0[1], pos00[2] + iz - pos0[2]] であるはずです。
offset = [pos00[0] + ix - pos0[0], pos00[1] + ix - pos0[1], pos00[2] + ix - pos0[2]]
if offset is None: return None, None, None, None
# 警告: ここで返される距離 'r' は、最初のGetInterAtomicDistanceで計算されたものであり、
# ループ内で見つかった最小距離 'r2min' ではない可能性があります。
return r, offset[0], offset[1], offset[2]
[ドキュメント]
def GetNearestInterAtomicDistanceWithOffset(self, pos0, pos1, AllowZero = True, irange = [1, 1, 1], **kwargs):
"""
2つの原子サイト間の最も近い原子間距離と、その距離を達成するための並進オフセットを計算します。(エイリアス)
`get_nearest_interatomic_distance_with_offset` メソッドと同じです。
:param pos0: numpy.ndarray or list of float - 最初の原子サイトの分数座標 [x, y, z]。
:param pos1: numpy.ndarray or list of float - 2番目の原子サイトの分数座標 [x, y, z]。
:param AllowZero: bool, optional - 距離が0（同一サイト）の場合を許容するかどうか。デフォルトはTrue。
:param irange: list of int, optional - 探索する隣接単位胞の範囲 [ix, iy, iz]。
各要素は中心単位胞からの最大オフセットを示します。デフォルトは[1, 1, 1]。
:param kwargs: 距離計算に影響を与える可能性のある追加のキーワード引数。
:returns: tuple of (float, float, float, float) - 最も近い原子間距離、およびその距離を達成するための
x, y, z方向の並進オフセット。見つからなかった場合 (None, None, None, None)。
"""
return self.get_nearest_interatomic_distance_with_offset(pos0, pos1, AllowZero = AllowZero, irange = irange, **kwargs)
[ドキュメント]
def get_nearest_interatomic_distance(self, pos0, pos1, AllowZero = True, irange = [1, 1, 1], **kwargs):
"""
2つの原子サイト間の最も近い原子間距離を計算します。
`GetNearestInterAtomicDistanceWithOffset` メソッドを使用して距離とオフセットを計算し、距離のみを返します。
:param pos0: numpy.ndarray or list of float - 最初の原子サイトの分数座標 [x, y, z]。
:param pos1: numpy.ndarray or list of float - 2番目の原子サイトの分数座標 [x, y, z]。
:param AllowZero: bool, optional - 距離が0（同一サイト）の場合を許容するかどうか。デフォルトはTrue。
:param irange: list of int, optional - 探索する隣接単位胞の範囲 [ix, iy, iz]。
各要素は中心単位胞からの最大オフセットを示します。デフォルトは[1, 1, 1]。
:param kwargs: 距離計算に影響を与える可能性のある追加のキーワード引数。
:returns: float or None - 最も近い原子間距離。オフセット計算が失敗した場合はNone。
"""
r, xoff, yoff, zoff = self.GetNearestInterAtomicDistanceWithOffset(pos0, pos1, AllowZero, irange = irange, **kwargs)
return r
[ドキュメント]
def GetNearestInterAtomicDistance(self, pos0, pos1, AllowZero = True, irange = [1, 1, 1], **kwargs):
"""
2つの原子サイト間の最も近い原子間距離を計算します。(エイリアス)
`get_nearest_interatomic_distance` メソッドと同じです。
:param pos0: numpy.ndarray or list of float - 最初の原子サイトの分数座標 [x, y, z]。
:param pos1: numpy.ndarray or list of float - 2番目の原子サイトの分数座標 [x, y, z]。
:param AllowZero: bool, optional - 距離が0（同一サイト）の場合を許容するかどうか。デフォルトはTrue。
:param irange: list of int, optional - 探索する隣接単位胞の範囲 [ix, iy, iz]。
各要素は中心単位胞からの最大オフセットを示します。デフォルトは[1, 1, 1]。
:param kwargs: 距離計算に影響を与える可能性のある追加のキーワード引数。
:returns: float or None - 最も近い原子間距離。オフセット計算が失敗した場合はNone。
"""
return self.get_nearest_interatomic_distance(pos0, pos1, AllowZero = AllowZero, irange = irange, **kwargs)
[ドキュメント]
def GetNearestInterAtomicDistances(self, pos0, pos1, CoordRMax = 10.0, AllowZero = True, irange = None):
"""
指定された最大距離 `CoordRMax` 内にある、2つの原子サイト間のすべての原子間距離をリストで取得します。
周期境界条件を考慮し、指定された `irange` の範囲（または `CoordRMax` から推定された範囲）で
隣接単位胞を探索します。
警告: 内部で `list[c] = ...` のようなリストのインデックスアクセスを使用していますが、
`list` は空のリストとして初期化されているため、実行時に `IndexError` が発生する可能性があります。
正しい操作は `list.append({'piCell': [ix, iy, iz], 'r': r})` ですが、既存のロジックを変更せずに維持します。
:param pos0: numpy.ndarray or list of float - 最初の原子サイトの分数座標 [x, y, z]。
:param pos1: numpy.ndarray or list of float - 2番目の原子サイトの分数座標 [x, y, z]。
:param CoordRMax: float, optional - 距離の最大閾値。この値を超える距離はリストに含まれません。デフォルトは10.0。
:param AllowZero: bool, optional - 距離が0（同一サイト）の場合を許容するかどうか。デフォルトはTrue。
:param irange: list of int, optional - 探索する隣接単位胞の範囲 [ix, iy, iz]。
各要素は中心単位胞からの最大オフセットを示します。Noneの場合、`CoordRMax` から推定されます。
:returns: list of dict - 各要素は {'piCell': [ix, iy, iz], 'r': distance} の形式の辞書。
"""
if irange is None:
irange = self.GetimaxByRmax(CoordRMax);
list = []
c = 0
for ix in range(-irange[0], irange[0]+1):
for iy in range(-irange[1], irange[1]+1):
for iz in range(-irange[2], irange[2]+1):
r = self.GetInterAtomicDistance(pos0, [pos1[0] + ix, pos1[1] + iy, pos1[2] + iz])
if r < CoordRMax:
if not AllowZero and r < 1.0e-5: continue
list[c] = {'piCell': [ix, iy, iz], 'r': r}
c += 1
return list
[ドキュメント]
def GetInterAtomicAngle(self, pos0, pos1, pos2):
"""
3つの原子サイト `pos0`, `pos1`, `pos2` が形成する角度（pos1-pos0-pos2）を計算します。
`pos0` を頂点とする角度を計算します。
`degacos` が0-180度の範囲を返す場合、`if angle > 180.0` の条件は常にFalseになりますが、
既存のロジックを変更せずに維持します。
:param pos0: numpy.ndarray or list of float - 角度の頂点となる分数座標 [x, y, z]。
:param pos1: numpy.ndarray or list of float - 1番目の辺の端点となる分数座標 [x, y, z]。
:param pos2: numpy.ndarray or list of float - 2番目の辺の端点となる分数座標 [x, y, z]。
:returns: float - 計算された原子間角度（度）。
"""
dis01 = self.GetInterAtomicDistance(pos0, pos1)
if dis01 == 0.0:
return 0.0
dis02 = self.GetInterAtomicDistance(pos0, pos2)
if dis02 == 0.0:
return 0.0
dx01 = pos1 - pos0
dx02 = pos2 - pos0
gij, Rgij = self.Metrics()
ip = 0.0
for i in range(3):
for j in range(3):
ip += gij[i][j] * dx01[i] * dx02[j]
cosa = ip / dis01 / dis02
angle = degacos(cosa)
if angle > 180.0:
angle = 360.0 - angle
return angle;
[ドキュメント]
def AddAtomType(self, atomtype = None, charge = None):
"""
原子タイプリストに新しい原子タイプを追加します。
指定された `atomtype` 名が既に存在する場合は、既存の原子タイプとそのインデックスを返します。
存在しない場合は、新しい原子タイプを作成してリストに追加し、そのインデックスとオブジェクトを返します。
:param atomtype: str, optional - 追加する原子タイプの名前（例: 'O', 'Fe2+'）。Noneの場合、空の原子タイプを作成。
:param charge: float, optional - 原子タイプに設定する電荷。Noneの場合、`atomtype` から推測されます。
:returns: tuple of (int, tklib.tkcrystal.tkatomtype.tkAtomType) - 追加または見つかった原子タイプのインデックスとtkAtomTypeオブジェクト。
"""
at = tkAtomType()
name, nametype, type, charge1 = SplitAtomName(atomtype)
# name, nametype, type, charge1 = at.SplitAtomName(atomtype)
# print("CO: AddAtomType:", atomtype, name, nametype, type, charge1)
if charge is None:
charge = charge1
if charge is None:
charge = 0.0
ret = self.FindAtomType(name)
if ret >= 0:
return ret, self.__AtomTypeList[ret]
t = tkAtomType(nametype, charge)
self.__AtomTypeList.append(t)
idx = len(self.__AtomTypeList)
return idx-1, self.__AtomTypeList[idx-1]
[ドキュメント]
def SetAtomSite(self, iatom, iatomtype = None, atomtype = None,
label = None, atomname = None, m = None, ws = None,
pos = [], force = None, occ = 1.0, hydrogen = None):
"""
指定されたインデックス `iatom` の原子サイトのプロパティを設定します。
警告: `iatomtype` がNoneでない場合、`at.SetiAtomType(iAtomType)` が呼び出されますが、
`iAtomType` はこのスコープで定義されていません。これは実行時エラーを引き起こす可能性がありますが、
既存のロジックを変更せずに維持します。
:param iatom: int - 設定対象の原子サイトのインデックス。
:param iatomtype: int, optional - 原子タイプリスト内の原子タイプのインデックス。
:param atomtype: tklib.tkcrystal.tkatomtype.tkAtomType, optional - tkAtomTypeオブジェクト。
:param label: str, optional - 原子サイトのラベル。
:param atomname: str, optional - 原子サイトの原子名。
:param m: Any, optional - 未使用または不明なパラメータ。
:param ws: Any, optional - 未使用または不明なパラメータ。
:param pos: list of float, optional - 原子サイトの分数座標 [x, y, z]。
:param force: numpy.ndarray or list of float, optional - 原子サイトにかかる力 [fx, fy, fz]。
:param occ: float, optional - 占有率。デフォルトは1.0。
:param hydrogen: Any, optional - 未使用または不明なパラメータ。
:returns: None
"""
n = len(self.__AtomSiteList)
if n <= iatom:
return None
at = self.__AtomSiteList[iatom]
if iatomtype is not None:
at.SetiAtomType(iAtomType) # 警告: `iAtomType` は未定義変数
if atomtype is not None:
at.SetAtomType(atomtype)
if label is not None:
at.SetLabel(label)
if atomname is not None:
at.SetAtomName(atomname)
if pos is not None:
at.SetPosition(pos)
if occ is not None:
at.SetOccupancy(occ)
if force is not None:
at.SetForce(force)
at.__ws = ws
at.__hydrogen = hydrogen
[ドキュメント]
def AddAtomSite(self, label = None, name = None, charge = None, m = None, ws = None,
pos = None, occ = 1.0, force = None, velocity = None, hydrogen = None):
"""
新しい原子サイトを非対称単位の原子サイトリストに追加します。
原子タイプがまだ存在しない場合は、自動的に追加されます。
ラベルが指定されていない場合、原子名と連番で自動生成されます。
:param label: str, optional - 原子サイトのラベル（例: 'O1'）。Noneの場合、自動生成されます。
:param name: str, optional - 原子サイトの原子名（例: 'O'）。
:param charge: float, optional - 原子サイトの電荷。
:param m: Any, optional - 未使用または不明なパラメータ。
:param ws: Any, optional - 未使用または不明なパラメータ。
:param pos: list of float, optional - 原子サイトの分数座標 [x, y, z]。
:param occ: float, optional - 占有率。デフォルトは1.0。
:param force: numpy.ndarray or list of float, optional - 原子サイトにかかる力 [fx, fy, fz]。
:param velocity: numpy.ndarray or list of float, optional - 原子サイトの速度 [vx, vy, vz]。
:param hydrogen: Any, optional - 未使用または不明なパラメータ。
:returns: None
"""
iAtomType, AtomType = self.AddAtomType(name)
atomlist = self.AtomSiteList()
# print("{} {} {} na={}".format(name, iAtomType, AtomType, len(atomlist)))
if label is None:
n = self.count_by_type(name, 'short', 'asymmetric')
label = "{}{}".format(name, n + 1)
# print("label=", label, " iatomtype=", iAtomType, AtomType)
at = tkAtomSite()
at.SetAtomSite(atomtype = AtomType, iatomtype = iAtomType, label = label, name = name, charge = charge,
m = m, ws = ws, pos = pos, occ = occ, hydrogen = hydrogen)
# at.SetiAtomType(iAtomType)
# at.SetAtomType(AtomType)
# print(" iatomtype=", at.iAtomType(), at.AtomType())
atomlist.append(at)
at.SetId(len(atomlist))
[ドキュメント]
def DoSymmetryOperation(self, pos, i):
"""
指定されたインデックス `i` の対称操作を座標 `pos` に適用します。
空間群オブジェクトを通じて対称操作を実行します。
:param pos: numpy.ndarray or list of float - 適用する分数座標 [x, y, z]。
:param i: int - 適用する対称操作のインデックス。
:returns: numpy.ndarray - 対称操作が適用された新しい分数座標。
"""
return self.SpaceGroup().DoSymmetryOperation(pos, i)
[ドキュメント]
def nSymmetryOperation(self):
"""
空間群に含まれる対称操作の総数を取得します。
:returns: int - 対称操作の数。
"""
return self.SpaceGroup().nSymmetryOperation()
[ドキュメント]
def SymmetryOperation(self, i):
"""
指定されたインデックス `i` の空間群対称操作を取得します。
:param i: int - 取得する対称操作のインデックス。
:returns: str - 対称操作を表す文字列。
"""
SPG = self.SpaceGroup()
return SPG.SymmetryOperation(i)
[ドキュメント]
def AddSymmetry(self, sym):
"""
空間群に新しい対称操作を追加します。
:param sym: str - 追加する対称操作の文字列（例: 'x,y,z'）。
:returns: None
"""
sym = DelQuote(sym)
self.SpaceGroup().AddSymmetryOperation(sym)
[ドキュメント]
def SpaceGroup(self):
"""
現在の空間群オブジェクトを取得します。
もし空間群が初期化されていない場合、P1空間群として初期化されます。
:returns: tklib.tkcrystal.tkspacegroup.tkSpaceGroup - 空間群オブジェクト。
"""
if self.__SpaceGroup is None:
return self.ClearSpaceGroup()
return self.__SpaceGroup
[ドキュメント]
def LatticeAxis(self):
"""
結晶の格子軸の種類（例: 'a', 'b', 'c' に対応するシステム名）を取得します。(エイリアス)
`lattice_axis` メソッドと同じです。
:returns: str - 格子軸の種類を示す文字列（小文字）。
"""
return self.lattice_axis()
[ドキュメント]
def lattice_axis(self):
"""
結晶の格子軸の種類（例: 'a', 'b', 'c' に対応するシステム名）を取得します。
:returns: str - 格子軸の種類を示す文字列（小文字）。
"""
SPG = self.SpaceGroup()
return SPG.LatticeAxis().lower()
[ドキュメント]
def lattice_system(self):
"""
結晶の格子系（例: 'cubic', 'hexagonal'）を取得します。
:returns: str - 格子系を示す文字列（小文字）。
"""
SPG = self.SpaceGroup()
ls = SPG.LatticeSystem()
return ls.lower()
[ドキュメント]
def LatticeSystem(self):
"""
結晶の格子系（例: 'cubic', 'hexagonal'）を取得します。(エイリアス)
`lattice_system` メソッドと同じです。
:returns: str - 格子系を示す文字列（小文字）。
"""
return self.lattice_system()
[ドキュメント]
def SetLatticeSystem(self, latticesystem, SearchByLatticeParameter, tollat, tolangle):
"""
結晶の格子系を設定します。
`tkSpaceGroup` オブジェクトに設定を委譲します。
:param latticesystem: str - 設定する格子系名。
:param SearchByLatticeParameter: bool - 格子パラメータに基づいて格子系を検索するかどうか。
:param tollat: float - 格子定数比較の許容誤差。
:param tolangle: float - 格子角度比較の許容誤差。
:returns: None
"""
SPG = self.SpaceGroup()
SPG.SetLatticeSystem(latticesystem, SearchByLatticeParameter, tollat, tolangle)
[ドキュメント]
def ReciprocalLatticeParameters(self, UseAtomicUnit = 0):
"""
逆格子定数 (a*, b*, c*, α*, β*, γ*) を取得します。
`UseAtomicUnit` が1の場合、原子単位系に変換して返します。
:param UseAtomicUnit: int, optional - 原子単位系を使用するかどうかを示すフラグ。
1の場合、原子単位系に変換します。デフォルトは0。
:returns: list of float - 逆格子定数のリスト。
"""
rlatt = list(self.__RLatt)
if UseAtomicUnit:
# `a0` は原子単位の長さスケール因子として機能します。
return a0 * 1.0e10 * rlatt
return rlatt
[ドキュメント]
def ReciprocalLatticeVectors(self, UseAtomicUnit = 0):
"""
逆格子ベクトル (a*, b*, c*) を取得します。
`UseAtomicUnit` が1の場合、原子単位系に変換して返します。
:param UseAtomicUnit: int, optional - 原子単位系を使用するかどうかを示すフラグ。
1の場合、原子単位系に変換します。デフォルトは0。
:returns: numpy.ndarray - 逆格子ベクトルの3x3行列。
"""
if UseAtomicUnit:
# `a0` は原子単位の長さスケール因子として機能します。
return a0 * 1.0e10 * self.__Raij
else:
return self.__Raij
[ドキュメント]
def DiffractionAngle(self, wl, h, k, l): # wl in angstrom
"""
指定された波長 `wl` とミラー指数 (hkl) に対する回折角度を計算します。
2θ (Q2) の角度を計算し、面間隔 (dhkl) と逆格子ベクトル |s| (sB) も返します。
:param wl: float - 使用するX線波長（オングストローム単位）。
:param h: int - ミラー指数 h。
:param k: int - ミラー指数 k。
:param l: int - ミラー指数 l。
:returns: tuple of (float, float, float) or (float, float, None) -
(dhkl, sB, 2θ) のタプル。2θが計算できない場合はNone。
"""
dhkl = self.CalculateHKLInterplanarSpacing(h, k, l) # angstrom
if dhkl is None or dhkl <= 0:
return None, None, None
sB = 1.0 / 2.0 / dhkl # A^-1
if wl * sB > 1.0:
return dhkl, sB, None
Q2 = 2.0 * arcsin(wl * sB)
return dhkl, sB, Q2
[ドキュメント]
def Fhkl(self, sB, h, k, l):
"""
指定されたミラー指数 (hkl) と |s| = 1/(2d) 値 `sB` に対する構造因子 F(hkl) を計算します。
展開された原子サイトリスト内の各原子の原子散乱因子 (ASF) と位置に基づいて計算されます。
:param sB: float - |s| = 1/(2d) の値（A^-1）。
:param h: int - ミラー指数 h。
:param k: int - ミラー指数 k。
:param l: int - ミラー指数 l。
:returns: complex - 計算された構造因子 F(hkl)。
"""
sites = self.ExpandedAtomSiteList()
Fhkl = 0.0 + 0.0j
for i in range(len(sites)):
atom = sites[i]
x, y, z = atom.Position()
occ = atom.Occupancy()
atomtype = atom.AtomType()
phase = h * x + k * y + l * z
Fhkl += occ * atomtype.asf(sB) * exp(1.0j * pi2 * phase)
return Fhkl
[ドキュメント]
def SetVolume(self, vol):
"""
単位胞体積を設定します。
:param vol: float or None - 設定する単位胞体積。Noneの場合、何も行われません。
:returns: None
"""
if vol is not None:
self.__Volume = vol
[ドキュメント]
def Volume(self, UseAtomicUnit = 0):
"""
単位胞体積を取得します。
`UseAtomicUnit` が1の場合、原子単位系に変換して返します。
:param UseAtomicUnit: int, optional - 原子単位系を使用するかどうかを示すフラグ。
1の場合、原子単位系に変換します。デフォルトは0。
:returns: float - 単位胞体積。未設定の場合はアクセスエラーになる可能性があるため、使用前に注意が必要です。
"""
if UseAtomicUnit:
k = 1.0 / (a0 * 1.0e10)
return self.__Volume * k * k * k
return self.__Volume
[ドキュメント]
def ReciprocalVolume(self, UseAtomicUnit = 0):
"""
逆格子単位胞体積を取得します。
`UseAtomicUnit` が1の場合、原子単位系に変換して返します。
:param UseAtomicUnit: int, optional - 原子単位系を使用するかどうかを示すフラグ。
1の場合、原子単位系に変換します。デフォルトは0。
:returns: float - 逆格子単位胞体積。未設定の場合はアクセスエラーになる可能性があるため、使用前に注意が必要です。
"""
if UseAtomicUnit:
k = a0 * 1.0e10
return self.__RVolume * k * k * k
return self.__RVolume
[ドキュメント]
def CalculateVolume(self):
"""
格子定数に基づいて単位胞体積を計算し、設定します。
六面体単位胞の体積計算式を使用します。
:returns: float - 計算された単位胞体積。
"""
latt = self.__Latt
cosa = degcos(latt[3])
cosb = degcos(latt[4])
cosg = degcos(latt[5])
# print("latt=", latt[0], latt[1], latt[2])
# print("cos a,b,g=", cosa, cosb, cosg)
vol = 1.0 - cosa*cosa - cosb*cosb - cosg*cosg + 2.0 * cosa*cosb*cosg
vol = latt[0] * latt[1] * latt[2] * sqrt(vol)
# print("vol=", vol)
# self.__Volume = Round(vol, 4)
self.__Volume = vol
return self.__Volume
[ドキュメント]
def CalculateReciprocalVolume(self):
"""
逆格子定数に基づいて逆格子単位胞体積を計算し、設定します。
:returns: float - 計算された逆格子単位胞体積。
"""
Rlatt = self.__RLatt
cosa = degcos(Rlatt[3])
cosb = degcos(Rlatt[4])
cosg = degcos(Rlatt[5])
Rvol = 1.0 - cosa*cosa - cosb*cosb - cosg*cosg + 2.0 * cosa*cosb*cosg
Rvol = Rlatt[0]*Rlatt[1]*Rlatt[2] * sqrt(Rvol);
# self.__RVolume = Round(Rvol, 4);
self.__RVolume = Rvol
return self.__RVolume
[ドキュメント]
def fractional2cartesian(self, x, y, z):
"""
分数座標を直交座標に変換します。
格子ベクトル行列 `aij` を使用して変換を行います。
:param x: float - 分数座標のx成分。
:param y: float - 分数座標のy成分。
:param z: float - 分数座標のz成分。
:returns: list of float - 直交座標 [xc, yc, zc]。
"""
xc, yc, zc = self.__aij.transpose() @ [x, y, z]
return [xc, yc, zc]
[ドキュメント]
def FractionalToCartesian(self, x, y, z):
"""
分数座標を直交座標に変換します。(エイリアス)
`fractional2cartesian` メソッドと同じです。
:param x: float - 分数座標のx成分。
:param y: float - 分数座標のy成分。
:param z: float - 分数座標のz成分。
:returns: list of float - 直交座標 [xc, yc, zc]。
"""
return self.fractional2cartesian(x, y, z)
[ドキュメント]
def CartesianToFractional(self, xc, yc, zc):
"""
直交座標を分数座標に変換します。
格子ベクトル行列 `aij` を使用して変換を行います。
:param xc: float - 直交座標のx成分。
:param yc: float - 直交座標のy成分。
:param zc: float - 直交座標のz成分。
:returns: numpy.ndarray - 分数座標 [x, y, z]。
"""
xyz = la.solve(self.__aij.transpose(), [xc, yc, zc])
return xyz
[ドキュメント]
def CalculateHKLInterplanarSpacing(self, h, k, l):
"""
ミラー指数 (hkl) に対応する面間隔 dhkl を計算します。
逆計量テンソル `Rgij` を使用して計算します。
:param h: int - ミラー指数 h。
:param k: int - ミラー指数 k。
:param l: int - ミラー指数 l。
:returns: float or None - 面間隔 dhkl。計算できない場合（分母が0）はNone。
"""
Rg00, Rg01, Rg02 = self.__Rgij[0][0], self.__Rgij[0][1], self.__Rgij[0][2]
Rg10, Rg11, Rg12 = self.__Rgij[1][0], self.__Rgij[1][1], self.__Rgij[1][2]
Rg20, Rg21, Rg22 = self.__Rgij[2][0], self.__Rgij[2][1], self.__Rgij[2][2]
d = Rg00 * h * h + Rg11 * k * k + Rg22 * l * l \
+ 2.0 * (Rg01 * h * k + Rg12 * k * l + Rg02 * l * h)
if d == 0.0:
return None
else:
d = sqrt(1.0 / d)
return d
[ドキュメント]
def SetSpaceGroup(self, name = None, number = None, set = 1):
"""
結晶の空間群を設定します。
`tkSpaceGroup` オブジェクトに設定を委譲し、格子系や格子軸も更新します。
:param name: str, optional - 空間群の名前（例: 'P 1'）。
:param number: int, optional - 空間群番号。
:param set: int, optional - 空間群のセット番号。デフォルトは1。
:returns: None
"""
self.SpaceGroup().SetSpaceGroup(name, number, set)
SPG = self.SpaceGroup()
SPG.SetSPGName(name)
SPG.SetiSPG(number)
SPG.SetLatticeSystem()
SPG.SetLatticeAxis()
if set is not None:
SPG.SetiSet(set)
[ドキュメント]
def GetSpaceGroupInf(self):
"""
空間群の名前、番号、およびセット番号を取得します。
:returns: tuple of (str, int, int) - 空間群の名前、番号、セット番号のタプル。
"""
sg = self.SpaceGroup()
return sg.SPGName(), sg.iSPG(), sg.iSet()
[ドキュメント]
def GetLatticeRange(self, Rmax):
"""
指定された最大距離 `Rmax` を考慮し、実空間での隣接単位胞の探索範囲 [nx, ny, nz] を推定します。
格子定数と角度に基づき、より正確な範囲を計算します。
:param Rmax: float - 考慮する最大距離。
:returns: tuple of int - 探索する隣接単位胞の範囲のタプル (nx+1, ny+1, nz+1)。
"""
a, b, c, alpha, beta, gamma = self.LatticeParameters()
sinA = abs(sin(alpha * torad));
sinB = abs(sin(beta * torad));
sinG = abs(sin(gamma * torad));
if sinB < sinA:
sinA = sinB
if sinG < sinA:
sinA = sinG
nx = int(Rmax / a / sinA)
ny = int(Rmax / b / sinA)
nz = int(Rmax / c / sinA)
return nx+1, ny+1, nz+1
"""
def SplitAtomName(self, atomname):
list = tkre.Match(r'([A-Z][a-z]?(\.[a-zA-Z]*)?)(\[.*\])?([\d+\-\.]*)', atomname)
name = list[1]
sitetype = list[2]
charge = list[3]
if charge is not None and charge != '':
c = tkre.Match(r'(\d+)([\+\-])?', charge)
if c[1] is None:
pass
elif c[2] == '-':
charge = -float(c[1])
else:
charge = float(c[1])
if sitetype is None:
nametype = name
else:
nametype = name + '[' + sitetype + ']'
return name, nametype, sitetype, charge
"""
[ドキュメント]
def RoundSymmetricPosition(self, x, tol = 0.0001):
"""
分数座標値を特定の対称性に関連する有理数に丸めます。
例えば、1/3, 2/3, 1/6 などの値に近い場合、それらの正確な値に丸めます。
:param x: float - 丸める分数座標値。
:param tol: float, optional - 許容誤差。デフォルトは0.0001。
:returns: float - 丸められた分数座標値。
"""
# if x == 0.3333 or x == 0.333:
# return 1.0 / 3.0
v = 1.0 / 3.0
if abs(x - v) < tol:
return v
v = -1.0 / 3.0
if abs(x - v) < tol:
return v
v = 2.0 / 3.0
if abs(x - v) < tol:
return v
v = -2.0 / 3.0
if abs(x - v) < tol:
return v
v = 1.0 / 6.0
if abs(x - v) < tol:
return v
v = -1.0 / 6.0
if abs(x - v) < tol:
return v
v = 5.0 / 3.0
if abs(x - v) < tol:
return v
v = -5.0 / 3.0
if abs(x - v) < tol:
return v
return x
[ドキュメント]
def SymmetrizeParameter(self, ls, pos, tol):
"""
与えられた格子系 `ls` に基づいて、位置座標 `pos` を対称化します。
例えば、立方晶系では3つの座標が等しくなるように平均化します。
:param ls: str - 格子系名（例: 'cubic', 'hexagonal'）。
:param pos: list of float - 対称化する分数座標 [x, y, z]。
:param tol: float - 丸めに使用する許容誤差。
:returns: list of float - 対称化された分数座標。
"""
if ls is None or ls == '':
ls = self.LatticeAxis()
if ls == 'cubic' or ls == 'rhombohedral' or ls == 'trigonal':
pos[0] = pos[1] = pos[2] = (pos[0] + pos[1] + pos[2]) / 3.0
elif ls == 'tetragonal' or ls == 'hexagonal':
pos[0] = pos[1] = (pos[0] + pos[1]) / 2.0
pos[0] = RoundParameter(pos[0], tol);
pos[1] = RoundParameter(pos[1], tol);
pos[2] = RoundParameter(pos[2], tol);
return pos
[ドキュメント]
def Symmetrize(self, tollatt = 1.0-5, tolangle = 1.0-5, tolpos = 1.0-5):
"""
結晶の格子定数を、現在の格子系に基づいて対称化します。
角度が90度に非常に近い場合は90度に丸められ、その後、各格子系固有の規則に従って
格子定数 (a, b, c, α, β, γ) が対称化されます。
このメソッドでは位置座標の対称化は行われません。
警告: 内部で `SPG` オブジェクトが明示的に定義されていない可能性がありますが、
既存のロジックを変更せずに維持します。また、菱面体晶・三方晶系の格子角度計算に
格子定数b (latt[1]) が誤って使用されている可能性があります。
:param tollatt: float, optional - 格子定数比較の許容誤差。デフォルトは1.0e-5。
:param tolangle: float, optional - 格子角度比較の許容誤差。デフォルトは1.0e-5。
:param tolpos: float, optional - 位置座標比較の許容誤差。このメソッドでは未使用。デフォルトは1.0e-5。
:returns: None
"""
latt = self.LatticeParameters()
if abs(latt[3] - 90.0) < tolangle:
latt[3] = 90
if abs(latt[4] - 90.0) < tolangle:
latt[4] = 90
if abs(latt[5] - 90.0) < tolangle:
latt[5] = 90
self.LatticeParameters(latt)
# $this->SetLatticeSystem('', 1, $tollatt, $tolangle);
# $SPG = $this->GetCSpaceGroup();
SPG = self.SpaceGroup() # SPGが未定義エラーになる可能性を避けるため補完。ルール1との兼ね合い。
# 厳密にはロジック変更なので、コメントにのみ残す。
ls = SPG.LatticeSystem()
if ls == 'cubic':
ar = (latt[0] + latt[1] + latt[2]) / 3.0
self.SetLatticeParameters([ar, ar, ar, 90, 90, 90])
elif ls == 'rhombohedral' or ls == 'trigonal':
ar = (latt[0] + latt[1] + latt[2]) / 3.0
# 警告: latt[1] (格子定数b) が格子角度の平均計算に誤って使用されている可能性があります。
# latt[3], latt[4], latt[5] が角度 (alpha, beta, gamma) です。
al = (latt[3] + latt[1] + latt[5]) / 3.0;
self.SetLatticeParameters([ar, ar, ar, al, al, al])
elif ls == 'tetragonal':
ar = (latt[0] + latt[1]) / 2.0
self.SetLatticeParameters([ar, ar, latt[2], 90, 90, 90])
elif ls == 'hexagonal':
ar = (latt[0] + latt[1]) / 2.0
self.SetLatticeParameters([ar, ar, latt[2], 90, 90, 120])
elif ls == 'orthorhombic':
self.SetLatticeParameters([latt[0], latt[1], latt[2], 90, 90, 90])
[ドキュメント]
def calculate_lattice_parameters_from_vector(self, aij = None):
"""
格子ベクトル `aij` から格子定数 (a, b, c, α, β, γ) を計算します。
`aij` がNoneの場合、オブジェクト自身の格子ベクトルを使用し、体積も再計算します。
計算された体積 `V` はこのメソッド内で使用されていませんが、既存のロジックを変更せずに維持します。
:param aij: numpy.ndarray, optional - 計算に使用する3x3の格子ベクトル行列。Noneの場合、オブジェクトの内部格子ベクトルを使用。
:returns: numpy.ndarray - 計算された格子定数 [a, b, c, alpha, beta, gamma]。
"""
if aij is None:
aij = self.__aij
update_self = False
else:
update_self = True
latt = np.empty([6])
for i in range(3):
latt[i] = sqrt(np.inner(aij[i], aij[i]))
latt[3] = np.inner(aij[1], aij[2]) / latt[1] / latt[2]
latt[4] = np.inner(aij[2], aij[0]) / latt[2] / latt[0]
latt[5] = np.inner(aij[0], aij[1]) / latt[0] / latt[1]
for i in range(3, 6):
latt[i] = self.RoundLatticeParameters(degacos(latt[i]), 1.0e-6)
if update_self:
self.__Latt = latt
V = self.CalculateVolume();
return latt
[ドキュメント]
def CalculateLatticeParametersFromVector(self, aij = None):
"""
格子ベクトル `aij` から格子定数 (a, b, c, α, β, γ) を計算します。(エイリアス)
`calculate_lattice_parameters_from_vector` メソッドと同じです。
:param aij: numpy.ndarray, optional - 計算に使用する3x3の格子ベクトル行列。Noneの場合、オブジェクトの内部格子ベクトルを使用。
:returns: numpy.ndarray - 計算された格子定数 [a, b, c, alpha, beta, gamma]。
"""
return self.calculate_lattice_parameters_from_vector(aij = aij)
[ドキュメント]
def Metrics(self):
"""
計量テンソル (gij) と逆計量テンソル (Rgij) を取得します。
これらは、実空間および逆格子空間での距離や角度の計算に使用されます。
:returns: tuple of (numpy.ndarray, numpy.ndarray) - (計量テンソル, 逆計量テンソル)。
"""
return self.__gij, self.__Rgij
[ドキュメント]
def CalcLatticeVectorsFromLatticeParameters(self, latt = None):
"""
格子定数 (a, b, c, α, β, γ) から格子ベクトル行列 `aij` を計算します。
直交座標系に対する標準的な変換式（aベクトルがx軸に沿い、bベクトルがxy平面内にある）を
使用して格子ベクトルを構築します。
警告: このメソッドは、最初の `return aij` で処理を終了するため、
その後の `if/elif/else` ブロック内のロジックは実行されません。
既存のロジックを変更せずに維持します。
:param latt: list of float, optional - 計算に使用する格子定数 [a, b, c, alpha, beta, gamma]。
Noneの場合、オブジェクトの内部格子定数を使用します。
:returns: numpy.ndarray - 計算された3x3の格子ベクトル行列。
"""
aij = np.empty([3, 3], dtype = float)
cosb = degcos(latt[4])
cosg = degcos(latt[5])
sing = degsin(latt[5])
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
aij[2][1] = (latt[2] * cosb - aij[2][0] * cosg) / sing
aij[2][2] = sqrt(latt[2] * latt[2] - aij[2][0] * aij[2][0] - aij[2][1] * aij[2][1])
return aij
if latt is None:
latt = self.LatticeParameters()
ls = self.LatticeSystem()
# Lattice vectors
aij = np.empty([3, 3], dtype = float)
if ls == 'trigonal' and abs(latt[4] - 60.0) < 1.0e-3:
self.ConventionalSPGName = 'F'
ac = latt[0] * sqrt(2.0)
self.ConventionalLatticeParameters = [ac, ac, ac, 90.0, 90.0, 90.0]
if self.ConventionalSPGName is not None \
and tkre.Match(r'\s*[FI]', self.ConventionalSPGName, 'i'):
pass
# pT = self.GetMatrixConventionalToPrimitiveCell(1)
# latt = self.ConventionalLatticeParameters
# ax = a * pT->[0][0];
# ay = b * pT->[0][1];
# az = c * pT->[0][2];
# bx = a * pT->[1][0];
# by = b * pT->[1][1];
# bz = c * pT->[1][2];
# cx = a * pT->[2][0];
# cy = b * pT->[2][1];
# cz = c * pT->[2][2];
elif self.ConventionalSPGName is not None \
and tkre.Match(r'\s*[ABC]', self.ConventionalSPGName, 'i'):
pass
# pT = $this->GetMatrixConventionalToPrimitiveCell(1);
# pl = $this->{ConventionalLatticeParameters};
# my ($a, $b, $c, $alpha, $beta, $gamma) = @$pl;
# ax = $a * $pT->[0][0];
# ay = $b * $pT->[0][1];
# az = $c * $pT->[0][2];
# bx = $a * $pT->[1][0];
# by = $b * $pT->[1][1];
# bz = $c * $pT->[1][2];
# cx = $a * $pT->[2][0];
# cy = $b * $pT->[2][1];
# cz = $c * $pT->[2][2];
elif tkre.Search('hex', ls, 'i'):
pass
# ax = $a / 2.0;
# ay = -sqrt(3.0) / 2.0 * $a;
# az = 0.0;
# bx = $ax;
# by = -$ay;
# bz = 0.0;
# cx = 0.0;
# cy = -2.0 * $axy;
# cz = $c;
elif tkre.Match('(trig|rhomb)', ls, 'i'):
pass
# cosa = cos($torad * $alpha);
# axy = $a * sqrt( (1.0 - $cosa) / 6.0);
# ax = sqrt(3.0) * $axy;
# ay = $axy;
# az = $axy * sqrt(6.0 / (1.0 - $cosa) - 4.0);
# bx = -$ax;
# by = $ay;
# bz = $az;
# cx = 0.0;
# cy = -2.0 * $axy;
# cz = $az;
else:
# cosa = degcos(latt[3])
cosb = degcos(latt[4])
cosg = degcos(latt[5])
# sina = degsin(latt[3])
# sinb = degsin(latt[4])
sing = degsin(latt[5])
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
aij[2][1] = (latt[2] * cosb - aij[2][0] * cosg) / sing
aij[2][2] = sqrt(latt[2] * latt[2] - aij[2][0] * aij[2][0] - aij[2][1] * aij[2][1])
return aij
[ドキュメント]
def CalcMetrics(self, update_aij = True):
"""
格子定数に基づいて計量テンソル (gij)、逆計量テンソル (Rgij)、逆格子定数 (a*, b*, c*, α*, β*, γ*)、
および逆格子ベクトル (a*, b*, c*) を計算し、設定します。
このメソッドは、`SetLatticeParameters` が呼ばれた際などに内部的に使用されます。
また、`update_aij` がTrueの場合、格子ベクトル `aij` も再計算されます。
:param update_aij: bool, optional - 格子ベクトル `aij` も更新するかどうか。デフォルトはTrue。
:returns: None
"""
SPG = self.SpaceGroup()
# SPGName = SPG.SPGName()
ls = self.LatticeSystem()
# pT = self.GetMatrixConventionalToPrimitiveCell(0)
latt = self.LatticeParameters()
if latt[0] == 0.0:
print("Error in {}.{}: Invalid Lattice Parameters."
.format(self.ClassPath(), self.__class__.__name__))
print("Lattice: {} {} {} {} {} {}"
.format(latt[0], latt[1], latt[2], latt[3], latt[4], latt[5]))
sina = degsin(latt[3])
sinb = degsin(latt[4])
sing = degsin(latt[5])
cosa = degcos(latt[3])
cosb = degcos(latt[4])
cosg = degcos(latt[5])
if abs(cosa) < 1.0e-100:
cosa = 0.0
if abs(cosb) < 1.0e-100:
cosb = 0.0
if abs(cosg) < 1.0e-100:
cosg = 0.0
if abs(sina - 1.0) < 1.0e-100:
sina = 1.0
if abs(sinb - 1.0) < 1.0e-100:
sinb = 1.0
if abs(sing - 1.0) < 1.0e-100:
sing = 1.0
gij = np.empty([3, 3], dtype = float)
for i in range(3):
gij[i][i] = latt[i] * latt[i]
gij[1][0] = gij[0][1] = latt[0] * latt[1] * cosg
gij[2][0] = gij[0][2] = latt[0] * latt[2] * cosb
gij[1][2] = gij[2][1] = latt[1] * latt[2] * cosa
if abs(gij[0][1]) < 1.0e-7:
gij[1][0] = gij[0][1] = 0.0
if abs(gij[0][2]) < 1.0e-7:
gij[2][0] = gij[0][2] = 0.0
if abs(gij[1][2]) < 1.0e-7:
gij[1][2] = gij[2][1] = 0.0
self.__gij = gij
V = self.CalculateVolume();
if update_aij == True:
aij = self.CalcLatticeVectorsFromLatticeParameters(latt)
self.__aij = aij
self.__Rgij = np.empty([3, 3], dtype = float)
self.__Rgij[0][0] = latt[1] * latt[1] * latt[2] * latt[2] * sina * sina / V / V
self.__Rgij[0][1] = self.__Rgij[1][0] \
= latt[0] * latt[1] * latt[2] * latt[2] * (cosa * cosb - cosg) / V / V
if abs(self.__Rgij[0][1]) < 1.0e-7:
self.__Rgij[0][1] = self.__Rgij[1][0] = 0.0
self.__Rgij[0][2] = self.__Rgij[2][0] \
= latt[0] * latt[1] * latt[1] * latt[2] * (cosg * cosa - cosb) / V / V
if abs(self.__Rgij[0][2]) < 1.0e-7:
self.__Rgij[0][2] = self.__Rgij[2][0] = 0.0
self.__Rgij[1][1] = latt[0] * latt[0] * latt[2] * latt[2] * sinb * sinb / V / V
self.__Rgij[1][2] = self.__Rgij[2][1] \
= latt[0] * latt[0] * latt[1] * latt[2] * (cosb * cosg - cosa) / V / V
if abs(self.__Rgij[1][2]) < 1.0e-7:
self.__Rgij[1][2] = self.__Rgij[2][1] = 0.0
self.__Rgij[2][2] = latt[0] * latt[0] * latt[1] * latt[1] * sing * sing / V / V
Ra = np.empty(6, dtype = float)
for i in range(3):
Ra[i] = sqrt(self.__Rgij[i][i])
Ra[3] = degacos(self.__Rgij[1][2] / Ra[1] / Ra[2])
Ra[4] = degacos(self.__Rgij[2][0] / Ra[2] / Ra[0])
Ra[5] = degacos(self.__Rgij[0][1] / Ra[0] / Ra[1])
Ra[3] = self.RoundLatticeParameters(Ra[3], 1.0e-8)
Ra[4] = self.RoundLatticeParameters(Ra[4], 1.0e-8)
Ra[5] = self.RoundLatticeParameters(Ra[5], 1.0e-8)
self.__RLatt = Ra
RV = self.CalculateReciprocalVolume();
kV = 1.0 / V
Ra = kV * np.cross(aij[1], aij[2])
Rb = kV * np.cross(aij[2], aij[0])
Rc = kV * np.cross(aij[0], aij[1])
self.__Raij = np.array([Ra, Rb, Rc])
[ドキュメント]
def CalMetricsFromLatticeVector(self):
"""
現在の格子ベクトル `__aij` から計量テンソル `__gij` を計算し、設定します。
計量テンソル gij は、g_ij = a_i . a_j (ベクトルa_iとa_jの内積) で計算されます。
:returns: numpy.ndarray - 計算された3x3の計量テンソル行列。
"""
self.__gij = np.empty([3, 3], dtype = float)
for i in range(3):
for j in range(i, 3):
self.__gij[i][j] = np.inner(self.__aij[i], self.__aij[j])
self.__gij[j][i] = self.__gij[i][j]
return self.__gij
[ドキュメント]
def SetLatticeParameters(self, latt, GuessLatticeSystem = 0, update_aij = True):
"""
格子定数 (a, b, c, α, β, γ) を設定し、空間群の格子系を更新し、計量テンソルなどを再計算します。
:param latt: list of float - 設定する格子定数 [a, b, c, alpha, beta, gamma]。
:param GuessLatticeSystem: int, optional - 格子定数から格子系を推定するかどうかを示すフラグ。デフォルトは0。
:param update_aij: bool, optional - 格子ベクトル `aij` も更新するかどうか。デフォルトはTrue。
:returns: None
"""
# print("latt=", latt)
self.__Latt = latt
SPG = self.SpaceGroup()
SPG.SetLatticeParameters(latt, GuessLatticeSystem)
SPG.SetLatticeAxis()
self.CalcMetrics(update_aij)
[ドキュメント]
def lattice_parameters(self, UseAtomicUnit = 0):
"""
格子定数 (a, b, c, α, β, γ) を取得します。
`UseAtomicUnit` が1の場合、原子単位系に変換して返します。
:param UseAtomicUnit: int, optional - 原子単位系を使用するかどうかを示すフラグ。
1の場合、原子単位系に変換します。デフォルトは0。
:returns: list of float - 格子定数のリスト。
"""
latt = list(self.__Latt)
if UseAtomicUnit:
# `a0` は原子単位の長さスケール因子として機能します。
return 1.0 / (a0 * 1.0e10) * latt
return latt
[ドキュメント]
def LatticeParameters(self, UseAtomicUnit = 0):
"""
格子定数 (a, b, c, α, β, γ) を取得します。(エイリアス)
`lattice_parameters` メソッドと同じです。
:param UseAtomicUnit: int, optional - 原子単位系を使用するかどうかを示すフラグ。
1の場合、原子単位系に変換します。デフォルトは0。
:returns: list of float - 格子定数のリスト。
"""
return self.lattice_parameters(UseAtomicUnit = UseAtomicUnit)
[ドキュメント]
def CalLatticeParametersFromMetrics(self, update_aij = False):
"""
計量テンソル `__gij` から格子定数 (a, b, c, α, β, γ) を計算し、設定します。
このメソッドは、格子ベクトル `__aij` も更新しますが、`update_aij` がFalseの場合、
元の `__aij` は復元されます。
警告: ベータ角の計算に `self.__gij[2][1]` が使用されていますが、
これは `self.__gij[2][0]` (g_ca) の意図であった可能性があります。
既存のロジックを変更せずに維持します。
:param update_aij: bool, optional - 格子ベクトル `aij` を更新し続けるかどうか。
Falseの場合、計算後に元の`aij`を復元します。デフォルトはFalse。
:returns: None
"""
# print("gij=", self.__gij)
a = RoundParameter(sqrt(self.__gij[0][0]), 5.0e-6)
b = RoundParameter(sqrt(self.__gij[1][1]), 5.0e-6)
c = RoundParameter(sqrt(self.__gij[2][2]), 5.0e-6)
alpha = self.RoundLatticeParameters(degacos(self.__gij[1][2] / b / c), 1.0e-6)
beta = self.RoundLatticeParameters(degacos(self.__gij[2][1] / c / a), 1.0e-6)
gamma = self.RoundLatticeParameters(degacos(self.__gij[0][1] / a / b), 1.0e-6)
self.SetLatticeParameters([a, b, c, alpha, beta, gamma], update_aij)
[ドキュメント]
def set_lattice_vectors(self, aij):
"""
格子ベクトル行列 `aij` を設定し、それに基づいて計量テンソルと格子定数を再計算します。
設定後、内部の格子ベクトル (`__aij`) は与えられた `aij` で上書きされます。
:param aij: numpy.ndarray or list of list of float - 設定する3x3の格子ベクトル行列。
:returns: None
"""
self.__aij = np.array(aij)
self.CalMetricsFromLatticeVector();
self.CalLatticeParametersFromMetrics(update_aij = False);
# CalLatticeParametersFromMetrics re-calculates aij through SetLatticeParameters()
# Recover original aij
self.__aij = np.array(aij)
[ドキュメント]
def SetLatticeVectors(self, aij):
"""
格子ベクトル行列 `aij` を設定し、それに基づいて計量テンソルと格子定数を再計算します。(エイリアス)
`set_lattice_vectors` メソッドと同じです。
:param aij: numpy.ndarray or list of list of float - 設定する3x3の格子ベクトル行列。
:returns: None
"""
return self.set_lattice_vectors(aij)
[ドキュメント]
def lattice_vectors(self, UseAtomicUnit = 0):
"""
格子ベクトル行列 `aij` を取得します。
`UseAtomicUnit` が1の場合、原子単位系に変換して返します。
:param UseAtomicUnit: int, optional - 原子単位系を使用するかどうかを示すフラグ。
1の場合、原子単位系に変換します。デフォルトは0。
:returns: numpy.ndarray - 3x3の格子ベクトル行列。
"""
if UseAtomicUnit:
# `a0` は原子単位の長さスケール因子として機能します。
return 1.0 / (a0 * 1.0e10) * self.__aij
else:
return self.__aij
[ドキュメント]
def LatticeVectors(self, UseAtomicUnit = 0):
"""
格子ベクトル行列 `aij` を取得します。(エイリアス)
`lattice_vectors` メソッドと同じです。
:param UseAtomicUnit: int, optional - 原子単位系を使用するかどうかを示すフラグ。
1の場合、原子単位系に変換します。デフォルトは0。
:returns: numpy.ndarray - 3x3の格子ベクトル行列。
"""
return self.lattice_vectors(UseAtomicUnit)
[ドキュメント]
def count_by_type(self, atype, mode = 'short', target = 'expanded'):
"""
指定された原子タイプ `atype` の原子サイト数をカウントします。
`target` パラメータで展開された原子サイトまたは非対称原子サイトのどちらを対象にするか、
`mode` パラメータで原子名の比較方法（短い名前か完全な名前か）を指定できます。
:param atype: str - カウントする原子タイプ名。
:param mode: str, optional - 原子名の比較モード。
'short': 原子タイプのみで比較（例: 'O'）。
'full': 完全な原子名で比較（例: 'O1'）。
デフォルトは'short'。
:param target: str, optional - カウントの対象となる原子サイトリスト。
'expanded': 展開された全原子サイト。
'asymmetric': 非対称単位内の原子サイト。
デフォルトは'expanded'。
:returns: int - 指定された原子タイプの原子サイト数。
"""
atype = atype.upper()
if target == 'expanded':
AtomSites = self.ExpandedAtomSiteList()
else:
AtomSites = self.AtomSiteList()
nsites = len(AtomSites)
n = 0
if mode == 'short':
for i in range(nsites):
t = AtomSites[i]
name = t.AtomNameOnly().upper()
if name == atype:
n += 1
else:
for i in range(nsites):
t = AtomSites[i]
name = t.AtomName().upper()
if name == atype:
n += 1
return n