Relativistic Hamiltonian issues - both DKH and X2C

In testing both the DKH and X2C modules in PSI4 against Molpro, I’ve found some odd behavior.

  1. For DKH, even for the Ne atom the HF energy with cc-pVDZ-DK differs from Molpro by a few tenths of a mEh. When one moves down to Zn, however, with cc-pVQZ-DK the energy differs by more than a Hartree. Has this module been tested recently? (I notice there is no test case for dkh)

  2. For X2C, the energies match those from Molpro very closely until one tries cc-pwCVQZ-DK for Zn atom. The calculation completes but the energy blows up, being too low by about 70 Hartrees. I"m guessing there are linear dependency issues with this set if it uses the fully uncontracted basis for the X2C integrals. Is there a way to specify a separate RI basis for these integrals?

I remember there was some discussion about DKH being an unknown factor as there were no test cases. You can use BASIS_RELATIVISTIC to select the X2C basis, although I am not quite sure I would call that an RI basis so hopefully we are referring to the same thing. Skimming through the code in src/libmints/x2cints.cc I see a few calls to general_invert which calls general LAPACK inversion routines, as this is not a protected inversion some of the numerical instabilities could be coming from here.

@francesco @jturney should probably comment more on this.

This is pretty much an echo of @dgasmith. I’ve not matched Molpro with DKH either. And you should believe the x2c test cases (like https://github.com/psi4/psi4/blob/master/tests/x2c3/input.dat) more than the documentation (which has at least one typo) for now.

We’ve confirmed now that if we use a properly decontracted basis set (i.e., manually removing obvious linear dependencies) and use this for basis_relativistic we can reproduce Molpro’s X2C results. In the DKH case it seems it uses the contracted basis set to calculate the DKH integrals so even if one uses a properly DKH contracted basis set the results really aren’t very accurate. We’re testing using a completely uncontracted basis set from the outset for DKH and my expectation is that this will then agree with Molpro (assuming linear dependencies are removed).

As a sidebar, the decontract option doesn’t seem to work for us when we append that to the basis set name. We get a basis set not found error.

Yep, using a completely uncontracted basis set for DKH calculations reproduces Molpro’s RHF energy nicely.

Tracked the numerical instabilities in X2C down to the general invert in X2CInt::form_X(). As this appears to be a non-symmetric matrix, general power routines will not work and a custom setup is required.

@Kirk As you can tell we do not use or test the relativistic capabilities in Psi4 much. Do you have a list of recommended changes that we can begin working on to make Psi4 a bit more canonical in this regard?

Thanks Daniel for looking so closely into this. Probably the biggest issue right now is that one currently can’t use X2C with correlated methods (see the recent post by my student Rulin). This is a show-stopper for sure. Otherwise getting the linearly dependent functions removed in the uncontracted relativistic basis set via some sort of SVD would be the next target. This latter one can at least be worked around by making these by hand where needed. Also, does the -decontract option work for you? It always crashes for us with a basis set not found error.

Hi Kirk,
I have fixed the first issue, see here. I’ll check with @francesco on the linear dependencies and see what we can do. For the decontract option I think this needed a revamp, maybe @loriab can comment on that one.

Hi Daniel,

I just did a git pull (my student Rulin gave me a heads up that you’d fixed this), but core.o now doesn’t build - perhaps due to a different commit:

[ 98%] Building CXX object src/CMakeFiles/core.dir/core.cc.o
/home/kipeters/psi4/psi4/src/core.cchttp://core.cc: In function ‘bool py_psi_set_local_option_array(const string&, const string&, const pybind11::list&, psi::DataType*)’:
/home/kipeters/psi4/psi4/src/core.cc:786http://core.cc:786:13: error: ‘isinstance’ is not a member of ‘py’
if (py::isinstancepy::list(values[n])) {
^~
/home/kipeters/psi4/psi4/src/core.cc:786http://core.cc:786:36: error: expected primary-expression before ‘>’ token
if (py::isinstancepy::list(values[n])) {
^
/home/kipeters/psi4/psi4/src/core.cchttp://core.cc: In function ‘bool py_psi_set_global_option_array(const string&, pybind11::list, psi::DataType*)’:
/home/kipeters/psi4/psi4/src/core.cc:834http://core.cc:834:13: error: ‘isinstance’ is not a member of ‘py’
if (py::isinstancepy::list(values[n])) {
^~
/home/kipeters/psi4/psi4/src/core.cc:834http://core.cc:834:36: error: expected primary-expression before ‘>’ token
if (py::isinstancepy::list(values[n])) {

Looks like this is due to a commit in core.cc by @loriab on 11/20

That looks like you need a full recompile as PyBind11 (a core library) was updated. Also, my changes will not be in the master branch for a few days yet. You can pull directly from my fork if you wish.

I reverted the small changes to core.cc to get that working. I then realized that I hadn’t really pulled your changes yet after all, so I manually made your changes based on the diffs posted on GitHub. I think it’s all working now. Once your changes clear, I’ll bring my release up to date and rebuild the whole shebang. Thanks for getting this sort out so quick!

The decontraction procedure should work insofar as this test case demonstrates (it’s “decon”, rather than “decontract” now). But decontraction fully tested was one of my tasks that’s half finished on a branch, outstripped by all the infrastructure changes, so I’m glad to collect any samples you think ought to be running that don’t.

1 Like

Thanks Lori, it does seem to be working now. Perhaps my last git pull did the trick (or I was incorrectly using -decontract). I"ll let you know if we can break it :slight_smile:

Kirk, I have not checked the forum lately and I saw this msg only now. I’ll be happy to work on the removal of linear dependencies in the X2C code during this break. For most part it should be straightforward but it is not clear to me what to do if the linear dependency is also present in the contracted basis. In that case I think it is mandatory that both at the X2C and SCF level linear dependencies are removed in a consistent way. If you have any suggestion/solution I’ll be happy to hear it.

For DKHn methods, like you said, the current implementation in psi4 does not decontract the computational basis. I am only maintaining the X2C code so I am not sure if and when this issue will be fixed.

Also, @dgasmith thanks for enabling X2C for post-SCF methods. I think this functionality existed in older version of psi4 but I did not check it in the newer refactored versions.

At this point I’m not worried about linear dependencies in the contracted basis set. In principle those should not exist, but in the worse case I haven’t seen eigenvalues of the overlap smaller than 1e-7 with correlation consistent basis sets and eigenvalues that small are rare (e.g., using diffuse augmented H-atom basis sets in cyclic systems). I know that Reiher worked hard in the Molcas and Molpro versions to make the linear dependency problem in the uncontracted basis more or less vanish, but it took a bit of work if I remember correctly. At least for my group, we’re switching over completely to X2C, so I’m not so concerned about DKH.