ChimeraX can color some surface, which could be a molecular surface, by the values in a 3D map, which could be electrostatic potential (ESP).  See "color electrostatic"

<https://rbvi.ucsf.edu/chimerax/docs/user/commands/color.html#map>

The ESP map can come from a Poisson-Boltzmann calculation (APBS, DelPhi, others) or from a more simplistic Coulombic calculation which can be done directly in ChimeraX.  

Poisson-Boltzmann calculations require atomic coordinates, partial charges, and radii, but the different programs differ in how you input these values. 

If I recall correctly, APBS takes as input a pqr file (q means charge, r means radii) which would usually be generated by PDB2PQR, but you said PDB2PQR throws out the nonstandard residues like cofactors, I guess because it doesn't have charges for them.  ChimeraX assigns radii and can assign atomic partial charges (via lookup table for standard residues and Antechamber calculation for nonstandard residues like cofactors).  However, ChimeraX does not write pqr format.  Also in ChimeraX, these partial atomic charges are for all atoms (i.e. all hydrogens explicit) and not a united-atom model (in which hydrogens on carbons are merged with those carbons so that they become implicit) which is often used for Poisson-Boltzmann calculations.  So if you wanted to go the APBS route,  in ChimeraX you would add hydrogens and charges, then save the cofactor radii and atomic partial charges (atom attributes "radius" and "charge") in the attribute file format, and if you were using a united-atom charge model, manually generate united atoms from CH, CH2, and CH3 groups in the cofactor by adding the hydrogen charges to those of the carbons where they are attached and figuring out the appropriate united-atom radii, then reformat the cofactor atom coordinates, radii and charge values into pqr format and then add these lines to your protein pqr file from PDB2PQR, then use the resulting combined pqr file as input to APBS to get the ESP map that includes cofactor contributions.

Actually I don't even know if APBS would accept that pqr file or not, but that is what I imagine you would try.  So let's assume the APBS calculation will run using that pqr file.

Then you would open the APBS ESP map that includes cofactor contributions in ChimeraX and color the molecular surface.  If the cofactor is covalently attached to the protein then the protein molecular surface may already include it.  However, if the cofactor is not covalently attached, you would need to tell ChimeraX specifically to lump it inside of the same surface as the protein with the "include" option of "surface" when you are calculating and displaying the molecular surface.

The "coulombic" approach wil do more of these steps for you, given the "include" trick as well as the "surfaces" option of the "coulombic" command, but calculating ESP with the relatively simplistic Coulomb's law method.  For example, 1qs0 chain A includes a noncovalently bound TPP residue and I could show the surface including the TPP and colored including its potential contribution:

open 1qs0

surf include :tpp
coulombic surfaces #1.2

<https://rbvi.ucsf.edu/chimerax/docs/user/commands/surface.html#options>

<https://rbvi.ucsf.edu/chimerax/docs/user/commands/coulombic.html>

See first image attached with :tpp selected (green outline surface patch)

If you don't use the "include" and "surfaces" options in these commands, the surface will be redrawn to include only the protein part with the TPP outside of it, like other ligands and the Coulombic ESP will not include the TPP contribution.    See second image attached.

I hope this helps,

Elaine