檔案建立_LAMMPS data檔案建立工具--moltemplate
技術標籤：檔案建立
一 介紹
長久以來，lammps的資料檔案構建就是一個大問題（對我來說）。一方面，LAMMPS的資料格式比較特殊，很多軟體不能直接匯出; 另一方面小分子到大分子的變換，大分子在系統中的排布和大分子中的拓撲結構都需要牽扯很多的精力。經過長時間的試錯，我算是終於找到一個比較合理且簡約的技術棧，可以直接生成lammps所需的in檔案和data檔案。
如圖所示，整個技術棧分為：
- 基本繪圖：繪製基本的單元，如高分子結構單元，單個水分子或者甲烷分子等;
- 資料檔案：匯出成常見格式的資料檔案，如PDB，XYZ;
- [空間排布]：packmol可以自動地根據所設定的限制向系統中填放單元。比如想向C60籠中填放原子，只需要一行命令將允許範圍限制在球內即可;
- 生成工具：moltemplate可以操作已有的資料檔案，在系統中複製，變換和隨機，同時根據已有的拓撲結構補充生成鍵角二面角等資訊，最後生成lammps所需的in檔案和data檔案。
由此，得到的所有檔案都可以直接使用lammps執行。此套技術棧的好處是，可以接受各種格式的資料檔案並且可以輕鬆地操縱它們，而且都是開源的，安裝方便快捷佔地小。開發moltemplate的人同樣也是lammps的主程，在手冊中他成這款軟體是lammps的前端（front-end）
唯一的缺點是，moltemplate（包括packmol）手冊有如16年以前的lammps手冊，比較 ……隨意灑脫。
二 檔案結構
moltemplate本身具有匯入（import）和繼承（inherit）功能，因此我們可以儘量地把檔案區分按照功能區分開。當然檔名都是可以自定義的，這裡只是介紹一個思路，告訴各位應該怎樣上手。
第一部分是forcefield.lt 這個檔案中儲存著lammps的系統和力場引數。
# -- ForceField -- # ForceField{ write_once("In Init"){ units lj boundary p p p atom_style full pair_style lj/cut 10.5 bond_style hybrid fene harmonic angle_style harmonic } write_once("Data Boundary") { 0 100.0 xlo xhi 0 100.0 ylo yhi 0 100.0 zlo zhi } write_once("Data Masses"){ @atom:M 1 @atom:G 1 @atom:R 1 } write_once("In Settings"){ pair_coeff @atom:M @atom:M 1 1 2.5 pair_coeff @atom:M @atom:G 1 1 2.5 pair_coeff @atom:M @atom:R 1 1 2.5 pair_coeff @atom:G @atom:G 1 1 2.5 pair_coeff @atom:G @atom:R 1 1 2.5 pair_coeff @atom:R @atom:R 1 1 2.5 } write_once("In Settings"){ bond_coeff @bond:MM fene 30 1.5 bond_coeff @bond:GR fene 30 1.5 bond_coeff @bond:RR harmonic 1000 0.66 } write_once("In Settings"){ angle_coeff @angle:RRR 200 180 } # write_once("Data Dihedrals By Type"){ # @dihedral: @atom: # } # write_once("Data Angles By Type"){ # @angle: @atom: # } }
第二部分是基本單元資料 這裡面儲存這一個片段的資訊
# -- file matrix.lt --
import "forcefield.lt"
Matrix inherits ForceField{
write("Data Atoms"){
$atom:1 $mol:. @atom:M 0 1 1 1
$atom:2 $mol:. @atom:M 0 2 1 1
$atom:3 $mol:. @atom:M 0 3 1 1
$atom:4 $mol:. @atom:M 0 4 1 1
$atom:5 $mol:. @atom:M 0 5 1 1
}
write("Data Bonds"){
$bond:1 @bond:MM $atom:1 $atom:2
$bond:2 @bond:MM $atom:2 $atom:3
$bond:3 @bond:MM $atom:3 $atom:4
$bond:4 @bond:MM $atom:4 $atom:5
}
}第三部分是操作部分 在這個檔案中我們可以對基本單元進行操作，同時也是程式執行的入口
# -- file system.lt -- #
import "forcefield.lt"
import "matrix.lt"
m1 = new Matrix [3].move(0,0,3)在最簡單的情況下，我們僅需要這三個檔案就可以描述一個系統的初始化狀態，然後執行命令列使用moltemplate生成實際檔案：
moltemplate.sh system.lt得到：
system.data
system.in
system.in.init
system.in.settings
# 如果出了bug還會出現生成樹output_ttree三 語法教程
1 forcefield.lt
# -- example of a forcefield.lt inscript -- #
ForceField{
write_once("In Init"){
units lj
boundary p p p
atom_style full
pair_style lj/cut 10.5
bond_style hybrid fene harmonic
angle_style harmonic
}
# write_once("Data Angles By Type"){
# @angle: @atom:
# }
}1.1 write_once() 和 write()
write_once("file_name"){
text_content
}每出現一次write_once()，就會新生成一個檔案，整個程式執行過程中此命令僅僅執行一次。注意，In 和 Data 開頭是系統保留關鍵詞,將會被合併到in檔案和data檔案中，而特別離譜的才會被單獨創立。比如，In Init就是個保留檔名，如果你輸入In init就會被強制要求改成大寫。接下來會有專門的專題來介紹保留檔名的特殊之處。還有一些檔名會被合併，比如Data Mass、Data boundary。花括號內的文字將會被轉移到檔案內，同時計數器都會被替換。
write則是可以在一個文字塊內出現多次，可以被重複執行的單元（通過new）
1.2 import 和 inherits
import 關鍵字要求系統去查詢需要引用的檔案的位置
inherits 表示本結構單元中的變數名稱（比如bond type名MM來自於繼承的單元）
1.3 forcefield 和 By Type
write_once("Data Angles By Type"){
@angle: type @atom:atom1 @atom:atom2 @atom:atom3 @bond：b1 @bond：b2
}這裡介紹moltemplate最終要的功能之一。眾所周知我們的邏輯結構不單單包括2-body的鍵接，還有3-body，4-body的angle和dihedral甚至還有improper。如果自己寫程式去搜索拓撲結構的話會非常費事。因此這裡給出了一個功能，只要你給出了結構的鍵接資訊，那麼就可以自動去查詢所有的angle和dihedral和improper。其中語法也十分清晰明瞭，第一，指明這個angle的三個原子型別，第二，指明兩個鍵接型別，僅此就可以遍歷所有原子，標記出所有的角。
2 unit.lt
# -- example of a unit.lt inscript -- #
import "forcefield.lt"
Matrix inherits ForceField{
write("Data Atoms"){
# atom id mol id type charge coordinate
$atom:1 $mol:. @atom:M 0 1 1 1
$atom:2 $mol:. @atom:M 0 2 1 1
$atom:3 $mol:. @atom:M 0 3 1 1
$atom:4 $mol:. @atom:M 0 4 1 1
$atom:5 $mol:. @atom:M 0 5 1 1
}
write("Data Bonds"){
$bond:1 @bond:MM $atom:1 $atom:2
$bond:2 @bond:MM $atom:2 $atom:3
$bond:3 @bond:MM $atom:3 $atom:4
$bond:4 @bond:MM $atom:4 $atom:5
}
}2.2 $和@計數器
如上面的$atom：1和@atom：1，都被成為變數。前面的符號是計數器，後面叫做類名（category），冒號後則是變數名，我們很快就會知道這個為什麼要這麼寫。
計數器分為兩種，$是動態計數器，每增加一個，在生成的時候都會賦給一個唯一的id，@是靜態計數器，每一個類名在全域性只會唯一給一個type。簡單說，就是靜態計數器和系統複雜度有關，動態計數器和系統規模有關。
2.3 變數作用域
上面提到的變數有兩種形式，一種是全名，另一種是簡名。簡名通常是在一個結構物件（molecule-object）內，不會產生歧義。而全名就厲害了：
@cpath/catname：lpath
cpath：category的作用域，通常省略以表示全域性。除了有極為特殊的要求之外不需要管
catname：類名
lpath：包括變數名和可選的指明那個分子的index這樣可以從一個結構物件內引用其他結構物件中的原子或者結構片段。但為了清晰起見，按照程式設計的基本原則，請不要交叉引用。
2.4 $mol:. 和 $mol:...
你所看到的點，斜線的用法和unix作業系統檔案路徑表示方法一樣，既可以用來表示相對路徑，又可以表示絕對路徑。所謂路徑就是指的是不同的結構物件之間的關係（molecule-object）。類似的，一個點意為“這個”，就是代指本結構，如果一個結構單元中所有的原子都是$mol:.，那每new一次，同一個結構片段內mol id相同，不同結構片段mol id遞增
而三個點則是告訴系統，我這個片段從屬於一個大分子，那麼這個片段在大分子的結構片段內不管new幾次，new一次大分子得到的所有原子從屬於一個mol id，不同大分子的mol id才不同。
3 system.lt
這一部分就是對之前已構建的所有單元片段的操縱了，不論是旋轉，平移還是隨機都是一行搞定。這些命令都在手冊的一開頭。
關於Debug
每次檔案生成結束，都會給出一個名為output_ttree的資料夾，裡面有各種檔案。其中最終要的是就是ttree_assignments，這裡面儲存了計數器被替換的數和該計數器第一次出現的位置。