In semiconductor heterostructures it is important to treat the case of position-dependent effective mass. The (mostly) accepted form of kinetic energy operator in this case is:
$$ \hat{T}=- \frac{\hbar^2}{2} \nabla \left( \frac{1}{m^*(\vec{r})} \right) \nabla $$
However, for electron under magnetic field with vector potential $\vec{A}$ and constant effective mass the correct kinetic energy operator is:
$$ \hat{T}=\frac{1}{2m^*} \left(-i \hbar \nabla+\frac{e}{c} \vec{A} \right)^2 $$
This operator is Hermitian, since matrix elements of the imaginary part are antisymmetric (at least in the bases I checked), because it takes a gradient of a wave function:
$$ -i \frac{\hbar e}{m^* c} \vec{A}\nabla \Psi $$
Now how do we treat the general case with both magnetic field and PDEM? It is usually done this way:
$$ \hat{T}=\frac{1}{2}\left(-i \hbar \nabla+\frac{e}{c} \vec{A} \right) \frac{1}{m^*(\vec{r})} \left(-i \hbar \nabla+\frac{e}{c} \vec{A} \right) $$
Or at least in most papers I've seen (such as here and here).
But this operator in general is not Hermitian.
Second update
Let me clarify:
I apply the above operator to a wavefunction. I omit the constants.
$$ \left(-i \nabla+\vec{A} \right) \frac{1}{m} \left(-i \nabla+\vec{A} \right) \Psi= $$
$$ =\left(-i \nabla+\vec{A} \right) \left(-i \frac{1}{m} \nabla \Psi+ \frac{1}{m} \vec{A} \Psi \right)= $$
$$ =-\nabla \left( \frac{1}{m} \nabla \Psi \right)-i \nabla \left( \frac{1}{m} \vec{A} \Psi \right)-i \frac{1}{m} \vec{A} \nabla \Psi+\frac{1}{m} A^2 \Psi= $$
$$ = \left[- \frac{1}{m} \Delta \Psi - \left( \nabla \frac{1}{m} \right) \left( \nabla\Psi \right) +\frac{1}{m} A^2 \Psi \right] - i \left[ 2 \frac{1}{m} \vec{A} \nabla \Psi + \vec{A} \left( \nabla \frac{1}{m} \right) \Psi + \frac{1}{m} \left( \nabla \vec{A} \right) \Psi\right] $$
Now the first two terms are identical to the case with no magnetic field, and they are Hermitian together. Third term is obviously Hermitian. The fourth term (the first term of the imaginary part) is also Hermitian.
I have the problem with last two terms:
$$ \hat{X}=-i \left( \vec{A} \left[ \nabla \frac{1}{m} \right]+\frac{1}{m} \left[ \nabla\vec{A} \right] \right) \Psi $$
The wavefunction is just multiplied by this operator. Meaning if we take:
$$ <\Psi_l | \hat{X}| \Psi_n> $$
It will be symmetric under $l \leftrightarrow n$ exchange, so the resulting matrix will not be Hermitian, unless:
$$ \vec{A}(\vec{r}) \left[ \nabla \frac{1}{m^*(\vec{r})} \right]+\frac{1}{m^*(\vec{r})} \left[ \nabla\vec{A}(\vec{r})\right]=0 $$
Note, that $\Psi$ is real in the bases I use.
Edit 3
Now I finally get it. I can use integration by parts this way:
$$ i \int^{+\infty}_{-\infty} \Psi_l \frac{d}{dx} \left( \frac{1}{m} A_x \Psi_n \right) dx =i \int^{+\infty}_{-\infty} \Psi_l d \left( \frac{1}{m} A_x \Psi_n \right) = $$
Integrating by parts and discarding boundary values (since they are zero anyway in the case relevant to me):
$$ =-i \int^{+\infty}_{-\infty} \Psi_n \frac{1}{m} A_x d\Psi_l = - i \int^{+\infty}_{-\infty} \Psi_n \frac{1}{m} A_x \frac{d\Psi_l}{dx} dx $$
Question is answered thanks to @ValterMoretti