Ammm:Lig ion params

From biowiki
Jump to navigation Jump to search

Introduction

Parameters in AMOEBA

  • Water first being parameterized
  • ions, organic molecules, proteins and nucleic acids (under improvement)

How to get parameters

  • Experiment
  • Ab initio QM
  • Other force field

Unique parameters in AMOEBA

  • Multipole (up to quadrupole)
  • Polarizability
  • Damping factor

Electrostatic Parameters of AMOEBA

Metal Ion Parameterization

Why model metal ions?

Ions play important roles in biological function

Nearly half of all proteins contain metal ions
  • Ion channel (Na+/K+)
  • Protein biosynthesis (Mg2+)
  • Coagulation cascade (Ca2+)
  • Pathogenesis of viruses (Zn2+)

What are the challenges

  • High charge density
  • Strong polarization effect
  • Charge transfer effect
    • Can’t be explained by coulumb’s law. Quantum effect.
  • Split orbital
    • Transition metal, e.g. Zn++

Polarizable vs non-polarizable force field in ion modeling

  • Both polarizable and nonpolarizable force fields handle monovalent ions reasonably.
  • In case of divalent ions, polarizable force fields give much better result. [1].

Above figures adapted from Warshel et al. JCTC 2007.

Parameterization

What parameters are needed?

  • Bonded parameters (bond, angle, torsion and out-of-plane)
  • Non-bonded parameters
    • Van der Waals
    • Electrostatics

How to parameterize

  • Charge is known (2+)
  • Dipole and quadrupole are 0
  • Polarizability determined from B3LYP/6-31G* calculation
  • Van der Waals and damping factor (R,ε,a) are determined by fitting to QM dimer energy
    1. Optimize (Cation-water) dimer structure with B3LYP/6-31G*
    2. Vary the separation between cation and water molecule with fixed optimized geometry.
    3. Compute the binding energy versus separation with both QM and AMOEBA. Binding energy = total energy - isolated energy of ion - isolated energy of water
    4. Tune the parameters to get the best fit.

HF: Hartree-Fock method is an approximate method for the determination of the ground-state wave function, sometimes called SCF (self-consistent field)

MP2: Møller–Plesset theory is based on HF method. It improves on the Hartree–Fock method by adding electron correlation effects to the second order. Mp3 (third), MP4 (fourth).

B3LYP: Hybrid functional. Approximation to density functional theory that incorporates a portion of exact exchange from Hartree-Fock theory with exchange and correlation from other sources (ab initio or empirical)

Case study

Example of gaussian input file for polarizability calculation

 %Mem=800MB 
 %Chk=polar.chk

 #B3LYP/6-311++G(2d,2p) SP Polar MaxDisk=10GB

   sp energy

   2,1
  Ca1   0.871756   -0.000003    0.000002  

