如何使用python中的量子化学软件包（包括psi4和pyscf） - How to use quantum chemistry software packages in Python (including psi4 and pyscf)

发表于 - Posted on 编辑于 - Edited on 字数 - Word count: 1.3k 阅读时间 - Reading time ≈ 5 mins.

虽说现有的量子化学软件如gaussian，orca，和gamess的功能非常强大，但有时为了方便，或者为了加深理解，也可以用python来做简单的量子化学计算。

While current quantum chemistry software like Gaussian, Orca, and GAMESS are incredibly robust, there are times when using Python for simple quantum chemistry calculations is more convenient or helps deepen our understanding.


python中的量子化学软件包有以下几个：

The following quantum chemistry software packages are available in Python:

psi4

pyscf

PyQuante（因为年代过于久远而安装失败 failed to install due to its age）


此外，还有几个类似的量子计算包，但不能做单分子量子化学计算，在本文中暂不考虑：

Additionally, there are several similar quantum computing packages that can't perform single-molecule quantum chemistry calculations, which are not considered in this article:

QuTiP（用于量子物理 used in quantum physics）

Quantum ESPRESSO（用于固体物理 used in solid-state physics）


本文将介绍如何使用psi4和pyscf。出于方便，本文将使用google colab安装并使用所有的软件包。

This article focuses on how to utilize psi4 and pyscf. For ease of use, all software packages will be installed and employed using Google Colab.


安装psi4的最简单的方法是使用conda，在google colab中使用conda的方法如下所示。

The simplest method to install Psi4 in Google Colab is using Conda, as demonstrated below:

1
2
3
4
5
6
!wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
!chmod +x Miniconda3-latest-Linux-x86_64.sh
!./Miniconda3-latest-Linux-x86_64.sh -b -f -p /usr/local

import sys
sys.path.append('/usr/local/lib/python3.10/site-packages/')


下载完conda之后，使用conda安装psi4：

After Conda is downloaded, psi4 can be installed using Conda.

1
!conda install psi4 python=3.10 -c conda-forge/label/libint_dev -c conda-forge -y


这里测试一下psi4是否安装正常：使用b3lyp/6-31G进行乙醇分子的单点能计算。乙醇分子的结构由rdkit产生（也是一个在python上可以安装并使用的软件包，rdkit使用代码如下所示）：

To test if psi4 is properly installed, we'll perform a single-point energy calculation on an ethanol molecule using B3LYP/6-31G. The molecular structure of ethanol is generated by RDKit, another Python package that can be installed and used.

1
2
3
4
5
6
7
8
9
10
11
12
13
!pip install rdkit

from rdkit import Chem
from rdkit.Chem import AllChem

ethanol_smiles = "CCO"
ethanol_molecule = Chem.MolFromSmiles(ethanol_smiles)

ethanol_molecule = Chem.AddHs(ethanol_molecule)
AllChem.EmbedMolecule(ethanol_molecule, AllChem.ETKDG())
AllChem.MMFFOptimizeMolecule(ethanol_molecule)

print(Chem.MolToXYZBlock(ethanol_molecule))


使用psi4计算乙醇的单点能的代码如下：

The code for calculating the single-point energy of ethanol using psi4 is as follows:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import psi4

molecule = psi4.geometry("""
C     -0.887089    0.175064   -0.012535
C      0.460489   -0.515516   -0.046535
O      1.442965    0.307267    0.565572
H     -0.847478    1.127768   -0.550817
H     -1.658782   -0.453327   -0.465842
H     -1.176944    0.403676    1.018306
H      0.768712   -0.724328   -1.075460
H      0.419486   -1.462073    0.500177
H      1.478640    1.141468    0.067135
""")

psi4.set_options({'basis': '6-31G'})

energy = psi4.energy('B3LYP', molecule=molecule)
print(f"B3LYP Energy: {energy} Hartree")


结果为：

The result is:

B3LYP Energy: -154.9882218123809 Hartree

其中量子化学计算所需时间大约为5s。

The quantum chemistry calculation takes about 5 seconds.

pyscf可以使用pip安装：

Pyscf can be simply installed by pip:

1
!pip install --prefer-binary pyscf


使用pyscf计算乙醇的单点能的代码如下：

The code for calculating the single-point energy of ethanol using pyscf is as follows:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
from pyscf import gto, dft

