3.1. 纤维素-聚丙烯复合材料#

教程来自芬兰 Aalto University 的wiki网站, Antti J. Karttunen 教授团队, 相关论文 《Molecular Dynamics Simulations on the Elastic Properties of Polypropylene Bionanocomposite Reinforced with Cellulose Nanofibrils》(纤维素纳米纤维增强聚丙烯纳米复合材料弹性性能的分子动力学模拟)，发表于 nanomaterials 期刊，DOI: 10.3390/nano12193379。

该论文通过平衡和非平衡分子动力学模拟，研究了由纤维素纳米纤维 (Cellulose Nanofibrils, CNF)、聚丙烯 (PP) 基体和马来酸酐 (Maleic Anhydride, MAH) 偶联剂组成的聚丙烯生物纳米复合体系的弹性性能。

使用工具: #GROMACS
本教程用到的文件打包在此: 百度网盘下载
教程链接: Cellulose nanofibril MD tutorial - Aalto University Wiki(可能需要梯子)

3.1.1. 力场参数#

本教程使用优化后的 Amber14 力场, 其中描述纤维二糖的参数基于 GLYCAM06 力场进行量化计算优化:

具体地, 使用 Acpype.py 工具将参数从 AMBER 转换为 Gromacs 兼容的格式，并利用两级限制静电势 (Restricted Electrostatic Potential, RESP) 方法重新推导了部分电荷，以再现与 AMBER 力场一致的电荷参数化。

将作者打包的力场文件解压后，会有一个力场目录 amber14sb_OL15.ff 和两个文件 residuetypes.dat、specbond.dat, 将它们放在你的工作目录下.

如果你希望 Gromacs 库能够从计算机的任何工作目录都能调用该力场，把解压的所有文件都放到一个合适的地方 <path-to-force-field-directory>，然后将路径添加到系统的 bashrc 文件中:

echo export GMXLIB=”<path-to-force-field-directory>” >> ~/.bashrc
# or
export GMXLIB=”<path-to-force-field-directory>”

查看 amber14sb_OL15.ff 目录中, 额外添加的参数文件有:

New-glc.rtp: 葡萄糖原子电荷和键描述的残基拓扑
New-glc.r2b: 葡萄糖残基名称数据库
New-glc.arn: 葡萄糖原子重命名数据库
pp.rtp: 聚丙烯原子电荷和键描述的残基拓扑
pp.r2b: 聚丙烯残基名称数据库
pp.arn: 聚丙烯原子重命名数据库

3.1.2. 水中的纤维素-平衡MD#

3.1.2.1. 纤维素模型构建#

使用 Cellulose-Builder 工具生成纤维素结构（cellulose fibril）, Cellulose-Builder 的安装和使用介绍见 必备技能 部分, 命令:

cd cellulose-builder_July_2013
./cellulose-builder fibril 20

如果脚本正确运行，你将看到一个名为 crystal 的新目录，其中将有一个名为 crystal.pdb 的 PDB 格式结构文件. 这是一个由36条葡萄糖链组成的, 每条链含10个纤维二糖的纤维素原纤维(cellulose elementary fibril).

https://2301032.xyz/cellulose-pp-1.jpg — 图. 3.1.1 纤维素纳米纤维的结构#

检查一下 Cellulose-builder 生成的 crystal.pdb 文件，将里面的残基命名 BGLC 改为 IBG（有时候也会命名成BGLCX）, 使其匹配力场中记录的残基名.

sed -i "s/BGLC/IBG /" crystal.pdb

3.1.2.2. 构建拓扑文件#

这里为了模拟具有无限长度的纤维素纳米纤维，需要每条链上的末端葡萄糖单体跨越盒子的周期性边界键合在一起, 以减少了链的波动，避免了松散末端链段带来的问题。

(1). 调用pdb2gmx生成一个初步的 topology file

gmx pdb2gmx -f crystal.pdb -o pdb2gmx.gro -p topol-pdb2gmx.top -water tip3p

输入 1 以选择力场 AMBER14sb-OL15 + beta-glucose GLYCAM06h
忽略警告 long bond warning，因为pdb2gmx不支持周期性共价键的生成, 将在后面使用 x2top 重新生成这些越过盒子的键参数。

(2). 设置盒子大小

gmx editconf -f pdb2gmx.gro -o mfc-20mono.gro -box 7 8 10.38

