pygamessがTD-DFTとNMR計算をサポートしました

なんとなくこの本の15章をやりたくなったので、朝からゴニョゴニョしてた。色々できたのでバージョン上げておいた。

TD-DFT

データはおなじみのPubchemQCからダウンロードしてきてopenbabelでmolファイルにコンバートしたものを利用。

>>> from pygamess import Gamess
>>> from rdkit import Chem
>>> m = Chem.MolFromMolFile("examples/methyl_yellow.mol", removeHs=False)
>>> g = Gamess()
>>> g.dft_type("b3lyp", tddft=True)
>>> g.basis_sets("6-31G*")
>>> r = g.run(m)
>>> r.uv_spectra # (exitation ev, oscillator strength)

NMR spectra

これは時間がかかった(2020 MBAで30分以上)

>>> from pygamess import Gamess
>>> from rdkit import Chem
>>> m = Chem.MolFromMolFile("examples/C=CCBr.mol", removeHs=False)
>>> g = Gamess(num_cores=1) # PARALLEL EXECUTION IS NOT ENABLED.
>>> g.basis_sets("6-31G*") 
>>> g.run_type("nmr")
>>> r = g.run(m)
>>> r.isotropic_shielding

独習量子化学計算について

わかりやすくて良い本だった。ケミストだったら本書に沿って一通り手を動かしておけば、後々役に立つのではないかと思った。

こっちはいわゆるクックブックみたいな感じ

Free Ligand 1D NMR Conformational Signatures To Enhance Structure Based Drug Design

This was a very interesting paper.

The correlation of kon with free ligand preorganization implies 1D NMR signatures can augment structure−kinetic relationships; when kon is optimized by highly populated bioactive states, potency will be driven by koff.

Storing GAMESS calculation results in RDKit Chem object

Japan is under a state of emergency for the third time, which forced me to do self-quarantine. After all, my five days of GW holidays were spent developing the pygamess library.

This is an example of storing GAMESS calculation results into RDKit's Chem object. I will calculate the two molecules, ethane, and ethanol, to export as an SDF file later.

>>> from pygamess import Gamess
>>> from pygamess.utils import rdkit_optimize
>>> from rdkit import Chem
>>> ethane = rdkit_optimize("CC")
>>> ethanol = rdkit_optimize("CCO")
>>> g = Gamess()
>>> g.run_type("optimize")
>>> ethane_result = g.run(ethane)
>>> optimized_ethane = ethane_result.mol
>>> ethanol_result = g.run(ethanol)
>>> optimized_ethanol = ethanol_result.mol

Orbital energies and dipole moments including HOMO and LUMO energies are stored in the object. Each atom is also assigned a mulliken charge/lowdin charge.

>>> optimized_ethane.GetProp("HOMO")
'-0.46029999999999999'
>>> optimized_ethane.GetProp("LUMO")
'0.65559999999999996'
>>> optimized_ethane.GetProp("dipole_moment")
'1.9999999999999999e-06'
>>> optimized_ethane.GetProp("orbital_energies")
'-11.0355 -11.0352 -0.981 -0.8121 -0.572 -0.572 -0.4717 -0.4603 -0.4603 0.6556 0.6556 0.687 0.737 0.7814 0.7893 0.7893'
>>> for a in optimized_ethane.GetAtoms():
...   print("{}:\t{:.4f}".format(a.GetSymbol(), float(a.GetProp("mulliken_charge"))))
... 
C:  -0.1748
C:  -0.1748
H:  0.0583
H:  0.0583
H:  0.0583
H:  0.0583
H:  0.0583
H:  0.0583

Finally, I'll export the object to SDF.

>>> optimized_mols = [optimized_ethane, optimized_ethanol]
>>> w = Chem.SDWriter('test.sdf')
>>> for m in optimized_mols:
...   w.write(m)
... 
>>> w.close()

Looking at the exported SDF, we can see these results are persisted.

$ cat test.sdf
     RDKit          3D

  8  7  0  0  0  0  0  0  0  0999 V2000
   -0.7687    0.0082   -0.0152 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7687   -0.0082    0.0152 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1742   -0.1016    0.9863 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1352    0.9432   -0.4285 H   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1499   -0.8048   -0.6260 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.1499    0.8048    0.6260 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.1742    0.1016   -0.9863 H   0  0  0  0  0  0  0  0  0  0  0  0
    1.1352   -0.9432    0.4285 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0
  1  3  1  0
  1  4  1  0
  1  5  1  0
  2  6  1  0
  2  7  1  0
  2  8  1  0
