Creating electrostatic potential with external charge

Hello,

I was trying to plot the electrostatic potential of a system embedded in external charges, and below is the input file

$rem
jobtype             sp
max_scf_cycles          250
geom_opt_max_cycles     500
method          wb97m-v
basis               def2-tzvp
QM_MM_INTERFACE         Zeolite
force_field         charmm27
user_connect        true
symmetry            off
sym_ignore          true
print_input         true
qmmm_print          false
qm_mm               true
qmmm_full_hessian       false
mem_total           12000
mem_static          450
geom_opt_dmax           80
pop_mulliken        false
XC_GRID         3
SKIP_CHARGE_SELF_INTERACT 1
esp_grid        -2  ! with esp of external charges
make_cube_files  1
$end

$plots
 plot esp
 70 -40 30
 70 -40 30
 70 -40 30
 0 0 0 0
 0
$end

$qm_atoms
1:19
$end

$molecule
0 1
H     7.525320    14.506700   12.920000   0  2   0   0   0
O     7.444290    13.808430   13.575020   0  1   3   0   0
Ti    7.650510    12.157530   14.094890   0  2   4   5   6
O     6.305990    11.048370   13.593010   0  3   7   0   0
O     8.192280    12.030100   15.840290   0  3   8   0   0
O     9.007710    11.321810   13.304180   0  9   3   0   0
Si    5.326610    9.749430    13.469470   0  12  33  11  4
Si    8.822730    10.748150   16.651220   0  10  14  5   13
Si    9.811420    9.976980    13.847360   0  10  6   31  32
O     9.500570    9.792950    15.469240   0  9   8   0   0
O     4.504500    9.620800    12.058970   0  7   15  0   0
O     4.284680    9.898860    14.719090   0  7   16  0   0
O     7.571970    10.049940   17.441580   0  17  8   0   0
O     9.961260    10.980290   17.821670   0  18  8   0   0
Si    3.202150    8.901300    11.373070   0  11  20  21  22
Si    3.151090    8.826050    15.244640   0  19  12  23  24
Si    6.548170    8.807490    17.725610   0  25  26  27  13
Si    10.467760   10.399850   19.294170   0  30  28  29  14
H     2.628330    9.257670    16.704750   0  16  0   0   0
O     3.655180    8.243890    9.904800   -1  1934    15  0   0
O     1.883800    9.864290    11.740140  -1  1882    15  0   0
O     2.563570    7.441520    11.893590  -1  1883    15  0   0
O     1.846410    8.917370    14.296250  -1  1882    16  0   0
O     3.830350    7.293790    15.141810  -1  1957    16  0   0
O     5.690230    8.446710    16.392960  -1  17  1957    0   0
O     5.316810    9.183850    18.743400  -1  17  1935    0   0
O     7.517080    7.487690    17.980940  -1  17  2005    0   0
O     11.837100   9.465600    19.144540  -1  18  2098    0   0
O     10.618340   11.685360   20.339020  -1  18  2058    0   0
O     9.357610    9.406280    19.991690  -1  18  2027    0   0
... (another ~2000 lines of external charges)
$end

However, the cube file I obtained has NAN inside at the beginning and can’t be plotted, which looks like this

Cube file for electrostatic potential
Generated using Q-Chem
   33     0.000000     0.000000         -nan
   80     0.000000     0.000000     0.000000
    0     0.000000    -0.000000     0.000000
-2091130759     0.000000     0.000000         -nan
    1     1.000000    14.220794    27.413690    24.415262
    8     8.000000    14.067669    26.094151    25.653070
   22    22.000000    14.457369    22.974402    26.635482
    8     8.000000    11.916594    20.878393    25.687066

Any ideas on how to fix this problem? Thank you.

can you please send a minimal example that exhibits this bug? (e.g., small molecule and small basis set so that it runs quickly). Please send full input for such an example.