python

from importlib.metadata import version
print(version('pywfn'))
from datetime import datetime
print(datetime.now())

1.0.18
2026-07-21 14:28:53.489048

原子性质

所有的原子性质的计算器都在pywfn.atomprop子包下，其中的每一个模块封装了一种类型的原子性质的计算器

其包含的模块有

charge 计算原子的各种电荷及自旋
activity 计算原子的活性指标
direction 计算以原子为中心的不同类型的方向向量。如：原子轨道方向、法向量方向、可能的反应方向等。
energy 将分子轨道的能量分布到每个原子上。

每个模块下都有一个Calculator类，实例化时传入要计算的分子即可

原子电荷

说是电荷，其实直接计算得到的是每个原子上的电子数，用原子的核电荷数减去电子数即可得到原子电荷。

该计算器实例有个form属性，用以控制打印的结果时电子数量的形式还是电荷的形式，可以为number或charge

Mulliken电荷

计算公式

q_{A} = Z_{A} - \sum_{μ \in A} \sum_{ν} P_{μ ν} S_{μ ν}

示例代码

python

from pywfn.base import Mole
from pywfn.atomprop import charge

path='mols/C6H6.out'
mol=Mole.from_file(path)
caler=charge.Calculator(mol)
caler.mulliken()

array([6.12839599, 6.1284028 , 6.1284028 , 6.12839599, 6.1284028 ,
       6.1284028 , 0.87152149, 0.87153433, 0.87153433, 0.87152149,
       0.87153433, 0.87153433])

lowdin电荷

计算公式

q_{I} = Z_{I} - \sum_{μ \in I} (S^{1 / 2} \cdot P \cdot S^{1 / 2})

示例代码

python

from pywfn.base import Mole
from pywfn.atomprop import charge

path='mols/C6H6.out'
mole=Mole.from_file(path)
caler=charge.Calculator(mole)
caler.lowdin()

array([6.15640809, 6.15976901, 6.15976901, 6.15640809, 6.15976901,
       6.15976901, 0.84287667, 0.84048448, 0.84048448, 0.84287667,
       0.84048448, 0.84048448])

hirshfeld电荷

计算公式

示例代码

python

from pywfn.base import Mole
from pywfn.atomprop import charge

path='mols/C6H6.out'
mole=Mole.from_file(path)
caler=charge.Calculator(mole)
caler.hirshfeld()

array([6.04257693, 6.04258752, 6.04258752, 6.04257693, 6.04258752,
       6.04258752, 0.95733417, 0.95734849, 0.95734849, 0.95733417,
       0.95734849, 0.95734849])

pocv (分子轨道系数投影法)

根据pocv方法，将原子的p轨道投影到不同方向得到投影系数矩阵，得到新的分子轨道系数矩阵，带入原子电荷计算公式

示例代码

将苯环的每个C原子的投影方向指定为分子平面的法向量，同时忽略其它原子和其它符号的轨道系数，计算的结果就是π电子布居

python

from pywfn.base import Mole
from pywfn.atomprop import charge

path='mols/C6H6.out'
mole=Mole.from_file(path)
caler=charge.Calculator(mole)
dirs={
    0:[0.,0.,1.],
    1:[0.,0.,1.],
    2:[0.,0.,1.],
    3:[0.,0.,1.],
    4:[0.,0.,1.],
    5:[0.,0.,1.],
}
caler.pocv(dirs,False,False,"mulliken") # 我们对四个原子的轨道投影到随机的方向

[0.9851325082883295,
 0.9851422336033631,
 0.9851422336033628,
 0.9851325082883295,
 0.9851422336033631,
 0.9851422336033632,
 0.0,
 0.0,
 0.0,
 0.0,
 0.0,
 0.0]

π 电子数 (pocv)

使用pocv方法计算 π 电子数，返回原子投影方向(法向量)和每个原子的π电子布居

示例代码

同样计算苯环的π电子布居

python

