Given a function $G\left(c_{1}, \ldots, c_{N}\right)$ \begin{aligned} G\left(c_{1}, \ldots, c_{N}\right) &=\sum_{i j} c_{i}^{*} c_{j} H_{i j}-\epsilon \sum_{i j} c_{i}^{*} c_{j} \delta_{i j} \\ &=\sum_{i j} c_{i}^{*} c_{j}\left(H_{i j}-\epsilon \delta_{i j}\right) \end{aligned}
Im not so sure how $$\frac{\partial G}{\partial c_{i}}=\sum_{j}\left(H_{i j}-\epsilon \delta_{i j}\right) c_{j}.$$ There is some explanation in this file (https://iramis.cea.fr/spcsi/cbarreteau/physique_du_solide/documents/variational_method.pdf) but an explicit derivation would greatly help me to understand better. In particularly, I'm puzzled by the absence of imaginary coefficients.
you can see this derivation as a "standard" partial derivative of G with respect to the variable $c_i^*$ An important fact to note is that when you do this derivation you should consider $c_i^* $ and $c_i$ as independent parameters, leading to tour remark (absence of imaginary complex conjugate coefficients)
The key formula to use is then
$\frac{\partial c_j}{\partial c_i^*}=0$ and $\frac{\partial c_j^*}{\partial c_i^*}=\delta_{ij}$