Estimation of the gradient of the logarithm of a probability density function is a versatile tool in statistical data analysis. A recent method for model-seeking clustering called the least-squares log-density gradient clustering (LSLDGC) [Sasaki et al., 2014] employs a sophisticated gradient estimator, which directly estimates the log-density gradients without going through density estimation. However, the typical implementation of LSLDGC is based on a spherical Gaussian function, which may not work well when the probability density function for data has highly correlated local structures. To cope with this problem, we propose a new gradient estimator for log-density gradients with Gaussian mixture models (GMMs). Covariance matrices in GMMs enable the new estimator to capture the highly correlated structures. Through the application of the new gradient estimator to mode-seeking clustering and hierarchical clustering, we experimentally demonstrate the usefulness of our clustering methods over existing methods.