from pywfn.base import Mole
from pywfn.atomprop import charge

path='mols/C6H6.out'
mole=Mole.from_file(path)
caler=charge.Calculator(mole)
dirs,eles=caler.pi_pocv('mulliken')
print(dirs) # 原子的投影方向(法向量)
print(eles) # 原子的π电子布居

{1: [0.0, 0.0, 1.0], 2: [0.0, 0.0, 1.0], 0: [-0.0, -0.0, 1.0], 3: [0.0, 0.0, 1.0], 4: [0.0, 0.0, 1.0], 5: [-0.0, -0.0, 1.0]}
[0.9851325082883295, 0.9851422336033631, 0.9851422336033628, 0.9851325082883295, 0.9851422336033631, 0.9851422336033632, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

mocv (分子轨道系数映射法)

根据mocv方法，为原子定义局部坐标，得到新的分子轨道系数矩阵，带入原子电荷计算公式

python

from pywfn.base import Mole
from pywfn.atomprop import charge

π 电子数 (mocv)

使用轨道分解法得到的 π 电子数，包含d轨道

示例代码

python

from pywfn.base import Mole
from pywfn.atomprop import charge

path='mols/C6H6.out'
mole=Mole.from_file(path)
caler=charge.Calculator(mole)
stms,atos,eles=caler.pi_mocv()
print(stms) # 每个原子的局部坐标系
print(atos) # 局部坐标系下保留的原子轨道
print(eles) # 每个原子的π电子布居

{4: T  ex: (  0.0000  0.0000  1.0000)  |  ey: ( -0.5000  0.8660  0.0000)  |  ez: (  0.8660  0.5000 -0.0000), 0: T  ex: ( -0.0000 -0.0000  1.0000)  |  ey: (  1.0000  0.0000  0.0000)  |  ez: (  0.0000 -1.0000  0.0000), 1: T  ex: (  0.0000  0.0000  1.0000)  |  ey: (  0.5000 -0.8660  0.0000)  |  ez: ( -0.8660 -0.5000  0.0000), 5: T  ex: ( -0.0000 -0.0000  1.0000)  |  ey: (  0.5000  0.8660  0.0000)  |  ez: (  0.8660 -0.5000  0.0000), 2: T  ex: (  0.0000  0.0000  1.0000)  |  ey: ( -0.5000 -0.8660  0.0000)  |  ez: ( -0.8660  0.5000  0.0000), 3: T  ex: (  0.0000  0.0000  1.0000)  |  ey: ( -1.0000  0.0000  0.0000)  |  ez: (  0.0000  1.0000 -0.0000)}
[2, 6, 12, 13, 17, 21, 27, 28, 32, 36, 42, 43, 47, 51, 57, 58, 62, 66, 72, 73, 77, 81, 87, 88]
[0.9999862455850268, 0.9999971240976323, 0.9999971240976323, 0.9999862455850267, 0.9999971240976322, 0.9999971240976322, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

原子活性

原子活性就是分子中每个原子发生化学反应能力。

福井函数&双描述符

计算公式

福井函数

f_{A}^{+} = q_{A}^{N} - q_{A}^{N + 1}

f_{A}^{-} = q_{A}^{N - 1} - q_{A}^{N}

f_{A}^{0} = (q_{A}^{N - 1} - q_{A}^{N + 1}) / 2

双描述符

Δ f_{A} = f_{A}^{+} - f_{A}^{-}

其中 $Q_{A}^{N + 1}$ 是具有N+1个电子的分子的原子电荷， $Q_{A}^{N}$ 是具有N个电子的分子的原子电荷， $Q_{A}^{N - 1}$ 是具有N-1个电子的分子的原子电荷。

$f^{+}$ 被亲核试剂进攻的能力
$f^{-}$ 被亲电试剂进攻的能力
$f^{0}$ 发生自由基反应的能力
$Δ f$ 越正表示该位点越容易受到亲核试剂进攻

其实就是计算三个分子的原子电荷之差，因此可以使用不同的电荷计算方式，可以使用的有：Mullien、Lowdin和Hirshfeld电荷

示例代码

python

from pywfn.base import Mole
from pywfn.atomprop import activity

mol0=Mole.from_file('./mols/C6H6.out')
molN=Mole.from_file('./mols/C6H6_N.out')
molP=Mole.from_file('./mols/C6H6_P.out')

caler=activity.Calculator(mol0)
caler.fukui(molN,molP,'hirshfeld')

array([[ 6.04257693,  6.10596772,  5.87627143,  0.16630551,  0.06339079,
         0.11484815, -0.10291472],
       [ 6.04258752,  6.17344958,  5.95639472,  0.0861928 ,  0.13086206,
         0.10852743,  0.04466926],
       [ 6.04258752,  6.17344958,  5.95639472,  0.0861928 ,  0.13086206,
         0.10852743,  0.04466926],
       [ 6.04257693,  6.10596772,  5.87627143,  0.16630551,  0.06339079,
         0.11484815, -0.10291472],
       [ 6.04258752,  6.17344958,  5.95639472,  0.0861928 ,  0.13086206,
         0.10852743,  0.04466926],
       [ 6.04258752,  6.17344958,  5.95639472,  0.0861928 ,  0.13086206,
         0.10852743,  0.04466926],
       [ 0.95733417,  1.00385036,  0.8959606 ,  0.06137357,  0.04651619,
         0.05394488, -0.01485737],
       [ 0.95734849,  1.02165415,  0.90749866,  0.04984983,  0.06430566,
         0.05707774,  0.01445583],
       [ 0.95734849,  1.02165415,  0.90749866,  0.04984983,  0.06430566,
         0.05707774,  0.01445583],
       [ 0.95733417,  1.00385036,  0.8959606 ,  0.06137357,  0.04651619,
         0.05394488, -0.01485737],
       [ 0.95734849,  1.02165415,  0.90749866,  0.04984983,  0.06430566,
         0.05707774,  0.01445583],
       [ 0.95734849,  1.02165415,  0.90749866,  0.04984983,  0.06430566,
         0.05707774,  0.01445583]])

该段代码实例化计算器时也是传入一个N个电子的分子，调用函数的时候传入N+1个电子的分子、N个电子的分子及计算电荷的类型（这里使用Mulliken电荷）。

输出了7列数据，分别对应着 $q_{(N)}$ ， $q_{(N + 1)}$ ， $q_{(N - 1)}$ ， $f^{-}$ ， $f^{+}$ ， $f^{0}$ ， $Δ f$

π电子福井函数与双描述符

使用 π 电子分布计算得到的福井函数与双描述符

示例代码

分别使用pocv和mocv方法计算 $π$ 福井函数

python

from pywfn.base import Mole
from pywfn.atomprop import activity

mol0=Mole.from_file('./mols/C6H6.out')
molN=Mole.from_file('./mols/C6H6_N.out')
molP=Mole.from_file('./mols/C6H6_P.out')

caler=activity.Calculator(mol0)
dirs,vals=caler.fukui_pi_pcov(molN,molP,'mulliken')
# 打印每个原子的法向量
for atm,dir in dirs.items():
    dx,dy,dz=dir
    print(f"{atm+1:>2}|{dx:>10.4f}{dy:>10.4f}{dz:>10.4f}")
# 打印计算数据
print(f"{'e0':>10}{'en':>10}{'ep':>10}{'fp':>10}{'fn':>10}{'f0':>10}{'df':>10}")
for e0,en,ep,fp,fn,f0,df in vals:
    print(f"{e0:>10.4f}{en:>10.4f}{ep:>10.4f}{fp:>10.4f}{fn:>10.4f}{f0:>10.4f}{df:>10.4f}")

 2|    0.0000    0.0000    1.0000
 4|    0.0000    0.0000    1.0000
 1|   -0.0000   -0.0000    1.0000
 3|    0.0000    0.0000    1.0000
 5|    0.0000    0.0000    1.0000
 6|   -0.0000   -0.0000    1.0000
        e0        en        ep        fp        fn        f0        df
    0.9851    1.0202    0.7060    0.2792    0.0351    0.1571   -0.2440
    0.9851    1.2101    0.8769    0.1082    0.2249    0.1666    0.1167
    0.9851    1.2101    0.8769    0.1082    0.2249    0.1666    0.1167
    0.9851    1.0202    0.7060    0.2792    0.0351    0.1571   -0.2440
    0.9851    1.2101    0.8769    0.1082    0.2249    0.1666    0.1167
    0.9851    1.2101    0.8769    0.1082    0.2249    0.1666    0.1167
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

python

from pywfn.base import Mole
from pywfn.atomprop import activity

mol0=Mole.from_file('./mols/C6H6.out')
molN=Mole.from_file('./mols/C6H6_N.out')
molP=Mole.from_file('./mols/C6H6_P.out')

caler=activity.Calculator(mol0)
stms,atos,vals=caler.fukui_pi_mocv(molN,molP,'mulliken')
for atm,stm in stms.items():
    print(f"{atm:>2} {stm}")
print(atos)
print(f"{'e0':>10}{'en':>10}{'ep':>10}{'fp':>10}{'fn':>10}{'f0':>10}{'df':>10}")
for e0,en,ep,fp,fn,f0,df in vals:
    print(f"{e0:>10.4f}{en:>10.4f}{ep:>10.4f}{fp:>10.4f}{fn:>10.4f}{f0:>10.4f}{df:>10.4f}")

 0 T  ex: ( -0.0000 -0.0000  1.0000)  |  ey: (  1.0000  0.0000  0.0000)  |  ez: (  0.0000 -1.0000  0.0000)
 4 T  ex: (  0.0000  0.0000  1.0000)  |  ey: ( -0.5000  0.8660  0.0000)  |  ez: (  0.8660  0.5000 -0.0000)
 2 T  ex: (  0.0000  0.0000  1.0000)  |  ey: ( -0.5000 -0.8660  0.0000)  |  ez: ( -0.8660  0.5000  0.0000)
 5 T  ex: ( -0.0000 -0.0000  1.0000)  |  ey: (  0.5000  0.8660  0.0000)  |  ez: (  0.8660 -0.5000  0.0000)
 1 T  ex: (  0.0000  0.0000  1.0000)  |  ey: (  0.5000 -0.8660  0.0000)  |  ez: ( -0.8660 -0.5000  0.0000)
 3 T  ex: (  0.0000  0.0000  1.0000)  |  ey: ( -1.0000  0.0000  0.0000)  |  ez: (  0.0000  1.0000 -0.0000)
[2, 6, 12, 13, 17, 21, 27, 28, 32, 36, 42, 43, 47, 51, 57, 58, 62, 66, 72, 73, 77, 81, 87, 88]
        e0        en        ep        fp        fn        f0        df
    1.0000    1.0400    0.7200    0.2800    0.0400    0.1600   -0.2400
    1.0000    1.2300    0.8900    0.1100    0.2300    0.1700    0.1200
    1.0000    1.2300    0.8900    0.1100    0.2300    0.1700    0.1200
    1.0000    1.0400    0.7200    0.2800    0.0400    0.1600   -0.2400
    1.0000    1.2300    0.8900    0.1100    0.2300    0.1700    0.1200
    1.0000    1.2300    0.8900    0.1100    0.2300    0.1700    0.1200
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000
    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000

自由价

原子的最大成键数减去原子的键级之和，使用mayer键级。自由价越大，原子的剩余成键能力越大，原子的反应能力就越大

示例代码

python

from pywfn.base import Mole
from pywfn.atomprop import activity

mol=Mole.from_file('./mols/C6H6.out')
caler=activity.Calculator(mol)
caler.freev()

[[4.0, 3.8370293597429224, 0.16297064025707764],
 [4.0, 3.837064737173867, 0.1629352628261329],
 [4.0, 3.837064737173867, 0.1629352628261329],
 [4.0, 3.8370293597429224, 0.16297064025707764],
 [4.0, 3.8370647371738675, 0.16293526282613247],
 [4.0, 3.837064737173867, 0.1629352628261329],
 [4.0, 0.930465601835863, 3.069534398164137],
 [4.0, 0.930484989485762, 3.069515010514238],
 [4.0, 0.930484989485762, 3.069515010514238],
 [4.0, 0.9304656018358632, 3.0695343981641368],
 [4.0, 0.9304849894857621, 3.069515010514238],
 [4.0, 0.9304849894857621, 3.069515010514238]]

返回两列数据，第一列为原子的最大成键数，第二列为原子的自由价

方向自由价（活性矢量）

基于pocv方法及原子指定的方向，对原子的p轨道进行投影，然后根据得到的分子轨道计算相关键级之和，随后用每种原子的标准值减去键级之和即可得到自由价（反应矢量）。

示例代码

python

from pywfn.base import Mole
from pywfn.atomprop import activity

mol=Mole.from_file('./mols/C6H6.out')
caler=activity.Calculator(mol)
caler.freev_dir(0,[0.0,0.0,1.0],[1,5])

(3.0, [0.7984412064219861, 0.7984412064219861], 1.4031175871560277)

原子方向

计算得到与原子相关的各种方向

法向量

计算指定原子所在局部平面的法向量

示例代码

计算苯环第一个原子的法向量

python

from pywfn.base import Mole
from pywfn.atomprop import direction

mole=Mole.from_file('./mols/C6H6.out')
caler=direction.Calculator(mole)
caler.normal_vector(1)

[0.0, 0.0, 1.0]

反应方向

计算原子左右可能发生反应的一系列方向，主要的原则是与原子的键之间的夹角要大于等于90°

示例代码

计算苯环的第一个原子的反应方向，结果是垂直于平面的上下两个方向

python

from pywfn.base import Mole
from pywfn.atomprop import direction

mol=Mole.from_file('./mols/C6H6.out')
caler=direction.Calculator(mol)
caler.reactions(0)

array([[-0., -0.,  1.],
       [ 0.,  0., -1.]])

局部坐标系

计算原子的局部坐标系

x轴是法向量方向，x、y轴是随机的。对于线性分子，z轴是键轴方向

示例代码

计算苯环的第一个原子的局部坐标系

python

from pywfn.base import Mole
from pywfn.atomprop import direction

mole=Mole.from_file('./mols/C6H6.out')
caler=direction.Calculator(mole)
caler.LCS(0)

T  ex: ( -0.0000 -0.0000  1.0000)  |  ey: (  1.0000  0.0000  0.0000)  |  ez: (  0.0000 -1.0000  0.0000)

原子性质 ​

原子电荷 ​

Mulliken电荷 ​

lowdin电荷 ​

hirshfeld电荷 ​

pocv (分子轨道系数投影法) ​

π 电子数 (pocv) ​

mocv (分子轨道系数映射法) ​

π 电子数 (mocv) ​

原子活性 ​

福井函数&双描述符 ​

π电子福井函数与双描述符 ​

自由价 ​

方向自由价（活性矢量） ​

原子方向 ​

法向量 ​

反应方向 ​

局部坐标系 ​

原子性质

原子电荷

Mulliken电荷

lowdin电荷

hirshfeld电荷

pocv (分子轨道系数投影法)

π 电子数 (pocv)

mocv (分子轨道系数映射法)

π 电子数 (mocv)

原子活性

福井函数&双描述符

π电子福井函数与双描述符

自由价

方向自由价（活性矢量）

原子方向

法向量

反应方向

局部坐标系