Hi @mbaczyk,

Thank you for your question.

I have looked into your code and I have the following comments.

First, I don’t see any problem in the way you are building the Hamiltonian. Perhaps you can save few lines of code by using the molecular_hamiltonian function available in PennyLane. In your example, this can be done as follows:

```
bohr_angs = 0.529177210903
symb = ["Li", "H"]
coords = np.array([0.0, 0.0, -0.39870062, 0.0, 0.0, 1.19046204])/bohr_angs
H, qubits = qml.qchem.molecular_hamiltonian(symb,
coords,
method='pyscf',
active_electrons=2,
active_orbitals=5,
mapping='bravyi_kitaev')
```

My second observation is probably more important. You have to be consistent with the chosen basis to encode the molecular wave function. If you use the occupation number basis, where the states of the qubits encode the occupation number of the spin molecular orbitals, you have to use the Jordan-Wigner transformation to map the fermionic to the qubit operator. Similarly, if you build your Hamiltonian using the Bravyi-Kitaev (BK) transformation you have to encode the molecular state using the BK basis. In the example code you provided these two representations are mixed as the defined quantum circuits rely on the occupation number basis encoding. In short, if you want to work with the BK Hamiltonian you would need to re-define the circuits to prepare the trial states in the BK basis. The transformations that can help you to do this can be found in this paper (e.g. Eq. (24):

The Bravyi-Kitaev transformation for quantum computation of electronic structure

I hope this helps. If you have further questions, do not hesitate to get back to us.

Thank you.