分子を表示させる雑記　その１

#前書き

学会やらなんやらで忙しい時期を抜け出しましたー

で、いい加減Pythonで何かしらの成果物を発表しようと考えました(イマココ)

過去の私の雑記を見た人は分かるかもですが、専攻は計算化学です

主に反応設計をしています

分かりやすく言えば、化学系のシミュレーター屋さんと考えていただければ

```
（　＾ω＾）・・・(全然わかりやすくないんですがそれは)

 ⊂彡☆))Д´）
```

量子力学齧ったり有機化学齧ったりと色々齧り虫しているわけですよ

~~(お尻は齧って)ないです~~

そんなわけ(?)で、世の中には色々と便利な分子モデリングがあるんですが、

ぶっちゃけ私的には使いにくいんですよ

特に二面角の操作とか、どことどこを面にしてんねん！とか

仕様書を読めばいいだろいい加減にしろ！と思われるかもしれませんが、

正直読む気が失せる・・・英文字の束とか英語初級者に無理言うなって

使いやすいようにいじればいいだろいい加減にしろ！という人もいるとは思いますが、

権利ガｰと煩い昨今、波風もできる限り抑えたい

というより、ソースコード読むのはMOPAC7ですらヒーヒー言っているのに、

無理ゲーと言わざるを得ない・・・

そーでした、MOPACのソースコード解析は別途Qiitaに上げようと思うので、

気になる人は**<font color="FireBrick">地獄に付き合ってね！</font>**


#本題

流れとしてはz-matrixを読み込んで3Dとして表示させることがしたい

なので、まず第一段階として、結合距離、角度、二面角の計算プログラムを作ります

初っ端からz-matrix(内部座標系)はイマイチ作り方を把握できてないので、

今回はカーテシアン座標、直交座標系からの計算にしました

```prot.py
import math
A=input()
with open(A,encoding='utf-8') as B:
    txt = B.read()

C=txt.split()

del C[:(C.index('Required')+3)]
number = [i*4 for i in range(int(len(C)/4))]

D=[]
s="No"

def point(a,b):
    p=float(a)-float(b)
    return p

def line(c,d,e):
    b=float((c**2+d**2+e**2)**0.5)
    return b

def inner(f,g):
    inn=float(f*g)
    return inn

def arg(h):
    degr=math.degrees(math.acos(h))
    return degr

print(str(len(number)) + "個の原子があります")

print("いくつ選択しますか？")
n=int(input())

while s=="No":
    for i in range(n):
        print("選択する分子はどれですか？")
        D.append(input())
    print("よろしいですか?")
    s=input()

if n<=1:
    print("計算できませんでした")

elif n>=3:
    x1=point(C[number[int(D[0])-1]+1],C[number[int(D[1])-1]+1])
    y1=point(C[number[int(D[0])-1]+2],C[number[int(D[1])-1]+2])
    z1=point(C[number[int(D[0])-1]+3],C[number[int(D[1])-1]+3])
    x2=point(C[number[int(D[2])-1]+1],C[number[int(D[1])-1]+1])
    y2=point(C[number[int(D[2])-1]+2],C[number[int(D[1])-1]+2])
    z2=point(C[number[int(D[2])-1]+3],C[number[int(D[1])-1]+3])
    
    if n==4:
        x3=point(C[number[int(D[3])-1]+1],C[number[int(D[1])-1]+1])
        y3=point(C[number[int(D[3])-1]+2],C[number[int(D[1])-1]+2])
        z3=point(C[number[int(D[2])-1]+3],C[number[int(D[1])-1]+3])
        bx=float(y1*z2)-float(z1*y2)
        cx=float(y1*z3)-float(z1*y3)
        by=float(z1*x2)-float(x1*z2)
        cy=float(z1*x3)-float(x1*z3)
        bz=float(x1*y2)-float(y1*x2)
        cz=float(x1*y3)-float(y1*x3)
        b3=line(bx,by,bz)
        b4=line(cx,cy,cz)
        i4=inner(bx,cx)
        i5=inner(by,cy)
        i6=inner(bz,cz)
        cs2=(i4+i5+i6)/(b3*b4)
        print(cs2)
        larg2=arg(cs2)
        print(D[0]+"番目の"+C[number[int(D[0])-1]]+\
              D[1]+"番目の"+C[number[int(D[1])-1]]+\
              D[2]+"番目の"+C[number[int(D[2])-1]]+"で作られる平面と"+\
              D[0]+"番目の"+C[number[int(D[0])-1]]+\
              D[1]+"番目の"+C[number[int(D[1])-1]]+\
              D[3]+"番目の"+C[number[int(D[3])-1]]+"で作られる平面のなす角度は"+\
              str(larg2)+"゜ です。")
    else:
        b1=line(x1,y1,z1)
        b2=line(x2,y2,z2)
        i1=inner(x1,x2)
        i2=inner(y1,y2)
        i3=inner(z1,z2)
        cs1=(i1+i2+i3)/(b1*b2)
        larg1=arg(cs1)
        print(D[1]+"番目の"+C[number[int(D[0])-1]]+"を頂点とした角度は"+str(larg1)+"゜ です。")

else:
    x1=point(C[number[int(D[0])-1]+1],C[number[int(D[1])-1]+1])
    y1=point(C[number[int(D[0])-1]+2],C[number[int(D[1])-1]+2])
    z1=point(C[number[int(D[0])-1]+3],C[number[int(D[1])-1]+3])
    Bond=line(x1,y1,z1)
    print(D[0]+"番目の"+C[number[int(D[0])-1]]+"と"+\
    D[1]+"番目の"+C[number[int(D[0])-1]]+"結合距離は "+str(Bond)+" です。")
```
なんだこのあたまわるわるプログラムは・・・たまげたなぁ

