LAMMPS Simulation — 聚乙烯玻璃化转变温度 (Tg) 计算

项目结构

├── Tg_simulation/              # 聚乙烯 (PE) Tg 模拟
│   ├── build_pe_melt.data      # PE 熔体初始构型 (4040 atoms)
│   ├── build_pe_melt.esh       # 构型构建脚本 (可能需要esh2input)
│   └── palce_polymer.lmp       # LAMMPS 输入脚本 (压缩 → 弛豫 → 逐级降温)
├── analysis.py                 # Tg 分析脚本 (分段线性拟合)
├── log.lammps                  # LAMMPS 输出日志 (含所有 thermo 数据)
├── density.dat                 # 逐温度点的密度记录
├── dump_pe_traj.lammpstrj      # 轨迹文件 (用于可视化)
├── Tg_analysis.png             # 分析结果图
│
├── Nanosheared_electrolyte/    # 纳米剪切电解液模拟
├── "PEG in water"/            # 水溶液中PEG模拟
├── Water_adsorption_in_silica/ # 二氧化硅吸水模拟
├── calculate_free_energy/      # 自由能计算
├── reactive_molecular_dynamics/ # 反应分子动力学
└── reactive_silicon_dioxide/   # 二氧化硅反应模拟

---

Tg 模拟流程

力场与模型

参数	值
力场	TraPPE-UA (联合原子)
原子类型	CH₂ (type 1), CH₃ (type 2)
原子数	4040
系综	NPT
压力	1 atm (降温阶段)

三步模拟协议 (palce_polymer.lmp)

1. 压缩阶段: NPT @ 500K, 100 atm, 50000 steps
   → 将盒子从 ~2.7e6 压缩至 ~1.3e5 ų

2. 弛豫阶段: NPT @ 500K, 1 atm, 100000 steps
   → 体系在目标密度下平衡

3. 降温阶段: NPT @ 1 atm, 逐级降温
   400K → 350K → 300K → 250K → 200K → 150K → 125K → 100K
   每温度 100000 steps (~100 ps)，每 1000 步输出密度

关键参数

• 时间步长: 1 fs (LAMMPS real units)
• 每温度点时长: 100 ps
• 等效冷却速率: 每 50K 降温 / 100 ps ≈ 500 K/ns
• 截断半径: 14.0 Å

---

Tg 分析 (analysis.py)

方法

经典的密度-温度曲线分段线性拟合法：

1. 从 log.lammps 提取降温阶段 thermo 数据 (Step Temp Density Volume)
2. 去重（每温度末态 = 下温度初态，取平均）
3. 自动遍历所有可能分割点，选择两段线性回归 R² 和最大的方案
4. 高温段（橡胶态）和低温段（玻璃态）分别做线性回归
5. Tg = 两条拟合直线的交点

使用方法

python analysis.py

输出: 终端显示拟合参数 + Tg_analysis.png (密度-T 和比容-T 双图)

分析结果

数据点 (已排除 500K 预平衡点)

温度 (K)	密度 (g/cm³)	相态
398.4	0.7759	橡胶态
355.1	0.8133	橡胶态
299.8	0.8355	橡胶态
251.8	0.8577	橡胶态
———— Tg ————	————	————
205.5	0.8716	玻璃态
151.3	0.8807	玻璃态
124.7	0.8852	玻璃态
97.8	0.8904	玻璃态

计算结果

方法	高温段斜率	低温段斜率	Tg
密度法	-1.73×10⁻⁴	-5.37×10⁻⁴	227–244 K
比容法	—	—	~247 K

注: 不同分割位置导致 Tg 在 201–332 K 间变化，主要原因是数据量有限且冷却速率太快，密度-T 曲线接近线性。

---

关于 Tg 的重要讨论

为什么 MD 算出的 Tg 比实验值高？

因素	实验	本模拟	影响
冷却速率	~1 K/min (1.7×10⁻⁵ K/ns)	~500 K/ns	速率差 ~10¹⁰ 倍
力场	真实化学	TraPPE-UA	联合原子降低灵活性
平衡时间	分钟-小时	100 ps/温度	远未充分弛豫

实验 PE 的 Tg ≈ 150 K，但快冷 MD 模拟得到 200–250 K 是完全正常的 —— 冷却速率每降低 10 倍，Tg 下降约 3–10 K。

改进方向

1. 降低冷却速率: 每温度点跑 500k–5M 步 (参见外推法)
2. 加密温度点: 在 150–250 K 区间每 10–15 K 取点
3. 多速率外推: 跑 3–4 组不同冷却速率，外推至实验速率

---

参考文献

• TraPPE-UA force field: J. Phys. Chem. B, 1998, 102, 2569–2577
• Glass transition in polymers: Macromolecules, 1995, 28, 500–510
• Cooling rate dependence of Tg in MD: J. Chem. Phys., 2013, 138, 12A508