实验四：介观动力学模拟_27396

《计算材料学》实验讲义

实验八：介观动力学模拟

一、前言

1、介观模拟简介

长期以来，化学家致力于从分子水平研究物质及其变化，而化学工程工作者主要研究物质在宏观体系的行为，介观层次的化学正是联系微观及宏观的桥梁，是从分子到材料的必由之路，同生命过程也有密切的关联。

由于介观模拟能够模拟的空间尺度（纳米到微米）、时间尺度（纳秒到微秒）更大，应用介观模拟方法可以模拟更加复杂的体系，例如：高分子熔体，高分子稀溶液自组装，表面活性剂溶液自组装，磷脂膜等胶体化学，高分子，生物大分子相关的内容。

目前介观模拟的方法很多，例如耗散颗粒动力学模拟方法(dissipative particle dynamics，DPD)，它是根据Hoogerbrugge和Koelman提出的一种针对柔性(soft)球模型流体动力学的模拟，并通过引入粒子间的谐振动势，来模拟聚合物的性质；元胞动力学方法（CDS），基于重整化群理论，对时间相关的Ginzburg-Landau方程直接用数值计算的方法在离散空间上进行描述。其中单个元胞的演化通常用双曲正切函数表示；动态密度泛函方法（DDFT或MesoDyn），应用于高分子体系，建立在粗粒化高斯链模型的基础上，实际上是一个动态的自洽场方法，使用了朗之万方程（Langevin’s equation）来描述体系演化的动力学。 （1）MS-Mesocite简介

MS Mesocite是一个基于粗粒度模拟方法的、可以对广泛体系进行模拟研究的分子力学工具集，模拟的对象大小尺寸在纳米到微米尺度范围，相应地，模拟变化的时间范围落在纳秒至微秒区间。MS Mesocite的模拟对象遍及多种工业领域，比如复合材料、涂料、化妆品以及药物控缓释等，它可以提供流体在平衡态下、在有剪切力存在下以及其它受条件下的结构与动力学性质。MS Mesocite的突出特点是使用完全区别于传统介观模拟技术，转而采用力场(Forcefield)方法—比如MS Martini力场—来描述粗粒度之间的相互作用，从而得到体系的结构、

1

和动力学特性，分析函数主要有角度分布，密度分布，径向分布函数，二面角分布，均方根位移等。同时，您还可以使用力场编辑工具对MS Mesocite的力场进行编辑，以获得满足特殊要求的力场，从而拓展了MS Mesocite的应用范围。

应用Mesocite进行动力学模拟时，最主要的是得到精确的力场。Martini力场，是由Marrink提出的，可以应用于生物分子体系。Martin力场中包括四种主要的力场类型：极性（polar-P）、非极性（apolar-C）、无极性（nonpolar-N）、带电（charged-Q）。每种力场类型又分为若干子类型，极性和非极性根据极性高低下分有五种类型（用下坐标1-5表示），无极性和带电的更具氢键结合能力分为四种类型（d-氢键供体，a氢键受体，da-两个都有，o-都没有），这样使得Martini力场能够更加精确的描述体系性质，应用于更多不同的有机分子体系。

二、实验目的

1、了解介观模拟方法及应用领域 2、了解Martini力场的

3、掌握Mesocite模块的基本操作

三、实验内容

以下以介观动力学模拟脂质双分子层为例，熟悉Mesocite的基本操作。 1、打开MS，选择created new project，键入CG-bilayer作为工程的名称，点击OK。本实例是在软件所有参数在默认的情况下进行的，选择Tools-Settings Organizer，选中CG-bilayer，点击Reset。

2、建脂质分子，建模过程要用到Mesostructure toolbar

，如在

工具栏中没有此建模工具，点击菜单栏中的view-toolbar-mesostructure，调出此建模工具。

（1）点击Bead Types按钮

，打开Bead Types 对话框。

点击Properties…按钮，打开 Bead Type Properties 对话框，点击Defaults…按钮，设置Mass为72，Radius为2.35，关闭Bead Type Defaults和Bead Type Properties对话框。

在Bead Types对话框中，定义一下珠子类型：C、GL、PO和NC，关闭对话框。

2

（2）点击Mesomolecule按钮，打开Build Mesomolecule对话框。