あるプログラムで生成されるcomファイルを読み込んで計算して！というやつです

多次元配列を作ろうとしたんですが、うまく作れなかったのでBUZAMAなプログラムになりました


![草草野草.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/293790/6ce2005f-1498-7432-00aa-02fe04e201e2.png)

計算自体は何とかなっているんで許してクレメンス

計算をスマートにするために、分子、X座標、Y座標、Z座標をそれぞれリスト化

```listp.py
M = []; X = []; Y = []; Z = []

for i in range(int(len(C)/4)):
    i = i*4
    M.append(C[i])

for i in range(int(len(C)/4)):
    i = i*4+1
    X.append(C[i])

for i in range(int(len(C)/4)):
    i = i*4+2
    Y.append(C[i])

for i in range(int(len(C)/4)):
    i = i*4+3
    Z.append(C[i])
```
としたんですど、多分もっとスマートにできるのでは？

で、関数を利用することに

```listpro2.py
Number=int(len(C)/4)

def abst(L,u,n):
    for i in range(n):
        i = i*4+u
        L.append(C[i])

M = []; X = []; Y = []; Z = []

abst(M,0,Number); abst(X,1,Number); abst(Y,2,Number); abst(Z,3,Number)
```
これで6行かな？短縮

これ以上短くは思いつかないので、思いついた人は教えてクレメンス

#後書き

次回は数字の書かれたtxtファイルを読み込んで、

自動的に距離、角度、二面角の計算をできるようにしたいと思います

<br>
P.S. 今MOPACのソースコードの解析を行っているんですが、

皆さんはソースコードの解析ってどんな感じでやっているんですかね？

頭からぼりぼりやろうとは決めたんですけど、

手書き癖と言うかそういうのがあって、コメントでワラワラ書くのに慣れないんですよね・・・

ほんとは全部印刷してガタガタ書きたいんですが、印刷すると枚数が多すぎるんですよ

今の所、OneNoteを利用してペンタブ片手にPDF化したソースコードへがりがり書いているんです

こんな感じでやってるよーとかあったら教えてください

![おねがいsm.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/293790/3512dc82-6395-c1e7-8476-a1af54d7d0a4.png)

<br><br><br>

#個人的な感想という名の愚痴
勘の良い人は気付いたと思います、Numpy使えよ！と

正直Numpyを使った方が数値の扱いは楽です

まぁ今回は基本的にPythonのみインストールした環境で作ると縛っています

つまり、標準モジュールで戦うってことです

一応例外として、自作のモジュールは可と言う感じにしてます

とはいえ、モジュールなんて作れないんですけど・・・

標準モジュールのみでFortranとPythonを繋げられたらなんですけどネ

Pythonの強みであるライブラリが豊富って言うのは私にとっては旨味じゃないんですよね

これをやったらこの結果が出たラッキー、で終わらないわけですよ

だってソコ突っ込まれたら終わりなわけじゃないですか

楽できるに越したことないとは思いますが、少なくとも私は自身を化学系だと考えているので、

何故そうなるのか、どうしてその結果が得られるのかを把握しなければならないと思っています

だから、他人に説明できもしない物を基本的には使いたくないわけですよ

import math ですら私的には屈辱でしたね使うの、だって中の動きを説明できないので

ただ今回はarccos、cosの値の計算をさせているんだ、というふわふわ理論で納得することにしました

角度の範囲が固定されている以上、cosの値の範囲も決まっているので、

単純な比例で角度を出すこともできるはず、とは思いますが

それに拘り過ぎても時間の無駄でしかないので、標準モジュールの使用は許可したという感じです

本当にPython用にFortranでモジュール作ってお手軽に流用できないものか・・・

最悪exeファイルの実行できれば、多分モジュール化しなくても済みそうではある

うごごごご、標準モジュールにそんなのないかなぁ