Adding two- and three-body operators

[farrenal]

It seems that whenever we add operators in the order 2b + 3b, the 3b is simply ignored whereas in the order 3b + 2b the correct addition is performed. The same issue persists when using += in python. See example below.

#!/usr/bin/env python3

from pyIMSRG import *
import numpy as np

hw =            16
cLS =           -2.00
emax =          3

imsrg2or3 =     '3n7'   # 2, 3n7 or 3

ref =           'O14'
val =           ref

############ Create model space ###############################################
ms = ModelSpace(emax,ref,val)
ms.SetHbarOmega(hw)
ut = UnitTest(ms)

############ Create operators #################################################
H = Operator(ms,0,0,0,2)

H = Operator(ms,0,0,0,2)
H = OperatorFromString(ms, "VMinnesota")
H += OperatorFromString(ms, 'Trel')
H += cLS * OperatorFromString(ms, "LdotS")
H += OperatorFromString(ms, 'VCoul')

T2HO = OperatorFromString(ms, 'Iso2')

HNO = H.DoNormalOrdering()
T2NO = T2HO.DoNormalOrdering()  
HHF = HNO
T2HF = T2NO

############ Decide IMSRG 2,3f2,3n7 or 3 ######################################
if imsrg2or3 == '3':
    Commutator.SetUseIMSRG3(True)
    HHF.ThreeBody.SetMode("pn")
    T2HF.ThreeBody.SetMode("pn")
    
if imsrg2or3 == '3n7':
    Commutator.SetUseIMSRG3(True)
    Commutator.SetUseIMSRG3N7(True)
    HHF.ThreeBody.SetMode("pn")
    T2HF.ThreeBody.SetMode("pn")
    
if imsrg2or3 == '3f2': 
    BCH.SetUseFactorizedCorrection(True)
    Commutator.FactorizedDoubleCommutator.SetUse_1b_Intermediates(True)
    Commutator.FactorizedDoubleCommutator.SetUse_2b_Intermediates(False)
    
############ Show Bug #########################################################
### First example
clean2b = 0*HHF
clean2b.SetAntiHermitian()
induced3b = 0*HHF
induced3b.SetAntiHermitian()
Commutator.comm223ss(HHF,T2HF,induced3b)

print('\n---\n')
print('[H,T2]_3b norm')
print(induced3b.ThreeBodyNorm())

print('\n2b + 3b')
supposed3b = clean2b + induced3b
print(supposed3b.ThreeBodyNorm())

print('\n3b + 2b')
fixed3b = induced3b + clean2b
print(fixed3b.ThreeBodyNorm())


print('\n---\n')
### Second example
print('2b += 3b')
clean2b += induced3b
print(clean2b.ThreeBodyNorm())

print('\n2b = 3b')
clean2b = induced3b
print(clean2b.ThreeBodyNorm())

[ragnarstroberg]

It looks like the issue is caused by calling HHF.ThreeBody.SetMode("pn"), which allocates the 3N storage, but doesn't set the particle rank to 3. If I explicitly set the particle rank to 3, the correct addition happens with either ordering.
Beyond that, I haven't tracked down the issue.