Dielectric constant of crystals: Feynman's mistake? Feynman mentioned in his lecture that we can, in some way, obtain dielectric constant of crystals by summing up the electric field due to dipoles. For example, equation (11.33) from above link shows the sum for $\rm BaTiO_3$. He also mentioned that the same calculation for simple cubic results to a factor
of 1/3 in simple cubic crystals.
However, if one calculates the exact field due to dipoles and integrate for all the neighbors in simple cubic lattice (e.g. using $(3(p.r)r-p)/r^3$ instead of the one he used $p/r^3$), the result would be zero field due to surrounding dipoles.
I cannot find anyone extending or using this approach to figure out what I am missing. Also equation 15 and 16 in these notes proves the sum $(3(p.r)r-p)/r^3$ is zero. Moreover, in a simple cubic lattice, the simplest example is to consider 6 nearest neighbor dipoles (6 neighbors: top, bottom, left, right, front and back) which gives zero field.
Is this a mistake in Feynman lectures?
To show the sum for 6 neighbors, one can evaluate equation 6 times for 6 neighbors:
        Rs = np.array([[ 0.,  0.,  1.],
                  [ 0.,  1.,  0.],
                  [ 1.,  0.,  0.],
                  [-1.,  0.,  0.],
                  [ 0., -1.,  0.],
                  [ 0.,  0., -1.]])
    pdir = np.array([0, 0, 1.0])
    for ii in range(len(Rs)):
        u = Rs[ii]
        d = 1.0
        Efield_i = (3 * pdir.dot(u) * u - pdir) / (d * d * d)
        print("Efield_i", Efield_i)

Which gives with zero sum:
Efield_i [0. 0. 2.]
Efield_i [ 0.  0. -1.]
Efield_i [ 0.  0. -1.]
Efield_i [-0.  0. -1.]
Efield_i [ 0. -0. -1.]
Efield_i [-0. -0.  2.]

Edit: Please notice this question is not about E-field of dipoles, but about the approach and assumptions Feynman used to calculate dielectric constant.
 A: Since you do not show us your 6 nearest neighbour calculation, I am not able to tell anything about possible flaws in it. However, Maxwell's equations (namely Gauß' law) generally tell us that
$$\vec \nabla\cdot \vec E=4\pi\rho$$
that is, if $\vec E=0$ (which is what you claim to be your result) for a certain distribution, then inevitably
$$\rho = \frac{1}{4\pi}\vec \nabla\cdot \vec E=0$$
which is not a dipole distribution, nor any distribution of any charges, at all. Please compare this to the known charge distribution of a set of dipoles:
$$\rho(\vec r) = \sum_i \vec p_i \cdot \vec\nabla \delta(\vec r-\vec r_i)$$
Also, in your second reference, after equation (14) you can read

Now  we  consider  the  field  that  acts  on  the  atom  at  the  center  of  the  sphere. If  all  dipoles  are parallel  to  the  z  axis  and  have  magnitude p,  the z component  of  the  field  at  the  center  due  to  all other dipoles is...

In other words, what follows is a very special calculation for a very special configuration: not only is a spherical crystal, i.e. charge/dipole distribution assumed, but also the value of the field strength is only calculated at the center of this symmetric distribution. Hence, after equation (15)

The  latter  equation  comes  from  the  fact  that  the x,  y,  z directions  are  equivalent  because  of  the symmetry of the lattice and of the sphere

That is, the field of a set of dipoles does NOT vanish in general (that would be ridiculous), but only in this special setting, which is dedicated to calculating the local field at one atom site, and not the macroscopic electric field (which is described by the dielectric permittivity).
The reason (the Lorentz theory of dielectrics) for this special assumption follows later on after equation (17), where an artifical distinction is made between the rapidly fluctuating near-field of an atom, and the smooth far-field. This distinction is made by "drawing" a more or less arbitrary sphere around the atom in question, which clearly shows the spherical symmetry used earlier.
To reinterate, what vanishes is not the dipole field in general, but only the fluctuating near-field inside an artificial sphere around a specific atom.
The general behavior of dielectrics, as already indicated in my comment, is also described after equation (4):

