How to plot wavefunction amplitude?

[sgr]

Hi Psi4 Community,

I have what feels like a dumb question, but I can’t figure it out, so maybe you can help me!

Things that work: I am trying to build a tutorial on CI. I am using the psi4API in python (which is amazing). I am motivating the need for correlation with the failure of HF to get singlet and triplet energies correct for H2 dissociation. I can plot singlet and triplet surfaces and triplet is lower than singlet, as expected. Then, CISD fixes the problem like it should. Yea! So far, so good.

I see how to get the updated wavefunction

e,wfn_out=psi4.energy(‘SCF’,molecule=mol,return_wfn=‘true’)

Where I am stuck: I want to show the uncorrelated nature of the wavefunction in HF and the correlation added by CI. I want to plot the two-electron wavefunction amplitude in real-space \Psi(r1,r2) for the HF and CI wavefunctions. I want to show that in HF the amplitude of two electrons to be on the same nucleus is high for all distances but CI fixes that. I don’t see how to get this!

What I have tried: I have looked through the tutorials and Psi4 lab experiments for some hints, but I don’t see anything helpful. I have looked through the documentation that I could find - traversing help(psi4), help(psi4.core), help(psi4.core.Wavefunction), help(psi4.core.MintsHelper), and a bunch of others. I have learned a lot but not the easy thing I am looking for! (At least I don’t recognize it as such!)

This thread is very similar to what I am looking for

Compute property at given 3D coordinates

Dear all, Is there a function which returns one of the properties accessible through cubeprop() (e.g. electron density), but at any arbitrary 3D coordinate, given a wavefunction? I assume the information itself could be easily computed, given that cubeprop() can return this property if you choose the grid correctly, but this of course comes along with the entire cubefile. I’m thinking of something like this: ```E, wfn = psi4.energy(level_of_theory, return_wfn=True) result = pointprop(wfn, [x, …

but it is for the density. I want to put in all the electron coordinates and get the total wavefunction amplitude (3N cartesian coordinates as inputs and one number as an output).

My current work around: I have made a cartoon by showing 1-D cuts through sigma and sigma* MOs (from 1s AOs) and show that the double excitation sigma* (times) sigma* added to the ground state reduces the chance for both electrons to be on the same atom. BUT I would love to show that happening with the CISD wavefunction as a function of r_AB.

Sorry if I missed something obvious! Any help is welcome! Thanks!

[hokru]

Why the wave function? The probability distribution is better shown with the density.
The HF wave function could be plotted I suppose, but a CI wave function is generally too large to store.

[jmisiewicz]

If it’s just H2, then I can store FCI wavefunction coefficients at cc-pVTZ even with my laptop. Unfortunately, there isn’t a way to access the CI coefficients Python-side at the moment. (Obligatory reminder ping of @loriab about making sure this is on CDS’s to-do list.)

Will plotting electron density work for your purposes? Also, which module are you using for this computation: Does your output file say DETCI or FNOCC?

[sgr]

Thanks for the replies!

@hokru, I aim to show the effect of mixing of the ground state wavefunction with the doubly excited wavefunction. At the HF level, both the bonding and doubly excited wavefunctions have equal probabilities for both electrons to be on each nucleus. But, at the wavefunction amplitude level, they have different signs. Adding the two wavefunctions then reduces the joint probability density that both electrons are on the same nucleus. So, at the HF level showing the wavefunction would be optimal.

At the CI level showing the joint probability density would still convey the essential insight. But this would not be the usual density! I am trying to show electron correlation, which depends on all the electron coordinates. I want to plot the amplitude (or probability) that electron two can be at position r2 given that electron one is at r1. That information is contained in Psi(r1,r2).

@jmisiewicz, thanks for the insights. As mentioned above, for the CI calculation yes the density might work – if I can get to the multi electron density | Psi(r1,r2) |^2 .

The output file says DETCI (I just went with the default). I don’t know if FNOCC is an option – I want to show that CI gets the convergence of the singlet and triplet energies at large rAB. There are some words in the documentation that say FNOCC is implemented for closed shell references, that would be a problem for the triplet state, is that right?

Thank you both for your time!

[jmisiewicz]

Yes, fnocc is unusable for the triplet state.

Unfortunately, I’m not aware of any way to compute a 2-electron real space density in Psi. You can get a 2-electron orbital space reduced density matrix, but that isn’t going to help unless you can transform that into a real space object, and I don’t think Psi has that functionality.

[sgr]

okay I will come up with some kind of hand-waving argument! Thanks for your time!

[hokru]

Making a python script to plot a density is not a big problem. Especially for H2 minimal basis.
FCI densities for H2 can be obtained from CCSD densities.

Amplitudes for HF should also be straightforward. Though spontaneously I am a little unsure about the formula/loops.

However, for CI one needs to know also all determinants to place the occupations.
We have python-based CI codes (https://github.com/psi4/psi4numpy/tree/master/Configuration-Interaction) and python-based determinant generation too, so one could cook something up from there perhaps. Not so familiar with that part.

[jmisiewicz]

But how would you go from FCI coefficients or an orbital-space density matrix to a two-electron real-space density in Psi? That’s the part that I can’t figure out.

[hokru]

For the density you do it like this:

where the phi are the basis functions, D the density matrix and r your real space point.

[jmisiewicz]

Well, you’d need the 2RDM for the desired pair density, not the 1RDM, but in principle, it’s the same trick. How would you evaluate the orbitals at points? I imagine it’s possible based on cubeprop, but I’m not sure how. If you can answer that, I can see how you could actually get this amplitude.

[hokru]

cubeprop can plot the (1-particle) density matrix, but needs to made available outside of cube files.

I wrote this python script, it works for s-functions. I can reproduce old H2 data I have. It should work for higher angular momenta, but it’s untested: https://gist.github.com/hokru/9b93fffc093dce0e8af8b9e54c120bb2

The current setup plots the (1-particle) density matrix for H2 along the bond vector (z-axis).
It plots an amplitude which looks qualitatively ok-ish.
python density.py > plot.dat

It’s not a particular nice script, it’s late here and I am watching Space Force on Netflix.

[crawdad]

I do this in ccdensity/densgrid_RHF.cc on the CC one-electron density, and I used the same basic code to compute values of the CC two-electron density in http://dx.doi.org/10.1063/1.1535440 .