This is best understood by approximating the dipole as a pair of finite charges $\pm q$ separated by a finite distance $d$. In a uniform electric field, the electrostatic forces on each of the charges will cancel out exactly, but in a non-uniform one the forces on the two will be slightly different, leading to a slight imbalance and therefore a non-zero net force. As you take the distance to zero, the difference in electric field goes to zero, but the charge also grows to exactly cancel it out.
To be more quantitative, suppose the negative charge is at $\mathbf r$ and the positive charge at $\mathbf r+d\mathbf n$. The total force is then
$$
\mathbf F=q\left[\mathbf E(\mathbf r+d\mathbf n)-\mathbf E(\mathbf r)\right].
$$
To get the correct form for the limit, change from the charge $q$ to the electric dipole $p=qd$, to get
$$
\mathbf F=p\frac{\mathbf E(\mathbf r+d\mathbf n)-\mathbf E(\mathbf r)}{d}.
$$
The true force on a point dipole is the limit of this as $d\to0$,
$$
\mathbf F=p\lim_{d\to0}\frac{\mathbf E(\mathbf r+d\mathbf n)-\mathbf E(\mathbf r)}{d},
$$
and this is exactly the directional derivative along $\mathbf n$, typically denoted $\mathbf n\cdot \nabla$, so
$$
\mathbf F=p\mathbf n\cdot \nabla\mathbf E=\mathbf p\cdot \nabla\mathbf E.
$$
The torque $ \tau $ on an electric dipole with dipole moment p in a uniform electric field E is given by $$ \tau = p \times E $$ where the "X" refers to the vector cross product.
Ref: Wikipedia article on electric dipole moment.
I will demonstrate that the torque on an ideal (point) dipole on a non-uniform field is given by the same expression.
I use bold to denote vectors.
Let us begin with an electric dipole of finite dimension, calculate the torque and then finally let the charge separation d go to zero with the product of charge q and d being constant.
We take the origin of the coordinate system to be the midpoint of the dipole, equidistant from each charge. The position of the positive charge is denoted by $\mathbf r_+ $ and the associated electric field and force by $\mathbf E_+$ and $ \mathbf F_+$, respectively. The notation for these same quantities for the negative charge are similarly denoted with a - sign replacing the + sign.
The torque about the midpoint of the dipole from the positive charge is given by
$$ \mathbf \tau_+ = \mathbf r_+ \times \mathbf F_+ $$
where
$$ \mathbf F_+ = q\mathbf r_+ \times \mathbf E_+(\mathbf r+) $$
Similarly for the negative charge contribution
$$ \mathbf \tau_- = \mathbf r_- \times \mathbf F_- $$
where
$$ \mathbf F_- = -q\mathbf r_- \times \mathbf E_-(\mathbf r-) $$
Note that
$$ \mathbf r_- = -\mathbf r_+ $$
We can now write the total torque as
$$ \mathbf \tau_{tot} = \mathbf \tau_- + \mathbf \tau_+ =q\mathbf r_+ \times (\mathbf E(\mathbf r_+)+\mathbf E(\mathbf r_-))$$
It is clear that in taking the limit as the charge separation d goes to zero, the sum of electric fields will only contain terms of even order in d.
Noting that $$ \mathbf |r_+| = \frac{d}{2} $$
and defining in the usual way $$ \mathbf p = q\mathbf d = q(\mathbf r_+ - \mathbf r_- ) $$
We can write that $$ \tau_{tot} = \mathbf p \times \mathbf E(0) + \ second \ order \ in \ d $$
As we take the limit in which d goes to zero and the product qd is constant, the second order term vanishes.
Thus, for an ideal (point) dipole in a non-uniform electric field, the torque is given by the same formula as that of a uniform field.
Note that it is not correct to start with the expression for a force on an ideal/point dipole in a non-uniform field and then calculate torque from this force. To derive this expression one ends up first taking the limit of a point dipole (on which there is zero force in a uniform field) and then one finds a torque of zero, which is incorrect. One must start with the case of a finite dipole, calculate torque and only then pass to the limit.
When p and E are parallel and anti-parallel, the torque is zero, so yes zero is possible. But the case in which p and E are anti-parallel is one of an unstable equilibrium, and a small angular perturbation will cause the dipole to experience a torque which attempts to align the dipole with the electric field.
Best Answer
Both formulas are equivalent, if you are in the electrostatic approximation and your dipole vector does not depend on the position $\mathbf{r}$.
Let's consider the expression $\mathbf{F}=\nabla_{\mathbf{r}}(\mathbf{p} \cdot \mathbf{E})$ which can be easily obtained from the potential energy function
$U=-\mathbf{p} \cdot \mathbf{E}$
and its relation with the force $\mathbf{F}=\nabla_\mathbf{r} U$. Now, recall the vector identity
$\nabla_\mathbf{r}(\mathbf{a}\cdot \mathbf{b})= (\mathbf{a} \cdot \nabla_\mathbf{r}) \mathbf{b}+(\mathbf{b} \cdot \nabla_\mathbf{r}) \mathbf{a} + \mathbf{a} \times (\nabla_\mathbf{r} \times \mathbf{b})+ \mathbf{b} \times (\nabla_\mathbf{r} \times \mathbf{a})$
for $\mathbf{a}=\mathbf{a}(\mathbf{r})$ and $\mathbf{b}=\mathbf{b}(\mathbf{r})$ two arbitrary vectors. For $\mathbf{p}=\mathbf{a} \neq \mathbf{p}(\mathbf{r})$ [independent of the position] and $\mathbf{b}=\mathbf{E}(\mathbf{r}$) we have
$\nabla_\mathbf{r}(\mathbf{p}\cdot \mathbf{E})= (\mathbf{p} \cdot \nabla_\mathbf{r}) \mathbf{E}+(\mathbf{E} \cdot \nabla_\mathbf{r}) \mathbf{p} + \mathbf{p} \times (\nabla_\mathbf{r} \times \mathbf{E})+ \mathbf{E} \times (\nabla_\mathbf{r} \times \mathbf{p})$
As the dipole vector does not depend on the position we can drop the second and the fourth terms. In the electrostatic approximation, Faraday's law reads $\partial_t \mathbf{B}=\mathbf{0}\Leftrightarrow \nabla_\mathbf{r} \times \mathbf{E}(\mathbf{r})=\mathbf{0} $ [this is known as ''Carn's law''] so that the electric field is irrotational and the curl vanishes. Then we can drop the third term and
$\nabla_\mathbf{r}(\mathbf{p}\cdot \mathbf{E})= (\mathbf{p} \cdot \nabla_\mathbf{r}) \mathbf{E}$
so that your definitions agree.