Problems with Wavefunction object and fchk() interface

It appears that post-HF methods are not properly updating the Wavefunction object that is being passed to the fchk() interface. Users of fchk() may be unwittingly passing the wrong information to external programs for analysis. The following input writes several .fchk files whose content may be compared to see the problems:

memory 2 gb

molecule {
0 1
f 1 1.4
symmetry c1

set {
basis aug-cc-pvtz
guess sad
scf_type df
mp2_type df
cc_type df

nrg, wfn = energy(‘scf’, return_wfn=True)
fchk_writer = psi4.FCHKWriter(wfn)

nrg, wfn = gradient(‘mp2’, return_wfn=True)
fchk_writer = psi4.FCHKWriter(wfn)

nrg, wfn = energy(‘mp2’, return_wfn=True)
fchk_writer = psi4.FCHKWriter(wfn)

nrg, wfn = gradient(‘ccsd’, return_wfn=True)
fchk_writer = psi4.FCHKWriter(wfn)

nrg, wfn = energy(‘ccsd’, return_wfn=True)
fchk_writer = psi4.FCHKWriter(wfn)

Whether one uses energy(‘ccsd’) or gradient(‘ccsd’) the information in the .fchk file appears to be that for the SCF, not CCSD, including the method written at the top of the file.

It is already documented that one must use gradient(‘mp2’) to get the MP2 density written to the .fchk file, but in this case the orbital energies and wave function are wrong. It looks like they may be for the perturbed part of the wave function rather than the total wave function.

The “Total Energy” written to the .fchk file is always the SCF energy, not that from the requested method.

If fchk() is just a formatter and the real problem is in the Wavefunction object, then this may affect the Molden interface and anything else that uses the Wavefunction object. It is probably related to the segfault issue with conventional MP2 that I reported on 21 April. I hope this can be fixed in the near term as the work I’m currently doing depends upon a valid .fchk file, at least for the use of DF-MP2 and conventional MP2.

I am currently using the binary resulting from a ‘conda update psi4’ on 28 May.


You are right! I can see that this problem does exist for conventional-MP2, df-MP2 and df-CCSD modules who are not updating their wave-functions properly. The conventional CC modules don’t have this issue though. I am going to create a pull request soon to fix this. Thanks for bringing this to our attention.

Did a pull request get filed to fix this? I didn’t see one…


Sorry for the delay as I was busy with other projects. I will try to open a pull request this weekend.

Hi, is there any good new about this fixing?
I recently enconter the similar problem when using 1.2.1 version. The attemp to output the df-ccsd(t) density in the fchk file just give the same density as SCF result.

@bozkaya, lots of appreciations for your great contribution to the dfocc module, this density update issue may also need your attention. The dfocc calculation doesn’t update the density outputted in fchk file.

Is there any workaround to print the “total df-ccsd(t)/df-occsd(t) density” in the format like those in fchk file? Then a manual replacement of the density part in fchk flie maybe achieved to facilitate analyses using third-party codes.

Best regards.

The problems mentioned by @abdale were fixed by this pull request:

Thanks for you efforts and contribution.
All the best.

I think the original problem I posted was fixed, but the FCHKWriter interface broke again afterwards, and is still broken in v1.2.1. I reported the new problems with it on this forum on Sept 25-- please see that email. I think the data being written to the fchk file may be mostly correct (but I haven’t checked df-ccsd), but the data labels are wrong, so 3rd party apps cannot parse it.

Dear wanJ,

Thank you for drawing my attention to this issue. I will try to consider your suggestions in future releases. I will take a look to fchk format, then I may add a piece of code to dot it.

Best regards,

I just found that in v1.2.1, there is still a bug when writing uhf mp2 densities: Sample code:
molecule {
0 1
H 0.000000 0.000000 0.000000
H 0.000000 0.000000 1.900000
symmetry c1

set basis sto-3g

set reference uhf
grad, wfn = gradient(‘mp2’, return_wfn=True)

grad, wfn = gradient(‘scf’, return_wfn=True)

set reference rhf
grad, wfn = gradient(‘mp2’, return_wfn=True)

grad, wfn = gradient(‘scf’, return_wfn=True)

The total densities are:
UHF MP2 8.77062824e-01 8.77062824e-01 8.77062824e-01
UHF SCF 8.77062824e-01 8.77062824e-01 8.77062824e-01
RHF MP2 9.04428342e-01 6.81830761e-01 9.04428342e-01
RHF SCF 8.77062824e-01 8.77062824e-01 8.77062824e-01

Thus, using RHF, the mp2 density is written correctly to the checkpoint file, but using UHF, the written density is actually the SCF density even though it is termed “Total DF-MP2 Density”

Could someone fix this (or any ideas how I could change the code to get the correct checkpoint files until the problem gets fixed later on)?

Best regards,

We literally started working on this yesterday, so expect a fix in the near future.

1 Like

A status update: The current 1.3 release candidate fixes the FCHK based issues, but there are some underlying issues with wavefunctions that will still need to be fixed. I think the discrepancy with UHF SCF density being written in place of UHF MP2 density remain.

Instructions to get the release candidate are here. Try it out and let us know if there are more bugs to squash before the official 1.3 release.

Still not working in this latest version.

MP2 and CI jobs both write sections labeled “CC Density”. Instead, MP2 jobs should write “MP2 Density” and CI jobs should write “CI Density”. “CC Density” is for coupled-cluster jobs only.

Also, all post-HF methods are writing the same density to the “SCF Density” section as to their post-HF section. In the case of MP2 and CC jobs, it looks like it’s just the SCF Density. In the case of CI jobs, I think it’s the CI density, but in that case the CI density shouldn’t be written to the “SCF Density” section.

Thanks. I reopened the issue. Feel free to add anything there you think would be helpful.

I suggest writing a job file that will loop over all basic levels of theory (HF, DFT, MP2, CI and CC), for both open and closed-shell (using H2O for closed shell, CH3 radical for open-shell) and then running an energy calculation and writing the fchk file for each. This will give you a single, multi-step job file that will quickly generate all the basic fchk file types that need to be checked to verify that this feature is working. Then use:

grep “[A-DF-Z]” <fchkfile>

to pull out the headings so that you can quickly verify that they are correct for each job. (omit ‘E’ in case it is used in the exponents, though it might be lower case)

In the HF jobs, set e_convergence and d_convergence to 1e-8, that is, to post-HF values. This will tell you what the density matrix elements should be for HF, and what they should not be for post-HF. If you see those same numbers in the post_HF density, then that job is writing the wrong data.

As you assumed, the UHF SCF density if still written instead of the UHF MP2 density. Would be very happy if that would be fixed in the near future. Thanks!

Trying to revive the thread here. What is the current status on these issues?

I haven’t tested since the last official release. I hope that this issue is resolved before the next release, so that there is no risk of anyone publishing bad data unwittingly. If one intends to analyze the CC density but Psi is actually passing the SCF density to fchk(), the resulting analysis won’t necessarily reveal that. In my email of Feb 19 I suggested a way of doing the testing that should be reasonably thorough.

I don’t think there have been any changes.

I’ll take a more detailed inventory of where things are at when I finish another Psi project I’m working on. Ping me in 72 hours, if I haven’t posted back here by then.

1 Like

I’m assuming that all details posted after Feb '19 have been resolved.

In that case, the only issues are:

  1. The correlated density needs to be labeled as MP2/CI/CC rather than the current catchall “CC Density”
  2. MP2/CI/CC are currently writing the same density to HF and post-HF. Investigation will probably be necessary to determine which is actually being written.

Fixing this doesn’t look bad. That said, I am very busy right now. Among other things, the previously mentioned code refactoring is taking longer than expected… I’ll try to get somebody else to deal with this. If that fails, I can deal with this myself after the current refactoring project. Either way, I make no guarantees about this getting in 1.4.

(Technical notes for other developers: Use to determine which, if any, label should be used. Pull the HF density from reference_wavefunction. If reference wavefunction and current wavefunction’s density match, assume the correlated density was not calculated and skip density writing to the file. We may also need to add name attributes to MP2/CC/CI wavefunctions, and get the density set on the wavefunction. I’m going to need the density set on the wavefunction anyways for another project anyways.)