注意这里的盒子尺寸 10.38nm 是为了匹配纤维素z轴方向上的长度, 使得接下来能够通过xtop命令生成越过盒子的键.

(3). 如前所述, pdb2gmx 不支持周期性共价键的生成, 我们要使用x2top生成一个具有周期性键连参数的拓扑。因此我们只需要 topol-pdb2gmx.top 的[atom]部分, 因为后面的键参数是不正确的(因此才会有 long bond warning), 后面的参数由x2top命令生成。

gmx x2top -f mfc-20mono.gro -o topol-x2top.top -ff select -pbc -name IBG -noparam -alldih

融合 topol-pdb2gmx.top（从头到[bonds]之间）和 topol-x2top.top（从bonds到system之间）, 生成一个新的拓扑文件 newtopol.top
修改 [ angles ] 或者[ dihedrals ] 的函数类型，和 topol-pdb2gmx.top 中的函数类型相同.

至此我们就搭建好了最基础的模拟体系，获得了用于GROMACS软件模拟所需的gro文件和top文件, 如果你还希望加点其他东西（比如蛋白或者其他分子的），像常规的蛋白质-配体体系搭建方法一样，整合 .gro 文件，在 .top 中添加 .itp 文件的引入即可

3.1.2.3. 溶剂化盒子#

(1). 添加水分子

gmx solvate -cp mfc-20mono.gro -cs spc216.gro -p newtopol.top -o box-sol.gro

3.1.2.4. 能量最小化#

在开始动力学之前，使用能量最小化来松弛系统并消除任何空间冲突。

gmx grompp -f em.mdp -p newtopol.top -c box-sol.gro -o em
gmx mdrun -deffnm em

可视化观察em.gro文件, 判断纤维素整体是否偏移出了盒子, 使用一些命令将结构作为一个整体放回盒子中心。

echo “Cellulose System” | gmx trjconv -f em.gro -s em.tpr -pbc mol -center -o em.gro

注意:

模拟周期性/无限（periodic/infinite）纤维素链参数所需的一个关键参数是在mdp文件中添加 “periodic_molecules = yes”。你可以在提供的 em.mdp 文件中看到这个参数.
结束后可以使用 gmx energy 来检查一下体系的能量变化, 选择“Potential”选项编号，然后选择“0”终止分析输入。使用 Xmgrace 或任何其他选择的程序来可视化能量图。

3.1.2.5. 平衡#

然后将对溶剂化的纤维进行两步平衡: 首先用NVT系综进行500ps温度平衡，然后进行500ps的恒温恒压NPT平衡。

(1). NVT平衡: 将系统的温度设定到 300K

gmx grompp -f nvt.mdp -p newtopol.top -c em.gro -r em.gro -o nvt 
gmx mdrun -deffnm nvt -dlb yes	
# -dlb也可以设置为 auto

在成功运行结束时，可以分析体系的温度演变: gmx_mpi energy -f nvt.edr -o nvt.xvg, 选择“Temperature”选项编号，然后输入“0”以结束。使用绘图工具（Xmgrace）可视化温度波动。
同样可以使用 trjconv 工具使结构完整：echo “Cellulose System” | gmx trjconv -f nvt.gro -s nvt.tpr -pbc mol -center -o nvt.gro

(2). NPT平衡: 使用 V-rescale 恒温器、弱耦合 Berendsen 压浴, 并对重原子进行位置约束.

注意这里采取了半各向同性压力耦合（pcoupltype = semiisotropic），以分离原纤维沿着链轴（Z轴）和其他两个轴感受压力的方式。压力耦合在X和Y方向上是各向同性的。可以在提供的 npt.mdp 文件中看到这个参数.

gmx grompp -f npt.mdp -p newtopol.top -c nvt.gro -r nvt.gro -t nvt.cpt -o npt 
gmx mdrun -deffnm npt -dlb yes 

结束后可以使用以下命令分析压力平衡：gmx_mpi energy -f npt.edr -o npt.xvg, 选择“ pressure“选项编号，然后输入“0“结束输入, xmgrace可视化输出文件 npt.xvg
同样可以使用 trjconv 工具使结构完整：echo “Cellulose System” | gmx trjconv -f npt.gro -s npt.tpr -pbc mol -center -o npt.gro

