とりあえず書く
ベンゼンの重心計算のプログラムを作成したので載せてみる
例えばこのようなデータが16128個連なっている場合、その重心を求めたいとする。ちなみに左から4つ目から6つ目がxyz座標のデータであり、最後は速度xyzデータとなっている。
BENZENE
16128
1BNZ C1 1 5.769 2.782 5.588 0.3603 0.2619 0.5938
1BNZ C2 2 5.718 2.819 5.710 -0.1587 -0.4771 0.6026
1BNZ C3 3 5.596 2.768 5.747 -0.2580 -0.3057 0.5127
1BNZ C4 4 5.521 2.693 5.659 -0.0603 -0.0666 0.1392
1BNZ C5 5 5.570 2.666 5.533 0.5944 0.2517 0.3192
1BNZ C6 6 5.696 2.708 5.498 0.3389 0.8589 0.1068
1BNZ H1 7 5.866 2.817 5.556 0.8850 -2.4481 -1.0496
1BNZ H2 8 5.773 2.880 5.780 0.9081 -3.0610 2.1750
1BNZ H3 9 5.560 2.794 5.845 -2.0099 0.0088 -0.1769
1BNZ H4 10 5.424 2.665 5.697 0.1623 0.4350 1.1151
1BNZ H5 11 5.512 2.607 5.464 -0.9770 -1.7258 3.0982
1BNZ H6 12 5.735 2.667 5.406 -0.1421 0.1094 0.2260
2BNZ C1 13 3.809 3.388 2.568 0.1936 0.6084 0.3144
2BNZ C2 14 3.710 3.320 2.635 0.4378 0.1620 0.2311
2BNZ C3 15 3.583 3.373 2.634 0.4549 0.2041 -0.2045
2BNZ C4 16 3.552 3.479 2.552 -0.7824 -0.5633 -0.7582
2BNZ C5 17 3.652 3.543 2.483 -0.3636 -0.3166 0.0660
2BNZ C6 18 3.783 3.504 2.499 -0.1683 0.2085 -0.2319
2BNZ H1 19 3.907 3.343 2.567 -0.5937 -1.3216 2.3472
2BNZ H2 20 3.738 3.231 2.689 0.5879 0.8724 1.3627
2BNZ H3 21 3.507 3.318 2.687 0.7589 -2.4117 -2.2964
2BNZ H4 22 3.451 3.518 2.551 -0.8241 -0.4911 -3.6901
2BNZ H5 23 3.633 3.635 2.430 -1.7266 -0.5897 0.0502
2BNZ H6 24 3.868 3.547 2.448 -0.8328 2.0720 0.1500
3BNZ C1 25 3.770 0.062 3.474 0.1817 0.0610 0.4697
3BNZ C2 26 3.796 -0.022 3.580 0.3213 0.1508 0.5058
3BNZ C3 27 3.704 -0.035 3.681 -0.0229 -0.2349 0.1411
3BNZ C4 28 3.590 0.041 3.678 0.5097 0.5785 0.4873
3BNZ C5 29 3.563 0.124 3.571 0.1231 0.0653 0.1832
3BNZ C6 30 3.651 0.131 3.466 0.0520 -0.1970 0.1035
3BNZ H1 31 3.851 0.084 3.407 -0.0536 0.3014 0.2611
3BNZ H2 32 3.897 -0.060 3.587 -0.4650 -2.7843 -1.8647
3BNZ H3 33 3.723 -0.098 3.766 0.8231 -1.1911 -0.7223
3BNZ H4 34 3.515 0.020 3.754 -0.8979 -1.0196 -1.2396
3BNZ H5 35 3.465 0.168 3.559 0.2341 -0.0289 -1.2420
3BNZ H6 36 3.638 0.192 3.378 1.7645 -0.0052 -0.0504
# coding: utf-8
"""----------------------------------------------------
Calculate the center of benzens
----------------------------------------------------"""
def cal_center(scalar_list):
# リストの重心計算
# return sum(scalar_list)/len(scalar_list)
return sum(scalar_list[0:6])/6 #炭素だけで計算だから6個だけ
with open('./melting.325K.gro', 'r') as f:
lines = f.readlines()
# print(lines)
# 必要な行抽出, 整形
# element_line_list = lines[2:]
# print(element_line_list)
# element_list = [element_line.split() for element_line in element_line_list]
# # なんか変な行があるので除外する 要素が9個ある行だけ読む
element_list = []
for element_line in lines:
if len(element_line) == 45:
element_list.append(element_line.split())
# print(element_list)
# print(element_list)
e_list = []
for e in element_list:
if len(e) == 6:
e_list.append(e)
elif len(e) == 5:
if e[0] != "BENZENE":
# スライス処理をかく
# ['1BNZ', 'C1', '1', '0.099', '-0.045', '-0.034']
e.insert(1,e[1][0:2]) #2個目の要素を2個スライスしたやつをinsert
e.insert(2,e[2][2:]) #
del e[3]
e_list.append(e)
# print(e)
# print(e_list)
# print(e_list[1])
# print(e_list)
# 重心を格納するリスト
center_list = []
bensen_size = 12
# 12行ごとにイテレート この場合はベンゼン
# for st_idx in range(0, len(lines), bensen_size): #３個目の引数は12行ごと
# bensen = lines[st_idx: st_idx+bensen_size]
# print(bensen)
# for bensen in e_list(0, len(lines), bensen_size):
x_list = []
y_list = []
z_list = []
# for element in range(0, len(lines), bensen_size):
for e in e_list:
# x, y, z = [float(i) for i in e_list[3:6]]# 座標を抽出, floatに変換
# print(element)
# print(e_list[3])
# print(e_list[4])
# print(e_list[5])
# print(e_list[element][3])
x = float(e[3])
y = float(e[4])
z = float(e[5])
# x = float(e_list[element][3])
# y = float(e_list[element][4])
# z = float(e_list[element][5])
# y = float(e_list[4])
# z = float(e_list[5])
x_list.append(x)
y_list.append(y)
z_list.append(z)
# print(x_list)
center = [cal_center(x_list), cal_center(y_list), cal_center(z_list)] #座標表示っぽくリスト作成
center_list.append(center) #リストの中にリスト
# 出力を文字列として変換
output_text = ""
output_text += 'frame\t t=0.000\n'
output_text += f'\t{len(center_list)}\n'
# print(center_list[0])
# >>> x,y,z
for i in range(len(center_list)):
# 重心
center = center_list[i]
center_text = "\t".join(['{:.4f}'.format(c) for c in center]) #重心の出力 一気に結合
bensen_index = i+1
# 出力
output_text += f'{bensen_index}BEZ\tBEZ\t{bensen_index}\t{center_text}\n'
# print(output_text[:200])
# テキストファイルへ保存
with open('./tmp_melting.325K.txt', 'w') as out:
out.write(output_text)
#close()