Output file

      ....skipping...
   Electronic spatial extent (au):  <R**2>=    15.4404
   Charge=     2.0000 electrons
   Dipole moment (field-independent basis, Debye):
      X=     0.0000    Y=     0.0000    Z=     0.0000  Tot=     0.0000
   Quadrupole moment (field-independent basis, Debye-Ang):
     XX=    -6.9226   YY=    -6.9226   ZZ=    -6.9226
     XY=     0.0000   XZ=     0.0000   YZ=     0.0000
   Traceless Quadrupole moment (field-independent basis, Debye-Ang):
     XX=     0.0000   YY=     0.0000   ZZ=     0.0000
     XY=     0.0000   XZ=     0.0000   YZ=     0.0000
   Octapole moment (field-independent basis, Debye-Ang**2):
    XXX=     0.0000  YYY=     0.0000  ZZZ=     0.0000  XYY=     0.0000
    XXY=     0.0000  XXZ=     0.0000  XZZ=     0.0000  YZZ=     0.0000
    YYZ=     0.0000  XYZ=     0.0000
   Hexadecapole moment (field-independent basis, Debye-Ang**3):
    XXXX=    -3.3848 YYYY=    -3.3848 ZZZZ=    -3.3848 XXXY=     0.0000
    XXXZ=     0.0000 YYYX=     0.0000 YYYZ=     0.0000 ZZZX=     0.0000
    ZZZY=     0.0000 XXYY=    -1.1283 XXZZ=    -1.1283 YYZZ=    -1.1283
    XXYZ=     0.0000 YYXZ=     0.0000 ZZXY=     0.0000
   N-N= 0.000000000000D+00 E-N=-1.591920174717D+03  KE= 6.759685000409D+02
   Symmetry AG   KE= 4.581170345941D+02
   Symmetry B1G  KE= 2.744378782701D-61
   Symmetry B2G  KE= 2.109222553139D-61
   Symmetry B3G  KE= 3.141249198819D-61
   Symmetry AU   KE= 0.000000000000D+00
   Symmetry B1U  KE= 7.261715514895D+01
   Symmetry B2U  KE= 7.261715514895D+01
   Symmetry B3U  KE= 7.261715514895D+01
   Exact polarizability:   3.216   0.000   3.216   0.000   0.000   3.216
  Approx polarizability:   4.622   0.000   4.622   0.000   0.000   4.622
  1\1\GINC-NODE12\SP\RB3LYP\6-311++G(2d,2p)\Ca1(2+)\OSCAR\05-Jan-2006\0\
  \#B3LYP/6-311++G(2D,2P) SP POLAR MAXDISK=10GB\\sp energy\\2,1\Ca,0,0.8
  71756,-0.000003,0.000002\\Version=IA32L-G03RevC.02\State=1-A1G\HF=-676
  .9057847\RMSD=6.491e-10\Dipole=0.,0.,0.\Polar=3.2162085,0.,3.2162085,0

....skipping...

Bohr Radius = 0.529 177 2083 Å

Basis set superposition error (BSSE). In quantum chemistry, calculations of interaction energies are susceptible to BSSE. As two molecules approach one another, their basis functions overlap. Counterpoise method to eliminate the error. To define basis set internally, use the keyword Gen.

EMSL Basis Set Exchange


Input file of Gaussian energy calculation

 #MP2/Gen SP counterpoise=2 MaxDisk=10GB

sp energy
 
 2,1,2,1,0,1
Ca1   -0.628244   -0.000003    0.000002    1
O2    1.471756   -0.000020   -0.000018   2
H3    2.069508   -0.764665    0.000058   2
H4    2.069178    0.764883    0.000058   2
 
 CA  0
 S  19  1.00
      7503960.00         0.224000000E-05
      1123486.00         0.174400000E-04
   
 ....skipping....
   P   1  1.00
     0.144000000          1.00000000    
P   1  1.00
     0.580000000E-01      1.00000000    
P   1  1.00
     0.231600000E-01      1.00000000    
D   4  1.00
      26.2708300         0.734600000E-02
      7.30782000         0.409010000E-01
      2.45498000         0.130065000    
     0.887140000         0.244151000    
D   1  1.00
     0.309250000          1.00000000    
D   1  1.00
     0.103160000          1.00000000    
D   1  1.00
     0.337400000E-01      1.00000000    
F   1  1.00
     0.488300000          1.00000000    
F   1  1.00
     0.116500000          1.00000000    
G   1  1.00
     0.146600000          1.00000000    
 ****
H O 0
 cc-pvqz
 ****

Output

