You are interested in finite Lie algebras, which is the same thing as finite-dimensional Lie algebras over finite fields. The usual way to understand such a problem is first to consider finite-dimensional Lie algebras algebraically closed fields of characteristic $p$.
So let's assume the ground field $k$ is algebraically closed.
Hogeweij has determined here the automorphism group of the following Lie algebras:
${\mathfrak g}={\rm Lie}(G)$ where $G$ is a simple algebraic group;
${\mathfrak f}={\mathfrak g}/{\mathfrak z}({\mathfrak g})$ where ${\mathfrak g}={\rm Lie}(G)$ for a simply-connected simple algebraic group $G$.
Slightly confusingly, he called ${\mathfrak f}$ the classical Lie algebra of type $X_l$ (where $X_l$ is the root system). (Hogeweij remarked that the automorphism groups of most of the classical Lie algebras were earlier determined by Steinberg.) He used the notation: $\widetilde{G}$, $\widetilde{\mathfrak g}$ for the simply-connected case (called universal by Hogeweij) and $\bar{G}$, $\bar{\mathfrak g}$ for the adjoint case. There are no other possibilities except in type $A_l$ with $(l+1)$ composite or type $D_l$ with $p=2$. Note that ${\mathfrak z}({\mathfrak g})$ is trivial and $\widetilde{\mathfrak g}\cong\bar{\mathfrak g}$ except:
when $p$ divides $(l+1)$ in type $A_l$;
when $p=2$ in types $B_l$, $C_l$, $D_l$, $E_7$;
when $p=3$ in type $E_6$.
The appearance of $G_2$ from your calculations is (surely) related to the following fact, well-known to specialists in positive characteristic Lie algebras: in characteristic 3 (and only in characteristic 3) a Lie algebra of type $G_2$ has an ideal, generated by the short root elements, isomorphic to $\mathfrak{psl}(3,k)$. (Note that this is a 7-dimensional Lie algebra which is not the Lie algebra of an algebraic group.) The simple algebraic group of type $G_2$ therefore acts as automorphisms of $\mathfrak{psl}(3,k)$. However, I don't understand quite how it arises in your set-up.
Indeed, Hogeweij's tables show that (in characteristic 3) the automorphism group of $\mathfrak{psl}(3,k)$ is simple of type $G_2$, but both $\mathfrak{sl}(3,k)$ and $\mathfrak{pgl}(3,k)$ have automorphism group equal to a semidirect product of ${\rm PGL}(3,k)$ and the group of diagram automorphisms (which is cyclic of order 2). So unless ${\rm PGL}(3,\overline{{\mathbb F}_3})$ somehow contains something like $G_2({\mathbb F}_3)$ (which seems unlikely to me, though I am not certain) then I can't see how your result can be correct.
This brings me to my second point. Your algorithm is doomed to fail for large $p$, since the exponential of ${\rm ad}\, x$ is only defined if $({\rm ad}\, x)^p = 0$. In general this holds for some but not all nilpotent elements - it holds for all nilpotent elements if and only if $p$ is greater than or equal to the Coxeter number. The usual construction of the Chevalley groups (if I remember correctly) uses only $\exp(\lambda\,{\rm ad}\, e_\alpha)$ where $e_\alpha$ is a simple (positive or negative) root element. The initial choice of a ${\mathbb Z}$-form ensures that these exponentials are well-defined. But I don't see any reason to expect that the group generated by the exponentials of these simple root elements will be equal to the group generated by the exponentials of all nilpotent ($p$-th power zero) elements. This probably (at least in part) explains the discrepancy between your group of order 9285337152 and the group you expected to obtain.