Strange (DF-)MP2 forces when point charges are included

[jxzou]

Dear Molpro developers and users,

Recently I'm comparing the (DF-)MP2 results in Molpro with the same method in Gaussian/ORCA/PySCF etc. These software packages are certainly in accord with each other in pure gas phase. But when background point charges are included in an MP2 or DF-MP2 analytical gradient calculation, the obtained atomic forces of Molpro seem strange.

Let's look at a simple example, the linear CO2 molecule which is located on the z-axis, and there are two point charges located in the x-axis:

C   0.0   0.0   0.0
O   0.0   0.0   1.2584
O   0.0   0.0  -1.2584

with point charges of

-4.0   0.0   0.0   -0.1
 4.0   0.0   0.0   -0.1

where the first three columns are positions of point charges, and the fourth column contains the charges.

The non-zero force components are the z-components of two oxygen atoms, and thus we only need to look at these two values in detail. The aug-cc-pVDZ basis set is used in all computations. Electronic energies will not be discussed here since they are always correct. The obtained two z-component forces are shown in the following picture

where the MP2 in Gaussian, MP2 in ORCA, and RI-MP2 in ORCA are extremely close to each other, but the MP2 and DF-MP2 forces of Molpro seem to be different with 3e-4 Hartree/Bohr. I have other large systems which show larger difference, but this example is supposed to be sufficient to illustrate the problem.

Here is my guess: is it possible that the point charges are not taken into consideration in the CP-HF step or during calculation derivatives of H_core?

Thanks for any suggestion!

Best,

Jingxiang

[andreas...@gmail.com]

Dear Jingxiang,

assumingly both the coordinates you have given for the molecule as well as for the point charges are in angstrom (?), and then all should be fine. One can adjust the unit for the lattice charges using:

lattice,<name of file>,unit=[bohr,angstrom]

in the input file (with unit=angstrom already being the default). But I made the following observation: when I use a lattice coordinate file that already is in bohr units (so Molpro multiplies them again with a factor of 1.8897261246257702 because it thinks it is given in angstrom) the result for the z-component of the gradient changes to +/- 0.117970414, much closer to what you get with the other programs (and without the lattice I get +/- 0.117705772). This leads me to think that maybe somewhere you might have used a wrong unit for the coordinates of the point charges. You can probably check if that is the case.

Best wishes,

Andreas

[jxzou]

Hi,  Andreas,

Thanks a lot for your reply. The following answer would be long, before that, here are my short conclusions:
(1) The unit for coordinates of point charges are correct (in A), and these coordinates are used in Gaussian/ORCA/Molpro. These 3 programs use Angstrom as default. These programs have identical HF and MP2 electronic energies, which indicates the point charges have identical effect.
(2) It is possible that the point charges are mistakenly viewed as unit Bohr in the MP2 analytical gradient step somewhere in the Molpro code. The HF energy, HF forces, MP2 energy with point charges included are identical to results from Gaussian/ORCA. Only the MP2 forces are strange/odd.

Now comes the long analysis:
(i) The unit is correct.
To illustrate this, we need some computational details. The content of the point file `Lattice,infile=co2_molpro_mp2.chg1` is
```
Molpro background charge file generated by AutoMR of MOKIT
2
        -4.00000000         0.00000000         0.00000000        -0.10000000   0
         4.00000000         0.00000000         0.00000000        -0.10000000   0
```
The keywords in the input file is `Lattice,infile=co2_molpro_mp2.chg1`. The HF energies calculated by all 3 programs are -187.61454513 a.u., and the MP2 energy is -188.15140634 a.u. Besides, I convert the converged HF MOs from Gaussian to Molpro using [MOKIT](https://gitlab.com/jxzou/mokit), the SCF in Molpro is converged in 1 cycle, which indicates MOs are equivalent and point charges are taken into considerations correctly. If we change the coordinates of point charges (like your lattice coordinates that already is in Bohr units) in Molpro, the SCF cannot be converged in 1 cycle, and the converged HF energy is totally different.

(ii) Only MP2 forces are strange/odd.
I've tested the HF/MP2 energies + forces without any point charge. And I've tested the HF energies + forces, and the MP2 energies with point charges included. The differences among 3 programs are extremely tiny. For example, the MP2 forces without any point charge are +/-0.1177058 Hartree/Bohr for 3 programs. Only the MP2 forces with point charges included are strange. A difference of 3e-4 Hartree/Bohr for MP2 and for this small molecule cannot be viewed as numerical error. So my guess is that the point charges are mistakenly viewed as unit Bohr in the MP2 analytical gradient step somewhere in the Molpro code (but viewed as Angstrom is HF and MP2 energy steps). I have larger examples which show larger force deviations, but introducing them will probably make the discussion complicated.

(iii) A similar force deviation can be found in this example for CCSD with point charges included. I think this problem is VERY important so I have to illustrate it in detail.

Thanks for your patience for reading here. Any further suggestions are welcomed.

Best,
Jingxiang