Hexadecapole moment (field-independent basis, Debye-Ang**3):
XXXX=   -20.3681 YYYY=    -5.8454 ZZZZ=    -6.0461 XXXY=    -0.0019
XXXZ=     0.0006 YYYX=     0.0008 YYYZ=     0.0000 ZZZX=     0.0001
ZZZY=     0.0000 XXYY=    -1.9087 XXZZ=    -5.9481 YYZZ=    -2.3208
XXYZ=     0.0000 YYXZ=     0.0002 ZZXY=     0.0000
N-N= 9.069616433324D+00 E-N=-1.988243128533D+02  KE= 7.595642684817D+01
Counterpoise: corrected energy =    -750.764481575760
Counterpoise: BSSE energy =       0.003336990775
1\1\GINC-NODE5\SP\RMP2-FC\Gen\Ca1H2O1(2+)\OSCAR\29-Dec-2005\0\\#MP2/GE
N SP COUNTERPOISE=2 MAXDISK=10GB\\sp energy\\2,1\Ca,0,0.471756,-0.0000
03,0.000002\O,0,1.471756,-0.00002,-0.000018\H,0,2.069508,-0.764665,0.0
00058\H,0,2.069178,0.764883,0.000058\\Version=IA32L-G03RevC.02\State=1
-A\HF=-76.0635232\MP2=-76.3473292\RMSD=2.782e-09\PG=C01 [X(Ca1H2O1)]\\
@

1 Hartree = 627.509 kcal/mol

TINKER input (xyz)

    4  Ca-wat
     1  Ca+   -0.628244   -0.000003    0.000002  18 
     2  O     1.471756   -0.000020   -0.000018    61     3     4
     3  H     2.069508   -0.764665    0.000058    62     2
     4  H     2.069178    0.764883    0.000058    62     2

TINKER output (command analyze)

Total Potential Energy :                -10.9037 Kcal/mole

Intermolecular Energy :                 -11.3330 Kcal/mole

Energy Component Breakdown :           Kcal/mole      Interactions

Bond Stretching                           0.1828              2
Angle Bending                             0.2243              1
Urey-Bradley                              0.0223              1
Van der Waals                            -0.0492              3
Atomic Multipoles                       -10.0915              3
Polarization                             -1.1923              3

Small ligand parameterization

Build your ligand

(1) Make your molecule with Chem3D Ultra or Arguslab (windows)

(2) Minimize the structure. Once the molecule is built, minimize energy with any of the methods available in Chem3D. There are other ways to minimize energy: MOPAC, MM2, Gaussian, Mechanics. You can use any of them as long as it works.

(3) Save the structure as xyz file.

Save as Tinker MM3 Input (*.xyz), benz.xyz

Optimize structure with QM

(1) Convert .xyz to .com

Just copy the atoms and their coordinates into the .com file with following format:

 %chk=benz.chk
 %mem=50MW
 %nproc=2
  # opt hf/6-31g*
  
   opt energy
 
   1 1
C    0.310792    0.313000   -0.224530
C   -0.125208    0.944000    0.941470
C   -0.473208    0.250000    2.081470
C   -0.397208   -1.110000    2.090470
C    0.025792   -1.764000    0.961470
C    0.376792   -1.071000   -0.177530
C    0.655792    1.023000   -1.380530
N    1.378792    0.421000   -2.305530
N    0.251792    2.268000   -1.515530
N   -0.666208   -1.665000    2.975470
H   -0.197208    2.021000    0.966470
H   -0.807208    0.776000    2.964470
H    0.078792   -2.841000    0.960470
H    0.706792   -1.627000   -1.043530
H    1.626117   -0.505205   -1.989348
H    1.663792    0.848000   -3.156530
H   -0.246151    2.554668   -0.678924
H    0.471792    2.841000   -2.290530
H   -0.548431   -2.662838    2.799895
H   -1.663792   -1.553551    3.156530
 <blank line>


(2) Run Gaussian program to optimize the structure.

Single Point calculation

(1) Make the single point file (.com)

Extract the optimized coordinates from opt output file (.log) and set up a single point calculation

 %Mem=500MB 
 %Nosave 
 %Chk=benzsp.chk
 %Nproc=2
  #MP2/6-311++G(2d,2p) Sp Density=MP2 MaxDisk=960MW

  sp energy
 
   1 1
