Dear PSI4 developers and users,
Recently I’m testing and comparing MP2 analytical gradients among different quantum chemistry software packages. I find that when there are point charges included in a MP2 analytical gradient calculation using PSI4, the obtained MP2 gradients seem strange/odd.
Here I use the linear CO2 molecule to illustrate the problem. cc-pVTZ is used in all calculations. No RI/DF approximation is used. No core orbital is frozen (i.e. all-electron calculation). The Cartesian coordinates (in Angstrom) are
C 0.0 0.0 0.0
O 0.0 0.0 1.2584
O 0.0 0.0 -1.2584
and the point charges (in Bohr) are
pt_chg = [
[ 0.100000,[ -5.66917837, 0.00000000, 0.00000000]],
[ 0.100000,[ 5.66917837, 0.00000000, 0.00000000]],
[ 0.100000,[ 0.00000000, -5.66917837, 0.00000000]],
[ 0.100000,[ 0.00000000, 5.66917837, 0.00000000]]
]
The 1st column contains the charges, and the 2nd~4th columns contain their coordinates in Bohr (since PSI4 external_potentials requires the usage of unit Bohr). For convenience, positions of point charges in Angstrom are shown below
-3.0 0.0 0.0
3.0 0.0 0.0
0.0 -3.0 0.0
0.0 3.0 0.0
The last few lines in the PSI4 input file are
...
set {
guess read
e_convergence 1e-8
d_convergence 1e-6
}
set num_frozen_docc 0
set mp2_type conv
gradient('mp2',external_potentials=pt_chg)
Now comes to the results. This linear CO2 molecule only has non-zero gradients on the z-components of two O atoms. To see it more clearly, here is the result of PSI4
-Total gradient:
Atom X Y Z
------ ----------------- ----------------- -----------------
1 -0.000000000000 -0.000000000000 -0.000000000102
2 0.000000000000 0.000000000000 0.135559687833
3 -0.000000000000 0.000000000000 -0.135559687730
And here are results (in Hatree/Bohr) collected from various programs
| program | gradient |
|---|---|
| Gaussian | ± 0.138602256 |
| ORCA | ± 0.138601671 |
| PySCF | ± 0.138602251 |
| PSI4 | ± 0.135559688 |
One can see that the PSI4 MP2 analytical gradient has 3e-3 Hatree/Bohr difference with results of other packages, which cannot be viewed as numerical error.
The RHF energy and MP2 energy of various programs are in excellent agreement (in Hartree)
| program | E_RHF | E_MP2 |
|---|---|---|
| Gaussian | -187.64278870 | -188.31236705 |
| ORCA | -187.64278870 | -188.31236705 |
| PySCF | -187.64278870 | -188.31236705 |
| PSI4 | -187.64278870 | -188.31236705 |
This means that the unit of positions for atoms and point charges are correctly taken into considerations. Note that PSI4 HF and MP2 energies in the output file do not include the self-interaction of point charges 0.00675305 a.u. After adding this term, one can compare the PSI4 energy with that of Gaussian.
Here are some PSI4 input/output files in case that anyone wants to reproduce the calculation (please remove .txt suffix and open it as .zip package)
strage_psi4_mp2_force.zip.txt (79.5 KB)
This problem does not exist for the RI-MP2 analytical gradients (in Hatree/Bohr), e.g.
| program | gradient |
|---|---|
| ORCA | ± 0.138515736 |
| PSI4 | ± 0.138475826 |
I only test PSI4 v1.8.2 and v1.9.1, and this problem occurs in both versions. Similar problems exist for CCSD gradients when point charges are included. Would you please be so kind as to figure it out? Thanks a lot.
Best,
Jingxiang