You don't need to include the XC potential. If you can ensure that the current POT.cube includes ionic (nuclear contribution) + Hartree (classical electronic contribution), in principle it should directly correspond to ESP (but the sign may need to be inverted since electron carries one unit of negative charge).

Since your system is too complex, I suggest you first using VASP to calculate a simple system such as a water molecule to check whether the ESP over the vdW surface can be normally produced via the aforementioned way.

So you think I should use the Total electrostatic potential, which also include the XC part ? Because vasp already include the ionic+Hartree in the output file.

Hartree potential is different to ESP, it only corresponds to the electrostatic potential due to electron density. In contrast, in the definiton of ESP, the nuclear charge contribution is also taken into account.

I updated Multiwfn today, this version is able to normally perform quantitative molecular surface analysis based on your cube file, the resulting map is shown below

Thank you very much for your reply. I am sure there is no problem of charge and potential.  CHG.cube is the charge density and POT.cube is the Local Hartree potential from LOCPOT output file. I think the abnormal may comes from the mismatch of VASP and Multiwfn. I will check more carefully.
   Thank you again for your suggestion and patience

I updated Multiwfn today, this version is able to normally perform quantitative molecular surface analysis based on your cube file, the resulting map is shown below

If the grid data contained in the POT.cub is ESP, the map seems to be weird. Since I am familiar with VASP, I cannot give you more advice.

Dear Lu
        I followed your suggestion and check each item for installing multiwfn and it remained the same. 
        So,  I find a new computer with win10 and tried the multiwfn_3.7_dev. And the program collapsed at   
0  // Start calculation
ESP.cub  // The cube file recording ESP
         I do not know what the problem. Can you check what happens?

On my computer, both Windows and Linux version of Multiwfn can normally load your file.
This problem of Linux version commonly implies that the Multiwfn has not been installed properly, please carefully follow the installation steps in Section 2.1.2 of Multiwfn maual.

Did you use the Linux version? I tried the Linux version on my server because my laptop is not strong enough.
Here is the output information

Differential element:   0.0021970507 Bohr^3
Summing up positive value in grid file:              1141005.3814225656
After multiplied by differential element:               2506.8466724032
Summing up negative value in grid file:             -1141005.3785599838
After multiplied by differential element:              -2506.8466661140
Summing up all value in grid file:                         0.0028625818
After multiplied by differential element:                  0.0000062892

Loaded POT.cube successfully! 
^Cforrtl: error (69): process interrupted (SIGINT)
Image              PC                Routine            Line        Source             
Multiwfn           0000000001AF4D64  Unknown               Unknown  Unknown
libpthread-2.17.s  00002B10BA0C05F0  Unknown               Unknown  Unknown
libc-2.17.so       00002B10BA3C0BED  __poll                Unknown  Unknown
libxcb.so.1.1.0    00002B10BBBB305A  Unknown               Unknown  Unknown
libxcb.so.1.1.0    00002B10BBBB0E43  xcb_connect_to_fd     Unknown  Unknown
libxcb.so.1.1.0    00002B10BBBB4A73  xcb_connect_to_di     Unknown  Unknown
libX11.so.6.3.0    00002B10B9826F0A  _XConnectXCB          Unknown  Unknown
libX11.so.6.3.0    00002B10B9817D92  XOpenDisplay          Unknown  Unknown
Multiwfn           0000000000B81EAC  Unknown               Unknown  Unknown
Multiwfn           0000000000B9F634  Unknown               Unknown  Unknown
Multiwfn           00000000007934FB  Unknown               Unknown  Unknown
Multiwfn           0000000000430922  Unknown               Unknown  Unknown
libc-2.17.so       00002B10BA2EF505  __libc_start_main     Unknown  Unknown
Multiwfn           0000000000430829  Unknown               Unknown  Unknown

