GAMESS を使って RESP 電荷を計算する


本記事はこちら（https://qiita.com/Ag_smith/items/430e9efb32a855d4c511）を参考にGAMESS を利用して分子動力学計算用の RESP 電荷を計算する方法についてのメモです。


想定している環境

今どきでそこそこ良いスペックの Mac（Catalina）

PyMOL（ファイルを眺めたり変換するだけなのでオープンソース版でも可）


Molby のインストール

Molby は分子構造を作成するためのソフト。Gaussian でいうところの GaussView 的に使える。GAMESS の入力ファイルをわかりやすく作ってくれる。Molby の公式サイト（http://molby.osdn.jp/doc/ja/index.html）からファイルをダウンロードしてインストールするだけ。

GAMESS のインストール

GAMESS (US) の公式サイト（https://www.msg.chem.iastate.edu/gamess/download.html）からメールアドレス等を登録すると、ユーザ名とパスワードが送られてくる。今回ほしいのは「Pre-compiled Binary Distributions」の中にある「GAMESS version September 30, 2019 R2 for Apple MacOS X」。ダウンロードしたらアプリケーションフォルダに展開しておく。 

ワークステーションやスパコンなどで計算する場合は適切なバージョンをインストールしてください。


実際の計算方法

計算したい分子の PDB ファイルを用意する。

Molby で PDB ファイルを読み込む。

必要であれば水素を付加する。

「MD/MM」→「Guess MM/MD Parameters」
「Guess atom types」にチェックを入れた状態で実行すると Antechamber がatom type を割り当ててくれる。

「MD/MM」→「GAMESS and RESP」
Step1の「Create GAMESS Input」を押すと設定画面が現れる。

設定ファイルとしては、以下の２項目が追加されている。
 $ELPOT  IEPOT=1 OUTPUT=PUNCH WHERE=PDC $END
 $PDC    CONSTR=NONE PTSEL=CONNOLLY $END

ローカルマシンで計算する場合

「Execute GAMESS on this machine」にチェックを付けて、実行ファイルのパスを指定する。（例：/Applications/gamess_2018R1/gamess.30Sep2018R3.x）
「OK」を押すと計算が始まる。それなりに負荷がかかるので注意。

別のマシンで計算する場合

「Edit GAMESS Input and Go」から現在の設定ファイル（.inp）が見れるので、コピーして ligand.inp などとして保存する。


RESP電荷計算後の処理

実行後はligand.datが作成されるので、Step2 の「Import GAMESS dat」から読み込む。Step3 の「Run RESP...」を押すと RESP 電荷が適用される。

この状態で、ligand_resp.pdb および ligand.psf として保存する。ちなみに電荷が含まれているのは .psf の方。PyMol で ligand_resp.pdb を読み込み、「Export molecule」から ligand_resp.mol2 として保存する。この時、「Original atom order」にチェックを入れておく。

ligand_resp.mol2、ligand_resp.psf を強いテキストエディタ（例：mi）で開く。矩形選択モードに変更して ligand_resp.psf の電荷をコピーする。Excel に貼り付けて  ROUND(元の電荷, 5) で桁を丸め込む。丸め込んだものについて総和を計算する。ゼロ（あるいは設定した電荷）にならず0.0001 など端数が生じることがあるので、適切に電荷を調節して総電荷が整数になるようにする。計算結果を ligand_resp.mol2 に貼り付ける（一番右のカラム）。atom type もligand_resp.psf から ligand_resp.mol2 にコピーする（座標の次のカラム）。

ligand_resp.mol2 に二重結合を反映させる。@<TRIPOS>BOND 行以下に結合の情報が書かれている（例：1 3 5 1）。例では結合番号１は原子３と５の間にあり、１重結合である、という形式です。この際、Molby の原子番号はゼロ始まりなのに対して、mol2 では１始まりなので注意してください。構造を見ながら、二重結合である部位については一番右の数字を２に変更してください。ちなみにアミドの場合は am、芳香環やCOO-については ar を指定します。

antechamber を使ってファイル形式を変換すると、AMBER（tleap） で読み込めるパラメータファイルができます。（http://ambermd.org/antechamber/ac.html#antechamber）

antechamber -i ./ligand_resp.mol2 -fi mol2 -o ligand.prep -fo prepi -nc 0


コメント

コメントを投稿

このブログの人気の投稿

Excelで近似直線の傾きと切片を得る