C1  0.293566  0.297615  -0.240256
C2  -0.088764  0.943618  0.941596
C3  -0.462820  0.230827  2.043653
C4  -0.490542  -1.177405  2.018026
C5  -0.112718  -1.822504  0.824061
C6  0.275038  -1.102026  -0.268251
C7  0.697450  1.057005  -1.402781
N8  1.606270  0.574513  -2.231120
N9  0.160802  2.239323  -1.645792
N10  -0.865410  -1.882805  3.098002
H11  -0.050701  2.013222  1.010215
H12  -0.729196  0.743799  2.944741
H13  -0.141810  -2.891482  0.774237
H14  0.525857  -1.628596  -1.168273
H15  2.146762  -0.228883  -1.984752
H16  1.818078  1.010838  -3.106430
H17  -0.641774  2.550071  -1.138264
H18  0.503042  2.844767  -2.365315
H19  -0.880774  -2.880211  3.088084
H20  -1.139941  -1.431611  3.944290
 <blank line>

We use higher energy level and basis set MP2/6-311++G(2d,2p) in order to get accurate energy result. Lower energy level has been used for optimization since the energy is sensitive to energy level while structure is not.

Multipoles calculation

Distributed Multipole Analysis (DMA) is a technique for describing a molecular charge distribution by using local multipoles at a number of sites within a molecule. It gives a much more accurate representation of the charge density than a single-point multipole expansion.

Generalized Distributed Multipole Analysis (GDMA) program was developed by Anthony Stone. It is a Fortran 90 program for performing Distributed Multipole Analysis of wavefunctions calculated using the Gaussian program system. That is, it calculates electric multipole moments at the atomic positions, or at other specified sites, that can give a very accurate representation of the electrostatic field of the molecule.

GDMA program

Multipole information can be extracted from formatted check file. Since check file is binary, need to run formchk to make it readable

(1)Do "formchk benzsp.chk". This will generate benzsp.fchk file.

(2) Make GDMA input file.

 Title 

 File benzsp.fchk density MP2
 Angstrom
 AU
 Multipoles
 Limit 2
 Radius H 0.31
 Punch benzsp.punch

 Start

 Finish

GDMA program sets all the atom radii to be equal by default, except for hydrogen. The use of equal radii for all sites is the most efficient choice for convergence of the resulting multipole expansion of the electrostatic potential. Default radius is 0.65 Å. Different value may be explicitly specified. 0.325 Å for hydrogen, 1.11Å for chloride is found to give more acceptable values.


(3) Run GDMA to calculate multipoles of this ligand.

  /opt/gdma-2.2/bin/gdma < gdmain > benzsp.gdmaout"

There are two output files benzsp.punch and benzsp.gdmaout.

Standard DMA

Positions and radii in angstrom
Multipole moments in atomic units, ea_0^k for rank k

C          x =  0.640305  y = -0.000002  z = -0.000115
          Maximum rank =  2   Relative radius =  0.650
                  Q00  =  -0.101884
|Q1| =   0.024901  Q10  =   0.000007  Q11c =  -0.024901  Q11s =  -0.000018
|Q2| =   1.365955  Q20  =  -1.283977  Q21c =   0.000260  Q21s =   0.375277
                  Q22c =   0.276408  Q22s =   0.000001

C          x = -0.077336  y =  1.181820  z =  0.220245
          Maximum rank =  2   Relative radius =  0.650
                  Q00  =  -0.095143
|Q1| =   0.204753  Q10  =   0.005562  Q11c =  -0.073651  Q11s =  -0.190967
|Q2| =   1.225313  Q20  =  -1.139120  Q21c =   0.041511  Q21s =   0.427129
                  Q22c =   0.084762  Q22s =  -0.111577

C          x = -1.442045  y =  1.184957  z =  0.229039
          Maximum rank =  2   Relative radius =  0.650
                  Q00  =  -0.091024
