Is there a quantum chemistry package which calculates a second-order wavefunction?
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.