ベンゼンのxyz座標から重心座標を求める

とりあえず書く

ベンゼンの重心計算のプログラムを作成したので載せてみる

例えばこのようなデータが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
```

```python
# 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()

```