CP2Kは、量子化学、固体物理学のシミュレーションを行うためのオープンソースのパッケージです。
WikipediaのList of quantum chemistry and solid-state physics softwareの中のリスト、Quantum chemistry and solid-state physics characteristicsでは、ちょっと前(2018年ころ?)までですが、すべて緑色のYESで埋まっていたのは、商用版のGaussianとこの無料版のCP2Kくらいでした。
(2021年5月23日現在では、機能を示す項目が増えたため、すべてYESのソフトウェアはなくなっています。)
ともあれ、現在でも商用版に劣らない、多くの機能を備えた量子化学、固体物理学のシミュレーションパッケージといえます。
(追記)
Wikipediaですが、日本語の「量子化学および固体物理計算ソフトの一覧」ですと2021年6月19日現在、以前のままでした。
Linux版のバイナリは、2021年5月23日現在、v5.1からv8.1まで公式githubで配布されています。
以前はMac用のバイナリの配布もあったようなのですが、現在は配布されていないようです。
Android用のバイナリは、Alan Liškaさんによるアプリが、Google Playで配布されています。
Windows版のバイナリは、私のgithubで配布しています。
ここではWindows版について説明します。
64ビットネイティブビルドは少なく、ほとんどが32ビットビルドですが、いずれも「2GBを超えるアドレスを処理する」(Visual Studioにおける/LARGEADDRESSAWAREオプション付きでのコンパイルに相当)設定でのビルドのため、32ビットビルドであっても、64ビットPCでは、2GB以上のメモリを使用することができます。
しかしながら、現在Windows PCをお使いの方の多くは、64ビットWindows10をお使いでしょう。
その場合は、タグ「0.0.220531」のリリースをお使いになるのがよいと思います。
インストーラは設けていません。いわゆるportable版のようなアプリケーションになっています。
以下の図はgithub releaseのタグ「0.0.200301」ですが、assetsの中のzipファイルをダウンロードして、zipファイルを展開すれば、インストール完了です。
また、展開したフォルダごと削除すれば、アンインストールされます。
CP2Kに計算をさせるための入力ファイルを作成します。
入力ファイルはテキストファイルです。セクション、サブセクション、キーワードで構成されます。
セクションの中にサブセクションがあり、サブセクションの中にキーワードが入ります。
サブセクションの中にさらにサブセクションが入ることもあります。
サブセクションやキーワードはどこに書いてもいいわけではなくて、それが属する親セクション、親サブセクションの中に書かないといけません。
セクション、サブセクションは、&SECTION_NAMEで始まり、&ENDで終わります。
実際書いてみると、&ENDだらけになりますから、どの(サブ)セクションの終わりか明示するために&ENDでなく&END SECTION_NAMEなどと書いたほうがよいでしょう。
キーワードはなんらかの値をもちます。
これだけだとなんだかわからないですから具体的に書いてみましょう。
「Avogadro 1.x に原子や分子を表示する 」の「マウスで一から作る 」では、ピリミジンの分子モデルを作成しました。
このとき分子構造を古典力場で最適化しましたが、ここではもう少し精度を上げ、量子力学(量子化学)的計算で構造最適化をする入力ファイルを作成してみましょう。
構造最適化とは、人間が、分子モデリングソフトウェアなどで適当につくった分子構造を、現実の構造に近づける作業のことです。
CP2Kの計算では、GLOBALセクションはどんな計算でも必ず必要です。
このGLOBALセクションの中で、プロジェクト名やランタイプ(エネルギーの一点計算、構造最適化、振動解析、動力学...などなど)などをキーワードの値で指定します。
&GLOBAL
PROJECT pyrimidine
RUN_TYPE GEO_OPT
&END GLOBAL
PROJECTやRUN_TYPEがキーワードで、それぞれpyrimidine、GEO_OPTという値を指定しています。
プロジェクト名はピリミジンの計算であることがわかるようにpyrimidineとしています。
構造最適化計算の場合はGEO_OPTを指定します。
これがたとえば指定した分子構造でのエネルギー一点計算でしたら
RUN_TYPE ENERGY
と指定します。
空行があってもいいです。インデントもあったほうがわかりやすいと思います。
次ですが、たいていの計算にはFORCE_EVALセクションが必要です。
FORCE_EVALとは力の評価ですが、力は原子間に働く力のことです。
すなわち電子状態のことを指しますが、古典力場も指定できます。
量子力学的計算では次のように、
・DFTサブセクションで、計算方法を
・SUBSYSサブセクションで、分子構造を
それぞれ指定します。
&FORCE_EVAL
&DFT
&QS
METHOD xTB
&END QS
&END DFT
&SUBSYS
&CELL
ABC 20.0 20.0 20.0 20.0
&END CELL
&COORD
C 1.2502607908 1.4433052443 0.0000000000
C -0.9403223106 -0.0606634697 0.0000000000
C -0.0141333425 2.0100033698 0.0000000000
C 1.3248470438 0.0597315438 0.0000000000
N -1.1079366402 1.2558521162 0.0000000000
N 0.2280679834 -0.6900854069 0.0000000000
H 2.1427366348 2.0560442778 0.0000000000
H -1.8361739357 -0.6757201742 0.0000000000
H -0.1584509209 3.0866295677 0.0000000000
H 2.2777356046 -0.4617614311 0.0000000000
&END COORD
&END SUBSYS
&END FORCE_EVAL
計算方法は、DFTサブセクションとQSサブセクションの中に、METHODキーワードを置いて指定します。
ここではxTBを指定します。
xTBは、計算負荷の少ない、半経験的な密度汎関数法計算である「DFTB法」の中でも、特に「分子構造(Geometry)」「分子振動数(Frequency)」「非共有結合性相互作用(Noncovalent interaction)」の再現性に優れた方法で、頭文字をとって「GFN-xTB」とも呼ばれています。
ここでは分子構造の予測が目的ですから、この方法で計算してみます。
半経験的な方法ですので、パソコンでも十分計算可能なくらい軽量な計算です。
次に、SUBSYSサブセクションの中に、
・CELLサブセクションで、シミュレーションボックスのサイズを、
・COORDサブセクションで、分子構造(構造最適化の場合は初期構造)を
それぞれ指定します。
シミュレーションボックスは、分子1個が入るサイズであればいいと思うかもしれないですが、あまりギリギリの小さなボックスですと、計算の収束が悪くなりますので、大きめに設定します。オングストローム(Å)単位です。
分子構造は
原子名 X座標 Y座標 Z座標
の値を、これもオングストローム(Å)単位で書きます。
これはAvogadroの値を使うこともできます。
Avogadroで分子を表示した状態で、「ビルド(B)」→「Cartesian Editor...」を選びます。
カルテシアンエディタが開きます。
左下のドロップダウンリストは、「Angstroms」と「XYZ」を選びます。次に右下のアイコンの「紙が二枚重なっている」ボタンの「Copy All」を押します。
クリップボードにコピーされますので、入力ファイルの&COORDと&END COORDの間にペーストします。
入力ファイルは、エネルギー一点計算、構造最適化などでは、これだけが必要最小限です。
これだけ指定すれば、とりあえずほかの設定はデフォルト値で、計算を実行できます。
CP2Kの全てのセクション、サブセクション、キーワードを網羅したリファレンスマニュアルは
にあります。
実行方法は「8.0まで」(タグ0.0.200331まで)と「8.1以降」(タグ0.0.200313、0.0.200313-2以降)とで少し違います。
ここでは「8.1以降」の実行方法からまず説明します。
ところで、入力ファイルをもう一度書くと
&FORCE_EVAL
&DFT
&QS
METHOD xTB
&END QS
&END DFT
&SUBSYS
&CELL
ABC 20.0 20.0 20.0
&END CELL
&COORD
C 1.2609286369 1.4506293911 0.0000000000
C -1.0074419935 -0.1067452373 0.0000000000
C -0.0058826261 2.0342516809 0.0000000000
C 1.3504406104 0.0587193314 0.0000000000
N -1.1141420092 1.2481700536 0.0000000000
N 0.2186693732 -0.6931164381 0.0000000000
H 2.1531665541 2.0632054576 0.0000000000
H -1.8988717038 -0.7187660865 0.0000000000
H -0.1087390082 3.1111638215 0.0000000000
H 2.3184921054 -0.4241838669 0.0000000000
&END COORD
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT pyrimidine
RUN_TYPE GEO_OPT
&END GLOBAL
のようになっていると思います。
これを「pyrimidine.txt」という名前のテキストファイルで保存し、ダウンロード展開したzipファイルから生成したフォルダの中に入れます。
次にコマンドプロンプトを開きます。まず「Windowsキー」と「R」キーを同時押しで「ファイル名を指定して実行」ダイアログを開きます。
「cmd」と入れて「OK」ボタンを押します。
コマンドプロンプトが開きます。
コマンドプロンプトのカレントディレクトリを、CP2Kのディレクトリに変更します。
どんなやり方でもいいのですが、私がよくやる方法は、プロパティダイアログの「場所」情報を使う方法です。
まずは、CP2Kのディレクトリの中のファイルどれでもいいので、右クリックでプロパティダイアログを開きます。
プロパティダイアログの「場所」をマウスで選択し、右クリックでメニューを出して「コピー」を選び、「CP2Kのフォルダの場所をクリップボードにコピー」します。
コマンドプロンプトに戻ります。「cd 」(cdの二文字とスペース一文字)を打ち込んでから「CP2Kのフォルダの場所」を貼り付け(ペースト)します。
エンターキーを押して実行すると、コマンドプロンプトのカレントディレクトリが「CP2Kのフォルダ」に変わっていると思います。
カレントディレクトリがCP2Kのフォルダ(ディレクトリ)に変わっていることを確認したら
startcp2k -i pyrimidine.txt -o pyrimidine.out
を打ち込んでエンターエンターキーを押し、計算を実行します。
コマンドの意味ですが、「命令「startcp2k」を、入力ファイル「pyrimidine.txt」で実行し、結果を「pyrimidine.out」というファイルに出力せよ」という意味です。
割とあっさり計算は終わったのではないでしょうか。
計算が終わると、コマンドプロンプトは入力を受け付ける状態に戻ります。
また、CP2Kフォルダの中に次のようにファイルがたくさん生成します。
pyrimidine.outは、こちらで計算時に指定した名前の出力ファイルで、テキストファイルです。どんな計算が行われたかの情報が書かれています。
肝心の最適化された分子構造は、pyrimidine-pos-1.xyzファイルに、xyz形式のテキストファイルで書かれています。
それ以外のBFGS.Hessian、restart(.bak)、wfn(.bak)ファイルは、バイナリファイルで、エディタなどで見ても「?」です。構造最適化計算の際のHessian情報、再計算をする際の情報などが書かれています。
ところで、CP2Kが出力するテキストファイルは、メモ帳、ワードパッドなどで見ることができるのですが、改行コードがUNIX式のLFになっており、ワードパッドで見る場合は問題ないのですが、メモ帳で見ると改行なしの一行に書かれているように見えてしまいます。メモ帳で見るには、Windows式のCR+LFにする必要がありますが、このための簡単なツールを私のgithubに準備しています。
unix2winディレクトリの中に、drag_and_drop_on_me.batというバッチスクリプトがあります。中身はすごく単純なものですが、これにUNIX式改行コードのテキストファイルをドラッグ&ドロップすると、もとのテキストファイルと同じディレクトリに「ファイル名_win.txt」というWindows式に改行コードが変換されたテキストファイルが生成します。便利だと思いますので、ぜひダウンロードの上、お試しあれ。
出力された「pyrimidine-pos-1.xyz」を(unix2win/drag_and_drop_on_me.batで変換した後)開いてみましょう。
私が計算させた結果では、6回目の計算(i=6)で最適化されました。
10
i = 6, E = -16.9539191122
C 1.2502607908 1.4433052443 -0.0000000000
C -0.9403223106 -0.0606634697 0.0000000000
C -0.0141333425 2.0100033698 0.0000000000
C 1.3248470438 0.0597315438 0.0000000000
N -1.1079366402 1.2558521162 -0.0000000000
N 0.2280679834 -0.6900854069 -0.0000000000
H 2.1427366348 2.0560442778 -0.0000000000
H -1.8361739357 -0.6757201742 -0.0000000000
H -0.1584509209 3.0866295677 0.0000000000
H 2.2777356046 -0.4617614311 -0.0000000000
これをAvogadroのカルテシアンエディタに貼り付けて構造を見ることができます。
カルテシアンエディタのペーストボタンか、右クリックのポップアップメニューから貼り付けなどで、カルテシアンエディタに座標値のテキストを貼り付けます。
貼り付けたら「Apply」ボタンを押して適用した後、右上のクローズボタンでカルテシアンエディタを閉じます。
UFFで最適化されたときとの構造の違いは、見た目ではわからないですけれどねえ・・・
以上が「CP2K 8.1以降」の使い方です。一番簡単な例です。
「8.0まで」もほとんど変わらないのですが、zipファイルを解凍して生成するフォルダ構造が結構違うので、とまどうかもしれません。
「8.0まで」を使う場合、注意が必要なのは、入力ファイルの書き方です。
上の図でも見えますが、「8.0まで」には、dataというフォルダがあって、計算に必要な各種パラメータファイルがこの中に格納されています。
このdataフォルダを、CP2Kがうまく見つけてくれない不具合があるので、たとえば、上で行ったpyrimidine.txtの場合、明示的にdataフォルダを指定するため、次の赤字の部分を追加する必要があります。
&FORCE_EVAL
&DFT
&QS
METHOD xTB
&xTB
&PARAMETER
PARAM_FILE_PATH ./data/
&END PARAMETER
&END xTB
&END QS
&END DFT
&SUBSYS
&CELL
ABC 20.0 20.0 20.0
&END CELL
&COORD
C 1.2609286369 1.4506293911 0.0000000000
C -1.0074419935 -0.1067452373 0.0000000000
C -0.0058826261 2.0342516809 0.0000000000
C 1.3504406104 0.0587193314 0.0000000000
N -1.1141420092 1.2481700536 0.0000000000
N 0.2186693732 -0.6931164381 0.0000000000
H 2.1531665541 2.0632054576 0.0000000000
H -1.8988717038 -0.7187660865 0.0000000000
H -0.1087390082 3.1111638215 0.0000000000
H 2.3184921054 -0.4241838669 0.0000000000
&END COORD
&END SUBSYS
&END FORCE_EVAL
&GLOBAL
PROJECT pyrimidine
RUN_TYPE GEO_OPT
&END GLOBAL
この入力ファイル「pyrimidine.txt」をCP2Kのフォルダに入れる点は、「8.1」と同じです。
「8.1」と同じように、コマンドプロンプトを開き、カレントディレクトリをCP2Kのフォルダ(ディレクトリ)に変更します。
計算実行のコマンドは、ちょっと違って
cp2k.sopt -i pyrimidine.txt -o pyrimidine.out
と「startcp2k」ではなく「cp2k.sopt」とします。
これで「8.1」と同じように計算することができます。
(ただし、xTBは7.1以降にしかありませんので、あしからず)