【PySCF学习1】-分子结构之定义分子

分子结构

定义一个分子
- 法一(使用Mole.build())
- 法二(先定义分子坐标信息及基组，再调用Mole.build())
- 法三(使用pyscf.M() 或者Mole.M())
分子几何结构位置信息
- 坐标单位(使用unit)
- 分子原子数量(使用natm)
- 使用内坐标(Z-矩阵)定义分子几何结构位置信息
- 元素表达方式
- 使用label标记某个原子
- 坐标信息输入方式
- 从文件中输入分子几何结构位置信息

一天，导师分享了一篇知乎的文章在群里，我点进去一看，感觉这文章应该就是分享给我看的。然后浏览了优秀答主的近期动态，刚好发现了基于python语言而开发的量化程序，这简直就是救命稻草。因为本人目前只有python语言用得上手一点，在学习量化计算的过程中，时常因为上手常用的量化程序非常慢而感到苦恼。因此看到pyscf后，我想如果我能从熟悉的语言出发，学习量化计算也许压力会小一点。因此基于官方教程，我复现一些我会用到的功能，并加上自己的理解，增强自己的学习记忆。所以接下来分享的内容大多数仅是单纯地重复官方教程，感兴趣的同学可以前往官网找适合自己的教程。

定义一个分子

以定义一个水分子为例

法一(使用Mole.build())

直接导入gto包，并引用gto中的Mole()函数，直接在使用mole.build()时定义分子坐标信息及基组信息

from pyscf import gto
mol = gto.Mole()
mol.build(atom = '''O 0 0 0; H  0 1 0; H 0 0 1''',basis = 'sto-3g')

这里的：O 0 0 0; H 0 1 0; H 0 0 1就是分子元素及笛卡尔坐标信息；sto-3g为所选用的基组
基组对元素定义了各自计算时要用的壳层类型及数目、收缩系数和高斯函数指数。

法二(先定义分子坐标信息及基组，再调用Mole.build())

from pyscf import gto
mol = gto.Mole()
mol.atom = '''O 0 0 0; H  0 1 0; H 0 0 1'''
mol.basis = 'sto-3g'
mol.build()

法三(使用pyscf.M() 或者Mole.M())

可以使用更为简短的函数pyscf.M() or Mole.M()对分子进行定义
使用pyscf.M()

import pyscf
mol = pyscf.M(atom = '''O 0 0 0; H  0 1 0; H 0 0 1''',basis = 'sto-3g')

使用Mole.M()


from pyscf import gto
mol = gto.M(atom = '''O 0 0 0; H  0 1 0; H 0 0 1''',basis = 'sto-3g')

分子几何结构位置信息

在前面定义分子的代码中可以发现分子元素及坐标信息是调用了关键字atom，上面都是使用了元素+笛卡尔坐标信息形式。下面详细介绍其他信息。

坐标单位(使用unit)

单位一般可选Angstrom (缩写AU)或者Bohr(缩写B)

from pyscf import gto
mol = gto.Mole()
mol.atom = '''O   0. 0. 0.H   0. 1. 0.H   0. 0. 1.
'''
mol.unit = 'B'#这里定义了单位为Bohr，如果单位选为Angstrom，可以改成AU

分子原子数量(使用natm)

检查所建立的分子中的原子数量

from pyscf import gto
mol = gto.Mole()
x=mol.build(atom = '''O 0 0 0; H  0 1 0; H 0 0 1''',basis = 'sto-3g')
print(x.natm)#输出结果是3

使用内坐标(Z-矩阵)定义分子几何结构位置信息

上面的例子在定义分子几何结构位置信息时均是使用了笛卡尔坐标，还有另一种方式是使用内坐标。
笛卡尔坐标中原子的位置是绝对的，相对于笛卡尔坐标轴；
内坐标是通过键长、键角及二面角信息来定义原子之间的相对位置。
Z-矩阵格式
元素符号原子1 键长原子2 键角原子3 二面角
(二面角可以使用右手法则来确定方向)

from pyscf import gto
mol = gto.Mole()
mol.atom = '''OH  1  1.2H  1  1.2  2 105
'''

元素表达方式

可以使用核电荷数、元素符号(大小写均支持)来表示某种元素

from pyscf import gto
mol = gto.Mole()
mol.atom = '''8 0. 0. 0.; h 0. 1. 0; H 0. 0. 1.'''

使用label标记某个原子

在研究中可能出现不同位置的原子需要采用不同基组来处理的情况，因此在定义分子几何结构位置信息时，可以标记某个原子，方便再后续为标记的原子指定规定的基组。
【标记方法】在原子符号前面或者后面加上数字1234567890 或者一些特殊符号，比如 ~!@#$%^&*()_+.?:<>[]{}| (不能用 ,与 ; )

from pyscf import gto
mol = gto.Mole()
mol.atom = '''8 0 0 0; h:1 0 1 0; H@2 0 0 1'''
mol.unit = 'B'#单位
mol.basis = {'O': 'sto-3g', 'H': 'cc-pvdz', 'H@2': '6-31G'}#对标记的原子分别指定使用不同的基组
mol.build()
print(mol._atom)

输出结果为：

[('O', [0.0, 0.0, 0.0]), ('H:1', [0.0, 1.0, 0.0]), ('H@2', [0.0, 0.0, 1.0])]

坐标信息输入方式

在定义分子几何结构位置信息时，我们可以灵活地结合循环、numpy等功能，利用信息存在的规律性，更省事地输入分子位置信息，这对于定义一些大分子时，对于节省工作量非常有效。
【原子位置信息输入格式】

atom = [[atom1, (x, y, z)],[atom2, (x, y, z)],...[atomN, (x, y, z)]]

例如：

mol.atom = [['O',(0, 0, 0)], ['H',(0, 1, 0)], ['H',(0, 0, 1)]]

可以换一种方式表示：

mol.atom = (('O',numpy.zeros(3)), ['H', 0, 1, 0], ['H',[0, 0, 1]])

从文件中输入分子几何结构位置信息

从文件中导入atom的信息，利用：

mol.atom = "my_molecule.xyz"

例如：

from pyscf import gto
mol = gto.Mole()
mol.atom = "atom.xyz"
mol.unit = 'B'#注意要定义单位
x=mol.build(basis = 'sto-3g')
print(x.natm)
print(mol._atom)

【输出信息为】

3
[('O', [0.0, 0.0, 0.0]), ('H', [0.0, 1.0, 0.0]), ('H', [0.0, 0.0, 1.0])]

其中文件atom.xyz内的信息为：

3O 0. 0. 0.
H 0. 1. 0.
H 0. 0. 1.

注意.xyz文件中的原子位置信息基本格式是：
原子数
空行
元素坐标