从零开始的磷脂双分子层模拟

事情起源于一篇讨论跨膜过程的文章：10.1021/acscentsci.0c00893。

研究磷脂双分子层的行为往小处说可以研究化合物跨膜特性，往大了说还可以研究膜蛋白识别过程嘛。 是个从一个小小硕士的视角看起来，是十分有前途的路子。下面简单记录一下首次磷脂双分子层模拟。

生成磷脂双分子层

磷脂双分子层的生成一般通过CHARMM-GUI或者VMD来生成，前者一般要学校邮箱才能注册， 而后者正常注册个账户就能下载。因此为了不失一般性选择使用VMD来生成磷脂双分子层。

VMD可以调用Extension->Modeling->Membrane Builder来生成：

该插件只能生成POPC（磷脂酰胆碱，卵磷脂）和POPE（磷脂酰乙醇胺，脑磷脂）两种类型的磷脂，虽说有点少， 而且也不能生成混合磷脂层，但又不是不能用。使用该插件，特别是在Windows上使用， 务必需要在工作目录中使用start命令启动VMD，因为该插件会在VMD运行的目录下创建生成好的psf和pdb文件。 如果直接启动，至少在1.9.3版的VMD上，会出现VMD生成的处在VMD安装目录的文件只有它自己能看见的BUG。

比如这里选择生成一个100x100的膜，于是就得到了Z轴两端薄薄一个水层的磷脂双分子模型：

哪怕作为一个不是那么专业的人，也会觉得这个水层也忒薄了，而且还没有加离子。因此还需要加个水盒， 比如说Z轴两端各加10埃水：

solvate membrane.psf membrane.pdb +z 10 -z 10

然后再用Extension->Modeling->Add ions随机加点NaCl，比如这边加到0.15M。

首先，可以看到一个明显的界面，当然这其实是个小问题；其次可以看到磷脂疏水部分也出现了水分子， 这个问题其实就相对比较讨人厌了，不过实际做下来，只要能升温到310K，稍微平衡几个ps就能正确地将水分子排出。 使用measure命令量一下周期性边界盒子的大小:

measure minmax [atomselect top all]
measure center [atomselect top all]

接下来，在这里采取了三次升温平衡，最终取得了还算令人满意的结构，参考配置文件如下：

# input
coordinates             ionized-minimization.coor
structure               ionized.psf
parameters              /home/hzy/lipid/toppar_c36_jul21/par_all22_prot.prm
parameters              /home/hzy/lipid/toppar_c36_jul21/par_all36_carb.prm
parameters              /home/hzy/lipid/toppar_c36_jul21/par_all35_ethers.prm
parameters              /home/hzy/lipid/toppar_c36_jul21/par_all36_cgenff.prm
parameters              /home/hzy/lipid/toppar_c36_jul21/par_all36_lipid.prm
parameters              /home/hzy/lipid/toppar_c36_jul21/par_all36_lipid_ljpme.prm
parameters              /home/hzy/lipid/toppar_c36_jul21/par_all36_na.prm
parameters              /home/hzy/lipid/toppar_c36_jul21/par_all36m_prot.prm
parameters              /home/hzy/lipid/toppar_c36_jul21/toppar_water_ions_namd.str
paratypecharmm          on
mergeCrossterms         yes

# output
set output              "/home/hzy/lipid/ionized-eq1"
outputname              $output
dcdfile                 ${output}.dcd
xstFile                 ${output}.xst
dcdfreq                 250
xstFreq                 250

binaryoutput            no
binaryrestart           no
outputEnergies          500
outputPressure          500
restartfreq             500

fixedAtoms              off

# Periodic Boundary Condition
cellBasisVector1        103.0    0.0    0.0
cellBasisVector2          0.0  102.0    0.0
cellBasisVector3          0.0    0.0   96.0
cellOrigin                0.86   1.59  10.62
wrapAll                 on

# Basic dynamics
exclude                 scaled1-4
1-4scaling              1
COMmotion               no
dielectric              1.0

# Simulation space partitioning
switching               on
switchdist              9
cutoff                  10
pairlistdist            12

