对生物大分子进行第一性原理计算一般采用什么code？

在生物大分子领域，利用第一性原理（Ab initio）方法进行计算，特别是涉及到电子结构和量子力学层面的模拟，是一个复杂但信息丰富的过程。这类计算的目的是从最基础的物理原理出发，不依赖于经验参数或实验拟合，来描述分子的性质，例如能量、力、电荷分布、光谱特征等。对于像蛋白质、核酸、多糖这样的大型分子，其复杂性对计算的精度和效率提出了极高的要求。

在生物大分子第一性原理计算中，最核心的代码通常是基于量子化学（Quantum Chemistry）或密度泛函理论（Density Functional Theory, DFT）的计算包。这些计算包实现了求解薛定谔方程（或其等价形式）的各种近似方法，并将这些方法应用于具体的分子体系。

下面我将详细介绍一些常用的计算代码，并尽可能深入地解释它们的原理和在生物大分子研究中的应用：

1. 密度泛函理论 (DFT) 计算包

DFT 是目前在计算化学和材料科学领域最广泛使用的方法之一，因为它在精度和计算成本之间取得了较好的平衡，尤其适合处理包含数百甚至数千个原子的生物大分子。

核心原理：

DFT 的核心思想是，系统的基态能量完全由其电子密度唯一确定。这意味着我们不需要直接求解多电子波函数，而是转而求解一系列单电子的（ KohnSham）方程。 KohnSham 方程的形式与 HartreeFock 方程相似，但其关键在于一个交换关联泛函 (ExchangeCorrelation Functional, XC Functional) 的存在。这个泛函包含了所有难处理的多体量子效应，其精确形式是未知的，因此 DFT 的精度很大程度上取决于所使用的 XC 泛函。

常用的 XC 泛函类别：

LDA (Local Density Approximation): 最早也是最简单的泛函，假设交换关联能量只取决于局部的电子密度。对于电子密度变化缓慢的体系较好，但对于电子密度变化剧烈的生物分子，精度有限。
GGA (Generalized Gradient Approximation): 在 LDA 的基础上引入了电子密度梯度，考虑了电子密度变化率的影响。例如：
PBE (PerdewBurkeErnzerhof): 一种广泛使用的、从头推导的 GGA 泛函，对许多体系表现良好。
BLYP (BeckeLeeYangParr): 结合了 Becke 的交换泛函和 LeeYangParr 的关联泛函，在描述分子键和官能团方面有不错的表现。
MetaGGA: 进一步引入了电子动能密度，提供了更丰富的电子密度信息，理论上能获得更高的精度。
Hybrid Functionals: 将一部分 HartreeFock 的精确交换能混合到 DFT 交换能中。这通常能显著提高描述电子激发、键解离能等性质的精度，尤其对于范德华力（如生物分子间的相互作用）的描述有一定改进。例如：
B3LYP: 最流行的杂化泛函之一，在许多化学反应和分子性质预测中表现出色。
PBE0: PBE 的杂化版本。
M06 系列 (M06L, M06, M062X): 由 Truhlar 团队开发的一系列泛函，针对不同类型的化学键和相互作用进行了参数化，在描述非共价相互作用方面表现突出，这对生物大分子研究非常重要。
Double Hybrid Functionals: 进一步包含了 MP2 (MøllerPlesset perturbation theory) 风格的关联能，精度更高，但计算成本也显著增加。

数值实现与近似：

为了解决 KohnSham 方程，通常需要将电子波函数展开成基函数 (Basis Functions) 的线性组合。基函数又可以分为：

