LAMMPS（Large-scale Atomic/Molecular Massively Parallel Simulator）是一種經典的分子動力學仿真代碼，專注於材料建模。它旨在在並行計算機上高效運行，並且易於擴展和修改。LAMMPS 最初由美國能源部機構桑迪亞國家實驗室開發，現在包括來自許多機構的許多研究小組和個人的貢獻。LAMMPS 的大部分資金來自美國能源部（DOE）。LAMMPS 是根據 GNU 公共許可證版本 2（GPLv2）的條款分發的開源軟件。
在本次教程中，我們通過改變材料的晶格常數，實現模擬對施加材料單軸應變的情況，後續再計算並繪製材料的應變應力曲線。通過本教程的學習您將學會：

1. 使用 npt 讓系統結構達到弛豫
2. 使用 fix 結合 variable 命令改變晶格常數
3. 使用 compute 命令計算系統應力

教程鏈接：https://go.openbayes.com/kuLKA

該教程將在雲平台 http://OpenBayes.com 上進行演示，使用下方邀請鏈接註冊即可獲得 4 小時 RTX 5090 免費使用時長：
https://openbayes.com/console/signup?r=Dennis9801_1ohB

一、輸入文件説明

本次教程將會用到以下輸入文件：

• Al99.eam.alloy 材料的 eam 勢
• in.txt 輸入文件
• pl.py 繪製圖像腳本

系統初始化與基礎設置

units metal               # 採用金屬單位制（長度：Å，時間：ps，質量：g/mol，能量：eV等）
dimension 3               # 三維空間模擬
boundary p p p            # 三個方向均為週期性邊界條件（模擬無限大晶體）
atom_style atomic         # 原子類型為「基本原子」，僅包含位置和類型信息
variable latparam equal 4.05  # 定義鋁的晶格常數為 4.05Å（FCC 鋁的典型值）

原子結構構建

lattice  fcc ${latparam} orient x 1 0 0 orient y 0 1 0 orient z 0 0 1  # 定義 FCC 晶格
# 晶格常數為 ${latparam}，晶向沿 x[100]、y[010]、z[001]（標準笛卡爾座標系）

region whole block 0 10 0 10 0 10  # 定義模擬區域：一個立方體，在 x、y、z 方向各包含 10 個晶格常數長度
create_box 1 whole                 # 創建容納原子的「盒子」，僅包含 1 種原子類型，邊界為上述 region 定義的範圍
create_atoms 1 region whole        # 在盒子內填充類型為 1 的原子，按 FCC 晶格排列

原子間相互作用勢設置

pair_style eam/alloy  # 採用嵌入原子法（EAM）合金勢描述原子間相互作用（適合金屬體系）
pair_coeff * * Al99.eam.alloy Al  # 配置勢函數參數：
                                  # - 第一個*：所有原子類型（此處僅 1 種）
                                  # - 第二個*：所有原子類型的相互作用
                                  # - Al99.eam.alloy：勢函數文件（預訓練的鋁勢）
                                  # - Al：指定原子類型 1 對應元素為鋁

計算量定義（用於後續分析）

compute csym all centro/atom fcc  # 計算每個原子的中心對稱參數（用於識別 FCC 結構完整性，缺陷處會偏離）
compute peratom all pe/atom       # 計算每個原子的勢能（用於分析能量分佈）

平衡模擬階段（消除初始應力，達到穩定狀態）

reset_timestep 0  # 重置時間步計數器為 0（平衡階段作為新起點）
timestep 0.001    # 時間步長為 0.001ps（即 1fs，金屬模擬常用精度）

#初始化原子速度：温度 300K（室温），隨機種子 12345，消除整體動量（mom yes）和旋轉（rot no）
velocity all create 300 12345 mom yes rot no
#施加 NPT 系綜（恆温、恆壓、恆粒子數）進行平衡：
#- temp 300 300 1：温度維持 300K，熱浴阻尼係數 1（單位：ps）
#- iso 0 0 1：各向同性壓力控制，目標壓力 0bar，壓力阻尼係數 1（單位：ps）
#- drag 1：原子運動阻尼係數（用於穩定系綜）
fix 1 all npt temp 300 300 1 iso 0 0 1 drag 1 
#熱力學輸出設置：
thermo 1000  #每 1000 步輸出一次熱力學信息
thermo_style custom step lx ly lz press pxx pyy pzz pe temp  #輸出內容：
                                                             #步數、盒子尺寸（x/y/z）、總壓力、壓力張量分量、勢能、温度