|Q1| =   0.176199  Q10  =  -0.027643  Q11c =   0.050291  Q11s =  -0.166592
|Q2| =   1.177708  Q20  =  -1.106636  Q21c =  -0.000588  Q21s =   0.379762
                  Q22c =   0.119889  Q22s =   0.061321

C          x = -2.168671  y =  0.000011  z = -0.000002
          Maximum rank =  2   Relative radius =  0.650
                  Q00  =   0.033606
|Q1| =   0.004956  Q10  =  -0.000020  Q11c =   0.004956  Q11s =  -0.000009
|Q2| =   1.025836  Q20  =  -0.964500  Q21c =  -0.000179  Q21s =   0.302654
                  Q22c =   0.174581  Q22s =   0.000026

   ....skipping.....
 
H          x = -4.032047  y = -0.837065  z = -0.154879
          Maximum rank =  2   Relative radius =  0.310
                  Q00  =   0.222435
|Q1| =   0.084255  Q10  =   0.009661  Q11c =   0.065456  Q11s =   0.052164
|Q2| =   0.181344  Q20  =  -0.142929  Q21c =   0.017003  Q21s =   0.053053
                  Q22c =  -0.029051  Q22s =   0.092245

H          x = -4.032046  y =  0.836924  z =  0.155699
          Maximum rank =  2   Relative radius =  0.310
                  Q00  =   0.222432
|Q1| =   0.084256  Q10  =  -0.009709  Q11c =   0.065457  Q11s =  -0.052156
|Q2| =   0.181345  Q20  =  -0.142895  Q21c =  -0.017130  Q21s =   0.053161
                  Q22c =  -0.029030  Q22s =  -0.092222

Total multipoles
referred to origin at x =   0.000000,  y =    0.000000,  z =    0.000000
                  Q00  =   0.999990
|Q1| =   1.655329  Q10  =   0.000477  Q11c =   1.655328  Q11s =  -0.000178
|Q2| =  39.203832  Q20  = -29.420772  Q21c =  -0.003923  Q21s =  -2.843583
                  Q22c =  25.754080  Q22s =   0.000629


(4) Edit the multipoles from the GDMA output with poledit.x.

  The Multipole Editing Facility can Provide :

   (1) Multipole Parameters from GDMA Output
   (2) Alter Local Coordinate Frame Definitions
   (3) Removal of Intramolecular Polarization

  Enter the Number of the Desired Choice :  1

Global Frame Cartesian Multipole Moments :

Site:       1         Name:   C         Atomic Number:       6

Coordinates:            0.640305      -0.000002      -0.000115

Charge:                 -0.10188
Dipole:                 -0.02490       -0.00002        0.00001
Quadrupole:              0.88136
                         0.00000        0.40261
                         0.00023        0.32500       -1.28398

Site:       2         Name:   C         Atomic Number:       6

Coordinates:           -0.077336       1.181820       0.220245

Charge:                 -0.09514
Dipole:                 -0.07365       -0.19097        0.00556
Quadrupole:              0.64297
                        -0.09663        0.49615
                         0.03595        0.36990       -1.13912

Site:       3         Name:   C         Atomic Number:       6

Coordinates:           -1.442045       1.184957       0.229039

Charge:                 -0.09102
Dipole:                  0.05029       -0.16659       -0.02764
Quadrupole:              0.65714
                         0.05311        0.44949
                        -0.00051        0.32888       -1.10664
  
  ...skipping....
  
