相较于之前做的单晶硅纳米磨削，聚合物的建模要复杂的多，用lammps自带密令通常无法解决。常用的建模方法有用python/C/MATLAB代码生成聚合物data文件，不过就算建模个简单的聚合物往往也需要写上上百行代码。moltemplate作为专门生成lammps的软件在代码量上比前者少很多。使用可视化软件material studio的话能够省去coding的过程。

一、模型建立


文中采用的联合原子模型的方法进行建模，将聚乙烯分子链中的CH2整体看成是一个原子，这样的话能够将模型中的原子、键角、二面角等类型都简化为一个，极大的减少了计算量。（节省了超算上跑仿真花的money。）
本次建模采用的是material studio，步骤如下：

  • 绘制乙烯单体，采用clean结构优化。
  • build>repeat unit。选择两个碳上各一个H作为聚合的头尾。因为采用的联合原子模型，所以这里需要删除其余未被选中的H。
  • build>hopolymer。建立单条链长为50的聚乙烯分子链，并删除两端的H。
  • 建立amorphous cell，包含十个之前建立的聚乙烯分子链，力场选择cvff，密度1g/cm3。最终建立的模型如下：

image.png

  • 导出car和mdf文件，打开终端输入以下密令生成data文件:

    1. ./msi2lmp PE - i -class 1 -frc cvff >data.pe

  • 导出的data文件部分内容： ```python LAMMPS data file. msi2lmp v3.9.8 / 06 Oct 2016 / CGCMM for PE

    1000 atoms 990 bonds 980 angles 970 dihedrals

    1. 0 impropers

    1 atom types 1 bond types 1 angle types 1 dihedral types

    3.286577762 63.286577762 xlo xhi 5.871187987 65.871187987 ylo yhi 5.606698939 65.606698939 zlo zhi

Masses

1 12.011150 # c

Pair Coeffs # lj/cut/coul/long

1 0.1599999990 3.4745050026 # c

Bond Coeffs # harmonic

1 322.7158 1.5260 # c-c

Angle Coeffs # harmonic

1 46.6000 110.5000 # c-c-c

Dihedral Coeffs # harmonic

1 1.4225 1 3 # c-c-c-c

Atoms # full

  1.   1      1   1  0.000000    12.856758660    13.246611987    23.346400350   0   0  -1 # c
  2.   2      1   1  0.000000    10.871831278    15.621789577    21.922500446   0   0  -1 # c
  3.   3      2   1  0.000000    15.801140500    11.904915192    27.872642417   0   0  -1 # c
  4.   4      2   1  0.000000    13.866863126    14.306208863    26.423071451   0   0  -1 # c
  5.   5      3   1  0.000000    13.038691021    11.402970592    32.676112667   0   0  -1 # c
  6.   6      3   1  0.000000    15.946417178    12.440044884    31.234411673   0   0  -1 # c
  7.   7      4   1  0.000000    10.852350394    10.817753706    37.758812525   0   0  -1 # c
  8.   8      4   1  0.000000    12.719850433    12.838692910    35.749547837   0   0  -1 # c
  9.   9      5   1  0.000000     8.568337190    10.414035534    42.816153454   0   0  -1 # c
  10.  10      5   1  0.000000     9.836755987    12.567439705    40.500368386   0   0  -1 # c
  1. <br />
  2.  注意因为采用的联合原子模型，需要将data文件中的mass由12（C原子）改成14（CH2）。data文件中势函数部分之间删除，之后再in文件中按照文献进行设置。
  3. <a name="jJ7PD"></a>
  4. ## 二、in文件部分
  5. 1.初始化和读取data文件
  6. ```python
  7. #-------聚乙烯拉伸in文件-------#
  9. #-------初始化设定-------#
  11. units real #单位制
  12. boundary p p p #周期性边界
  13. atom_style full #分子+电荷
  15. #--------模型读取-------#
  17. read_data PE.data #读取MS建立的模型

2.势函数设置

  1. #-------势函数设置  文献：doi:10.1016/j.polymer.2010.10.009-------#
  3. bond_style harmonic #键势势类型
  4. bond_coeff 1 350 1.53 #CH2-CH2键 Kb=350 r0=1.53
  6. angle_style harmonic #键角势类型
  7. angle_coeff 1 60 109.5 #CH2-CH2键角 Kθ=60 键角109.5度
  9. dihedral_style multi/harmonic #二面角势类型
  10. dihedral_coeff 1 1.736 -4.490 0.776 6.990 0 #CH2-CH2二面角（四种） C0=1.736 C1=-4.490 C2=0.776 C3=6.990 C4=0
  12. pair_style lj/cut 10 #对相互作用势
  13. pair_coeff 1 1 0.122 4.01 #CH2 CH2 ε=0.112 σ=4.0

