Hi Jin,
Thanks for the question. So this is the intended, and expected behavior of the ITensor * operator, namely that first A and B are contracted, then the result of that is contracted with C.
One reason for this is that this is just how operators work in C++. There is no notation in C++ that would allow simultaneous contraction of three or more ITensors in that language using operator overloading (*,+, and similar), apart from defining a function call such as multicontract(A,B,C,D,...)
(which we might do in the future).
Another reason is that this behavior is quite often the desired behavior that one does want, so it is not in any way a bad thing. It just may not have been what you were expecting. The ITensor interface is based on tensor diagram notation which also uses a pairwise contraction convention. Contracting over an index shared between three or more tensors requires the notion of a "hyperedge" in a tensor network graph, which is perfectly fine to introduce but I'm just mentioning it to say it's not the standard thing that diagrams express.
We do offer something like a hyperedge in ITensor though: it is the delta(m,n,p)
ITensor, which is a special ITensor with all its (multi-)diagonal elements set to 1.0. Not all routines involving delta
tensors are as optimized as they could be, though some are and in general it's a useful tool.
Hope that helps -
Miles