#author("2024-06-21T10:51:04+09:00;2024-03-08T20:40:20+09:00","default:masami","masami")
#author("2024-08-02T13:58:53+09:00;2024-03-08T20:40:20+09:00","default:masami","masami")
*lammpsメモ(2023/06/15：更新2024/6/21更新)
#image(https://mkanzaki.sakura.ne.jp/images/OVITO_cristobalite.png,center,60%)
#contents
**概要
lammpsは非常にポピュラーな分子動力学法(MD)のプログラムです。最近必要に迫られてちょっと使い始めたのですが、大規模なシステムでマニュアルは2000ページ以上。なかなか初心者には取っ付きにくいところがあります。色々調べて何とか結晶のシミュレーションができるようになったので、そのメモを残しておきます（自分もすぐ忘れるので）。無機結晶のMDシミュレーションが対象です。lammps以外にはOVITOという可視化ソフトを使っています（一番上の図はOVITOの画面）。OVITOはlammpsの入力ファイルとなる初期構造も作ることができるので、最低限OVITOだけあれば取り敢えず計算ができます。それ以外にVESTAも使ってます。以下MD計算自体は知っていることを想定。MD計算については最近だと「分子動力学法と原子間ポテンシャル」（渡邉孝信著・森北出版）が出版されてます。lammpsのインストール方法やReaxFF、GAPについても説明があります。
2023年8月にlammpsの「Virtual Workshop and Symposium 2023」が開催され、その時の発表は今後個々の講演毎に整理されて2021年の時と同様に見られるようになる予定ですが、現時点ではyoutubeで各日の講演全体が見られます。ただ、整理されてないそのままの動画なので、冒頭２０分以上始まらなかったり、昼休み１時間もそのまま流しているので、適当に飛ばして見ないといけません。2024年6月21日時点でもまだ整理されていませんでした…
2023年8月にlammpsの「Virtual Workshop and Symposium 2023」が開催され、その時の発表がやっとまとめられて、個々の講演（ハンズオンは除く）が見られるようになってました。また、2024年6月にMaster Classが開催されたようです。現時点ではyoutubeで内容が見られるようです（まだちゃんと見てません）。その時に使った内容はダウンロードできるようです。
**インストール
lammpsは、[[lammpsのサイト:https://www.lammps.org/]]からダウンロードできます。今回はMacBook Airに実行版をインストールしました。pythonの関係でanacondaを使っているので、私はcondaを使った方法でインストールを行なっています。brewでインストールすることもできます。ただ、現在のところMacの場合インストールされるのはIntel CPU版だけのようです。特に指定しなくても１core用(lmp_serial)とmpiを使った複数core用もの(lmp_mpi)の２種のlammps実行版がインストールされます。Linuxでソースからインストールする場合は最後をご覧ください。
OVITOのことを知ったのはlammpsの「Virtual Workshop and Symposium 2021」の内容を見ていた時でした。OVITOの発表があって、構造など綺麗に可視化できるようだったので、インストールして使ってます。2023でもOVITOの発表はありました（2021とそれほど変わらず）。OVOTOは[[こちら:https://www.ovito.org/]]からダウンロードできます。BasicとProがあります。Proはライセンスが必要です。Proのアカデミックシングルライセンスは年間¥72,428となってました。毎年払うのはちょっときついですね。Mac版のBasicをダウンロードしました。Linux版も使ってます。
**結晶構造ファイルの作成
金属やダイヤモンドなどの簡単な構造については結晶構造ファイルを用意する必要はありませんが、ほとんどの無機構造では結晶構造ファイルをまず用意する必要があります。最初cifなどから構造を変換してくれるAtomsk、lammpsのpythonコードなどウェブ上で紹介されているソフトを試してみたのですが、Macでは動かないなどうまく行きません(WindowsではAtomskは問題なかった）。ところが可視化ソフトのOVITOでlammpsで使う構造を作成することができると分かったので、今はOVITOで作成しています。
まずは必要な構造のcifファイルをCODなどから探してダウンロードします。cifファイルとは結晶構造情報をまとめたファイルのことで、結晶構造データベースや結晶構造を報告した雑誌のsupplemental materialsなどとして入手できます。CODは結晶構造データベースで、完全freeなのでお勧めです。CODについては別のページで説明しています。
cifファイルを得たら、OVITOからLoad Fileでcifファイルを読み込ませます。ちゃんと構造が見えていることを確認します。そしてFileメニューのExport fileで、LAMMPS Data Fileを選んで出力します。ここで出力オプションの選択が出ますが、通常はatomicでいいのですが、電荷を計算中に変化させるようなポテンシャルの場合（例えばreaxFF）はchargeを選択する必要があります。これで必要な結晶構造ファイルができました。直接supercellを作りたい時は一度VESTAでcifを読んで、そちらでsupercellにして、cifでexportすればできそうです。ただ、結晶構造ファイルをベースにlammps側でsupercellを作成することができるので、普通は必要はありません（欠陥とか導入する時は必要になりますが）。
なお、実際に計算を試すにはポテンシャルが必要なので、lammpsのポテンシャルのディレクトリを見て、ポテンシャルが使える原子群を含む構造を選ぶ必要があります（本末転倒してますが）。今回の例では鉱物であるcristobalite (SiOSUB{SIZE(9){2}})の計算を行ってます。SiOSUB{SIZE(9){2}}については、lammpsプログラムに付属のTersoffタイプやVashishtaのポテンシャル（３種）が利用可能です。
**計算のスクリプト（構造最適化）
lammpsの計算では、計算の条件などを書いた入力ファイル（テキストファイル）を用意する必要があります。まずは構造最適化から。適当な構造から始めるとMD計算が不安定になりかねないので、最初に構造最適化をしておく方が無難です。以下にcristobalite構造を4x4x4のsupercellにして、Tersoffタイプのポテンシャルで実行する場合の入力ファイルの例を示します。なお、このポテンシャルは体積を実際よりも大きく評価する傾向があります（それはベースとしたDFT計算の影響が大きいようですが）。
# test Tesoff SiO2 potential for cristobalite
units metal
boundary p p p
atom_style atomic
read_data cristobalite.lmp
replicate 4 4 4
pair_style tersoff
pair_coeff * * SiO.tersoff Si O
thermo_style custom step etotal vol
thermo 10
fix 1 all box/relax aniso 0.0
minimize 0.0 1.0e-12 1000 10000
write_data cristobalite-4x4x4-min.lmp
最初の行はコメントです。
2行目のunitsは使う単位系の指定で、metalって変ですが、これを指定すると普段使いなれている単位系になるので、無機物の場合にはこれを普通使います。なお、使っているポテンシャルの単位系とも関係しているので、ポテンシャルがrealを使っているなら、realに変更する必要があります。これらのコマンドはマニュアル(マニュアルのhtmlファイルがlammpsのDocファイルにあるので、それをウェブブラウザーで開いておくと便利）のコマンドのところで詳しい情報を得ることができます。metalの使う単位系もそこに書いてあります。例えば圧力単位はbar, 長さはÅ、時間はpsなどになります。
3行目は境界条件で、pは周期境界条件で、ここではxyz全ての方向で周期境界条件を付けてます（３次元結晶の場合はこれ）。
4行目は原子スタイル。これも無機結晶構造のMDの場合はatomicにする。ただcharge可変型のポテンシャルの時はchargeにする。
5行目は使う結晶構造ファイルの指定です。この場合はcristobalite.lmpというファイルを先に述べた手順で予めOVITOを使って作成しています。
6行目はsupercellの指定で、上の行で指定した構造の4x4x4のsupercellを使うことになります。3x4x4x4x4なので、これで作られるMDセル内には768個の原子があることになります。
7、8行目はポテンシャルの指定で、今回はTersoffタイプを使っていて、ポテンシャルファイルはSiO.tersoffです。このスクリプトと同じディレクトリに存在する必要があります。lammpsのpotentialsディレクトリにあるので、コピーするかリンクを作っておきます。Si Oは使う原子の指定です。
9行目のthermo_styleは計算時の出力指定で、ここで指定したパラメータが標準出力に出ますが、log.lammpsファイルにも書き込まれます。ここではカスタム指定で、ステップ数、全エネルギー、体積を出力させてます。ここはMDの前段階の構造最適化目的なので、全エネルギー以外はあまり必要ないと思います。様々な出力が可能ですが、何を出力できるかはマニュアルを見る必要があります。
10行目のthermoは出力を何ステップ毎にするかの指定です。ここでは10ステップ毎にthermo_styleで指定したパラメータを出力させます。
11行目はfixで、計算中の制約や計算方法を指定するものです。ここではbox(cell)の体積を異方的に緩和させる指定になってます。なかなか複雑で、まだあまりよく理解できていません。
12行目は構造最適化で使うコマンドminimizeで、収束条件と最大のステップ数を指定しています。これが計算を実施させるコマンドになります。
13行目は計算終了後の構造を出力させるコマンドです。4x4x4のsupercellで構造最適化された構造が出力されます。次のMD計算の出発構造として使いたいので、ここで出力させます。このファイルもOVITOで見ることができます。
以上はミニマムなものですが、今回は構造最適化自体が主目的ではないので、これで十分だと思います。また、上記をテキストファイルとしてコピーするとそのまま使えると思いますが、wiki表示のために先頭にスペースが入る場合はそれらを削除してください。
このスクリプトをin.cristobaliteとしてセーブして（examplesなど見ていると、入力ファイルにin.*とつけるのがlammps流のようです）、ターミナルから以下のようにしてlammpsを実行します。最初はシリアル版で、２つ目は並列版の場合で、この例では２つのコアを指定しています。使っているMacが２コアのためですが、もっとある場合はその数を指定すると計算が速くなります。ただ、古いMac mini (Intel CPU）で使ってみたところ、コア数多くすると逆に遅くなることがありました。一度、コア数を変えて、計算時間をチェックした方がよさそうです。
% lmp_serial < in.cristobalite
% mpirun -np 2 lmp_mpi < in.cristobalite
この計算はすぐ終わります。体積やエネルギーが収束していることを確認します。構造最適化されたcristobalite-4x4x4-min.lmpをOVITOからオープンして眺めます。構造が壊れていたりしたら、ポテンシャルか、初期構造がおかしいのでしょう。なお上の計算ではredirect (<)を使ってますが、-inを使うようにwarningが出ます（並列の場合特に）。Mac, Linuxでredirectを使って特に問題になったことはありませんが、問題になる場合は上の例では、
% mpirun -np 2 lmp_mpi -in in.cristobalite
とすればOKです。
**計算のスクリプト（NVT MD計算）
次はNVTアンサンブルでMD計算をさせてみます。NVTは体積と温度を一定に保つMD計算です。
# test Tesoff SiO2 potential for cristobalite
units metal
boundary p p p
atom_style atomic
read_data cristobalite-4x4x4-min.lmp
#replicate 4 4 4
pair_style tersoff
pair_coeff * * SiO.tersoff Si O
thermo_style custom step etotal temp press
thermo 10
velocity all create 300 345678 dist gaussian
timestep 0.001
dump 1 all atom 50 cristobalite-4x4x4-nvt.lammpstrj
#dump 2 all cfg 100 *.cfg mass type xs ys zs vx vy vz fx fy fz
fix 1 all nvt temp 300.0 300.0 0.1
run 10000
write_data cristobalite-4x4x4-nvt.lmp
構造最適化の最初の方とほぼ同じです。ただ、今回は構造最適化した構造を読み込んでいて、それは既に4x4x4のsupercellになっているので、replicateは不要になります（残しておくとsupercellをさらに4x4x4にするので、コメントアウトしてます）。
thermo_styleは今回はNVTなので、圧力を表示するように指定しています。出力はbar単位です。
ここで出てくる新しいコマンドは、１つ目がvelocityです。MDなので各原子に初期速度を指定する必要があります。構造最適化では速度や温度は必要ありませんでした。その指定を行うのがこのコマンドです。300 Kになるように指定しています。345678としているのは乱数のための種で、整数であれば何でもかまいません。速度分布をガウシアンと指定しています。
その次は時間ステップで、metalの場合はps（ピコ秒）単位なので、0.001は1 fs（フェムト秒）と指定していることになります。普通は1 fsくらいが適当ですが、原子が軽い時（水素など）や超高温時には原子数がおかしくなるエラーを生じます。そういう場合には時間ステップをもっと短くします。
dumpは途中の計算結果をファイルに書き込ませるコマンドで、この例では原子位置を50ステップ毎にtrajectoryファイルに落としています。このファイルは後でOVITOで読んで、アニメーション表示することができます。
その次のコメントアウトしているdumpは別の条件で構造などを書き込むものです。この場合、各指定ステップ毎に*.cfgファイルができるので、沢山のファイルが生成されることになります（ちょっと鬱陶しい）。このファイルもOVITOで一括して読めます（動画的に見ることができる）。速度や力も入っているので、OVITOでそれらを矢印で示すことができますが、粒子数が多いので煩雑すぎて見てもよく分からないものになります。
fixは今回はNVTアンサンブル（体積と温度一定）を指定しています。温度を指定する必要があるので、temp 300.0 300.0 0.1としてます。tempは温度指定、最初の300.0は計算最初の温度、次の300.0は計算最後の温度で、加熱や冷却したい場合は後者の値を変えれば可能になります。最後はダンピングの値で、これは温度制御で使うパラメーター（Nose-Hoover thermostatで使う)で、100ステップ程度が適当とマニュアルには書いてありますので、100 step * 0.001で0.1としています。ここは0.1と書く代わりに$(100.0*dt)としてもOKです。この場合、timestepを変えてもこちら側は自動的に変更されるので、timestepをよく変える場合には便利でしょう。
runが計算を実施させるコマンドで、数値は計算する総ステップ数です。
最後は計算最後の構造を出力するコマンドです。
計算は以下のようにして実施します。上のファイルをin.cristobalite-nvtと保存しているとして。
% lmp_serial < in.cristobalite-nvt
% mpirun -np 2 lmp_mpi < in.cristobalite-nvt
ターミナルに出てくる温度や圧力が収束していることを確認します。構造もOVITOで見てみます。
**計算のスクリプト（NPT MD計算）
最後にNPTアンサンブルでMD計算をさせてみます。NPTは温度と圧力を一定に制御するMD計算なので、実験（現実）と対応させやすいアンサンブルです。
# test Tersoff SiO2 potential for cristobalite
units metal
boundary p p p
atom_style atomic
read_data cristobalite-4x4x4-nvt.lmp
pair_style tersoff
pair_coeff * * SiO.tersoff Si O
thermo_style custom step etotal temp density press lx ly lz
thermo 10
#velocity all create 300 456789 dist gaussian
timestep 0.001
dump 1 all atom 20 cristobalite-4x4x4-npt.lammpstrj
#dump 2 all cfg 10 *.cfg mass type xs ys zs vx vy vz fx fy fz
fix 1 all npt temp 300.0 300.0 $(100.0*dt) aniso 1 1 $(1000.0*dt)
run 20000
write_data cristobalite-4x4x4-npt.lmp
NVT計算とほぼ同じです。最初に前のNVT計算結果の構造を読み込んでいます。
thermo_styleは今回はNPTなので、温度以外に密度や格子定数の一部(lx,ly,lz)を表示するように指定しています。
前回のNVTはMD計算なので、その出力には原子毎の速度も含まれています。そのため、velocityコマンドは今回は不要です。そのためコメントアウトしています。
fixは今回はNPTアンサンブル（圧力と温度一定）を指定しています。NPTなので、温度に加えて、圧力を指定する必要があります。圧力も温度と似ていますが、いくつか制御の仕方が選択できて、iso, aniso, tri, 直接応力指定が可能です。isoは等方的に、anisoはx,y,zは独立に、triは軸間が90度以外になることも許容します。さらに応力で直接指定すると剪断応力をかけたりできるようです（変形や粘性計算）。今回は正方晶系で単斜や三斜に歪むことも予想されないので、anisoを選択しています。立方晶系の結晶やガラス・液体の場合はisoの指定でいいでしょう。圧力の指定は最初と最後を1 barとしてます。温度同様に最後のパラメーターはダンピングで、マニュアルによると1000 stepくらいが適当となっているので、それを今回は数式で指定しています。この例では20000ステップの計算してます。
計算は以下のようにして実施します。in.cristobalite-nptとして保存されているとして、
% lmp_serial < in.cristobalite-npt
% mpirun -np 2 lmp_mpi < in.cristobalite-npt
この計算は前よりは少し時間がかかります。温度や圧力が制御されていることを確認します。log.lammpsというファイルにthermo_styleで指定した内容が保存されています。もし、保存したい場合は名前を変えておきます（次の計算で上書きされるため）。各パラメータの時間変化をプロットして見ることに使えます。
今の計算はテスト計算ですが、実際の計算ではもっとステップ数を長くしたり、平衡を確認するために繰り返したり、supercellをもう少し大きくしたり、温度や圧力を変えたりします。なお、今回各計算を別々にやってますが、１つのファイルで複数の計算をまとめることもできます。pythonでスクリプトを作れば、条件を変えた一連の計算なども自動でできるはずです。
msd(平均二乗変位)は以下のようにすれば出力させることができます(これ以外にもOVITO側で処理することもできます）。上記のfix前後を以下で置き換えます。
compute 1 Si msd
compute 2 O msd
fix 1 all npt temp 300.0 300.0 $(100.0*dt) aniso 1 1 $(1000.0*dt)
fix 2 all ave/time 1 1 5 c_1[4] file Si.msd
fix 3 all ave/time 1 1 5 c_2[4] file O.msd
そうすると計算実施後にSi.msd, O.msdのファイルが出来ているはずです。結晶の場合は最初少し上がって、一定値に落ち着きます。融けている場合や拡散が起きている場合は、時間とともに単調に増加するはずです。その傾きから拡散係数を求めることができます。
**OVITOによるMD trajectoryファイルの表示
MD計算で作られた*.lammpstrjをOVITOで読み込みます。中央下部にあるPlayボタンをクリックすると原子を時間変化アニメーションとして見ることができます。最初原子だけが表示されます。結合を棒で表示したい場合は、右側上のAdd modification...からCreate bondsを選択します。Cutoff distanceを適当に調整することでSi-O間の結合を示すことができます。それ以外にも原子の大きさを変えたり色々できます。多面体表示するには、まずAdd modification...でSelect Typeをクリック。そこで多面体の中心原子、今の場合Siを選びます。次にCreate bonds、その次にCoordination polyhedraを選びます。これで多面体表示になります。動径分布関数の表示などもできます（これはlammps側で計算させることもできますが、OVITO側でも計算できます）。軌跡を表示することも可能です。色々な機能がありますが、まだまだ使いこなせていません。なお、(Pro)と書いてある機能はOVITO Proのライセンスが必要です。
右側のカメラタブをクリックすると動画を出力できます。Complete animationのラジオボタンをクリックします。Save to fileをチェックして、Choose...でフォーマットを選んで、ファイル名を指定して、下の方のRender active viewportで動画ファイルの作成を開始します。そうやって作ったcristobaliteのNPT計算のgifファイルを下に示しています（これは2000 stepだけでbondを表示）。赤がSiで青が酸素。
#image(https://mkanzaki.sakura.ne.jp/images/cristobalite_lammps.gif,center,60%)
**格子常数
いまだdumpで格子常数をうまく出力する方法が分かりません。lx, ly, lzで一応各軸の長さは出せるのですが、角度の出し方が分かりません（cristobaliteの例では必要ありませんが）。もっともdump *.cfgで出力すると格子のマトリックスが書いてあるので、そこから計算することができそうでした。ただ、各出力ステップ毎でファイルが出力されるので処理が面倒そうです。trajectoryファイルにも格子常数が入っているのですが、フォーマットがよく分かりません。２つ目の値は格子の長さに対応しているようですが。
OVITOで見えている構造を出力して、VESTAで眺めようとxyz, POSCARフォーマットで出力したのですがうまくいきません。いくつか試して、FHI-aimsで出力するとVESTAで読めることが分かりましたが、拡張子によっては読んでくれません。*.inはOKなので、読めない時は.inにrenameしてみてください。POSCARで読めないのも拡張子の問題なのかも。
**ポテンシャル
現状、lammpsでは金属系ポテンシャル（EAM)が多く、無機化合物のポテンシャルは多くないようです。[[OpenKIM:https://openkim.org/]]で探すともう少し出てきます。OpenKIMのポテンシャルをlammpsで使うことは可能です。いくつかのポテンシャルでは追加のライブラリを入れて、lammpsをコンパイルし直す必要があります（実行版をダウンロードしている場合にはこれができませんが）。さらにNISTもポテンシャルのデータベースを作っているそうです（2023年のWorkshopでNISTの発表があった）。別のページで紹介しているGULPの方が使える無機物のポテンシャルが多いと思います。ただ、GULPは元々MDとして開発された訳ではなく、作者もMDは得意でないとされてます。ポテンシャルを使えるように変換すればいいのでしょうが、まだlammps用のポテンシャルのフォーマットを十分理解できていません。
**kimを使う
実行版ではkimパッケージは最初から入っているようで、特に何もしなくても使えました（MacBook Airで）。kimを使う場合はopenKIMにアクセスして、必要なポテンシャルを探して（周期律表で元素をクリックする）、そしてkimIDを見つける必要があります。例えばSiでStillinger-Weberのポテンシャルは探すと、"SW_StillingerWeber_1985_Si__MO_405512056662_006"がありました。右上の方のfileをクリックすると、ポテンシャルのファイルをダウンロードできます。実際には既にプログラム中に内臓されているので、このStillinger-Weberポテンシャルをダウンロードする必要はありません。
使い方はexamples/kimの入力ファイルを見れば分かりますが、普通pair_style, pair_coeffで指定するポテンシャルの代わりに、kim initとkim interactionsの２つのコマンドで使うポテンシャルを指定するだけです。kim initでkimIDを指定して、kim interactionsで元素を指定するようです。
ちなみにexamples/kimにあるin.kim-pm.meltはこのStilinger-Weberポテンシャルを使った例ですが、中身は"SW_StillingerWeber_1985_Si__MO_405512056662_005"とバージョンが１つ前のものになっているために、そのまま実行するとエラーが出ました。openKIMで探すとバージョンが古いことが分かって、最後の5を6にすると正常に計算できました。同様な対応がin.kim-sm.meltでも必要でした。
なお、実行版に付属しているkim用のファイル(*.soになっている）は400あまりでした（これ自体は2021に作られていた）。そこにない場合はOpenKIMからダウンロードしてくるしかないようです。ダウンロードした場合、ポテンシャルによってはkimを使わないでそのまま使えるものもあって（SiOSUB{SIZE(9){2}}のtersoffポテンシャルなど）、通常のpair_style, pair_coeffで指定することで使えました。
**NISTのポテンシャルリポジトリ
[[NIST:https://www.ctcms.nist.gov/potentials/]]もポテンシャルのリポジトリーを作ってます。こちらも周期律表の元素をクリックすることで必要な元素のポテンシャルを探すことができます。Si-OだとopenKIMとほぼ同じですが、kimにはないものもありました。こちらの方は各ポテンシャルの説明のところに必要なファイルのリンクがあるので、ダウンロード可能です。
**reaxFFを使う
実行版ではreaxFFは有効になっているので、すぐ使えます。examples/reaxffの中にHSUB{SIZE(9){2}}O, ZnOH, Fe(OH)SUB{SIZE(9){3}}などの例があります。「分子動力学法と原子間ポテンシャル」（渡邉孝信著・森北出版）にはその中のHSUB{SIZE(9){2}}Oを実行する例があります。今回試したのはCa-Si-O-H系のもので、これはNISTのリポジトリーの方にあるので、そこからダウンロードしました(ffield-CaSiOH-lammps.reax.txt)。このポテンシャルはセメント材料の水和プロセス解明などに使われています。ここでは含水カルシウムケイ酸塩のpoldervaartite (CaSUB{SIZE(9){2}}SiOSUB{SIZE(9){3}}(OH)SUB{SIZE(9){2}})構造でまずは試してみました。結晶構造はCODの9011378を使いました（水素位置が決まっているものが必要です）。DLしたcifをOVITOで読んで、lammps形式で出力します。ここで１つ注意は、reaxFFは以下でも見るようにatom_styleがchargeであるので、lammps形式で出力する際にデフォルトのatomicではなく、chargeを選択する必要があります（そうしないと実行時にエラーが出るし、逆にatomicの場合にchargeを選択してもエラーになる）。それをpoldervaartite.lmpとして保存しています。前と同様にまずは構造の最適化(minimize)から。reaxFFは時間かかるので、今回は比較的小さい2x2x2のスーパーセルで実施しています。
# test reaxxff potential for poldervaartite
units metal
boundary p p p
atom_style charge
read_data poldervaartite.lmp
replicate 2 2 2
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * ffield-CaSiOH-lammps.reax.txt Ca Si O H
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes
thermo_style custom step etotal vol
thermo 10
fix 1 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff maxiter 400
minimize 0.0 1.0e-10 1000 10000
write_data poldervaartite-2x2x2-min.lmp
前と違うところはatom_styleがchargeになっていることで、reaxFFでは各原子のチャージは固定ではなく、その環境に応じて変化します。pair_style, pair_coeffも前とちょっと違いがあります。neighbor, neigh_modifyはデフォルトでもいいようですが、examples/reaxffの例から適当に設定してます（問題ある可能性がありますが）。fixもちょっと違いがあります。qeq/reaxffが必要です。それ以外は前のcristobaliteの例と同じです。この計算はすぐ終わります。次にNVTアンサンブルの計算。前の計算で出力されたpoldervaartite-2x2x2-min.lmpを使ってます。なお、この系では水素が含まれています。軽い元素は速度が速いために、cristobaliteの時と同じtimestepで計算すると問題が起こります。そのために0.0002 (ps)にしています。
# test reaxff CaSiOH potential for poldervaartite
units metal
boundary p p p
atom_style charge
read_data poldervaartite-2x2x2-min.lmp
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * ffield-CaSiOH-lammps.reax.txt Ca Si O H
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
thermo_style custom step etotal temp press
thermo 10
velocity all create 300 345678 dist gaussian
timestep 0.0002
dump 1 all atom 50 poldervaartite-2x2x2-nvt.lammpstrj
fix 1 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff maxiter 400
fix 2 all nvt temp 300.0 300.0 1
run 1000
write_data poldervaartite-2x2x2-nvt.lmp
次にNPTアンサンブルの計算へ。
# test reaxff CaSiOH potential for poldervaartite
units metal
boundary p p p
atom_style charge
read_data poldervaartite-2x2x2-nvt.lmp
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * ffield-CaSiOH-lammps.reax.txt Ca Si O H
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
thermo_style custom step etotal temp density press lx ly lz
thermo 10
timestep 0.0002
#compute 1 Si msd
dump 1 all atom 50 poldervaartite-2x2x2-300K-1bar.lammpstrj
#dump 1 all cfg 100 *.cfg mass type xs ys zs vx vy vz fx fy fz
fix 1 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff maxiter 400
fix 2 all npt temp 300.0 300.0 $(100.0*dt) aniso 1 1 $(1000.0*dt)
run 10000
write_data poldervaartite-2x2x2-300K-1bar.lmp
これはNVTからNPTにして、圧力指定(1 bar)が追加されているだけです。これらを実行してみて構造が保たれていて問題なさそうでした（下の図）。なお、この系ではtimestepを短くしているので、平衡になるためのステップ数がより多く必要になります。
#image(https://mkanzaki.sakura.ne.jp/images/poldervaartite-300K-1bar.png,center,33%)
次に同じReaxFFポテンシャルを使って、SiOSUB{SIZE(9){2}}-HSUB{SIZE(9){2}}Oメルトを作ってみた。5000 KでNPT計算。1200個の原子で計算。OHとHSUB{SIZE(9){2}}O両方が存在。まあ大体いいと思うが、いくらか変な欠陥があったり、OSUB{SIZE(9){2}}分子があったりする。まだ最初の構造作成とその後のアニーリングに課題があるかも。原子数ももう少し増やしたいところ。下の図は一度OVITOから出力して、それをVESTAで読んで作ってます。
ちなみにSiOSUB{SIZE(9){2}}-HSUB{SIZE(9){2}}O系メルトの初期構造を今回どう作ったかというと、実は別のMDプログラムで河村雄行先生のmxdシステムがあるのですが、それの入力ファイル作成用プログラムのmxdinputを使いました。そこで出力されるfile07.datをVESTAが読めるので、それをcifとして出力して、後はOVITOで同様に処理しました。昔file07.datがVESTAの前身のソフトで読めるように泉富士夫先生に少し協力しました。
その後、もう少し初期構造を自分でカスタマイズしたいので、pythonでCaSiOH系メルト（フルイド？）の初期構造を作るコードを自作しました（原子数と密度を指定するようにした）。基本的には乱数で原子を配置していくのですが、既に置いている原子と接近しすぎる場合は乱数で別の位置を選びます。pythonなので結構時間かかります。ファイルをxtl形式にしたので（簡単なので）、やはりVESTAで読んで、cifで出力します。もっともlammpsくらいになると、そういうコードは誰かが既に作っていると思いますが。
#image(https://mkanzaki.sakura.ne.jp/images/SiO2-H2O-5000K.png,center,50%)
**Linux版ソースからインストール（2023/12/26)
既にコンパイルされている実行版は便利なのですが、色々と試そうとすると別途モジュールやらライブラリーのインストールが必要な場合があり、実行版では対応できなくなります。そこでLinux (Fedora39)でソースからインストールしてみました。最近出版された「分子動力学法と原子間ポテンシャル」（渡邉孝信著・森北出版）にもlammpsのインストール法が紹介されていました。私の環境ではQuantum Espressoを導入している関係で、gcc (gfortran), openmpiは既にソースからインストール済みなので、特に追加で何かする必要はありませんでした。
ソース（tarボール）は同じく[[lammpsページ:https://www.lammps.org/]]からダウンロードできます。最新版は2Aug2023となってました。これを解凍して、srcディレクトリでmake psまたはmake package-statusを実行します。追加パッケージの情報が出てきます。実行版はこれらを適当に選択してコンパイルされているようですが、全てが入っている訳ではないので、入ってないパッケージは使えません。そういう場合にソースからのコンパイルが必要となります。なお、パッケージによっては別途事前にインストールが必要なものもあります。
デフォルトでは全てNOになっているので、必要なパッケージをyesにしてから、makeする必要があります。make時間はそれほどかからないので、必要になった時に変更して、makeし直せばいいと思います。可能なオプションはmakeとだけタイプすると出てきます。Tersoffポテンシャルを使う時は、これは多体ポテンシャル(MANYBODY)なので、make yes-MANYBODYとします。ReaxFFポテンシャルを使う場合は、make yes-REAXFF。なお、make yes-mostとすると、大半のモジュールをyesにしてくれます。再度make psすると何がyesになっているか分かります。make yes-basicもあります。なお、実行版に何が入っているかはlmp_mpi -hで分かります（出力の最初の方）。
パッケージを選択したら、コンパイルします。私はopenmpiをインストールしているので、mpi版（並列版）を作ります。これには以下を実行します。特に問題なければsrcディレクトリにlmp_mpi実行ファイルができているはずです。
make mpi
シリアル版を作る場合は以下を実行します。lmp_serialができているはずです。
make serial
.bashrcにsrcへのパスを書き込んでおきます。lammpsの実行は上記と同じで、
mpirun -np 8 lmp_mpi < in.cristobalite-npt
**Macの実行版
M3のMacBookAirを買ったので、lammpsをAnacondaでupdateしたのですが、Intel版のままでした。rosetta2を使ったエミュレーションで走らせることになります。Apple Silicon用のソースからコンパイルするのは結構大変そうなので、今のところそのまま。普段、計算にはLinux PC使っているので必要性はあまりないのです。古いIntel CPUのMacBookAirと計算速度を比べてみました。計算にはreaxffポテンシャルを使ってます。
-MacBookAir Intel core i5 1.1 GHz CPU 2core使用：184 s
-MacBookAir M3 CPU 2core使用：119 s
-MacBookAir M3 CPU 4core使用：74 s
-Linux PC Intel i7 13700 8core使用：40 s
M3で走らせても多少早くなってました。
OVITOはApple Silicon対応版があったのでインストール。VESTAもuniversal binary版が最近出てたので、それもインストールしました。ただし、これはdevelopment versionで、windowサイズ変えると表示がおかしくなることがありました。
**注意
最初、石英を扱おうとしたのですが、途中でエラーが出ました。軸間の角度が120度以上のものは計算開始時にエラーが出ることがあります（近接原子のbook keepingの関係らしい）。対処法としてはsupercellを使うといいようです。マニュアルに記述があります。1x1x2のsupercellをVESTAで作って、cifで出力して（格子外の原子は除いておく）、OVITOで*.lmpへ変換したら、エラーが出なくなりました。または軸を取り替えて120度以下にすることも有効です（VESTAでやっておく）。
単斜晶系等でbeta角度が120度近いと計算中にセルが変換されることがあります。計算自体は全く問題ないのですが、OVITOでムービー見ている時に突然セルの取り方が変わることがありました。
atom_styleがatomicの場合は、先に述べたようにOVITOから出力するファイルが初期構造として使えるのですが、chargeの場合はエラーがでます。これはどうもファイルの座標データのところにchargeが入ってないためのようでした。対応として、OVITOから出力する時にatomicではなくchargeを選べば自動的に0.0を挿入してくれます。例えばReaxFFの時にはatom_styleがchargeとなり、tersoffはatomicです。
unitsがrealとmetalでは単位が色々と異なります。特にtimestepはrealがfs、metalがps単位で３桁違うので注意が必要です。
水素のような軽い元素は移動が速いので（温度が非常に高い時も）、timestepは0.001(unitsがmetalの時）だとエラーが出るので、短くする必要があります。
ポテンシャルを指定するところ(pair_coeff)に書き込む元素の順番は、使う構造ファイル中で出てくる順番と合わせる必要があるようです。間違えてもエラーはでないのですが、結果がおかしくなります。
COD等で拾ってくる構造には固溶体、原子が部分的にしか席を占有してないケースや位置がいくつかの等価な席にdisorder分布している場合などがあります。そういう場合はcifをそのまま使えません。それをベースにモデルを作る必要があります。
最近のIntelのCPUではパーフォーマンスコアとエフィシェンシーコアの両方が搭載されています。lammpsの計算には当然パーフォーマンスコアの方が速いので、指定する並列数はパーフォーマンスコアの数にした方がよいようです。パーフォーマンスコア数以上に設定した場合、むしろ計算が遅くなる場合があります。ただし、実際にコア数を変えて、計算時間がどうなるかチェックした方がいいと思います。
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")