高斯型基函数 (Gaussian Type Orbitals, GTOs):
Slater Type Orbitals (STOs): 更接近氢原子轨道，理论上更精确，但积分计算困难。
Cartesian Gaussians: 容易计算积分，但会引入“伪键”和线性相关问题，需要进行正交化。
Spherical Gaussians: 避免了伪键问题，是目前主流高斯基组实现方式。
基组大小 (Basis Set Size): 基组的大小直接影响计算的精度和成本。常用的基组系列包括：
STOnG: 基于 STOnG 的近似（如 STO3G, STO6G）。
Pople 系列 (e.g., 631G, 631G, 6311++G): 包含不同数量的基函数（高斯函数）来描述原子轨道，星号 () 表示极化函数（允许原子轨道变形），加号 (+) 表示弥散函数（扩展轨道以更好地描述弱束缚电子）。
Dunning 系列 (e.g., ccpVDZ, ccpVTZ, augccpVnZ): 称作完全一致的基组，它们是通过外推到无穷大基组极限的系统性方法构建的，在能量和性质预测上通常提供更高的精度。`aug` 前缀表示加入了弥散函数。
平面波基函数 (Plane Wave Basis Sets, PWs):
原理: 将电子波函数展开成平面波的傅里叶级数。
优点: 具有优良的收敛性，且与周期性边界条件（如晶体或模拟盒中的分子）非常兼容。
缺点: 对于局域化的电子（如生物分子中的孤立分子），需要非常大的平面波截止能量（Energy Cutoff）才能达到与高斯基组相当的精度，计算成本会很高。
适用性: 在模拟周期性体系（如蛋白质在溶液中的周期性盒子模拟）时，PWs 是更优的选择。
原子中心基函数 (AtomicCentered Basis Functions, ACBFs) / 有限元方法 (Finite Element Method, FEM):
原理: 使用 piecewise polynomial 样式的函数来表示波函数。
优点: 能够自适应地在电子密度变化剧烈的区域（如化学键附近）使用更精细的网格，而在电子密度平坦的区域使用粗糙的网格，从而实现计算效率和精度的良好平衡。
缺点: 实现和优化相对复杂。

关键计算步骤：

1. 几何构型构建: 首先需要一个合理的初始原子坐标。这可以通过实验数据（如 X 射线晶体衍射、NMR）获得，或者通过分子力学（Molecular Mechanics, MM）方法进行粗略优化。
2. 电子密度和波函数迭代求解:
初始化: 猜测一个初始的电子密度和 KohnSham 轨道。
构造 Fock 矩阵: 根据当前的电子密度和选定的 XC 泛函，计算 Fock 矩阵（包含动能、核吸引、电子电子排斥、交换和关联贡献）。
对角化 Fock 矩阵: 求解 Fock 矩阵的本征值和本征向量，得到新的 KohnSham 轨道和能量。
更新电子密度: 根据新的轨道计算新的电子密度。
收敛性检查: 将新计算的电子密度与前一步的电子密度进行比较。如果差异小于预设的阈值，则认为计算收敛。否则，重复上述步骤。
3. 性质计算: 一旦计算收敛，就可以计算各种性质，如总能量、原子力的梯度（用于几何优化）、电荷分布、偶极矩、振动频率（用于红外、拉曼光谱）等。

常用的 DFT 计算包 (Code):

Gaussian:
特点: 历史悠久，功能强大，支持多种 DFT 泛函、基组，以及各种化学反应（如过渡态搜索）、光谱计算、分子动力学（QM/MM）、含时 DFT (TDDFT) 等。是学术界最常用的计算化学软件之一。
生物大分子应用: 尽管 Gaussian 本身并非为大规模生物分子专门设计，但通过 QM/MM 方法，可以将活性中心（如酶的催化位点）用 DFT 计算，而周围的大分子环境用 MM 方法描述，这大大降低了计算成本，使得研究酶促反应机理、药物靶点相互作用等成为可能。
ORCA:
特点: 免费且功能强大，特别适合计算大型分子。提供了多种 DFT 泛函、STL (Scaled Total Energy) 方法（用于处理范德华力），以及用于光谱学、稀土元素计算等的高级功能。其在并行计算和内存管理方面做了很多优化，对生物大分子研究友好。
生物大分子应用: ORCA 在处理具有多个活性位点、金属配位、大共轭体系的生物分子方面表现突出。
QChem:
特点: 专注于高精度量子化学计算，支持多种前沿的 DFT 方法、杂化泛函，以及高效的 TDDFT 算法，特别适合计算激发态性质，如光吸收、荧光。
生物大分子应用: 用于研究光合作用中的光捕获蛋白、荧光蛋白的发色团性质、以及生物传感器等。
NWChem:
特点: 一个开源的、高性能的量子化学计算软件，支持从头算（ab initio）和 DFT 计算，以及分子动力学模拟。能够处理大型体系，并且在并行计算方面有很好的扩展性。
生物大分子应用: 适合用于研究蛋白质折叠、DNA 结构、以及大型复合物的电子结构。
VASP (Vienna Ab initio Simulation Package):
特点: 主要面向固体材料和表面科学，但其强大的平面波基组和 PAW (Projector Augmented Wave) 赝势方法，也使其能够模拟周期性盒子中的生物分子，例如蛋白质在晶体中的状态，或者蛋白质吸附在表面上的情况。
生物大分子应用: 蛋白质晶体学（电子密度拟合）、蛋白质表面相互作用、以及某些在固相或膜环境中具有特殊性质的生物分子。
CP2K:
特点: 一个开源的、高性能的分子模拟软件，结合了 DFT 和 MM 方法，特别适合模拟大量粒子（如溶剂分子）的体系。其特殊的基组方法（如 GTH 赝势和基组）以及稀疏矩阵技术，使其在模拟大体系时具有较高的效率。
生物大分子应用: 蛋白质在溶液中的动力学行为、溶剂效应、水合物壳结构等。
TERACHEM:
特点: 专注于高效的 DFT 计算，特别是在处理大量原子和长程相互作用方面有优势。其采用了多种技术来降低计算成本，如快速多极方法 (FMM) 和基于网格的近似。
生物大分子应用: 模拟长距离的非共价相互作用，如长链脂肪酸与蛋白质的结合，或者离子在生物通道中的传输。

