We may take $p=(0,0)$ without loss of generality since any optimizer with $p$ nonzero can be reflected and translated. Now we want to find $(a,b),(c,d)\in\{0,1,\ldots,n\}^2$ such that the angle between the lines spanned by these vectors is as small as possible. Since $p=(0,0)$, it is equivalent to minimize the sine of the angle. To this end, recall that $$ \det\left(\begin{array}{cc}a&b\\c&d\end{array}\right)=\|(a,b)\|\|(c,d)\|\sin\theta. $$ Furthermore, the determinant is integer. The accepted answer on MSE suggests taking $(a,b)=(n-1,n)$ and $(c,d)=(n-2,n-1)$. This leads to a determinant of 1 (the smallest we could ask for since we want the angle to be nonzero). Also, the norms of these vectors are both $(\sqrt{2}-o(1))n$ (almost the largest we can ask for, since we are confined to a grid).

It remains to show that this particular choice of $(a,b)$ and $(c,d)$ is, indeed, optimal. To this end, an improvement can only come from selecting $(a,b)$ and $(c,d)$ such that the product of their norms is larger. This reduces the search space considerably: For example, for sufficiently large $n$, the optimal points necessarily lie in the set $\{(n-x,n-y):x+y\leq 3\}$. Checking these remaining cases then gives the result.

We may take $p=(0,0)$ without loss of generality since any optimizer with $p$ nonzero can be reflected and translated. Now we want to find $(a,b),(c,d)\in\{0,1,\ldots,n\}^2$ such that the angle between the lines spanned by these vectors is as small as possible. Since $p=(0,0)$, it is equivalent to minimize the sine of the angle. To this end, recall that $$ \det\left(\begin{array}{cc}a&b\\c&d\end{array}\right)=\|(a,b)\|\|(c,d)\|\sin\theta. $$ Furthermore, the determinant is integer. The accepted answer on MSE suggests taking $(a,b)=(n-1,n)$ and $(c,d)=(n-2,n-1)$. This leads to a determinant of 1 (the smallest we could ask for since we want the angle to be nonzero). Also, the norms of these vectors are both $(\sqrt{2}-o(1))n$ (almost the largest we can ask for, since we are confined to a grid).

It remains to show that this particular choice of $(a,b)$ and $(c,d)$ is, indeed, optimal. To this end, an improvement can only come from selecting $(a,b)$ and $(c,d)$ such that the product of their norms is larger. This reduces the search space considerably: For example, for sufficiently large $n$, the optimal points necessarily lie in the set $\{(n-x,n-y):x+y\leq 3\}$. Checking these remaining cases then gives the result.