あるケミストの独り言(winchemwinの日記)

ケミスト(化学者)の視点で、面白そうな情報(シミュレーション関係など)を発信

Psi4による量子化学計算ー構造最適化(化学構造の描画)

 前回までの記事で水分子の構造最適化と座標データの表示について紹介させていただきました。今回は最適化構造の可視化について、紹介したいと思います。可視化にあたっては、jupyter notebook 上で可視化できるpy3Dmolを使用しています。py3Dmolについては「化学新しいカタチ」さんのブログで詳細に紹介されています(私も参考にさせて頂きました)ので、参考にして貰えればと思います。

future-chem.com

 

 以下はアセトニトリルの構造最適化とpy3Dmolを利用した可視化のコードになります。座標データについては、前回の記事で記載したようにxyz形式のファイルで保存しています。

 py3DMolですが、conda もしくはpipコマンドで予めインストールしておく必要があります。以下のコードではpy3Dmol のライブラリーをインポートし、さらに表示用のウィジェット(ipywidgets)もインポートしています。

アセトニトリルの構造最適化ー構造の描画

 いろんな化合物を表示するコードを前提にdef で関数設定していますが

、1分子だけの表示であればdef , return の部分は省いてもらってもいいいかと思います。表示は読み込みの座標(xyz), 表示範囲、表示スタイルなどを設定した上でview.show( )で行っています。

 以下は実際の描画結果になります。

最適化後のアセトニトリルの構造

 炭素原子が灰色、窒素原子が青色、水素原子が白色で表現されています。視覚化されることで最適化構造が非常ににわかりやすくイメージできるかと思います。計算で本当に妥当な立体構造に最適化されているかのチェックにも非常に有効です。

 Psi4はpythonコード上で利用可能ですので、上記のようにpsi4だけの機能では難しいことも、python で利用できる他のライブラリーをうまく組み合わせることでいろいろできるようになります。この辺が実験化学者にとってPsi4の大きなメリットの一つと

言えるのではないでしょうか。

 次回は構造最適化後に実施する振動計算について紹介したいと思います。

機械学習による予測モデル作成(我が家の電気使用量 その1-2)

 

 この前の記事で我が家の電気使用量と気象データとの相関について、データの前処理について紹介しました。今回は実際の解析と解析結果の表示について紹介したいと思います。

 

sklearn による解析準備


 まず、機械学習を行うにあたり、目的変数と説明変数を設定しています。先の記事にも記載したようにウィークデイと週末では電気使用量に大きな差がありますので、それぞれのケースで機械学習のモデル作成を検討しています。

 続いてsklearn からモデルのインポートを行っています。今回はシンプルな'linear model'を選択しています。その後、データを学習用の訓練データとテストデータに分割して(train_test_split)準備は完了です。

  以下は休日データを使った学習とその結果の表示になります。

機械学習の実施 1

学習ですがlinear model で回帰モードで学習を行っています(fit)。また学習結果の確認ですが、回帰分析の精度評価指標となる決定係数(R2)と関係係数で確認を行い、プリントアウトしています。

 

 

 

学習結果の表示

 学習結果ですが、見ていただいたらわかるように(上段が訓練データの決定係数、下段が学習データの決定係数)、決定係数としては非常に悪い値となっています。また説明変数の寄与を示す関係係数も算出していますが、あまり良いモデルではないので、意味はないですね。

 次回以降、より良いモデルの作成の検討について、いくつか実施した内容について紹介してゆきたいと思います。

