#author("2024-03-08T07:37:31+09:00","","")
#author("2024-03-08T07:43:18+09:00","","")
*GULPの最初の一歩
#contents
----
**GULPについて
(ここで説明しているGULPは格子力学などの計算を行うもので、Javascriptとは無関係です）
GULP(General Utility Lattice Program)はJulian D. Gale教授（Curtin University） が開発したプログラムで格子力学計算などを行うものです。沢山の機能があり、分子動力学も可能ですが、作者自身はメインの機能だとは見做してないようです。様々な経験的ポテンシャルが使えますが、個人的にはイオン結晶でshell modelが使えるところが特徴だと思います。ポテンシャルの最適化も可能です。物性研のMateriAppsでも紹介されています。以下は作者の説明文(gulp.introから）。
>GULP is a program designed for performing lattice dynamical calculations on solids, defects, surfaces, interfaces, polymers and molecules. For 3D solids calculations can use symmetry to generate structures and to save CPU time where this is useful. The emphasis on the program is on analytical solutions rather than on molecular dynamics.
**GULPの機能（これもGULPのIntroから）
-系：3-Dバルク、0-Dクラスター、3-D中の欠陥、1-Dのポリマー、表面(2-D)、界面と粒界(2-D)
-エネルギー最小化：定圧、定積、シェルだけの緩和（光学的）、基準配置に対する有限ひずみの最適化、 四元数を用いた剛体分子、ブリーシングのみの緩和、対称性を考慮した緩和、非拘束緩和、内部およびセル座標の拘束、Newton/RaphsonまたはRational Function最適化、HessianまたはそのインバースのDFPまたはBFGS更新、限定メモリBFGS
-自由エネルギー最小化：２体、3体、４体ポテンシャルに対する解析的導関数
-遷移状態：n次定常点の位置、モードフォローイング
-結晶物性
-結晶物性（以下の２項目）
-弾性定数：バルク、ヤング、シアーモデュラス、静的誘電率
-フォノン周波数：フォノン状態密度、射影フォノン状態密度、分散曲線、０点振動エネルギー、エントロピー、比熱、ヘルムホルツ自由エネルギー、ボルン有効電荷、バルク、ヤング、シアーモデュラス、ポアッソン比、S,P波速度、フォノン平均運動エネルギー、熱伝導率
-欠陥物性：欠陥エネルギー、欠陥移動の遷移状態、欠陥周波数
-クラスター：全ての計算がクラスターにもできる（意味があれば）
-分子動力学：有限質量または断熱法のshell model MD、NVE, NVT,NPTアンサンブル、PLUMED2プラグインを使った自由エネルギー計算
-モンテカルロ：固定Nまたはグランドカノニカル、剛体分子のオプション
-ポテンシャルへのフィット：計算された全ての物性への経験的フィット、shell位置と半径の同時緩和、勾配でなく変異にフィット、複数の構造に同時フィット、core/shell電荷分割を変化させる、全電荷を変化させる
-ポテンシャルライブラリー：標準ポテンシャルのライブラリーを持てるオプション
-shell model：ダイポーラーとスフェリカルブリーシングshellが可能
-Embedded atom method：合金スケーリングを含む金属に対する方法が使える
-エネルギー超曲面へのフィット：QM由来エネルギー曲面のフィット、ポテンシャルフィット
-電気陰性度平衡法：Mortiersの方法を使ってケイ酸塩と有機物系の電荷分布を求める、Rappe and Goddard III の QEq スキームを使用する、Streitz-Mintmireスキームを使用する
-構造解析：結合長、結合角、捻れ角、密度、体積を計算する
-構造の操作：スーパーセルの作成
-静電サイトポテンシャル
-電場勾配と非対称性パラメータの計算
-他プログラム用の入力ファイル生成：GDIS, marvin, Cerius2, Insight (.xtl, .arc), Materials Studio, Alamode
**GULPの入手方法
[[GULPのページ:https://gulp.curtin.edu.au/gulp/]]でまずレジストレーションする必要がありますが、私の場合はここで問題が起こりました（2021年11月）。大学のメールアドレスを使う必要がありますが（アカデミックかどうかのチェックのため？）、アドレスにハイフンが含まれていると自動処理がうまくいかないそうです(okayama-uとかダメ)。そのために結局、Gale先生に直接連絡してレジストレーションして頂きました。この状況がまだ続いているかどうかは不明です。レジストレーションされると後はメールアドレスを入れると、そのメールアドレスにダウンロードリンクが送られてきます。2021年11月に私がダウンロードした時は6.0でした。最新は6.1.2です（20230626)。tgzがダウンロードされるので解凍します。
**GULPのコンパイル
Windows版はコンパイルしたものがあるようですが、Mac, Linuxはコンパイルする必要があります。Fortran90コンパイラが必要です。gfortranでコンパイルできます。MPIで並列化もできますが、私の場合はシリアルでコンパイルしています。XtalOptのoptimizerとして使う予定があったのですが、XtalOptは同時にGULPに複数のジョブを投げるので、結局GULPが複数並行して起動されることになり、GULP自体を並列化する意味があまりないのが理由です。しかしGULP単独で重い計算をする場合は並列化した方がいいかもしれません。
Srcディレクトリに移動して、ここにF90ソースとコンパイル用のスクリプトmkgulpがあるので、ここで
./mkgulp
を実行する。シリアル版だとこれだけでOK。MPIを使う場合などのオプションはREADMEに書かれている（マニュアルにも）。
**GULPの使用方法
マニュアルもありますし、サンプルも多数あるので、まずはそれらを読むのがいいと思います。実際の入力ファイルのフォーマットにはhelp.txtかhelp.htmlが役に立ちます。また、Librariesの中にポテンシャル例が沢山あります。SiO2だとTsuneyuki, Vashishtaと懐かしい名前が。shell modelもあります。もちろん自分の計算したい原子のポテンシャルがないことが多いと思います。酸化物だと酸素のパラメータはLibrariesの中から探してきて、陽イオン側は結晶構造からポテンシャルパラメーターを適当に決めることが可能です。ただ注意が必要なのは、構造からポテンシャルパラメーターを決めると、ポテンシャルミニマム位置は求まり易くなりますが、ポテンシャルミニマム付近の傾きが制約されません。傾き（微分）が重要となる、弾性定数、格子力学などではあまり精度がよくないことになります。GULPは弾性定数など入れてもポテンシャルパラメーターを決めることができたはずです。ただし、弾性定数のよいデータはなかなか入手難しいかもしれません。鉱物関係だと結構測られてますが。
プログラムは
> gulp < test.gin > test.got
で実行できますが、実際ターミナルで実行してみると１つ問題がありました。実行結果を上記のようにリダイレクトした場合に計算の最後の方が書き込まれていませんでした（Macで）。ただ、標準出力でターミナル上に書き出すと問題ありません。そこでファイルに書き込む時は
> gulp < test.gin | cat > test.got
と一度パイプでcatに渡してからリダイレクトしたら全て書き込まれてました。
**forsteriteのポテンシャル最適化と構造最適化の入力ファイル
例としてforsterite (Mg2SiO4)のポテンシャル最適化と構造最適化の入力ファイル例を下に示します。これは酸素はshell modelで扱っていて、O-Si-O３体ポテンシャルも入ってます。catlow.libにあるポテンシャルからMg, Si, Oのポテンシャル初期値を取って来てます。
最初の行にどういう計算をするかのキーワードが入ってます。optiは構造最適化、fitはポテンシャル最適化、conpはconstant pressure、simulは同時にフィットすることを意味します。この例のようにoptとfitが同時にある場合は、ポテンシャルをまず最適化して、それを使った構造最適化が行われます。指定したポテンシャルで構造最適化するだけの場合にはfitはいりません。
２〜３行は格子定数、その次行から結晶中原子の内部座標が入ります。空間群を使うので格子内の全ての原子位置は必要ありません。この例の場合は酸素にshell modelを使っているために、同じ酸素の位置が２つあって、片方はcore、他方がshellになります。最初位置は同じですが、最適化するとcoreとshellは一般には少しズレます。shell model使わない場合でも元素記号の後にcoreと書いておきます。１つの原子に対してcoreとshellとほぼ同じ位置の座標があるので、VESTAなどで構造を描かせたときにちょっと変なことになりますし、XtalOptなどはshell modelはサポートしてないようです。
座標値の後には空間群の指定があります。番号かHM記号で指定。この例では標準軸設定ではないPbnmを指定しても問題ありませんでした。なので非標準の設定でも大丈夫なようですが（その場合は番号ではダメ）、origin選択ある場合どうするのかはよく分かりません。標準に直しておくのが無難なようです。
その次のspeciesは原子種で価数を定義します。ここでも酸素の場合にはcore, shell(shelとする）があることに注意。このcatlow.libでは形式電荷を使っていて、この酸素の場合もcoreとshellを足すと形式電荷（-2）となります。
次はポテンシャルの指定で、ここではよく使われるバッキンガムタイプの２体間ポテンシャルを使ってます。shell modelを使ってますので、酸素のshellの方がこのポテンシャルの対象となってます。Mg-O, Si-O, O-Oの３つの２体間ポテンシャルが定義されてます。最初のパラメーターがpre-exponentialのAパラメーターで、次がexpの中のρパラメーター、３つ目が分散のパラメーターCになります。４つ目は近い側のカットオフ距離、５つ目が遠い側のカットオフ距離。その後の３つはfitの時にどのパラメーター（最初の３つ）を最適化するかのフラグで、最適化する場合は1にします。fit自体しない場合は全て0にします。あまり一度に多くを最適化するとうまくいきません。うまく収束しても、あるパラメーターが大きく変わることがあります。それらの場合は、少しづつ最適化します。
その次のspringはshell modelで使われていて、酸素のcore-shell間がspringでつながっていると考えているので、そのバネ定数がここに入ります。次の0/1はそれを最適化するかどうかのフラグです。なおポテンシャルはLibrariesのファイルを指定して読み込むこともできたはずです（まだ試したことはありません）。
threeは3体間ポテンシャルのことで、この例ではO-Si-Oについて適用してます（酸素はやはりshellに対して）。当然4配位Siでのみ有効（といっても場合によっては配位数変化が起こります）。これもバネ的なポテンシャルで、最初が変角に対するバネ定数、次が平衡なO-Si-O角度、カットオフ距離、最初２つのフラグとなります。
圧力を指定したい場合はPressureで出来ます。それに応じて格子定数が変化しますので、一連の圧力を変えた計算から圧縮率が計算できます。指定しないと1 barとなります。
最適化した構造はcifで出力できます。一番最後のoutputがそれの指定。shellの位置は出力されません。cifファイルはVESTAで読むことができます（なぜかP1で出力される場合もあります）。出力されたcifはたまに変な感じになる時がありました。
fit simul opti conp
cell
4.756000 10.207000 5.980000 90.000000 90.000000 90.000000
fractional
Mg core 0.000000 0.000000 0.000000
Mg core 0.991500 0.277400 0.250000
Si core 0.426200 0.094000 0.250000
O core 0.765700 0.091300 0.250000
O shel 0.765700 0.091300 0.250000
O core 0.221500 0.447400 0.250000
O shel 0.221500 0.447400 0.250000
O core 0.277700 0.162800 0.033100
O shel 0.277700 0.162800 0.033100
space
P b n m
species
Mg core 2.00000
Si core 4.00000
O core 0.86902
O shel -2.86902
buckingham
Mg core O shel 946.627 0.31813 0.00000 0.0 10.0 1 0 0
Si core O shel 1283.907 0.32052 10.66158 0.0 10.0 1 0 0
O shel O shel 22764.000 0.14900 27.87900 0.0 12.0 0 0 0
spring
O 74.92 0
three
Si core O shel O shel 2.09724 109.47 1.9 1.9 3.5 0 0
pressure 0 GPa
output cif forsterite-fit
出力ファイルや構造をチェックして、おかしい場合は入力パラメーターを変えてみる必要があります。この例の場合だと最初３つのAバラメーターを最適化したら、O-OのAがトンデモなく大きな値になりました。多分おかしいので、パラメーターは少しづづ変えていく必要があります。
この例では構造は１つしかありませんが、複数の構造を入れることもできます。その場合はcellのところからspaceのところまでを構造毎に入れていけばいいだけです。その場合は同時にフィットすることになります。最高７構造まで入れてフィットしたことがあります。ただし、出力されるcifファイルは複数の構造が入ったものになります。VESTAは最初の構造だけしか読みません。必要ならエディタで構造毎に切り分けてやる必要があります。このcifはちょっと変なところがあって、単純格子でない場合には単純格子に変換して計算されますが、出力されるcifで格子定数と座標値の格子が一致しないケースがありました。仕方ないので、出力ファイル(got)から複数の構造をcifにして取り出すプログラムをPythonで作ってます。ちょっとまだ不具合があるので、公開まで至ってません。
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")