I don't personally have experience with this but a colleague told me a few things.
There are different ways to set up a bilayer with a potential across it and the way to measure the potential across the bilayer depends on the way your system is set up. As I understood it, if you use a single bilayer with an external electrostatic potential applied to it, you have to measure the electrostatic potential at many points in your solvent and apply statistics on those measurements. However, if your periodic box contains two bilayers or a bilayer and a vacuum gap, separating two bodies of solvent with different ion concentrations, then the electrostatic potential across the bilayer apparently can be inferred from the positions of the ions (if I understood that correctly). Here are 3 references:
- From the gating charge response to pore domain movement: initial motions of Kv1.2 dynamics under physiological voltage changes. Denning EJ, Crozier PS, Sachs JN, Woolf TB., Mol Membr Biol. 2009 Dec;26(8):397-421.
- Double bilayers and transmembrane gradients: a molecular dynamics study of a highly charged peptide. Denning EJ, Woolf TB. Biophys J. 2008 Oct;95(7):3161-73. Epub 2008 Jun 27.
- Atomistic simulations of biologically realistic transmembrane potential gradients. Sachs JN, Crozier PS, Woolf TB. J Chem Phys. 2004 Dec 8;121(22):10847-51.
Of course, that doesn't tell you how to do it in CHARMM. I hope someone on this forum can fill in the blanks. If not, your best bet is to take a look at the papers, the script archive and the CHARMM documentation.