babbush
commented
7 years ago @Spaceenter: I describe this task some more below. @jarrodmcc: Take a look at what I write here.
The MolecularOperator class is currently a special purpose class for manipulating fermionic operators which come from a molecule. In particular, we use this class to represent molecular Hamiltonians of the form: H = constant + \sum{pq} h{pq} a^\dagger_p aq + 0.5 * \sum{pqrs} h_{pqrs} a^\dagger_p a^\dagger_q a_r as. We could just store this as a FermionOperator but we could also specify it more compactly by just storing the constant, h{pq} and h_{pqrs}. Indeed, this is what the MolecularOperator class does now. In particular, it has the attribute "constant", "one_body_coefficients" and "two_body_coefficients". one_bodycoefficients is a numpy array that is N by N and stores the h{pq} and two_bodycoefficients is a numpy array that is N by N by N by N and stores the h{pqrs} values.
This class supports some convenient things such as myham = MolecularOperator(stuff) h{12} = myham[1,2] h{1234} = my_ham[1,2,3,4]. You can also set elements like this.
A good example of why we need a class completely different from FermionOperator can be seen in the function rotatebasis. Basically, one can rotate the orbital basis of the system with the unitary matrix R = e^{i W} where W = W^\dagger = \sum{pq} w_{pq} a^\dagger_p a_q. So that H' = R^{-1} H R. Naively, one might need N^8 operations to perform this basis rotation. It would be terrible to perform directly on a FermionOperator. But if you supply R to MolecularOperator.rotate_basis(R) it performs the basis transformation in time N^5 using tricks from this website: http://sirius.chem.vt.edu/wiki/doku.php?id=crawdad:programming:project4 By the way, that website has a lot of useful information about this sort of stuff.
Likewise, we have other routines such as the very complicated .jordan_wigner_transform() which are specialized to MolecularOperators for performance. For instance, one can do: fermion_operator = my_operator.get_fermion_operator() qubit_operator = fermion_operator.jordan_wigner_transform() but it is much slower than the direct qubit_operator = my_operator.jordan_wigner_transform(). Another crucial function on the MolecularOperator helps restrict one to an active space. It would also be nice to have a function on MolecularOperator which iterates through just the unique terms (see http://sirius.chem.vt.edu/wiki/doku.php?id=crawdad:programming:project3:hint3.1)
Currently, we are "abusing" the MolecularOperator class and using it for something else in addition to molecular Hamiltonians. We are also using it to store molecular RDM (reduced density matrices): in particular, the 1-RDM and the 2-RDM: D_{pq} = <psi| a^\dagger_p aq |psi> D{pqrs} = <psi| a^\dagger_p a^\dagger_q a_r a_s |psi> Clearly, the RDMs are related to the Hamiltonians and it makes sense for them to be very similar data structures. However, there are some things about these classes that are quite different.
For instance, you can't just apply the Jordan-Wigner transform in the normal way to the RDM. If you call jordan_wigner_transform() on the RDM, you get something that doesn't make sense. Also, you should not try to map a RDM to a fermion operator. So the RDM should not have the method .get_fermion_operator(). But one thing you do want to be able to do is evaluate the expectation value of a QubitTerm / QubitOperator with the molecular RDM. So the molecular RDM should have the basic functionality of the MolecularOperator including all get/set, math and string methods and rotate_basis. There also needs to be a function which gets the RDM in an active space, but we need @jarrodmcc to help with that. The RDM should not have jordan_wigner_transform() on it either. Once the class is made I will also add functionality to the MolecularRDM class.
I want to point out that I'm not 100% happy with MolecularOperators being a class that only works for molecules. We could probably generalize it to something that will work for all fermion systems that have, at most, two-body interactions. But maybe that is not worth the effort. I suggest we have MolecularCoefficients be a fairly general class and then have MolecularRDM and MolecularOperator be classes that inherit from MolecularCoefficients. Does this make sense? Jarrod, do you have any thoughts?
Spaceenter
We need to split MolecularOperator into two classes: one that is specialized for RDMs and the other which is specialized for Hamiltonians. We could have both of these inherit from a common parent which is similar to the current MolecularOperator class.