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:
ff_gen.reparameterize_dihedrals(rotatable_bond=(21,22), scan_file="../../data/md/27-22-21-16.xyz", visualize=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 MM energy
* Info * --------------------------------
* Info * 0.0 deg 38.292 kJ/mol
* Info * 10.0 deg 37.834 kJ/mol
* Info * 20.0 deg 36.768 kJ/mol
* Info * 30.0 deg 35.415 kJ/mol
* Info * 40.0 deg 34.037 kJ/mol
* Info * 50.0 deg 32.794 kJ/mol
* Info * 60.0 deg 31.758 kJ/mol
* Info * 70.0 deg 30.933 kJ/mol
* Info * 80.0 deg 30.287 kJ/mol
* Info * 90.0 deg 29.781 kJ/mol
* Info * 100.0 deg 29.391 kJ/mol
* Info * 110.0 deg 29.110 kJ/mol
* Info * 120.0 deg 28.962 kJ/mol
* Info * 130.0 deg 28.984 kJ/mol
* Info * 140.0 deg 29.218 kJ/mol
* Info * 150.0 deg 29.667 kJ/mol
* Info * 160.0 deg 30.253 kJ/mol
* Info * 170.0 deg 30.795 kJ/mol
* Info * 180.0 deg 31.076 kJ/mol
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 0.0 deg 9.329 kJ/mol 3.836 kJ/mol 5.494
* Info * 10.0 deg 8.872 kJ/mol 3.392 kJ/mol 5.480
* Info * 20.0 deg 7.806 kJ/mol 3.074 kJ/mol 4.731
* Info * 30.0 deg 6.453 kJ/mol 3.059 kJ/mol 3.394
* Info * 40.0 deg 5.074 kJ/mol 3.796 kJ/mol 1.278
* Info * 50.0 deg 3.832 kJ/mol 5.416 kJ/mol -1.584
* Info * 60.0 deg 2.796 kJ/mol 7.761 kJ/mol -4.965
* Info * 70.0 deg 1.971 kJ/mol 10.365 kJ/mol -8.395
* Info * 80.0 deg 1.324 kJ/mol 12.461 kJ/mol -11.136
* Info * 90.0 deg 0.819 kJ/mol 13.214 kJ/mol -12.395
* Info * 100.0 deg 0.428 kJ/mol 12.264 kJ/mol -11.835
* Info * 110.0 deg 0.148 kJ/mol 9.969 kJ/mol -9.821
* Info * 120.0 deg 0.000 kJ/mol 7.136 kJ/mol -7.136
* Info * 130.0 deg 0.022 kJ/mol 4.450 kJ/mol -4.428
* Info * 140.0 deg 0.256 kJ/mol 2.316 kJ/mol -2.060
* Info * 150.0 deg 0.705 kJ/mol 0.943 kJ/mol -0.238
* Info * 160.0 deg 1.291 kJ/mol 0.294 kJ/mol 0.997
* Info * 170.0 deg 1.833 kJ/mol 0.087 kJ/mol 1.746
* Info * 180.0 deg 2.114 kJ/mol 0.000 kJ/mol 2.114
* Info * Dihedral barriers [4.184 4.184 4.184 4.184] will be used as initial guess.
* Info * Validating the initial force field...
* 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 * 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. We can also skip the initial validation step by setting initial_validation to False.
ff_gen.add_dihedral((16, 21, 22, 27), barrier=0.0, phase=180.0, periodicity=1)
ff_gen.add_dihedral((16, 21, 22, 27), barrier=0.0, phase=0.0, periodicity=4)
* Info * Added dihedral 16-21-22-27
* 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, 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 MM energy
* Info * --------------------------------
* Info * 0.0 deg 38.292 kJ/mol
* Info * 10.0 deg 37.834 kJ/mol
* Info * 20.0 deg 36.768 kJ/mol
* Info * 30.0 deg 35.415 kJ/mol
* Info * 40.0 deg 34.037 kJ/mol
* Info * 50.0 deg 32.794 kJ/mol
* Info * 60.0 deg 31.758 kJ/mol
* Info * 70.0 deg 30.933 kJ/mol
* Info * 80.0 deg 30.287 kJ/mol
* Info * 90.0 deg 29.781 kJ/mol
* Info * 100.0 deg 29.391 kJ/mol
* Info * 110.0 deg 29.110 kJ/mol
* Info * 120.0 deg 28.962 kJ/mol
* Info * 130.0 deg 28.984 kJ/mol
* Info * 140.0 deg 29.218 kJ/mol
* Info * 150.0 deg 29.667 kJ/mol
* Info * 160.0 deg 30.253 kJ/mol
* Info * 170.0 deg 30.795 kJ/mol
* Info * 180.0 deg 31.076 kJ/mol
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 0.0 deg 9.329 kJ/mol 3.836 kJ/mol 5.494
* Info * 10.0 deg 8.872 kJ/mol 3.392 kJ/mol 5.480
* Info * 20.0 deg 7.806 kJ/mol 3.074 kJ/mol 4.731
* Info * 30.0 deg 6.453 kJ/mol 3.059 kJ/mol 3.394
* Info * 40.0 deg 5.074 kJ/mol 3.796 kJ/mol 1.278
* Info * 50.0 deg 3.832 kJ/mol 5.416 kJ/mol -1.584
* Info * 60.0 deg 2.796 kJ/mol 7.761 kJ/mol -4.965
* Info * 70.0 deg 1.971 kJ/mol 10.365 kJ/mol -8.395
* Info * 80.0 deg 1.324 kJ/mol 12.461 kJ/mol -11.136
* Info * 90.0 deg 0.819 kJ/mol 13.214 kJ/mol -12.395
* Info * 100.0 deg 0.428 kJ/mol 12.264 kJ/mol -11.835
* Info * 110.0 deg 0.148 kJ/mol 9.969 kJ/mol -9.821
* Info * 120.0 deg 0.000 kJ/mol 7.136 kJ/mol -7.136
* Info * 130.0 deg 0.022 kJ/mol 4.450 kJ/mol -4.428
* Info * 140.0 deg 0.256 kJ/mol 2.316 kJ/mol -2.060
* Info * 150.0 deg 0.705 kJ/mol 0.943 kJ/mol -0.238
* Info * 160.0 deg 1.291 kJ/mol 0.294 kJ/mol 0.997
* Info * 170.0 deg 1.833 kJ/mol 0.087 kJ/mol 1.746
* Info * 180.0 deg 2.114 kJ/mol 0.000 kJ/mol 2.114
* Info * Dihedral barriers [1.85973955 1.85973955 0. 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: [2.00263454 2.00263454 1.90578873 1.23978177 2.00263454 2.00263454]
* Info * Validating the fitted force field...
* 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')