Local Frame Definition for Multipole Sites :

    Site     Name      Axis Type     Z Axis  X Axis  Y Axis

      1      C         Z-then-X         7       2       0
      2      C         Bisector         1       3       0
      3      C         Z-then-X         4       2       0
      4      C         Z-then-X        10       3       0
      5      C         Z-then-X         4       6       0
      6      C         Bisector         1       5       0
      7      C         Bisector         8       9       0
      8      N         Z-then-X         7      15       0
      9      N         Z-then-X         7      17       0
     10      N         Z-then-X         4      19       0
     11      H         Z-then-X         2       1       0
     12      H         Z-then-X         3       2       0
     13      H         Z-then-X         5       4       0
 ....skipping...
 
 Enter Altered Local Frame Definition [<CR>=Exit] :
  
 Atomic Polarizabilities for Multipole Sites :
    
 Enter Atom Number & Polarizability Values [<CR>=Exit] :
  
 Enter a Bond between Polarization Groups [<CR>=Exit] :  (if the whole mole is one group - recommended for ligands - use default/Enter)
  
 Average the Multipole Moments of Equivalent Atoms [N] :   Y

 Remove Multipole Components Zeroed by Symmetry [N] :  Y
 
Site:       1         Name:   C         Atomic Number:       6

Local Frame:            Z-then-X             7       2       0

Charge:                 -0.10188
Dipole:                 -0.00002        0.00000       -0.02490
Quadrupole:              0.46307
                         0.00000       -1.34443
                         0.00004        0.00000        0.88136

Site:       2         Name:   C         Atomic Number:       6

Local Frame:            Bisector             1       3       0

Charge:                 -0.09514
Dipole:                 -0.02833        0.00000        0.19892
Quadrupole:              0.70208
                         0.00000       -1.21642
                        -0.01601        0.00000        0.51434

Site:       3         Name:   C         Atomic Number:       6

Local Frame:            Z-then-X             4       2       0

Charge:                 -0.09102
Dipole:                  0.13012        0.00000        0.11868
Quadrupole:              0.57375
                         0.00000       -1.17330
                        -0.08666        0.00000        0.59955

...skip...

This will generate benzsp.xyz and benzsp.key

benzsp.key has the electrostatic parameters

parameters ../amoeba.prm
(this comes from tinker.key in the current dir)


atom          1    1    C      "ttt                 "     6     12.011    3
atom          2    2    C      "ttt                 "     6     12.011    3
atom          3    3    C      "ttt                 "     6     12.011    3
atom          4    4    C      "ttt                 "     6     12.011    3
...
multipole     1    7    2             -0.10188
                                     -0.00002    0.00000   -0.02490
                                      0.46307
                                      0.00000   -1.34443
                                      0.00004    0.00000    0.88136
multipole     2    1   -3             -0.09514
                                     -0.02833    0.00000    0.19892
                                      0.70208
                                      0.00000   -1.21642
                                     -0.01601    0.00000    0.51434
...skipping...
polarize      1               1.334      0.390       2    6    7
polarize      2               1.334      0.390       1    3   11
polarize      3               1.334      0.390       2    4   12
polarize      4               1.334      0.390       3    5   10
polarize      5               1.334      0.390       4    6   13
 ...skipping...

Average the multipoles


Due to the symmetry of the ligand, we average the multipoles of these symmetric atoms.

  • Change atom index
  • Average multipoles
  • Shift the index and move it to the parameter file

Introduction to Multipoles

Refit the electrostatics with potential.x

This is useful to optimize the multipoles (after averaging over symmetric atoms,from one basis-set to another...)

Make sure your molecule is always oriented the same from QM to Tinker (always use standard orinetation, all "nosymm" in Gaussian potential calculation)

  1. Create grid points for molecule
  2. Compute electrostatic potential from QM formatted checkpoint file with cubegen
  3. Create potential file for grid points
  4. Fit Electrostatic parameters to target grid.
The TINKER Electrostatic Potential Facility Can :

   (1) Get QM Potential from a Gaussian CUBE File
   (2) Calculate the Model Potential for a System
   (3) Compare the Model Potentials of Two Systems
   (4) Compare a Model Potential to a Target Grid
   (5) Fit Electrostatic Parameters to Target Grid

 Enter the Number of the Desired Choice :  2
 
 Output Potential Value at Each Grid Point [N] :Y
 
  Average Electrostatic Potential over Atoms :
     (Kcal/mole per unit charge)

   Atom         Points        Potential

      1           135          64.2195
      2           371          59.3085
      3           228          54.8976
      4           302          52.9730
      5           235          54.9886
      6           366          59.4935
      7           286          66.5887
      8           523          65.5094
      9           526          65.6609
     10           654          49.3012
     11           635          58.6693
     12          3555          45.6758
     13          3555          45.6982
     14           624          58.6372
     15          1243          62.0622
     16          1637          64.5757
     17          1230          62.0704
     18          1645          64.6416
     19          1071          48.1094
     20          1074          48.0477
  
