$$\begin{array}{rcl} \sum_{i=1}^n (y_i-\hat{y_i})(\hat{y_i} - \bar{y}) = \sum_{i=1}^n \epsilon_i(\hat{y_i} - \bar{y})
= \underbrace{\sum_{i=1}^n \epsilon_i \hat{y_i}}_{\text{ equals zero if } \epsilon \perp \hat{y} } &-& \bar{y_i}\underbrace{ \sum_{i=1}^n \epsilon_i}_{ \text{equals zero} }
\end{array}$$
I do not believe that you have $\epsilon \perp \hat{y}$ for GLM, this is only the case if you have the residual vector represent the shortest distance between $y$ and the solution space of all possible $\hat{y}$, which is the least squares solution and relates to the assumption of Gaussian distributed errors.
Neither is $ \sum_{i=1}^n \epsilon_i = 0 $ general. This means $\epsilon$ is perpendicular to the intercept term. Even in the case of a simple linear model this may not be true (simple linear regression without the intercept term).
(Computational) counter-examples:
> # some data
> x <- c(1.0, 2.0, 3.0, 4.0)
> y <- c(1.1, 2.1, 3.2, 4.0)
>
> # OLS (minimize squared errors)
> hat_y <- predict(lm(y~x))
> sum((y-hat_y)*(hat_y-mean(y)))
[1] 6.366435e-16
>
> # GLM (using Gamma)
> hat_y <- predict(glm(y~x, family = Gamma(link="identity")),type="response")
> sum((y-hat_y)*(hat_y-mean(y)))
[1] -0.1042444
>
> # GLM (using Gaussian but different link)
> hat_y <- predict(glm(y~x, family = gaussian(link="log")),type="response")
> sum((y-hat_y)*(hat_y-mean(y)))
[1] 0.1504989
>
> # OLS without intercept term
> hat_y <- predict(lm(y~0+x))
> sum((y-hat_y)*(hat_y-mean(y)))
[1] -0.26