機械学習による予測モデル作成(我が家の電気使用量 その1-1)

 様々なデータの解析に機械学習の手法が取り入れられてきています。python機械学習のライブラリが比較的充実していますので、今回から機械学習によるいくつかのデータの予測モデル作成と検証についての内容を書いてゆきたいと思っています。

 まずは、我が家の電気使用量について行ってみました。我が家は関西電力さんと契約をしているのですが、関西電力さんの「はぴeみる電」のホームページから各自日々の電気使用量をダウンロートできます。そのデータのうち電気使用量を目的変数として、説明変数との関係について機械学習のモデルの作成を検討しました。

 説明変数ですが、なかなか関係が深そうなデータで系統的に集積されているものはないのですが、今回は気象庁で公開されている日々の気象情報との関係について、機械学習で整理できないか検討してみました。気温や湿度はエアコンなどの使用量と関係するので、それなりの相関関係があるのではないかと期待した訳です。

 

 以下はそのコードですが、まずはデータの読み込みと前処理です。それぞれのデータですが、電気使用量は、1ヶ月単位で「はぴeみる電」のホームページからダウンロードできます。

データの読み込みと前処理

 

データは毎日、毎時間の使用量と簡単な気象情報からなっていますが、今回は毎日の電気使用量を目的変数として設定することにしました。データは1ヶ月ごとのデータなので、年間のデータとして予め1年間分のデータをダウンロードして、

集約したファイルを「Year2020_EV.csv」としてエクセルで編集しています。Pandas のデータフレームとしてデータを読みこんでいますが、カラム名が割当られていないのでそれぞれのカラム名を df.renameでわかりやすい名前に割り当てています。

 元のデータの日付ですが、”4/1”のように年のデータが含まれ

ていないので’2020/’を加えた上でpd.to_datatime で日時データに変換しています。また電気使用量はウィークデイと週末とで大きく異なるので、曜日による分類ができるようdt.strftimeで曜日情報を取得した上で、新たにDay/Weekとしてデータフレームに列を追加しています。以上で電気使用量データの前処理は終了です。

 続いて気象データです。気象データも各地の気温、風速、気圧、湿度などがダウンロードできますが、一度にすべてのデータをダウンロードできるわけではないのでそれぞれの個別にダウンロードしたものをエクセルで集約、編集しています。気象データも毎日、毎時間のデータがありますが、今回は毎日のデータを説明変数として利用しています。

 加工した気象データ(Weather2020.csv)はPandas のデータフレームとして読み込んでいます。ここで、電気使用量データと気象データの日

付の並びが異なるので、気象データの方の順番をdf.iloc[::-1]で逆に並び替えて、電気使用量データをあわせる作業を行っています。続いて電気使用量と関係しそうな必要データを取り出して、気象データの前処理作業は終了です。

 最後に電気使用量データと気象データを統合して、機械学習用データ(dfsum)としています。また電気使用量ですが、ウィークデイと週末では大きく異なるので、それぞれ分けて機械学習モデルの作成を行うため、データを分割しています。

  

 以上、説明が長くなりましたが、電気使用量のデータ、気象データの読み込みと前処理までの説明でした。

 次回以降、sklearn による機械学習の実施と検証について紹介したいと思います。

Psi4による量子化学計算ー構造最適化(座標の確認)

前回までの記事で水分子の構造最適化の実施について、紹介しました。

通常は最適化計算後、座標データや構造の確認を行うかと思います。

以下は前回紹介したコードですが、計算のアウトプットファイルは"h2o-opt.txt"に保存されていて、このファイルを確認することで計算の詳細、過程や座標データを確認することができます。

水分子の構造最適化ー(前回のもの)

 以下はoutput ファイル内に記載された内容の一部になります。最適化

された座標データは一番したの欄で確認できます。ただ、いちいち

アウトプットファイルを確認するのも手間がかかりますので、座標データ

を確認できるコードを加えてみました。

 

水分子の構造最適化ーアウトプットファイルの内容

 以下は前回のコードにxyz座標データを取り出すコード(save_string_xyz_file( ) ) を加えたものになります。座標データはxyz形式としてファイル保存するとともに、print 文で表示もしています。座標データのみのファイルがあると別の計算を実施したりする際にも便利ですよね。

 

水分子の構造最適化ー座標データの表示

 上記のコードをjupyter notebook上で実施すると以下のような結果が表示され、最適化後の座標が確認できます。座標だけではどのように最適化されているかはイメージしにくいかもしれませんが、小さな分子であればなんとなくそれぞれの原子がどのように位置しているかは想像できるかもしれませんね。

 

