Is there a quantum chemistry package which calculates a second-order wavefunction?

I am trying to generate a change in density which corresponds to a second-order wave-function in standard Rayleigh-Schrodinger perturbation theory. I would like to try out a simple test on an atomic dimer. However, I am not sure if such a method is implemented in any of the standard quantum chemistry software packages. Does anyone know of such a package?

EDIT: What I am after is the following wavefunction:

$\Psi = \psi^{(0)} + \psi^{(1)} + \psi^{(2)}$

• What do you mean? Second-order Møller–Plesset perturbation theory (MP2) is available in a wide variety of packages. Be careful, though, since in some packages by default the HF density is used for post-HF methods unless you explicitly tell the program to do otherwise. Oct 14 '15 at 9:04
• And the last sentence is just a complete nonsense. The wave function that corresponds to the second-order corrected energy is the second-order corrected wave function. Thus, it is even more unclear what you are asking. Please, clarify. Oct 14 '15 at 9:07
• Wouldn't the second-order energy correspond to the first-order wavefunction in standard perturbation theory? In that sense, I need the second-order wavefunction corresponding to the third-order energy. Does that make sense? Refer to: en.wikipedia.org/wiki/… Oct 14 '15 at 10:57
• No, an $n$th-order corrected energy corresponds to a wave function corrected to the same order, i.e. to an $n$th-order corrected wave function . The fact that you don't need $n$th-order corrected wave functions (rather just $(n-2)$th-order ones) to calculate $n$th-order corrected energy has nothing to do with that. Oct 14 '15 at 12:24
• And this is why I said that you must careful. For instance, for MP2 calculations in Gaussian, the population analysis, etc. by default is done using the HF density that arises from the HF wave function which is the uncorrected wave function. If you want to use MP2 density that arises from the 2nd-order corrected wave function for analysis, you have to tell Gaussian about it. The keywords are Density=MP2. Oct 14 '15 at 12:40

As @Wildcat already mentioned, here is a Gaussian 09 input that will generate the MP3 density, which should correspond to the full second-order wavefunction, thanks to Wigner's $2n+1$ rule:

%chk=g09_mp3.chk
%mem=2gb
#p mp3(full)/sto-3g density=current pop=full output=wfn

wavefunction

-1 1
O     -2.0481500778    0.0447984581   -0.0601937089
H     -1.1514551262    0.4243541138   -0.0806303110
H     -1.8059176361   -0.8898508193    0.0217239645
Cl     1.2314805228   -0.2273431057   -0.0015462666

g09_mp3.wfn

The MO coefficients corresponding to the fully relaxed wavefunction $\Psi$ are also dumped to the *.wfn file in addition to the normal output, which I think can be parsed with cclib.

• But is not that $E^{(4)}$ that requires $\{ \psi_{i}^{(2)} \}$ in accordance with Wigner's rule? Since to calculate $E^{(3)}$ it should be enough to have $\{ \psi_{i}^{(1)} \}$ ($n=1 \to 2n+1=3$). Oct 15 '15 at 7:48
• Sorry for such a delayed response. You're absolutely correct. Nov 20 '15 at 22:16
• I'm actually a bit hesitant to say that this is the answer, though. If this truly gave the MPn wavefunction, the printed densities would be consistent with Wigner's rule, which they aren't. Nov 20 '15 at 22:21