Is it possible to multiply AutoMPO objects before converting to IQMPOs?

Question

I want to be able to do something like

auto A = AutoMPO(sites);
auto B = AutoMPO(sites);
//Add operators to A and B
auto C = A*B;
auto myIQMPO = IQMPO(C);

I want to be able to do this because I have operators A and B that don't individually conserve quantum numbers, but their product does. So I can't do something like

auto A = AutoMPO(sites);
auto B = AutoMPO(sites);
//Add operators to A and B
auto myIQMPOA = IQMPO(A);
auto myIQMPOB = IQMPO(B);
IQMPO C;
nmultMPO(A,B,C);

because I can't construct the IQMPO of A and B individually. Additionally, experimenting with nmultMPO, I find that nmultMPO gives errors like

scale too big, omitting from reported eigenvalues

even when the product does have an efficient representation as a MPO. So I'd like to be able to construct products of AutoMPO operators directly.

Ultimately, I want to be able to take products like A*B*C*D, where A,B,C,D, are sums of a large number of local operators, so just figuring out what these products are in terms of Cup, Cdn, etc quickly gets prohibitively messy.

Is there any way to do this? Or is this feature something planned? Thanks!

Answer 1

Hi, so this is a good idea for a feature that we could create with AutoMPO. However it's not a feature that we currently have. I'd like to add it at some point in the future but I can't promise I can get to it very soon.

Here are some things we could do:

(1) You could hack together your own workaround code that computes the product in your code before feeding the operators to AutoMPO. So like if the first MPO is an operator:

$$\sum_{i} h_{i} $$

and the second is

$$\sum_{j} k_{j} $$

then you could make a for loop that feeds operators into AutoMPO of the form

$$\sum_{ij} h_{i} k_{j} $$

It's not yet documented on the website, but AutoMPO now accepts any length string of operators, so it's no longer limited to just 2-site terms.

(2) If you're feeling ambitious, you could start diving into how AutoMPO works and create a function that builds on top of the AutoMPO system to create the product in the way you are asking. Actually this may not be too tough. Basically the AutoMPO just accumulates an array of string and integer pairs, together with a coefficient that also includes keeping track of any fermionic minus signs. So one could create a fancier function that does the double loop over the operators stored in two AutoMPO objects and merges them into the data stored in a single AutoMPO defined as their product (by data stored I mean the intermediate representation stored in the "AutoMPO::terms_" data member before finally converting everything to an actual IQMPO, which is a part you wouldn't have to change as long as the intermediate representation is valid).

But since as I just mentioned it should be somewhat straightforward for someone who understands the system, maybe it's a feature I can add over the next few weeks.

Finally, please let me know if you think there's a bug in nmultMPO as I should fix that. If the IQMPOs you are putting into it are invalid or ill-formed somehow (as you hinted they might be in your question) then it's probably not a bug if nmultMPO fails. But if they are valid IQMPOs and it still failed could you email me a minimal code that reproduces this bug? Thanks -

Miles

Comments

Thanks, I did the first option and it seems to be working fine. I'll email you about the bug. It's possible I'm doing something wrong, but you can judge for yourself!

---

Great. I'm glad the first option is working because if I made a function to multiply two AutoMPOs, what it would do internally is basically just what you did by hand in your code. So since it works I should go ahead and make that into a feature of AutoMPO.

Ok yes let me know about the MPO multiplication. In general numerical operations on MPOs are tricky and can become ill conditioned in some cases. So there's no "one true algorithm" I know of that always works perfectly for things like multiplying or adding MPOs etc. But we want to make sure nmultMPO works as best as it can.