Changes to Information Passing in Psi4 (Jan 2016)

Currently Psi4 heavily uses the global space to define how methods run. In order to facilitate more complex methodology sequences and, in general, have more control over a given calculation we are proposing a series of fundamental changes. It should be noted that this is not the best possible system that can be implemented, but a compromise of what we can do before the release while still stepping in the correct direction.

To begin, let us first observe how a single energy('SCF') call works.

energy('SCF')   -> lib/python/driver.py: def energy
def energy      -> lib/python/proc.py: def scf_helper
def scf_helper  -> psi4.scf(precallback, postcallback) (C++ signature `py_psi_scf_callbacks`)

src/bin/psi4/python.cc

double py_psi_scf_callbacks(PyObject *precallback, PyObject *postcallback)
{
    py_psi_prepare_options_for_module("SCF");
    if (scf::scf(Process::environment.options, precallback, postcallback) == Success) {
        return Process::environment.globals["CURRENT ENERGY"];
    }
    else
        return 0.0;
}

src/bin/scf/main.cc

scf = boost::shared_ptr<Wavefunction>(new RHF(options, psio));

// Skipping a few layers for clarity.

HF::HF(Options& options, boost::shared_ptr<PSIO> psio)
    : Wavefunction(options, psio),   // <--- Wavefunction is first created here.
      nuclear_dipole_contribution_(3),
      nuclear_quadrupole_contribution_(6)
{
    common_init();
}

This is shown to make the following points: 1) this sequence is relatively complex, 2) we obtain global options at several disparate places, and 3) we alternate between passing data objects (options) and obtaining from globals. In addition, when we create a molecule with HF (and HF derived classes) we receive the current global molecule object, further, many Wavefunction attributes that are SCF-agnostic are also set here. This creates some confusion in exactly what is contributing to a given calculation.

Let us also look at the last step of a correlated wavefunction:

src/bin/fnocc/fnocc.cc

boost::shared_ptr<FrozenNO> fno(new FrozenNO(Process::environment.wavefunction(),options));

// Setup and run FNOCC

wfn = (boost::shared_ptr<Wavefunction>)fno;
Process::environment.set_wavefunction(wfn);

return Success;

Here we grab the global Wavefunction again, build a new Wavefunction, and then set the new Wavefunction as global. The original global Wavefunction is completely discarded and we are losing data that is both relatively small and potentially useful for future methods.

---

Proposal

• Wavefunctions are immutable to all functions besides the code that generates them.
• All methods that use wavefunction information should take a Wavefunction and return a Wavefunction. There are notable exceptions to this rule (SAPT, Deriv, etc).
• The returned Wavefunction can (and probably should) be of a derived class.
• Wavefunctions should first be initialized and have a separate compute call that runs the method. Many wavefunctions already do this C++ side, we are moving more of these calls python side.
• Python-side energy function should always return a python float:
  • psi4.energy('SCF', ref_wfn=None, options=None, return_wfn=False)
    • ref_wfn=None means not initialized (i.e., use the current base Wavefunction)
    • options=None means grab Psi::enviroment.options since we're not messing (much) w/options at present (future change)
    • return_wfn=True means return the tuple (E, Wfn), and not the conventional E.
  • The energy function can accept wavefunctions and options objects and additionally return the wavefunction itself if requested.
  • If the wavefunction is not returned in this particular example it will be garbage collected by python.
• Similar to energy; gradient, Hessian, frequency, opt, etc calls should be similarly adapted.
• The Wavefunction base class is not a general storage area. We believe only the following parts are necessary:
  • Molecule
  • Primary Basis (Auxiliary basis?)
  • Orbitals
  • Orbital spaces (Should either be greatly enhanced or removed)
  • Density Matrices
  • Gradient
  • Hessian(future)
  • Value dictionary, e.g., (UHF -> "UHF One-Particle Energy", "UHF Two-Particle Energy", "UHF Spin Contamination Metric", ..., "UHF Total Energy")
• The global options object is only mutable by the input psithon script, local copies should be passed to methods. (future)

---

