MP2 Geometry Optimization with Forced Cartesian Functions

[wavefunction91]

Pyscf 2.4.0

I'm attempting to perform a geometry optimization with forced Cartesian GTOs. (EDIT: add script)

import sys,os
import pyscf
from pyscf import mp

geom_fname = sys.argv[1]
basis_name = sys.argv[2]

mol = pyscf.M(atom = geom_fname, output = 'coro_mp2_opt.log', basis = basis_name, cart = True)
mol.max_memory = 20*1024
mf = mol.RHF()
mf.chkfile = 'coro_mp2_opt.chk'
mf.verbose = 4
mf.run()

mp2 = mp.MP2(mf)
mp2.verbose = 4
mp2.run()
opt = mp2.Gradients().optimizer(solver='berny')
opt.verbose = 4

mol_eq = opt.kernel()
print(mol_eq.atom_coords())

HF/MP2 run fine with forced Cartesians, but the MP2 gradients error out with the following:

Traceback (most recent call last):
  File "../coro.py", line 21, in <module>
    mol_eq = opt.kernel()
  File "/home/dbwy/Documents/Projects/TransitionMetalDimers/venv/lib/python3.8/site-packages/pyscf/geomopt/berny_solver.py", line 203, in kernel
    kernel(self.method, callback=self.callback, **params)
  File "/home/dbwy/Documents/Projects/TransitionMetalDimers/venv/lib/python3.8/site-packages/pyscf/geomopt/berny_solver.py", line 143, in kernel
    energy, gradients = g_scanner(mol)
  File "/home/dbwy/Documents/Projects/TransitionMetalDimers/venv/lib/python3.8/site-packages/pyscf/grad/mp2.py", line 236, in __call__
    de = self.kernel(mp_scanner.t2)
  File "/home/dbwy/Documents/Projects/TransitionMetalDimers/venv/lib/python3.8/site-packages/pyscf/grad/mp2.py", line 293, in kernel
    de = self.grad_elec(t2, atmlst, verbose=log)
  File "/home/dbwy/Documents/Projects/TransitionMetalDimers/venv/lib/python3.8/site-packages/pyscf/grad/mp2.py", line 97, in grad_elec
    Imat += lib.einsum('ipx,iqx->pq', eri0.reshape(nf,nao,-1), dm2buf)
  File "/home/dbwy/Documents/Projects/TransitionMetalDimers/venv/lib/python3.8/site-packages/pyscf/lib/numpy_helper.py", line 239, in einsum
    out = _contract(subscripts, *tensors, **kwargs)
  File "/home/dbwy/Documents/Projects/TransitionMetalDimers/venv/lib/python3.8/site-packages/pyscf/lib/numpy_helper.py", line 166, in _contract
    raise ValueError(err)
ValueError: ERROR: In index string ipx,iqx->pq, the range of index x is different in A (88704) and B (73920)

This error goes away if I allow for auto-deduction of the basis (6-31G*) to spherical GTOs. Am I missing something obvious here?

[susilehtola]

Wait, 6-31G* is defined with cartesian D functions (and spherical F functions!)

[wavefunction91]

Correct, but if I don't specify anything, it's generated as spherical... Verified this against NIST CCCBDB MP2 reference values (plus printing nao_nr/nao_cart to discover they're different) - with cart=True I get the right thing, without I don't.

[verena-neufeld]

Hi @wavefunction91 , thank you for reporting this issue. Has it been resolved now via #2001 ?

[wavefunction91]

@verena-neufeld Yep! All resolved on my end. Thanks for addressing