3.1.2.6. 执行MD模拟#

接下来使用与 NPT 相同的设置进行 MD，但不施加任何位置限制，同时使用 “Parrionello Rahman” 进行更紧密的压力耦合。可以在提供的 md.mdp 文件中看到这个参数.

gmx grompp -f md.mdp -p newtopol.top -c npt.gro -o topol 
gmx mdrun -s topol -dlb yes 

至此我们完成了水中纤维素纳米纤维的平衡分子动力学模拟, 使系统达到了热力学平衡状态，使密度、温度等参数稳定，为后续通过非平衡分子动力学模拟来计算力学性能提供稳定的初始结构.

3.1.2.7. 分析#

(1). 使用 trjconv 工具再次使文件 _confout.gro 中_的输出结构完整:

echo “Cellulose System” | gmx trjconv -f confout.gro -s topol.tpr -pbc mol -center -o confout-10ns.gro`

(2). 评估纤维素原纤维的 RMSD 波动

echo “Cellulose Cellulose” | gmx rms -f traj.trr -s topol.tpr -o rmsd

均方根波动可能会有初始上升，但它应该在 2 ns 内稳定下来，只有很小的波动.

(3). 检查势能的演化

gmx energy -f ener.edr -o energy.xvg

选择“Potential”、“Temperature” 和 “Density” 选项的编号，然后选择“0“结束提示。打开绘图文件，并评估是否需要根据能量景观和密度波动续算。

3.1.3. 聚丙烯基体-退火模拟#

3.1.3.1. 聚丙烯模型构建#

(1). 使用 CHARMM-GUI 界面构建具有 10 个单体的聚丙烯(PP)链结构, 提供的文件名为pp-10mono_raw.pdb

(2). 修该该聚丙烯的残基和原子名称以匹配力场描述

10 个单体的聚丙烯组成: 起始残基（SPP）+ 8 个内部残基（IPP）+ 末端残基（EPP）

(3). 为 PP 链生成拓扑和 Gromacs 结构文件

echo " 1" | gmx pdb2gmx -water tip3p -f pp-10mono.pdb -o single-chain-10mono.gro -p topol-pp.top -i single-chain-10mono-posre.itp

3.1.3.2. 构建盒子#

像在上一节中在纤维素原纤维盒中加水一样, 接下来把聚丙烯填充到一个空的盒子里.

(1). 首先将 “top” 文件转换为“itp“格式的拓扑文件, 并创建一个新的拓扑文件“newtopol.top“, 将 “single-chain-10mono.itp” 包含其中.

sed -n '/moleculetype/,/endif/p' topol-pp.top | sed 's/Composite_chain_X/Composite        /' >& single-chain-10mono.itp 

sed '/moleculetype/,/endif/d' topol-pp.top | sed "/Include water topology.*/i ; Include PolyPropylene chain topology\n#include \"single-chain-10mono.itp\"\n " >& newtopol.top

(2). 根据最终盒子体积和聚丙烯基体的密度来计算需要添加多少个10单体的 PP单链.

Atomic mass of 10-monomer chain: 422.796 amu or 702.0692136431 x 10⁻²⁴ g   
Unit cell volume: 512 nm³ or 5.12 * 10-19 cm³  
Target chain count X (10-monomer chain): (Density*Volume)/(Chain mass)   
For our system, X: (0.85 g/cm³ * 5.12* 10-19 (cm³))/ (702.0692136431 * 10⁻²⁴ g) = ~620 chains

(3). 使用以下 Gromacs 命令创建具有约 620 个聚合物链的大盒子 , 后面会使用模拟退火方法将盒子压缩到目标体积/密度.

gmx insert-molecules -ci single-chain-10mono.gro -box 12 12 12 -nmol 620 -scale 0.48 -o box-large.gro >& inmol-620chains.log

注意:

默认 vdw 半径值“0.57”表示水的密度为 1 g/cm 3，但在这里我们将使用“0.48”半径值来获得约 0.85 g/cm 3 的密度。
在实践中，大盒子的大小是通过试错法确定的，直到我们找到适合目标聚合物链数的盒子大小

https://2301032.xyz/cellulose-pp-2.jpg — 图. 3.1.2 向盒子中填充纤聚丙烯链#

(4). 编辑拓扑文件“newtopol.top“的[ molecules ]部分，添加 620 条 PP 链。

[ molecules ]   
; Compound        #mols   
Polypropylene_chain_X  620

3.1.3.3. 能量最小化#

接下来，运行具有 1000 kJ/mol 阈值的能量最小化。我们可以使用相同的 em.mdp 文件，但需要注释掉 periodic_molecules 行，因为这里有有限的链长度。

gmx grompp -f em.mdp -p newtopol.top -c box-large.gro -o em
gmx mdrun -deffnm em

最小化运行的输出结构会有一些跨过盒子边界而“断开”的链。使用 trjconv 工具重建它：

echo “Polypropylene System” | gmx trjconv -f em.gro -s em.tpr -pbc mol -center -o em.gro

可以使用 “gmx energy” 工具分析最小化运行的势能曲线。

3.1.3.4. 平衡#

使用参数文件 nvt.mdp 以与之前相同的方式开始 1 ns 温度平衡运行，但对参数进行了少量修改：

;periodic_molecues = yes   
; Temperature coupling   
tcoupl                   = V-rescale    ; V-rescale Thermostat   
tc-grps                  = System   
tau-t                    = 0.1          ; coupling time constant [ps]   
ref-t                    = 300          ; reference temperature [K]

gmx grompp -f nvt.mdp -p newtopol.top -c em.gro -r em.gro -o nvt
gmx mdrun -deffnm nvt -dlb yes
echo “Polypropylene System” | gmx trjconv -f nvt.gro -s nvt.tpr -pbc mol -center -o nvt.gro

结束后可以使用能量工具检查温度演变
在输出文件 nvt.gro 可视化, 将看到 PP 链开始弯曲并形成分子间相互作用。

3.1.3.5. 执行退火模拟#

(1). 现在温度达到参考值，接下来平衡压力并执行退火模拟过程, 在等温等压系综 (NPT)(300 K→450 K→300 K) 中对系统进行退火模拟, 该过程允许模拟聚合物链的自然构象和相互纠缠的状态。这里使用各向同性耦合来允许三个盒子向量的均匀缩放，否则，我们可能最终得到一个非立方体（长方体）盒子。使用提供的mdp文件npt-anneal-1.mdp

gmx grompp -f npt-anneal-1.mdp -p newtopol.top -c nvt.gro -r nvt.gro -t nvt.cpt -o npt-anneal-1
gmx mdrun -deffnm npt-anneal-1 -dlb yes -dds 0.6
echo “Polypropylene System” | gmx trjconv -f npt-anneal-1.gro -s npt-anneal-1.tpr -pbc mol -center -o npt-anneal-1.gro

注意 mdrun 运行输入的“dds”值的附加指定。该值指示在运行期间单元尺寸可以增加/减少的分数的倒数，同时保持系统单元的区域分解。如果没有这个参数，退火运行很可能会由于盒子大小的压缩而崩溃。
使用 PyMol 打开结构文件 npt-anneal-1.gro 和 nvt.gro。在退火 NPT 运行后，基质单元显著收缩，并且链完全缠结，如从实际 PP 基质所预期的。

gmx energy -f npt-anneal-1.edr -o npt-anneal-1.xvg

在输入提示符上选择以下选项：Potential, Temperature, Pressure, Box-X, Box-Y, Box-Z, Volume, Density, 0。打开输出的 xvg 图文件并检查每个变量随时间的演变。

(2). 接下来，用新的盒子开始另一个 5 ns 时间的退火模拟。通过将系统从 300 K -> 600 K 加热约 5-10 次来平衡聚合物，每次 300K→600K→300K 的循环在 500 ps 的时间内完成。这些设置可以在提供的npt-anneal-2.mdp中看到.

gmx grompp -f npt-anneal-2.mdp -p newtopol.top -c npt-anneal-1.gro -t npt-anneal-1.cpt -r npt-anneal-1.gro -o npt-anneal-2 
gmx mdrun -deffnm npt-anneal-2 -dlb yes -dds 0.6
echo “Polypropylene System” | gmx trjconv -f npt-anneal-2.gro -s npt-anneal-2.tpr -pbc mol -center -o npt-anneal-2.gro
gmx energy -f npt-anneal-2.edr -o npt-anneal-2.xvg

在输入提示符上选择以下选项：Potential, Temperature, Pressure, Box-X, Box-Y, Box-Z, Volume, Density, 0。打开输出 xvg 图文件并检查每个变量的时间演变。

https://2301032.xyz/cellulose-pp-3.png — 图. 3.1.3 退火模拟过程中体系温度, 势能, 体积, 密度的波动#

3.1.3.6. MD#

接下来，我们运行额外的 10 ns NPT 平衡（没有退火和位置限制）以稳定新电池。通过最终的平衡运行，目的是稳定细胞大小和原子相互作用。使用提供的参数文件 npt-3.mdp

gmx grompp -f npt-3.mdp -p newtopol.top -c npt-anneal-2.gro -t npt-anneal-2.cpt -r npt-anneal-2.gro -o npt-3 
gmx mdrun -deffnm npt-3 -dlb yes

注意:

如果模拟崩溃，则以 1-fs 的较小时间步长重新启动运行。
运行成功结束后，使用适当的参数运行 trjconv 和 energy 命令，使输出结构（npt-3.gro）再次完整，并分析体系盒子的矢量、体积和 PP 密度的稳定性。

3.1.4. 复合材料的模拟#

3.1.4.1. 模型和拓扑构建#

我们保持纳米原纤维模型的盒子尺寸不变（7 * 8 * 10.38 nm3），并用 PP 链填充盒子。为了找到填充密度为 0.85 g/cm3 的盒子所需的链的数量，我们需要计算盒子中的空空间。使用“gmx sasa” 工具计算纤维素纳米纤丝所占据的体积.

这部分使用提供的动力学参数文件mdp-biocomp.tgz

# 解压提供的压缩包 mdp-biocomp.tgz
tar -xzf mdp-biocomp.tgz

gmx grompp -f em.mdp -p newtopol-mfc20.top -c mfc-20mono.gro -o sasa-mfc20
# newtopol-mfc20.top 即上述纤维素的 newtopol.top
echo "Cellulose" | gmx_mpi sasa -f mfc-20mono.gro -s sasa-mfc20.tpr -o area-cellulose -tv volume-cellulose

使用文本编辑器打开输出文件“volume-cellulose.xvg” ，文件的最后一行包含纤维素纳米原纤维占据的体积（见下图）

我们知道盒体积为 581.28 nm3，因此 PP 基质可用的体积将为 487.943 nm3。现在，我们可以计算要添加到盒子中的 PP 链的数量：

PP 10-monomer chain's mass = 422.796 amu or 702.0692136431 x 10⁻²⁴ g   
Number of PP chains, X = (Density*Volume)/(Chain mass)  
X = ~590 chains

为了获得 0.85 g/cm³的目标密度，将需要 590 个 10 单体长度的 PP 链。为了容纳足够的 PP 链的添加，我们现在通过延长 X 和 Y 轴但保持盒的 Z 轴与之前相同来产生用于纳米原纤系统的新盒，保证无限长度（周期性）的纤维素原纤并允许跨越盒边界的共价连接。

gmx editconf –f mfc-20mono.gro -o mfc-20mono-largebox.gro -c –box 18.0 18.0 10.38

打开结构输出“mfc-20 mono-largebox.gro” ，确保原纤维位于中心，且盒子的 z 轴不变。接下来，我们将 PP 链插入到盒子中：

gmx insert-molecules -f mfc-20mono-largebox.gro -ci single-chain-10mono.gro -o biocomp-mfc20-pp10.gro -nmol 590 -scale 0.48 >& im-nmol590.log

https://2301032.xyz/cellulose-pp-5.jpg — 图. 3.1.4 *将聚丙烯分子插入到纳米纤维素盒子中*#

注意：新盒的尺寸必须足够大，以容纳所需的 PP 链。在 “insert-molecules” 步骤检查文件 im-nmol590.log ，确保添加了足够(590条)的链。如果没有，则增加盒子 X 和 Y 维度上的尺寸, 并重复插入分子的步骤。

创建一个新的拓扑文件，以包含 cellulose 和 PP 拓扑。

sed "/Include water topology.*/i ; Include PolyPropylene chain topology\n#include \"single-chain-10mono.itp\"\n " newtopol-mfc20.top >& newtopol.top
# newtopol-mfc20.top 即上述纤维素的 newtopol.top
echo "Polypropylene_chain_X 590" >> newtopol.top

注意区分三个部分中的newtopol.top

3.1.4.2. 能量最小化#

gmx grompp -f em.mdp -p newtopol.top -c biocomp-mfc20-pp10.gro -o em
gmx mdrun -deffnm em  gmx_mpi mdrun -defnm em
echo "Cellulose System" | gmx trjconv -f em.gro -s em.tpr -pbc mol -center -o em.gro

分析最小化运行期间系统的潜在演变，并评价是否适合继续进行平衡步骤。

3.1.4.3. 平衡#

3.1.4.3.1. NVT#

gmx grompp -f nvt.mdp -p newtopol.top -c em.gro -r em.gro -o nvt
gmx mdrun -deffnm nvt -dlb yes
echo “Cellulose System” | gmx trjconv -f nvt.gro -s nvt.tpr -pbc mol -center -o nvt.gro

3.1.4.3.2. annealing#

接下来，我们实施 2 步模拟退火 MD 模拟（类似于聚丙烯教程），以使盒子压缩并实现聚丙烯基质的 0.85 g/cm 3 密度。

第一次运行是 500 ps 的（300 K → 450 K → 300 K）的单次退火循环。
第二次, PP 基体将经历 10 个退火循环（300 K -> 600 K），其中 300 K → 600 K → 300 K 的每次重复在 500 ps 内完成。

# 第一次annealing
gmx grompp -f npt-anneal-1.mdp -p newtopol.top -c nvt.gro -r nvt.gro -o npt-anneal-1
gmx mdrun -deffnm npt-anneal-1
echo “Cellulose System” | gmx trjconv -f npt-anneal-1.gro -s npt-anneal-1.tpr -pbc mol -center -o npt-anneal-1.gro
gmx energy -f npt-anneal-1.mdp -o npt-anneal-1.xvg

使用 PyMol 打开结构文件 npt-anneal-1.gro 和 nvt.gro。请注意，盒子沿着 X 和 Y 轴发生了变化，PP 链也相互缠绕。
在输入提示符上选择以下选项： Potential, Temperature, Pressure, Box-X, Box-Y, Box-Z, Volume, 0. 打开输出的 xvg 图文件并检查每个变量的时间演变。

# 第二次annealing
gmx grompp -f npt-anneal-2.mdp -p newtopol.top -r npt-anneal-1.gro -t npt-anneal-1.cpt -r npt-anneal-1.gro -o npt-anneal-2
gmx mdrun -deffnm npt-anneal-2 -dlb yes -dds 0.6
echo “Cellulose System” | gmx trjconv -f npt-anneal-2.gro -s npt-anneal-2.tpr -pbc mol -center -o npt-anneal-2.gro
gmx energy -f npt-anneal-2.edr -o npt-anneal-2.xvg

这里，再次对纤维素原纤维使用单独的温度耦合，以将其保持在 300 K 的参考温度，以避免任何结构变形
在输入提示符上选择以下选项： Potential, Temperature, Pressure, Box-X, Box-Y, Box-Z, Volume, 0. 打开输出的 xvg 图文件并检查每个变量的时间演变。

https://2301032.xyz/cellulose-pp-6.jpg — 图. 3.1.5 复合材料退火模拟过程中体系温度, 势能, 体积, 密度的波动#

https://2301032.xyz/cellulose-pp-7.jpg — 图. 3.1.6 复合材料盒子退火前后的尺寸变化和模型变化#

3.1.4.4. MD#

为了验证 PP 填充是稳定的，并且盒子矢量在 NPT 条件下不会进一步偏离，在 NPT 系综条件温度 300 K 和压力 1.0 bar 下进行另一个 10 ns 平衡运行。该步骤还允许 PP 基体在参考温度下在较长时间内平衡并稳定目标密度。此运行的参数文件应已位于当前目录中，名称为 “npt-3.mdp”。

gmx grompp -f npt-3.mdp -p newtopol.top -c npt-anneal-2.gro -o npt-3
gmx mdrun -deffnm npt-3 -dlb yes

在运行结束时，使用 “gmx 能量” 分析盒子矢量，以确保稳定的体积和 PP 密度。另外，请记住使用 trjconv 工具再次使结构完整。

3.1.5. 弹性性质研究-非平衡MD#

在对所有上述各种模型在 10 ns 平衡MD运行后, 进行了单轴压缩模拟.

#待更新