【代码分享】利用Gaussian与自写Python脚本计算分子间电荷转移积分

在有机电子器件的研究中，分子间的电荷转移积分（Transfer Integral, JJ 或 VV）是衡量电荷传输效率的核心参数。今天我们分享一段Python代码，它通过读取Gaussian的计算结果（Fchk文件和Log文件），利用“投影法”（Projection Hamiltonian）来计算这一重要物理量。

1. 导入库与辅助工具函数

首先，我们需要导入处理数值运算的numpy库和数学函数sqrt。

#!bin/env python3import numpy as npfrom math import sqrt

2. 读取基组数量 (find_nbasis)

这是计算的第一步。我们需要知道计算中使用了多少个基函数，这将决定所有矩阵（如系数矩阵、重叠矩阵）的大小。该函数遍历.fchk文件，找到 “Number of basis functions” 关键字并返回其数值。

deffind_nbasis(fchk_dir):    txt = open(fchk_dir)for line in txt:if line.split()[0] == "Number"and line.split()[1] == "of"and line.split()[2] == "basis":returnint(line.split()[5])else:continuereturn0

3. 读取分子轨道能量 (find_moene)

此函数用于读取分子的轨道能级。.fchk文件中的数据通常是固定格式（每行5个数值）。函数通过循环读取这些浮点数，并将其转换为Numpy数组以便后续处理。

deffind_moene(fchk_dir, num):    txt = open(fchk_dir)    txt.seek(0)    line = txt.readline()while1:if line.split()[0] == "Alpha"and line.split()[1] == "Orbital"and line.split()[2] == "Energies":            line = next(txt)breakelse:            line = next(txt)continue    moene = []for i inrange(int(num / 5)):        moene.append(float(line.split()[0]))        moene.append(float(line.split()[1]))        moene.append(float(line.split()[2]))        moene.append(float(line.split()[3]))        moene.append(float(line.split()[4]))        line = next(txt)if num % 5 == 0:        txt.close()else:for i inrange(num % 5):            moene.append(float(line.split()[i]))        txt.close()return np.array(moene)

4. 读取分子轨道系数 (find_mocoeff)

分子轨道系数矩阵 CC 是连接基函数与分子轨道的桥梁。代码逻辑与读取能量类似，但由于数据量是基函数数量的平方（N×NN×N），读取后我们需要使用 .reshape 和 .T 将其还原为矩阵形式并进行转置，以符合标准定义。

deffind_mocoeff(fchk_dir, numb):    mocoeff = []    txt = open(fchk_dir)    txt.seek(0)    line = txt.readline()while1:if line.split()[0] == "Alpha"and line.split()[1] == "MO"and line.split()[2] == "coefficients":            line = next(txt)breakelse:            line = next(txt)continuefor i inrange(int(numb * numb / 5)):        mocoeff.append(float(line.split()[0]))        mocoeff.append(float(line.split()[1]))        mocoeff.append(float(line.split()[2]))        mocoeff.append(float(line.split()[3]))        mocoeff.append(float(line.split()[4]))        line = next(txt)if numb * numb % 5 == 0:        txt.close()else:for i inrange(numb * numb % 5):            mocoeff.append(float(line.split()[i]))        txt.close()    moc=np.array(mocoeff).reshape(numb, numb).Treturn moc

5. 读取重叠矩阵 (find_overlap)

这是一个比较复杂的函数。Gaussian 默认不在 .fchk 文件中输出重叠矩阵 SS，因此我们需要从 .log 文件中提取。

• 难点
  ：Log文件中矩阵是以“下三角”格式打印的，且每5个元素换行，空行处理较繁琐。
• 格式转换
  ：Gaussian使用 D (如 1.0D-05) 表示科学计数法，Python不识别，因此需要用 replace('D', 'E') 进行转换。
• 矩阵补全
  ：代码先填充下三角部分，最后通过双重循环利用对称性 overlap[i, j] = overlap[j, i] 补全整个矩阵。

deffind_overlap(log_dir, num):  # 寻找重叠矩阵    txt = open(log_dir)    txt.seek(0)    line = txt.readline().strip()for i inrange(10000000):if line.strip() == '':  # 如果此行为空，则换行            line = next(txt)elif line.split()[0] == "***"and line.split()[1] == "Overlap":  # 如果找到关键词，则换行，退出循环            line = next(txt)breakelse:            line = next(txt)  # 不匹配则换行，继续循环continue    overlap = np.zeros((num, num), float)  # 建立Nbasis*Nbasis大小的矩阵，用来盛放重叠矩阵# ---（省略部分中间繁琐的逐行读取赋值代码，逻辑是按Log文件格式解析下三角矩阵）---# ...# 处理最后剩余不足5个元素的行if k == 1:  # ...         overlap[num - 1, num - 1] = float(line.split()[1].replace('D', 'E'))elif k == 2:# ... (省略具体赋值逻辑，原理同上)pass# ... (省略 k==3 和 k==4 的逻辑)    txt.close()# 利用对称性补全矩阵for i inrange(num):for j inrange(num):            overlap[i, j] = overlap[j, i]return overlap

6. 主程序：数据读取与Fock矩阵构建

在主程序中，我们首先初始化文件路径。这里硬编码了文件路径，实际使用时可以取消注释 input 部分进行交互式输入。核心逻辑如下：

1. 读取二聚体的基组数、轨道能量、系数矩阵和重叠矩阵。
2. 构造Fock矩阵
   ：根据量子化学原理，F=S⋅C⋅ε⋅C−1F=S⋅C⋅ε⋅C−1。我们利用二聚体计算得到的系数矩阵 CC 和能量 εε，反推出二聚体在基组下的Fock算符矩阵 FF。