run 20000  #運行 20000 步平衡模擬（總時間20000×0.001ps=20ps）
unfix 1    #平衡結束，移除 NPT 約束
#記錄平衡後的 x 方向盒子長度（用於後續計算應變）
variable tmp equal "lx"  #臨時變量存儲當前 x 方向長度
variable L0 equal ${tmp}  #定義 L0 為初始長度（拉伸前的平衡長度）
print "Initial Length, L0: ${L0}"  #輸出初始長度，便於驗證

單軸拉伸變形階段

reset_timestep 0  #重置時間步計數器為 0（拉伸階段作為新起點）

#重新施加 NPT 系綜，但僅控制温度和橫向壓力：
#- temp 300 300 1：維持温度300K
#- y 0 0 1 z 0 0 1：y 和 z 方向壓力保持 0bar（橫向自由，模擬單軸拉伸）
fix 1 all npt temp 300 300 1 y 0 0 1 z 0 0 1 drag 1
#定義拉伸速率：
variable srate equal 1.0e10  #應變速率為 1e10 /s（分子動力學常用高應變速率，遠高於實驗）
variable srate1 equal "v_srate / 1.0e12"  #轉換為LAMMPS單位（應變速率單位：ps⁻¹，1e10/s = 1e-2 ps⁻¹）
#施加 x 方向拉伸變形：
#- deform 1 x：沿 x 方向變形，變形組為 1（所有原子）
#- erate ${srate1}：按上述應變速率拉伸
#- units box：變形基於盒子尺寸
#- remap x：原子座標隨盒子拉伸同步更新（避免原子「跑出」盒子）
fix 2 all deform 1 x erate ${srate1} units box remap x

拉伸過程數據輸出

#定義應變和應力變量：
variable strain equal "(lx - v_L0)/v_L0"  #工程應變 =（當前x長度 - 初始長度）/ 初始長度
variable p1 equal "v_strain"              #p1 對應工程應變（無量綱）
variable p2 equal "-pxx/10000"            #x 方向應力（GPa）：pxx 為 LAMMPS 內置壓力張量（單位 bar），負號轉為拉應力，除以 10000 轉換為 GPa
variable p3 equal "-pyy/10000"            #y 方向應力（GPa）
variable p4 equal "-pzz/10000"            #z 方向應力（GPa）

#輸出應變-應力數據到文件：
#- 每 100 步輸出一次 p1（應變）、p2（x 應力）、p3（y 應力）、p4（z 應力）
#- 保存到 Al_tens_100.def1.txt，不在屏幕顯示
fix def1 all print 100 "${p1} ${p2} ${p3} ${p4}" file Al_tens_100.def1.txt screen no
#輸出原子軌跡文件（用於可視化）：
#- 每 100 步保存一次，文件名為 md.lammpstrj
#- 包含原子 ID、類型、x/y/z 座標
dump trace all custom 100 md.lammpstrj id type x y z
#熱力學信息輸出設置：
thermo 1000  #每 1000 步輸出一次
thermo_style custom step v_strain temp v_p2 v_p3 v_p4 ke pe press  #輸出步數、應變、温度、三方向應力、動能、勢能、總壓力
run 30000  #拉伸階段運行 30000 步（總時間 30ps，累積應變約0.3）

模擬結束

print "All done"  # 輸出完成信息，提示模擬結束

二、操作步驟

1. 克隆並啓動容器

登錄 OpenBayes.com ，在「公共教程」頁面，選擇「LAMMPS：以單晶鋁為例，模擬材料單軸拉伸」教程。

頁面跳轉後，點擊右上角「克隆」，將該教程克隆至自己的容器中。

選擇「NVIDIA GeForce RTX 4090」以及「lammps」鏡像，OpenBayes 平台提供了 4 種計費方式，大家可以按照需求選擇「按量付費」或「包日/周/月」，點擊「繼續執行」。可以使用文章開頭的邀請鏈接，獲得 RTX 4090 使用時長！

待系統分配好資源，當狀態變為「運行中」後，點擊「打開工作空間」。

2. 運行 lammps

打開「終端」，輸入以下命令運行 lammps。

mpirun --allow-run-as-root -np 2 lmp -sf gpu -pk gpu 1 < in.txt 

3. 文件輸出

待模型處理完成後，會得到兩個輸出文件：
Al_tens_100.def1.txt：應變 - 應力曲線數據，可用於分析彈性模量、屈服強度等力學性能。
md.lammpstrj：原子軌跡文件，可通過 VMD、Ovito 等軟件可視化拉伸過程中的晶體結構演變

4. 繪製應變-應力圖像

輸入以下命令即可得到繪製的應變-應力曲線圖。

python pl.py

此外，我們還可以將「md.lammpstrj」文件下載到本地並載入到 VMD、Ovito 等軟件，得到實時模型。