Reparameterize force fields#
The MMForceFieldGenerator class in VeloxChem provides a number of functions to manipulate topology files. As an example, we consider the optimization of the force field parameters for a dihedral angle of a thiophene-based optical ligand named HS-276.
Load the required modules and the equilibrium geometry:
import veloxchem as vlx
import numpy as np
# load B3LYP optimized geometry in the xyz format
molecule = vlx.Molecule.read_xyz_file("../../data/md/hs276.xyz")
molecule.show(atom_indices=True, width=600, height=450)
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
We here consider the dihedral angle between the two thiophene rings, which corresponds to the rotation around bond 21-22.
The dihedral potential curves will be plotted using force field parameters from the general amber forcefield (GAFF) database, as compared to QM results. The force field will then be changed to better fit this potential.
Create the initial force field#
Initialize an MMForceFieldGenerator object and generate initial force field with pre-calculated RESP charges:
resp_chg = np.array([
0.014967, 0.224409, 0.089489, -0.551618, 0.153050, 0.106847, -0.315042,
0.165854, 0.206441, -0.425163, 0.225816, -0.039099, 0.179377, 0.002119,
0.014483, -0.043666, -0.230202, 0.165960, -0.071404, 0.188598, -0.044661,
0.088786, -0.173586, 0.175730, -0.089376, 0.164515, -0.118573, -0.150541,
0.736497, -0.570770, -0.273724, -0.144051, 0.112846, 0.112846, 0.112846
])
ff_gen = vlx.MMForceFieldGenerator()
ff_gen.partial_charges = resp_chg
ff_gen.create_topology(molecule)
* Info * Using GAFF (v2.11) parameters.
Reference: J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A. Case, J. Comput. Chem. 2004,
25, 1157-1174.
* Info * Updated bond angle 1-9-10 (ca-ca-cc) to 107.217 deg
* Info * Updated bond angle 1-14-15 (ca-na-cc) to 126.641 deg
* Info * Updated bond angle 2-1-14 (ca-ca-na) to 131.687 deg
* Info * Updated bond angle 7-9-10 (ca-ca-cc) to 135.297 deg
* Info * Updated bond angle 9-1-14 (ca-ca-na) to 107.608 deg
* Info * Updated bond angle 14-15-17 (na-cc-cd) to 128.324 deg
* Info * Updated bond angle 19-21-22 (cd-cc-cc) to 129.220 deg
* Info * Updated bond angle 21-22-23 (cc-cc-cd) to 128.571 deg
Reparameterize the force field#
Use the reparameterize_dihedrals method to fit dihedral parameters to QM scan. If you would like to see more output during the fitting, set verbose to True.
ff_gen.reparameterize_dihedrals(rotatable_bond=(21,22),
scan_file="../../data/md/27-22-21-16.xyz",
visualize=True,
verbose=False)
VeloxChem Dihedral Reparameterization
=======================================
* Info * Rotatable bond selected: 21-22
* Info * Dihedrals involved:[[16, 21, 22, 23], [16, 21, 22, 27], [19, 21, 22, 23], [19, 21, 22, 27]]
* Info * Reading QM scan from file...
* Info * ../../data/md/27-22-21-16.xyz
* Info * Performing dihedral scan for MM baseline by excluding the involved dihedral barriers...
* Info * Target dihedral angle: 16-21-22-27
* Info * Dihedral barriers [4.184 4.184 4.184 4.184] will be used as initial guess.
* Info * Validating the initial force field...
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 0.0 deg 7.215 kJ/mol 3.836 kJ/mol 3.379
* Info * 10.0 deg 7.767 kJ/mol 3.392 kJ/mol 4.375
* Info * 20.0 deg 9.607 kJ/mol 3.074 kJ/mol 6.532
* Info * 30.0 deg 12.706 kJ/mol 3.059 kJ/mol 9.648
* Info * 40.0 deg 16.790 kJ/mol 3.796 kJ/mol 12.994
* Info * 50.0 deg 21.360 kJ/mol 5.416 kJ/mol 15.944
* Info * 60.0 deg 25.786 kJ/mol 7.761 kJ/mol 18.025
* Info * 70.0 deg 29.413 kJ/mol 10.365 kJ/mol 19.047
* Info * 80.0 deg 31.673 kJ/mol 12.461 kJ/mol 19.212
* Info * 90.0 deg 32.177 kJ/mol 13.214 kJ/mol 18.963
* Info * 100.0 deg 30.777 kJ/mol 12.264 kJ/mol 18.513
* Info * 110.0 deg 27.591 kJ/mol 9.969 kJ/mol 17.622
* Info * 120.0 deg 22.990 kJ/mol 7.136 kJ/mol 15.854
* Info * 130.0 deg 17.550 kJ/mol 4.450 kJ/mol 13.100
* Info * 140.0 deg 11.971 kJ/mol 2.316 kJ/mol 9.656
* Info * 150.0 deg 6.959 kJ/mol 0.943 kJ/mol 6.016
* Info * 160.0 deg 3.092 kJ/mol 0.294 kJ/mol 2.798
* Info * 170.0 deg 0.728 kJ/mol 0.087 kJ/mol 0.641
* Info * 180.0 deg 0.000 kJ/mol 0.000 kJ/mol 0.000
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 19.212 kJ/mol
* Info * Standard deviation: 6.677 kJ/mol
* Info * Fitting the dihedral parameters...
* Info * Optimizing dihedral via least squares fitting...
* Info * New fitted barriers: [1.85973955 1.85973955 1.85973955 1.85973955]
* Info * Validating the fitted force field...
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 0.0 deg 7.215 kJ/mol 3.836 kJ/mol 3.379
* Info * 10.0 deg 7.206 kJ/mol 3.392 kJ/mol 3.814
* Info * 20.0 deg 7.432 kJ/mol 3.074 kJ/mol 4.357
* Info * 30.0 deg 8.058 kJ/mol 3.059 kJ/mol 4.999
* Info * 40.0 deg 9.107 kJ/mol 3.796 kJ/mol 5.311
* Info * 50.0 deg 10.449 kJ/mol 5.416 kJ/mol 5.032
* Info * 60.0 deg 11.841 kJ/mol 7.761 kJ/mol 4.080
* Info * 70.0 deg 12.994 kJ/mol 10.365 kJ/mol 2.628
* Info * 80.0 deg 13.639 kJ/mol 12.461 kJ/mol 1.179
* Info * 90.0 deg 13.583 kJ/mol 13.214 kJ/mol 0.369
* Info * 100.0 deg 12.743 kJ/mol 12.264 kJ/mol 0.480
* Info * 110.0 deg 11.172 kJ/mol 9.969 kJ/mol 1.203
* Info * 120.0 deg 9.044 kJ/mol 7.136 kJ/mol 1.908
* Info * 130.0 deg 6.639 kJ/mol 4.450 kJ/mol 2.188
* Info * 140.0 deg 4.289 kJ/mol 2.316 kJ/mol 1.973
* Info * 150.0 deg 2.310 kJ/mol 0.943 kJ/mol 1.368
* Info * 160.0 deg 0.917 kJ/mol 0.294 kJ/mol 0.623
* Info * 170.0 deg 0.167 kJ/mol 0.087 kJ/mol 0.081
* Info * 180.0 deg 0.000 kJ/mol 0.000 kJ/mol 0.000
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 5.311 kJ/mol
* Info * Standard deviation: 1.753 kJ/mol
* Info * Dihedral MM parameters have been reparameterized and updated in the topology.
It is not uncommon that the existing dihedral terms in the force field are not sufficient to reproduce the QM potential. Here we can try adding more dihedral terms to the force field, and then redo the fitting. In this particular case we would like to lower the MM energy at 0\(^{\circ}\) so that it matches the QM energy better. This corresponds to adding a dihedral term with periodicity 1 and phase 180\(^{\circ}\). The barrier can be set to 0 since it will be fitted.
ff_gen.add_dihedral((16, 21, 22, 27), barrier=0.0, phase=180.0, periodicity=1)
* Info * Added dihedral 16-21-22-27
After adding the dihedral term we can redo the fitting. This time we can skip the initial validation step by setting initial_validation to False, since we already know how the MM potential energy curve looks. We can also set show_diff to True so that the difference between QM and MM energies will be plotted after the fitting.
ff_gen.reparameterize_dihedrals(rotatable_bond=(21,22),
scan_file="../../data/md/27-22-21-16.xyz",
visualize=True,
verbose=False,
initial_validation=False,
show_diff=True)
VeloxChem Dihedral Reparameterization
=======================================
* Info * Rotatable bond selected: 21-22
* Info * Dihedrals involved:[[16, 21, 22, 23], [16, 21, 22, 27], [19, 21, 22, 23], [19, 21, 22, 27]]
* Info * Reading QM scan from file...
* Info * ../../data/md/27-22-21-16.xyz
* Info * Performing dihedral scan for MM baseline by excluding the involved dihedral barriers...
* Info * Target dihedral angle: 16-21-22-27
* Info * Dihedral barriers [1.85973955 1.85973955 0. 1.85973955 1.85973955] will be used as initial guess.
* Info * Fitting the dihedral parameters...
* Info * Optimizing dihedral via least squares fitting...
* Info * New fitted barriers: [1.9819721 1.9819721 1.97102933 1.9819721 1.9819721 ]
* Info * Validating the fitted force field...
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 0.0 deg 3.273 kJ/mol 3.836 kJ/mol -0.563
* Info * 10.0 deg 3.323 kJ/mol 3.392 kJ/mol -0.069
* Info * 20.0 deg 3.723 kJ/mol 3.074 kJ/mol 0.649
* Info * 30.0 deg 4.624 kJ/mol 3.059 kJ/mol 1.566
* Info * 40.0 deg 6.031 kJ/mol 3.796 kJ/mol 2.234
* Info * 50.0 deg 7.785 kJ/mol 5.416 kJ/mol 2.368
* Info * 60.0 deg 9.617 kJ/mol 7.761 kJ/mol 1.856
* Info * 70.0 deg 11.212 kJ/mol 10.365 kJ/mol 0.847
* Info * 80.0 deg 12.275 kJ/mol 12.461 kJ/mol -0.186
* Info * 90.0 deg 12.590 kJ/mol 13.214 kJ/mol -0.625
* Info * 100.0 deg 12.063 kJ/mol 12.264 kJ/mol -0.201
* Info * 110.0 deg 10.738 kJ/mol 9.969 kJ/mol 0.769
* Info * 120.0 deg 8.792 kJ/mol 7.136 kJ/mol 1.656
* Info * 130.0 deg 6.508 kJ/mol 4.450 kJ/mol 2.058
* Info * 140.0 deg 4.232 kJ/mol 2.316 kJ/mol 1.916
* Info * 150.0 deg 2.291 kJ/mol 0.943 kJ/mol 1.348
* Info * 160.0 deg 0.912 kJ/mol 0.294 kJ/mol 0.618
* Info * 170.0 deg 0.167 kJ/mol 0.087 kJ/mol 0.080
* Info * 180.0 deg 0.000 kJ/mol 0.000 kJ/mol 0.000
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 2.368 kJ/mol
* Info * Standard deviation: 0.969 kJ/mol
* Info * Dihedral MM parameters have been reparameterized and updated in the topology.
Now we get better agreement between QM and MM energies. Since the difference between QM and MM is also plotted, we can tell that the MM energies can be brought closer to QM energies by adding one more dihedral term with periodicity 4 and phase 0\(^{\circ}\).
ff_gen.add_dihedral((16, 21, 22, 27), barrier=0.0, phase=0.0, periodicity=4)
* Info * Added dihedral 16-21-22-27
ff_gen.reparameterize_dihedrals(rotatable_bond=(21,22),
scan_file="../../data/md/27-22-21-16.xyz",
visualize=True,
verbose=False,
initial_validation=False)
VeloxChem Dihedral Reparameterization
=======================================
* Info * Rotatable bond selected: 21-22
* Info * Dihedrals involved:[[16, 21, 22, 23], [16, 21, 22, 27], [19, 21, 22, 23], [19, 21, 22, 27]]
* Info * Reading QM scan from file...
* Info * ../../data/md/27-22-21-16.xyz
* Info * Performing dihedral scan for MM baseline by excluding the involved dihedral barriers...
* Info * Target dihedral angle: 16-21-22-27
* Info * Dihedral barriers [1.9819721 1.9819721 1.97102933 0. 1.9819721 1.9819721 ] will be used as initial guess.
* Info * Fitting the dihedral parameters...
* Info * Optimizing dihedral via least squares fitting...
* Info * New fitted barriers: [2.00263406 2.00263406 1.90579111 1.23978221 2.00263406 2.00263406]
* Info * Validating the fitted force field...
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 0.0 deg 3.521 kJ/mol 3.836 kJ/mol -0.315
* Info * 10.0 deg 3.285 kJ/mol 3.392 kJ/mol -0.107
* Info * 20.0 deg 2.962 kJ/mol 3.074 kJ/mol -0.113
* Info * 30.0 deg 3.045 kJ/mol 3.059 kJ/mol -0.014
* Info * 40.0 deg 3.926 kJ/mol 3.796 kJ/mol 0.130
* Info * 50.0 deg 5.701 kJ/mol 5.416 kJ/mol 0.285
* Info * 60.0 deg 8.097 kJ/mol 7.761 kJ/mol 0.336
* Info * 70.0 deg 10.539 kJ/mol 10.365 kJ/mol 0.173
* Info * 80.0 deg 12.339 kJ/mol 12.461 kJ/mol -0.122
* Info * 90.0 deg 12.937 kJ/mol 13.214 kJ/mol -0.277
* Info * 100.0 deg 12.105 kJ/mol 12.264 kJ/mol -0.159
* Info * 110.0 deg 10.020 kJ/mol 9.969 kJ/mol 0.051
* Info * 120.0 deg 7.206 kJ/mol 7.136 kJ/mol 0.070
* Info * 130.0 deg 4.341 kJ/mol 4.450 kJ/mol -0.109
* Info * 140.0 deg 2.028 kJ/mol 2.316 kJ/mol -0.288
* Info * 150.0 deg 0.598 kJ/mol 0.943 kJ/mol -0.344
* Info * 160.0 deg 0.028 kJ/mol 0.294 kJ/mol -0.266
* Info * 170.0 deg 0.000 kJ/mol 0.087 kJ/mol -0.087
* Info * 180.0 deg 0.117 kJ/mol 0.000 kJ/mol 0.117
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 0.344 kJ/mol
* Info * Standard deviation: 0.197 kJ/mol
* Info * Dihedral MM parameters have been reparameterized and updated in the topology.
Save the force field to file#
Use write_gromacs_files or write_openmm_files to write the force field to files.
ff_gen.write_gromacs_files('hs276')