Hamiltonian with ext. vector potential - complex kinetic energy in a given (TD)DFT code with an atomic basis set, i.e.
$$
\psi(\mathbf{r},t) = \sum_i c_i(t) \phi_i(\mathbf{r})
$$
(where the non-on-site basis functions $\phi_i$ aren't necessarily orthogonal), the kinetic energy matrix is calculated via
$$
T_{ij}=\langle \phi_i|\hat{\mathbf{T}}|\phi_j\rangle \sim\langle\phi_i|\hat{\mathbf{p}}^2|\phi_j\rangle \sim -\langle\phi_i|\nabla^2|\phi_j\rangle.
$$
Now in the case where an external time-dependent and spacial homogeneous vector potential $\mathbf{A}(t)$ is considered, the momentum operator reads
$$
\hat{\mathbf{p}}'=\hat{\mathbf{p}}+e\mathbf{A}
$$
and the kinetic energy matrix elements will be like
$$
T_{ij}'\sim\langle\phi_i|-\nabla^2-2ie\mathbf{A}\cdot\nabla+e^2\mathbf{A}^2|\phi_j\rangle
$$
since $\nabla\cdot\mathbf{A}=0$ is used. 
Now my problem: obviously, there is a complex contribution to the kinetic energy that won't vanish. How can this be resolved, what am I missing?
Thanks and best regards!
 A: The operator $\vec{p} \vec{A}+\vec{A} \vec{p}$ is Hermitian which means that its eigenvalues are real. So if you evaluate the matrix element between identical states (it corresponds to the mean value of the observable) which are real in some basis, you will just get $0$. The same story as with the operator of momentum.
For the general matrix element not only 2nd term may be complex, but also the 1st and 3rd ones (just sandwich the Laplacian between two arbitrary complex functions and you'll see it).
A: Thanks for your answers. I do understand that, formally, the expectation value of a hermitian operator must be real. This is more about the implementation I think. 
The basis functions $\phi_i$ (not eigenfunctions of the hamiltonian) are given real functions and the program gives you $\phi_i(\mathbf{r})$, $\nabla \phi_i(\mathbf{r})$ and so on. Then, the matrix elements are calculated via a numerical integration as defined above:
$$
T_{ij}'\sim-\int_V\phi_i(\mathbf{r})[\nabla^2\phi_j](\mathbf{r})-2ie\mathbf{A}\cdot\int_V\phi_i(\mathbf{r})[\nabla\phi_j](\mathbf{r})+e^2\mathbf{A}^2\int_V\phi_i(\mathbf{r})\phi_j(\mathbf{r})
$$
So this is a complex number. Maybe I am missing the point in the sense, that, finally, the solution to the generalized eigenvalue problem,
$$
\left(\mathbf{H}-E\mathbf{S}\right)\mathbf{c}=0
$$
must deliver real coefficients $c_{ij}$ while $\mathbf{H}$ can be complex (thus, the eigensolver has to be chosen appropriately)?
EDIT: apparently, in this formulation, the resulting hamiltonian isn't hermitian which is another problem...
