Use frozen core with orbital optimized methods?

I would like to know if it is possible to evaluate the frozen core correlation energy when using orbital optimized methods (oo-mp2, etc.) after the orbital optimization has completed. It’s customary to do this with non-op methods (mp2, ccsd, etc.) when using valence basis sets and it seems like this might make sense with oo methods too. Is there something I can set in my input to try this out?

Is set freeze_core true what you’re after?

Not exactly. I want to optimize the orbitals including the core, then evaluate the frozen core mp2 energy with the optimized orbitals.

Ah, I see. This actually sounds like a good use-case for the grand information-passing refactoring we just did. We did the top-level refactoring, and I’m not sure one of the modules you need is fully bought-in. Will make inqueries.

Due to some under the hood option passing issues this may not be possible to do easily. However, you have alerted us to a significant hole in some of our option capabilities.

@bozkaya Do you have any tricks to accomplish this.

I can implement it in the DFOCC module, but I will need to get the number of frozen core orbitals for frozen-core energy from user, since it will be read as zero from wfn object in case of all-electron OO computations. I can create an option for it. Alternatively, I may compute it just like the wfn object. However, I did not understand why you need it.

It seems like a somewhat strange thing to do, to do orbital optimization of all orbitals, and then freeze the core and do another computation with frozen core. On the other hand, it should be possible. I think what we should do is the all-electron computation, retrieve the wavefunction for this (with the new wavefunction-passing mechanisms just added to Psi4), and then use that wavefunction for input into the frozen-core computation. We would need a way to set no frozen core in the first part of the computation, and then turn on frozen core in the second part of the computation. Frozen core should always be accepted from user input when present, and that should override any defaults that might be in the wavefunction object.

BTW we might also want to allow orbital optimization for all the orbitals but the frozen core. This is occasionally useful in, say, MCSCF, where sometimes core orbitals mix with valence orbitals in strange and undesirable ways. (Although this kind of some-optimized, some-frozen scheme messes up the orbital response terms if you want to do analytic gradients).

My original question resulted from a conversation with a student about the possible usefulness of frozen core orbital optimized MP2 (or CCD) and how best to perform such a calculation. This led to my earlier question, as we were interested in looking at some test cases with different choices for the core orbitals.

This question was motivated in part by an apparent change in how frozen core Brueckner doubles calculations are performed in the Gaussian programs. Older versions left the Hartree-Fock core unchanged, whereas the current version rotates them to conform to the Brueckner condition. I don’t have access to the current version of Gaussian, so I haven’t been able to see if this matters much.

I’ve noticed this option in the Gaussian keyword list, but they don’t indicate any advantage this gives unless one is already correlating the core. Do you know if there’s a reference for this in the literature?

-TDC

@bozkaya, I think you can use the existing option FROZEN_CORE. Borrowing @dgasmith’s code snippet

    set freeze_core false 
    omp2_e, omp2_wfn = energy('omp2', return_wfn=True)

    set freeze_core true
    new_wfn = psi4.new_wavefunction(omp2_wfn.molecule(), psi4.get_global_option('BASIS')
    new_wfn.Ca().copy(omp2_wfn.Ca())
    new_wfn.Cb().copy(omp2_wfn.Cb())
    fzc_omp2_e = energy('mp2', ref_wfn=omp2_wfn)

is the way we’d ideally want this to work. But I think your code (and SCF?) would need some adaptation for this to really work.

@crawdad, This is discussed in JCTC 5 (2009) 2687-2693 in the context of using Brueckner doubles in a modified version of W1 theory.

Yea, I don’t think this will work at the moment as I think frozen core is only updated within SCF.

@loriab, I don’t know why it would need to be as complicated as in your example input. I’d prefer something like this:

set freeze_core false
omp2_e, omp2_wfn = energy(‘omp2’, return_wfn=True)
set freeze_core true
fzc_omp2_e = energy(‘mp2’, ref_wfn=omp2_wfn)

and post-HF modules should always let user input override any default frozen core behavior that might otherwise be inferred from the wavefunction that’s passed in.