python

# 导入pywfn所在文件路径
import sys;sys.path.append('d:/code/pywfn')

原子性质

所有的原子性质的计算器都在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.reader import LogReader
from pywfn.atomprop import charge

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

array([-0.1284, -0.1284, -0.1284, -0.1284, -0.1284, -0.1284,  0.1285,
        0.1285,  0.1285,  0.1285,  0.1285,  0.1285])

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.reader import LogReader
from pywfn.atomprop import charge

path='mols/C6H6.out'
reader=LogReader(path)
mol=Mole(reader)
caler=charge.Calculator(mol)
caler.lowdin()

array([-0.1564, -0.1598, -0.1598, -0.1564, -0.1598, -0.1598,  0.1571,
        0.1595,  0.1595,  0.1571,  0.1595,  0.1595])

hirshfeld电荷

计算公式

示例代码

python

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

path='mols/C6H6.out'
reader=LogReader(path)
mol=Mole(reader)
caler=charge.Calculator(mol)
caler.hirshfeld()

array([-0.0426, -0.0426, -0.0426, -0.0426, -0.0426, -0.0426,  0.0426,
        0.0427,  0.0427,  0.0426,  0.0427,  0.0427])

方向电子数

根据pocv方法，将原子的p轨道投影到不同方向得到投影系数矩阵，再以此带入原子电荷计算公式（Mulliken和Lowdin）得到的结果即为方向电子数，这里的方向需要自己指定。

示例代码pywfn.atomprop.Calculator.dirElects包含三个参数

atms:list[int] # 需要进行投影操作的原子索引列表（从1开始）
dirs:np.ndarray # 对应原子需要投影到的方向，形状为[n,3]，n==len(atms)
ctype:str='mulliken' # 可以为 mulliken 或 lowdin

python

from pywfn.base import Mole
from pywfn.reader import LogReader
from pywfn.atomprop import charge
import numpy as np

path='mols/C6H6.out'
reader=LogReader(path)
mol=Mole(reader)
caler=charge.Calculator(mol)
dirs=np.random.rand(4,3)
dirs=dirs/np.linalg.norm(dirs,axis=1)[:,None] # 将方向向量归一化
atms=[1,2,3,4] # 
caler.dirElects(atms,dirs,'mulliken') # 我们对四个原子的轨道投影到随机的方向

array([0.654 , 0.6438, 0.6531, 0.7069, 0.    , 0.    , 0.    , 0.    ,
       0.    , 0.    , 0.    , 0.    ])

π 电子数 (轨道投影法)

指定每个原子的方向为其所在局部平面的法向量方向，没有法向量的则忽略（通常为饱和C原子），然后计算方向电子数，得到的结果即为 π 电子数

示例代码

python

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

path='mols/C6H6_N.out'
reader=LogReader(path)
mol=Mole(reader)
caler=charge.Calculator(mol)
caler.piElects('mulliken') # 我们对四个原子的轨道投影到随机的方向

array([1.0202, 1.2101, 1.2101, 1.0202, 1.2101, 1.2101, 0.    , 0.    ,
       0.    , 0.    , 0.    , 0.    ])

π 电子数 (轨道分解法)

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

示例代码

python

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

path='mols/C6H6.out'
reader=LogReader(path)
mol=Mole(reader)
caler=charge.Calculator(mol)
caler.piDecomElects()