mol = gto.Mole()
mol.atom = '''
C     -0.887089    0.175064   -0.012535
C      0.460489   -0.515516   -0.046535
O      1.442965    0.307267    0.565572
H     -0.847478    1.127768   -0.550817
H     -1.658782   -0.453327   -0.465842
H     -1.176944    0.403676    1.018306
H      0.768712   -0.724328   -1.075460
H      0.419486   -1.462073    0.500177
H      1.478640    1.141468    0.067135
'''
mol.basis = '6-31g'
mol.build()

mf = dft.RKS(mol)
mf.xc = 'b3lyp'

energy = mf.kernel()
print(f"B3LYP Energy: {energy} Hartree")


其结果为：

The result is:

B3LYP Energy: -154.98821705367368 Hartree

量子化学计算所需时间大约为8s，比psi4慢上一些。

The quantum chemistry calculation takes about 8 seconds, which is slightly slower than psi4.

这边让这两个软件包进行一个更耗时的任务：优化苯甲酸的几何结构，来判断这两个软件包是否足够实用。苯甲酸的初始结构来源于此。

We then tasked both software packages with a more time-consuming job: optimizing the geometric structure of benzoic acid to determine their practicality. The initial structure of benzoic acid was sourced here.


对于psi4，代码如下：

For psi4, the code is as follows:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import psi4

molecule = psi4.geometry("""
  H      2.2291     -0.1517      2.9720
  C      1.5782     -0.1115      2.0923
  C      0.1964     -0.0808      2.2505
  H     -0.2402     -0.0983      3.2547
  C     -0.6352     -0.0273      1.1374
  H     -1.7228     -0.0023      1.2721
  C     -0.0788     -0.0060     -0.1439
  C     -0.9318      0.0432     -1.3534
  O     -0.6482     -0.2725     -2.4972
  O     -2.1954      0.5046     -1.1851
  H     -2.6479      0.4922     -2.0230
  C      1.3094     -0.0420     -0.3022
  H      1.7464     -0.0340     -1.3085
  C      2.1329     -0.0923      0.8167
  H      3.2205     -0.1183      0.6913
""")

psi4.set_options({'basis': '6-31G'})

energy, wfn = psi4.optimize('B3LYP', molecule=molecule, return_wfn=True)

print(energy)
geometry = wfn.molecule()
print(geometry.save_string_xyz())

计算耗时3分钟，最终优化后的能量为-420.7053888500819。

This calculation took about 3 minutes, resulting in a final optimized energy of -420.7053888500819.


pyscf不能做几何优化，但其可以引用一些几何优化包的代码。这里使用geometric包，其全部代码如下：

PySCF is unable to perform geometry optimization on its own, but it can leverage code from geometry optimization packages. Here, we use the geometric package.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
from pyscf import gto, dft

mol = gto.Mole()
mol.atom = '''
  H      2.2291     -0.1517      2.9720
  C      1.5782     -0.1115      2.0923
  C      0.1964     -0.0808      2.2505
  H     -0.2402     -0.0983      3.2547
  C     -0.6352     -0.0273      1.1374
  H     -1.7228     -0.0023      1.2721
  C     -0.0788     -0.0060     -0.1439
  C     -0.9318      0.0432     -1.3534
  O     -0.6482     -0.2725     -2.4972
  O     -2.1954      0.5046     -1.1851
  H     -2.6479      0.4922     -2.0230
  C      1.3094     -0.0420     -0.3022
  H      1.7464     -0.0340     -1.3085
  C      2.1329     -0.0923      0.8167
  H      3.2205     -0.1183      0.6913
'''
mol.basis = '6-31g'

!pip install geometric
from pyscf.geomopt.geometric_solver import optimize

mf = dft.RKS(mol)
mf.xc = 'b3lyp'

mol_eq = optimize(mf, maxsteps=100)
print(mol_eq.atom_coords())

它大约计算了11分钟，能量为-420.705454579。

This process took roughly 11 minutes, with an energy of -420.705454579.


在CCCBDB数据库中，苯甲酸的b3lyp/6-31G的能量为-420.705447。两个软件包的能量误差均较小，均可以使用。

In the CCCBDB database, the B3LYP/6-31G energy for benzoic acid is -420.705447. Both software packages show a small energy error and are considered usable.


综上，psi4似乎执行速度更快，但是pyscf的输出比较全，而且更灵活。

In summary, psi4 appears to be faster in execution, but PySCF offers more comprehensive outputs and greater flexibility.