Boundary integral methods have long been used to simulate fully nonlinear, time-dependent water waves. For surface waves in 3D it appears that choices of discretization can strongly affect the numerical stability as well as accuracy. We will discuss the formulation, accuracy, and stability in the case of a doubly periodic surface with Lagrangian markers. One version of the method, not too different from ones in use, has been shown to converge to the actual solution as long as it is smooth. In the method to be described, the normal velocity at the surface is determined from the potential by solving an integral equation. The singular integrals are computed using a simple quadrature rule with a regularization and then adding a correction. A sign condition for the discrete single layer potential seems to be important in maintaining wavelike behavior. The error analysis leading to convergence depends on preserving a structure in the linearized discrete equations which appears in the original system.