t_as_lambda option control

[cjcscott]

In some recent work on the solvers I've come across some confusing aspects of the code with respect to CC solver calculations with or without solution of the lambda equations. It seems like we need to properly decide what functionality we support and how to control it, as at the moment the interaction of various options are definitely confusing.
Much of these changes stem from 085372b, where we removed t_as_lambda as a CCSD solver option. IIRC there was some discussion of what behaviour we wanted with these options between automatically applying this approximation if required (with warnings) vs failing with a NotCalculatedError.

I can't remember what the conclusion of this was- we currently do the former, so I'm going to assume that in the rest of this. The CCSD solver has the form

        if self.opts.solve_lambda:
...
        else:
            self.log.info("Using Lambda=T approximation for Lambda-amplitudes.")
            l1, l2 = self.solver.t1, self.solver.t2

so in the absence of solving the lambda equations the l1 and l2 values will be set to t1 and t2 regardless of anything else.

The confusion comes in when then looking at routines to calculate global wavefunction and density matrix values, where we seem to have a fair amount of code trying to implement the other approach, and allow calculation of the t_as_lambda=True expectation values when solve_lambda=True, but with inconsistent approaches to control this option.

Looking for instance at get_global_t1_rhf in ewf/amplitudes.py we have code expecting control of this approximation to be in the fragment option t_as_lambda (can be inherited from emb):

    for x in emb.get_fragments(contributes=True, mpi_rank=mpi.rank, sym_parent=None):
        pwf = x.results.pwf.restore().as_ccsd()
        if for_dm2 and x.solver == 'MP2':
            continue
        t1x = pwf.l1 if (get_lambda and not x.opts.t_as_lambda) else pwf.t1
        if t1x is None:
            raise NotCalculatedError("Fragment %s" % x)

