How to compute the formal group law of a Shimura variety (using its invariant differentials)?
I don't know how to do this purely with invariant differentials. One of the main reasons is that, in order to solve this problem, you need to construct an "algebraic" lift of this variety from the complex numbers to a much smaller number field. This lift isn't necessarily unique and different choices will give you different formal group laws (though they'll become isomorphic over extension fields). What I'll describe, then, is a method for determining -- as best I know how -- an algebraic lift of $A$ and the $p$-series in terms of a coordinary on the formal group law.
Let's pick a finite field $\Bbb F_q$ where $q = p^k$ and discuss things in terms of the Tate-Honda classification of abelian varieties over $\Bbb F_q$.
Your description in terms of the Shimura-Taniyama formula means that our lift of $A$ is going to give us an abelian variety over $\Bbb F_q$ which is $3$-dimensional and whose Newton polygon breaks up into a component of slope $1/3$ and a component of slope $2/3$. This forces $A$ to be simple because the $p$-divisible group can't break up into simpler pieces that still satisfy self-duality.
The Frobenius endomorphism $\pi$ of $A$ is an element of the ring of integers ${\cal O}_K$, and it's a Weil $q$-integer: it has absolute value $q^{1/2}$ in any complex embedding. Moreover, the Shimura-Taniyama formula tells you exactly which primes of ${\cal O}_K$ divide $\pi$ and how much. The only indeterminacy left is possibly multiplying $\pi$ by a unit $u$ which has absolute value $1$ in any complex embedding, which forces $u$ to be a root of unity (which we've got a goodly stock of here). Replacing $\pi$ with $u \pi$ replaces your abelian variety by a Galois twist: extending fields from $\Bbb F_q$ to $\Bbb F_{q^k}$ replaces $\pi$ by $\pi^k$, and so this root-of-unity indeterminacy goes away after a sufficient extension of finite fields.
We should also note for a moment that while, in theory, the (Galois conjugacy class of) Frobenius $\pi$ only determines $A$ up to isogeny, the constraint that $A$ should have an action of ${\cal O}_K$ actually means that $\pi$ determines $A$ up to isomorphism. This is because any isogenous abelian variety that still has ${\cal O}_K$-multiplication is determined by an ${\cal O}_K$-stable lattice in $\Bbb Q \otimes {\cal O}_K$ -- i.e., a fractional ideal -- and because $K$ has class number one, any such fractional ideal is principal. Multiplication by a generator of the ideal gives an isomorphism between the isogenous variety and $A$.
This means that, up to a Galois twist, we can determine $A$ completely over a finite field just from this slope data. We can't get any farther without making some choices, so let's make some. I'll say our prime is $p=2$ for the purposes of illustration.
Fix $\zeta \in K$ to be our chosen seventh root of unity, and let $w = \zeta + \zeta^2 + \zeta^4$. Then $w^2 + w + 2 = 0$ so $w = \tfrac{-1\pm \sqrt{-7}}{2}$. The element $w$ has absolute value $2^{1/2}$ in any complex embedding, and $w \bar w = 2$; in $K$, this reflects the splitting of prime ideals $(2) = (w) (\bar w)$.
Let $\pi = 2 \bar w$. Then $\pi$ has minimal polynomial $t^2 + 2t + 8$ and has absolute value $8^{1/2}$ in any complex embedding. It's also divisible by $w$ once and by $\bar w$ twice. This means that $\pi$ is a Weil $8$-integer and an abelian variety associated to $\pi$ over $\Bbb F_8$ has an endomorphism ring $End(B) \otimes \Bbb Q$ a central division algebra over $\Bbb Q(\sqrt{-7})$ with invariant $1/3$ at $w$ and $2/3$ at $\bar w$, with trivial invariant at all other places. The field $K$ embeds into this one, so $B$ has complex multiplication by $K$. By replacing $B$ with an isogenous abelian variety $A$, we can get an abelian variety $A$ over $\Bbb F_8$ with an action of ${\cal O}_K$ whose Frobenius (the $8$'th power Frobenius, since we're working over $\Bbb F_8$) satisfies the minimal polynomial $F^2 + 2F + 8$. This abelian variety, having the correct slope data, will be an algebraic lift of the abelian variety you constructed over the complex numbers because of the argument we sketched above.
Now let's look at this 1-dimensional summand of the formal group law, which corresponds to a $1$-dimensional summand of the $2$-divisible group. If we look at the $2$-divisible group, the Frobenius acts on it too, and it still satisfies the same minimal polynomial -- but now over $\Bbb Z_2$ instead of integrally. The polynomial $t^2 + 2t + 8$ has two solutions over $\Bbb Z_2$ (which we might call $-1 \pm \sqrt{-7}$): one divisible by $2$ and the other divisible by $4$. Let's call the first root $\alpha$ and the other $8/\alpha$. Our minimal polynomial factors now as $F^2 + 2F + 8 = (F - \alpha)(F - 8/\alpha)$, and the splitting of the $2$-divisible group into a piece of slope $1/3$ and one of slope $2/3$ corresponds to the factors where $F$ acts by $\alpha$ and by $8/\alpha$ respectively.
This was a bit of a slog to get to this point, but let's reflect on what we've found: we've found our $1$-dimensional piece of the $2$-divisible group, and we've found that the Frobenius acts on it as multiplication by this element $\alpha = -1 - \sqrt{-7}$. Or, equivalently, we get an identity $$ \tfrac{2}{-1-\sqrt{-7}} F = 2 $$ as endomorphisms of this $1$-dimensional formal group law $G$. But if we write that in a power series expansion, we get $$ \left[\tfrac{2}{-1-\sqrt{-7}}\right]_G(x^8) = [2]_G(x) $$ If we want to push this a little farther, we have to start doing a little more gruntwork. Namely, we've got to write this element $\tfrac{2}{-1-\sqrt{-7}}$ in binary expansion as a $2$-adic integer: $$ \tfrac{2}{-1-\sqrt{-7}} = a_0 + 2 a_1 + 4 a_2 + \dots $$ where $a_i \in \{0,1\}$. Then our identity for the Frobenius becomes $$ \sum_G a_i [2]_G^{(i)}(x^8) = [2]_G(x). $$ And then we can start solving for the 2-series iteratively. It starts with $x^8$. Let's write $f(x)$ for $[2]_G(x)$. I believe the $a_i$ starts with $a_0 = 1, a_1 = 0, a_2 = 1, a_3 = 1, a_4 = 0$, etc. So we find saying $$ f(x) = x^8 +_G f(f(x^8)) +_G f(f(f(x^8))) + \dots $$ and so we can start recursively solving to find $$ [2]_G(x) = x^8 + x^{8^3} + \dots $$ I would imagine that there are ways to be a little more clever about this by finding a series expansion for $2$ in terms of $F$ instead of $F$ in terms of $2$. If we can figure out how to write $$ 2 = \sum b_i \alpha^i \equiv \sum b_i F^{(i)} $$ for $b_i \in \{0,1\}$, this becomes $$ [2]_G(x) = \sum_G b_i x^{8^i}, $$ which tells you directly that the Araki generators satisfy $b_i = v_{3i}$ (and the $v_j$ are zero when $j$ is not a multiple of $3$).
While there are a lot of steps in here, the endgame takeaway is not that bad: find an element in ${\cal O}_K$ with appropriate prime divisibility, split it $p$-adically, and then use a $p$-adic power expansion to determine how Frobenius is acting on the appropriate $1$-dimensional summand. Choosing $\pi$ in $\Bbb Q(\sqrt{-7})$ will have a tendency to make your life a little easier.
Hope this helps.