Df-lccd and df-ccd broken

1

I updated my binary distribution today and find that df-lccd and df-ccd, which previously worked fine, now fall over. I hope this is just a matter of a few lines of python and that this functionality can be restored in future updates.

2

Yes, I got rid of df- and cd- prefixes to energy(), etc. Just use

set cc_type df
energy('ccd')

http://psicode.org/psi4manual/master/proc_py.html#id4 shows what’s available for methods that one can get through Psi4 through multiple codes. ccd follows the same pattern.

3

Thanks, Lori! That works fine for CCD, is there a way to do LCCD? energy(‘lccd’) falls over now, and energy(‘cepa(0)’) does a QCISD(T) calculation (!)

4

Oh right, the codes were using lccd for non-orbital-optimized and ocepa0 for orbital-optimized. So I standardized to cepa(0) all around. For non-OO that calculation is also do-able by another module (fnocc; DePrince; QCISD(T) header) and defaults to that other module. So two items:

  • adding set qc_module occ will switch to the non-default module for doing the lccd calc and will get you the output you’re accustomed to seeing
  • the header of the calculation you got may say QCISD(T), but is the cepa(0) energy not in there somewhere? or is the problem that the cepa(0) energy is not what’s being returned? if either case, thanks for bringing it to my attention.
5

Using set qc_module occ gives a result that is almost, but not quite correct:

On a small test case, CH4 using cc-pVTZ:

psi4 previous version, using
energy(‘df-lccd’)

=========================================================
================ LCCD FINAL RESULTS =====================

Nuclear Repulsion Energy (a.u.) : 13.51348362938091
SCF Energy (a.u.) : -40.21346131953246
REF Energy (a.u.) : -40.21346131953246
DF-LCCD Correlation Energy (a.u.) : -0.22364568317991
DF-LCCD Total Energy (a.u.) : -40.43710700271236

psi4 current version, using
cc_type df
qc_module occ
energy(‘cepa(0)’)

=========================================================
================ CEPA FINAL RESULTS =====================

Nuclear Repulsion Energy (a.u.) : 13.51348362938091
SCF Energy (a.u.) : -40.21346131953202
REF Energy (a.u.) : -40.21346131953202
CEPA(0) Correlation Energy (a.u.) : -0.22349651055890
CEPA(0) Total Energy (a.u.) : -40.43695783009092

which differs from my DF-LCCD reference value

Total Correlation Energy -0.223645587040
Total energy -40.437106906610

by a small amount, but enough to make me scratch my head about it.

In the case where I have not used set qc_module occ, the results are
the same if I replace cepa(0) with qcisd(t) in the input.

I don’t mean to create a lot of work for you, I can fall back to the previous
version (0.3.594) and get the results I want.

6

Thanks very much for the numbers. Will investigate and report back.

7

Many thanks! For reference, my test case used this geometry

C 0.000000 0.000000 0.000000
H 0.625675 0.625675 0.625675
H -0.625675 -0.625675 0.625675
H -0.625675 0.625675 -0.625675
H 0.625675 -0.625675 -0.625675

8

I wasn’t able to reproduce your SCF energy, but I was able to match v0.3 series energies with current code.

  • I should have advised you to use set cepa_type df, not set cc_type df. Because of that, you were getting conventional integral results, not the DF-LCCD labels you were accustomed to.
  • The QCISD calc was a mistake. Should have gone to different func w/i the same module.
  • For cepa only, orbital-optimization was not getting set right due to old names (either a multi-directional merge problem or I’m losing it b/c I remember changing those).

Anyway, switch to cepa_type and apply the below to proc.py (or wait for https://github.com/psi4/psi4/pull/270 to go through), and current psi4 should be useable to you again. Thanks very much for bringing up this issue.

diff --git a/share/python/proc.py b/share/python/proc.py
index 2453274..020e438 100644
--- a/share/python/proc.py
+++ b/share/python/proc.py
@@ -452,7 +452,7 @@ def select_cepa_0_(name, **kwargs):
             if module == 'OCC':
                 func = run_occ
             elif module in ['', 'FNOCC']:
-                func = run_fnocc
+                func = run_cepa
         elif mtd_type == 'DF':
             if module in ['', 'OCC']:
                 func = run_dfocc
@@ -1121,12 +1121,10 @@ def run_dfocc(name, **kwargs):
         raise ValidationError('Unidentified method %s' % (lowername))
 
     # conventional vs. optimized orbitals
-    if lowername in ['ugur-mp2', 'mp2', 'mp2.5', 'mp3',
-                     'lccd', 'ccd', 'ccsd', 'ccsd(t)', 'ccsd(at)',
-                     'ccdl', 'ccsdl']:
+    if lowername in ['mp2', 'mp2.5', 'mp3', 'cepa(0)',
+                     'ccd', 'ccsd', 'ccsd(t)', 'ccsd(at)']:
         psi4.set_local_option('DFOCC', 'ORB_OPT', 'FALSE')
-    elif lowername in ['omp2', 'omp2.5', 'omp3',
-                     'olccd']:
+    elif lowername in ['omp2', 'omp2.5', 'omp3', 'ocepa(0)']:
         psi4.set_local_option('DFOCC', 'ORB_OPT', 'TRUE')
9

Thanks again, Lori! I applied your patch and things seem to work properly now.

10

Just FYI, your queries in this thread led to the finding that one module was returning cepa(0) values that included singles and another excluded them. (Disscussion Naming considerations for cepa(0) if curious.) The upshot of it all is that once this code is merged in, energy() calls to lccd will return linear coupled-cluster doubles (no singles) and lccsd will return linear coupled-cluster singles and doubles (conventional integrals, RHF only). cepa(0) will just be an alias to lccsd, and higher-level cepa methods will continue to include singles.

So for your purposes, you’ll be back to lccd calls where you started. Sorry for the detour, though there really are many things better arranged under the hood because of it.

11

This seems very reasonable to me, thanks again for your assistance!

12

Are there any plans to update the manual to reflect these changes? One finds examples of df-mp2 throughout the current online manual even though that syntax no longer works.

13

There is a ticket out for this. It is one of the things holding back the 1.0 release.