array([1.0011, 1.0013, 1.0013, 1.0011, 1.0013, 1.0013, 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.reader import LogReader
from pywfn.atomprop import activity

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

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

array([[-0.0426, -0.106 ,  0.1237,  0.1663,  0.0634,  0.1148, -0.1029],
       [-0.0426, -0.1734,  0.0436,  0.0862,  0.1309,  0.1085,  0.0447],
       [-0.0426, -0.1734,  0.0436,  0.0862,  0.1309,  0.1085,  0.0447],
       [-0.0426, -0.106 ,  0.1237,  0.1663,  0.0634,  0.1148, -0.1029],
       [-0.0426, -0.1734,  0.0436,  0.0862,  0.1309,  0.1085,  0.0447],
       [-0.0426, -0.1734,  0.0436,  0.0862,  0.1309,  0.1085,  0.0447],
       [ 0.0426, -0.0039,  0.104 ,  0.0614,  0.0465,  0.0539, -0.0149],
       [ 0.0427, -0.0217,  0.0925,  0.0498,  0.0643,  0.0571,  0.0145],
       [ 0.0427, -0.0217,  0.0925,  0.0498,  0.0643,  0.0571,  0.0145],
       [ 0.0426, -0.0039,  0.104 ,  0.0614,  0.0465,  0.0539, -0.0149],
       [ 0.0427, -0.0217,  0.0925,  0.0498,  0.0643,  0.0571,  0.0145],
       [ 0.0427, -0.0217,  0.0925,  0.0498,  0.0643,  0.0571,  0.0145]])

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

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

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

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

示例代码

python

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

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

caler=activity.Calculator(mol0)
caler.piFukui(molN,molP,'mulliken')

array([[ 0.9851,  1.0202,  0.706 ,  0.2792,  0.0351,  0.1571,  0.244 ],
       [ 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.706 ,  0.2792,  0.0351,  0.1571,  0.244 ],
       [ 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.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ],
       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ,  0.    ]])

输出的结果有7列，分别为：

N 电子分子π电子数
N+1 电子分子π电子数
N-1 电子分子π电子数
f-
f+
f0
df

化合价

原子的键级之和，使用mayer键级。化合价越大，原子的剩余成键能力越小，原子的反应能力就越小。

示例代码

python

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

mol=Mole(LogReader('./mols/C6H6.out'))
caler=activity.Calculator(mol)
caler.valence()

array([3.9412, 3.9412, 3.9412, 3.9412, 3.9412, 3.9412, 0.9334, 0.9334,
       0.9334, 0.9334, 0.9334, 0.9334])

自由价（反应矢量）

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

示例代码

python

import sys;sys.path.append('d:/code/pywfn')

from pywfn.base import Mole
from pywfn.reader import LogReader
from pywfn.atomprop import activity
import numpy as np

mol=Mole(LogReader('./mols/C6H6.out'))
caler=activity.Calculator(mol)
atm=1
dir_=np.array([0.0,0.0,1.0])
caler.vector(atm,dir_)

1.147283692768253

原子方向

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

法向量

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

示例代码

python

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

mol=Mole(LogReader('./mols/C6H6.out'))
caler=direction.Calculator(mol)
caler.normal(1)

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

最大波函数方向

通过迭代搜索的方法获取原子波函数最大值的方向，需要指定分子轨道索引及原子轨道符号（可以使用正则表达式）

示例代码

获取苯环第一个原子10号分子轨道中p原子轨道最大值的方向

python

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

mol=Mole(LogReader('./mols/C6H6.out'))
caler=direction.Calculator(mol)
caler.maxWeave(1,10,'P[XYZ]')

array([ 0.9964, -0.0842,  0.    ])

化学反应方向

pocv算法文章中定义的原子的可能的一些反应方向 [n,3]

示例代码

python

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

mol=Mole(LogReader('./mols/C6H6.out'))
caler=direction.Calculator(mol)
caler.reactions(1)

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

基座标

分子中每个原子的局部基座标。以字典的形式返回，并不是每个原子都有局部坐标系，主要用于轨道分解法

z轴是法向量，x、y轴是随机的。

示例代码

python

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

mol=Mole(LogReader('./mols/C6H6.out'))
caler=direction.Calculator(mol)
caler.bases()

{1: array([[-0.6763, -0.7367, -0.    ],
        [ 0.7367, -0.6763, -0.    ],
        [ 0.    , -0.    ,  1.    ]]),
 2: array([[-0.3156, -0.9489,  0.    ],
        [ 0.9489, -0.3156,  0.    ],
        [ 0.    ,  0.    ,  1.    ]]),
 3: array([[-0.9735, -0.2286,  0.    ],
        [ 0.2286, -0.9735,  0.    ],
        [ 0.    ,  0.    ,  1.    ]]),
 4: array([[-0.7928, -0.6095,  0.    ],
        [ 0.6095, -0.7928, -0.    ],
        [ 0.    ,  0.    ,  1.    ]]),
 5: array([[-0.4287, -0.9035, -0.    ],
        [ 0.9035, -0.4287,  0.    ],
        [-0.    ,  0.    ,  1.    ]]),
 6: array([[-0.7199, -0.6941,  0.    ],
        [ 0.6941, -0.7199,  0.    ],
        [ 0.    ,  0.    ,  1.    ]])}

原子能量

计算分子轨道能量在原子上的分布

原子电子能量

其实就是将每个分子轨道的能量分配到所有的原子的基函数

示例代码

python

from pywfn.base import Mole
from pywfn.reader import LogReader
from pywfn.atomprop import energy

mol=Mole(LogReader('./mols/C6H6.out'))
caler=energy.Calculator(mol)
caler.atmEngs()

array([-22.3995, -22.3993, -22.3993, -22.3995, -22.3993, -22.3993,
        -0.4095,  -0.4095,  -0.4095,  -0.4095,  -0.4095,  -0.4095])

原子π电子能量能量

根据分子轨道每个原子的π电子的分布，计算π电子的能量

示例代码

python

from pywfn.base import Mole
from pywfn.reader import LogReader
from pywfn.atomprop import energy

mol=Mole(LogReader('./mols/C6H6.out'))
caler=energy.Calculator(mol)
caler.atmPiEngs()

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

原子性质 ​

原子电荷 ​

Mulliken电荷 ​

lowdin电荷 ​

hirshfeld电荷 ​

方向电子数 ​

π 电子数 (轨道投影法) ​

π 电子数 (轨道分解法) ​

原子活性 ​

福井函数&双描述符 ​

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

化合价 ​

自由价（反应矢量） ​

原子方向 ​

法向量 ​

最大波函数方向 ​

化学反应方向 ​

基座标 ​

原子能量 ​

原子电子能量 ​

原子π电子能量能量 ​

原子性质

原子电荷

Mulliken电荷

lowdin电荷

hirshfeld电荷

方向电子数

π 电子数 (轨道投影法)

π 电子数 (轨道分解法)

原子活性

福井函数&双描述符

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

化合价

自由价（反应矢量）

原子方向

法向量

最大波函数方向

化学反应方向

基座标

原子能量

原子电子能量

原子π电子能量能量