M  END
>  <total_energy>  (1) 
-78.306179626800002

>  <HOMO>  (1) 
-0.46029999999999999

>  <LUMO>  (1) 
0.65559999999999996

>  <nHOMO>  (1) 
-0.46029999999999999

>  <nLUMO>  (1) 
0.65559999999999996

>  <dipole_moment>  (1) 
1.9999999999999999e-06

>  <dx>  (1) 
-0

>  <dy>  (1) 
-1.9999999999999999e-06

>  <dz>  (1) 
0

>  <orbital_energies>  (1) 
-11.0355 -11.0352 -0.981 -0.8121 -0.572 -0.572 -0.4717 -0.4603 -0.4603 0.6556 0.6556 0.687 0.737 0.7814 0.7893 0.7893

>  <atom.dprop.mulliken_charge>  (1) 
-0.174764 -0.174764 0.058255000000000001 0.058255000000000001 0.058254 0.058255000000000001 0.058255000000000001 0.058255000000000001

TODO: Add GAMESS calculation conditions as attributes, like "HF/6-31G**".

フロンティア軌道で化学を考える

私の一番大きな興味は蛋白質と低分子化合物の相互作用を如何に深く理解するかってことで、最近は量子化学のパラメータを記述子として使ったりしてるんだけど、ガチ量子系の研究者の言ってることがにわかの私には理解できないことがたまにある。

最近であったconfusionは励起状態への遷移に伴うインデュースドフィットの描像が古典力学の解釈でのそれとはなんか違うような気がして時間的な変化どうなってんのかねーという疑問が湧いたので、この本を読み直している。

レチナールとか視覚系は量子化学的な取り扱いしないといけないとは思っているけど、普通の低分子化合物の場合はどこまで非占有有軌道が効いてくるのか興味があるところ。

Cation-π相互作用をSAPT0してみた。

別件で調べ物をしていたら10年くらい前に検討したエントリを見つけたので、これをpsi4のSAPTでやったらどういう結果になるかな?と。

今回psi4コマンドで流すのでinputはこんな感じです(いずれはpsikitでサクッと計算できるようにしたい)。EmacsやめちゃったからVimでgamessのインプットいじってつくった。繰り返しの技を覚えたけど、こういうの手になじまないとめんどくさいですね。

molecule cation_pi {
     1 1
     N 26.56861  22.88472  72.47206
     C 25.05591  23.05763  72.38527
     C 27.18276  23.04159  71.08205
     C 27.14794  23.94295  73.40648
     C 26.89124  21.49601  73.01604
     H 24.65761  22.29324  71.71973
     H 24.84300  24.05067  71.99229
     H 26.95089  24.04166  70.71972
     H 26.73892  22.29197  70.42917
     H 28.26387  22.89551  71.15590
     H 26.89847  24.92550  73.00913
     H 28.22884  23.81270  73.44995
     H 26.70726  23.81001  74.39332
     H 26.46855  20.75492  72.33939
     H 26.44825  21.40230  74.00640
     H 27.97398  21.38820  73.06968
     H 24.63622  22.94695  73.38408
     --
     0 1
     C 30.61324  20.49644  72.40614 
     C 30.47183  21.03406  71.19420 
     C 30.79531  21.59051  73.36029 
     C 30.75216  22.77934  72.62888 
     C 31.00908  21.63948  74.74382 
     N 30.43745  22.45766  71.27884 
     C 30.92671  24.02761  73.23103 
     C 31.18435  22.86580  75.34238 
     C 31.14962  24.05129  74.58925 
     H 30.33993  20.55500  70.23268
     H 31.04481  20.72825  75.32842
     H 30.86231  23.00634  70.52231
     H 30.90815  24.93911  72.64846
     H 31.36065  22.92555  76.40870
     H 31.30378  24.99933  75.08996
     H 30.61446  19.44515  72.64676

     units angstrom
     no_reorient
     symmetry c1
}

set basis 6-31g_d_p_

energy('sapt0')

10年も経つと6-31G**でもサクッと終わるのね。subaracy!

結果はこのあたりを見ればいいんですかね。ESも効いてるけどInduction(CT+mixみたいな項?)とDispersionも効いてるから、次はカチオンの有無でトリプトファンのそれぞれの原子のRESP chargeがどのくらい変化するかとか双極子どうなっているかとか確認すればいいんかな。

Special recipe for scaled SAPT0 (see Manual):
  Electrostatics sSAPT0         -13.53632077 [kcal/mol]     
  Exchange sSAPT0                10.77916633 [kcal/mol]      
  Induction sSAPT0               -5.12922298 [kcal/mol]     
  Dispersion sSAPT0              -5.69289464 [kcal/mol]     