print('===========================================================')# ... (打印标题信息) ...# 这里定义输入文件路径 (使用时可用 input() 替代)dimerlog = "..."dimerfchk = "..."monomer1fchk = "..."monomer1log = "..."monomer2fchk = "..."monomer2log = "..."# 读取二聚体数据print("reading the number of dimer basis functions...")NBasis = find_nbasis(dimerfchk)print("reading the dimer molecular orbital energies...")MOene = find_moene(dimerfchk, NBasis)print("reading the dimer molelular orbital coefficients...")MOcoeff = find_mocoeff(dimerfchk, NBasis)print("reading the dimer overlap matrix...")Overlap = find_overlap(dimerlog, NBasis)# 计算二聚体 Fock 矩阵: F = S * C * E * C^(-1)print("calculating the dimer Fock matrix...\n")cfrag = np.zeros((NBasis, NBasis), dtype=np.double)tmp1 = np.dot(Overlap, MOcoeff)  # S * Ce = np.zeros((NBasis, NBasis), dtype=np.double)for i inrange(NBasis):    e[i, i] = MOene[i]          # 构造对角能量矩阵tmp = np.dot(tmp1, e)          # (S * C) * EMOcoeff_inv = np.linalg.inv(MOcoeff)Fock = np.dot(tmp, MOcoeff_inv) # 得到 Fock 矩阵

7. 组装片段轨道与计算转移积分

计算电荷转移积分通常采用“片段轨道”方法。

1. 片段轨道组装
   ：我们将单体1和单体2的分子轨道系数（MO系数）填充到与二聚体基组相同大小的矩阵 cfrag 中。这相当于在二聚体的基组下描述单体的轨道。
2. 对角化计算
   ：利用公式 Jeff=J12−0.5(e1+e2)S121−S122Jeff=1−S122J12−0.5(e1+e2)S12，其中 J12=C1TFC2J12=C1TFC2。
• e1,e2e1,e2
  : 单体轨道在二聚体势场下的有效轨道能。
• J12J12
  : 原始耦合矩阵元。
• S12S12
  : 轨道重叠积分。
• 最终输出考虑了重叠修正后的有效转移积分 JeJe，单位换算为 meV 和 kcal/mol。

# 读取单体1数据并填入 cfrag 左上角NBasis1 = find_nbasis(monomer1fchk)MOene1 = find_moene(monomer1fchk, NBasis1)MOcoeff1 = find_mocoeff(monomer1fchk, NBasis1)for i inrange(NBasis1):for j inrange(NBasis1):        cfrag[i, j] = MOcoeff1[i, j]# 读取单体2数据并填入 cfrag 右下角NBasis2 = find_nbasis(monomer2fchk)MOene2 = find_moene(monomer2fchk, NBasis2)MOcoeff2 = find_mocoeff(monomer2fchk, NBasis2)for i inrange(NBasis2):for j inrange(NBasis2):        cfrag[i + NBasis1, j + NBasis1] = MOcoeff2[i, j]# 设置需要计算的轨道范围 (HOMO/LUMO索引)ihomo1, ilumo1 = 53, 55ihomo2, ilumo2 = 165, 167print("calculating the charge transfer integral...")# 循环计算每对轨道之间的耦合for i inrange(ihomo1 - 1, ilumo1):for j inrange(ihomo2 + NBasis1 - 1, ilumo2 + NBasis1):# 计算 e1: C1^T * F * C1        ctmp = np.dot(Fock, cfrag[:, i])        e1 = np.dot(cfrag[:, i], ctmp) * 27.21138505# Hartry转eV# 计算 e2: C2^T * F * C2        ctmp = np.dot(Fock, cfrag[:, j])        e2 = np.dot(cfrag[:, j], ctmp) * 27.21138505# 计算 J12: C1^T * F * C2        ctmp = np.dot(Fock, cfrag[:, j])        J12 = np.dot(cfrag[:, i], ctmp) * 27.21138505# 计算 S12: C1^T * S * C2        ctmp = np.dot(Overlap, cfrag[:, j])        S12 = np.dot(cfrag[:, i], ctmp)# 计算有效轨道能量和有效转移积分 (Generalized 2-level model)        ee1 = 0.5 * ((e1 + e2) - 2.0 * J12 * S12 + (e1 - e2) * sqrt(1.0 - S12 ** 2)) / (1 - S12 ** 2)        ee2 = 0.5 * ((e1 + e2) - 2.0 * J12 * S12 - (e1 - e2) * sqrt(1.0 - S12 ** 2)) / (1 - S12 ** 2)        Je12 = (J12 - 0.5 * (e1 + e2) * S12) / (1 - S12 ** 2)# 打印结果print('{:5d}{:5d}{:10.3f}{:10.3f}{:10.3f}{:10.3f}{:10.3f}{:10.3f}{:10.3f}{:10.3f}'.format(i + 1, j - NBasis1 + 1,                                                                                                e1, e2, J12 * 1000, S12,                                                                                                ee1, ee2, Je12 * 1000,                                                                                                Je12 * 23.05))

总结：这段代码通过Python实现了从Gaussian原始数据到物理量计算的完整流程。虽然Log文件的解析部分稍显繁琐，但核心思想——即利用二聚体Fock矩阵和单体轨道系数的投影——是计算电子耦合最准确的方法之一。希望这段代码对你的科研工作有所帮助！

(注：代码中的 find_overlap 部分处理了剩余逻辑，实际运行时请确保Log文件路径正确且包含完整的Overlap矩阵输出)