  The present paper studies concentration phenomena of semiclassical
approximation of a massive Dirac equation with general nonlinear self-coupling:
\[ -i\hbar\alpha\cdot\nabla w+a\beta w+V(x)w=g(|w|)w \,. \] Compared with some
existing issues, the most interesting results obtained here are twofold: the
solutions concentrating around local minima of the external potential; and the
nonlinearities assumed to be either super-linear or asymptotically linear at
the infinity. As a consequence one sees that, if there are $k$ bounded domains
$\Lambda_j\subset\mathbb{R}^3$ such that $-a<\min_{\Lambda_j}
V=V(x_j)<\min_{\partial\Lambda_j}V$, $x_j\in\Lambda_j$, then the $k$-families
of solutions $w_\hbar^j$ concentrates around $x_j$ as $\hbar\to 0$,
respectively. The proof relies on variational arguments: the solutions are
found as critical points of an energy functional. The Dirac operator has a
continuous spectrum which is not bounded from below and above, hence the energy
functional is strongly indefinite. A penalization technique is developed here
to obtain the desired solutions.
