Is there a way to make the identity of a gradient of a product of matrix and vector, similar to divergence identity, that would go something like this:
$\nabla(\textbf{M}.\textbf{c})= \nabla(\textbf{M}).\textbf{c}\ +\ ... (\text{not necessarily like this}) $,
where M is a $n\times n$ matrix and $c$ is a $n\times 1$ matrix (column vector)?
I think you mean chain rule rather than divergence identity (since this is what it seems from your equation). It is convenient to use Einstein's notation in this case:
$ \nabla ( \textbf{M} \cdot \textbf{c}) = \dfrac{\partial}{\partial x_k}(M_{ij}c_j) = \dfrac{\partial M_{ij}}{\partial x_k} c_j + \dfrac{\partial c_{j}}{\partial x_k} M_{ij} = \nabla ( \textbf{M} ) \textbf{c} + \nabla ( \textbf{c} ) \textbf{M}$
where the gradient of a matrix M is defined as a tensor of order three: $ \dfrac{\partial M_{ij}}{\partial x_k} $