小分子动力学模拟：溶剂效应与构象采样的计算策略

小分子（分子量<500 Da）的MD模拟和生物大分子的MD相比，物理本质更纯粹但也更挑剔。小分子没有蛋白质折叠那样的层次架构来约束自己的构象空间，缺少这种”内部约束”意味着构象采样收敛成了一个更棘手的问题。

溶剂模型的选择直接影响构象分布。显式溶剂（TIP3P、SPC/E等水模型）给出的溶剂化自由能面和隐式溶剂（GBSA、PBSA）的差异在某些分子上可以大到2~3 kcal/mol。做过一个含多个可旋转键的柔性药物分子，显式溶剂模拟发现它在水中存在三个不同构象簇，而用隐式溶剂模拟时这三个簇合并为一个过于光滑的自由能面——精细结构被抹平了。这组对比直观地说明：如果你关心的是构象之间的精确布居比，显式溶剂是绕不开的。

模拟时长不是信手拈来随便设的。小分子的构象转化时间尺度在纳秒级别——单键旋转的能垒一般为2~6 kcal/mol，对应的弛豫时间在室温下是几十皮秒到几纳秒。但多个可旋转键的协同构象转变需要更长的时间窗。经验是：可旋转键数≤5时，100~200ns的显式溶剂模拟通常能给出收敛的构象分布；可旋转键数≥10时，相同时间窗的构象采样可能严重不足——你会看到模拟反复访问某几个局部极小态而从未越过某个较高的旋转势垒。这种情况下需要用增强采样方法——REMD或MetaD——在不同温度窗口之间交换构象来加速势垒穿越。

REMD的参数设置也值得展开说一下。对于含50个重原子以下的小分子，10~16个温度窗口通常够用，最高温度建议设在400~500K之间（低于溶剂的沸点，否则盒子压力会失控）。相邻窗口之间的交换接受率控制在20%~40%之间是比较理想的——太低说明窗口间隔太大，太高则意味着窗口数过多浪费计算资源。这个接受率可以直接从模拟开始的几纳秒中提取，然后调整窗口分布。

力场对于小分子的参数化质量可以通过与QM计算结果做基准对比来验证。具体做法是：对小分子的关键可旋转键做一维或二维的QM（ωB97X-D/6-311++G(d,p)级别）势能面扫描，把QM得到的旋转势垒与MD力场下同一势垒做对比。如果偏差大于1 kcal/mol，就说明力场对这个分子的构象能量有定性级别的偏差——用它跑出的构象布居比没有物理意义。这种情况下要么换用专门针对该分子类型做了参数化的力场（如CGenFF或GAFF2），要么对问题参数做手动调整并在论文中详细说明。
自由能计算是小分子MD模拟的一个自然延伸。从构象分布到构象间的自由能差，中间需要经过玻尔兹曼反演——ΔG=-RT ln(P_A/P_B)，其中P_A和P_B是两个构象态的布居概率。这里有一个容易被忽略的数值陷阱：当某个构象态的布居概率很低时（比如不到1%），模拟中统计到的该态出现次数只有个位数，由此计算出的自由能差的统计误差可以大到2~3 kcal/mol——对于药物设计中要比较的亚μM级别的结合自由能来说，这个误差是完全不可接受的。解决方案是：如果关心的是两个低布居态之间的自由能差，不要从总轨迹中直接数布居——用伞形采样沿着一个低维度的反应坐标来做约束采样，WHAM重建自由能面。

https://www.keyanxueshu.com/category/md/