Using ITensor for computational not-quantum chemistry

Question

This is an amazing piece of code and I'm delighted to have found it! I think it might be useful for an application I am working on. The goal is to solve the chemical master equation

$$
\frac{\mathrm{d} P}{\mathrm{d} t} =\mathbf{A} P
$$

where @@\mathbf{A}@@ is a matrix of connections and @@P@@ represents the vector of probabilities of being in each state. (See this Wikipedia page for more details on master equations.) I got the idea for solving the master equation using the tensor formalism from this paper and its companion paper. Our system is similar in form to the CO oxidation model in the papers, but it has different chemical species and more complicated reactions.

To model this, I think I can use ITensor by applying an MPO (which may act on more than two sites!) to an MPS repeatedly and, at each timestep, measure the properties we are interested in (say, total number of atoms/molecules of each type). I haven't seen an example of this being done in ITensor yet, but this Wikipedia page suggests that a reaction-diffusion system could be modeled as a time-evolving set of non-quantum states.

Has anyone done something like this in ITensor before? Are there any fundamental problems with using the library this way, particularly since this isn't exactly the canonical use case?

I also understand that toExpH doesn't work for MPOs that operate on more than two sites -- which, sadly, mine probably does. Barring the obvious numerical instabilities, is there any reason why I couldn't do something like an explicit Euler integration instead of using the matrix exponential form to solve the above equation?

Comments

Hi emprice, I may not have time this week for a longer answer, but just wanted to say that a project along the lines you describe would certainly be implementable in ITensor. I would argue it would be the best library to implement it in. However, many of the specific steps you'd need to do aren't implemented in an automatic way in ITensor (nor would they be in any other library). So you would have to do significant algorithm development to carry out some of these steps.

For example, time evolving a density matrix (represented as an MPO) using Trotter gates is not significantly different from time evolving an MPS wavefunction *if* the interactions are nearest-neighbor. If they are long-ranged and involve fermions, then you'd need to use a method to compress the interactions and deal carefully with the Jordan-Wigner string needed to capture the fermion minus signs.