Biomolecular interactions determine how transcription factors recognize their DNA binding sites, how proteins interact with each other, and consequently how a biological system functions. Since many biological molecules bear considerable electric charge, electrostatic interactions are among the most important when studying biomolecular interactions. However, electrostatic interactions in biological systems are difficult to calculate accurately in practice. Aside from the significant charges carried by biomolecules such as DNA and proteins, the solvent itself, namely, water produces considerable electrostatic effects. Furthermore, hydrogen bonds, known to be involved in helix formation in both DNA and proteins, are essentially electrostatic in origin. Indeed, it seems that electrostatic effects often drive the physical-chemical processes in biological systems and, thereby, determine biological function. Therefore, any attempt to perform molecular dynamics (MD) simulations of biological systems will require an adequate description of these electrostatic forces. Previously, we have developed a rigorous method to calculated the crucial electrostatic forces in a biomolecular system. Exact analytical results have been obtained for systems with sufficient symmetry. For example we have completely analyzed a system consisting of a pair of objects with planar surfaces containing embedded charges, which is an idealized model of molecular contact. We have also developed methods to compute the electrostatic forces for a biomolecular system in which the atoms are represented by spheres. The method is rigorous in the context of the model and the accuracy can be tuned to any desired level. However, the underlying idea we developed to compute the electrostatic forces in a biomolecular system can, with some modification, be applied to a more realistic and flexible model in which a biomolecule is represented by an arbitrary surface which is decomposed into a set of small triangular patches. Once again, the electrostatic forces can be calculated to any desired accuracy, in this case by adjusting the number of triangular patches used to represent the surface. For the past year, the main effort was to implement this method in a computer code which, after optimization, will be used to demonstrate the utility of the method in a biomolecular system.