Hello, I am trying to calculate RESP charges using RESP plugin. For this, I changed and added PD Radii (as 1.55, got from AMBER ANTECHAMBER dat file) to vdw_surface.py (resp/resp/vdw_surface.py at master · cdsgroup/resp · GitHub) :
vdw_r = {'H': 1.20, 'HE': 1.20, 'LI': 1.37, 'BE': 1.45, 'B': 1.45, 'C': 1.50, 'N': 1.50, 'O': 1.40, 'F': 1.35, 'NE': 1.30, 'NA': 1.57, 'MG': 1.36, 'AL': 1.24, 'SI': 1.17, 'P': 1.80, 'S': 1.75, 'CL': 1.70, 'PD': 1.55 }
Then I tried to calculate resp but how can I change basis set
I have PD atom and I want to calculate RESP charges #HF/CEP-31G level but resp pluginis trying to calculate 6-31G* which not include PD Aatom
so here is the expected error:
BasisSet::construct: Unable to find a basis set for atom 1 for key BASIS among: !
! Shell Entries: [‘PD’] !
! Basis Sets: [(‘6-31G*’, ‘6-31G*’, None)]
NOTE: I also added ‘BASIS_ESP’ : ‘HF-CEP-31G’, (in resp/resp/driver.py at master · cdsgroup/resp · GitHub driver.py line:80 BASIS_ESP was default g-31g so I thought if I added this variable it may change the result):
options = {'VDW_SCALE_FACTORS' : [1.4, 1.6, 1.8, 2.0],
'VDW_POINT_DENSITY' : 1.0,
'RESP_A' : 0.0005,
'RESP_B' : 0.1,
'BASIS_ESP' : 'HF-CEP-31G',
}
here is myinput:
import psi4
import resp
mol = psi4.geometry("""PD 0.0920 -2.2970 -0.3020
CL 1.8320 -3.7270 -0.8090
CL -1.4390 -3.8110 -1.1660
C -2.0400 -1.5570 1.4980
N -1.4200 -0.9740 0.2680
C -0.7060 0.2940 0.6050
C 0.5750 0.0210 1.3360
N 1.4160 -0.9330 0.5500
C 2.3630 -1.5910 1.5180
C -2.4540 -0.7180 -0.8120
C -3.4900 0.3390 -0.4610
C -4.6910 -0.0140 0.2120
C -5.0770 -1.3720 0.4840
C -6.2240 -1.6670 1.1560
C -7.0790 -0.6440 1.6300
C -6.7960 0.6350 1.3630
C -5.6180 1.0160 0.6440
C -5.3240 2.3270 0.3760
C -4.1890 2.7200 -0.3410
C -3.9260 4.0850 -0.6320
C -2.8420 4.4660 -1.3430
C -1.9720 3.4780 -1.8580
C -2.1760 2.1500 -1.6200
C -3.2650 1.6920 -0.8030
C 2.1830 -0.2720 -0.5930
C 3.2790 0.6910 -0.1920
C 3.0330 2.1000 -0.0660
C 1.7530 2.6970 -0.3790
C 1.5650 4.0480 -0.2050
C 2.5950 4.8730 0.2650
C 3.8120 4.3580 0.5480
C 4.0730 2.9610 0.3860
C 5.3280 2.4470 0.6570
C 5.6250 1.0910 0.4710
C 6.9270 0.5620 0.7320
C 7.2430 -0.7330 0.5260
C 6.2600 -1.5940 0.0290
C 4.9840 -1.1600 -0.2210
C 4.5910 0.2090 0.0120
H -2.4890 -2.3600 1.2860
H -1.3460 -1.7280 2.1500
H -2.6590 -0.9300 1.8720
H -0.4950 0.7700 -0.2020
H -1.2700 0.8350 1.1570
H 0.3630 -0.3620 2.1990
H 1.0490 0.8380 1.4570
H 2.9240 -2.2040 1.0340
H 2.9220 -0.9180 1.9360
H 1.8660 -2.0580 2.1810
H -2.0020 -0.4430 -1.6070
H -2.9170 -1.5400 -0.9790
H -4.5070 -2.0700 0.1800
H -6.4520 -2.5860 1.3250
H -7.8650 -0.8580 2.1400
H -7.3980 1.3260 1.6680
H -5.9060 3.0130 0.7120
H -4.5360 4.7540 -0.3090
H -2.6760 5.3980 -1.5090
H -1.2190 3.7440 -2.3880
H -1.5780 1.4980 -2.0030
H 1.5500 0.2110 -1.1210
H 2.5700 -0.9700 -1.1280
H 1.0360 2.1390 -0.6980
H 0.7090 4.4320 -0.4090
H 2.4420 5.8070 0.3970
H 4.5070 4.9270 0.8430
H 6.0100 3.0410 0.9660
H 7.6090 1.1620 1.0700
H 8.1270 -1.0550 0.7150
H 6.4950 -2.5060 -0.1410
H 4.3390 -1.7780 -0.5730 """)
mol.update_geometry()
options = {'VDW_SCALE_FACTORS' : [1.4, 1.6, 1.8, 2.0],
'VDW_POINT_DENSITY' : 1.0,
'RESP_A' : 0.0005,
'RESP_B' : 0.1,
'BASIS_ESP' : 'HF-CEP-31G',
}
# Call for first stage fit
charges1 = resp.resp([mol], options)
print('Electrostatic Potential Charges')
print(charges1[0])
print('Restrained Electrostatic Potential Charges')
print(charges1[1])
# Change the value of the RESP parameter A
options['RESP_A'] = 0.001
# Add constraint for atoms fixed in second stage fit
constraint_charge = []
for i in range(4, 8):
constraint_charge.append([charges1[1][i], [i+1]])
options['constraint_charge'] = constraint_charge
options['constraint_group'] = [[2, 3, 4]]
options['grid'] = ['1_%s_grid.dat' %mol.name()]
options['esp'] = ['1_%s_grid_esp.dat' %mol.name()]
# Call for second stage fit
charges2 = resp.resp([mol], options)
# Get RESP charges
print("\nStage Two:\n")
print('RESP Charges')
print(charges2[1])