Implement VibrationalOp.from_polynomial_tensor

[mrossinek]

What should we add?

#665 proposes the addition of a PolynomialTensor class to handle coefficient storage.

This issue tracks the implementation of the method which can construct a VibrationalOp from those coefficients.
Below, I summarize the proposed structure of the PolynomialTensor.data field:

VibrationalOp

Finally, the vibrational structure use case is a bit more special. Let me first explain briefly, what the vibrational operator does:

• it has m modes
• each mode can have an individual number of modals,

$$ n_i \forall i \in m $$

• this means, the total length of the VibrationalOp is given by

$$ \sum_{i \in m} n_m $$

• the splitting of each mode into its modals is determined by the basis (e.g. HarmonicBasis) which we use to map to the second-quantization space
• thus, the coefficients which we extract from a classical code (e.g. via the GaussianLogResult) actually only index the modes and, thus, have no concept of modals

What this means, is that the coefficients stored in a PolynomialTensor must be combined with a basis in order to produce a full VibrationalOp. We see the same in the current stack, because the VibrationalEnergy requires a basis in combination with the VibrationalIntegrals.

That said, what should we actually store in a PolynomialTensor in the vibrational scenario? I would argue, that the splitting into n-body terms which we do in the VibrationalIntegrals is not very suitable here, because of the following:

1-Body Terms:
        <sparse integral list with 13 entries>
        (2, 2) = 352.3005875
        (-2, -2) = -352.3005875
        (1, 1) = 631.6153975
        (-1, -1) = -631.6153975
        (4, 4) = 115.653915
        (-4, -4) = -115.653915
        (3, 3) = 115.653915
        (-3, -3) = -115.653915
        (2, 2, 2) = -15.341901966295344
        (2, 2, 2, 2) = 0.4207357291666667
        (1, 1, 1, 1) = 1.6122932291666665
        (4, 4, 4, 4) = 2.2973803125
        (3, 3, 3, 3) = 2.2973803125

A simple 1-body term, can have very different lengths. The same holds for n-body terms in general.
Instead, I suggest we store the coefficients in the same format as done in the GaussianLogResult, where we split them into quadratic, cubic and quartic force terms. These all have length 2, 3 and 4, respectively.
The kinetic terms (negative indices in the example above) can be inferred from the quadratic terms by a simple negation of the coefficient, so they do not need to be stored separately.
This results in a structure like this:

{
  "??": [[...]],  # quadratic
  "???": [[[...]]],  # cubic
  "????": [[[[...]]]],  # quartic
}

The strings here do not carry useful information in their characters, as only their length is relevant. I am open to suggestions on what to do here.