定义粗粒化分子，依次选择4个C、1个GL、1个PO、1个GL和4个C，确定不选Randomize order within repeat unit，点击Build按钮。

在Mesomolecule.xsd文件中左击PO珠子，删除Build Mesomolecule对话框中所有的珠子。

选中Add to branch points，点击more…按钮，打开Mesomolecule Branches对话框。设置Number of branches to attach为1，关闭对话框。

在Build Mesomolecule对话框中选择1个NC。点击Build按钮。

（在显示面板中右击，选择Label，打开label对话框，在properties一栏中选择BeadTypeName，点击Apply，可以检测建立的粗粒化分子是不是正确，可以对比下图。

（3）关闭Build Mesomolecule对话框。在Project Explorer，把Mesomolecule.xsd文件名改为DPPC.xsd。我们得到以下粗粒化分子结构：

3、更改Martini力场，分配力场，优化脂质分子。

（1）选择Modules -Mesocite - Forcefield Manager或点击Mesocite tools

，选择

Forcefield Manage，选择MS Martini，点击>>，打开力场文件。在Project Explorer中，把文件名改为MSMartiniCIS.off。

（2）打开MSMartiniCIS.off文件，点击Interactions。在Show interaction下拉选项中选择Angle Bend。在空白框中，设置Fi 和Fk 到Na 以及 Fj 到Qa。改变 Functional Form 为Cosine Harmonic设置TO为120，KO为10.8。关闭力场文件并保存。

（3）选择Modules | Mesocite | Calculation或点击Mesocite toolsCalculation;

打开Mesocite Calculation对话框，点击Energy，在Forcefield的下拉选项中选择Browse...，在Choose Forcefield对话框中选择MSMartiniCIS.off。

3

选择

（4）打开DPPC.xsd文件。按下ALT键，双击任意C类型珠子，选中所有的C类型珠子。在Mesocite Calculation对话框中，点击More...打开Mesocite Proparation options对话框，选择C1，点击Assign按钮。

重复此步，为GL、PO、NC分配力场，分配类型如下表所示：

BeadTypeName C GL PO NC

MS Martini Forcefield Type C1 Na Qa Q0

Charge 0 0 -1.0 1.0

选择PO珠子，在Properties Explorer中，设置Charge为-1，同样把NC设置为1。

（5）在Mesocite Calculation对话框中，点击Setup，改变Task为Geometry Optimization。点击Run按钮。得到以下结构：

（6）在工具栏中，选择Measure/Change按钮，下拉选项中点击Angel，

依次点击左边的C-GL-PO,同样选择右边的PO-GL-C。此时会显示出两个接近156.50的角度，选在两个角度，在Properties Explorer中，设置Angels为230。按下ALT键，双击角度，按下Delete。得到以下分子结构：

4

（7）参照第二步，定义珠子W，用Build Mesomolecule建模工具，建立一个仅包含W的粗粒化分子。更改文件名为solvent.xsd。 4、建立双分子层结构。

（1）选择Build | Build Mesostructure | Mesostructure Template或点击Mesostructure toolbar

中的Mesostructure Template

，打开Build

Mesostructure Template对话框。

改变X、YExtents为，Z Extent为100。在Filler中，键入solvent。点击Build按钮。

在Build Mesostructure Template对话框中，改变Former type为Slab。改变Depth为44.15，Orientation为Along Z。

选中Enable surface packing;

在Filler中键入lipid。点击Add，关闭对话框。

（2）选择Build | Build Mesostructure | Mesostructure或点击Mesostructure toolbar

中的Mesostructure

5

，打开Build Mesostructure对话

框。

solvent filler 中的Mesoscale Molecule，选择solvent.xsd。 lipid filler选择优化的DPPC.xsd。

点击Packing，设置Length scale (L)为1，Density为0.00836； 不选Randomize conformations。

在Packing中，点击More...按钮，打开Bead Packing Options对话框； 双击打开已经优化过的DPPC.xsd。选择NC，点击Create bead Head set from selection按钮。

按下CTRL + D取消选定，之后按下CTRL键，选择尾部的两个C珠子。在Bead Packing Options对话框中，改变Bead tag为Tail，点击Create bead Tail set from selection。

