看到“schrodinger 薛定谔”这个关键词,很多人的第一反应是那只既死又活的猫。但我作为一个常年跑计算模拟的人,脑子里闪过的其实是另一个东西:Schrödinger Suite,业内直接叫“薛定谔”,一套做分子模拟、药物设计、材料科学的高通量计算平台。薛定谔方程是量子力学的基石,而 Schrödinger 这家公司把这套物理模型做成了工程化软件,让人可以直接在图形界面里做分子对接、虚拟筛选、自由能微扰计算,算得上是药物发现领域绕不开的“重型装备”。
这篇文章不聊量子力学科普,就围绕我在实际项目里用 Schrödinger 套件做分子对接和虚拟筛选的经验展开。我会从最初的软件选型逻辑讲起,一步一步把结构准备、Glide 对接参数、MM-GBSA 重打分、FEP+ 边界这些核心环节拆开,再把实操中踩过的坑、总结的排查方法一并放出来。适合正在学计算化学、生物信息学,或者刚拿到薛定谔 license 准备跑第一个真实课题的人参考。
1. 先说清楚:此薛定谔非彼薛定谔
1.1 从量子力学的薛定谔到计算平台
Schrödinger 这个名字源于物理学家 Erwin Schrödinger,他提出的薛定谔方程描述了微观粒子的波函数演化。现代分子模拟的底层逻辑,本质上还是解这个方程——不过不是解析解,而是通过各种近似方法(分子力学、半经验方法、从头算、DFT)在计算机上逼近真实体系的能量和运动状态。Schrödinger 这家公司成立于1990年,核心就是用物理模型去预测分子的性质和行为,把学术圈里散落的算法工具包装成一套工业级工作流。
很多人第一次接触“薛定谔”是在药企实习或者计算化学课程里,打开 Maestro 界面,看到一堆模块菜单。整套套件里最常见的几个模块是:Glide(分子对接)、Prime(蛋白结构预测与 MM-GBSA 重打分)、LigPrep(配体准备)、Desmond(分子动力学模拟)、FEP+(自由能微扰)、Canvas(化学信息学与构效关系分析)。它们不是独立的小工具,而是共享同一套力场参数、同一套文件格式、同一个图形界面,这意味着从受体准备到结果分析,整个流程都是串起来的,不太需要在不同软件之间玩“文件搬家”。
1.2 谁在什么场景下会用到它
我见过三类人会真正用到这套软件。第一类是药物研发团队的计算化学家,他们拿它做靶点的虚拟筛选、先导化合物优化,配合湿实验验证;第二类是结构生物学出身的研究者,手里有晶体结构或冷冻电镜结构,想快速看这个靶点能装下什么形状的分子;第三类是跟计算沾边的学生和跨界研究者,比如做农药设计、材料筛选、甚至酶催化机理的人,也会借用对接和动力学模拟来补充实验解释。
这些场景有一个共同点:核心需求都是“预测”。结晶一个蛋白可能花掉几个月,合成一个化合物又需要时间和经费,而计算可以通过打分函数先筛一遍,把候选范围从百万级压到几十个。所以 Schrödinger 解决的不是“算出真相”,而是“用较少资源找到最可能正确的方向”。这点非常重要,后面所有参数选择、流程设计,都是围绕它来的。
2. 用 Schrödinger 搭一条分子对接流水线
2.1 结构准备:Protein Prep Wizard 和 LigPrep,别跳过这一步
我第一次用薛定谔跑对接的时候,以为拿一个 PDB 结构直接 Glide 就能出结果,结果对出来的打分烂得离谱。后来才明白,PDB 文件里存的是实验解析的原子坐标,但这个坐标系里其实缺了很多化学上必须的信息:氢原子基本都没有,组氨酸的质子化状态不明确,Asp、Glu 等残基该带几个质子也没定,活性口袋里可能还混着结晶缓冲液分子、金属离子、二聚体界面上的水。
Protein Prep Wizard(蛋白准备向导)做的就是这些“补齐”工作。实际操作里我会依次做三件事。第一,给蛋白加氢并分配正确的质子化状态;第二,用 Prime 补缺失的残基和侧链,同时处理掉不合理的原子位置冲突;第三,做一轮能量最小化,让整个结构在力场下达到一个稳定构象。这里要特别提醒:最小化的原子位移限制不要设太大,默认 0.3 Å 就够了,否则会把活性位点压变形,后面的对接全废。
配体端的准备用 LigPrep。它会生成配体的 3D 构象,计算特定 pH 下的质子化状态,还能产生合理的立体异构体和环构象。这一步很多人不重视,直接拿 2D 结构的 smiles 去对接,结果发现 Glide 报错或者产出一堆高能构象。常见做法是设定 pH 7.0 ± 2.0,让 Epik 模块按生理条件计算所有可电离基团的质子化状态,再对每个输入结构生成最多 32 个低能 3D 构象。听起来简单,但这一步直接决定后续对接能不能找到“对的姿势”。
提示:结构准备是整个流程的地基。模板结构里如果带着共晶配体,Prep Wizard 默认会保留它,不要盲目删掉。这个配体往往是网格生成的定位点,也是判断对接结果是否合理的锚定参照。
2.2 Glide 分子对接参数怎么选
Glide 是 Schrödinger 的招牌对接工具,它把配体放进受体结合口袋,搜索所有可能的结合构象,然后用 GlideScore 打分。具体用哪种精度,取决于你要解决的问题。
Glide 提供了三档模式:HTVS(高通量虚拟筛选)、SP(标准精度)、XP(超高精度)。HTVS 最粗糙也最快,适合动辄百万级的超大库初筛;SP 是日常最常用的精度,兼顾速度和可靠性;XP 会对配体-受体之间的形状互补、疏水作用做更精细的惩罚和奖励,但极其耗时,适合对几十个命中做深度分析。
网格生成这一步最容易被忽略。Glide 需要先定义一个受体网格(Receptor Grid),网格的中心通常放在共晶配体的质心上,这样能保证结合口袋被完整覆盖。网格盒子的大小不需要贪大,默认的 15 Å × 15 Å × 15 Å 一般够用;盒子里还可以指定某些残基为柔性残基,让侧链在对接时微调,但每加一个柔性残基,计算量会指数增加,所以一般只对活性位点里最关键的几个残基开启。
如果做一个小规模的筛选,我的习惯是先用 SP 跑一遍,把打分前 20% 的配体拿出来再用 XP 重新对接,这两轮的结果综合起来决定下一步实验。很多人上来直接跑 XP,一个配体跑十几分钟,换来的是被噪声淹没的分数,纯属浪费时间。
2.3 打分函数怎么选、怎么看
对接完成后,你会在 Maestro 的结果表里看到一列 GlideScore。它表示的是配体结合到受体上时的一种“经验性评分”,不是真实的结合自由能,但数值越负,说明预测的结合越强。它由好几项加在一起:静电相互作用、范德华力、氢键、疏水接触、溶剂化效应、以及配体内部的应变能。
理解“打分函数不是能量”这一点很重要。真实结合自由能要考虑熵、去溶剂化、蛋白构象变化,这些在对接打分里要么是近似项,要么根本没有。所以分数差 0.3 分根本不算差距,我自己一般以 1 分以上作为筛选阈值。另一点是,分数高(更负)不代表活性好,因为打分函数可能被某些“表面友好”的分子误导,比如分子太大、埋在口袋里强行填满空隙导致范德华项虚高。
我见过不少人直接按 GlideScore 排序取前 50 个去做活性测试,结果命中率惨淡。我的做法是:先按分数看排名,再逐个目视检查结合姿势——看配体有没有完全离开口袋、有没有严重的原子碰撞、极性基团有没有和水或缺电子区域形成合理相互作用。视觉检查能过滤掉 30% 到 50% 的“分数好看但姿势离谱”的候选。
3. 进阶:从虚拟筛选到结合自由能计算
3.1 虚拟筛选流程设计与富集率
如果你面对的化合物库有几万甚至上百万个分子,那就不是一个个跑 Glide 的事了,而是一个系统流程设计问题。我习惯的虚拟筛选管线分四步。
第一步是预过滤。所有配体先过一遍物理化学性质过滤器,比如分子量、AlogP、氢键供体/受体数量、可旋转键数,参考类药五规则。再跑一遍 ADMET 预测,把有明显毒性风险或透膜性差的分子提前踢掉。这一步能把库缩到原来的三分之一甚至更少。
第二步是构象和异构体准备。用 LigPrep 批量生成 3D 结构。注意控制立体异构体数量,没有手性中心的分子不要生成无谓的异构体,否则后面对接时间直接翻倍。
第三步是分级对接。所有化合物先跑 HTVS,按分数保留前 20%~30%;剩下的跑 SP;SP 结果里取前 1000 到 2000 个跑 XP。每一级都没有必要把上一级的全部结果都拿过来,筛的就是“低分基本不值得细看”的思路。
第四步是重打分和聚类。用 Prime MM-GBSA 对 XP 命中的前 200 个分子重新打分,再对化学结构聚类,选每一类里分数最高的一两个。这样最终留下的几十个化合物在结构和理化性质上都有代表性和差异,避免 50 个命中全是同一个母核的不同尾巴。
富集率是评价筛选流程优劣的关键指标。操作上可以准备一组已知活性化合物和一组诱饵分子,混入候选库中跑一遍完整流程,看看活性化合物是否比诱饵排在更前面。如果活性化合物大部分排在库的头部,说明流程的富集能力靠谱;如果活性化合物排名跟随机差不多,那即使表面分数好看,也说明你的受体结构或准备步骤有问题。
3.2 Prime MM-GBSA 的批量重打分
GlideScore 是快速筛选的机枪,但到了决赛圈,你会需要精度更高一点的评价。Prime MM-GBSA 计算的是配体结合前后的能量差,综合了分子力学能量(MM)、连续溶剂化模型(GBSA)里的极性项和非极性项。它的核心优势在于部分考虑了溶剂化效应,这是对接打分里处理得比较粗的部分。
实际操作不复杂:在 Glide 的结果列表里选中一批配体,右键选择 Prime MM-GBSA,设置力场用 OPLS4、溶剂模型用 VSGB,其他保持默认,然后提交任务。它会自动对每个配体做一个小规模的构象采样,并给出一个 delta G bind 的估计值。通常可以把 Glide 的排名和 MM-GBSA 的排名做一个交叉验证,两者都靠前的分子优先做实验。
但这个手段也有明显的边界。它仍是一个“端态”计算,忽略了配体结合过程中的熵效应与通路上的中间态,因此它输出的数值适合做排序,不适合当真实自由能解读。比如对比 A、B 两个类似物,说“B 比 A 预估强 5 kcal/mol”是合理的;但说“A 的结合自由能就是 -48.6 kcal/mol”就没有物理意义,因为绝对数值受力场参数和模型简化影响太大。在我自己的流程里,MM-GBSA 永远是“相对比较”工具,绝不单独用它拍板一个化合物能不能进合成清单。
3.3 FEP+ 的思路和适用边界
如果项目推进到先导化合物优化阶段,常见的问题是:在某个核心骨架上,我想把这个苯环换成吡啶,把甲基换成乙基,或者在这里加一个氟原子,活性会变好还是变坏?这时 MM-GBSA 的精度已经不够了,就需要引入自由能微扰(FEP+)。
FEP+ 的核心思路是用一系列中间态连接两个结构高度相似的配体,计算从一个配体“演化”到另一个配体过程中的自由能差值。因为两个分子结构接近,很多误差会在求差的过程中相互抵消,所以它预测的相对结合自由能精度可以做得非常高,误差通常在 1 kcal/mol 以内,远好于 MM-GBSA,也更适合指导化学家往哪个方向修饰。
但 FEP+ 不是拿来随便跑的。它的第一个硬性要求是配体之间结构足够相似,最好是只在局部官能团上有差异;如果你拿一个完全不同的骨架去做 FEP+,中间态很难收敛,结果就失去了意义。第二个要求是计算资源充足,每个配体需要搭一张合理的图(哪些配体之间建微扰边),每个边都要跑分子动力学模拟,动辄需要 GPU 集群跑几天。第三个要求是你得有一个可靠的共晶结构或对接结构作为起点,起点蛋白质构象如果不对,后面一切都白搭。
我自己项目里的分工是这样的:初筛用 Glide,中筛用 MM-GBSA,最后的十来个类似物决定具体合成顺序时,才上 FEP+。这个三级递进既保证了效率,也把最贵的计算留给了最重要的决策点,不会出现“GPU 集群跑了一礼拜,结果发现模拟的体系根本不是活性构象”这种惨剧。
4. 实操中我踩过的坑
4.1 别拿原始 PDB 直接对接
我刚开始做激酶项目时,直接从 PDB 下载了一个分辨率不错的晶体结构,没有用 Prep Wizard 处理,只是删掉了水分子和其他链,就生成了网格开始对接。结果筛出来一个打分极高的化合物,合成出来却完全没活性。事后回头检查才发现,原始结构里 DFG 基序附近有一个关键残基的侧链密度不完整,PDB 里存的坐标是扭曲的,Prep Wizard 能自动检测并修正,而我跳过了这一步。
正确操作是:每次下载 PDB 后,先看序列信息和配体信息,再用 Protein Prep Wizard 按默认流程走一遍,重点观察 Log 窗口里有没有“missing residues”或“alternate conformations”提示。如果有缺失长度超过五六个残基的 loop,看它离活性位点远不远,远的话直接保持缺口,不要强行建模;近的话要考虑换一个更完整的结构,因为模型补出来的 loop 不确定性很大,对接结果可信度低。
4.2 质子化状态和异构体的问题
配体端最容易翻车的是异构体数量失控。有一次我想筛一个含有两个手性中心的母核,LigPrep 跑完提醒我生成了 4 个立体异构体各 32 个构象,总共 128 个文件,对接时间直接翻了 16 倍。后来我学会了在有明确实验数据的前提下,用 chiral 标记指定只有目标手性构型参与计算,或者用 LigPrep 的“retain specified stereoisomers”功能把不需要的异构体删掉。
另一件常被忽略的事是水分子的处理。有些活性位点里会有一个水分子,它同时在蛋白和配体之间搭桥,就好像两个陌生人之间的传话人,处于关键的位置。准备受体时如果无条件删除所有水分子,某些配体可能就少了一个重要的极性相互作用,打分被低估。我的做法是:保留所有太阳下看起来协调的水分子(与蛋白极性残基距离在 3.0 Å 左右,且不与非极性残基冲突的),生成网格时再把它们当作受体的一部分参与打分。
4.3 资源消耗与并行调优
Schrödinger 的任务调度逻辑跟普通程序不太一样,很多新手会误以为把所有 CPU 核心都填满就能加速。实际上,单个 Glide 对接任务本身对多核利用有限,提速的关键是“同时跑多个任务”。我习惯在 Linux 服务器上把配体库拆成几个子文件,每个子文件单独提交一个 Glide 任务,利用批处理脚本把它们并行调度起来。
Desmond 和 FEP+ 则刚好相反,它们强烈依赖 GPU 加速,没有 GPU 的话跑一个 100 ns 的分子动力学模拟可能要等一周,而用单张中端 GPU 能缩短到一天以内。所以如果你准备长期跑薛定谔的自由能计算,别省 GPU;如果只是跑对接和虚拟筛选,CPU 就够用,优先保证内存不要爆掉就行。
许可证连接也是一个隐蔽的坑。Schrödinger 的 license 是网络浮动式的,如果服务器网络不稳定或者 license 服务没启动,启动 Maestro 或命令行工具时会一直卡在“could not checkout license”的报错。排查时先看 license 服务器日志,再确认本机环境变量 SCHROD_LICENSE_FILE 是否指向正确的端口和地址,通常都能解决。
5. 常见问题速查
| 现象 | 可能原因 | 处理方法 |
|---|---|---|
| Prep Wizard 报错缺少大量残基 | PDB 文件本身分辨率差或截断 | 换用分辨率更高的结构,或删除无序区域,尽力补全关键残基 |
| Glide 网格生成失败 | 共晶配体坐标不完整或原子类型冲突 | 检查配体是否含金属离子或非标准残基编号,必要时重新用 LigPrep 处理 |
| LigPrep 产出的构象太多 | 立体异构体数量未限制 | 指定手性构型,控制每个结构的构象生成上限 |
| 对接分数很好但结合姿势奇怪 | 配体没进入口袋,或原子碰撞严重 | 用 Maestro 检查受力分布,关掉不合理的柔性残基设置 |
| MM-GBSA 结果不收敛 | 配体电荷设置不合理或结构构象过少 | 重新检查配体质子化状态,生成更多的初始构象 |
| FEP+ 自由能结果方差很大 | 中间态不够密,或者模拟时间不够 | 增加两个相邻配体间的中间态数量,延长平衡时长,检查是否选了结构差异过大的配体对 |
还有一个值得单独说的技巧:批量任务报错时,优先看 job 目录下的 write.log 文件。Schrödinger 的图形界面经常把错误吞掉,只显示一个红色的 “Job failed”,但具体的失败原因一定写在日志里。学会从日志的第一条 error 开始处理,比反复重跑要有效得多。实际项目里,我第一次排查一个批量对接任务失败的问题,就是因为某个配体包含未被力场定义的卤素原子类型,日志里一行 “cannot type atom” 瞬间让我定位到问题,换了参数版本后整个队列就顺了。
另外建议定期备份 Maestro 中构建好的自定义网格和结果视图。薛定谔的项目文件结构里除了 .mae 主文件,还有配套的 .maegz、.log、.grd 等文件。只拷贝主文件回家换个机器继续看结果时,经常发现结构缺失,就是因为没有把同一目录下的附属文件一起打包。把这些文件统一放进一个项目文件夹里,算是一个做了就会感谢自己的好习惯。
这套软件不便宜,学习曲线也不平缓,但只要理顺了“准备结构—生成网格—分级对接—重打分—自由能计算”这条主线,它就能成为决策效率的强杠杆。我个人这几年最深的体会是:薛定谔的价值不在于帮你跑出更高的分数,而在于逼着你想清楚每一步到底在算什么、误差在哪里、结论能支撑多大范围的判断。想明白这层,软件选哪家、参数怎么调,反而都是次要问题了。