# Multiple timestepping
firsttimestep           0
timestep                1
stepspercycle           10
nonbondedFreq           2
fullElectFrequency      2

# Temperature control

set temperature         310
temperature             $temperature;  # initial temperature

# Constant Temperature Control
langevin                on    ;# do langevin dynamics
langevinDamping         1     ;# damping coefficient (gamma) of 1/ps
langevinTemp            $temperature
langevinHydrogen        off    ;# don't couple langevin bath to hydrogens

# PME (for full-system periodic electrostatics)
PME                     yes
PMEGridSpacing          1.0

# Constant Pressure Control (variable volume)
useGroupPressure        yes ;# needed for rigidBonds
useFlexibleCell         yes  ;# seems needed for lipid
useConstantArea         yes

langevinPiston          on
langevinPistonTarget    1.01325 ;#  in bar -> 1 atm
langevinPistonPeriod    100.0
langevinPistonDecay     50.0
langevinPistonTemp      $temperature

# Scripting
minimize                20000
for { set TEMP 0 } { $TEMP < 100 } { incr TEMP 20 } {
    reinitvels          $TEMP
    langevinTemp        $TEMP
    run                 5000
}

minimize                20000

for { set TEMP 0 } { $TEMP < $temperature } { incr TEMP 10 } {
    reinitvels          $TEMP
    langevinTemp        $TEMP
    run                 5000
}

reinitvels              $temperature
langevinTemp            $temperature
langevinPistonTemp      $temperature

run                     100000;

minimize                20000

for { set TEMP 200 } { $TEMP < $temperature } { incr TEMP 5 } {
    reinitvels          $TEMP
    langevinTemp        $TEMP
    run                 5000
}

reinitvels              $temperature
langevinTemp            $temperature
langevinPistonTemp      $temperature

run                     500000;

水分子对磷脂双分子层的两面包夹芝士俨然已成：

加入研究对象

所研究的对象分子可以通过CGenFF得到一个 与CHARMM力场相容的小分子力场文件，不妨称之为target.str文件。 然后打开一个空白VMD，读入双分子层和对象分子，移动对象分子到适当位置， 导出此时对象分子的坐标为pdb文件，不妨称之为target.pdb文件。

接着重新打开一个空白VMD，接下来就在终端里操作了，先用topology引入拓扑定义。 接着用readpsf把双分子层结构定义读入，同时用coordpdb读入相应坐标。 最后创建一个不重名的segment，这里取名叫XCUM，然后把坐标也读入：

topology  "C:/Users/hzy/Desktop/membrane demo/toppar_c36_jul21/top_all22_prot.rtf"
topology  "C:/Users/hzy/Desktop/membrane demo/toppar_c36_jul21/top_all36_carb.rtf"
topology  "C:/Users/hzy/Desktop/membrane demo/toppar_c36_jul21/top_all35_ethers.rtf"
topology  "C:/Users/hzy/Desktop/membrane demo/toppar_c36_jul21/top_all36_cgenff.rtf"
topology  "C:/Users/hzy/Desktop/membrane demo/toppar_c36_jul21/top_all36_lipid.rtf"
topology  "C:/Users/hzy/Desktop/membrane demo/toppar_c36_jul21/top_all36_lipid_ljpme.rtf"
topology  "C:/Users/hzy/Desktop/membrane demo/toppar_c36_jul21/top_all36_na.rtf"
topology  "C:/Users/hzy/Desktop/membrane demo/toppar_c36_jul21/toppar_water_ions_namd.str"
topology  "C:/Users/hzy/Desktop/membrane demo/target.str"

readpsf "Z:/lipid/ionized.psf"
coordpdb "Z:/lipid/ionized-eq.pdb"

segment XCUM { pdb "C:/Users/hzy/Desktop/membrane demo/target.pdb" }
coordpdb "C:/Users/hzy/Desktop/membrane demo/target.pdb" XCUM

writepsf "C:/Users/hzy/Desktop/membrane demo/test2.psf"
writepdb "C:/Users/hzy/Desktop/membrane demo/test2.pdb"

至此就得到了复合模型，就差最后加水盒了：