Martini 3.0教程：参数化一个新的小分子
0. Introduction
本教程将讨论如何为新的小分子构建Martini 3.0拓扑。该教程基于Martini 2教程，但加入了Martini 3.0中的一些更改。建议阅读文献[1]
原文用的例子是1-乙基萘，并使用Gromacs 2019.x或更高版本。这里我们以2-萘酰胺为例，因为我要做的体系是2-萘酰胺（叉腰）
注意：本教程基于Martini 3的公开测试版（v.3.0.b.3.2），将在Martini 3的最终参数可用时立即进行更新。
1）生成原子参考数据
我们将需要原子参考数据来提取CG模型的成键参数。请注意，在本教程中，我们将需要所有氢原子来提取键长，因此请确保您的原子结构包含所有氢。
原文使用了LigParGen服务器[2]来获得原子级（全原子，AA）结构和力场拓扑的方法，当然也可以随意使用自己喜欢的原子力场。这里我改用我习惯的acpype程序（对应Amber力场）或ATB（对应GROMOS力场）来生成拓扑。
通过acpype程序，我们已经获得了2-萘酰胺的结构文件（gro）和Amber拓扑，分别名为NAP_GMX.gro，NAP_GMX.top和NAP_GMX.itp。习惯上将NAP_GMX.top重命名为topol.top。
进入该文件夹，将刚才生成的itp、top、gro文件复 制过来，并执行：
将NAP_GMX.itp的[ atomtypes ]字段剪切到amber99s b-ildn.ff/ffnonbond.itp的原子类型中，并删除topol.top中的[ default ]字段，添加两行：
运行如下命令：
gmx insert-molecules -ci NAP_GMX.gro -nmol 1 -box 3.8 3.8 3.8 -rot xyz -seed 0 -o NAP_box.gro
gmx solvate -cp NAP_box.gro -p topol.top -o NAP_sol.gro
gmx grompp -f em.mdp -c NAP_sol.gro -p topol.top -o em
gmx mdrun -v -deffnm em
gmx grompp -f eq.mdp -c em.gro -p topol.top -o eq
gmx mdrun -v -deffnm eq
gmx grompp -f run.mdp -c eq.gro -p topol.top -o md
gmx mdrun -v -deffnm md上述命令运行能量最小化的命令，然后执行250 ps的NPT平衡，并执行50 ns的MD运行（检查脚本和各mdp文件以了解更多信息）。在这种情况下，使用的溶剂是水；但是，也可输入其它平衡好的溶剂盒子。应该选择一种代表分子将花费大部分时间的环境的溶剂。
2）原子到CG粒子的映射
映射（即，将分子拆分为由CG粒子描述的构建单元）是粗粒化的核心，并且依赖于经验、化学知识和反复的试验。将分子映射到Martini 3.0模型时，应遵循以下准则：
仅考虑非氢原子来定义映射；
避免在两个小粒子之间切断特定的化学基团（例如酰胺或羧酸根）;
保证分子的对称性；此外，希望尽可能保留AA结构的体积和形状；
4:1、3:1和2:1映射的默认选项是Regular（R），Small（S）和Tiny（T）粒子；它们是线性片段的默认选项，例如，辛烷包含的两个4:1片段；
就计算性能而言，R型粒子是最佳选择，其粒子大小可以合理地代表4:1线性分子。
T型粒子特别适合代表芳环的平面度；
S型粒子通常更好地模拟脂族环的“大块”形状。
应优化粒子的数量，以使映射中的最大错配为原子结构中每10个非氢原子±1个非氢原子；
完全分支的碎片通常应使用较小尺寸的珠子（合理的理由是分支基团的中心原子被掩埋，也就是说，它不暴露于环境中，从而减少了其对相互作用的影响）；例如，一个新戊烷基团包含5个非氢原子，但是由于它是完全支链的，因此可以安全地将其建模为规则的珠子。
注意，在Martini 3.0中，共轭的、厚原子的结构最好由Tiny（T）粒子描述。这样可以确保与堆积相关的特性与全原子数据一致[1]。因此，在这种情况下，萘部分的10个碳原子映射到5个T珠，如下图所示：
剩下的酰胺基团，映射为S粒子，因为S粒子的大小适合描述3个非氢原子，且酰胺不宜被切断。粒子在图中也已编号，以供进一步参考。
3）从原子模拟生成CG映射轨迹
使用刚才创建的映射，现在将在1)中进行的模拟转换为CG分辨率。一种方法是创建一个Gromacs（AA到CG）索引文件，其中每个索引组代表一个珠子，并包含映射的原子数。
不用从头开始手动创建索引文件，而是可以使用CGbuilder工具[3]获得AA到CG的索引文件。只需加载分子的原子pdb / gro文件，单击要第一个粒子包含的原子。如果改变主意，则再次单击以将其删除，单击“新珠子”按钮创建下一个珠子；完成后下载文件。实际上，该工具还允许根据所选的映射获得珠子的初始CG配置（gro文件）和CG-AA映射文件（map文件）。
注意，CGbuilder要求原子不能以不同于1的权重贡献给特定的粒子，而在将原子结构映射到Martini时有时会需要这样做。在这种情况下，索引和/或映射文件应随后手动完善。
在开始之前：关于Martini 2.x的一个重要变化是，在将原子结构映射到CG分辨率时，现在要考虑氢原子以确定珠子的相对位置。这应该反映在AA到CG索引文件中，也就是说，索引还应该包含氢（对CGbuilder而言，就是也要单击氢！）。一般规则是将某个氢原子映射到含有与其相连的非氢原子的珠子上。
现在，可以尝试通过CGbuilder映射NAP.gro。
映射后的结构文件如下，重命名为NAP_CG.gro。
Generated with cgbuilder
6
1NAP BB 1 -0.358 -0.002 -0.046
1NAP SC1 2 -0.088 -0.095 -0.003
1NAP SC2 3 -0.052 0.206 0.027
1NAP SC3 4 0.111 0.012 0.002
1NAP SC4 5 0.273 -0.183 -0.021
1NAP SC5 6 0.349 0.099 0.004
10 10 10索引文件如下，重命名为mapping.ndx。
[ BB ]
21 22 20 17 18
[ SC1 ]
10 5 19
[ SC2 ]
2 1 11 12
[ SC3 ]
4 3
[ SC4 ]
14 6 7 15
[ SC5 ]
16 8 9 13将mapping.ndx文件复 制到当前目录。
现在，我们考虑了氢原子，因为对AA结构的基于几何中心（COG）的映射，是在Martini 3.0中获得成键参数的标准程序[1]。因此，在将AA结构映射到CG分辨率时，我们需要考虑氢原子。由于gmx traj的异常行为（潜在的错误，请参见注释[4]），如果我们要坚持使用gmx traj（替代方法包括使用MDAnalysis Python库），我们需要在运行gmx traj之前进行一些修改，即我们首先需要创建一个原子结构相同，但原子质量全部相等（e.g.，全为1）的全原子的tpr文件。
复 制NAP_GMX.itp为NAP_COG.itp，并将原子质量全部改为1。
NAP_COG.itp文件内容（部分）：
[ atoms ]
; nr type resi res atom cgnr charge mass ; qtot bond_type
1 ca 1 UNK C4 1 -0.110000 1.00000 ; qtot -0.110
2 ca 1 UNK C5 2 -0.122000 1.00000 ; qtot -0.232
3 ca 1 UNK C6 3 -0.023000 1.00000 ; qtot -0.255
4 ca 1 UNK C11 4 -0.050000 1.00000 ; qtot -0.305
5 ca 1 UNK C12 5 -0.054000 1.00000 ; qtot -0.359
6 ca 1 UNK C10 6 -0.108000 1.00000 ; qtot -0.467
7 ca 1 UNK C9 7 -0.129000 1.00000 ; qtot -0.596
8 ca 1 UNK C8 8 -0.117000 1.00000 ; qtot -0.713
9 ca 1 UNK C7 9 -0.124000 1.00000 ; qtot -0.837
10 ha 1 UNK H12 10 0.157000 1.00000 ; qtot -0.680
11 ha 1 UNK H4 11 0.134000 1.00000 ; qtot -0.546
12 ha 1 UNK H5 12 0.136000 1.00000 ; qtot -0.410
13 ha 1 UNK H7 13 0.135000 1.00000 ; qtot -0.275
14 ha 1 UNK H10 14 0.138000 1.00000 ; qtot -0.137
15 ha 1 UNK H9 15 0.136000 1.00000 ; qtot -0.001
16 ha 1 UNK H8 16 0.134000 1.00000 ; qtot 0.133
17 c 1 UNK C 17 0.670701 1.00000 ; qtot 0.804
18 o 1 UNK O 18 -0.611101 1.00000 ; qtot 0.193
19 ca 1 UNK C3 19 -0.144600 1.00000 ; qtot 0.048
20 n 1 UNK N 20 -0.673001 1.00000 ; qtot -0.625
21 hn 1 UNK H1 21 0.312500 1.00000 ; qtot -0.313
22 hn 1 UNK H2 22 0.312500 1.00000 ; qtot 0.000复 制topol.itp为topol_COG.itp，并改为引入NAP_COG.itp。
topol_COG.itp文件内容（部分）：
; NAP_GMX.top created by acpype (Rev: 0) on Fri Jan 8 22:17:10 2021
#include "amber99s b-ildn.ff/forcefield.itp"
#include "amber99s b-ildn.ff/spc.itp"
; Include NAP_GMX.itp topology
#include "NAP_COG.itp"
[ system ]
NAP in water
[ molecules ]
; Compound nmols
NAP 1
SOL 1692准备好上述文件后，执行命令：
echo 0 | gmx trjconv -f md.xtc -o AA-traj.whole.xtc -s md.tpr -pbc whole
gmx grompp -f run.mdp -c eq.gro -p topol_COG.top -o AA-COG.tpr
rm mdout.mdp # clean-up
# 因为要划分为6个粒子, 所以用seq 0 5生成0到5的序列, 且指定-ng 6
seq 0 5 | gmx traj -f AA-traj.whole.xtc -s AA-COG.tpr -oxt mapped.xtc -n mapping.ndx -ng 6 -com这将保证：
首先，确保AA轨迹是完整的，即目标分子在轨迹文件的一个或多个帧中未被周期性边界条件分割（
gmx trjconv -pbc whole ...命令）；随后创建一个
AA-COG.tpr，将在接下来的步骤（gmx grompp -p ...命令）中用于COG映射；最后，将AA轨迹映射到CG分辨率：
gmx traj -f ...命令将进行COG映射，因为它使用的是AA-COG.tpr。
4）创建初始CGitp和tpr文件
CGitp文件的创建必须手动完成，尽管从现有itp文件中进行一些复 制粘贴可能有助于正确设置格式。您（可能）需要的不同指令是：
[ moleculetype ]：一行包含分子名称和许多排除项。对于Martini，排除的标准数量为1。
[ atoms ]每个原子一行，含有原子编号，粒子类型，残基编号，残基名，原子名，电荷组，电荷，质量。小分子的残基编号和残基名称将相同。原子名称可以自由选择。在Martini，每个粒子都有自己的电荷组。通常不指定质量，在这种情况下，质量取自粒子定义。
[ bonds ]每个键一行，包含原子1，原子2，函数类型，平均长度，作用力常数。Martini的函数形式为1。可以将键长设置为在步骤5中获得的平均长度。应将力常数根据所获得的直方图的宽度来调整：较小的宽度表示较高的力常数，反之亦然（另请参见以下各节）。
[ constraints ]，[ angles ]，[ dihedrals ]，[ exclusions ] 的正确格式见Gromacs手册。
构建初步的CG itp文件，名为NAP_initial.itp。实际上，您甚至不需要真正关心成键项。查看第5步的开头，以了解需要哪些成键项。
NAP-initial.itp的文件内容如下：
;;;;;; 2-NAPHTHAMIDE
[ moleculetype ]
; molname nrexcl
NAP 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 SC2 0 NAP BB 1 0
2 TC4 0 NAP SC1 2 0 45
3 TC4 0 NAP SC2 3 0 45
4 TC4 0 NAP SC3 4 0 0
5 TC4 0 NAP SC4 5 0 45
6 TC4 0 NAP SC5 6 0 45
[ virtual_sitesn ]
; site funct constructing atom indices
4 1 2 3 5 6在继续下一步之前，我们需要一个CG tpr文件来从映射的轨迹生成键，角度和二面角的分布。为此，应再准备一个topol_CG.top，文件内容如下：
#include "martini_v3.0.b.3.2.itp"
#include "NAP_initial.itp"
[ system ]
One molecule
[ molecules ]
NAP 1编写用于CG的mdp文件，名为martini.mdp：
;
; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
; Updated 15 Jul 2015 by DdJ
; Adaptation to em/NVT/etc. by RA
;
; for use with GROMACS 5
; For a thorough comparison of different mdp options in combination with the Martini force field, see:
; D.H. de Jong et al., Martini straight: boosting performance using a shorter cutoff and GPUs, submitted.
title = Martini ('new' parameters for) NPT run
; TIMESTEP IN MARTINI
; Most simulations are numerically stable with dt=40 fs,
; however better energy conservation is achieved using a
; 20-30 fs timestep.
; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.
integrator = md
dt = 0.02
nsteps = 1000000 ; 30 ns
nstcomm = 100
comm-grps =
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 100
nstxout-compressed = 1000
compressed-x-precision = 100
compressed-x-grps = System
energygrps = System
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic interactions in the Martini model are
; not considered to be very accurate to begin with, especially as the
; screening in the system is set to be uniform across the system with
; a screening constant of 15. When using PME, please make sure your
; system properties are still reasonable.
;
; With the polarizable water model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable water itself will perform the explicit screening
; in aqueous environment. In this case, the use of PME is more realistic.
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
pbc = xyz
verlet-buffer-tolerance = 0.005
coulombtype = cutoff
coulomb-modifier = Potential-shift-verlet
rcoulomb = 1.1
epsilon_r = 15; 2.5 (with polarizable water)
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale (V-rescale)
; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling
; constant to 0.1 ps, although with such tight coupling (approaching
; the time step) one can no longer speak of a weak-coupling scheme.
; We therefore recommend a coupling time constant of at least 0.5 ps.
; The Berendsen thermostat is less suited since it does not give
; a well described thermodynamic ensemble.
;
; Pressure can be controlled with the Parrinello-Rahman barostat,
; with a coupling constant in the range 4-8 ps and typical compressibility
; in the order of 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,
; the Berendsen barostat probably gives better results, as the Parrinello-
; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
; coupling should be done semiisotropic.
tcoupl = v-rescale
tc-grps = System
tau_t = 1.0
ref_t = 298
Pcoupl = parrinello-rahman
Pcoupltype = isotropic
tau_p = 12.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 3e-4
ref_p = 1
gen_vel = no
gen_temp = 298
gen_seed = 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = none
constraint_algorithm = Lincs最后在官网下载martini.itp（当前名为martini_v3.0.b.3.2.itp），然后运行脚本：
seq 0 5 | gmx traj -f AA-traj.whole.xtc -s AA-COG.tpr -n mapping.ndx -oxt molecule.gro -ng 6 -com -e 0
gmx grompp -f martini_v2.x_new_run.mdp -c molecule.gro -p topol_CG.top -o CG.tpr -maxwarn 1这将：
从轨迹中提取一帧（当然将其映射）；
使用该帧，以及top和mdp（请参阅网站上的示例）文件，为分子创建CG.tpr文件。
遵循《Martini 3.0 圣经》[1]，您还可以根据化学构造单元指定粒子类型。有关粒子选择的详细讨论，请参阅本教程的最后一节。
5）从CG映射轨迹生成目标CG分布
我们需要从步骤3映射到CG的原子模拟中，测量我们想要在CG模型中使用的成键相互作用的长度/角度（bonds, constraints, angles, proper and improper dihedrals）。但是，我们需要定义哪些成键项？让我们回到绘图表中，确定应该在哪些粒子之间定义成键相互作用。
5.1）关于CG模型的成键项的选择
连接2-萘酰胺内的T粒子的键很可能非常刚性，也就是说，它们的分布将非常集中。这要求使用constraints[1]。将粒子模型保持在一起的“幼稚”方法是把所有粒子都约束起来（下图A）：
但是，由于约束算法要满足越来越多的连接约束，这种模型容易出现数值不稳定性。另一种选择是构建“铰链”模型，其中通过5个约束将4个外部粒子（图中的粒子2、3、5和6）连接起来，以形成“铰链”构造，而中心粒子（编号 4）被描述为虚拟位点（Virtual Site, VS），即位置完全由其构造粒子定义的粒子（上图B）。使用虚拟粒子不仅可以提高模型的数值稳定性，还可以提高性能。由于VS没有质量，4号粒子的质量应分摊给四个构造粒子，这样，粒子2、3、5和6的质量应分别为45（= 36 + 36/4）。
然后可以通过两个键即1-2和1-3连接粒子1。需要两个improper项（1-3-6-2和2-1-3-4），以将粒子1保持在萘环上的正确位置。最后还需要使用improper项 2-3-6-5来保持萘环平坦。在这种情况下，所有T粒子之间的排斥也应适用，因为分子非常硬并且在这种情况下不需要分子内相互作用。
5.2）索引文件和目标分布的生成
一旦确定了成键项，就为每个键创建一个带有索引[bondX]的键索引文件，其中包含成对的CG粒子，例如：
[bond1]
1 2
[bond2]
1 4
...角度（带有三个CG粒子）和二面角（带有四个）的情况类似。编写脚本，为感兴趣的所有键，角度和二面体生成分布。对于2-萘酰胺，有七个键（5个约束和2个键）和三个二面角，如前所述。
教程提供了5_generate_target_distr.sh脚本，可通过运行
bash 5_generate_target_distr.sh自动调用轨迹，生成分布。如果要用于新的体系，记得修改其中的gmx命令。
检查文件夹bond_mapped和dihedrals_mapped的结果，会发现每个成键分布为bonds_mapped / distr_bond_X.xvg，而bonds_mapped / data_bonds.txt 给出了每个键的均值和标准偏差。
对于每个键，脚本使用以下命令（在此示例中，该命令应用于索引为0的第一个键）：
echo 0 | gmx distance -f mapped.xtc -n bonds.ndx -s CG.tpr -oall bonds_mapped/bond_0.xvg -xvg none
gmx analyze -f bonds_mapped/bond_0.xvg -dist bonds_mapped/distr_bond_0.xvg -xvg none -bw 0.001对二面角也类似：
echo 0 | gmx angle -type dihedral -f mapped.xtc -n dihedrals.ndx -ov dihedrals_mapped/dih_0.xvg
gmx analyze -f dihedrals_mapped/dih_0.xvg -dist dihedrals_mapped/distr_dih_0.xvg -xvg none -bw 1.06）创建CG模拟
现在，我们可以完成CG模型的第一阶段NAP_take1.itp的确定，在该模型中，我们可以使用data_bonds.txt和data_dihedrals.txt文件中包含的信息来更好地猜测结合参数。
复 制NAP_initial.itp为NAP_take1.itp。使用上一步的输入来调整ENAP_take1.itp。
;;;;;; 2-NAPHTHAMIDE
[ moleculetype ]
; molname nrexcl
NAP 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 SC2 0 NAP BB 1 0
2 TC4 0 NAP SC1 2 0 45
3 TC4 0 NAP SC2 3 0 45
4 TC4 0 NAP SC3 4 0 0
5 TC4 0 NAP SC4 5 0 45
6 TC4 0 NAP SC5 6 0 45
[bonds]
; i j funct length
1 2 1 0.28079 25000 ; cog
1 3 1 0.38768 25000 ; cog
#ifndef FLEXIBLE
[constraints] ; For minimizations, you may use define=-DFLEXIBLE to use a stiff-bond version of the topology that is more amenable to minimization.
#endif
; i j funct length
2 3 1 0.30060 1000000 ; cog
2 5 1 0.36661 1000000 ; cog
2 6 1 0.47223 1000000 ; cog
3 6 1 0.41003 1000000 ; cog
5 6 1 0.29058 1000000 ; cog
[ dihedrals ]
; improper
; i j k l funct ref.angle force_k
1 3 6 2 2 0 50 ; controls planarity of amide group
2 1 3 4 2 0 50 ; keeps amide on position 1 of naphthalene
2 3 6 5 2 0 200 ; keeps naphthalene planar
[ virtual_sitesn ]
; site funct constructing atom indices
4 1 2 3 5 6
[ exclusions ]
1 2 3 4 5 6
2 3 4 5 6
3 4 5 6
4 5 6
5 6同时修改topol_CG.top，如下：
#include "martini_v3.0.b.3.2.itp"
#include "martini_v3.0_solvents.itp"
#include "NAP_takel.itp"
[ system ]
One molecule
[ molecules ]
NAP 1执行：
gmx solvate -cp NAP_CG.gro -cs box_CG_WN_eq.gro -p topol_CG.top -o NAP_sol.gro -box 4.5 4.5 4.5
gmx grompp -f martini_em.mdp -c NAP_sol.gro -p topol_CG.top -o 1-min.tpr -po 1-min.mdp -maxwarn 1
gmx mdrun -v -deffnm 1-min
gmx grompp -p topol_CG.top -c 1-min.gro -f martini_eq.mdp -o 2-eq.tpr -po 2-eq.mdp
gmx mdrun -v -deffnm 2-eq
gmx grompp -p topol_CG.top -c 2-eq.gro -f martini_run.mdp -o 3-run.tpr -po 3-run.mdp
gmx mdrun -v -deffnm 3-run将对水中的Martini系统进行能量最小化，然后进行NPT平衡，并执行50 ns的MD产生相模拟（检查脚本和各种mdp文件以了解更多信息）。
检查所生成的结构文件和轨迹文件。
一旦运行了MD，就可以使用为映射轨迹生成的索引文件来生成CG轨迹的分布：
cp ../5_target-distr/bonds.ndx .
cp ../5_target-distr/dihedrals.ndx .
bash 6_generate_CG_distr.sh它将生成一系列文件，与上一步中的5_generate_target_distr.sh所做的相似，但现在针对CG轨迹。
现在，可以绘制使用脚本draw_bond_distribution.py和draw_dihedral_distribution.py绘制二者的分布并进行比较。
发现萘环上的键长分布与全原子吻合较好，而1-3和1-4（前两张子图）分布不太合理。
观察模拟结构，发现酰胺基团和3号原子（SC2）距离较近，不合理。因此在NAP_take1.itp中添加键角约束。
添加键角约束后重新跑一次能量最小化、平衡相和产生相模拟，再计算键长分布，如下图。
可以看到两者分布比上一次模拟合理一些。
同理计算二面角分布，如下图。
注意，有时分布会出现两个峰，而CG模型无法考虑双峰性。如果在第一次迭代中不令人满意（很可能会发生），则应使用CG itp中的平衡值和力常数进行迭代，直到达成令人满意的协议。
8）与实验结果进行比较，进一步完善和最终考虑
粗粒化后的模拟结果应从划分自由能、纯液体或有机晶体的密度等进行合理性验证，以及其它更具体的参考材料。
（注：本文章来源于网络,如有侵权,请联系我们删除。）
