Very different converged charges with different basis sets. about resp

[alenaizan]

Thanks Peter for this report.
First, you are using an old version of the code, likely from Psi4Numpy. The code that you are using did not correctly handle geometries defined in Bohr. Update the geometry to Angstrom and you will get different answers. Please update the code to the current master in this repository. Once you update the code, use the following script to generate the charges:

import psi4
from resp import resp
import numpy as np

mol = psi4.geometry('''
    O  0.000000000000 -0.143225816552 0.000000000000
    H  1.638036840407 1.136548822547 -0.000000000000
    H  -1.638036840407 1.136548822547 -0.000000000000
    noreorient
    units Bohr
''')

mol.update_geometry()

options = {
           'VDW_SCALE_FACTORS'  : [1.0],
           'VDW_POINT_DENSITY'  : 1.0,
           'CONSTRAINT_GROUP'   : [[2,3]], 
           'BASIS_ESP'          : '6-31G*' ,
}

charges1 = resp([mol], options)

options["BASIS_ESP"] = "STO-3G"
charges2 = resp([mol], options)

options["BASIS_ESP"] = "aug-cc-pVDZ"
charges3 = resp([mol], options)

options["BASIS_ESP"] = "DZ"
charges4 = resp([mol], options)

print('Basis, ESP, RESP')
print("6-31G*      ", charges1[0], charges1[1])
print("STO-3G      ", charges2[0], charges2[1])
print("aug-cc-pVDZ ", charges3[0], charges3[1])
print("DZ          ", charges4[0], charges4[1])

You will get the following output:

Basis, ESP, RESP
6-31G*       [-0.68816353  0.34408177  0.34408177] [-0.68547359  0.3427368   0.3427368 ]
STO-3G       [-0.48345216  0.24172608  0.24172608] [-0.48079071  0.24039535  0.24039535]
aug-cc-pVDZ  [-0.61675067  0.30837533  0.30837533] [-0.6140676  0.3070338  0.3070338]
DZ           [-0.78604978  0.39302489  0.39302489] [-0.78335328  0.39167664  0.39167664]

These values are better but they still do not agree with the data that you generated with other codes.

There are two possible sources of errors: (1) the fitting procedure, and (2) the grid generation. To test (1), I used the original RESP implementation by Chris Bayly (available in ambertools) with my grid. The results are the following:

6-31G*  -0.690039  0.345020  0.345020
STO-3G  -0.484770  0.242385  0.242385
aug-cc-pVDZ  -0.618432  0.309216  0.309216
DZ  -0.788192  0.394096  0.394096

As you can see, they are very similar to my values.

Next, I tested the grid. Please note that your options for the grid are not recommended. Please use the following values:

'VDW_SCALE_FACTORS'  : [1.4, 1.6, 1.8, 2.0],

This ensures that grid points are outside the van der Waals surface of the molecule and that multiple layers are used. When I do this, I get grid points that are identical with the Connolly surface algorithm available in GAMESS and generated using the following input file:

 $CONTRL ICHARG=0 MULT=1 RUNTYP=ENERGY MOLPLT=.T.
         MPLEVL=0 UNITS=BOHR MAXIT=200 EXETYP=RUN 
         SCFTYP=RHF 
         COORD=UNIQUE                                  $END
 $SCF    DIRSCF=.T. CONV=1.0E-06                       $END
 $SYSTEM TIMLIM=5000 MWORDS=32 MEMDDI=0                $END
 $BASIS  GBASIS=N31 NGAUSS=6 NDFUNC=1                  $END
 $GUESS  GUESS=HUCKEL                                  $END
 $ELPOT  IEPOT=1 WHERE=PDC OUTPUT=BOTH                 $END 
 $PDC    PTSEL=CONNOLLY CONSTR=NONE                    $END
 $DATA
 MOLECULE
 C1
 O           8   0.000000000000  -0.143225816552   0.00000000
 H           1   1.638036840407   1.136548822547   0.00000000
 H           1  -1.638036840407   1.136548822547   0.00000000
 $END

So in summary, I think my code is correct. If you want to debug the issue further, please attach files containing the grids and ESP potentials generated by the other codes.

[pwborthwick]

Hi @alenaizan ,
Thank you. The code I used was from psi4 conda download version 1.3.2. I won't persue this further but note that the sto-3g result (even though it's a trivial basis and I was only using one shell - for testing) still doesn't agree with https://github.com/erikkjellgren/SlowQuant/blob/master/data/testfiles/outQFIT.txt or my code. (I also get results closer to Slowquant with my own code).
But I do see the problem with units eg

    data['invr'].append(invr*bohr_to_angstrom) # convert to atomic units
    data['coordinates'][-1] /= bohr_to_angstrom # convert to angstroms

invr is surely in angstroms being derived from points and coordinates, both in angstroms. So if your comment is right you should be dividing by 'bohr->angstrom'. I thought the coordinates were in angstroms at this point and anyway dividing by bohr->angstrom is converting them to bohr.
But if you're happy then I'm happy and thanks for your time.
Peter

[alenaizan]

Peter, I have a couple of points. First, the SlowQuant code you referred to seems to constrain the dipoles, which I don't do in my code. You should disable this if you want to compare my code to SlowQuant. Second, the fitting procedure is sensitive to the choice of grid points. I see that SlowQuant uses a different set of van der Waals radii and a different algorithm for generating the grid points. The first point is more important, though.