关闭对话框。标记后的DPPC结构如下：

（3）双击mesostructure template.msd。在Build Mesostructure对话框中，点击Build按钮。得到下图所示结构：

6

（4）在菜单栏中选择File | Export...，打开Export对话框，在保存类型下拉选项中选择Materials Studio 3D Atomistic Files (*.xsd)，点击Options...按钮，打开MSD/MTD Export Options对话框，设置Length scale为1，点击OK。

改变文件名为bilayer.xsd，保存在（I）：选择当前工程的根目录下的CG-bilayer Files/Documents。点击保存（S）。此时在project explorer会出现一个名为bialyer.xsd的文件。

（5）在菜单栏中选择File | Save Project，选择Window | Close All。 5、体系优化及动力学过程。

在Project Explorer中，双击bilayer.xsd，打开文件。

（1）分配力场：如第三步中的第四小步，为每种粗粒子珠子分配力场，分配电荷。分配类型如下表所示：

BeadTypeName

C GL PO NC W

（2）第一次构型优化

打开Mesocite Calculation对话框；

点击Energy按钮，在summation method中的Electrostatic的下拉选项中选择Bead based。

确保Mesocite Calculation/Setup中的Task为Geometry Optimization； 选中Mesocite Calculation/Jop Control中的Run inparallel on[ ]of i processors，把可用的CPU调到最大值（此后在几何优化过程，还是动力学过程，为了充分利用服务器，CPU都调到最大值）。

MS Martini Forcefield Type

C1 Na Qa Q0 P4

Charge 0 0 -1.0 1.0 0

7

点击Run。 （3）第二次构型优化

双击打开优化过的bilayer.xsd

在Mesocite Calculation对话框中选择Setup按钮；

点击More...打开Mesocite Geometry Optimization对话框，选中Optimize cell； 关闭Mesocite Geometry Optimization对话框。 点击Run。 （4）动力学优化

双击打开第二次优化过的文件bilayer.xsd

在Setup中，选择Task为Dynamics，点击More...按钮，打开Mesocite Dynamics对话框。

设置Time step为20fs，Dynamic time 50ps，改变Ensemble为NPT。 选择Thermostat按钮，设置Thermostat为Velocity Scale。 点击Barostat按钮，设置Barostat为Andersen。 在Mesocite Calculation对话框中点击Run。 （5）第二次动力学优化

双击打开bilayer Mesocite Dynamics文件夹下的bilayer.xtd件；

在Mesocite Dynamics对话框中选择Thermostat按钮，设置Thermostat为Nose。

设置Q ratio为1600。

设置Time step为40fs，Dynamic time 200ps，点击Dynamics按钮，设置Frame output every为100steps。

在Mesocite Calculation对话框中，选中Restart； 点击Run。

弹出警告对话框，点击Yes。

（6）选择File | Save Project，选择Window | Close All。 6、结果分析，以角度分布和沿Z轴浓度分布为例。 （1）角度分布：

8

文

① 双击打开bilayer Mesocite Restart文件夹下的bilayer.xtd文件；

双击打开DPPC Mesocite GeomOpt文件夹下的DPPC.xsd文件。

在DPPC.xsd下，用Measure/change工具，选择下图所示两个角度。选择GL-PO-GL键角。

② 在菜单栏中选择Edit | Find Patterns，打开Find Patterns对话框。

定义优化过的DPPC.xsd文件作为Pattern document，并且确定键角GL-PO-GL仍然被选中。

改变Match property为BeadTypeName。 打开轨迹文件bilayer.xtd，点击Find。

点击New Sets...按钮，打开Define New Set对话框，键入GL-PO-GL Angles，点击OK按钮。

在bilayer.xtd文件中取消选定。 同样定义sets为 C-PO-C Angles。

③ 选择Modules | Mesocite | Analysis，或点击mesocite tools

，选择Analysis；

打开Mesocite Analysis对话框，在Analysis选项中选择Angle distribution； 在Sets下选项中选择GL-PO-GL Angles，点击Analyze。同理，分析键角C-PO-C Angles。把数据拷贝到excel中，作图可得：

9

（2）Z方向浓度分布

① 双击打开bilayer Mesocite Restart文件夹下的bilayer.xtd文件；