I can normally load your cube file using latest version of Multiwfn, as you can see from below output information, the loading is successful.

 Input file path, for example E:\iDOLM@STER\Makoto_Kikuchi.wfn (Supported: .mwfn/wfn/wfx/fch/molden/31/chg/pdb/xyz/mol/mol2/cub/grd, etc.) Hint: Press ENTER button directly can select file in a GUI window. To reload t e file last time used, simply input the letter "o". Input such as ?miku.fch can open the miku.fch in the same folder as the file last time used. H:\POT.cube Please wait... Title line of this file: CUBE FILE CONVERTED FROM VASP format OUTTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z Total number of atoms: 196 Translation vector: X Y Z (Bohr) Vector 1: 0.130001 0.000000 0.000000 Vector 2: 0.000000 0.130001 0.000000 Vector 3: 0.000000 0.000000 0.130001 The range of x is from 0.000000 to 27.950215 Bohr, 216 points The range of y is from 0.000000 to 27.950215 Bohr, 216 points The range of z is from 0.000000 to 27.950215 Bohr, 216 points Total number of grid points: 10077696 This grid data will take up at least 76 MB memory Loading grid data, please wait... The minimum value: -0.41183300E+01 at 3.640028 19.890153 21.190163 Bohr The maximum value: 0.22817000E+00 at 16.380126 0.390003 26.650205 Bohr Differential element: 0.0021970507 Bohr^3 Summing up positive value in grid file: 1141005.3814225656 After multiplied by differential element: 2506.8466724032 Summing up negative value in grid file: -1141005.3785599838 After multiplied by differential element: -2506.8466661140 Summing up all value in grid file: 0.0028625818 After multiplied by differential element: 0.0000062892 Loaded H:\POT.cube successfully! Formula: H28 Li7 C56 N7 O42 F42 S14 Molecule weight: 2766.27302 Point group: C1

Then the main function menu is shown, you can use proper options to start analysis.

Dear Sobereva
        Thank for your reply and suggestions.
        I tried as your suggestion with Multiwfn 3.7, But it remained the same.  Here are the files of CHG.cube and POT.cube
https://pan.baidu.com/s/1QgmZluAKMfeg2vD1duVzWg     ysdk

I didn't test compatibility with cube file of VASP. Since some codes related to grid data reading has been changed during development of Multiwfn 3.7(dev), I suggest you try 3.7(dev) first, if the latest version is still unable to load your cube file, please upload the cube file to a netdisk and show me link, I will download and try to figure out the reason.

Hi
  I also tried this method by generation charge density and potential cube file from VASP calculation by files format conversion. And when I followed your suggestion and Multiwfn just stopped at:  Loaded CHG.cube successfully!  I used the Multiwfn version 3.6 and The cell in my calculation is cubic cell. I also checked the CHG.cube and POT.cube by Vesta and they all show normally.   I do not know what happens and can you give me some tips?

Your Multiwfn version is too old. Please update to the latest version first, then if this problem still exists, please compress the two cube files and then send it to my mail box (the address can be found in the first page of Multiwfn manual)

Very thanks. I hope that in the future Multiwfn could read and manipulate periodic jobs. My another question is that when I use density.cube and esp.cube (according to the above method), the summary analysis in Multiwfn says that positive surface area is 0.000 and positive average value is also NaN.
My charge data in FHI-aims show that in Al-doped graphene, Al has positive charge and carbons have near to zero but negative charges. In VMD, the color scale bar shows that Al has negative ESP, and C atoms have positive ESP. I have attached both figures. I think if this problem is solved, I could use the ESP figures in the paper.
Very thanks

I think the reason is that currently Multiwfn only supports rectangle grid, while the grid in your cube file is evidently not rectangular. The only workaround is properly building a supercell and transform the cell to rectangle type, then use FHI-aims to perform calculation and generate cube files. (Note that Multiwfn never automatically moves origin to 0,0,0)

Currently Multiwfn is mainly designed for analyzing isolated system, there is no option that can remove the upper layer. You should manually erase it using image editor.