Inconsistency in naming convention of components of the quadrupole and octupole moments in TCE module

[Mith13]

Hi,

I am currently working on modifying the TCE code to calculate the quadrupole and octupole polarizabilities and met an inconsistency in name ordering of the quadrupole/octupole components.

For example "axisname" array with names of multipole moments is defined as

tce_energy.F:287-288

      data axisname/'X','Y','Z','XX','YY','ZZ','XY','XZ','YZ'
     & 'XXX','YYY','ZZZ','XXY','XXZ','YYX','YYZ','ZZX','ZZY','XYZ'/

and similarly the name of filename stubs for corresponding components is the same (tce_energy.F:278,284-285). However, when the calculated expectation values of multipole moments are printed (tce_energy.F:1010:1039) , the order is changed to XX XY XZ YY YZ ZZ and XXX XXY XXZ XYY XYZ XZZ YYY YYZ YZZ ZZZ.

This inconsistency then occurs in multiple files, for example tce_ints.fh uses the former ordering, while tce_aod1.F works with latter ordering. Based on the fact that tce_aod1.F handles the multipole matrices I presume the correct ordering is the latter. Moreover, I have done some tests comparing the calculation of multipole polarizabilites with DALTON code and figured out that the correct ordering is the latter. It looks like there are two parallel implementations of CC properties and only one is currently working, perhaps the ordering has got changed during the development of the currently non-functional part (tce_props.F) and was never fixed.

The fix seems relatively easy but I would like to ask the developers whether the latter ordering (i.e. XX XY XZ YY YZ ZZ and XXX XXY XXZ XYY XYZ XZZ YYY YYZ YZZ ZZZ) is really correct and that there is no "hidden" stuff that I missed in my assessment of the correct ordering.

Expected behavior The naming order of multipole moments/integrals should be consistent across the TCE module

[jeffhammond]

I wrote that code more than 10 years ago. It was never properly finished and I should have deleted the printout.

Feel free to make it work, and correct in comparison to e.g. Dalton, anyway you like.

To the extent that I remember it now, I was trying to match TCE to the multipole code in other parts of NWChem. This was somewhat tricky because there are two different ways to get the dipole operator. The original TCE code used the one that only works for dipoles. I then tried to integrate with the more general multipole module. However, none of this was ever seriously validated.

Anyways, I'll do my best to help here if I can, but if you want to walk through the code together, we can meet online and I can try to explain what I was thinking when this stuff was written.

The other thing is that I never finished the full refactoring of TCE, so there is a bunch of unused code that might be confusing. Figure out the call tree from tce_energy as best you can.

[Mith13]

I see. I think I got, it just wanted to make sure that the ordering I think is correct is really correct.

Based on my tests with DALTON I think I managed to correctly hack tce_energy.F and ccsd_lr_alpha.F to use higher multipoles than a dipole, at least when compared to helium, hlium dimer and neon calculations.

What I did was that I expanded rr1filename, tr1filename, ... data blocks
and then added proper offsets to their handles when running over axes from 1...19 in tce_energy.F and ccsd_lr_alpha.F, moreover dynaxis had to be modified with offset as well.

I have also modified CCSDT and CCSDTQ code but to validate it I am finishing FCI code as I am not aware of any other public code to compare it with (DALTON goes only up to CC3) and do calculation of 3 and 4 electron system.

When I am finished I can do PR.

I have an additional question, I discovered that probably "you" wrote a part for imaginary frequency polarizabilities code. As I want to later calculate the dispersion coefficients I want to ask if that code ever worked or not so I can reintroduce it and update it for higher poles to get higher coefficients than C6.

[jeffhammond]

You can validate quadrupole and octupole with finite field but it's a huge pain. I don't know if MRCC can do it but that's the other CCSDT and CCSDTQ response code out there, although it's not open-source.

However, I would expect that if you get stuff working for CCSD-LR with quadrupoles and octupoles, then CCSDT-LR and CCSDTQ-LR should be fine, because the nature of the multipole operator is orthogonal to the CC excitation level. Just make sure you validate that the general multipole integrals version of dipole matches the dipole-only version.

The C6 code was never validated and Alex Tkatchenko found some issues with it, which is why it was not published. I went through the code quite a few times and never figured it out. The equations should be right, so I think the solver is the problem. There is a chapter in my dissertation on it if you want the details of the implementation. The numbers seemed reasonable so the error is subtle. I think Alex compared 2-electron systems, so you could start there.

[Mith13]

You are right about the finite field approach, good point. Afaik, MRCC can only do dipole polarizability.

Thank you for the info on imaginary C6 code, I found your phd thesis and I am currently going through it.