dataRange = Map[MinMax, {xg, yg, zg}]; (* the 3D range *)
espRange = {-0.03, 0.01}; (* the range of values for the electrostatic potential*)
esp = ListInterpolation[Transpose[espdata, {3, 2, 1}], dataRange]; (* Interpolate the electrostatic potential so that we can use it to color the electron density isosurface via the ColorFunction directive *)
vp = ViewPoint -> {1.3, -2.4, 2.0};
labels = Directive[FontSize -> 24, Black, FontFamily -> "Helvetica"];
densityPlot =
Legended[ListContourPlot3D[Transpose[cubedata, {3, 2, 1}],
   DataRange -> dataRange, Contours -> {0.001}, PlotRange -> Automatic,
   BoxRatios -> Automatic, Axes -> False, vp, Boxed -> False,
   MeshStyle -> Directive[Thickness[0.0075], Blue]
   , ColorFunction ->
    Function[{x, y, z},
     ColorData["NeonColors", Rescale[esp[x, y, z], espRange]]],
   ColorFunctionScaling -> False, MeshFunctions -> {esp[#, #2, #3] &},
Mesh -> {{-0.03, -0.02, -0.01, 0.00, 0.01}}],
BarLegend[{"NeonColors", espRange}, LabelStyle -> labels]]