LAMMPSを使ってNPD分子100分子系のMD計算を行う
今更ではありますが、LAMMPSを使って有機分子集合体の状態を計算してみます。第一原理計算を用いて複数の分子を計算する場合、かなり時間がかかるのに対して、
分子動力学計算(MD計算)は普通のノートPCなどでも数時間で計算ができます。
もちろん精度の問題はありますが、セットアップなどの障壁がなくなれば多くの人が手軽にMD計算できるようになると思ったので、この記事でまとめておきます。
具体的には、本記事では、
Packmol + AmberTools + Moltemplate + LAMMPS を用いて有機分子（NPD）100分子系のMD計算を行います。そして最終的に、NPD 100 分子からなる系を作成し、エネルギー最小化 → NVT → NPT計算まで実行することを目標にします。
まずはMD計算を行うための環境を整えていきましょう。
Windowsだと全てのソフトのセットアップができなかったので、Linux推奨です。今回はWSL(Ubuntu)を使いました。WSLの導入は別の記事などを参考にしてください。
1. ソフトウェア環境の構築（Miniforge + LAMMPS周辺ツール）
Miniforge 
Ubuntu 環境で conda ベースの仮想環境を作るために Miniforge を使用します。Miniforgeが良い理由は以下の記事に書かれています。
以下のコマンドで Ubuntu (WSL) に Miniforge をインストールします。Ubuntuを開いて、以下のコードを打ってください。
wget https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-x86_64.sh
bash Miniforge3-Linux-x86_64.sh
途中でライセンスや初期化の確認が出てきますが、基本的にyesで進めてください。
インストールが終わったら、bashにcondaを有効化させてください。
~/miniforge3/bin/conda init bash
source ~/.bashrc
次にMD計算を行う仮想環境を作りましょう。
conda create -n lammps_md python=3.10
conda activate lammps_md
以下のように、プロンプトが(lammps_md)となっていればOKです。
(lammps_md) XXXX:~$
以降のインストール作業はすべてこの環境内で行ってください。
LAMMPS 
LAMMPSはオープンソースの古典分子動力学(MD)ソルバの一つです。様々な力場（相互作用）が使用可能で、材料シミュレーションに適しています。計算速度が遅いなどの短所があるみたいですが、今回扱う系くらいであれば問題なく処理できます。
様々なインストール方法がありますが、Miniforgeを入れたのでcondaを使ってインストールしましょう。
conda install -c conda-forge lammps
lmp -h
lmp -hと打ってヘルプが表示されれば、インストール完了です。
Large-scale Atomic/Molecular Massively Parallel Simulator - 29 Aug 2024
Usage example: lmp -var t 300 -echo screen -in in.alloy
List of command line options supported by this LAMMPS executable:
-echo none/screen/log/both : echoing of input script (-e)
-help : print this help message (-h)
-in none/filename : read input from file or stdin (default) (-i)
-kokkos on/off ... : turn KOKKOS mode on or off (-k)
-log none/filename : where to send log output (-l)
-mdi '<mdi flags>' : pass flags to the MolSSI Driver Interface
-mpicolor color : which exe in a multi-exe mpirun cmd (-m)
-cite : select citation reminder style (-c)
-nocite : disable citation reminder (-nc)
-nonbuf : disable screen/logfile buffering (-nb)
-package style ... : invoke package command (-pk)
-partition size1 size2 ... : assign partition sizes (-p)
-plog basename : basename for partition logs (-pl)
-pscreen basename : basename for partition screens (-ps)
-restart2data rfile dfile ... : convert restart to data file (-r2data)
-restart2dump rfile dgroup dstyle dfile ...
: convert restart to dump file (-r2dump)
-restart2info rfile : print info about restart rfile (-r2info)
-reorder topology-specs : processor reordering (-r)
-screen none/filename : where to send screen output (-sc)
-skiprun : skip loops in run and minimize (-sr)
-suffix gpu/intel/kk/opt/omp: style suffix to apply (-sf)
Packmol 
初期配座を作るために、複数分子をランダムに配置させるためのフリーソフトウェアです。
こちらもcondaを使って、以下のコードでインストールしましょう。
conda install -c conda-forge packmol
packmol < /dev/null
エラーメッセージと使用方法が表示されれば、正常に動いています。
################################################################################
PACKMOL - Packing optimization for the automated generation of
starting configurations for molecular dynamics simulations.
Version 21.1.3
################################################################################
Packmol must be run with: packmol < inputfile.inp
Userguide at: http://m3g.iqm.unicamp.br/packmol
Reading input file... (Control-C aborts)
WARNING: File type not (correctly?) specified, using PDB
Seed for random number generator: 1234567
ERROR: Output file not (correctly?) specified.
STOP 171
Ambertools 
力場パラメータ生成のために使用します。力場とは、原子間力を表現する関数形・パラメータの組み合わせのことです。簡単に言うと、原子と原子の間に働く相互作用のことです。
実際に使うのは、Ambertoolsの中のantechamberになります。こちらがWindowsではインストール・ビルドできないため、Linux(今回はWSL)を使用しないといけませんでした（できる方法もあるらしいが、かなり面倒）。公式でもWSLを推奨しています。
conda install -c conda-forge ambertools
antechamber -h
ヘルプが表示されればOKです。
Welcome to antechamber 22.0: molecular input file processor.
Usage: antechamber -i input file name
-fi input file format
-o output file name
-fo output file format
-c charge method
-cf charge file name
-nc net molecular charge (int)
-a additional file name
-fa additional file format
-ao additional file operation
crd : only read in coordinate
crg : only read in charge
radius: only read in radius
name : only read in atom name
type : only read in atom type
bond : only read in bond type
-m multiplicity (2S+1), default is 1
-rn residue name, overrides input file, default is MOL
-rf residue toplogy file name in prep input file,
default is molecule.res
-ch check file name for gaussian, default is 'molecule'
-ek mopac or sqm keyword, inside quotes; overwrites previous ones
-gk gaussian job keyword, inside quotes, is ignored when both -gopt and -gsp are used
-gopt gaussian job keyword for optimization, inside quotes
-gsp gaussian job keyword for single point calculation, inside quotes
-gm gaussian memory keyword, inside quotes, such as "%mem=1000MB"
-gn gaussian number of processors keyword, inside quotes, such as "%nproc=8"
-gdsk gaussian maximum disk usage keyword, inside quotes, such as "%maxdisk=50GB"
-gv add keyword to generate gesp file (for Gaussian 09 only)
1 : yes
0 : no, the default
-ge gaussian esp file generated by iop(6/50=1), default is g09.gesp
-tor torsional angle list, inside a pair of quotes, such as "1-2-3-4:0,5-6-7-8"
':1' or ':0' indicates the torsional angle is frozen or not
-df am1-bcc precharge flag, 2 - use sqm(default); 0 - use mopac
-at atom type
gaff : the default
gaff2: for gaff2 (beta-version)
amber: for PARM94/99/99SB
bcc : bcc
sybyl: sybyl
-du fix duplicate atom names: yes(y)[default] or no(n)
-bk component/block Id, for ccif
-an adjust atom names: yes(y) or no(n)
the default is 'y' for 'mol2' and 'ac' and 'n' for the other formats
-j atom type and bond type prediction index, default is 4
0 : no assignment
1 : atom type
2 : full bond types
3 : part bond types
4 : atom and full bond type
5 : atom and part bond type
-s status information: 0(brief), 1(default) or 2(verbose)
-eq equalizing atomic charge, default is 1 for '-c resp' and '-c bcc' and 0 for the other charge methods
0 : no use
1 : by atomic paths
2 : by atomic paths and structural information, i.e. E/Z configurations
-pf remove intermediate files: yes(y) or no(n)[default]
-pl maximum path length to determin equivalence of atomic charges for resp and bcc,
the smaller the value, the faster the algorithm, default is -1 (use full length),
set this parameter to 10 to 30 if your molecule is big (# atoms >= 100)
-seq atomic sequence order changable: yes(y)[default] or no(n)
-dr acdoctor mode: yes(y)[default] or no(n)
-i -o -fi and -fo must appear; others are optional
Use 'antechamber -L' to list the supported file formats and charge methods
Moltemplate 
LAMMPSで計算するインプットファイル生成に使います。
このソフトもcondaでインストールしましょう。
conda install -c conda-forge moltemplate
moltemplate.sh -h
moltemplate.sh v2.22.5 2025-11-16
lttree_check.py v0.81.2 2021-5-24
########################################################
## WARNING: atom_style unspecified ##
## --> "Data Atoms" column data has an unknown format ##
## Assuming atom_style = "full" ##
########################################################
Error: unable to open file
"-h"
for reading.
moltemplateのversionが表記され、正しくインストールされていればOKです。
Openbabel 
分子形式変換や2D→3D変換に使うソフトです。
conda install -c conda-forge openbabel
以上で、必要なソフト群をインストールできました。これにて環境構築は終了になります
2. LAMMPS用入力ファイル作成
以下の流れで、LAMMPS用のインプットファイルを作っていきましょう。改めて下記サイトから引用させていただきました。
流れとして、
- 分子描写ソフトで分子を描く
-
openbabelで2次元構造を3次元構造に変換（*_3d.molやpdbの作成） -
antechamber(Ambertoolsに入っている)でmolをmol2に変換（原子タイプ・電荷付与） - pdbファイルを使って
packmolにて100分子を格子中にランダム配置（初期配座） - mol2ファイルに対して、
mol22ltを使ってltに変換、力場を設定する -
moltemplateを使って、lammps用のインプットであるsystem.data、system.in.init、設定ファイルであるsystem.in.settingsを作成 - LAMMPS計算
となります。まず、作業用ディレクトリを作成しましょう。今回は ~/lammps/md_npdsim を作業ディレクトリとします。
mkdir -p ~/lammps/md_npdsim
cd ~/lammps/md_npdsim
pwd
以降の操作はこのディレクトリ内で行います。
分子描写
今回は有機ELで正孔輸送性材料としてよく使われているNPD分子を扱っていきます。
SMILESはC=1(C=2C(C=CC1)=CC=CC2)N(C=3C=CC=CC3)C=4C=CC(C=5C=CC(N(C=6C=7C(C=CC6)=CC=CC7)C=8C=CC=CC8)=CC5)=CC4となります。
こちらのサイトを参考にして、ChemSketchを使いましょう。どうやら水素原子も描写する必要があるみたいです。
ちなみに、私はJChemPaintというフリーソフトを使っていますが、このソフトで描写した構造をmol fileで保存して開きなおすと、ChemSketchが起動します。何故？
こんな感じで、構造を描写したら、ACD/Labsをクリックして、3D Viewerを起動します。起動しただけでは何も表示されていないので、左下に表示されるリストから 1-ChemSketch を選択して、一度 ChemSketch に戻り作成した分子構造を表示させます。
次に左下のリストに表示される 2-Copy to 3D を選択すると、 ChemSketch で作成した分子構造が 3D Viewer に表示されます。
そして、Tools メニューから 3Doptimization を実行して、省略されていた水素原子が付加され、同時に分子構造を三次元化することができます。
3D構造ができたら、左下の2-Copy to ChemSk → Structure → Clear current page を選んで、ChemSketch に戻ります。
すると、先ほどまではなかった水素原子が付加されていることが分かります。save asでmol形式として保存してください。
mol→mol2変換
そして、先のディレクトリ（~/lammps/md_npdsim）にnpd.molを格納し、openbabelを使って2D構造のmolファイルを3D構造に変更します。通常、2D描写ツールで書いた構造にはz軸の座標が存在しないためです。下記コードを実行しましょう。
obabel -imol npd.mol -omol -O npd_3d.mol --gen3d
head -n 20 npd_3d.mol
1 molecule converted
OpenBabel02082622013D
78 85 0 0 0 0 0 0 0 0999 V2000
1.1334 -0.3702 -1.8360 C 0 0 0 0 0 0 0 0 0 0 0 0
0.8494 0.8411 -1.2142 C 0 0 0 0 0 0 0 0 0 0 0 0
-0.4364 1.0939 -0.7191 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.4727 0.1406 -0.8578 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.1663 -1.0718 -1.5016 C 0 0 0 0 0 0 0 0 0 0 0 0
0.1277 -1.3258 -1.9789 C 0 0 0 0 0 0 0 0 0 0 0 0
-2.7967 0.4214 -0.3301 N 0 0 0 0 0 0 0 0 0 0 0 0
-3.3531 1.7557 -0.4279 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.0012 2.6195 -1.4712 C 0 0 0 0 0 0 0 0 0 0 0 0
-3.4868 3.9286 -1.5176 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.4066 4.4179 -0.5702 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.8056 3.5312 0.4377 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.2703 2.2372 0.5213 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.9010 5.8124 -0.6510 C 0 0 0 0 0 0 0 0 0 0 0 0
-4.9750 6.4609 -1.8932 C 0 0 0 0 0 0 0 0 0 0 0 0
-5.3059 7.8054 -1.9999 C 0 0 0 0 0 0 0 0 0 0 0 0
xyzの3次元座標に変換されていることが分かります。
次にmol→mol2へ変更し、原子タイプ・電荷の付与を行います。LAMMPS で有機分子を扱うためには、力場に対応した原子タイプと部分電荷が必要になります。
ここでは力場として GAFF(General AMBER Force Field)を使用します。
GAFF は有機分子全般を幅広くカバーする汎用力場で、AmberTools の antechamber を使うことで、分子構造から原子タイプと部分電荷を自動的に割り当てることができます。
今回は「まず MD を回せる系を作る」ことを目的として、手軽に使える GAFF + Gasteiger 電荷の組み合わせを採用します。
antechamber -fi mdl -i npd_3d.mol -fo mol2 -o npd.mol2 -at gaff -c gas -rn npd -dr yes
コマンドの説明は以下の通りです。
-
-fi: 入力ファイル形式の指定、ここではmol形式なので、-fi mdl、一方、Gaussianの出力を使う場合は-fi goutとなる -
-i: 入力ファイル名の指定 -
-fo: 出力ファイル形式の指定、ここではmol2形式として出力 -
-o: 出力ファイル名の指定 -
at: 原子タイプの指定、ここではGAFFを指定、gaff2という新しいものもあるみたい -
-c: 電荷を付与、ここでは、簡単なGasteigerを指定、他には、AM1-BCC(半経験的＋経験的補正)のbccや量子化学計算のESPからフィットするrespなども選べる -
-rn: 残基名(rename)の指定、mol2やAmber/LAMMPS変換時に使われる識別子 -
-dr: 構造チェックをするかどうか（変な構造になっていないかのチェック)
そして、中身を見てみましょう。
head -n 30 npd.mol2
@<TRIPOS>ATOM
1 C1 1.1330 -0.3700 -1.8360 ca 1 npd -0.062132
2 C2 0.8490 0.8410 -1.2140 ca 1 npd -0.060383
3 C3 -0.4360 1.0940 -0.7190 ca 1 npd -0.040234
4 C4 -1.4730 0.1410 -0.8580 ca 1 npd 0.037828
5 C5 -1.1660 -1.0720 -1.5020 ca 1 npd -0.040234
6 C6 0.1280 -1.3260 -1.9790 ca 1 npd -0.060383
7 N1 -2.7970 0.4210 -0.3300 nh 1 npd -0.269176
8 C7 -3.3530 1.7560 -0.4280 ca 1 npd 0.037871
9 C8 -3.0010 2.6190 -1.4710 ca 1 npd -0.039617
10 C9 -3.4870 3.9290 -1.5180 ca 1 npd -0.052498
11 C10 -4.4070 4.4180 -0.5700 cp 1 npd -0.018265
12 C11 -4.8060 3.5310 0.4380 ca 1 npd -0.052498
13 C12 -4.2700 2.2370 0.5210 ca 1 npd -0.039617
14 C13 -4.9010 5.8120 -0.6510 cp 1 npd -0.018265
15 C14 -4.9750 6.4610 -1.8930 ca 1 npd -0.052498
16 C15 -5.3060 7.8050 -2.0000 ca 1 npd -0.039617
17 C16 -5.6470 8.5680 -0.8760 ca 1 npd 0.037871
18 C17 -5.6440 7.9220 0.3710 ca 1 npd -0.039617
19 C18 -5.2680 6.5660 0.4800 ca 1 npd -0.052498
20 N2 -5.9440 9.9770 -1.0280 nh 1 npd -0.269176
21 C19 -6.6800 10.4450 -2.1920 ca 1 npd 0.037828
22 C20 -7.6240 9.6310 -2.8490 ca 1 npd -0.040234
このように、各原子に GAFF の原子タイプ（ca, cp, nh など）が割り当てられ、さらに各原子の部分電荷（行の末尾の小数値）も付与されていることが分かります。これにより、LAMMPS で力場パラメータを割り当てるための準備が整いました。
packmolにてbox.pdbを作る
packmolの入力用に、以下コマンドでmolファイルをPDB に変換します。
obabel -imol npd_3d.mol -opdb -O npd.pdb
head -n 5 npd.pdb
COMPND UNNAMED
AUTHOR GENERATED BY OPEN BABEL 3.1.0
HETATM 1 C UNL 1 1.133 -0.370 -1.836 1.00 0.00 C
HETATM 2 C UNL 1 0.849 0.841 -1.214 1.00 0.00 C
HETATM 3 C UNL 1 -0.436 1.094 -0.719 1.00 0.00 C
そして、packmol用のinputファイルを作ります。ファイル名は、packmol_100.inpとします。コマンドで作る場合は、以下の通り。
cat > packmol_100.inp << 'EOF'
tolerance 2.0
add_box_sides 2.0
filetype pdb
output box.pdb
structure npd.pdb
number 100
inside box 0. 0. 0. 60. 60. 60.
end structure
EOF
中身の説明として、
-
tolerance: 分子間最小距離の目安（Å）、今回は2Å -
add_box_sides: 配置領域に対して余白を足す（packmolが出力座標の端を少し広げる）、今回は2Å -
inside box: 分子の配置領域（初期座標の制約）
一片の長さが60Åの箱を作り、そこにNPD分子100個をランダムに配置しています。後程圧縮して密度を計算するので、最初は少しスカスカの状態で設定しています。
出力PDBに 座標の最小最大（extrema）と tolerance から推定したCRYST1を自動で書いてくれます。
それでは、下のコマンドでpackmol実行しましょう。
packmol < packmol_100.inp
無事に完了すると、Success!!と最後に表示され、box.pdbファイルが生成されます。中身を見てみましょう。
head -n 20 box.pdb | cat -A
HEADER $
TITLE Built with Packmol $
REMARK Packmol generated pdb file $
REMARK Home-Page: http://m3g.iqm.unicamp.br/packmol$
REMARK$
REMARK CRYST1 info below is (extrema(coordinates) +/- 1.1*tolerance) because no explicit$
REMARK PBCs were defined. To apply PBCs, use the `pbc` keyword.$
CRYST1 64.46 64.49 64.53 90.00 90.00 90.00 P 1 1$
HETATM 1 C UNL A 1 1.492 36.591 31.344 1.00 0.00 C $
HETATM 2 C UNL A 1 2.108 35.529 30.690 1.00 0.00 C $
HETATM 3 C UNL A 1 3.481 35.571 30.417 1.00 0.00 C $
HETATM 4 C UNL A 1 4.263 36.695 30.776 1.00 0.00 C $
HETATM 5 C UNL A 1 3.616 37.763 31.425 1.00 0.00 C $
HETATM 6 C UNL A 1 2.244 37.706 31.712 1.00 0.00 C $
HETATM 7 N UNL A 1 5.687 36.714 30.491 1.00 0.00 N $
HETATM 8 C UNL A 1 6.189 36.154 29.252 1.00 0.00 C $
HETATM 9 C UNL A 1 5.411 36.139 28.089 1.00 0.00 C $
HETATM 10 C UNL A 1 5.877 35.533 26.918 1.00 0.00 C $
HETATM 11 C UNL A 1 7.169 34.980 26.827 1.00 0.00 C $
HETATM 12 C UNL A 1 7.969 35.052 27.975 1.00 0.00 C $
生成したこのPDBファイルをそのままmoltemplateで使用すると、余計な行がある影響でエラーが起こってしまうので、下記コマンドで修正します。
CR=$(grep '^CRYST1' box.pdb | head -n 1)
awk '/^(ATOM|HETATM)/{print}' box.pdb > box_atoms.tmp
{ echo "$CR"; cat box_atoms.tmp; } > box.pdb.new
mv box.pdb.new box.pdb
rm -f box_atoms.tmp
何をしているかを解説すると、
-
CR=$(grep '^CRYST1' box.pdb | head -n 1)
box.pdbの中から行頭がCRYSTの行だけを抜き出し、そのうち最初の1行のみを抽出して、シェル変数CRに保存する。
→こうすることで、このPDBに含まれる「セルサイズ定義(CRYST1)」を1行だけ確保できる -
awk '/^(ATOM|HETATM)/{print}' box.pdb > box_atoms.tmp
box.pdbの中から行頭がATOMまたはHETATMの行だけを抜き出して、それを一時ファイルに保存
→ヘッダ類を捨てて、原子座標行だけを抽出している -
{ echo "$CR"; cat box_atoms.tmp; } > box.pdb.new
先ほど抽出した「CRYST1行」と「原子座標PDB」を出力し、それらを連結した結果をbox.pdb.newに書き出して、それを元のbox.pdbに置き換える
このような処理をすることで余計な行を消去でき、moltemplateで使用できるようになります。再度中身を見てみましょう。
head -n 20 box.pdb | cat -A
CRYST1 64.46 64.49 64.53 90.00 90.00 90.00 P 1 1$
HETATM 1 C UNL A 1 1.492 36.591 31.344 1.00 0.00 C $
HETATM 2 C UNL A 1 2.108 35.529 30.690 1.00 0.00 C $
HETATM 3 C UNL A 1 3.481 35.571 30.417 1.00 0.00 C $
HETATM 4 C UNL A 1 4.263 36.695 30.776 1.00 0.00 C $
HETATM 5 C UNL A 1 3.616 37.763 31.425 1.00 0.00 C $
HETATM 6 C UNL A 1 2.244 37.706 31.712 1.00 0.00 C $
HETATM 7 N UNL A 1 5.687 36.714 30.491 1.00 0.00 N $
HETATM 8 C UNL A 1 6.189 36.154 29.252 1.00 0.00 C $
HETATM 9 C UNL A 1 5.411 36.139 28.089 1.00 0.00 C $
HETATM 10 C UNL A 1 5.877 35.533 26.918 1.00 0.00 C $
HETATM 11 C UNL A 1 7.169 34.980 26.827 1.00 0.00 C $
HETATM 12 C UNL A 1 7.969 35.052 27.975 1.00 0.00 C $
HETATM 13 C UNL A 1 7.479 35.604 29.168 1.00 0.00 C $
HETATM 14 C UNL A 1 7.628 34.345 25.570 1.00 0.00 C $
HETATM 15 C UNL A 1 7.111 34.771 24.337 1.00 0.00 C $
HETATM 16 C UNL A 1 7.414 34.119 23.150 1.00 0.00 C $
HETATM 17 C UNL A 1 8.305 33.039 23.113 1.00 0.00 C $
HETATM 18 C UNL A 1 8.887 32.637 24.326 1.00 0.00 C $
HETATM 19 C UNL A 1 8.545 33.277 25.537 1.00 0.00 C $
余計なヘッダー部分が消去されて、きちんと整形されていることが分かります。
moltemplate用のltテンプレを作る
moltemplateをインストールしたときに同梱されているmol22ltを使って、mol2をltに変換します。ただ、このpythonファイルを使っただけでは、moltemplate.shを実行したときにGAFFを継承していない影響でエラーが出るので、sedのコマンドで、inherits GAFFというのを付け加える必要があります。下記コマンドを実行します。
python "$CONDA_PREFIX/bin/mol22lt.py" < npd.mol2 > npd.lt
name=$(head -n 1 npd.lt | awk '{print $1}')
cp npd.lt npd.lt.orig
sed -i "0,/^${name}[[:space:]]*{/s//${name} inherits GAFF {/" npd.lt
これも解説すると、
-
python "$CONDA_PREFIX/bin/mol22lt.py" < molecule.mol2 > molecule.lt
mol22lt.pyを使ってmol2をlt形式に変換、これにより原子、型（GAFFのatom type : ca, cp, nh...)、部分電荷、結合などの情報が.ltに継承されます。 -
name=$(head -n 1 npd.lt | awk '{print $1}')
LTファイルの１行目を抜き取り、その行の第１トークン（空白区切りの1個目）も取って、それをnameに保存。このとき、$1はnpdなので、name=npdとなっています。 -
sed -i "0,/^${name}[[:space:]]*{/s//${name} inherits GAFF {/" npd.lt
inherits GAFFを先頭のクラス定義に差し込む。具体的には
npd {
を
npd inherits GAFF {
に変換しています。inherits GAFFを付けることで、GAFF側で定義されているpair_coeff / bond_coeff / angle_coeffなどの係数設定をnpdクラスが参照できるようになります。
中身を見てみます。
head -n 30 npd.lt
npd inherits GAFF {
# atomId molId atomType charge X Y Z
write("Data Atoms"){
$atom:C1 $mol:m @atom:ca -0.062132 1.133 -0.37 -1.836
$atom:C2 $mol:m @atom:ca -0.060383 0.849 0.841 -1.214
$atom:C3 $mol:m @atom:ca -0.040234 -0.436 1.094 -0.719
$atom:C4 $mol:m @atom:ca 0.037828 -1.473 0.141 -0.858
$atom:C5 $mol:m @atom:ca -0.040234 -1.166 -1.072 -1.502
$atom:C6 $mol:m @atom:ca -0.060383 0.128 -1.326 -1.979
$atom:N1 $mol:m @atom:nh -0.269176 -2.797 0.421 -0.33
$atom:C7 $mol:m @atom:ca 0.037871 -3.353 1.756 -0.428
$atom:C8 $mol:m @atom:ca -0.039617 -3.001 2.619 -1.471
$atom:C9 $mol:m @atom:ca -0.052498 -3.487 3.929 -1.518
$atom:C10 $mol:m @atom:cp -0.018265 -4.407 4.418 -0.57
$atom:C11 $mol:m @atom:ca -0.052498 -4.806 3.531 0.438
$atom:C12 $mol:m @atom:ca -0.039617 -4.27 2.237 0.521
$atom:C13 $mol:m @atom:cp -0.018265 -4.901 5.812 -0.651
$atom:C14 $mol:m @atom:ca -0.052498 -4.975 6.461 -1.893
$atom:C15 $mol:m @atom:ca -0.039617 -5.306 7.805 -2.0
$atom:C16 $mol:m @atom:ca 0.037871 -5.647 8.568 -0.876
$atom:C17 $mol:m @atom:ca -0.039617 -5.644 7.922 0.371
$atom:C18 $mol:m @atom:ca -0.052498 -5.268 6.566 0.48
$atom:N2 $mol:m @atom:nh -0.269176 -5.944 9.977 -1.028
$atom:C19 $mol:m @atom:ca 0.037828 -6.68 10.445 -2.192
$atom:C20 $mol:m @atom:ca -0.040234 -7.624 9.631 -2.849
$atom:C21 $mol:m @atom:ca -0.060383 -8.277 10.065 -4.006
$atom:C22 $mol:m @atom:ca -0.062132 -8.044 11.345 -4.499
$atom:C23 $mol:m @atom:ca -0.060383 -7.149 12.186 -3.846
問題なく処理が行われていることが分かります。
次に、system.ltを作成します。
cat > system.lt << 'EOF'
import "gaff.lt"
import "npd.lt"
NPD = new npd [100]
EOF
system.ltのimport部分は、「使う力場」と「使う分子テンプレ」に必ず一致させる必要があります。そして、[100]の部分もpackmolで配置させた分子数と揃える必要があります。
system.ltにおける NPD = new npd のnpdの部分と、npd.lt におけるnpd inherits GAFFのnpdは一致させる必要があります。違う分子を扱う際は、この部分を揃えないとmoltemplate実行時にエラーが出ます。分子名が長いものだとmol22lt.pyでltファイルを作ったときに、名前が省略されていることがあるので、ちゃんと中身を確認してください。
今回はnpdという名前で実行していますが、一般名で統一してしまえばこのようなエラーを回避できるかと思います(M001とか)。
そして、moltemplate.shを実行しましょう。
moltemplate.sh -atomstyle full -pdb box.pdb system.lt
無事に操作が完了すると、同じディレクトリに、system.data, system.in.init, system.in.settingが新たに生成されます。これらがLAMMPSのインプットファイルになります。
3. LAMMPS計算実行
lammpsの計算はいくつかのstepsに分かれており、今回は「エネルギー最小化→NVT計算→NPT計算」とします。各計算の簡単な説明としては、
-
エネルギー最小化
速度や時間の概念はなく、単にポテンシャルエネルギーが小さくなるように、分子の反発が小さくなるように原子座標を移動させる。
-
NVT計算
速度を与えることで温度が定義される。サーモスタット(熱浴)と系でエネルギー交換しながら、カノニカル分布に緩和していく。NPT計算を行う際に圧力揺らぎを小さくするための前処理として利用できる。
-
NPT計算現実の環境(1atm, 300K)に最も近い
サーモスタット＋バロスタットを使って、温度を保ちつつ、箱サイズを変えて圧力を合わせる。結果として、セルが膨張・収縮して、密度、体積、構造が物理的に妥当な値に緩和させるのが目的。
となります。最後のNPT計算により、現実の環境に近い状態でシミュレーションを行うことで、得られる物性値が現実的に妥当な値に収束していきます。
エネルギー最小化
それでは、まずエネルギー最小化をやっていきます。エネルギー最小化用のインプットファイル in.minを以下の様に作ります。
cat > in.min << 'EOF'
# --- Minimization ---
units real
atom_style full
boundary p p p
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes
# Force-field styles + kspace + special_bonds etc.
include system.in.init
# Topology + coefficients
read_data system.data
include system.in.settings
# Output
thermo_style custom step time temp pe ke etotal enthalpy press vol density lx ly lz pxx pyy pzz pxy pxz pyz
thermo_modify flush yes
thermo 10
dump 1 all custom 100 lmp.dump id mol type xs ys zs ix iy iz xu yu zu
dump_modify 1 sort id
restart 1000 restart.1 restart.2
# Minimization
minimize 1e-4 1e-6 5000 1000000
write_restart lmp.start
write_data lmp_final.data
EOF
中身を一つずつ解説していきます。
-
units: 単位系の設定、今回はrealとしている。エネルギーはkcal/mol, 距離はÅ, 圧力はatmになる。
-
atom_style: 原子スタイルの設定。全てのタイプの相互作用を使用するfullに設定しており、電荷、分子ID、角度など分子トポロジー情報を全部扱える形式。
-
boundary: 境界条件の設定、x, y, z 全方向で周期境界条件を指定する。バルク物質内部の物性を見ることができ、表面効果を回避できる。
-
neighbor: 近傍リストの距離を2.0Åにしている。binはセル分割方式で近傍リストを作る方式（高速）
-
neighbor_modify: 近傍リストの更新タイミングの指定、delay 0で最初から更新を許可し、every 1で1ステップごとに更新をチェックし、check yesで原子が十分動いたら自動で再構築させる。最小化の場合、原子が大きく動く可能性があるため、このような設定にする。
-
include: 力場、クーロン計算（kspace）、special_bondsなどの基本設定を別ファイルから読み込む。ここでは、system.in.initとsystem.in.settingを読み込んでおり、moltemplateが自動生成している部分となる。system.in.initの中身を覗いてみると、
bond_style hybrid harmonic
angle_style hybrid harmonic
dihedral_style hybrid fourier
improper_style hybrid cvff
pair_style hybrid lj/charmm/coul/long 9.0 10.0 10.0
kspace_style pppm 0.0001
pair_modify mix arithmetic
special_bonds amber
となっている。これも一つずつ解説していくと、
-
bond_style: 結合相互作用の設定、harmonicは単純な調和ポテンシャルという意味で、分子内の結合長は「バネ」で拘束するということ。hybridは複数のbond_styleを混在させるための指定だが、今回は実質harmonicしか使用しない（moltemplateにより自動生成されているだけ）
-
angle_style: 結合角相互作用の設定、こちらも調和ポテンシャルで拘束している
-
dihedral_style: 二面角相互作用の設定、今回はfourier型にしている
-
improper_style: improper(二面角の1種、平面性やキラリティ保持用）のポテンシャル、cvffはCVFF力場由来のimproper形式で、芳香環の平面性やsp2中心の平面構造を保つために使われる
-
pair_style: 非結合相互作用（van der Waals + クーロン）の指定、lj/charmm/cou/longはLennard-Jones (12-6)ポテンシャル + クーロン相互作用（長距離部分はkspaceで処理）という意味で、引数として、9.0がLJのカットオフ距離、10.0はクーロン相互作用のカットオフ距離、最後の10.0は近傍リスト用の距離となる
-
kspace_style: 長距離クーロン相互作用の計算方法を指定、pppm(Particle-Particle Particle-Mesh)はEwald系の高速アルゴリズムで、周期境界系でのクーロン相互作用を効率よく正確に計算できる
-
pair_modify: 異種原子間のLJパラメータの混合則を指定、arithmeticは算術平均という意味、他にはgeometricもある
-
special_bonds: 分子内の1-4相互作用の設定、今回はamberルールとしており、1-2（結合原子）は非結合相互作用を除外、1-3（角度でつながる原子）も除外、1-4（二面角でつながる原子）では、LJとCoulombに係数を描けてスケールして入れている
ということです。system.in.initの内容はin.minの中に直接記載しても構いませんが、moltemplateが自動で作ってくれるので、今回はそのまま使います。
また、system.in.settingには、各原子タイプ・結合タイプに対する力場パラメータ（係数）が入っており、pair_coeff, bond_coeff, angle_coeff, dihedral_coeffなどがそれに該当します。
-
read_data: 原子座標、分子構造、結合情報などを含むsystem.dataを読み込む。Packmol + moltemplate で作った系の本体がここで入る
-
thermo_style:thermo出力（画面に出てくるログ）の項目をカスタム指定している- step : ステップ数
- time : 時間
- temp : 温度（※最小化では意味は薄いが出力可能）
- pe, ke, etotal : ポテンシャル、運動、全エネルギー
- enthalpy : エンタルピー
- press : 圧力
- vol : 体積
- density : 密度
- lx ly lz : ボックスサイズ
- pxx pyy pzz pxy pxz pyz : 応力テンソル成分
-
thermo_modify:flush yesとすることで出力をバッファせず、毎回すぐファイル/画面に書き出す。ジョブが落ちたときでも、途中までのログが残りやすい。
-
thermo: 10ステップごとにthermo情報を出力させる。最小化の進行を細かく設定したいので、やや高頻度。
-
dump 1 all custom 100 lmp.dump id mol type xs ys zs ix iy iz xu yu zu:1はこのdumpのID(ラベル)、allで出力対象の原子グループを全原子に設定し、100ステップごとに座標を lmp.dump に出力。出力項目として、- id, type, mol : LAMMPSが各原子に割り振っている一意な番号(原子ID), どの分子に属している原子か(分子ID), 力場が定義した種類(type)
- xs ys zs : 0–1 に正規化されたスケール座標（scaled）
- ix iy iz : 周期境界のイメージフラグ（unwrap用）
- xu yu zu : 周期境界をまたいだ分も含めたアンラップ座標
-
dump_modify: dump 出力を 原子ID順にソート。可視化や後処理で原子順が安定して扱いやすくなる。
-
restart: 1000ステップごとにリスタートファイルを書き出す。2ファイルを交互に使う方式（安全対策）。万一途中で計算が落ちても、ここから再開可能。
-
minimize: エネルギー最小化の実行。引数の意味は、- 1e-4 : エネルギー収束条件（|ΔE|）
- 1e-6 : 力の収束条件（|F|）
- 5000 : 最大イテレーション数
- 1000000 : 最大評価回数
目的はPackmol 初期構造の過大な反発や不自然な接触を解消するため。
-
write_restart lmp.start: 最小化後の状態を リスタートファイルとして保存。次の NVT/NPT 計算の初期構造としてそのまま使える。
-
write_data lmp_final.data: 最小化後の構造を data ファイル形式でも保存。別条件の計算や検証用に再利用しやすい。
といった具合です。長くなりましたが、それでは計算してみます。下記コマンドで計算開始です。
lmp -in in.min | tee log.min
途中経過が出力され、Total wall timeが表示されると、無事に計算が終了します。
実行後に、下記コマンドを入力して、出力されたlog.minの中身を見てみると計算終了を再確認できます。
tail -n 40 log.min
Comm | 0.047933 | 0.047933 | 0.047933 | 0.0 | 0.11
Output | 0.14013 | 0.14013 | 0.14013 | 0.0 | 0.32
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.1744 | | | 0.40
Nlocal: 7800 ave 7800 max 7800 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 4947 ave 4947 max 4947 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 1.45197e+06 ave 1.45197e+06 max 1.45197e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 1451974
Ave neighs/atom = 186.15051
Ave special neighs/atom = 10.538462
Neighbor list builds = 23
Dangerous builds = 0
System init for write_restart ...
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.14161198
grid = 9 9 9
stencil order = 5
estimated absolute RMS force accuracy = 0.032665926
estimated relative force accuracy = 9.8372466e-05
using double precision FFTW3
3d grid and FFT values/proc = 2744 729
Generated 3403 of 3403 mixed pair_coeff terms from arithmetic mixing rule
System init for write_data ...
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.14161198
grid = 9 9 9
stencil order = 5
estimated absolute RMS force accuracy = 0.032665926
estimated relative force accuracy = 9.8372466e-05
using double precision FFTW3
3d grid and FFT values/proc = 2744 729
Generated 3403 of 3403 mixed pair_coeff terms from arithmetic mixing rule
Total wall time: 0:00:43
エネルギー最小化では、初期配座から分子の反発が小さくなるように原子座標を移動させるのが目的でした。ですので、ポテンシャルエネルギーが小さくなっているか確認する必要があります。log.minをエクセルにドラッグして、Time vs PotEng をプロットしてみてください。
こんな感じでエネルギーが減少し、収束気味ならOKです。NVT計算用に配座を安定化させることに成功しました。
NVT計算
次にNVT計算用のinputファイルを以下の様に作ります。
cat > in.nvt << 'EOF'
# --- NVT ---
units real
atom_style full
boundary p p p
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes
include system.in.init
read_data lmp_final.data
include system.in.settings
velocity all create 300. 12345
fix 1 all shake 1e-5 1000 0 m 1.0
fix 2 all nvt temp 300.0 300.0 100.0 tchain 3
fix 3 all momentum 50 linear 1 1 1
timestep 2.0 # fs
thermo_style custom step time temp pe ke etotal enthalpy press vol density lx ly lz
thermo_modify flush yes
thermo 10
dump 1 all custom 200 nvt.dump id mol type xs ys zs ix iy iz xu yu zu
dump_modify 1 sort id
run 5000 upto
unfix 3
unfix 2
unfix 1
write_restart nvt.end
write_data nvt_final.data
EOF
エネルギー最小化と共通部分の説明は省略し、新しく表れた入力スクリプトを解説します。
-
velocity all create 300. 12345: 全原子に300Kに相当するランダムな初期速度を与える、12345は乱数シード
-
fix 1 all shake 1e-5 1000 0 m 1.0: 質量1の原子＝水素原子を含む結合をSHAKEアルゴリズムで固定、C-H結合などの高速振動を凍結し、タイムステップを大きくする
-
fix 2 all nvt temp 300.0 300.0 100.0 tchain 3: NVTアンサンブルでの温度制御、開始温度300K, 目標温度300K, 緩和時間100fsとしている
-
fix 3 all momentum 50 linear 1 1 1: 50ステップごとに系全体の重心速度をゼロにする、これにより系全体が箱の中を流れるドリフトを除去できる
それでは、下記コマンドでNVT計算をしましょう。
lmp -in in.nvt | tee log.nvt
計算が無事に完了したら、出力されたlog.nvtをエクセルで開いて、Time vs Tempをプロットしてみます。NVT計算は温度一定にすることが目的でした。
時間が進むにつれて、温度が300Kで安定していることが分かります。これで、本計算(NPT)前の前準備が終わりました。
NPT計算
最後に、このようなインプットファイルを作りましょう。
cat > in.npt << 'EOF'
# --- NPT ---
units real
atom_style full
boundary p p p
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes
include system.in.init
# start from NVT result
read_data nvt_final.data
include system.in.settings
fix 1 all shake 1e-5 1000 0 m 1.0
fix 2 all npt temp 300.0 300.0 100.0 tchain 3 iso 1.0 1.0 100.0 pchain 3
fix 3 all momentum 50 linear 1 1 1
timestep 2.0 # fs
thermo_style custom step time temp press density vol lx ly lz pe ke etotal
thermo_modify flush yes
thermo 10
dump 1 all custom 200 npt.dump id mol type xs ys zs ix iy iz xu yu zu
dump_modify 1 sort id
run 100000 upto
write_restart npt.end
write_data npt_final.data
EOF
NVT計算と大きくは変わらないのですが、一番大きな違いは、NPT計算用にfixを書き換えている点です。今回は、温度300K, 圧力1atmで等方的に体積を変化させるNPT計算を行いたいので、以下の様なスクリプトになります。
-
fix 2 all npt temp 300.0 300.0 100.0 tchain 3 iso 1.0 1.0 100.0 pchain 3:- Nosé–Hoover サーモスタット + バロスタット となり、粒子数N、圧力P、温度Tが一定のNPT計算
- 開始温度と目標温度が300Kで、温度緩和時間が100fs
-
tchainはNosé–Hooverチェーン長で熱浴の自由度を表す、1だと単純で3が標準 - 等方的に体積変化させたいので
isoとして、開始圧力と目標圧力を1atomとする、圧力緩和時間も100fs -
pchainは圧力側のチェーン長でこれも3としている
ちなみに、200000fsくらいまで見ないと収束しないので、time(fs) = timestep(2.0) x step(100000) としています。
それでは計算実行します。
lmp -in in.npt | tee log.npt
NPT計算は現実の環境に最も近い状態でのシミュレーションですので(300K, 1atm)、得られる物性値は物理的に妥当な値に近づきます。
よく見られる指標が密度です。有機物の密度はだいたい1くらいですので、NPT計算が問題なく完了すれば、密度がこの値に収束しているはずです。それでは、得られたlog.nptをエクセルに張り付け、time vs densityをプロットしてみます。
密度が1に収束していることが分かりました！！これにてNPT計算完了で、MD計算一連の流れを実行することができました。
4. 可視化
NPT計算が無事に完了しましたが、本当に系が収縮して密度が高くなっているかどうか、可視化ソフトを使って確認しようと思います。フリーの可視化ソフトは調べてみると、VMDとOVITOというものがありました。個人的にOVITOの方が使いやすかったので、簡単な操作方法をまとめておきます。
下記サイトに行くと、OVITOがダウンロードできます。無料版はbasicです（proはおそらく有料）。
ソフトのダウンロードが完了した後に開いてみると、こんな感じの初期画面に行きます。ここにNPT計算後のnpt.dumpをドラッグして読み込ませてください。File→Load Fileからでもいいです。
すると、こんな感じで立方体と詰まっている分子が表示されます。
4つのwindowsが表示されているので、1つに注目してみたい場合は、そのwindowをクリックして拡大のボタンを押します。今回はPerspectiveをクリックして拡大しました。
NPT計算で密度が向上していく様子を見る場合は、下の再生ポタン▶をクリックしてください。すると、このようにアニメーションで確認できます。
まだまだいろんな機能があるみたいですが、色々遊んでみてください。
まとめ
いかがでしたか？MD計算の一連の流れを実際にやってみました。これが基本的な流れで、さらに別のparameterが必要な時は追加の計算をします。ただ、最初のSet-Upは大変ですが、一度できるようになると、意外と簡単にMD計算ってできるんだなと実感できます。ぜひぜひ興味ある方やってみてください。
今回も長文読んでいただき有難うございました！！
引用
・「LAMMPSと連携ソフトウエアによる有機材料の分子動力学計算」
https://makoto-yoneya.github.io/LAMMPS-organics/
・有機分子の分子動力学シミュレーション方法
https://qiita.com/awaawa/items/8d84b46218777504f1d2
・antechamberを使ったMD用小分子化合物の電荷パラメータの作り方
https://qiita.com/Ag_smith/items/430e9efb32a855d4c511