From the macroscopic point of view the dielectric can be considered as a material with no net charges in the interior of the material and induced negative and positive charges on the left and right surfaces of the dielectric.

This clearly disproves any general statement like "the sum of those polarizations are zero" because then there could not be surface charges.
A: Despite your formula being wrong you get the correct result in this case. But all you did was calculate the electric field at one single point $\mathbf{r} = (0, 0, 0)^T$ which certainly doesn't say much about the field as a whole. The electric field of a single dipole $\mathbf{d}$ at $\mathbf{r}_0$ is given by
$$
\mathbf{E}(\mathbf{r}) = \frac{1}{4\pi\varepsilon_0}\frac{3(\mathbf{r}-\mathbf{r}_0) \mathbf{d}\cdot (\mathbf{r}-\mathbf{r}_0) - \mathbf{d} |\mathbf{r}-\mathbf{r}_0|^2}{|\mathbf{r}-\mathbf{r}_0|^5}
$$
If you redo your code like that an use a variable location at which to evaluate the electric field instead of choosing $(0,0,0)$ you will find a non zero field. It would look like this in the xy plane for example 
In the lecture Feynman is talking about an infinite chain along the z-axis of dipoles in the z direction and says that for a single dipole if you're a distance $r=z$ away from the dipole along the z-axis the field would be
$$
\mathbf{E}(0,0,z) \propto \frac{3 z d_z z \mathbf{e}_z - d_z\mathbf{e}_z z^2}{z^5} = \frac{2z^2 d_z}{z^5}\mathbf{e}_z = \frac{2d}{z^3}\mathbf{e}_z
$$
so the formula is not approximate, it's calculated for this special case. He then sums over the dipoles in this one chain and reasons that the other chains only contribute very little to the electric field.
A: It's quite simple. Consider just two neighbouring dipoles in the $+z$ and $-z$ direction, where $z$ also happens to be the direction of polarisation.
With your formula, they cancel. The field due to a dipole a distance $r$ along its axis of polarisation is equal and opposite to the field a distance $r$ in the opposite direction.
But if you think about it and draw the picture, they clearly add. The field at $(0,0,0)$ due to a $z-$direction dipole at $(0,0,a)$ is the same as the field at $(0,0,0)$ due to a similar dipole at $(0,0,-a)$.
In your formula you need to distinguish carefully between $\vec r$ as the vector from the dipole to the point of interest, and the vector from the point of interest to the dipole. There is an important minus sign between them. If you incorporate this correctly in your algebra or your program you will get a non-zero answer, as you must.
A: Here is first the calculation which I guess the OP did, and then the reason why it is the wrong calculation, and then a better calculation.
First let's put dipoles at $(\pm r,0,0)$,
$(0,\pm r,0)$ and $(0,0,\pm r)$ and consider the field at the origin.
We suppose the dipoles are all aligned in the $x$ direction, each of size $p$.
Let's drop the $1/4\pi \epsilon_0$, so we say the field of each dipole is
$$
{\bf E}_1 = \frac{3 ({\bf p} \cdot \hat{\bf r})\hat{\bf r} - {\bf p}}{r^3}
$$
At the origin the pair on the x axis give a net electric field
$2 \times 2p / r^3 = 4p/ r^3$, in the $+x$ direction. Each of the others gives an electric field $p/ r^3$, in the $-x$ direction. Hence the total field from the 6 dipoles under consideration, at the origin, is zero. Interesting!
However, if we imagine these 6 dipoles are part of a cubic lattice, then we have calculated their effect at the location exactly on top of one of the other dipoles in the lattice! So this is not a great place to do the calculation!
For a more useful insight, consider a location at the centre of a square of dipoles,  or at the centre of a cube with dipoles at the corners. In this case one has to track the signs both in ${\bf p} \cdot \hat{\bf r}$ and in $\hat{\bf r}$. One finds that the field components in the $y$ and $z$ directions cancel out (as could be seen by symmetry) but the $x$ component does not.
Finally, the case of a uniform polarization is well known (and easier to calculate).