Electrostatic Potential over all Grid Points :

Average Magnitude for Potential :             54.3511

benzsp.pot and benzsp.grid will be generated.

benzsp.grid has the coordinates of grid

   -1.34113200    -0.00020800    -3.20007300
   -1.38600600     0.19213500    -3.19397200
   -1.61966900     0.01881500    -3.18787100
   -1.53376000    -0.28251900    -3.18177000
   -1.22361200    -0.37674500    -3.17566900
   -0.94649500    -0.19658800    -3.16956800
   -0.87661000     0.13078800    -3.16346700
   -1.04108100     0.42578800    -3.15736600
   -1.35499000     0.55638800    -3.15126500
   -1.68754100     0.47771100    -3.14516400
   -1.92050000     0.22579500    -3.13906300
   -1.25966200    -0.73008500    -3.11465800
   -0.92533300    -0.63616200    -3.10855700
   -0.66777100    -0.40245600    -3.10245600
  ...skipping...
 

benzsp.pot has the potential at each grid points (from tinker)

    19895  
      1      -1.341132   -0.000208   -3.200073       70.7873
      2      -1.386006    0.192135   -3.193972       70.3143
      3      -1.619669    0.018815   -3.187871       72.7352
      4      -1.533760   -0.282519   -3.181770       73.6782
      5      -1.223612   -0.376745   -3.175669       71.7168
      6      -0.946495   -0.196588   -3.169568       69.0251
      7      -0.876610    0.130788   -3.163467       67.5038
      8      -1.041081    0.425788   -3.157366       67.5183
      9      -1.354990    0.556388   -3.151265       68.9467
     10      -1.687541    0.477711   -3.145164       71.3468
     11      -1.920500    0.225795   -3.139063       74.0949
     12      -1.259662   -0.730085   -3.114658       74.0606
   
  ...skipping...
  

/opt/g03/bin/cubegen 0 potential=MP2 benzsp.fchk benzsp.cube -5 h < benzsp.grid

This will compute the potential of the tinker grid from QM

cubegen

-1.341132000000 -0.000208000000 -3.200073000000 0.112800371358
-1.386006000000 0.192135000000 -3.193972000000 0.112181402714
-1.619669000000 0.018815000000 -3.187871000000 0.115556191658
-1.533760000000 -0.282519000000 -3.181770000000 0.116842473623
-1.223612000000 -0.376745000000 -3.175669000000 0.114038580816
-0.946495000000 -0.196588000000 -3.169568000000 0.110351219561
-0.876610000000 0.130788000000 -3.163467000000 0.108378055681
-1.041081000000 0.425788000000 -3.157366000000 0.108460613005
-1.354990000000 0.556388000000 -3.151265000000 0.110406541871
-1.687541000000 0.477711000000 -3.145164000000 0.113700818337
-1.920500000000 0.225795000000 -3.139063000000 0.117550862437
-1.259662000000 -0.730085000000 -3.114658000000 0.117304867056
-0.925333000000 -0.636162000000 -3.108557000000 0.112718952139
-0.667771000000 -0.402456000000 -3.102456000000 0.109223962529
-0.536946000000 -0.079718000000 -3.096355000000 0.107108578011

...skipping...

potential.x to take QM potential from cube file (1) and then fit the parameters to the grid (5)

