Example 2: RESP charge fitting of aniline¶
Fitting of charges¶
In order to generate the electrostatic potential of aniline, first we will optimise the geometry of aniline at the B3LYP/6-31G* level. This is the usual level of theory used for small molecules and transition-metal complexes accurate calculations), whereas calculations at the AM1 level of theory are usually performed for large organic molecules (relatively crude calculations). We will use Gaussian to perform the optimization. The Gaussian input file reads:
Geometry optimization Gaussian input file.¶
%chk=aniline.chk
# b3lyp/6-31g* opt=tight
Title Card Required
0 1
C -0.53956834 1.85851314 0.00000000
C 0.85559166 1.85851314 0.00000000
C 1.55312966 3.06626414 0.00000000
C 0.85547566 4.27477314 -0.00119900
C -0.53934934 4.27469514 -0.00167800
C -1.23695034 3.06648914 -0.00068200
H -1.08932734 0.90619614 0.00045000
H 1.40509966 0.90600014 0.00131500
H 1.40567566 5.22691614 -0.00125800
H -1.08947134 5.22697614 -0.00263100
H -2.33655434 3.06667214 -0.00086200
N 3.02312941 3.06637108 0.00084750
H 3.35602581 2.59534128 0.81773851
H 3.35696733 2.59464054 -0.81525454
Furthermore, in order too calculate the electrostatic potential, we will run a Gaussian single-point energy calculation, at the HF/6-31G* level (without transition metals) or at the B3LYP/DZpdf/6-31G* level. The following keywods have to be used:
IOp(6/33=2) makes Gaussian write out the potential points and potentials (do not change)
IOp(6/41=10) specifies that 10 concentric layers of points are used for each atom (do not change)
IOp(6/42) gives the density of points in each layer. A value of 17 gives about 2500 points/atom. Lower values may be needed for large molecules, since the programs cannot normally handle more than 100 000 potential points. A value of 10 gives about 1000 points/atom.
The Gaussian input file for this calculation reads:
ESP calculation Gaussian input file.¶
%chk=aniline_opt.chk
# HF/6-31G* SCF=Tight Pop=MK IOp(6/33=2,6/41=10,6/42=17)
Title Card Required
0 1
C 1.17182300 1.20260700 -0.00337800
C -0.22125800 1.20825000 0.00517500
C -0.93874400 0.00000500 0.01013900
C -0.22126400 -1.20824500 0.00517400
C 1.17181400 -1.20261300 -0.00337600
C 1.88169200 -0.00000400 -0.00833000
H 1.70568500 2.14975100 -0.00879700
H -0.76331600 2.15175300 0.01327800
H -0.76333600 -2.15174000 0.01327300
H 1.70567000 -2.14976000 -0.00879700
H 2.96765300 -0.00000700 -0.01660800
N -2.33738000 0.00000100 0.07891000
H -2.77753500 0.83490500 -0.28857500
H -2.77753600 -0.83490600 -0.28856700
Finally, in order to read the ESP data into ParaMol and perform RESP charge fitting, the following
# ParaMol imports
from ParaMol.System.system import *
from ParaMol.MM_engines.openmm import *
# ParaMol Tasks imports
from ParaMol.Tasks.resp_fitting import *
from ParaMol.Utils.settings import *
from ParaMol.Utils.gaussian_esp import *
from ParaMol.Utils.Symmetrizers.amber_symmetrizer import *
# --------------------------------------------------------- #
# Preparation #
# --------------------------------------------------------- #
# Create the OpenMM engine for caffeine
openmm_system = OpenMMEngine(init_openmm=True, topology_format='AMBER', top_file='aniline.prmtop', crd_format='AMBER', crd_file='aniline.inpcrd')
# Create Molecular System
aniline = ParaMolSystem(name="aniline", engine=openmm_system, n_atoms=14)
# Create ParaMol's force field representation and ask to optimize charges
aniline.force_field.create_force_field(opt_charges=True)
# Create ParaMol settings instance
paramol_settings = Settings()
paramol_settings.properties["include_energies"] = False
paramol_settings.properties["include_forces"] = False
paramol_settings.properties["include_esp"] = True
# --------------------------------------------------------- #
# Read ESP Data into ParaMol #
# --------------------------------------------------------- #
gaussian_esp = GaussianESP()
aniline.ref_coordinates, aniline.ref_esp_grid, aniline.ref_esp = gaussian_esp.read_log_files(["path_to_guassian_log_file"])
# --------------------------------------------------------- #
# Symmetrize ParaMol ForceField #
# --------------------------------------------------------- #
# Symmetry ParaMol ForceField so that it respects symmetries
# In this example, we are not setting any symmetry, but we still need to do this step as we want to save a .mol2 file
amber_symmetrizer = AmberSymmetrizer(top_file="aniline.prmtop")
amber_symmetrizer.get_symmetries(aniline.force_field)
amber_symmetrizer.symmetrize_force_field(aniline.force_field)
# Set number of structures
aniline.n_structures = len(aniline.ref_coordinates)
# --------------------------------------------------------- #
# RESP Charge Fitting #
# --------------------------------------------------------- #
resp_fitting = RESPFitting()
systems, parameter_space, objective_function, optimizer = resp_fitting.run_task(paramol_settings, [aniline], solver="scipy", total_charge=0)
#systems, parameter_space, objective_function, optimizer = resp_fitting.run_task(paramol_settings, [aniline], solver="explicit", total_charge=0)
# Write ParaMol Force Field file with final parameters
aniline.force_field.write_ff_file("aniline_resp.ff")
# Update amber symmetrizer and save .mol2 file
amber_symmetrizer.update_term_types_parameters(parameter_space.optimizable_parameters)
amber_symmetrizer.save("aniline_resp.mol2")
Comparision of the ESP charges¶
We may now check the results obtained using the solvers available in ParaMol (‘scipy’ and ‘explicit’).
ESP charges obtained when using the ‘SLSQP’ SciPy optimizer (‘scipy’ solver).¶
NonbondedForce 3
0 0 -0.06152112 0.33996695 0.35982400 1 0 0 X
1 1 -0.37029657 0.33996695 0.35982400 1 0 0 X
2 2 0.50425830 0.33996695 0.35982400 1 0 0 X
3 3 -0.37030046 0.33996695 0.35982400 1 0 0 X
4 4 -0.06152352 0.33996695 0.35982400 1 0 0 X
5 5 -0.24057025 0.33996695 0.35982400 1 0 0 X
6 6 0.13601036 0.25996425 0.06276000 1 0 0 X
7 7 0.18160558 0.25996425 0.06276000 1 0 0 X
8 8 0.18160579 0.25996425 0.06276000 1 0 0 X
9 9 0.13601494 0.25996425 0.06276000 1 0 0 X
10 10 0.14580236 0.25996425 0.06276000 1 0 0 X
11 11 -0.92021579 0.32499985 0.71128000 1 0 0 X
12 12 0.36956553 0.10690785 0.06568880 1 0 0 X
13 13 0.36956488 0.10690785 0.06568880 1 0 0 X
ESP charges obtained when using the ‘explicit’ solver.¶
NonbondedForce 3
0 0 -0.06286283 0.33996695 0.35982400 1 0 0 X
1 1 -0.36913364 0.33996695 0.35982400 1 0 0 X
2 2 0.50319411 0.33996695 0.35982400 1 0 0 X
3 3 -0.36913297 0.33996695 0.35982400 1 0 0 X
4 4 -0.06285939 0.33996695 0.35982400 1 0 0 X
5 5 -0.23913232 0.33996695 0.35982400 1 0 0 X
6 6 0.13628034 0.25996425 0.06276000 1 0 0 X
7 7 0.18140134 0.25996425 0.06276000 1 0 0 X
8 8 0.18140069 0.25996425 0.06276000 1 0 0 X
9 9 0.13627925 0.25996425 0.06276000 1 0 0 X
10 10 0.14549336 0.25996425 0.06276000 1 0 0 X
11 11 -0.92005714 0.32499985 0.71128000 1 0 0 X
12 12 0.36956380 0.10690785 0.06568880 1 0 0 X
13 13 0.36956540 0.10690785 0.06568880 1 0 0 X
ESP charges as calculated by Gaussian.¶
Fitting point charges to electrostatic potential
Charges from ESP fit, RMS= 0.00123 RRMS= 0.10456:
ESP charges:
1
1 C -0.062863
2 C -0.369134
3 C 0.503195
4 C -0.369134
5 C -0.062859
6 C -0.239133
7 H 0.136280
8 H 0.181401
9 H 0.181401
10 H 0.136279
11 H 0.145494
12 N -0.920057
13 H 0.369564
14 H 0.369565