2. 局域轨道/多中心展开方法 (Localized Orbitals / Multicenter Expansion Methods)

对于非常大的生物分子，全电子 DFT 计算的成本会变得极其高昂（通常与分子大小的四次方到五次方成正比）。为了克服这个问题，研究人员开发了许多近似方法，其中一种重要思路是利用电子的局域性。

核心原理：

生物大分子在空间上是庞大的，但其化学反应和相互作用往往发生在局部的区域，例如酶的活性位点、分子识别的结合位点等。这意味着，在计算这些局域性质时，远距离的电子对局部区域的影响相对较小，或者可以通过更粗糙的近似来处理。

常用方法/代码：

QM/MM (Quantum Mechanics/Molecular Mechanics):
原理: 将体系划分为两个区域：量子力学 (QM) 区域，包含对化学性质至关重要的原子（如反应中心），用 DFT 或其他量子化学方法描述；分子力学 (MM) 区域，包含其余的原子，用更快速但精度较低的分子力学力场描述。两个区域之间的相互作用则需要特别处理（如基于电荷的相互作用，可能需要 QM 区域的电荷分布来驱动 MM 区域）。
代码: 许多主流的量子化学和分子动力学软件都支持 QM/MM 接口，例如：
Amber: 广泛用于分子动力学模拟，可以与 Gaussian, ORCA, GAMESS 等量子化学包结合。
GROMACS: 另一款流行的分子动力学软件，同样支持 QM/MM。
CHARMM: 也是一个经典的分子模拟软件。
ONETEP (OneTime Electronic Structure Calculation):
原理: 采用延展的、非重叠的局域轨道 (Extended, Nonoverlapping Local Orbitals, NONO) 作为基函数，并结合快速多极方法 (Fast Multipole Method, FMM) 来高效计算长程的库仑相互作用。它能够实现伪O(N) 的计算复杂度，其中 N 是原子数。
特点: 专为处理大规模体系设计，对生物大分子和周期性体系都非常适用。
生物大分子应用: 模拟大型蛋白质、DNA 链，甚至病毒颗粒的电子结构。
Linear Scaling DFT (LSDFT) 方法:
原理: 通过对 Fock 矩阵进行稀疏化 (Sparsification) 或者利用局域性的核（Kernels），将计算复杂度降低到 O(N)。通常需要将波函数展开在局域化的基函数（如小波、Bsplines）上，并利用“电池小电池”或“区域分解”等技术。
代码: 许多研究组开发了自己的 LSDFT 代码，或者将 LSDFT 模块集成到现有软件包中。一些知名的 LSDFT 研究代码包括：
DC (Divide and Conquer): 将体系划分为多个小单元，每个单元内部进行 DFT 计算，并利用缓冲区域来处理单元间的相互作用。
FBP (FragmentBased Prescription): 基于分子碎片化和相互作用的分析。
Embedding Methods (嵌入法):
原理: 类似于 QM/MM，但更侧重于将一个 QM 区域嵌入到一个更大、更简单的环境（如点电荷、连续介质模型、或者低精度的 QM 模型）中。
代码: COSMO (Conductorlike Screening Model) 是一种常用的连续介质模型，常用于模拟溶剂效应，许多 DFT 代码都支持。