Example CASSCF computation with CISD NO starting guess

rhf_energy, rhf_wfn = energy('SCF', return_wfn=True)

# Here we could do OMP2 or many other orbital guesses.
cisd_energy, cisd_wfn = energy('CISD', ref_wfn=rhf_wfn, return_wfn=True)
# 8Jan2016 note: getting back a Wfn so most things that act on it should be generic Wfn-taking but through Python duck-typing, advanced applications can still use it as a CIWfn
cisd_wfn.compute_nos()
# 8Jan2016 note: line above is perhaps in future compute_nos(cisd_wfn) which preserves immutability of wfn

casscf_energy, casscf_wfn = energy('CASSCF', ref_wfn=cisd_wfn, return_wfn=True)
casscf_wfn.opdm("A").print_out() # Print the alpha OPDM

# Run a RAS calculation using the same CISD NO's
rasscf_energy, rasscf_wfn = energy('RASSCF', ref_wfn=cisd_wfn, return_wfn=True)

# Current behavior is retained
casscf_energy = energy('CASSCF')

From the energy function the call is:

energy('CISD')  -> lib/python/driver.py: def energy
def energy      -> lib/python/proc.py: def run_detci

lib/python/proc.py

def run_detci(name, **kwargs):

    # Reference/options setup
    # Alternatively this could be moved to the `energy` call itself
    ref_wfn = kwargs.get('ref_wfn', None)
    if ref_wfn is None:
        ref_wfn = scf_helper(name)            # Check if we need to first run SCF

    ciwfn = psi4.CIWavefunction(ref_wfn)
    ciwfn.compute_energy()
    return ciwfn

Obviously not all methods are currently set up this way or, in some cases, cannot be set up to run this simply.

---

Possible Naming Tweak to Above (CDS)

DGAS comment: This is ultimately where we want to head, but are not quite there yet on the backend. The currently proposal can be transparently moved to the following at a later date.

The naming and syntax might be even clearer if we get wavefunctions from a function called something like "wavefunction" (or "wfn"). I think this is more intuitive than getting a wavefunction from a function called "energy". The CISD example would then look like this:

rhf_wfn = wavefunction('SCF')

cisd_wfn = wavefunction('CISD', rhf_wfn)
cisd_wfn.compute_nos()

# Run a CASSCF with CISD NO's
casscf_wfn = wavefunction('CASSCF', cisd_wfn)
casscf_wfn.opdm("A").print_out() # Print the alpha OPDM
casscf_energy = casscf_wfn.energy() 

# If user only wants energy, do something like this
casscf_energy = energy('CASSCF', cisd_wfn)

# if no special guess orbitals to pass in, can just use current behavior:
casscf_energy = energy('CASSCF')

Then we'd just need a little wrapper function called "energy" that would call wavefunction() and return its energy.

---

Questions

• How do we pass JK, DFERI, LibTrans, etc objects?
• We apparently do not use the Success/Failure tags C-side. Do they have a use?
• How do we check if SCF has been run for post-SCF methods.
• Who actually makes the base Wavefunction? Probably should not be SCF.
• Should psivars and psiarrays be in Wavefunction?
• If modules take in and emit Wavefunctions, need some rules on absorbing into emitted Wavefunction any info from input wfn that is (a) still valid and (b) strongly influences emitted Wavefunction. (I.e., for ret_wfn = detci(wfn=SCFwfn) (equiv to detci w/bypass_scf=True at present) and ret_wfn = hf(wfn=GUESSwfn), former ret_wfn should have SCF info from input arg while latter needn't.)
• Should Wavefunctions contain orbital spaces? We can either remove this from the base class so space handling is Wavefunction specific (RHF, UHF, CCWavefunction, CIWavefunction, etc) or make it more general by introduce a dictionary of spaces. See here for current orbital specifications.
• Should a derived class of Wavefunction be named in a specific way, e.g., (RHF -> RHFWavefunction).

More unrelated questions for today

• FAE mentioned a script/code that would convert different representations of orbital spaces. Can we use this to create a module for unified MR space parsing.