Calculating electrostatic potential using "ESPGrid" file

Questions
#1

I want to calculate electrostatic potential for my systems using a manually provided grid information in an external file. According to the manual this can be done using setting the ESP_GRID =n, where n is the number of grid points in the ESPGrid file.

I encountered a possible bug where the grid printed in the plot.esp file does not match the grid provided in the ESPGrid file. Here is the input and output grid:

ESPGrid:

-1.87308530    0.54366160    2.60058023
-1.34808530   -0.36566506    2.60058023
-1.08190665    2.34380368    1.83192688
-2.01404996    1.91810810    1.83192688
-2.56807041    1.05603546    1.83192688
-2.56807041    0.03128774    1.83192688
-1.56775556    2.50719571    0.78192688
-2.39495787    1.93621918    0.78192688
-2.86206311    1.04622450    0.78192688
-2.86206311    0.04109871    0.78192688

plot.esp

Grid positions and esp values from SCF density (all in a.u.)

       X            Y            Z           GS

  -0.17718237E+01   0.52919838E+01  -0.77511187E+00  -0.61382639E-01
  -0.23984743E+01   0.44931714E+01  -0.50732899E+00  -0.49403268E-01
   0.32499777E+00   0.51178931E+01  -0.94928234E+00  -0.97358514E-01
  -0.39227236E+00   0.56349953E+01  -0.14671982E+01  -0.65959505E-01
  -0.13648564E+01   0.56733032E+01  -0.17876988E+01  -0.48467693E-01
  -0.22839658E+01   0.52206545E+01  -0.18090276E+01  -0.38633580E-01
   0.56652072E+00   0.50511552E+01  -0.20905355E+01  -0.60713087E-01
  -0.24328775E+00   0.54246994E+01  -0.25541656E+01  -0.45553206E-01
  -0.12096365E+01   0.53849249E+01  -0.28277819E+01  -0.35759615E-01
  -0.21111467E+01   0.49409435E+01  -0.28487022E+01  -0.28941210E-01

Input file:

$molecule
0 1
    N       0.84399998       0.45600000       0.92299998
    H       1.37199998       0.49500000       1.79799998
    C       1.64300001       1.12500000      -0.11300000
    H       2.00000000       2.09699988       0.24300000
    H       0.98699999       1.30700004      -0.97100002
    C       2.77900004       0.15300000      -0.50400001
    H       3.68899989       0.38400000       0.05900000
    H       3.03099990       0.22200000      -1.56599998
    C       2.23099995      -1.24000001      -0.11200000
    H       2.86599994      -1.70700002       0.64700001
    H       2.18499994      -1.92799997      -0.96200001
    C       0.82999998      -0.94099998       0.46900001
    H       0.54200000      -1.60899997       1.28600001
    H       0.06800000      -1.02699995      -0.31400001
$end

$rem
   METHOD     hf
   BASIS      6-31g**
   IANLTY     200
   ESP_GRID   10        !643
$end

$plots
   plot the electrostatic potential on a line
   1   0.0   0.0
   1   0.0   0.0
  1  0.0   0.0
   0  0  0  0
   0
$end

Please let me know if there is any mistakes in the input file or in the formatting of the ESPGrid file.

#2

Please add the following to your $rem:

$rem
no_reorient true
sym_ignore true   ! or point_group_symmetry = false
....other stuff...
$end

What NO_REORIENT=TRUE will do is to make sure that Q-Chem uses the user’s coordinate system (in $molecule), rather than rotating into the standard nuclear orientation (which is the default behavior). Since your ESPGrid is presumably made in the user’s frame, that’s the only choice that makes sense. Turning off symmetry avoids an occasional problem when you turn off reorientation for a molecule with nontrivial point group symmetry.

This is pretty subtle and I will add a cautionary note to the manual.

#3

That solved the problem, thank you so much.

#4

For future reference, I checked in a fix for Q-Chem 6.2.2 so that one can set NO_REORIENT = TRUE without messing around with the symmetry keywords. And the cautionary note in the manual should also appear with v. 6.2.2.