此处的参数设置完全按照文献中所给：
image.png
3.neighbor参数

  1. #-------设置邻居参数-------#
  3. neighbor 0.4 bin #截断距离0.4
  4. neigh_modify every 10 one 10000 #每十步重新构建列表，一个原子最大邻居数为10000

4.驰豫过程
按照文中所述的四阶段驰豫：500K下NVT驰豫，500K下NPT驰豫，500K-期望温度NPT驰豫，期望温度NPT驰豫。

  1. #-------平衡态分子模拟-------#
  3. #-------nvt条件下弛豫-------#
  5. dump 1 all custom 100 1.xyz id type x y z #输出nvt弛豫下的原子序号，种类和坐标到nvt.xyz。
  7. fix 1 all nvt temp 500 500 50 #500k下弛豫
  9. thermo 100 #100步输出一次
  11. reset_timestep 0 #起始步数归零
  13. timestep 1
  15. run 10000 #弛豫10000步
  17. unfix 1 #取消 fix 1
  18. undump 1 #取消 dunmp 1
  20. #-------npt条件下弛豫-------#
  22. fix 2 all npt temp 500 500 50 iso 0 0 1000 #npt弛豫。起始终止温度500，波动温度50。起始终止压力0，波动压力1000。
  24. dump 1 all custom 1000 2.xyz id type x y z #输出npt弛豫下所有原子的xyz坐标。
  26. thermo 100 #100步输出一次
  28. reset_timestep 0 #起始步数归零
  30. run 50000 #跑5000步
  32. unfix 2 #取消 fix 1
  33. undump 1 #取消 dump 1
  35. #-------npt条件下弛豫-------#
  37. fix 3 all npt temp 500 200 50 iso 0 0 1000 #npt弛豫。起始终止温度500-200，波动温度50。起始终止压力0，波动压力1000。
  39. dump 1 all custom 1000 3.xyz id type x y z #输出npt弛豫下所有原子的xyz坐标。
  41. thermo 100 #100步输出一次
  43. reset_timestep 0 #起始步数归零
  45. run 50000 #跑5000步
  47. unfix 3 #取消 fix 1
  48. undump 1 #取消 dump 1
  50. #-------npt条件下弛豫-------#
  52. fix 4 all npt temp 200 200 50 iso 0 0 1000 #npt弛豫。起始终止温度200，波动温度50。起始终止压力0，波动压力1000。
  54. dump 1 all custom 1000 4.xyz id type x y z #输出npt弛豫下所有原子的xyz坐标。
  56. thermo 100 #100步输出一次
  58. reset_timestep 0 #起始步数归零
  60. run 50000 #跑5000步
  62. unfix 4 #取消 fix 1
  63. undump 1 #取消 dump 1

驰豫过程的ovito轨迹：
lammps学习（二）联合原子模型聚乙烯拉伸 - 图3
5.拉伸过程

  1. #-------恒温恒压拉伸-------#
  3. variable temp equal "lx" #设定参数temp等于体系长度
  4. variable L0 equal ${temp} #设定L0等于temp固定值
  5. variable strain equal "(lx - v_L0)/v_L0" #计算应变大小
  7. variable p1 equal "v_strain" #应变
  8. variable p2 equal "-pxx/10000*1.01325" #x方向应力（单位转换）
  9. variable p3 equal "-pyy/10000*1.01325" #y方向应力（单位转换）
  10. variable p4 equal "-pzz/10000*1.01325" #z方向应力（单位转换）
  11. variable step0 equal step
  13. reset_timestep 0  
  15. dump 1 all custom 1000 dump.pe id type x y z #输出拉伸状态下原子坐标
  17. fix 1 all npt temp 200 200 50 y 0 0 1000 z 0 0 1000 drag 2 #恒温恒压，摩擦系数2。
  18. fix 2 all deform 1 x erate 1e-5 units box remap x #x方向拉伸，速率1e-5，box每步变形，原子重新映射到变形后的box。
  19. fix 3 all ave/time 10 5 100 v_p1 v_p2 v_p3 v_p4 file stress.txt #平均化采样
  21. thermo_style custom step temp v_strain pxx pyy pzz lx ly lz epair ebond eangle edihed #自定义热力学输出
  22. thermo 100 #每100步输出一次
  24. reset_timestep 0 #起始步数归零
  26. run 200000 #跑20000步

拉伸过程的ovito轨迹：
lammps学习（二）联合原子模型聚乙烯拉伸 - 图4

三、数据分析

利用输出文件stress.txt中的v_p1和v_p2降噪绘制应力应变曲线：
image.png
各种相互作用能随应变的变化：
image.png
同理可以得到文中不同链长、不同温度下的拉伸情况。