Symmetric part of hamilton matrix and eigenvalues I'm working with a Density Functional Theory package and I am not sure with a procedure in the code. At some point, after the Hamilton matrix $H_{ij}=\langle\phi_i|\hat{H}|\phi_j\rangle\in\mathbb{R}$ (where $\{\phi_i\}$ is the basis set used) is calculated, the following step is done before the eigenvalues and -vectors of $\mathbf{H}$ are calculated:
$$
\mathbf{H} := \frac{1}{2}\left(\mathbf{H} + \mathbf{H}^T\right),
$$
i.e. only the symmetric part of the Hamilton matrix is used for further computation. Since $\mathbf{H}$ should be symmetric (or self-adjoint in general, but here $H_{ij}$ are purely real) by definition, I don't understand why this 'symmetrization' is done explicitly. Don't we lose information here? What about the eigenvalue spectrum?
 A: If you want to describe inversion symmetric systems the Hamiltonian matrix is a symmetric matrix and otherwise it is a Hermitian matrix. This is an important property to obtain real-valued energy eigenvalues.
Depending on the approach to density functional theory this property of the Hamiltonian matrix is not mathematically guaranteed in the first place but has to be ensured by postprocessing measures like the one you describe in the question. 
Let us assume that you are using a (linearized) augmented-plane-wave basis set as it is used in software packages like Wien2k, Fleur, Elk, Exciting and others. With such a representation the unit cell is separated into non-overlapping but nearly touching so-called muffin-tin (MT) spheres around each atom $\alpha$ and an interstitial region (IR) in between. In the IR such a basis function is a plane wave with wave vector $\mathbf{k} + \mathbf{G}$, where $\mathbf{k}$ is a Bloch vector and $\mathbf{G}$ a reciprocal lattice vector. In the MT spheres each basis function consists of radial functions times spherical harmonics up to a certain angular momentum cutoff parameter $l_\text{max}$. The radial functions are scaled to ensure continuity in value (and slope) at the MT spheres boundaries. But this matching of course only happens for $l \le l_\text{max}$. For higher $l$ the basis functions feature discontinuities at the MT sphere boundaries.
Let us test whether we obtain a Hermitian matrix with such a basis. Hermiticity means
\begin{equation}
H_{\mathbf{G}\mathbf{G}'}^{\mathbf{k}} - \left(H_{\mathbf{G}'\mathbf{G}}^{\mathbf{k}}\right)^\ast = 0.
\end{equation}
With the sketched basis consisting of IR representations $\phi_{\mathbf{k}\mathbf{G}}^\text{IR}$ and MT representations $\phi_{\mathbf{k}\mathbf{G}}^\alpha$ this expression becomes
\begin{align}
H_{\mathbf{G}\mathbf{G}'}^{\mathbf{k}} - \left(H_{\mathbf{G}'\mathbf{G}}^{\mathbf{k}}\right)^\ast & = \left\langle \phi_{\mathbf{k}\mathbf{G}}^\text{IR} \middle\vert \hat{H} \middle\vert \phi_{\mathbf{k}\mathbf{G}'}^\text{IR} \right\rangle_\text{IR} - \left(\left\langle \phi_{\mathbf{k}\mathbf{G}'}^\text{IR} \middle\vert \hat{H} \middle\vert \phi_{\mathbf{k}\mathbf{G}}^\text{IR} \right\rangle_\text{IR} \right)^\ast \nonumber \\ & \phantom{=} + \sum\limits_\alpha \left\langle \phi_{\mathbf{k}\mathbf{G}}^\alpha \middle\vert \hat{H} \middle\vert \phi_{\mathbf{k}\mathbf{G}'}^\alpha \right\rangle_{\text{MT}^\alpha} - \left(\left\langle \phi_{\mathbf{k}\mathbf{G}'}^\alpha \middle\vert \hat{H} \middle\vert \phi_{\mathbf{k}\mathbf{G}}^\alpha \right\rangle_{\text{MT}^\alpha} \right)^\ast \nonumber \\ & = -\frac{1}{2} \left\lbrace\int\limits_\text{IR} \left(\phi_{\mathbf{k}\mathbf{G}}^\text{IR}(\mathbf{r})\right)^\ast \nabla^2 \phi_{\mathbf{k}\mathbf{G}'}^\text{IR}(\mathbf{r}) - \phi_{\mathbf{k}\mathbf{G}'}^\text{IR}(\mathbf{r}) \nabla^2 \left(\phi_{\mathbf{k}\mathbf{G}}^\text{IR}(\mathbf{r})\right)^\ast d^3r \right. \nonumber \\ & \phantom{=} + \left. \sum\limits_\alpha \int\limits_{\text{MT}^\alpha} \left(\phi_{\mathbf{k}\mathbf{G}}^\alpha(\mathbf{r})\right)^\ast \nabla^2 \phi_{\mathbf{k}\mathbf{G}'}^\alpha(\mathbf{r}) - \phi_{\mathbf{k}\mathbf{G}'}^\alpha(\mathbf{r}) \nabla^2 \left(\phi_{\mathbf{k}\mathbf{G}}^\alpha(\mathbf{r})\right)^\ast d^3r \right\rbrace,
\end{align}
where we already used the fact that the contributions related to the potential cancel each other. With Greens's second identity we can turn these volume integrals into surface integrals over the MT sphere boundaries. We obtain
\begin{align}
H_{\mathbf{G}\mathbf{G}'}^{\mathbf{k}} - \left(H_{\mathbf{G}'\mathbf{G}}^{\mathbf{k}}\right)^\ast & = -\frac{1}{2}\left\lbrace \sum\limits_\alpha \int\limits_{\partial \text{MT}^\alpha} -\left[\left(\phi_{\mathbf{k}\mathbf{G}}^\text{IR}(\mathbf{r})\right)^\ast \nabla \phi_{\mathbf{k}\mathbf{G}'}^\text{IR}(\mathbf{r})\right] \mathbf{n} + \left[\phi_{\mathbf{k}\mathbf{G}'}^\text{IR}(\mathbf{r}) \nabla \left(\phi_{\mathbf{k}\mathbf{G}}^\text{IR}(\mathbf{r})\right)^\ast\right] \mathbf{n} \right. \nonumber \\ & \phantom{=} \left. \vphantom{\sum\limits_\alpha \int\limits_{\partial \text{MT}^\alpha}} + \left[\left(\phi_{\mathbf{k}\mathbf{G}}^\alpha(\mathbf{r})\right)^\ast \nabla \phi_{\mathbf{k}\mathbf{G}'}^\alpha(\mathbf{r})\right] \mathbf{n} - \left[\phi_{\mathbf{k}\mathbf{G}'}^\alpha(\mathbf{r}) \nabla \left(\phi_{\mathbf{k}\mathbf{G}}^\alpha(\mathbf{r})\right)^\ast\right] \mathbf{n} dS \right\rbrace.
\end{align}
This expression clearly becomes $0$ if the basis functions are continuous in value and slope. But as mentioned above this is only the case for $l \le l_\text{max}$. The cutoff parameter therefore forces us to perform postprocessing steps like the one sketched in the question. 
The final results should not depend on such tricks as one typically increases the cutoff parameters until the results are converged. Of course, other approaches to DFT may have similar reasons to perform such measures.
