One way to do this is to compute the homography (planar perspective transformation) that maps the square board onto the detected quadrangle in the image, then use that transformation to map grid line endpoints onto the image. There are some advantages to this method. First, this computation is available as a library function in many graphics APIs, so you can do it with basically no work on your part. Second, it allows you to map any point in the plane of the board onto the image, which can be handy depending on what else you’re going to be doing. Third, it might allow you to blast the grid out as a texture with suitable alpha for transparency. If this function isn’t available in the APIs that you’re using, you can do this computation yourself. You can find a description of the process here. For this application, computing the homography is particularly simple: you can use the unit square for the board, in which case you don’t have to compute a second intermediate matrix when using the linked method—the homography is simply $A^{-1}$.
Since you’ve computed the vanishing points of the board edges, another way to go is to use cross-ratios to compute the grid line endpoints: if $A$, $B$, $C$, $D$ are distinct colinear points, then the ratio $$(A,B;C,D) = {AC \cdot BD \over AB \cdot CD}$$ is invariant under projective transformations. If you have three fixed pairs of corresponding points on a line and its image, you can use the equality of the cross-ratios to find the image of any other point on the line.
For each of the edges you’ve got the images of the endpoints $P_1$ and $P_2$ and can compute the vanishing point $P_\infty$, which gives the necessary three pairs. The latter is the image of a point at infinity in the overhead view of the board, which presents a slight difficulty as its distance from any of the other points on the line is infinite, but that can be dealt with by introducing homogeneous coordinates on the line and using determinants of pairs of points. Taking the unit square as the original chess board, we have $\mathbf 0=[0:1]$, $\mathbf 1=[1:1]$, $\mathbf\infty=[1:0]$ and $\mathbf\lambda=[\lambda:1]$, with cross-ratio $$(\mathbf 0,\mathbf\lambda;\mathbf 1,\mathbf\infty) = {\begin{vmatrix} 0 & 1 \\ 1 & 1 \end{vmatrix} \cdot \begin{vmatrix} \lambda & 1 \\ 1 & 0 \end{vmatrix} \over \begin{vmatrix} 0 & \lambda \\ 1 & 1 \end{vmatrix} \cdot \begin{vmatrix} 1 & 1 \\ 1 & 0 \end{vmatrix}} = {(-1)\cdot(-1) \over (-\lambda)\cdot(-1)} = \frac1\lambda.$$ On the image side of things, we have the endpoints $P_0=[x_0:y_0:1]$, $P_1=[x_1:y_1:1]$ and the vanishing point $P_\infty=[x_\infty:y_\infty:1]$ (I’m assuming for now that the vanishing point is finite). We want the value of $\mu$ for which $P_\mu = (1-\mu)P_0+\mu P_1$ is the image of $\mathbf\lambda$, that is, for which $(P_0,P_\mu;P_1,P_\infty)=(\mathbf 0,\mathbf\lambda;\mathbf 1,\mathbf\infty)$. We can compute distances between these various image points and then get the cross-ratio, but it can also be computed by choosing a point not on the line and using determinants of $3\times3$ matrices formed from various combinations of the points. A convenient choice is the origin $O=[0:0:1]$, because the determinants then become wedge products of the Cartesian coordinates of pairs of points. (If the linear extension of the edge happens to pass through the origin, just translate the points so that it doesn’t.) The cross-ratio of the image points is then $$(P_0,P_\mu;P_1,P_\infty) = {[O,P_0,P_1]\cdot[O,P_\mu,P_\infty] \over [O,P_0,P_\mu]\cdot[O,P_1,P_\infty]} = {(x_0y_1-y_0x_1) \left((1-\mu) (x_0y_\infty-y_0x_\infty)+\mu (x_1y_\infty-y_1x_\infty)\right) \over \mu (x_0y_1-y_0x_1) (x_1y_\infty-y_1x_\infty)}.$$ Setting the two cross-ratios equal to each other and solving for $\mu$ produces $$\mu = {\lambda (x_0y_\infty-y_0x_\infty) \over (1-\lambda) (x_1y_\infty - y_1x_\infty) + \lambda (x_0y_\infty - y_0x_\infty)}.$$ This is an instance of a linear fractional transformation. Note that even though the cross ratio is undefined for $\lambda=0$, this formula produces $\mu=0$, as we’d like.
If a pair of edges is parallel, or nearly so, with the vanishing point far away from the image of the board, then you can just subdivide them evenly instead of going through the above mapping. (In this case, $P_\infty=[x_1-x_0:y_1-y_0:0]$, and if you work through the cross-ratio computation, you end up with $1/\mu$, i.e., $\mu=\lambda$.)
