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 functions. 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 calculate 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. There is, of course, fundamental effect that require more than electrostatics to understand. For the past few years, we have been intensively testing the accuracy and adequacy of surface triagularization by comparing the exact analytical result with those numerically obtained via surface triangularization. We found that the usual triangularization method (using flat triangles) is generally pathological. This pathology is somewhat expected as we know from classical electrostatics that pointy vertices and sharp edges will acquire infinite electric field. This pathological situation, plaguing all related numerical methods, disappears only if all the edges across all ajoining pathces are smooth. Using curved pathes such as spherical triangles and saddle surface pathes can in principle alleviate this problem. However, coming up with an exact analytic mathematical expression is very challenging. This year, we have finally met this challenge by providing for the first time a systematic method in using curved pathes. We are in the process of writing this up. In addition, we have also tested the validity of the classical method at atomic scale. For systems composed of atoms in the S state, our formulation for dielectric spheres works amazingly well, ofen better than crude quantum mechanical methods such as the DFT with moderate basis set. The results of this discovery has been written and the manuscript was recently published in EuroPhysics Letters as one of the editors choice articles of year 2016.