イメージ
近似直線・近似曲線を求めるには SLOPE 関数と INTERCEPT 関数を使おう。 学生実習のTAをしていて、気になった話。 １年生の実習ということで、酵素活性を測定してミカエリス・メンテンの  K m 値を求めるというものでした。目的の１つとして、Excelを使ってセルに計算式を入れて必要な値を求める方法を押えるというのがあります。その時、グラフの機能から近似直線の傾きと切片を得るのですが、どうしても傾きと y 切片は手入力で別のセルに書き込む必要がありました。これはとても不便で、何のための表計算ソフトなんだと思ったので少し調べてみました。 グラフから作成した近似式の係数を取得する汎用性ある関数はないようです。幸いなことに、直線近似については関数を使って別途求めることで係数を計算できます[1]。 SLOPE関数: SLOPE(yの値, xの値) INTERCEPT関数: INTERCEPT(yの値, xの値) SLOPE関数では近似直線の傾きが、INTERCEPT関数ではy切片が得られます。これを使うと表示された関数から手入力を介すること無く計算を進められます。 他の曲線についてはどうするんだということが気になりますが、直線にしてやれば良いようです。 累乗近似: y = ax^b 自然対数を求めるLN関数を使うと、LN(y) = LN(a) + b LN(x) となるので SLOPE(LN(yの値), LN(xの値)) で 指数の b が傾きとして求められ、 EXP(INTERCEPT(LN(yの値), LN(xの値))) で係数の a が切片として求められます。 指数近似: y = a e^bx a = SLOPE(yの値, EXP(xの値)) LN(y) = LN(a) + bx より b = SLOPE(LN(yの値), xの値) 対数近似: y = a log(x) + b これは直線近似そのまま xの値 を LN(xの値)にすれば問題ない。 また、FORECAST関数を使うと近似直線を使った yの値を直接求められます。この値と元の yの値から計算すると決定係数 R2 を求められます。RSQ関数を使えばもっと簡単。 さて多項式はどうしようか言うヒトのために LINEST関数なるものがあ
続きを読む

塩基配列のデータを扱うときは専用エディタ ApE を使おう

イメージ
塩基配列のデータを扱うときは MS Word じゃなくて ApE を使おう。という話 ApE を使おう ApEとは その名も A plasmid Editor という事で、プラスミドのためのテキストエディタです。とはいえプラスミド以外の塩基配列も扱うことができる優れたソフトです。現在は無料で公開されており、Windows でも Mac でも使うことが出来ます。 MS Word ではダメなのか？ 「餅は餅屋」ということで、ApE をはじめとした専用のエディタには塩基配列の編集に適した機能が豊富に備わっています。MS Word でも同じようなことは出来ますが、無駄な時間と変な工夫をしなくても済むならそちらを使った方が良いでしょう。間違いも少なくなります。 ApEのインストール 早速、インストールしましょう。 http://biologylabs.utah.edu/jorgensen/wayned/ape/ ここから使っているOSのロゴをクリックしてダウンロードしてインストールする。 ApEの初期設定 最初にこれだけは設定しておくと良いです。 メニューから設定画面を開く。Mac だと [ApE] ->  [Preference]（⌘ + ,）。Windows だとたぶん [Edit] か [Tools] から [Preference]。 [Files]タブから、ファイル形式を GenBank、拡張子を .gb にしておく。（もとからこの状態かもしれないけど念のため確認） Save Files in Format: GenBank Default File Extension: .gb Genbank 形式のファイルを編集してみる Addgene のサイトから pUC18 ベクターの .gb ファイルをダウンロードしてみましょう。適当なファイルがある方はそれで試して下さい。 https://www.addgene.org/browse/sequence_vdb/4490/ [Sequence] タブから [Genbank] を選択すると、GenBank 形式の配列が表示されます。配列をクリックしてから、Ctrl + A、Ctrl + C ですべて選択してコピーしてください。 ApE を起
続きを読む

Illustrator で稲妻マークを描く

イメージ
１．長方形を描く（Mキー） ２．一番下の長方形の角を１つ削除して三角形にする（ダイレクト選択ツール：Aキー） 切りっぱなしになるので一応、右上と左下のアンカーポイントを選択して「右クリック」→「連結」しておく。 ３．上辺のアンカーポイントを選択して右にずらす Shift キーを押しながら矢印キーを使うと再現性が良い。ここでは２回分動かした。 ４．右上のアンカーポイントを選択してさらに右にずらす。 ここでは Shift + → １回分。 ５．適当に重ねる ６．「パスファインダー」パネルから「形状モード」の「合体」を選択してまとめる おまけ Illustrator では選択したオブジェクトのみをファイルとして書き出す機能はない（少なくともCS6時点では）。 アートボードを指定したりスライスを作ったりすれば良いけどとてもめんどくさい。 そこで一旦、Powerpoint に貼り付けて配置し、右クリックメニューから「図として保存」する。 これで選択していたオブジェクトのみをPNGやPDFとして簡単に書き出すことができる。
続きを読む