Total sSAPT0                    -13.57927206 [kcal/mol]     
-------------------------------------------------------

Scanning the torsional potential in Psikit(RDKit+Psi4)

Considering the conformational effects of the compound is important in Structure Based Drug Design, this paper discussed about it, in terms Protein-Ligand binding using torsional scan of each ligands(PDB:2JH0,2JH5,2JH6). They calculated torsional energy, and explained the relation between inhibitory activity and torsional energy.

Torsional scanning is the task of the quantum chemistry rather than that of chemoinformatics. But I wanted to conduct quantum chemical calculation as an extention of chemoinformatics way, so I implemented it in Psikit

日本語訳

jupyterでRDKitからのB3LYP/6-31G*でのDFT計算さっくりうごくの素晴らしい。

Psi4 install error

I got a this error. Conda always irritates me ;-)

>>> import psi4
Traceback (most recent call last):
  File "/Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/__init__.py", line 55, in <module>
    from . import core
ImportError: dlopen(/Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so, 2): Library not loaded: @rpath/libiomp5.dylib
  Referenced from: /Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so
  Reason: image not found

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/__init__.py", line 60, in <module>
    raise ImportError("{0}".format(err))
ImportError: dlopen(/Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so, 2): Library not loaded: @rpath/libiomp5.dylib
  Referenced from: /Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so
  Reason: image not found

@rpathという変数がMacでは何らかの理由により正しく設定されないらしい

朝からイラツキすぎてまともに解決する気力もないので、当面放置するかも。これから買い物いかないといけないし。 それにしてもCondaよくわからんエラー多すぎるわ。RDKitもそうだったしさー。あんまよくないんじゃないの?職場では使ってないしー

forum見たけど報告されてないっぽいし、あとで気力が戻ったら投げておきます(未定)。

追記 0800くらい

エントリ書いて朝ごはん食べたら冷製になったのでなおした。

$ install_name_tool -change @rpath/libiomp5.dylib /Users/kzfm/anaconda3/pkgs/intel-openmp-2018.0.0-h8158457_8/lib/libiomp5.dylib /Users/kzfm/anaconda3/lib/python3.6/site-packages/psi4/core.so

Mishima.syk #6 やりました

参加した皆様お疲れ様でした。色々と面白い話が聞けて楽しかったです。

僕は分子設計を軽くディスってきましたw 本当にもっと頑張って欲しいですね。リードオプティマイゼーションのスタンスからFMOを使っているんだけど、もっとアーリーな段階での使い方の議論は新鮮だったので今度試してみたい。

あとは告知の仕方に問題があったので改善したい。事前に知っていれば参加したのにと何人かのヒトに言われたので。twitterやれとしか言えないが、それも難しいでしょうからメールを送るとか考えないとあかんのかな…

懇親会で興味を伺ったらディプラーニングに注目しているヒトが多いようなので次回はそれ系のネタを集めようかなと思っています。

反省会は「鈴木屋」か「やごみ」でやりましょうw

pygamessの便利なところ

懇親会でpygamessじゃなくてもFacioで良かろうと言われたのだけど、対話環境よりも、スクリプトでやれたほうがいい場面が多いので書いておく。ちなみに僕はAvogadroユーザーです。

適当なsdfが手元になかったので

C
CC
CCC

というSMILESのリストを

babel -ismi carbons.smi -osdf carbons.sdf

とやってsdfにしたものを使いました。

sdfの分子群に対して計算したい場合、対話環境だとチマチマと作業しないといけないが、スクリプトだとforループを回せばよいし、結果がmoleculeオブジェクトで返ってくるので、そのあとのケモインフォマティクス的な取り回しが楽です。Gamessのアウトプット睨みながらコピペとか苦行だし、アウトプットの処理用のパーサー書くならinputから全てスクリプトで処理したほうがいいんじゃないかなと。

>>> import pybel
>>> from pygamess import Gamess
>>> g = Gamess()
>>> mols = pybel.readfile('sdf', 'carbons.sdf')
>>> for mol in mols:
...   nmol = g.run(mol)
...   print nmol, nmol.energy
... 
C   
 -39.7265813363
CC  
 -78.3052918313
CCC 
 -116.885136991

ケミストのように化合物を個別のものとして一つ一つ丁寧に対応していくのか、情報論的にまとめて取り扱いたいのかとかそういうあたりで使い分ければいいのだと思う。オービタル見たい場合にはGUI必須だけど、エネルギーみたいな数字だけに関心があるのであればGUIは要らないかな。