3. 从头算 (Ab Initio) 方法（HartreeFock, MP2, CCSD(T) 等）

虽然 DFT 是目前处理生物大分子最常用的第一性原理方法，但一些更精确的从头算方法（如 MP2, CCSD(T)）在描述特定性质时仍然不可或缺，尤其是在需要精确捕捉电子相关性时。

核心原理：

这些方法直接求解多电子薛定谔方程，并且不依赖于 XC 泛函。它们的精度通常随着方法级别的提高而增加，但计算成本也呈指数级增长。

HartreeFock (HF): 最基础的从头算方法，将多电子波函数近似为单电子轨道的斯莱特行列式，忽略了电子相关性。
MP2 (MøllerPlesset perturbation theory of second order): 在 HF 的基础上引入了二阶微扰理论，考虑了电子相关的部分效应，能较好地描述范德华力。
CCSD(T) (Coupled Cluster Singles Doubles with perturbative Triples): 被认为是“金标准”方法，能非常精确地描述电子相关性，但计算成本极高。

适用性与局限性：

由于计算成本原因，将 HF、MP2、CCSD(T) 等方法直接应用于完整的生物大分子（如一个完整的蛋白质）通常是不可行的。它们通常仅用于：

小分子模型: 研究生物大分子中的关键基团或反应中心的简化模型。
QM/MM 的 QM 部分: 在 QM/MM 方法中，如果活性中心非常小，且需要非常高的精度，可以选择这些方法。
基准测试: 用于验证 DFT 泛函或 MM 力场的准确性。

常用代码：

mesmos (MultiPlatform Quantum Chemistry Package):
特点: 一个强大的、开源的量子化学计算包，支持各种从头算方法 (HF, MP2, CCSD, CCSD(T)) 和 DFT 方法。
生物大分子应用: 主要用于对小分子模型或 QM/MM 方法中的 QM 区域进行高精度计算。
DALTON:
特点: 专注于计算光谱性质（如 IR, Raman, NMR, UVVis）和各种响应性质。支持多种高精度从头算方法。
生物大分子应用: 用于预测小分子模型的光谱特征，帮助理解生物分子在光谱实验中的信号。
PSI4:
特点: 开源、高性能的量子化学软件，支持广泛的从头算和 DFT 方法，尤其在激发态计算和强关联体系方面有特色。
生物大分子应用: 类似 DALTON，用于高精度计算小分子模型。

4. 关键考虑因素与挑战

计算资源: 生物大分子第一性原理计算需要大量的计算时间和内存。高性能计算（HPC）集群、GPU 加速是必需的。
精度与成本的权衡: 选择合适的理论方法（DFT 泛函、基组）和近似方法（QM/MM, LSDFT）是在可接受的计算时间内获得足够精度信息的核心问题。
体系的复杂性: 生物大分子常常是柔性的、动态的，并处于复杂的溶剂和离子环境中。这些因素都会影响计算结果的解释。
处理周期性边界条件 vs. 孤立分子: 蛋白质在晶体或膜中是周期性的，而研究其在溶液中的行为则常采用周期性盒子模拟。选择适合模拟体系的计算方法（如平面波 vs. 高斯基组）很重要。
泛函的选择: 对于生物大分子，尤其是涉及非共价相互作用（如氢键、范德华力、ππ 堆积）的体系，泛函的选择至关重要。研究表明，某些杂化泛函或专门针对非共价相互作用优化的泛函（如 M062X, ωB97XD）通常表现更好。
溶剂效应: 溶剂（主要是水）在生物过程中起着至关重要的作用。如何准确地将溶剂效应包含在内是一个持续的研究课题，常用的方法包括连续介质模型（PCM, COSMO）和显式溶剂分子模拟（QM/MM 的一部分）。

总结一下，在生物大分子领域，第一性原理计算的主要代码集中在基于密度泛函理论 (DFT) 的计算包。其中，Gaussian、ORCA、QChem、NWChem 是在学术界最常用的代表。

