lammps学习（二）联合原子模型聚乙烯拉伸

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

本次的参考文献为：doi:10.1016/j.polymer.2010.10.009 。

一、模型建立

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

本次建模采用的是material studio，步骤如下：

1. 绘制乙烯单体，采用clean结构优化。
2. build>repeat unit。选择两个碳上各一个H作为聚合的头尾。因为采用的联合原子模型，所以这里需要删除其余未被选中的H。
3. build>hopolymer。建立单条链长为50的聚乙烯分子链，并删除两端的H。
4. 建立amorphous cell，包含十个之前建立的聚乙烯分子链，力场选择cvff，密度1g/cm3。最终建立的模型如下：
5. 导出car和mdf文件，打开终端输入以下密令生成data文件:

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

 可执行文件msi2lmp下载链接：msi2lmp.exe_msi2lmp.exe-教育文档类资源-CSDN下载

     6.导出的data文件部分内容：

LAMMPS data file. msi2lmp v3.9.8 / 06 Oct 2016 / CGCMM for PE

   1000 atoms
    990 bonds
    980 angles
    970 dihedrals
      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  0.000000    12.856758660    13.246611987    23.346400350   0   0  -1 # c
      2      1   1  0.000000    10.871831278    15.621789577    21.922500446   0   0  -1 # c
      3      2   1  0.000000    15.801140500    11.904915192    27.872642417   0   0  -1 # c
      4      2   1  0.000000    13.866863126    14.306208863    26.423071451   0   0  -1 # c
      5      3   1  0.000000    13.038691021    11.402970592    32.676112667   0   0  -1 # c
      6      3   1  0.000000    15.946417178    12.440044884    31.234411673   0   0  -1 # c
      7      4   1  0.000000    10.852350394    10.817753706    37.758812525   0   0  -1 # c
      8      4   1  0.000000    12.719850433    12.838692910    35.749547837   0   0  -1 # c
      9      5   1  0.000000     8.568337190    10.414035534    42.816153454   0   0  -1 # c
     10      5   1  0.000000     9.836755987    12.567439705    40.500368386   0   0  -1 # c

注意因为采用的联合原子模型，需要将data文件中的mass由12（C原子）改成14（CH2）。data文件中势函数部分之间删除，之后再in文件中按照文献进行设置。

二、in文件部分

1.初始化和读取data文件。

#-------聚乙烯拉伸in文件-------#

#-------初始化设定-------#

units real #单位制
boundary p p p #周期性边界
atom_style full #分子+电荷

#--------模型读取-------#

read_data PE.data #读取MS建立的模型

2.势函数设置

#-------势函数设置  文献：doi:10.1016/j.polymer.2010.10.009-------#

bond_style harmonic #键势势类型
bond_coeff 1 350 1.53 #CH2-CH2键 Kb=350 r0=1.53

angle_style harmonic #键角势类型
angle_coeff 1 60 109.5 #CH2-CH2键角 Kθ=60 键角109.5度

dihedral_style multi/harmonic #二面角势类型
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

pair_style lj/cut 10 #对相互作用势
pair_coeff 1 1 0.122 4.01 #CH2 CH2 ε=0.112 σ=4.0

此处的参数设置完全按照文献中所给：

 3.neighbor参数

#-------设置邻居参数-------#

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

4.驰豫过程

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

#-------平衡态分子模拟-------#

#-------nvt条件下弛豫-------#

dump 1 all custom 100 1.xyz id type x y z #输出nvt弛豫下的原子序号，种类和坐标到nvt.xyz。

fix 1 all nvt temp 500 500 50 #500k下弛豫

thermo 100 #100步输出一次

reset_timestep 0 #起始步数归零

timestep 1

run 10000 #弛豫10000步

unfix 1 #取消 fix 1
undump 1 #取消 dunmp 1

#-------npt条件下弛豫-------#

fix 2 all npt temp 500 500 50 iso 0 0 1000 #npt弛豫。起始终止温度500，波动温度50。起始终止压力0，波动压力1000。

dump 1 all custom 1000 2.xyz id type x y z #输出npt弛豫下所有原子的xyz坐标。

thermo 100 #100步输出一次

reset_timestep 0 #起始步数归零

run 50000 #跑5000步

unfix 2 #取消 fix 1
undump 1 #取消 dump 1

#-------npt条件下弛豫-------#

fix 3 all npt temp 500 200 50 iso 0 0 1000 #npt弛豫。起始终止温度500-200，波动温度50。起始终止压力0，波动压力1000。

dump 1 all custom 1000 3.xyz id type x y z #输出npt弛豫下所有原子的xyz坐标。

thermo 100 #100步输出一次

reset_timestep 0 #起始步数归零

run 50000 #跑5000步

unfix 3 #取消 fix 1
undump 1 #取消 dump 1

#-------npt条件下弛豫-------#

fix 4 all npt temp 200 200 50 iso 0 0 1000 #npt弛豫。起始终止温度200，波动温度50。起始终止压力0，波动压力1000。

dump 1 all custom 1000 4.xyz id type x y z #输出npt弛豫下所有原子的xyz坐标。

thermo 100 #100步输出一次

reset_timestep 0 #起始步数归零

run 50000 #跑5000步

unfix 4 #取消 fix 1
undump 1 #取消 dump 1

驰豫过程的ovito轨迹：

 5.拉伸过程

#-------恒温恒压拉伸-------#

variable temp equal "lx" #设定参数temp等于体系长度
variable L0 equal ${temp} #设定L0等于temp固定值
variable strain equal "(lx - v_L0)/v_L0" #计算应变大小

variable p1 equal "v_strain" #应变
variable p2 equal "-pxx/10000*1.01325" #x方向应力（单位转换）
variable p3 equal "-pyy/10000*1.01325" #y方向应力（单位转换）
variable p4 equal "-pzz/10000*1.01325" #z方向应力（单位转换）
variable step0 equal step

reset_timestep 0  

dump 1 all custom 1000 dump.pe id type x y z #输出拉伸状态下原子坐标

fix 1 all npt temp 200 200 50 y 0 0 1000 z 0 0 1000 drag 2 #恒温恒压，摩擦系数2。
fix 2 all deform 1 x erate 1e-5 units box remap x #x方向拉伸，速率1e-5，box每步变形，原子重新映射到变形后的box。
fix 3 all ave/time 10 5 100 v_p1 v_p2 v_p3 v_p4 file stress.txt #平均化采样

thermo_style custom step temp v_strain pxx pyy pzz lx ly lz epair ebond eangle edihed #自定义热力学输出
thermo 100 #每100步输出一次

reset_timestep 0 #起始步数归零

run 200000 #跑20000步

 拉伸过程的ovito轨迹：

 三、数据分析

利用输出文件stress.txt中的v_p1和v_p2降噪绘制应力应变曲线：

 各种相互作用能随应变的变化：

 自由体积变化：

同理可以得到文中不同链长、不同温度下的拉伸情况。

键长键角分布可以使用VMD统计。取向、缠结参数的计算见原文的参考文献。