Ghost atoms use ECP

[agoetz]

Ghost atoms should only carry basis functions (and in case of DFT xc quadrature grid points) but not ECPs. Psi4 does not remove the ECP. Tested for Psi4 Version 1.4a2.dev839+e273d32 installed via anaconda.

Example: HF/def2-SVP for He-Xe dimer at 1 Angstrom separation, Xe as ghost atom.

Expected energy: -2.855863 Hartree Psi4: -2.818142 Hartree

Expected energy is obtained when removing ECP information from basis set file.

Input file used for calculation:

import psi4

geo = """
0 1
he  0.0  0.0 0.0
@xe  1.0  0.0  0.0
"""

method = "hf/def2-svp"
mol = psi4.geometry(geo)

psi4.set_options({
    'scf_type': 'direct'
})

e = psi4.energy(method, molecule=mol)

Output file is attached. he-xe_ghost_hf_def2-svp_sp.log

[loriab]

Thanks for the report -- will investigate.

[JonathonMisiewicz]

Moved to a 1.5 target. This will probably be a good thing to re-examine when we're overhauling ECPs.

[tallakahath]

Is anyone presently working on this? I see it's gotten moved from milestone to milestone, just curious what priority level this is (or isn't). I've gotten some dimer interaction energy mismatches of >0.1 kcal/mol due to this bug in CP-corrected dimer interaction energies, though these errors are still << the error of DFT vs a CCSD(T)-level method, so it's not a big deal, just a little obnoxious.

[JonathonMisiewicz]

I'm not aware of anybody working on this, though it should be a priority... I've recently moved research groups, so my time is spoken for right now. Maybe this could be a good group debugging project for GA Tech? @loriab

[loriab]

Naively, this looks within the scope of a group programming project for GT. Will try to emphasize it for v1.9. Thanks for poking the issue.

[tallakahath]

Adding purely for the sake of search results for anyone else who goes down this rabbit hole:

If you are getting a mismatch between Orca and Psi4 on, say, a dimer interaction energy, and your system contains Xe or I, and you've been tearing your hair out trying to reconcile the difference, this is probably why. Orca removes the ECPs on ghost sites.

[philipmnel]

git diff to nominally fix. Further testing required.

--- a/psi4/src/export_mints.cc
+++ b/psi4/src/export_mints.cc
@@ -151,6 +151,7 @@ std::shared_ptr<BasisSet> construct_basisset_from_pydict(const std::shared_ptr<M
             std::string atomlabel = atominfo[0].cast<std::string>();
             std::string hash = atominfo[1].cast<std::string>();
             int ncore = atominfo[2].cast<int>();
+            printf("%s %s %d is ghost %8.4f\n", atomlabel.c_str(), hash.c_str(), ncore, mol->Z(atom));
             for (int atomshells = 3; atomshells < py::len(atominfo); ++atomshells) {
                 // Each shell entry has p primitives that look like
                 // [ angmom, [ [ e1, c1, r1 ], [ e2, c2, r2 ], ...., [ ep, cp, rp ] ] ]
@@ -166,7 +167,9 @@ std::shared_ptr<BasisSet> construct_basisset_from_pydict(const std::shared_ptr<M
                     coefficients.push_back(primitiveinfo[1].cast<double>());
                     ns.push_back(primitiveinfo[2].cast<int>());
                 }
+                if (mol->Z(atom)) {
                 vec_shellinfo.push_back(ShellInfo(am, coefficients, exponents, ns));
+                }
             }
             basis_atom_ncore[name][atomlabel] = ncore;
             basis_atom_ecpshell[name][atomlabel] = vec_shellinfo;

[JGrantHill]

Is it likely that this would also affect SAPT calculations where an ECP is used? I'm seeing some very odd errors and trying to track down what the issue is.

I'm using 1.8.2.

[loriab]

@JGrantHill yes, I think it's definitely affected. The quicker solution will probably slow down the delta Hartree--Fock step b/c it won't reuse integrals. But eventually that could be restored.

I'm still working on the gradients issue #3066 but I hope to look into this soon.