选择Edit |Edit sets，打开Edit sets对话框，按下ALT键，双击任意W珠子，选中了所有的W珠子。

在Edit sets对话框中，点击New…，打开Define New Set对话框对话框，键入W，点击OK。

同理，定义Sets NC、PO、GL、C。

② 选择Modules | Mesocite | Analysis，或点击mesocite tools

，选择Analysis；

打开Mesocite Analysis对话框，在Analysis选项中选择Concentration profile； Sets选择W，选中Specified direction (hkl)，改为0 0 1； 点击Analyze；

同理分析NC、PO、GL、C。把数据拷贝到excel中，作图可得：

10

本实例为软件帮助中的实例教程，参数设置原因可参考Help帮助文件。

参考文献：S.J. Marrink, H.J. Risselada, S. Yefimov, D.P. Tieleman, A.H. de Vries., \"The MARTINI forcefield: coarse grained model for biomolecular simulations.\J. Phys. Chem. B, 111:7812-7824, 2007.

实验步骤及注意的问题

1、 构建5种珠子C、GL、PO、NC、W，设置Mass为72，Radium为2.35 2、 构建DPPC脂质分子 3、 构建力场

4、 给DPPC分子分配力场，优化分子结构，调整角度，获得DPPC分子的最终构型

5、 构建水分子构型 6、 构建盒子 7、 填充盒子 8、 导出.xsd构型文件 9、 给盒子分配力场 10、

初步优化盒子

11

11、 选中optimized cell 进一步优化盒子（注意能量变化曲线，如太高，需进

一步优化，一般需要优化2-3次） 12、

对优化后的构型进行初步分子动力学模拟（time step 20fs，Dynamic time

50ps） 13、

改变参数设置，再次进行分子动力学模拟（time step 40fs，Dynamic time

200ps） 14、

在DPPC分子中选中GL-PO-GL以及C-PO-C两个角度，选定find pattern，

在.xtd轨迹文件中find所有的角度 15、 16、 17、

对角度分布进行分析

在.xtd轨迹文件edit sets，选定5种原子，对其浓度分布进行分析。 将两个角度分布导入EXCEL，将5中原子浓度分布导入EXCEL，分别作图。

四、作业

1、模拟油水混合溶液的分层构型，油选择癸烷作为油相代表

（1）构建癸烷分子并对其进行优化获得如下结构

（2）选择粗粒化珠子 点击粗粒度转化分子工具

，弹出Coarse grain对话框，在对话框中选中Motion

groups，点击more，选中分子中前面个碳原子作为一个粗粒化珠子，点击Motion groups对话框中的create，产生第一个类型的粗粒化珠子；同理，定义中间的两个碳原子作为第二个类型的粗粒化珠子。关闭Motion groups对话框。

12

（2）分配珠子，获得粗粒化珠子结构

点击Coarse grain对话框中的bead typing中的create，获得Bead Typing.std的珠子类型文件。选择patterns，从下拉框中选中该珠子类型文件；点击Coarse grain对话框下面的built文件。获得粗粒化之后的珠子结构。

（3）对珠子分配力场（MS Martini，该力场为软件自带力场，不需要修改参数），两种珠子均属于C1类型，电荷为零，优化珠子构型，获得稳定结构。 （4）构建水珠子，Mass 72，Radius 2.35 （5）构建32*32*32的盒子，在Filtter中键入一个名字，如solvent，点击built，获得一个空的立方体盒子。

（6）在对话框中点击Filtters，点击Add，键入另外一个名字，如Water。关闭对话框。 （7）点击

，填充盒子，选中优化后的粗

13

粒化硅烷分子和水珠子，比例（Relative amount）为1：1；点击Packing，将Length scale设为1，Density设为0.00836，点击built构建填充后的盒子。 （8）将构建后的后缀为.msd的盒子导出为xsd构建文件（如上面所属）。

（9）分配力场（MS Martini），水分子属于P4力场，优化分子构型

（10）进行分子动力学优化（10fs，500ps，500step）获得构型后，对其显示方式（display style）进行设置，扩展构型，观察油水分离现象。

14

（10）分别显示水珠子和两种类型碳珠子的密度分布。

W

C2

C1

15