Gaussian 因其功能全面和用户友好性而受到青睐，尤其适合中小型体系或 QM/MM 计算。
ORCA 则以其对大型体系的优化和强大的功能在处理复杂生物分子时表现出色。
QChem 在激发态性质计算方面具有优势，适合研究光化学过程。
NWChem 是一个强大的开源选项，尤其在并行计算方面有不错表现。

对于非常大的体系，ONETEP 等采用局域轨道和 FMM 的线性标度 DFT 方法，以及 CP2K 这种结合了 DFT 和 MM 的方法，是更具竞争力的选择。同时，QM/MM 仍然是连接量子精度和模拟大规模体系的关键桥梁，其接口支持众多计算包。

选择哪种代码，最终取决于研究的具体问题、对精度和计算成本的要求，以及可用的计算资源。通常，一个项目的成功也离不开对理论方法的深入理解和对计算结果的仔细验证。

对生物大分子进行第一性原理计算通常有两个思路：

1. 暴力硬算，代表者是Stanford大学，Todd Martinez开发的TeraChem。Stanford是目前世界上已知各专业的学生都能写出良好代码的大学（估计这也是为什么没听说Stanford大学毕业生劝退什么的，个人读Stanford几个课题组的工作，感到在Stanford，CS已经成了大学数学一样的基础知识了）。TeraChem是目前已知最好用的可GPU加速的量子化学软件。其他的比如VASP什么的据说也能支持GPU，但是都不那么好使。我自己博士期间参与过编写分块计算蛋白质电子结构的程序，比较了之后发现，TeraChem居然比分块计算还快。。。。。。但是由于DFT、HF之类对体系的计算量成O(N^4)增长，硬算是有上限的。我3年前测试的体系大约在200个残基以下，是小型蛋白质。我用MFCC分块，每一个CPU只计算一个残基，包含30个左右的原子；居然没有一个节点，两个GPU的TeraChem速度快。不过我不知道目前TeraChem的计算上限在哪。顺带说一下，TeraChem是按节点卖的。他们定义一个标准节点，里边只能包含两块GPU，5万元人民币一个节点。如果是一个节点装了4块GPU，价格就要增加了。
   刚才看到，基于平面波的PWmat也蹦出来票了一把。
   第一性原理计算2019-nCoV病毒分子3CL水解酶结构
   不过，个人对这东西用在生物大分子量化计算中的价值持保留态度。等以后开会的时候咨询好了。PWmat算材料好像还是很有帮助的。知乎上也有他们的团队成员

分块计算。这又有许多流派。比较流行的有日本人开发的FMO、我熟悉的MFCC等。去年正好出了一篇综述：John M. Herbert，J. Chem. Phys. 151, 170901 (2019); https://doi.org/10.1063/1.5126216

1.5 目前DeepMD发展迅猛，以后估计可以做机器学习加速的第一性原理分子动力学来模拟蛋白质了。分子力场有可能成为历史，或者只占小半壁江山了。 @深势科技 DP Technology 应该给我加鸡腿。

1. 上边有 @翟羽 说得好，只用分子力场做MD足够了。这是因为，对蛋白质进行理论计算，主要关心的是药物和蛋白质的结合能，这是一个自由能计算问题。且实验误差在1 kcal/mol数量级，完全不需要第一性原理计算。用只基于数据的化学信息学、生物信息学就可以做药物开发，结果又快又好。关于复杂体系的分子动力学模拟，我课题组近期做了一系列工作，开发了一套用于分子晶体模拟的极化（polarized，不是polarizable）力场，目前最好的结果是成功模拟了生物铁电体的相变，发表在一区杂志J. Phys. Chem. Letters上，欢迎引用：J. Phys. Chem. Letters 10, 1319 (2019).
2. 但是针对某特定的药物和特定的蛋白质相结合，详细研究机理，就要上量子化学。不过通常的做法也不需要把整个体系都用量子力学计算，一般只需要ab initio QM/MM MD，针对药物作用位点进行高的精度计算即可。而更进一步研究蛋白质催化化学反应机理、生物体系中的量子效应如光合作用、视觉机理、生物发光机理的话，量子动力学也可能要用上，如MCTDH。甚至，如果研究生物大分子中的量子统计问题，对于经典环境、量子系统，就需要上Ehrenfest MD、surface hopping等，如果环境也是量子的，就需要HEOM、QUAPI等等。