氧合血红素cpdI分子动力学模拟（Amber + Gromacs）
一、配体结构获取：
- 配体结构由iqmol或者高斯view手动绘制；
- 绘制好后另存为pdb或mol2文件待用；
- 可根据RCSB网站直接下载；
二、血红素-配体复合物分子对接：
1. 对接软件选择：
- 薛定谔(收费，无授权，自行解决)
- chai_lab（免费，需要科学上网，要排队，在线版，普通邮箱注册即可，需要自己绘制SMILES字符）
- SMILES字符可以通过iqmol另存为
- swisstargetprediction网站绘制（更推荐）
- chai_lab预测结果为cif文件，需要转为pdb
- chai_lab预测结果原子名非常不对，需要手动修改，尤其是配体！！！
- autodock(群里有对接教程，请自行参阅)
- CBDOCK（在线网页版，结果不太可靠，仅供参考）
- 以上仅为个人推荐使用的软件，其余诸如HDOCK,ZDOCK等需要教育邮箱的不做推荐（本人没有）
2. 结果整理：
- 选择合适的理想构象作为输出结果；
- 检查输出结果有无化学键的断裂或者残基缺失等问题；
三、模拟体系构建
1. 配体处理:配体分子名必须始终保持一致！！！
- 配体需要使用amber工具antechamber 将小分子pdb文件转为mol2文件（必须用该工具！！！）
- 启动ambertools:
source ~/scripts/evamber.sh -
antechamber -i Y.pdb -fi pdb -o Y.mol2 -fo mol2 -c bcc -s 2
- 获得mol2文件后，使用iqmol打开，添加H原子，并且使用MMFF94s力场优化之后保存为xyz文件用于计算RESP2电荷
- ORCA和Multiwfn计算RESP2电荷：
- 使用文本阅读器(notepad++,emeditor,vscode)打开，复制坐标文件；
- 打开lig_opt.inp文件，将坐标文件复制在**之间，同时注意检查电荷，如果不带电，将 xyz 1 1 改为 xyz 0 1
- shift + 右键工作文件夹，点击Cmder Here,输入
E:\software\OCRA\orca.exe lig_opt.inp > lig_out.inp - 上述计算完成后会输出lig_out.xyz文件，计算未完成会有一堆临时文件，请不要做任何给改动
- 等待计算结果完成后，点开lig_sp1.inp和lig_sp2.inp，将lig_out.xyz文件坐标以同样方式粘贴近上述两个inp文件中
- 在Cmder界面依次输入
E:\software\OCRA\orca.exe lig_sp1.inp > lig_sp1.inpE:\software\OCRA\orca.exe lig_sp2.inp > lig_sp2.inp一个计算完再输入另外一个！！！计算完成后请继续输入以下命令：E:\software\OCRA\orca_2mkl.exe lig_sp1 -moldenE:\software\OCRA\orca_2mkl.exe lig_sp2 -molden - 上述命令结束后会输出.input文件；同样在cmder界面输入
E:\software\Multwfn\Multiwfn_3.8_dev_bin_Win64\Multiwfn_3.8_dev_bin_Win64\Multiwfn.exe将.input文件拖入cmder界面，按下回车 - 随后依次输入：7 #Population analysis and caculation of atomic charges18 #Restrained ElectroStatic Poteintial atomic charge1 #Start standard two-stage RESPfitting caculationy #yes
- 将得到的.chg文件最后一列拖入excel中，计算总电荷数，是不是和.inp文件初始设定保持一致(1E-10也算一致，gmx和amber单精度只有1e-6),随后每个原子电荷取两个.chg文件的平均值，替换mol2文件的电荷即可
- 得到lig.mol2文件后，使用parmchk2 产生小分子参数文件：
parmchk2 -i lig.mol2 -f mol2 -o lig.frcmod
- 使用tleap程序生成小分子lib文件：
1.vim lig.tleap ##复制以下内容 source leaprc.gaff lig = loadmol2 lig.mol2 check lig loadamberparams lig.frcmod check lig saveoff lig lig.lib saveamberparm lig lig.prmtop lig.inpcrd quit #按下esc，再输入:wq!按下回车保存并退出 :wq! #运行tleap tleap -f lig.tleap #没有报错会输出lig.lib文件
- 启动ambertools:
2. (cpdI)HEM的处理：
- 从对接结果中提取cpdI的坐标文件，另存为cpdi.pdb(注意分子名改为HEM);
- 使用ambertools:
antechamber -i hem.pdb -fi pdb -o hem.mol2 -fo mol2 -c bcc -s 2 - 得到mol2文件之后，请根据原子名将hem.mol2的坐标复制到附件一对应的原子坐标，另存为hem_ref.mol2；
- 此时你得文件夹内应该有:hem_ref.mol2hem.pdbsub.mol2sub.pdbsub.frcmodsub.lib等文件！！！如果没有，请仔细检查！！！
3. 受体蛋白处理：
- 使用文本编辑器打开X.pdb文件，删掉除原子信息外的所有行，只保留ATOM字段，另存为receptor.pdb
- 封端处理：
vim fengduan.tleap
#按下a，输入以下内容
source leaprc.protein.ff19SB # 或者根据你的蛋白质选择适合的力场
source leaprc.water.tip3p # 如果你还需要考虑水分子
protein = loadpdb "receptor.pdb" #加载受体蛋白文件
addions protein ACE #为N-端添加ACE(乙酰基)
addions protein NME #为C-端添加NME(甲基氨基乙酰基)
savepdb protein "receptor_termini.pdb" #保存修改后结构
quit #退出
:wq!保存并退出vim编辑器
tleap -f fengduan.tleap
- 使用pdb4amber对受体蛋白重新编码加H：
pdb4amber -i receptor_tgermini.pdb -o receptor_tgermini_reduce.pdb -y --reduce
#-i：输入蛋白质的pdb格式文件
#-y：该参数是指去除所有氢原
#--reduce ： 运行Reduce来加氢
#之所以用了-y和–reduce选项，是因为pdb4amber会检查每个组氨酸（HIE、HIP、HID）的类型
#如果用其他软件加氢则不一定是amber所支持的原子类型，所以还是用pdb4amber一劳永逸吧
#最后，对保存的receptor_tgermini_reduce.pdbpdb手动修改组氨酸质子化状态的分子名称
- HEM链接的CYS处理：
- 将HEM链接的CYS的分子名改为CYP；
- 将自己的cyp坐标替换附件二的CYP.mol2；
- 复合物体系搭建：
- 复制HEM.pdb,sub.pdb到receptor_tgermini_reduce.pdb末尾；
- 将原子序号和分子序号修改，按序顺延！！！
- 请务必仔细检查hem和sub的分子名是否和pdb，mol2文件一一对应！！！
- 另存为receptor_tgermini_reduce_complex.pdb
- 创建复合物体系拓扑参数文件：
vim compplex.tleap #请务必使用该命令，否则tleap程序不识别dos文件格式！！！
###以下为需要黏贴的内容，上面为启动编辑器的命令
source leaprc.protein.ff19SB #加载19SB力场（amber20以上）
source leaprc.gaff2 #加载小分子gaff2力场（amber20以上）
source leaprc.water.tip3p #加载水分子模型TIP3P力场
CYP = loadmol2 cyp.mol2 #加载和HEM链接的CYS分子
HEM = loadmol2 heme_ref.mol2 #加载氧合血红素分子
SUB = loadmol2 sub_h.mol2 #加载配体分子
loadamberparams frcmod.ions1lm_126_tip3p #加载水分子参数
loadamberparams CPDI.frcmod #加载氧合cpdI力场参数（包括FE-S成键项，不用了另外加载CYP.frcmod）
loadamberparams sub_h.frcmod #加载配体参数
mol = loadpdb receptor_tgermini_reduce_complex.pdb #加载受体分子
bond mol.43.SG mol.242.FE #链接HEM的FE原子和CYP的SG原子
bond mol.42.C mol.43.N #链接CYP的上游残基
bond mol.43.C mol.44.N #链接CYP下游残基
savepdb mol complex4amber.pdb #保存为complex4amber.pdb
saveamberparm mol test.prmtop test.inpcrd #输出复合物单独的力场参数文件
solvatebox mol TIP3PBOX 12.0 #创建模拟盒子，蛋白距盒边12.0埃大小
addions mol Na+ 0 #自动添加Na+中和电性
addions mol Cl- 0 #自动添加Cl-中和电性
savepdb mol solv.pdb #保存为solv.pdb
check mol #检查体系是否正常
saveamberparm mol solv.prmtop solv.inpcrd #输出体系力场参数文件
quit
###上述为complex.tleap文件内容
###运行命令
tleap -f complex.tleap
### 运行成功即可出现solv.prmtop solv.inpcrd
###amber转gromacs
###方法1：parmed方法
pthon3
import parmed as pmd
amber=pmd.load_file('complex_ions.prmtop','complex_ions.inpcrd')
amber.save('complex_ions.top')
amber.save('complex_ions.gro')
###方法2：amb2gro_top_gro.py
amb2gro_top_gro.py -p solv.prmtop -c solv.inpcrd -t solv.top -g solv.gro -b solv_gromacs.pdb
###如果上述两个方法出现下列报错：
### RuntimeError:Could not find <Atom NC[3371];In HEM 216>
###建议使用方法3
###方法3：使用acpepe转化
acpype -p solv.prmtop -x solv.inpcrd -b solv
###最后，如果您在 Amber 中添加 Na+ 和 Cl- 离子
###则它们可以在 solv_GMX.gro 文件中按随机顺序排列
###检查您的 .gro 文件，看看 Na+ 和 Cl- 是否乱序
###GROMACS 需要将 Na+ 和 Cl- 全部分组到坐标文件中。如果需要对离子重新排序
###可以使用脚本reorder.py详见附件四
四、氧合血红素结构，参数，脚本获取：
CPDI对接结构手绘，软件选择高斯view或者iqmol，不想手绘可以参考以下结构：
不建议直接使用文献的mol2文件对接，原子类型不对，会直接报错！！！
HETATM 3787 O2D HEM A 247 10.103 -15.394 -0.015 1.00 12.46 O1-
HETATM 3788 CGD HEM A 247 9.192 -14.965 0.074 1.00 10.92 C
HETATM 3789 O1D HEM A 247 8.866 -14.797 1.265 1.00 11.23 O
HETATM 3790 CBD HEM A 247 8.225 -14.433 -0.977 1.00 9.51 C
HETATM 3791 CAD HEM A 247 8.532 -15.002 -2.313 1.00 8.86 C
HETATM 3792 C3D HEM A 247 7.615 -14.498 -3.449 1.00 8.32 C
HETATM 3793 C2D HEM A 247 6.447 -15.029 -3.902 1.00 8.44 C
HETATM 3794 CMD HEM A 247 5.755 -16.244 -3.360 1.00 8.65 C
HETATM 3795 C4D HEM A 247 7.856 -13.294 -4.287 1.00 7.77 C
HETATM 3796 CHA HEM A 247 8.966 -12.476 -4.106 1.00 7.75 C
HETATM 3797 ND HEM A 247 6.916 -13.152 -5.248 1.00 8.16 N
HETATM 3798 C1D HEM A 247 6.025 -14.169 -5.040 1.00 8.79 C
HETATM 3799 CHD HEM A 247 4.872 -14.398 -5.727 1.00 9.01 C
HETATM 3800 FE HEM A 247 6.852 -11.835 -6.664 1.00 7.90 Fe2+
HETATM 3801 NB HEM A 247 6.688 -10.348 -8.077 1.00 7.99 N
HETATM 3802 C4B HEM A 247 5.728 -10.336 -9.010 1.00 8.02 C
HETATM 3803 C3B HEM A 247 6.055 -9.190 -9.894 1.00 8.47 C
HETATM 3804 CAB HEM A 247 5.208 -8.828 -11.054 1.00 9.41 C
HETATM 3805 CBB HEM A 247 5.595 -8.072 -12.099 1.00 11.75 C
HETATM 3806 C2B HEM A 247 7.235 -8.637 -9.396 1.00 8.22 C
HETATM 3807 CMB HEM A 247 7.951 -7.439 -9.984 1.00 8.95 C
HETATM 3808 C1B HEM A 247 7.617 -9.402 -8.269 1.00 7.67 C
HETATM 3809 CHB HEM A 247 8.720 -9.106 -7.505 1.00 8.13 C
HETATM 3810 NC HEM A 247 5.031 -12.558 -7.344 1.00 8.11 N1-
HETATM 3811 C4C HEM A 247 4.425 -13.628 -6.757 1.00 8.73 C
HETATM 3812 C3C HEM A 247 3.216 -13.875 -7.445 1.00 9.41 C
HETATM 3813 CAC HEM A 247 2.362 -15.046 -7.114 1.00 10.67 C
HETATM 3814 CBC HEM A 247 1.544 -15.542 -8.006 1.00 14.62 C
HETATM 3815 C2C HEM A 247 3.126 -12.917 -8.450 1.00 8.86 C
HETATM 3816 CMC HEM A 247 2.065 -12.752 -9.527 1.00 9.63 C
HETATM 3817 C1C HEM A 247 4.284 -12.134 -8.360 1.00 8.64 C
HETATM 3818 CHC HEM A 247 4.622 -11.077 -9.169 1.00 7.77 C
HETATM 3819 NA HEM A 247 8.434 -10.907 -5.938 1.00 7.63 N1-
HETATM 3820 C1A HEM A 247 9.191 -11.386 -4.883 1.00 7.92 C
HETATM 3821 C4A HEM A 247 9.099 -9.832 -6.402 1.00 8.04 C
HETATM 3822 C3A HEM A 247 10.232 -9.558 -5.584 1.00 7.80 C
HETATM 3823 CMA HEM A 247 11.217 -8.439 -5.757 1.00 7.11 C
HETATM 3824 C2A HEM A 247 10.326 -10.487 -4.658 1.00 7.87 C
HETATM 3825 CAA HEM A 247 11.351 -10.653 -3.545 1.00 9.00 C
HETATM 3826 CBA HEM A 247 10.779 -10.241 -2.144 1.00 8.86 C
HETATM 3827 CGA HEM A 247 10.399 -8.791 -2.081 1.00 10.30 C
HETATM 3828 O1A HEM A 247 9.216 -8.413 -1.667 1.00 11.32 O
HETATM 3829 O2A HEM A 247 11.235 -7.949 -2.496 1.00 10.09 O1-
HETATM 3830 O1 HEM A 247 5.972 -10.851 -5.459 1.00 0.00 O
HETATM 3831 HBD1 HEM A 247 8.301 -13.344 -1.004 1.00 0.00 H
HETATM 3832 HBD2 HEM A 247 7.197 -14.675 -0.699 1.00 0.00 H
HETATM 3833 HAD1 HEM A 247 8.457 -16.091 -2.290 1.00 0.00 H
HETATM 3834 HAD2 HEM A 247 9.563 -14.797 -2.603 1.00 0.00 H
HETATM 3835 HMD1 HEM A 247 4.689 -16.054 -3.238 1.00 0.00 H
HETATM 3836 HMD2 HEM A 247 5.878 -17.084 -4.043 1.00 0.00 H
HETATM 3837 HMD3 HEM A 247 6.150 -16.535 -2.387 1.00 0.00 H
HETATM 3838 HHA HEM A 247 9.650 -12.704 -3.304 1.00 0.00 H
HETATM 3839 HHD HEM A 247 4.296 -15.250 -5.408 1.00 0.00 H
HETATM 3840 HAB HEM A 247 4.200 -9.192 -11.076 1.00 0.00 H
HETATM 3841 HBB1 HEM A 247 4.903 -7.863 -12.900 1.00 0.00 H
HETATM 3842 HBB2 HEM A 247 6.588 -7.663 -12.193 1.00 0.00 H
HETATM 3843 HMB1 HEM A 247 8.654 -6.978 -9.294 1.00 0.00 H
HETATM 3844 HMB2 HEM A 247 8.509 -7.733 -10.873 1.00 0.00 H
HETATM 3845 HMB3 HEM A 247 7.237 -6.662 -10.257 1.00 0.00 H
HETATM 3846 HHB HEM A 247 9.346 -8.273 -7.777 1.00 0.00 H
HETATM 3847 HHC HEM A 247 2.399 -15.483 -6.130 1.00 0.00 H
HETATM 3848 HBC1 HEM A 247 0.955 -16.406 -7.760 1.00 0.00 H
HETATM 3849 HBC2 HEM A 247 1.444 -15.152 -9.004 1.00 0.00 H
HETATM 3850 HMC1 HEM A 247 1.992 -11.726 -9.866 1.00 0.00 H
HETATM 3851 HMC2 HEM A 247 2.287 -13.384 -10.387 1.00 0.00 H
HETATM 3852 HMC3 HEM A 247 1.078 -13.006 -9.142 1.00 0.00 H
HETATM 3853 HHC HEM A 247 3.955 -10.867 -9.984 1.00 0.00 H
HETATM 3854 HMA1 HEM A 247 11.981 -8.381 -4.986 1.00 0.00 H
HETATM 3855 HMA2 HEM A 247 11.739 -8.564 -6.704 1.00 0.00 H
HETATM 3856 HMA3 HEM A 247 10.697 -7.481 -5.779 1.00 0.00 H
HETATM 3857 HAA1 HEM A 247 12.275 -10.120 -3.769 1.00 0.00 H
HETATM 3858 HAA2 HEM A 247 11.656 -11.699 -3.514 1.00 0.00 H
HETATM 3859 HBA1 HEM A 247 11.552 -10.418 -1.398 1.00 0.00 H
HETATM 3860 HBA2 HEM A 247 9.931 -10.873 -1.877 1.00 0.00 H
附件一：cpdI参考结构及电荷文件cpdI.mol2
# Shahrokh,K; Orendt,A; Yost, G.S; and Cheatham III, T.E. Jour Comp Chem (2011)
# Coumpound I High-spin heme coordinates and atomic charges
@<TRIPOS>MOLECULE
HEM-CPDI
74 81 1 0 1
SMALL
USER_CHARGES
@<TRIPOS>ATOM
1 NC 2.870000 -0.982000 -0.754000 nc 1 HEM -0.0333 ****
2 C1C 4.206000 -0.781000 -0.528000 cc 1 HEM -0.0836 ****
3 C4C 2.714000 -2.307000 -1.022000 cc 1 HEM -0.0148 ****
4 C2C 4.934000 -2.026000 -0.645000 cc 1 HEM -0.0025 ****
5 C3C 3.998000 -2.997000 -0.943000 cc 1 HEM -0.1636 ****
6 CHD 1.504000 -2.923000 -1.333000 cg 1 HEM 0.0470 ****
7 HHD 1.518000 -3.994000 -1.504000 ha 1 HEM 0.0773 ****
8 C1D 0.275000 -2.292000 -1.459000 cd 1 HEM -0.1256 ****
9 ND 0.033000 -0.964000 -1.253000 nd 1 HEM 0.0685 ****
10 C4D -1.296000 -0.768000 -1.518000 cd 1 HEM -0.0340 ****
11 C3D -1.935000 -2.029000 -1.895000 cd 1 HEM -0.0325 ****
12 C2D -0.948000 -2.976000 -1.850000 cd 1 HEM 0.0288 ****
13 CHA -1.955000 0.450000 -1.421000 cg 1 HEM -0.0124 ****
14 HHA -3.015000 0.459000 -1.641000 ha 1 HEM 0.1920 ****
15 C1A -1.380000 1.648000 -1.011000 cc 1 HEM -0.0365 ****
16 C2A -2.100000 2.914000 -0.894000 cc 1 HEM -0.0626 ****
17 C4A 0.084000 3.125000 -0.329000 cc 1 HEM -0.0543 ****
18 NA -0.069000 1.806000 -0.659000 nc 1 HEM -0.0015 ****
19 C3A -1.181000 3.829000 -0.460000 cc 1 HEM 0.0097 ****
20 CHB 1.271000 3.739000 0.037000 cg 1 HEM 0.0210 ****
21 HHB 1.216000 4.793000 0.285000 ha 1 HEM 0.1004 ****
22 C1B 2.526000 3.138000 0.092000 cd 1 HEM -0.1156 ****
23 C2B 3.750000 3.834000 0.427000 cd 1 HEM 0.0805 ****
24 NB 2.778000 1.824000 -0.184000 nd 1 HEM 0.0339 ****
25 C4B 4.122000 1.646000 -0.047000 cd 1 HEM -0.0866 ****
26 CHC 4.792000 0.439000 -0.209000 cg 1 HEM 0.0569 ****
27 HHC 5.861000 0.438000 -0.029000 ha 1 HEM 0.0740 ****
28 FE 1.406000 0.441000 -0.796000 fe 1 HEM 0.2620 ****
29 C3B 4.763000 2.899000 0.343000 cd 1 HEM -0.0237 ****
30 CAB 6.194000 3.060000 0.540000 cc 1 HEM -0.0963 ****
31 HAB 6.817000 2.326000 0.027000 ha 1 HEM 0.1169 ****
32 CBB 6.836000 3.995000 1.265000 cd 1 HEM -0.3967 ****
33 HBB1 6.314000 4.743000 1.853000 ha 1 HEM 0.1410 ****
34 HBB2 7.922000 4.014000 1.310000 ha 1 HEM 0.1410 ****
35 CAC 4.181000 -4.419000 -1.186000 cd 1 HEM 0.0915 ****
36 HAC 3.404000 -4.891000 -1.788000 ha 1 HEM 0.0767 ****
37 CBC 5.182000 -5.213000 -0.760000 cc 1 HEM -0.4715 ****
38 HBC1 5.981000 -4.859000 -0.117000 ha 1 HEM 0.1517 ****
39 HBC2 5.207000 -6.266000 -1.032000 ha 1 HEM 0.1517 ****
40 CMB 3.863000 5.294000 0.744000 c3 1 HEM -0.1258 ****
41 HMB1 3.798000 5.494000 1.824000 hc 1 HEM 0.0397 ****
42 HMB2 4.826000 5.693000 0.402000 hc 1 HEM 0.0397 ****
43 HMB3 3.069000 5.874000 0.261000 hc 1 HEM 0.0397 ****
44 CMC 6.418000 -2.183000 -0.508000 c3 1 HEM 0.0353 ****
45 HMC1 6.720000 -2.440000 0.518000 hc 1 HEM 0.0049 ****
46 HMC2 6.786000 -2.983000 -1.160000 hc 1 HEM 0.0049 ****
47 HMC3 6.946000 -1.261000 -0.777000 hc 1 HEM 0.0049 ****
48 CMD -1.057000 -4.446000 -2.127000 c3 1 HEM -0.2123 ****
49 HMD1 -0.334000 -4.774000 -2.886000 hc 1 HEM 0.0653 ****
50 HMD2 -0.870000 -5.047000 -1.226000 hc 1 HEM 0.0653 ****
51 HMD3 -2.057000 -4.704000 -2.489000 hc 1 HEM 0.0653 ****
52 CMA -1.392000 5.286000 -0.173000 c3 1 HEM -0.0747 ****
53 HMA1 -1.304000 5.510000 0.900000 hc 1 HEM 0.0306 ****
54 HMA2 -0.657000 5.915000 -0.693000 hc 1 HEM 0.0306 ****
55 HMA3 -2.388000 5.608000 -0.491000 hc 1 HEM 0.0306 ****
56 CAA -3.554000 3.168000 -1.187000 c3 1 HEM -0.0115 ****
57 HAA1 -3.660000 4.127000 -1.709000 hc 1 HEM 0.0253 ****
58 HAA2 -3.949000 2.419000 -1.880000 hc 1 HEM 0.0253 ****
59 CAD -3.395000 -2.247000 -2.184000 c3 1 HEM -0.0659 ****
60 HAD1 -3.818000 -1.387000 -2.716000 hc 1 HEM 0.0364 ****
61 HAD2 -3.526000 -3.100000 -2.860000 hc 1 HEM 0.0364 ****
62 CBA -4.444000 3.191000 0.074000 c3 1 HEM -0.0304 ****
63 HBA1 -4.048000 3.901000 0.813000 h1 1 HEM -0.0023 ****
64 HBA2 -4.428000 2.209000 0.564000 h1 1 HEM -0.0023 ****
65 CBD -4.217000 -2.482000 -0.900000 c3 1 HEM -0.0022 ****
66 HBD1 -4.053000 -1.651000 -0.207000 h1 1 HEM 0.0140 ****
67 HBD2 -3.864000 -3.386000 -0.384000 h1 1 HEM 0.0140 ****
68 CGD -5.729000 -2.615000 -1.120000 c 1 HEM 0.5749 ****
69 CGA -5.891000 3.581000 -0.216000 c 1 HEM 0.6402 ****
70 O1A -6.681000 3.739000 0.763000 o 1 HEM -0.6508 *****
71 O1D -6.463000 -2.631000 -0.087000 o 1 HEM -0.6552 *****
72 O2A -6.295000 3.767000 -1.392000 o 1 HEM -0.6508 *****
73 O2D -6.192000 -2.733000 -2.281000 o 1 HEM -0.6552 *****
74 O1 1.689000 0.783000 -2.374000 oa 1 HEM -0.3729 ****
@<TRIPOS>BOND
1 1 2 1
2 1 3 1
3 1 28 1
4 2 4 1
5 2 26 1
6 3 5 1
7 3 6 1
8 4 5 1
9 4 44 1
10 5 35 1
11 6 7 1
12 6 8 1
13 8 9 1
14 8 12 1
15 9 10 1
16 9 28 1
17 10 11 1
18 10 13 1
19 11 12 1
20 11 59 1
21 12 48 1
22 13 14 1
23 13 15 1
24 15 16 1
25 15 18 1
26 16 19 1
27 16 56 1
28 17 18 1
29 17 19 1
30 17 20 1
31 18 28 1
32 19 52 1
33 20 21 1
34 20 22 1
35 22 23 1
36 22 24 1
37 23 29 1
38 23 40 1
39 24 25 1
40 24 28 1
41 25 26 1
42 25 29 1
43 26 27 1
44 28 74 1
45 29 30 1
46 30 31 1
47 30 32 1
48 32 33 1
49 32 34 1
50 35 36 1
51 35 37 1
52 37 38 1
53 37 39 1
54 40 41 1
55 40 42 1
56 40 43 1
57 44 45 1
58 44 46 1
59 44 47 1
60 48 49 1
61 48 50 1
62 48 51 1
63 52 53 1
64 52 54 1
65 52 55 1
66 56 57 1
67 56 58 1
68 56 62 1
69 59 60 1
70 59 61 1
71 59 65 1
72 62 63 1
73 62 64 1
74 62 69 1
75 65 66 1
76 65 67 1
77 65 68 1
78 68 71 1
79 68 73 1
80 69 70 1
81 69 72 1
@<TRIPOS>SUBSTRUCTURE
1 HEM 1 **** 0 **** ****
3. 附件二：同HEM链接的CYS参数——CYP.mol2
# Shahrokh,K; Orendt,A; Yost, G.S; and Cheatham III, T.E. Jour Comp Chem (2011)
# Compound I cysteine coordinates and atomic charges
@<TRIPOS>MOLECULE
CYP-CPDI
10 9 1 0 0
SMALL
USER_CHARGES
@<TRIPOS>ATOM
1 N -2.6330 -1.3520 2.7840 N 1 CYP -0.4157
2 H -2.8330 -2.0920 3.4500 H 1 CYP 0.2719
3 CA -1.2870 -1.4470 2.2240 CT 1 CYP 0.0725
4 HA -1.3570 -1.6240 1.1450 H1 1 CYP -0.0095
5 CB -0.4780 -0.1320 2.4210 CT 1 CYP 0.0020
6 HB2 -0.4010 0.0680 3.4970 H1 1 CYP 0.0508
7 HB3 -1.0820 0.6620 1.9770 H1 1 CYP 0.0508
8 SG 1.2250 -0.0810 1.7090 SH 1 CYP -0.4381
9 C -0.6660 -2.7000 2.8810 C 1 CYP 0.5973
10 O -1.2840 -3.3630 3.7280 O 1 CYP -0.5679
@<TRIPOS>BOND
1 1 3 1
2 1 2 1
3 3 4 1
4 3 5 1
5 3 9 1
6 5 6 1
7 5 7 1
8 5 8 1
9 9 10 2
@<TRIPOS>SUBSTRUCTURE
1 CYP 1 RESIDUE 4 A CYP 0 ROOT
附件三：cys和氧合血红素链接的amber力场参数文件——CPDI-HEM.frcmod
CPDI.frcmod- Shahrokh,K; Orendt,A; Yost, G.S; and Cheatham III, T.E. Jour Comp Chem (2011)
MASS
fe 55.85
oa 16.00
BOND
fe-nc 114.00 2.029
fe-nd 114.00 2.029
fe-SH 39.00 2.565
cg-ha 341.50 1.089 # same as ce-ha 341.5 1.089 SOURCE3
fe-oa 572.00 1.639
ANGLE
nc-fe-oa 65.000 92.406 # average angle
nd-fe-oa 65.000 92.406 # average angle
SH-fe-oa 0.000 174.087 #
nc-fe-nd 239.000 89.900 # average angle
fe-nc-cc 146.000 126.651 # average angle
fe-nd-cd 146.000 126.651 # average angle
nc-fe-nc 0.000 174.731
nd-fe-nd 0.000 175.636
SH-fe-nc 48.00 87.595 # average angle
SH-fe-nd 48.00 87.595 # average angle
CT-SH-fe 39.00 105.885
cc-cc-cg 65.6 124.539 # same as cc-cc-cf 65.6 123.92 SOURCE3
cd-cd-cg 65.6 124.539 # same as cc-cc-cf 65.6 123.92 SOURCE3
nc-cc-cg 68.5 125.044 # same as cf-cf-n2 68.5 123.00 SOURCE3
nd-cd-cg 68.5 125.044 # same as cf-cf-n2 68.5 123.00 SOURCE3
ha-cc-ha 38.0 117.65 # same as ha-c2-ha 38.0 117.65 SOURCE3
ha-cd-ha 38.0 117.65 # same as ha-c2-ha 38.0 117.65 SOURCE3
cc-cg-ha 46.6 116.969 # same as cd-cd-ha 46.6 123.74 SOURCE3
cd-cg-ha 46.6 116.969 # same as cd-cd-ha 46.6 123.74 SOURCE3
cd-cg-cc 63.8 126.057 # same as ce-ce-cf 63.8 130.92 SOURCE3
DIHEDRAL
cd-nd-fe-oa 1 0.000 180.000 2.000
cc-nc-fe-oa 1 0.000 180.000 2.000
CT-SH-fe-nc 1 0.00 180.000 2.000
CT-SH-fe-nd 1 0.00 180.000 2.000
CT-SH-fe-oa 1 0.00 180.000 2.000
cc-cc-c3-c3 3 0.000 180.000 2.000
cc-cc-c3-hc 3 0.000 180.000 2.000
cd-cd-c3-c3 3 0.000 180.000 2.000
cd-cd-c3-hc 3 0.000 180.000 2.000
X -cg-cd-X 4 16.00 180.000 2.000
X -cg-cc-X 4 16.00 180.000 2.000
cg-cc-nc-fe 1 0.000 180.000 2.000
cg-cd-nd-fe 1 0.000 180.000 2.000
cc-cc-nc-fe 1 0.000 180.000 2.000
cd-cd-nd-fe 1 0.000 180.000 2.000
nd-fe-nc-cc 1 0.000 180.000 2.000
nc-fe-nd-cd 1 0.000 180.000 2.000
cc-nc-fe-nc 1 0.000 180.000 2.000
cd-nd-fe-nd 1 0.000 180.000 2.000
cc-nc-fe-SH 1 0.000 180.000 2.000
cd-nd-fe-SH 1 0.000 180.000 2.000
IMPROPER
NONBON
fe 1.3 0.01
oa 1.6612 0.21
请注意！！！原始文献的附件中frcmod文件的DIHEDRAL项中：cd-nd-fe-SH为cd-nd-fe-sh，
需要将小写的sh改为大写的SH！！！
否则报错二面角参数不对！！！
附件四：离子重新排序脚本：reorder.py
#!/usr/bin/env python
import sys
import os
import shutil
import re
import getopt
# Current directory
cwd = os.getcwd()
# infile = "ref_GMX.gro"
# outfile = "ref_GMX_reordered.gro"
infile = ""
outfile = ""
try:
opts, args = getopt.getopt(sys.argv[1:],"hi:o:",["infile=","outfile="])
except getopt.GetoptError:
print 'reorder.py -i <inputfile> -o <outputfile>'
sys.exit(2)
for opt, arg in opts:
if opt == '-h':
print 'reorder.py -i <inputfile> -o <outputfile>'
sys.exit()
elif opt in ("-i", "--infile"):
infile = arg
elif opt in ("-o", "--outfile"):
outfile = arg
print 'Input file is "', infile
print 'Output file is "', outfile
# Open input file
filerun=os.path.join(cwd, infile)
f = open(filerun,'r')
# Extract the number of atoms
header = f.readline()
natoms_raw = f.readline()
natoms = int(natoms_raw.strip())
# Extract all atom information
atoms = {}
residues = {}
resnum = []
resname = []
atname = []
satnum = []
atnum = []
coords = []
for iatom in range(0,natoms):
line = f.readline()
resnum.append(line[0:5])
resname.append(line[5:10])
atname.append(line[10:15])
satnum.append(line[15:20])
atnum.append(int(satnum[iatom].strip()))
coords.append(line[20:80])
# Extract cell information
box = f.readline()
# BEWARE : standard list copy is by pointer, not by value
# The newlist=list(oldlist) systax is needed to copy by value
newresnum=list(resnum)
newresname=list(resname)
newatname=list(atname)
newatnum=list(atnum)
newsatnum=list(satnum)
newcoords=list(coords)
# Indexes for Na+ and Cl- need to be corrected so that they are grouped together
def update_indexes():
newresnum[new_index]=resnum[new_index]
newresname[new_index]=resname[old_index]
newatname[new_index]=atname[old_index]
newatnum[new_index]=atnum[new_index]
newsatnum[new_index]=satnum[new_index]
newcoords[new_index]=coords[old_index]
na_indexes = []
cl_indexes = []
for iatom in range(0,natoms):
if resname[iatom] == " Na+":
na_indexes.append(iatom)
else :
if resname[iatom] == " Cl-":
cl_indexes.append(iatom)
nacl_indexes=na_indexes+cl_indexes
nacl_indexes.sort()
nacl_indexdiff=nacl_indexes[-1]-nacl_indexes[0]
if len(nacl_indexes)!= nacl_indexdiff+1:
print "Na+ and Cl- ions are scattered at different places, cannnot continue"
print "Na+/Cl- First index : ", nacl_indexes[0]
print "Na+/Cl- Last-first index difference :", nacl_indexdiff
print "Na+/Cl- list length", len(nacl_indexes)
quit()
if min(na_indexes)<min(cl_indexes):
first="Na"
else :
first="Cl"
if first == "Na":
firstna=min(na_indexes)
new_na_indexes=range(min(na_indexes),min(na_indexes)+len(na_indexes))
for i in range(0,len(na_indexes)):
old_index=na_indexes[i]
new_index=new_na_indexes[i]
update_indexes()
new_cl_indexes=range(min(na_indexes)+len(na_indexes),min(na_indexes)+len(na_indexes)+len(cl_indexes))
for i in range(0,len(cl_indexes)):
old_index=cl_indexes[i]
new_index=new_cl_indexes[i]
update_indexes()
else :
firstcl=min(cl_indexes)
new_cl_indexes=range(min(cl_indexes),min(cl_indexes)+len(cl_indexes))
for i in range(0,len(cl_indexes)):
old_index=cl_indexes[i]
new_index=new_cl_indexes[i]
update_indexes()
new_na_indexes=range(min(cl_indexes)+len(cl_indexes),min(cl_indexes)+len(cl_indexes)+len(na_indexes))
for i in range(0,len(na_indexes)):
old_index=na_indexes[i]
new_index=new_na_indexes[i]
update_indexes()
# Write vthe output file
fileout=os.path.join(cwd, outfile)
f = open(fileout,'w')
f.write(header,)
f.write(natoms_raw,)
for iatom in range(0,natoms):
f.write('%s%s%s%s%s' % (newresnum[iatom],newresname[iatom],newatname[iatom],newsatnum[iatom],newcoords[iatom] ),)
f.write(box)
f.close()
Reference
- Shahrokh, K., Orendt, A., Yost, G.S. and Cheatham, T.E., III (2012), Quantum mechanically derived AMBER-compatible heme parameters for various states of the cytochrome P450 catalytic cycle. J. Comput. Chem., 33: 119-133. https://doi.org/10.1002/jcc.21922
- Dai, M., Yu, X., Peng, X., Guo, Z., Yu, H., & Li, A. (2024). Modulating regioselectivity of CYP107J3-catalyzed isophorone hydroxylation by disrupting the hydrophobic balance of the substrate binding pocket. Green Synthesis and Catalysis.
致谢
感谢“凯算-生物医药交流”微信群各位大佬的帮助
感谢“凯算-生物医药交流”微信群大佬"华东理工大学王龙兴"大佬的倾力帮助
感谢湖北大学酶分子工程及合成生物技术研究团队郭志勇老师答疑解惑
祝以上同志身体健康，科研顺利~~~
扫一扫
843
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
