I only read until the Gâteaux derivative, that I would do this way:
func[p_] :=
Inactive[Sum][
Inactive[D][p[k],(*function*)
Inactive[Part][x, k]],(*\[PartialD]/\[PartialD]x_k kept inert*){k,
1, 3}];
MapAt[D[#, \[Epsilon]] &,
func[Function[i,
p[i] @@ x + \[Epsilon] \[Psi] @@ x KroneckerDelta[k, 1]]],
1]
