Prof. Tian, In case someone wants to do the same thing you mentioned in the 2nd point of your answer, the manual is not clear (v3.8). I want a set of orbitals suppose homo, lumo, homo+-2 and lumo+-2 into different cube files. how can one proceed?

You should input indices of orbitals of interest together, and let Multiwfn to export cube individually, for example:

...
17 Calculate Coulomb and exchange integrals between two orbitals
18 Calculate bond length/order alternation (BLA/BOA)
20 Bond order density (BOD) and natural adaptive orbital (NAdO) analyses
3
Input orbital index. e.g. 1,3-6,8,10-11 denotes 1,3,4,5,6,8,10,11
You can also input for example: (quotation marks should not be inputted)
"h" for HOMO, "h-5" for HOMO-5, "l" for LUMO, "l+3" for LUMO+3
Input q can return
3,9,12

Please select a method to set up grid
-10 Set extension distance of grid range for mode 1~4, current:  6.000 Bohr
1 Low quality grid,    covering whole system, about 125000 points in total
2 Medium quality grid, covering whole system, about 512000 points in total
3 High quality grid,   covering whole system, about 1728000 points in total
4 Input the number of points or grid spacing in X,Y,Z, covering whole system
5 Input original point, grid spacings, and the number of points
6 Input center coordinate, number of points and extension distance
7 The same as 6, but input two atoms, the midpoint will be defined as center
8 Use grid setting of another cube file
10 Set box of grid data visually using a GUI window
11 Select a set of atoms, set extension distance around them and grid spacing
2
Coordinate of origin in X,Y,Z is      -6.000000   -7.434213   -6.901838 Bohr
Coordinate of end point in X,Y,Z is    5.952251    7.506100    6.212438 Bohr
Grid spacing in X,Y,Z is    0.166003    0.166003    0.166003 Bohr
Number of points in X,Y,Z is   73   91   80   Total:      531440

1 Output the grid data of these orbitals as separate cube files
2 Output the grid data of these orbitals as a single cube file

1 Using main function 5 to calculate grid data, in the interface of selecting real space function, choose "4 Value of orbital wavefunction" and input orbital index. After calculation of grid data is finished, choose "2 Export data to a Gaussian-type cube file in current folder".
This way can only export cube file for only one orbital each time.

2 Using the function described in Section 3.200.3 of Multiwfn manual (corresponding to subfunction 3 of main function 200).
This way can simultaneously export a batch of orbitals to a single cube file or to respective cube files.