I directly use Multiwfn to generate cube file of ESP here, you can also use other code. Boot up Multiwfn and input
examples\H2O.fch
5  // Calculate grid data
12  // ESP
2  // Medium quality grid
2  // Export grid data to cube file

Now you have totesp.cub in current folder.

Next, we calculate ESP charges based on the totesp.cub. You should set "iuserfunc" in settings.ini to -1 to make user-defined function equivalent to linearly interpolated function based on the loaded grid data. Then boot up Multiwfn and input

totesp.cub
100  // Other function (Part 1)
23  // Fit function distribution to atomic value
5  // Set constraint for total value
0
1  // Select a real space function and start calculation
100  // User-defined real space function

You will see

 Center       X           Y           Z             Value
    1O     0.000000    0.000000    0.225257       -0.725089
    2H     0.000000    1.434977   -0.901029        0.362542
    3H     0.000000   -1.434977   -0.901029        0.362548
Sum of values:    0.000000
RMSE:    0.002612   RRMSE:    0.124550

Note that the subfunction 23 in main function 100 is a general module for fitting a function as atomic values using the same algorithm as ESP fitting, see Section 3.100.23 of Multiwfn manual for detail. Since this module employs MK-type of fitting points, the "Value" shown above corresponds to Merz-Kollman ESP fitting charge.

Be aware that when generating the ESP cube file, the buffer distance between the boundary atom and the box should be large enough. Under the default setting of fitting points (standard MK fitting points), the buffer should be no less than two times of vdW radius of boundary atoms, otherwise some fitting points will be out-of-box, and ESP at those points cannot be evaluated by interpolation.

Let me ask another question to find an alternative.
Is a text file including grid points and ESP values enough to calculate the ESP charges in Multiwfn?

If yes, I will prepare grid points and ESP values using density cube file.

I’m wondering if we could calculate ESP fitting charge using density cube file.

I understood that ESP charges by such as CHELPG and MK were easily calculated using wfn or fchk files.
But I couldn’t find the way to do the same thing using only density cube file.

Thanks in advance for you help!