本稿は計算化学等で周期境界条件を用いた計算を行っている方向けのメモです。
並進ベクトル（Translation Vectors: TV）を含む3次元座標からcif形式のファイルを生成する方法を述べます。
以下のような直交座標系（Cartesian coordinate system）で記述された構造をcif形式に変換するには分率座標系（fractional coordinate system）に変換する必要があります。
Si 1.367182000000 4.101546000000 1.367182000000
Si 0.000000000000 0.000000000000 2.734364000000
Si 1.367182000000 1.367182000000 4.101546000000
Si 0.000000000000 2.734364000000 0.000000000000
Si 4.101546000000 4.101546000000 4.101546000000
Si 2.734364000000 0.000000000000 0.000000000000
Si 4.101546000000 1.367182000000 1.367182000000
Si 2.734364000000 2.734364000000 2.734364000000
TV 5.468730000000 0.000000000000 0.000000000000
TV 0.000000000000 5.468730000000 0.000000000000
TV 0.000000000000 0.000000000000 5.468730000000
直交座標系の点$(x,y,z)$を分率座標系の点$(x_{\text{frac}},y_{\text{frac}},z_{\text{frac}})$に変換する公式が知られており、
\begin{equation}
v=\sqrt{1-\cos ^2 \alpha-\cos ^2 \beta-\cos ^2 \gamma+2 \cos \alpha \cos \beta \cos \gamma} \quad
\end{equation}
として
\left[{\begin{array}{l}
x_{\text{frac}} \\
y_{\text{frac}} \\
z_{\text{frac}}
\end{array}}\right]=\left[\begin{array}{ccc}
1 / a & -\cos \gamma / a \sin \gamma & (\cos \alpha \cos \gamma-\cos \beta) / a v \sin \gamma \\
0 & 1 / b \sin \gamma & (\cos \beta \cos \gamma-\cos \alpha) / b v \sin \gamma \\
0 & 0 & \sin \gamma / c v
\end{array}\right]\left[\begin{array}{l}
x \\
y \\
z
\end{array}\right]
とすればよいことが知られています。ここで、$a$、$b$、$c$は各並進ベクトルの長さ、$\alpha$、$\beta$、$\gamma$は各並進ベクトルの夾角です。
これにしたがってxyz形式からcif形式に変換するコードの例を示します。
xyz2cif.py
import math
import os
print("Convert xyz cordinates into cif format. (by xyz2cif.py v1.0)")
filepath = 'C:/Users/'+os.getlogin()+'/Downloads/'
filename_list = ["Si.xyz", "CuO.xyz", "ZnO.xyz"]
for i in range(len(filename_list)):
target_file = filepath + filename_list[i]
with open(target_file, encoding='utf-8') as f:
tmp_atomnum = 0; tmp_nodenum = 0; atomcount = 0
atomnamelist=[]
coordlist = []
TVlist = []
TVcount = 0
while line := f.readline():
elem = line.split()
if len(elem) >= 4:
if elem[0] == "TV":
TVlist.append([float(elem[1]), float(elem[2]), float(elem[3])])
TVcount += 1
if TVcount == 4:
print("Warning: only three TVs are allowed.")
else:
atomnamelist.append(str(elem[0]))
coordlist.append([float(elem[1]), float(elem[2]), float(elem[3])])
if not line:
break
# angström
a_length = TVlist[0][0]
b_length = math.sqrt(TVlist[1][0]**2.0 + TVlist[1][1]**2.0)
c_length = math.sqrt(TVlist[2][0]**2.0 + TVlist[2][1]**2.0 + TVlist[2][2]**2.0)
# radian
alpha = math.acos(sum([i*j for (i, j) in zip(TVlist[1], TVlist[2])])/(b_length * c_length)) # b-c
beta = math.acos(sum([i*j for (i, j) in zip(TVlist[2], TVlist[0])])/(c_length * a_length)) # c-a
gamma = math.acos(sum([i*j for (i, j) in zip(TVlist[0], TVlist[1])])/(a_length * b_length)) # a-b
# degree
degree_alpha = math.acos(sum([i*j for (i, j) in zip(TVlist[1], TVlist[2])])/(b_length * c_length)) * (180/math.pi) # b-c
degree_beta = math.acos(sum([i*j for (i, j) in zip(TVlist[2], TVlist[0])])/(c_length * a_length)) * (180/math.pi) # c-a
degree_gamma = math.acos(sum([i*j for (i, j) in zip(TVlist[0], TVlist[1])])/(a_length * b_length)) * (180/math.pi) # a-b
print("\n\nTV coordinates to cell parameters < " + filename_list[i]) # TVs as a lower triangular matrix
print('{:>11.8f}'.format(a_length), '{:>11.8f}'.format(b_length), '{:>11.8f}'.format(c_length), sep=", ")
print('{:>11.8f}'.format(degree_alpha), '{:>11.8f}'.format(degree_beta), '{:>11.8f}'.format(degree_gamma), sep=", ")
print("Surface area (ang^2) =", '{:>14.10f}'.format(TVlist[0][0] * TVlist[1][1]))
print("Cell volume (ang^3) =", '{:>14.10f}'.format(TVlist[0][0] * TVlist[1][1] * TVlist[2][2]))
# create a cif file
outputall_file = filepath + filename_list[i] + ".cif"
with open(outputall_file, mode='w', encoding='utf-8') as outf:
outf.write("# generated using xyz2cif.py (by HN)\n" + filename_list[i] + "\n")
outf.write("_symmetry_space_group_name_H-M \'P 1\'\n_chemical_name_common ?\n_chemical_formula_sum ?\n")
outf.write("_cell_length_a "+'{:>11.8f}'.format(a_length) + "\n")
outf.write("_cell_length_b "+'{:>11.8f}'.format(b_length) + "\n")
outf.write("_cell_length_c "+'{:>11.8f}'.format(c_length) + "\n")
outf.write("_cell_angle_alpha "+'{:>11.8f}'.format(degree_alpha) + "\n")
outf.write("_cell_angle_beta "+'{:>11.8f}'.format(degree_beta) + "\n")
outf.write("_cell_angle_gamma "+'{:>11.8f}'.format(degree_gamma) + "\n")
outf.write("_symmetry_Int_Tables_number 1\n_cell_volume " + '{:>11.8f}'.format(TVlist[0][0] * TVlist[1][1] * TVlist[2][2]) + "\n_cell_formula_units_Z " + str(len(coordlist)) + "\n")
outf.write("loop_\n _symmetry_equiv_pos_site_id\n _symmetry_equiv_pos_as_xyz\n 1 'x, y, z'\n")
outf.write("loop_\n _atom_site_type_symbol\n _atom_site_label\n _atom_site_symmetry_multiplicity\n _atom_site_fract_x\n _atom_site_fract_y\n _atom_site_fract_z\n _atom_site_occupancy\n")
for j in range(len(atomnamelist)):
v = math.sqrt(1 - math.cos(alpha)**2 - math.cos(beta)**2 - math.cos(gamma)**2 + 2*math.cos(alpha)*math.cos(beta)*math.cos(gamma))
a_frac = coordlist[j][0] - coordlist[j][1] / math.tan(gamma) + coordlist[j][2] * ((math.cos(alpha)*math.cos(gamma)-math.cos(beta))/(v*math.sin(gamma)))
b_frac = coordlist[j][1] / math.sin(gamma) + coordlist[j][2] * ((math.cos(beta)*math.cos(gamma)-math.cos(alpha))/(v*math.sin(gamma)))
c_frac = coordlist[j][2] * math.sin(gamma) / v
outf.write(' {:>2s}'.format(atomnamelist[j]))
outf.write(' {:>2s}'.format(atomnamelist[j]) + str(j) + " 1")
outf.write(' {:>9.8f}'.format(a_frac / a_length)) # x componet only
outf.write(' {:>9.8f}'.format(b_frac / b_length)) # y componet
outf.write(' {:>9.8f}'.format(c_frac / c_length)) # z componet
outf.write(" 1\n")