By "simulating" do you mean MD simulations? If so, I'd suggest MM/PBSA or MM/GBSA analysis of the trajectory. In particular the GB model in MM/GBSA can give you per-residue contributions to the electrostatic binding energy, including both the direct charge--charge Coulomb interactions and the effects of desolvation. Typically the former is favorable to binding, the latter unfavorable and of approximately equal magnitude. The AmberTools sute has both perl and python scripts to do this.
By "simulating" do you mean MD simulations? If so, I'd suggest MM/PBSA or MM/GBSA analysis of the trajectory. In particular the GB model in MM/GBSA can give you per-residue contributions to the electrostatic binding energy, including both the direct charge--charge Coulomb interactions and the effects of desolvation. Typically the former is favorable to binding, the latter unfavorable and of approximately equal magnitude. The AmberTools sute has both perl and python scripts to do this.
Protein - ligand interactions can be simulated directly using self-assembled monolayers (SAMs) coupled to the protein of interest on the gold surface of various types of biosensors such as Surface Plasmon Resonance and Quartz Crystal Microbalance systems commercially available. Use of of a mixed SAM consisting of a thiolalkanePEG background/backfill with a thiolalkanePEG-coupled functionality provides a reasonable simulation for membrane-bound and solution-phase protein-ligand interactions.
Regarding Ewald and MM/GBSA, there's no conflict. It is quite common to do an MD simulation with explicit water, and therefore periodic boundary conditions, and therefore Ewald, and then to do MM/GBSA analysis of the resulting trajectory. The advantage in terms of *understanding* is that the implicit-solvent GB or PB analysis gives you a nice separation of direct Coulomb interactions and solvent effects.
In a more practical point of view, there is a plugin in VMD called "APBS electrostatics" that allows you to solve exactly the electrostatic potential and visualizes it (see tutorial for vmd: http://www.poissonboltzmann.org/apbs/examples/visualization/apbs-electrostatics-in-vmd). it requires that APBS (solver) is already installed on your computer. In case you just want to calculate the electrostatic contribution of the interaction energy between your ligand and your protein from a trajectory, you may simply use the NAMD energy plugin present in VMD and load the output log file of your simulation (selection of the ligand and specific residues is possible). PME (long range electrostatic) is also included but in the end this an approximation at 0°K. If you compare several ligand it's useful. It mainly depends if your are running explicit or implicit simulations.
I would suggest you to try DelPhi, which calculates electrostatic potentials and energies of systems comprised of biological macromolecule (protein, DNA, protein-ligand etc.).
DelPhi uses the finite-difference Poisson-Boltzmann equation to calculate electrostatic potentials in and around molecules. It incorporates the effects of ionic strength, and PDB files are used as input for electrostatic potentials calculations.
Maybe I can reccommend using the new PM7 semi-empirical method if you want to assess electrostatic interactions between molecules (after MD or minimization). Use a single SCF to calculated the partial charges, and make electrostatic potential map. It seems the partial charges compare quite well with those of B3LYP-DFT methods. see: http://openmopac.net/manual/