New potential file converted from cube output which is consistent with tinker pot file in format

    19895  sp energy potential=MP2
      1      -1.341132   -0.000208   -3.200073       70.7833
      2      -1.386006    0.192135   -3.193972       70.3949
      3      -1.619669    0.018815   -3.187871       72.5126
      4      -1.533760   -0.282519   -3.181770       73.3198
      5      -1.223612   -0.376745   -3.175669       71.5603
      6      -0.946495   -0.196588   -3.169568       69.2464
      7      -0.876610    0.130788   -3.163467       68.0083
      8      -1.041081    0.425788   -3.157366       68.0601
      9      -1.354990    0.556388   -3.151265       69.2812
     10      -1.687541    0.477711   -3.145164       71.3483
     11      -1.920500    0.225795   -3.139063       73.7643
     12      -1.259662   -0.730085   -3.114658       73.6099
     13      -0.925333   -0.636162   -3.108557       70.7322
     14      -0.667771   -0.402456   -3.102456       68.5391
     15      -0.536946   -0.079718   -3.096355       67.2116
     16      -1.668816    0.856317   -3.065850       70.4540
   
  ...skipping...
  
  Enter RMS Gradient Termination Criterion [0.001] :  0.1
  
Average Electrostatic Potential over Atoms :
     (Kcal/mole per unit charge)
 
   Atom         Points        Potential        Target        RMS Diff
 
      1           135          64.2195        64.4504         0.3508
      2           371          59.3085        60.0373         0.8938
      3           228          54.8976        56.0802         1.2937
      4           302          52.9730        53.8936         0.9659
      5           235          54.9886        56.0704         1.1716
      6           366          59.4935        60.1774         0.8223
      7           286          66.5887        66.4829         0.2056
      8           523          65.5094        65.3357         0.3199
 ...skipping....
    
 Electrostatic Potential over all Grid Points :

Average Magnitude for Potential :             54.3511
Average Magnitude for Target :                54.4149
Average Signed Potential Difference :         -0.0639
Average Unsigned Potential Difference :        0.5162
Root Mean Square Potential Difference :        0.6756

Potential Fitting of Electrostatic Parameters :

Parameter      Atom Type         Category            Value        Fixed

   --              1             Monopole         -0.05287          X
    1              1             X-Dipole         -0.00001
   --              1             Y-Dipole          0.00000          X
    2              1             Z-Dipole         -0.02333
    3              1             XX-Quad           0.16559
   --              1             XY-Quad           0.00000          X
    4              1             XZ-Quad          -0.00011
    5              1             YY-Quad          -0.51058
   --              1             YZ-Quad           0.00000          X
   --              1             ZZ-Quad           0.34499
   --              2             Monopole         -0.38438          X

  ...skipping...
  Optimally Conditioned Variable Metric Optimization :
 
VM Iter    F Value       G RMS     F Move    X Move      Angle   FG Call
 
    0    2250.6724  20301.4582                                        1
    1    1857.3986  18441.9488   393.2738    0.0002     0.0000        2
    2       0.8428     63.0015  1856.5558    0.0020     0.0000        3
    3       0.8388     62.2482     0.0039    0.0000     0.4472        4
    4       0.6737     61.3801     0.1652    0.0001     0.0001        5
    5       0.6699     57.8278     0.0038    0.0000    12.0351        6
  ...skipping...
     16          1637          63.3915        63.7087         0.3363
     17          1230          61.7354        61.8744         0.1606
     18          1645          63.4573        63.7545         0.3150
     19          1071          48.1338        47.8362         0.3167
     20          1074          48.0967        47.7956         0.3196
 
Electrostatic Potential over all Grid Points :
 
Average Magnitude for Potential :             54.4084
Average Magnitude for Target :                54.4149
Average Signed Potential Difference :         -0.0065
Average Unsigned Potential Difference :        0.1797
Root Mean Square Potential Difference :        0.2223
   

References

  1. Arieh Warshel, Mitsunori Kato, and Andrei V. Pisliakov. JCTC. 2007