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 surface charge method (SCM) to calculate the crucial electrostatic forces in a biomolecular system. Exact analytical results have been obtained for systems with sufficient symmetry. For example we 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. Our efforts in the past years have been to extend this exact formalism to a more realistic and flexible model in which a biomolecule is represented by an arbitrary surface. One way to do so is to incorporate intrinsic multipole moments other than just point charges. This should provide a reasonably good description of the system when the separation between biomolecules are large enough. The other way is to decompose the molecular surface into a set of small triangular patches and solve the problem fully numerically under our SCM. 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 patches such as spherical triangles and saddle surface patches can in principle alleviate this problem. However, coming up with an exact analytic mathematical expression is very challenging. This year, we have made progress in both ends: coming up with a formalism accommodating the intrinsic multipoles and providing for the first time a systematic method in using curved patches. The formalism for inclusion of intrinsic multipoles was recently published in Physical Review E; the first part of the formalism using curved surface patches was recently published in Journal of Physics Condensed Matter. We are in the process of writing up the second part of the curved surface patch formalism.