Finding a basis which is orthonormal with respect to symmetric bilinear form

Question

Let $X=\{a\in M_4(\mathbb{C}): a^T=-a\}$ be the $6$-dimensional vector space of matrices. Define the Pfaffian as $$ \mathrm{Pf}(a)=a_{12}a_{34}-a_{13}a_{24}+a_{14}a_{23} $$ There is an associated bilinear form defined by polarisation $B=\frac{1}{4}(\mathrm{Pf}(a+b)-\mathrm{Pf}(a-b))$. Suppose I want to find a basis for $X$ which is orthonormal relative to $B$; that is, a basis $\{a_1,a_2,a_3,a_4,a_5,a_6\}$ such that $B(a_i,a_j)=\begin{cases} \pm1\mbox{ if } i=j \\ 0 \mbox{ if }i\neq j\end{cases}$. Now, it is clear to me what the first three elements of the most obvious basis are: $$ a_1=\begin{pmatrix} 0 & 1 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -1 & 0 \end{pmatrix}, a_2=\begin{pmatrix} 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0\end{pmatrix}, a_3=\begin{pmatrix} 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & -1 & 0 & 0 \\ -1 & 0 & 0 & 0 \end{pmatrix} $$ But I've got no idea how to complete this to the full basis. It is possible to extend this to any basis and then perform a modified version of the Gram-Schmidt process, but is it maybe possible to "read off" the basis elements from the definitions of $\mathrm{Pf}$ and $B$?