One way to see “Dieudonné‐modules are really cohomology” is to observe that Messing’s Dieudonné‐crystal is by definition an Ext–group in the crystalline topos of the base. I will sketch two (equivalent) points of view:

1. The Ext–definition (Messing, BBM).
   Let S be a p‐adic base (say S=Spec k with k perfect of char p, or more generally any p‐nilpotent formal scheme) and let X→S be a p‐divisible group. On the crystalline site (S/W){cris}, X defines an fppf‐sheaf of groups X{cris}. Messing (and later Berthelot–Breen–Messing) show that the covariant Dieudonné‐crystal
   D(X)
   is naturally D(X)(T,J,γ) ≔ Ext^1_{(T/W){cris},,fppf}!\bigl(X_T,,,O{T/W}\bigr),,
   for every PD‐thickening (T,J,γ) of S. One checks that as (T,J,γ) varies this is indeed a crystal in finite locally free O_{T/W}–modules, endowed with the usual Frobenius and Verschiebung. Finally, for S=Spec k with its canonical W(k)–thickening one sets M(X) := D(X)\bigl(,W(k)\bigr), recovering the usual Dieudonné module of X.

   In particular M(X) is literally an Ext^1–group in the crystalline topos of the base. Modulo p, the same construction on the infinitesimal site gives
   D(X) mod p ≅ Ext^1_{(S/k){inf}}(X, O) ≅ M(X)/p,
   which in the abelian‐variety case recovers H^1{dR}(A/k).

2. Finite‐level approximation (BBM).
   If X[p^n] is the finite flat n–torsion of the p-divisible group, one shows (BBM, Thm. 3.3.6) that for each n H^i_{!cris}\bigl(X[p^n]/W(k)\bigr) vanishes for i≠0,1, and that H^1_{!cris}\bigl(X[p^n]/W(k)\bigr) is canonically the Dieudonné–module of X[p^n]. Passing to the inverse limit in n one obtains M(X) ≅ \varprojlim_n H^1_{!cris}\bigl(X[p^n]/W(k)\bigr),. When X=A[p^∞] this simply becomes H^1_{!cris}(A/W(k)), recovering the Mazur–Messing–Oda theorem.

In neither approach does one really invent a brand‐new “site of p-divisible groups over X” – rather, one uses the usual crystalline (or infinitesimal) site of the base (or the finite torsion‐schemes X[p^n]), and regards X itself as a sheaf there. The key point is:

• Dieudonné–crystals = Ext^1(X, O) in the crystalline topos of the base.
• M(X) is its evaluation on the canonical PD‐thickening (e.g. W(k)),
and in finite‐level language one really has
M(X) ≅ limₙ H^1_{cris}(X[pⁿ]/W(k)).

References.
– M. Messing, “The crystals associated to Barsotti–Tate groups: with applications to abelian schemes,” Lecture Notes in Math. 264 (1972).
– P. Berthelot, L. Breen, W. Messing, “Théorie de Dieudonné cristalline,” Lecture Notes in Math. 930 (1982).
– T. T. Liu, “Six operations and Lefschetz–Verdier formula for arithmetic D-modules,” especially the discussion of crystalline Ext.

So in this precise sense the Dieudonné module is nothing mysterious: it is a first Ext–group (or a first crystalline cohomology) of the Barsotti–Tate group itself.