1. 為什么說(shuō)in文件才是LAMMPS真正的門(mén)檻每次有新人來(lái)問(wèn)LAMMPS怎么入門(mén)我第一句話都是先把安裝這關(guān)過(guò)了然后立刻把重心放到in文件上。LAMMPS的安裝本身并不復(fù)雜網(wǎng)上有大量編譯好的二進(jìn)制包Windows、Linux、macOS都有真正讓大批初學(xué)者卡住兩周、甚至一個(gè)月的從來(lái)不是可執(zhí)行文件能不能跑起來(lái)而是那個(gè)決定模擬全過(guò)程的輸入文件——in文件。很多人對(duì)分子動(dòng)力學(xué)模擬的理解是“把原子坐標(biāo)丟進(jìn)去然后等著出結(jié)果”這個(gè)認(rèn)知偏差很致命。LAMMPS本身是一個(gè)沒(méi)有任何內(nèi)置物理模型的框架它不知道你要模擬水分子、銅納米線還是石墨烯剪切它只知道按照你在in文件里寫(xiě)的命令去讀數(shù)據(jù)、計(jì)算受力、更新坐標(biāo)。換句話說(shuō)in文件不是一份“配置文件”而是一份完整的、逐行執(zhí)行的模擬程序源碼。你用什么單位制讀入數(shù)據(jù)選哪個(gè)勢(shì)函數(shù)計(jì)算原子間作用力用哪種系綜控制溫度壓力每一步輸出什么信息全部由in文件里的每一條命令決定。這也是我寫(xiě)這篇文章的原因。市面上雖然有很多LAMMPS教程但大多分成兩種極端一種是官方手冊(cè)的翻譯式教學(xué)把每個(gè)命令的語(yǔ)法列一遍看了等于沒(méi)看另一種是論文復(fù)現(xiàn)式的命令堆疊直接給你一份能跑的in文件但完全不解釋為什么這樣寫(xiě)換個(gè)體系就抓瞎。我的目標(biāo)是走中間路線——以一份完整的拉伸模擬in文件為載體把文件的邏輯區(qū)塊、每個(gè)關(guān)鍵命令背后的物理邏輯、以及我在實(shí)際調(diào)試中踩過(guò)的坑全部講透讓真正需要用它來(lái)干活的人能夠看完之后自己寫(xiě)出適合自己體系的in文件而不是只能復(fù)制粘貼。這篇文章適合哪些人零基礎(chǔ)、剛剛裝好LAMMPS、正被第一份in文件折磨的人有一定經(jīng)驗(yàn)但只會(huì)在網(wǎng)上抄作業(yè)、遇到報(bào)錯(cuò)不知道怎么排查的人以及想把自己的模擬流程從“手工改文件”提升到“結(jié)構(gòu)化組織”的人。我會(huì)盡量把底層邏輯講清楚同時(shí)保留足夠的實(shí)操細(xì)節(jié)保證你能直接跟著步驟把案例跑通。先說(shuō)一個(gè)貫穿全文的核心觀點(diǎn)in文件的本質(zhì)是“組織計(jì)算流程”,而不是“描述物理模型”。理解了這句話你寫(xiě)完一行命令就會(huì)多想一步——它到底是在初始化環(huán)境、構(gòu)建體系、定義計(jì)算還是在控制運(yùn)行節(jié)奏這個(gè)分類(lèi)意識(shí)是后期排查錯(cuò)誤最快的利器。2. 一張in文件的結(jié)構(gòu)地圖從初始化到結(jié)果輸出的五個(gè)邏輯區(qū)塊我見(jiàn)過(guò)太多人拿到一份in文件上來(lái)就逐行去查命令是什么意思這是效率最低的方式。正確的讀法是把整個(gè) in 文件按功能切成區(qū)塊因?yàn)?LAMMPS 的命令執(zhí)行順序?qū)φ_性有決定性影響區(qū)塊之間天然存在嚴(yán)格的先后依賴(lài)關(guān)系。一份標(biāo)準(zhǔn)、完整的in文件無(wú)論體系是液體、固體、納米線還是多相復(fù)合材料骨架通常是這幾大塊初始化區(qū)塊設(shè)置單位制、維度、邊界條件、原子類(lèi)型、勢(shì)函數(shù)種類(lèi)等“全局默認(rèn)值”。建模區(qū)塊讀入或生成原子坐標(biāo)定義區(qū)域、創(chuàng)建晶格、設(shè)置原子種類(lèi)和質(zhì)量、建立近鄰列表參數(shù)。計(jì)算設(shè)置區(qū)塊分配系綜、施加溫度/壓力控制、定義需要統(tǒng)計(jì)的物理量、設(shè)置輸出頻率和格式。弛豫與預(yù)平衡區(qū)塊讓體系達(dá)到穩(wěn)定狀態(tài)消除初始構(gòu)型中的不合理接觸和應(yīng)力集中。生產(chǎn)運(yùn)行區(qū)塊執(zhí)行正式的模擬過(guò)程并在過(guò)程中持續(xù)輸出軌跡和熱力學(xué)數(shù)據(jù)。2.1 初始化區(qū)塊最先執(zhí)行的命令決定了后面所有數(shù)值的“語(yǔ)言”初始化區(qū)塊是整個(gè)文件的地基同時(shí)也是初學(xué)者最容易忽略的地方。很多人裝上LAMMPS后第一件事就是去網(wǎng)上找一份例子復(fù)制下來(lái)改改坐標(biāo)文件就跑結(jié)果把units real和units metal混著用最后算出來(lái)的溫度離譜到幾千K也沒(méi)有意識(shí)到是單位制出了問(wèn)題。常見(jiàn)的初始化命令包括units real dimension 3 boundary p p p atom_style atomic pair_style lj/cut 10.0我先解釋一下為什么這個(gè)區(qū)塊的順序本身就有講究。units必須放在最前面因?yàn)樗x了所有后續(xù)命令中數(shù)值的默認(rèn)單位。units real下能量單位是 kcal/mol距離單位是?溫度直接就是K時(shí)間單位是 fs飛秒而units metal下能量單位變?yōu)?eV距離仍然是?時(shí)間單位變成了 ps皮秒。同一個(gè)數(shù)值 0.5在不同單位制下代表了完全不同的物理場(chǎng)景。最坑的是很多勢(shì)函數(shù)參數(shù)文件里有默認(rèn)的單位制假設(shè)如果你自己的units設(shè)置和勢(shì)文件期望的不一致LAMMPS 不會(huì)報(bào)錯(cuò)但算出來(lái)的一切都是錯(cuò)的——這是同類(lèi)錯(cuò)誤里最難排查的一種。boundary命令設(shè)置的是三個(gè)方向上的邊界條件p表示周期性邊界f表示固定邊界s表示收縮邊界shrink-wrap。絕大多數(shù)晶體塑性、納米壓痕、拉伸模擬用的是p p p三方向周期性理由很簡(jiǎn)單周期性邊界條件下表面效應(yīng)被消除模擬盒子相當(dāng)于無(wú)限大體系的一個(gè)原胞。但如果你要模擬納米線或者薄膜表面必須真實(shí)暴露出來(lái)那至少某一個(gè)方向要用f或者s。我見(jiàn)過(guò)有人想模擬單軸拉伸把兩個(gè)方向設(shè)成p一個(gè)方向設(shè)成f結(jié)果固定邊界方向在原子的熱運(yùn)動(dòng)下產(chǎn)生非物理的應(yīng)力集中整個(gè)模擬早期就直接爆掉了。atom_style這個(gè)命令也很容易被一帶而過(guò)。它定義的是每個(gè)原子存儲(chǔ)哪些屬性。例如atomic只存坐標(biāo)和原子編號(hào)適合簡(jiǎn)單的 LJ 流體charge額外存儲(chǔ)電荷用于帶電位體系molecular存儲(chǔ)分子拓?fù)湫畔⒎肿泳幪?hào)、鍵角連接適合聚合物或分子晶體bond、angle、dihedral等更多樣式對(duì)應(yīng)復(fù)雜力場(chǎng)的全原子模型。初學(xué)者最常見(jiàn)的錯(cuò)誤是在read_data之前沒(méi)有正確設(shè)置atom_style導(dǎo)致 LAMMPS 讀取數(shù)據(jù)文件時(shí)發(fā)現(xiàn)原子屬性數(shù)目對(duì)不上。這個(gè)問(wèn)題我在后文案例里會(huì)再展開(kāi)講一次。pair_style定義了非鍵相互作用的計(jì)算方法。這一條嚴(yán)格來(lái)說(shuō)屬于物理模型的范疇但因?yàn)樗诔跏蓟A段就必須聲明而且與后續(xù)建模區(qū)塊中原子類(lèi)型分配直接關(guān)聯(lián)我習(xí)慣把它歸入初始化區(qū)塊。要注意的是pair_style必須在建模命令前指定。因?yàn)?LAMMPS 在構(gòu)建粒子列表和計(jì)算鄰居關(guān)系時(shí)需要提前知道應(yīng)該為哪些原子對(duì)計(jì)算力而不同勢(shì)函數(shù)對(duì)近鄰列表的存儲(chǔ)方式、截?cái)喟霃降哪J(rèn)處理都有差異。2.2 建模區(qū)塊坐標(biāo)從哪里來(lái)決定了你后面能走多遠(yuǎn)建模區(qū)塊有兩種截然不同的路線外部讀入和內(nèi)部生成。對(duì)應(yīng)兩個(gè)核心命令read_data和create_atoms。路線一從外部文件讀入。這是最常用、也最接近真實(shí)應(yīng)用的方式。你從 Materials Studio、OVITO 的待導(dǎo)出結(jié)構(gòu)或者實(shí)驗(yàn)晶體結(jié)構(gòu)數(shù)據(jù)生成一個(gè)data文件然后在 in 文件里用read_data命令讀入read_data my_system.data此時(shí)要特別注意的是 data 文件內(nèi)部結(jié)構(gòu)必須與 in 文件的atom_style嚴(yán)格匹配。比如你聲明了atom_style charge那 data 文件里的 Atoms 區(qū)域就必須包含電荷這一列如果你的體系有分子、鍵、角那么 data 文件的 Molecule 區(qū)域、Bond 區(qū)域、Angle 區(qū)域也要完整。LAMMPS 查這類(lèi)錯(cuò)誤時(shí)通常會(huì)給出類(lèi)似 “Inconsistent atom style” 的報(bào)錯(cuò)但很多時(shí)候報(bào)錯(cuò)信息出現(xiàn)的位置和真正原因離了十萬(wàn)八千里因?yàn)?data 文件里多個(gè)區(qū)域的字段是連續(xù)解析的一個(gè)錯(cuò)位后面全亂。路線二在 in 文件內(nèi)部生成。適合簡(jiǎn)單晶格、需要快速構(gòu)造規(guī)則晶體的情況。這時(shí)用到lattice、region、create_atoms三兄弟lattice fcc 3.615 region box block 0 10 0 10 0 10 create_atoms 1 boxlattice fcc 3.615定義了一個(gè)晶格常數(shù)為 3.615? 的面心立方晶格。region box block 0 10 0 10 0 10定義了一個(gè)邊長(zhǎng)為10個(gè)晶胞的立方體區(qū)域。create_atoms 1 box則是在這個(gè)區(qū)域內(nèi)按晶格格點(diǎn)生成類(lèi)型為1的原子。這三條命令配合起來(lái)是生成理想單晶最快捷的路徑。但如果你要建模多晶、非晶或帶有缺陷的結(jié)構(gòu)內(nèi)部生成就不夠用了必須依賴(lài)外部建模工具。無(wú)論走哪條路建模后必須做一件事檢查體系里有沒(méi)有原子重疊或原子過(guò)近。很多人嫌麻煩跳過(guò)直接進(jìn)入下一步計(jì)算結(jié)果就是后面fix nvt一開(kāi)體系能量爆到幾萬(wàn) kcal/mol溫度直接崩到幾十萬(wàn) K。所以建模后我通常緊跟一個(gè)簡(jiǎn)單的能量最小化做初步結(jié)構(gòu)優(yōu)化min_style cg minimize 1.0e-4 1.0e-6 1000 10000minimize四個(gè)參數(shù)分別代表能量收斂容差、力收斂容差、最大迭代步數(shù)和最大力計(jì)算次數(shù)。這一步的目的不是徹底的幾何優(yōu)化而是快速排除初始構(gòu)型中的不合理接觸。如果最小化過(guò)程中能量無(wú)法收斂基本可以判定模型構(gòu)建有問(wèn)題這時(shí)回頭看坐標(biāo)文件遠(yuǎn)比等生產(chǎn)模擬跑了一半再排查要高效得多。2.3 計(jì)算設(shè)置區(qū)塊系綜、溫控、輸運(yùn)參數(shù)的關(guān)系體系建好之后next就要決定如何推進(jìn)模擬。這一步對(duì)模擬結(jié)果的物理正確性至關(guān)重要而且也是in文件中最容易“照著別人的模板改但不適配自己體系”的部分。先解決一個(gè)最常見(jiàn)的困惑fix nvt、fix npt、fix nve分別什么時(shí)候用fix nve是牛頓運(yùn)動(dòng)方程的本征積分器不對(duì)溫度和壓力做任何控制能量在統(tǒng)計(jì)意義上守恒適合微正則系綜NVE下的動(dòng)力學(xué)采樣。fix nvt在 nve 的基礎(chǔ)上疊加了一個(gè)溫度耦合項(xiàng)thermostat保持原子數(shù)、體積、溫度不變適合在目標(biāo)溫度下弛豫體系至平衡態(tài)。fix npt同時(shí)控制溫度和壓力允許盒子體積變化適合弛豫到環(huán)境壓強(qiáng)下、或者在恒壓條件下做生產(chǎn)模擬。實(shí)際項(xiàng)目中最常見(jiàn)的組合是先用fix nvt做升溫/恒溫弛豫再用fix npt做等溫等壓平衡最后生產(chǎn)階段根據(jù)研究問(wèn)題換成nvt或者nve。典型的兩段式計(jì)算設(shè)置如下velocity all create 300.0 4928459 loop geom fix 1 all nvt temp 300.0 300.0 0.1 run 20000velocity命令根據(jù)目標(biāo)溫度生成初始速度分布其中l(wèi)oop geom是按原子順序分配隨機(jī)速度以確保質(zhì)心動(dòng)量為零。fix 1 all nvt中的1是 fix 的ID之后的時(shí)間常數(shù)0.1是溫控耦合時(shí)間單位取決于units設(shè)置。新手在這里常犯的錯(cuò)是直接把velocity里的隨機(jī)數(shù)種子設(shè)為固定值每次運(yùn)行得到的初始速度完全相同。重復(fù)性在某些時(shí)候是優(yōu)點(diǎn)但如果你要做統(tǒng)計(jì)分析或者需要多次獨(dú)立采樣保持種子一致會(huì)嚴(yán)重削弱樣本獨(dú)立性。建議隨機(jī)種子每次都改或者用時(shí)間相關(guān)的種子。另外dump和thermo的輸出頻率設(shè)置也是這一區(qū)塊的重要組成部分。thermo 100表示每100步在屏幕上輸出一次熱力學(xué)量溫度、壓能、總能量等dump 1 all custom 1000 dump.lammpstrj id type x y z vx vy vz則是每1000步輸出一個(gè) LAMMPS 軌跡幀包含原子坐標(biāo)和速度。輸出頻率怎么設(shè)牽涉到模擬的“時(shí)間成本”與“數(shù)據(jù)量”的平衡。thermo頻率太高會(huì)拖慢速度雖然現(xiàn)代機(jī)器上影響很小dump頻率太高則會(huì)讓軌跡文件膨脹到幾個(gè)TB。一般經(jīng)驗(yàn)是平衡階段thermo 100dump每1000步一幀生產(chǎn)階段如果復(fù)核能量變化thermo 1000即可dump頻率根據(jù)你想捕捉的物理過(guò)程的特征時(shí)間尺度來(lái)定取樣間隔至少要小于特征時(shí)間一個(gè)量級(jí)。2.4 弛豫、生產(chǎn)運(yùn)行與輸出為什么“跑完”不等于“算完”前面所有區(qū)塊準(zhǔn)備的鋪墊都是為了最后能夠穩(wěn)當(dāng)?shù)匕焉a(chǎn)階段的模擬跑通。但“run 跑完”絕不等于“計(jì)算完成”。在我的工作流里模擬完成后至少還要做三件事第一檢查能量和溫度軌跡是否平穩(wěn)。如果溫度曲線在平衡階段一直在漂移說(shuō)明體系沒(méi)有充分弛豫這種情況下生產(chǎn)階段的結(jié)果是不可信的。第二檢查dump出的軌跡文件在可視化軟件如 OVITO中是否正常有沒(méi)有原子飛出盒子、有沒(méi)有斷鍵后原子漂移到異常位置。第三根據(jù)研究目標(biāo)確認(rèn)是否需要延長(zhǎng)生產(chǎn)時(shí)間。很多初學(xué)者在生產(chǎn)階段只跑 1 ns就急著提取力學(xué)曲線實(shí)際上體系在微正則系綜下可能根本還沒(méi)有達(dá)到穩(wěn)態(tài)。關(guān)于rerun和write_data我覺(jué)得是很多人沒(méi)有用起來(lái)的好命令。rerun允許你在已有軌跡文件上重新進(jìn)行后處理計(jì)算而不必重新跑一遍動(dòng)力學(xué)write_data則可以把當(dāng)前構(gòu)型寫(xiě)到 data 文件里方便以當(dāng)前狀態(tài)為起點(diǎn)做不同的后續(xù)模擬分支。3. 一個(gè)完整案例銅納米線單軸拉伸的in文件逐段拆解光講命令肯定不夠我直接把一份可以用來(lái)做單軸拉伸模擬的完整 in 文件拿出來(lái)逐段拆給你看。這個(gè)案例的核心任務(wù)是對(duì)一根銅納米線施加應(yīng)變計(jì)算應(yīng)力-應(yīng)變響應(yīng)研究其彈性模量和塑性變形機(jī)制。3.1 完整in文件展示可以直接復(fù)制# 銅納米線拉伸模擬 # 單位與全局設(shè)置 units metal dimension 3 boundary f f p atom_style atomic neighbor 0.3 bin neigh_modify delay 0 every 1 check yes # 勢(shì)函數(shù) pair_style eam pair_coeff * * Cu_u3.eam Cu # 建模生成fcc銅納米線 lattice fcc 3.615 region box block 0 10 0 10 0 20 create_box 1 box create_atoms 1 box mass 1 63.546 # 設(shè)置區(qū)域與計(jì)算輸出 thermo 100 thermo_style custom step temp press pe ke etotal pxx pyy pzz dump 1 all custom 2000 wire_nvt.dump id type x y z reset_timestep 0 # 最小化初始構(gòu)型 min_style cg minimize 1.0e-6 1.0e-8 1000 10000 # 溫度初始化與NVT弛豫 velocity all create 300.0 1234567 loop geom fix 1 all nvt temp 300.0 300.0 0.1 run 5000 unfix 1 # 拉伸加載對(duì)盒子施加恒應(yīng)變率 fix 2 all deform 1 z erate 1.0e-4 units box fix 3 all nvt temp 300.0 300.0 0.1 dump 2 all custom 2000 wire_stretch.dump id type x y z vx vy vz run 20000 # 保存最終構(gòu)型 write_data wire_final.data3.2 逐段拆解每一個(gè)選擇背后的理由先看初始化段。我在這個(gè)案例里用了units metal因?yàn)?EAM 勢(shì)函數(shù)的參數(shù)通常以 eV 和 ? 為基準(zhǔn)配合metal單位制后面所有能量和力的數(shù)值都不需要額外換算。boundary f f p的設(shè)定是這樣的納米線在 x 和 y 方向是自由表面所以用f固定邊界讓原子在表面處天然形成真空層z 方向是拉伸方向用p周期性邊界確保納米線在長(zhǎng)度方向上可以維持連續(xù)周期性、避免端部效應(yīng)。pair_style eam對(duì)應(yīng)的是嵌入原子勢(shì)方法非常適合金屬銅。EAM 勢(shì)不僅僅計(jì)算兩兩原子間的對(duì)勢(shì)還額外考慮了每個(gè)原子嵌入在周?chē)娮用芏缺尘爸械哪芰繉?duì)描述金屬鍵合、表面重構(gòu)、位錯(cuò)產(chǎn)生都很關(guān)鍵。pair_coeff * * Cu_u3.eam Cu告訴 LAMMPS 從Cu_u3.eam文件中讀取 EAM 參數(shù)并把這個(gè)勢(shì)文件應(yīng)用到所有原子種類(lèi)上。這里的* *是通配符表示所有可能的原子對(duì)Cu則是勢(shì)文件中元素名稱(chēng)的映射。如果你有多個(gè)元素這里要逐個(gè)列出比如pair_coeff * * FeCu.eam.alloy Fe Cu。建模部分我用的是lattice fcc 3.615生成晶體格點(diǎn)。注意 3.615 是室溫附近銅的晶格常數(shù)?。用create_box 1 box創(chuàng)建盒子然后create_atoms 1 box在盒子區(qū)域內(nèi)填充原子。由于 x 和 y 方向是固定邊界這些方向原子按照 fcc 晶格排布后表面的原子自然形成裸露的納米線表面——這比人為定義一個(gè)圓柱形區(qū)域再裁剪要簡(jiǎn)單得多。mass 1 63.546是銅的原子質(zhì)量單位制為 metal 時(shí)這里用的是 g/mol。到了計(jì)算設(shè)置區(qū)塊我要特別解釋reset_timestep 0這一行。LAMMPS 的時(shí)間步是從 0 還是從接續(xù)前一個(gè) run 的步數(shù)開(kāi)始取決于前面是否執(zhí)行過(guò)run命令。最小化的運(yùn)行不會(huì)影響 timestep但為了確保后面所有輸出文件中的第二步編號(hào)對(duì)應(yīng)一致我習(xí)慣在正式動(dòng)力學(xué)之前重置一次。再看最小化。minimize 1.0e-6 1.0e-8 1000 10000的兩個(gè)容差分別控制能量和力的相對(duì)收斂。這里我特意把能量容差設(shè)為 1e-6、力容差設(shè)為 1e-8是因?yàn)榻饘袤w系存在長(zhǎng)程應(yīng)力場(chǎng)時(shí)太寬松的收斂標(biāo)準(zhǔn)會(huì)讓表面原子在后續(xù) NVT 弛豫中出現(xiàn)很強(qiáng)的初始應(yīng)力波動(dòng)。隨后velocity all create 300.0 1234567 loop geom賦予體系 300K 的初始 Maxwell-Boltzmann 速度分布。注意這里的隨機(jī)數(shù)種子 1234567如果你需要多個(gè)獨(dú)立樣本記得改掉。弛豫階段用fix 1 all nvt temp 300.0 300.0 0.1溫度上下限都設(shè)為 300K時(shí)間常數(shù) 0.1 ps。總共運(yùn)行 5000 步因?yàn)橛玫氖莡nits metal時(shí)間步默認(rèn)是 1 fs所以相當(dāng)于 5 ps 的弛豫。對(duì)一根邊長(zhǎng)僅幾納米的納米線來(lái)說(shuō)5 ps 足夠讓原子位置弛豫到合理狀態(tài)。跑完unfix 1是因?yàn)榻酉聛?lái)要進(jìn)入加載階段舊的 fix 如果不刪掉會(huì)和新的 fix 疊加造成意外約束。拉伸加載的核心是fix 2 all deform 1 z erate 1.0e-4 units box。deform表示盒子在 z 方向隨時(shí)間以恒定應(yīng)變率變形erate 1.0e-4意思是每 psmetal 單位制下應(yīng)變?cè)黾?1e-4。units box表示應(yīng)變率的單位是盒子尺寸的比值。由于 box 的 z 方向初始長(zhǎng)度是 20 個(gè)晶胞 ×3.615? ≈ 72.3?在 1 ps 內(nèi)盒子長(zhǎng)度增加約 0.00723?這個(gè)速度對(duì)于金屬納米線的準(zhǔn)靜態(tài)拉伸來(lái)說(shuō)比較合理。與此同時(shí)fix 3 all nvt保持溫度穩(wěn)定。生產(chǎn)階段共運(yùn)行 20000 步即 20 ps對(duì)應(yīng)的總應(yīng)變約為 2%。如果你想模擬更大的塑性變形把步數(shù)調(diào)大到 50000 甚至 100000 即可。最后write_data wire_final.data保存終態(tài)構(gòu)型可以用于后續(xù)的繼續(xù)模擬或者結(jié)構(gòu)分析。3.3 從拉伸結(jié)果里拿到什么提取應(yīng)力-應(yīng)變曲線的思路跑完上一段模擬后你會(huì)得到一個(gè)wire_stretch.dump軌跡文件和一個(gè)日志文件log.lammps。從這些文件里提取應(yīng)力-應(yīng)變曲線是分子模擬最常用的分析動(dòng)作。應(yīng)力數(shù)據(jù)在 log 文件里thermo_style custom step temp press pe ke etotal pxx pyy pzz這一行已經(jīng)把六個(gè)應(yīng)力分量pxx pyy pzz pxy pxz pyz 我這里只輸出了三個(gè)對(duì)角線分量記錄在案。工程應(yīng)變可以通過(guò) dump 文件里盒子的 z 方向長(zhǎng)度變化計(jì)算。如果你用了 OVITO可以直接在軌跡上讀每一幀的盒子尺寸再用公式應(yīng)變 (Lz - Lz0) / Lz0。需要特別注意LAMMPS 輸出的應(yīng)力單位制在units metal下是 bar。很多繪圖腳本直接拿來(lái)當(dāng) MPa 用差了好幾個(gè)數(shù)量級(jí)。正確的換算關(guān)系是1 bar 0.1 MPa 1e-1 MPa。所以如果看到了 pxx 數(shù)值在幾萬(wàn) bar 附近那是非常正常的——銅的理想強(qiáng)度大概在幾千 MPa對(duì)應(yīng)的 bar 數(shù)值是幾萬(wàn)。提取數(shù)據(jù)后通常還要做一步平滑處理原始應(yīng)力-應(yīng)變曲線中存在高頻的熱振動(dòng)噪聲如果不加處理直接畫(huà)圖曲線會(huì)像鋸齒一樣密密麻麻?？梢杂玫钠交椒ㄊ侨∫粋€(gè)滑動(dòng)窗口比如每 200 個(gè)數(shù)據(jù)點(diǎn)平均一次或者對(duì)整段曲線做一個(gè)低通濾波。我自己的習(xí)慣是把 log 文件用 Python 腳本讀進(jìn)來(lái)先按應(yīng)變分箱再求平均應(yīng)力這樣得到的曲線既保留了物理趨勢(shì)又干凈。4. 我踩過(guò)的in文件深坑錯(cuò)誤排查的完整思路寫(xiě) in 文件這件事出問(wèn)題幾乎是必然事件。我調(diào)到現(xiàn)在的經(jīng)驗(yàn)是——真正有價(jià)值的能力不是你永遠(yuǎn)不出錯(cuò)而是你能夠在 30 分鐘內(nèi)定位到錯(cuò)誤的根源。下面這幾個(gè)坑我全都親手踩過(guò)每一個(gè)都花了我一下午甚至一整天。4.1 “Lost atoms”不是原子真的丟了這是 LAMMPS 用戶(hù)最經(jīng)典的報(bào)錯(cuò)ERROR: Lost atoms。報(bào)錯(cuò)信息很?chē)樔说鎸?shí)原因往往是體系中出現(xiàn)了一個(gè)原子被其他原子推開(kāi)到極其遙遠(yuǎn)的距離——它沒(méi)有真的憑空消失只是飛出了近鄰列表能夠追蹤的范圍。最常見(jiàn)的觸發(fā)場(chǎng)景有三類(lèi)第一初始構(gòu)型有原子重疊導(dǎo)致局部力過(guò)大某個(gè)原子被瞬間彈出。第二時(shí)間步長(zhǎng)過(guò)大。金屬體系建議時(shí)間步不超過(guò) 1 fsunits metal默認(rèn) 1 fs 通常沒(méi)問(wèn)題但如果溫度很高或者勢(shì)函數(shù)特別硬可能需要降到 0.5 fs。第三fix nvt的耦合時(shí)間常數(shù)過(guò)短速度更新過(guò)于劇烈導(dǎo)致溫度瞬間波動(dòng)過(guò)大。排查方法也有固定套路。先把dump輸出頻率調(diào)高比如每 10 步輸出一幀然后可視化觀察到底是哪一步哪個(gè)原子開(kāi)始被彈出。如果在很早期比如前 100 步就出現(xiàn)那基本是初始構(gòu)型問(wèn)題回到建模階段檢查坐標(biāo)如果在后期才出現(xiàn)優(yōu)先檢查時(shí)間步長(zhǎng)和勢(shì)函數(shù)參數(shù)。還有一個(gè)通用技巧把thermo輸出調(diào)密一點(diǎn)觀察溫度的變化趨勢(shì)。如果溫度在某個(gè)時(shí)間點(diǎn)突然出現(xiàn)脈沖式的尖峰那一瞬間就是原子被彈出的時(shí)刻。4.2 溫度失控NVE 體系熱量為什么越積越多有一次我模擬一個(gè)高分子體系跑了大概 200 ps 之后溫度一路飆升從 300 K 漲到了 420 K。查了很久才發(fā)現(xiàn)原因是體系內(nèi)部一直在產(chǎn)生熱量但 NVE 系綜下這些熱量無(wú)路可走只能轉(zhuǎn)化為原子動(dòng)能的增加。產(chǎn)生熱量的物理過(guò)程可能來(lái)自粘性形變、化學(xué)反應(yīng)雖然我還沒(méi)開(kāi) bond breaking、或非物理的數(shù)值耗散。對(duì)于我的情況罪魁禍?zhǔn)灼鋵?shí)是時(shí)間步長(zhǎng)太大——LAMMPS 的 Verlet 積分在時(shí)間步過(guò)大時(shí)會(huì)產(chǎn)生額外的數(shù)值能量漂移宏觀表現(xiàn)就是體系“變熱”。解決辦法有三種按可操作性強(qiáng)弱排序第一把時(shí)間步長(zhǎng)減半觀察溫度是否恢復(fù)平穩(wěn)第二確認(rèn)體系的初始溫度是不是過(guò)高或過(guò)低如果在低溫下直接跑體系內(nèi)部應(yīng)力會(huì)導(dǎo)致局部能量集中第三生產(chǎn)模擬階段改用fix nvt而不是fix nve讓溫度被恒溫器控制住。當(dāng)然如果你要嚴(yán)格研究微正則系綜nve 是必須的但前提是已經(jīng)通過(guò) nvt 弛豫到了一個(gè)穩(wěn)定的初態(tài)并且時(shí)間步長(zhǎng)足夠小。4.3 數(shù)據(jù)文件字段錯(cuò)位最隱蔽的靜默錯(cuò)誤前面提到過(guò)一個(gè)經(jīng)典坑atom_style與 data 文件字段不匹配。如果read_data之后 LAMMPS 沒(méi)有報(bào)錯(cuò)但計(jì)算結(jié)果明顯異常比如密度憑空少了 20%、勢(shì)能曲線形態(tài)完全不對(duì)、原子分布出現(xiàn)奇怪的帶狀那就要高度懷疑 fields 錯(cuò)位問(wèn)題。最常見(jiàn)的場(chǎng)景是你的 data 文件是某個(gè)舊版本軟件導(dǎo)出的Atom 區(qū)域包含 id type x y z 五列但你的 in 文件聲明了atom_style charge于是 LAMMPS 會(huì)認(rèn)為前五列分別是 id type q x y原本的原子坐標(biāo)被當(dāng)成了電荷量坐標(biāo)位置整個(gè)錯(cuò)亂。這種錯(cuò)誤在早期 bug 排查階段很容易被忽略因?yàn)閞ead_data通常只會(huì)在文件格式完全不匹配時(shí)報(bào)錯(cuò)而字段語(yǔ)義錯(cuò)位往往是靜默的。檢查方法也很簡(jiǎn)單建一個(gè)小體系讀入后用write_data重新寫(xiě)出來(lái)人工比對(duì)坐標(biāo)是否保持原值。如果坐標(biāo)發(fā)生了系統(tǒng)性偏移那就是原子數(shù)據(jù)字段解釋有誤。4.4 體積不守恒NPT 弛豫后為什么盒子收縮到離譜另一個(gè)高頻坑是模擬液體或高分子跑 NPT 時(shí)盒子體積越縮越小直到壓強(qiáng)歸零但體系明顯還剩下大量空隙。很多人第一反應(yīng)是壓力耦合常數(shù)設(shè)錯(cuò)了。實(shí)際上這通常是初始構(gòu)型的“真空”導(dǎo)致的。如果你用內(nèi)部生成的方法做了低密度初始構(gòu)型然后在 NPT 下弛豫盒子為了達(dá)到設(shè)定的 1 atm 壓力會(huì)不斷壓縮把空隙擠掉。這個(gè)過(guò)程在物理上是正確的但如果初始密度太低或者 NPT 溫控/壓控時(shí)間常數(shù)設(shè)置不合理盒子可能在很短時(shí)間被壓成一個(gè)極度畸形的形狀后續(xù)一切數(shù)據(jù)都不可用。解決辦法是初始構(gòu)型盡量接近目標(biāo)密度。如果是因?yàn)榻９ぞ弋a(chǎn)生的真空可以用delete_atoms overlap先刪除掉重疊原子再小心地跑了 MVT 壓縮到目標(biāo)密度最后繼續(xù)生產(chǎn)。另外壓力耦合時(shí)間常數(shù)比如tparam 1000.0過(guò)短是常見(jiàn)誤區(qū)一般建議設(shè)定在 1000 fs 的量級(jí)即 1 ps而不要給到幾十 fs。5. 從“能跑”到“跑得好”in文件的中級(jí)進(jìn)階技巧如果你已經(jīng)能把一份 in 文件跑通并且能熟練改參數(shù)那么恭喜你接下來(lái)可以進(jìn)入讓文件本身變得“可維護(hù)、可復(fù)用、可擴(kuò)展”的階段。這個(gè)階段的目標(biāo)不再是臨時(shí)湊一份能用的文件而是建立一個(gè)可以承載各種實(shí)驗(yàn)變化的工程框架。5.1 用include和變量把in文件拆成模塊一個(gè)常見(jiàn)的誤區(qū)是所有內(nèi)容都堆在一個(gè) in 文件里換體系就整體改一遍。等到項(xiàng)目多了之后你會(huì)發(fā)現(xiàn)改動(dòng)越來(lái)越容易出錯(cuò)因?yàn)橥粋€(gè)參數(shù)可能在多個(gè)位置出現(xiàn)而你不可能每次都記得全部同步。更好的做法是模塊化拆分。把一套模擬拆成幾個(gè)文件# 主文件 main.in units metal dimension 3 boundary f f p atom_style atomic # 建模部分單獨(dú)放一個(gè)文件 include model_build.in # 勢(shì)函數(shù)單獨(dú)放一個(gè)文件 include potential.in # 計(jì)算與輸出單獨(dú)放一個(gè)文件 include output_settings.in # 弛豫生產(chǎn)運(yùn)行 include run_production.in每個(gè)子文件內(nèi)部還可以通過(guò)variable命令定義參數(shù)方便批量修改。例如variable temp equal 300.0 variable timestep equal 0.001 variable total_steps equal 50000然后在命令里直接用${temp}來(lái)引用變量。這樣每次調(diào)整溫度只需要改一個(gè)地方而不是 grep 整份文件。如果要在同一服務(wù)器上批量跑不同溫度的模擬還可以配合 shell 腳本生成多個(gè)替換變量的變體文件——這在做相變溫度掃描、應(yīng)變率掃描時(shí)簡(jiǎn)直救命。5.2 用label和jump實(shí)現(xiàn)簡(jiǎn)單的循環(huán)與控制邏輯LAMMPS 雖然是一門(mén)“腳本語(yǔ)言”但它的控制流能力非常有限。不過(guò)labeljump的組合可以模擬一個(gè)簡(jiǎn)單的循環(huán)。典型用途是在同一個(gè) in 文件中把一個(gè)溫度范圍內(nèi)所有溫度點(diǎn)依次跑一遍而不用生成多個(gè)文件。比如做一個(gè)從 300K 到 600K 升溫掃描variable t index 300 350 400 450 500 550 600 label loop fix 1 all nvt temp ${t} ${t} 0.1 run 10000 unfix 1 variable t next jump main.in loop這里variable t index定義了一個(gè)可以取多值的索引變量label loop聲明循環(huán)起點(diǎn)末尾variable t next讓 t 取下一個(gè)值然后jump main.in loop跳回循環(huán)體開(kāi)頭。執(zhí)行到 t 的所有取值用完后循環(huán)自動(dòng)結(jié)束。這個(gè)技巧對(duì)于自動(dòng)化參數(shù)掃描非常高效而且不需要 Python 腳本介入。注意jump的語(yǔ)法在多個(gè) in 文件之間跳轉(zhuǎn)時(shí)要注意文件名匹配。如果你直接把主文件命名為 main.in那么jump main.in loop沒(méi)問(wèn)題如果文件被改名這里也要同步改否則 LAMMPS 會(huì)找不到文件。5.3 讓后處理管線化從in文件到分析的“一條龍”思路模擬本身只占了工作量的 50%后處理數(shù)據(jù)分析是另一半。我自己的項(xiàng)目習(xí)慣是每個(gè)模擬算例都配套一個(gè)分析腳本目錄in 文件之外還有提取力-應(yīng)變曲線的 Python 腳本、計(jì)算徑向分布函數(shù)的腳本、以及可視化用的 OVITO 狀態(tài)文件。以拉伸模擬為例我通常會(huì)在提交模擬算例的同時(shí)寫(xiě)好提取曲線的腳本。關(guān)鍵代碼思路是讀入 log 文件篩選出Step、Pxx列計(jì)算應(yīng)變?nèi)缓笞龌瑒?dòng)窗口平滑。如果你算的是多個(gè)應(yīng)變率下的響應(yīng)腳本里再套一層文件循環(huán)最后把所有曲線畫(huà)在同一張圖上。這里延伸一點(diǎn)很多人忽略了一個(gè)細(xì)節(jié)——thermo_style custom輸出應(yīng)力時(shí)LAMMPS 輸出的壓強(qiáng)是體系瞬時(shí)壓強(qiáng)的平均包含了動(dòng)能貢獻(xiàn)和維里貢獻(xiàn)。在做應(yīng)力-應(yīng)變分析時(shí)建議用維里應(yīng)力部分pzz而不是總壓強(qiáng)的直接值尤其在高速加載條件下動(dòng)能項(xiàng)的瞬時(shí)波動(dòng)會(huì)掩蓋材料的真實(shí)應(yīng)力響應(yīng)。LAMMPS 輸出的pzz已經(jīng)是包含動(dòng)能項(xiàng)的完整壓強(qiáng)不需要額外去除但在解讀時(shí)要知道它不等于材料內(nèi)部經(jīng)驗(yàn)意義上的“工程應(yīng)力”兩者差了約一個(gè)負(fù)號(hào)或者單位換算精度需要和你的分析方式匹配。6. 寫(xiě)在最后把in文件當(dāng)程序?qū)懚皇钱?dāng)配置改如果你讀到這里說(shuō)明你已經(jīng)開(kāi)始把 in 文件當(dāng)作一個(gè)需要認(rèn)真設(shè)計(jì)的計(jì)算流程來(lái)看待了。我個(gè)人這幾年最大的體會(huì)就是真正的高手不靠記住命令而是靠建立一套“先設(shè)計(jì)再實(shí)現(xiàn)”的工作習(xí)慣。拿到一個(gè)研究課題第一步不是打開(kāi)文本編輯器開(kāi)始寫(xiě)命令而是先在紙上列出研究問(wèn)題需要哪些輸出量、體系需要什么初始構(gòu)型、采用什么系綜和勢(shì)函數(shù)、生產(chǎn)階段跑多長(zhǎng)、掃描哪些參數(shù)。想清楚了這些問(wèn)題再打開(kāi) in 文件你會(huì)發(fā)現(xiàn)每一行命令的落點(diǎn)都清清楚楚。最后再分享一個(gè)小技巧如果你第一次跑一套全新體系的模擬不要直接上生產(chǎn)規(guī)模。先用最小體系比如 1000 個(gè)原子、短跑 5000 步驗(yàn)證 in 文件沒(méi)有報(bào)錯(cuò)、能量趨勢(shì)合理、軌跡文件能正常打開(kāi)然后再放大體系、延長(zhǎng)周期。這個(gè)“兩步走”看起來(lái)多花了一個(gè)小時(shí)實(shí)際上能幫你省掉后面花幾天時(shí)間排查一個(gè)只有在大體系上才會(huì)顯現(xiàn)的隱性錯(cuò)誤。如果你在調(diào)試過(guò)程中遇到了無(wú)法解決的問(wèn)題把手頭的 in 文件和 log 文件保留好逐行比對(duì)官方文檔的關(guān)鍵詞絕大部分問(wèn)題都能在文檔的 “Restrictions” 段落找到答案。祝你們都能跑出漂亮的數(shù)據(jù)。