While our model is designed to be inexpensive and flexible, aiming at the description of medium to large sized molecular crystals, it can potentially be used to converge towards the many-body limit description of crystals with correlated wave function theory. We are not currently aiming at the same level of accuracy demonstrated in Ref. , for example, but were interested in observing how the quality of the results changes when improving the QM description, and if the method is suited for such types of calculations. There are certain advantages in the current scheme when compared to other many-body expansions. The use of embedding helps reduce the order of the expansion required to reach a certain accuracy level. Furthermore, gradients are readily available.