while in ewf/rdm.py various functions have a similar keyword argument for control (that doesn't default to emb.opts.t_as_lambda in actual calls):

def make_rdm1_ccsd_global_wf(emb, t_as_lambda=False, with_t1=True, svd_tol=1e-3, ovlp_tol=None, use_sym=True,
                             late_t2_sym=True, mpi_target=None, slow=False):
...
               l2 = wfy.t2 if (t_as_lambda or fy_parent.solver == 'MP2') else wfy.l2
                if l2 is None:
                    raise RuntimeError("No L2-amplitudes found for %s!" % fy)

In calculation of correlation energies this t_as_lambda option is the difference between the 'dm-t2only' and 'dm' energy functionals.
As far as I can tell thanks to the code in the CCSD solver it's impossible for either of these values to be None and result in an error?

However, the main concern comes when looking at code (thankfully usually in the the slow implementations of various functions) which uses

        l1 = emb.get_global_l1() if not t_as_lambda else t1
        l2 = emb.get_global_l2() if not t_as_lambda else t2

These calls actually become calls to get_global_t1_rhf in ewf/amplitudes.py above, so depend on fragment.opts.t_as_lambda. To obtain the actual l1 value here then requires t_as_lambda=False, emb.opts.t_as_lambda=False, and solve_lambda=True. Someone using this function could be forgiven for thinking that t_as_lambda=False should be sufficient without setting the other parameters (which all default to appropriate values).

I'd say the first thing decide is when in a calculation we want to apply this approximation (in the solver or expectation value calculation).
The next thing would be taking into account that we're automatically applying this approximation where necessary in functions to calculate expectation values. Something like force_t_as_lambda or no_t_as_lambda to make it clear that the result can sometimes already incorporate this approximation regardless of this parameter value might be a little more transparent.
This obviously goes hand-in-hand with ensuring that we have a consistent approach to this between different functions.

(Sorry for the wall of text- I was getting test failures in ewf/test_tailoring.py which I traced back to this ambiguity after some effort and thought it best to raise while I had my head around it!)

[cjcscott]

Looking through more intermediate functions, I think that the two-body energy contribution does actually set t_as_lambda=emb.opts.t_as_lambda if no value is provided but not the one-body contribution, and dm-t2only just sets t_as_lambda=True for the t2 contribution.
This seems to be bourne out by some quick tests on STO-3G benzene (I can share a quick script if needed, atomic fragmentation, dmet bath and everything default except solve_lambda and t_as_lambda) giving correlation energies

Results without solving lambda equations:                         'dm' energy=-0.3481416644704918, 'dm-t2only' energy=-0.3481416644704918
Results with solving lambda equations, emb.opts.t_as_lambda=True: 'dm' energy=-0.35111340688875636, 'dm-t2only' energy=-0.35111340688875636
Results with default settings:                                    'dm' energy=-0.3472626167117462, 'dm-t2only' energy=-0.35111280807688394

I don't know what the intended functionality is here, but I'm slightly surprised that there's no option to reproduce the results of not solving the lambda equations when they've been solved. I think this is because only the two-body energy contribution defaults to emb.opts.t_as_lambda, so the one-body energy will always be computed using the true l1/l2 values if available.

[ghb24]

@maxnus would obviously be the best person to comment on this. My feeling is that t_as_lambda should be something enforced on the expectation value level, and to try to keep things consistent from there. I'll trust your judgement as to what makes sense, and then to try and make it all consistent.

[maxnus]

Hi Both.
This whole setting is a mess, partly since we made some changes in implied behaviour, without updating it everywhere rigourously.

We used to have both solve_lambda and t_as_lambda, both as global options, which can be overruled by fragment-specific, local options. The possible combinations of the options, for each fragment, are then

• (solve_lambda=True, t_as_lambda=False) : calculate and use Lambda-amps
• (solve_lambda=True, t_as_lambda=True) : calculate Lambda-amps but use T-amps
• (solve_lambda=False, t_as_lambda=True) : do not calculate Lambda-amps and use T-amps
• (solve_lambda=False, t_as_lambda=False) : Error

Lastly, for MP2-fragments there are no Lambda-amplitudes, and t_as_lambda should be implied
(the same may be true for CI fragments, i.e. T-amps, that come from variational C-amps).

IIRC, we then decided that the second combination is more of a development setting, where we want to reuse the same calculation to build DMs with or without T=Lambda approximation. In a "real-world" application, only the first and third combinations are useful and we got rid of solve_lambda. The Lambda-amplitudes are calculated if t_as_lambda is False.
In hindsight, I'm not sure this was a good idea, there are two problems with this:

1. we often allow t_as_lambda settings on a per make_rdm call level - but when we run the solver we cannot know if we later on want to build a DM with Lambda-amplitudes. Obviously, we could just build the Lambda-amplitudes late, but this would require either to store or rebuild ERIs.

2. The 3rd combination can be useful in a real application setting - perhaps to verify that the T-as-Lambda approximation gives similar results to solved Lambda-amplitudes, before you perform your larger calculation with t_as_lambda=True, exclusively.

Anyway, I suggest that we actually bring back solve_lambda as a separate, CCSD-specific, option.
If it's False, we do not simply set the L-amplitudes equal to the T-amplitudes within the solver itself; this is efficient in terms of lines of code, but just leads to errors down the line, which are difficult to track.

I would then suggest that we make sure we can calculate Lambda-amplitudes late - solve_lambda = True will then just serve to make the calculation more performant (by avoiding recalculating ERIs), but solve_lambda = False doesn't stop us from getting all the same results. t_as_lambda could become exclusively an expectation value option (i.e. a make_rdm function argument), which is not stored in the fragment options at all. However, I do think it's useful to set this on a per-fragment level, which would be more difficult as a function argument (it could be a dict, but it would be less intuitive). Secondly, t_as_lambda also influences the RDM-energy functional, so might want to keep this as a fragment option for this reason. If we allow both, we should have the function argument overrule the fragment option.