表示された座標データ

 次回は、最適化構造の描画について紹介したいと思います。

Psi4による量子化学計算ーインプットファイルの作成 その2

 

 前回までの記事で、座標データの取得まで説明してきましたので、ここからは、実際の計算用コードの作成に移ってゆきたいと思います。以下は水分子の構造最適化のコードの例です。

f:id:winchemwin:20220107161852j:plain

H2Oの構造最適化の参考コード1


 最初にpsi4のライブラリーをインポート(import psi4)した後、計算に使用するスレッド数とメモリーを設定します(psi4.set_num_threads(nthread=##), psi4.set_memory('###MB'))。私の場合は低スペックのPCで行っていたので上記のコードに記載のスレッド、メモリーですが、PCのスペックに合わせて適宜設定して頂いたら良いかと思います。

 次に化合物の構造情報を記載していきます。座標は前回の記事で記載したxyz形式の座標データを使っています(h2o=psi4.geometry(""" 座標データ"""))。

 続いて、計算結果を保存するためのファイルとしてアウトプットファイルの設定も行った上で(psi4.set_output_file('ファイル名’)、計算手法(今回は’hf/3-21g')、対象分子を設定した上で構造最適化を行うコード(psi4.optimize(計算手法、対象分子)となっています。

f:id:winchemwin:20220107164151j:plain

H2Oの構造最適化結果その1

 

 計算を実行すると最初に設定されたスレッド数とメモリーが表示されます。しばらくして構造最適化に成功するとOptimization completeの表示が現れるとともに、安定化された構造のエネルギーが表示されて計算が終了します。別の化合物を計算する場合は化合物名、座標データ、ファイル名、計算手法などを変えてコードを作成し、計算させることで同様に構造最適化計算を実施できます。

 以上、シンプルに構造最適化する手法を紹介させていただきましたが、次回以降

は最適化構造の座標の表示や化学構造の描画について、紹介したいと思います。

 

 

pythonの活用例について(fit ファイルの利用 その15:Cateye社 Avventuraのデータからーアプリの作成)

 前回までの記事で、fit ファイルをもとにpython で様々な解析を行ってきた例を紹介してきました。それぞれの解析は紹介したコードを参考に処理してもらえればできるかと思いますが、やはりいちいちコードを書いて実行するのは・・・ということでtkinter を使って簡単なアプリを作ってみました。下記が作成したアプリを立ち上げた画面になります。

 使い方ですが、「解析用のファイルを選択」と書かれたところの下にある「選択」ボタンを押して、解析用のfit ファイルを選択します。選択したファイルが上の空欄のところにディレクトリも含め表示されたらOKです。

 次に選択したfitファイルの解析範囲を設定するため(解析範囲を距離で設定)、その下の「解析範囲を設定」ボタンをおします。

アプリ メイン画面

 次に選択したfitファイルの解析範囲を設定するため(解析範囲を距離で設定)、その下の「解析範囲を設定」ボタンをおします。すると別画面に標高/距離のグラフが表示されます。(下図)グラフを参考にしながら、「解析開始点」、「解析終了点」を入力して解析したい範囲を指定します。数値を入力し終えたら、右下の「実行」ボタンを押して、解析範囲の設定を実行します。

解析範囲の設定画面

 解析範囲の設定が終わると、下図のように解析範囲がメイン画面に表示されます。

 

設定範囲の表示(解析範囲設定の結果)

 次にメイン画面の「作成するグラフ、データ等を選択してください」書かれている下のプルダウンメニューから行いたい処理を選択します。選択した処理はその下の「実行」ボタンを押すことで実行されます。解析結果は別画面で例えば下図(いくつかあります)のように表示されます。解析結果の表示画面で表示画面の拡大縮小や画面の保存も行うこともできます。

距離/速度、勾配のグラフ

 

ヒストグラム(勾配)

 

バイオリンプロット(勾配)

激坂ファクター

 読み込んだファイルはメイン画面の一番下の「データ保存」のボタンを押すことでcsv形式のファイルに保存することも可能です。

 シンプルなアプリですが、よかったら試してみてください。

 

 尚、このアプリのコード(英語版)は下記のgithub にMITライセンスで公開していますので、そちらから利用してみてもらったらと思います。

 

github.com

 

Psi4による量子化学計算ーインプットファイルの作成 1

デモのプログラムでPsi4の動作が確認できたら、いよいよ実際の掲載に移っていきます。

Psi4では化合物の座標情報として、デカルト座標, Z-matrix 系 のインプットの入力が可能です。ここではデカルト座標系での入力で計算を実施してゆきたいと思います。化合物の構造から座標データを得るには、モデリングソフトなどを利用するのが簡単です。市販の有料ソフトウェアを使って作成することは可能ですが、無料で簡単に行うとなるとWeb上で構造式描画を行えるMolviewがひとつの選択肢になります。

 Molvewについては下記のリンクからアクセスすることができます。

 

 MolView

 

 リンク先にアクセスすると下記の初期画面が現れます。利用規約などを確認した上で”Close"をクリックすると次の初期入力画面(カフェインの構造がデフォルトで入力されている)が現れます。

f:id:winchemwin:20211212154453j:plain

Mol View初期画面

 

 

 

f:id:winchemwin:20211212155543j:plain

初期入力画面

 MolViewの使い方ですが、初期入力画面からカフェインの構造を消去した後に、作成したい構造を書いていきます。下記は入力画面ですが、左側に必要な構造のパーツがあり、右側に主要元素が並んでいます。化学系の描画ソフトを使ったことのある方なら、それほど違和感なく書きたい構造をかけるのではないかと思います。構造の書き方はこのブログでは詳しくは説明はしませんが、それほど難しいソフトではないので触ってなれてみてもらえたらと思います。構造がかけたら、3Dデータを作成するため、上部のメニューバーにある2D to 3Dをクリックします(これが大事なポイントです)。すると右側の画面に描画した分子の3Dの分子モデルが表示されます。

f:id:winchemwin:20211218040435j:plain

入力画面

  下記は「酢酸エチル」を入力して、3D表記させた例になります。

f:id:winchemwin:20211218042559j:plain

入力例(酢酸エチル)

 構造を作成できたら、Mol ファイルへの保存を行います。ToolのプルダウンメニューからExport - MOL file を選択します。

 

f:id:winchemwin:20211218043037j:plain

Mol file での保存 1

 保存を行うとソフトウェアの下部に保存したファイル名が表示されます。表示されたファイルをダウンロードして自分のPCにファイルを保存することで、作成した構造のMol fileを取得することができます。

 

f:id:winchemwin:20211218043628j:plain

Mol file での保存 2

 以上でMol file を作成することができましたので、次にインプットファイルの形式に対応したデカルト座標のデータを得るため別のソフトウェアで処理を行います。単純にファイル形式の変更ですので、mol ファイルからxyz形式に変更できる手法であればここで紹介する手法以外(openbabelなどの利用)でも問題はないかと思います。ここでは無料利用可能なAvogadroを使った方法を紹介します。Avogadroについては下記リンクからアクセスすることができます。

 

Avogadro - Free cross-platform molecular editor - Avogadro

 

内容確認の上、使用OSにあったプログラムをダウンロードしてもらえればと思います。

f:id:winchemwin:20211226035936j:plain

Avogadro ホームページ画面

導入後、ソフトウェアを立ち上げると下記のような画面が現れます。

 

f:id:winchemwin:20211226062805j:plain

Avogadro 立ち上げ画面

立ち上げ後、作成したmolファイルを読み込み、改めてxyzとして保存します。

f:id:winchemwin:20211226063314j:plain

保存画面

保存したファイルを確認すると下記のようにデカルト座標系のデータとなっていますので、このデータを用いてPsi4のインプットファイルを作成してゆきます。

f:id:winchemwin:20211226063711j:plain

座標データ(H2Oの場合) 

実際のPsi4のインプットの作成は次の記事で紹介します。