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...                          
                                                                                                                          
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[3], line 1
----> 1 ff_gen.reparameterize_dihedrals(rotatable_bond=(21,22),
      2                                 scan_file="../../data/md/27-22-21-16.xyz",
      3                                 visualize=True,
      4                                 verbose=False)

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/mmforcefieldgenerator.py:600, in MMForceFieldGenerator.reparameterize_dihedrals(self, rotatable_bond, scan_results, scan_file, visualize, fit_extrema, initial_validation, show_diff, verbose)
    597 self.ostream.print_blank()
    598 self.ostream.flush()
--> 600 initial_data = self.validate_force_field(0, verbose=verbose)
    602 qm_energies = np.array(initial_data['qm_scan_kJpermol'])
    603 mm_baseline = np.array(initial_data['mm_scan_kJpermol'])

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/mmforcefieldgenerator.py:3138, in MMForceFieldGenerator.validate_force_field(self, i, verbose)
   3136 geom = self.scan_geometries[i]
   3137 angles = self.scan_dih_angles[i]
-> 3138 mm_scan = self.perform_mm_scan(dih, geom, angles, verbose=verbose)
   3139 mm_scan = np.array(mm_scan) - min(mm_scan)
   3141 qm_scan = np.array(self.scan_energies[i]) - min(self.scan_energies[i])

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/mmforcefieldgenerator.py:3184, in MMForceFieldGenerator.perform_mm_scan(self, dihedral, geometries, angles, verbose)
   3176 for i, (geom, angle) in enumerate(zip(geometries, angles)):
   3177     # energy minimization with dihedral constraint
   3178     constraints = [
   3179         'set dihedral {} {} {} {} {}'.format(dihedral[0] + 1,
   3180                                              dihedral[1] + 1,
   3181                                              dihedral[2] + 1,
   3182                                              dihedral[3] + 1, angle)
   3183     ]
-> 3184     pot_energy = self.minimize_mm_energy(geom, constraints)
   3185     # convert to kJ/mol
   3186     pot_energy *= hartree_in_kjpermol()

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/mmforcefieldgenerator.py:3215, in MMForceFieldGenerator.minimize_mm_energy(self, molecule, constraints)
   3213 opt_drv = OptimizationDriver(grad_drv)
   3214 opt_drv.constraints = constraints
-> 3215 opt_results = opt_drv.compute(molecule)
   3216 final_mol = Molecule.read_xyz_string(opt_results['final_geometry'])
   3217 self.ostream.unmute()

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/optimizationdriver.py:380, in OptimizationDriver.compute(self, molecule, *args)
    377 with redirect_stdout(StringIO()) as fg_out, redirect_stderr(
    378         StringIO()) as fg_err:
    379     try:
--> 380         m = geometric.optimize.run_optimizer(
    381             customengine=opt_engine,
    382             coordsys=self.coordsys,
    383             check=self.check_interval,
    384             trust=default_trust,
    385             tmax=default_tmax,
    386             maxiter=self.max_iter,
    387             converge=self.conv_flags(),
    388             constraints=constr_filename,
    389             transition=self.transition,
    390             hessian=self.hessian,
    391             input=optinp_filename)
    392     except geometric.errors.HessianExit:
    393         hessian_exit = True

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/geometric/optimize.py:1318, in run_optimizer(**kwargs)
   1316     params.xyzout = prefix+"_optim.xyz"
   1317 if ic == 0:
-> 1318     progress = Optimize(coords, M, IC, engine, dirname, params)
   1319 else:
   1320     progress += Optimize(coords, M, IC, engine, dirname, params, print_info=False)

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/geometric/optimize.py:1116, in Optimize(coords, molecule, IC, engine, dirname, params, print_info)
   1088 """
   1089 Optimize the geometry of a molecule. This function used to contain the whole
   1090 optimization loop, which has since been moved to the Optimizer() class;
   (...)   1113     A molecule object for opt trajectory and energies
   1114 """
   1115 optimizer = Optimizer(coords, molecule, IC, engine, dirname, params, print_info)
-> 1116 return optimizer.optimizeGeometry()

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/geometric/optimize.py:1035, in Optimizer.optimizeGeometry(self)
   1033 if self.state == OPT_STATE.NEEDS_EVALUATION:
   1034     self.calcEnergyForce()
-> 1035     self.evaluateStep()
   1036 if self.recalcHess:
   1037     # If a Hessian recalculation is needed at this point, we need to
   1038     # call calcEnergyForce() which will compute the cartesian Hessian,
   1039     # convert it to IC, and then store it.
   1040     self.calcEnergyForce()

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/geometric/optimize.py:841, in Optimizer.evaluateStep(self)
    839 Quality = del_E/self.expect
    840 # The internal coordinate gradient (actually not really used in this function)
--> 841 self.G = self.IC.calcGrad(self.X, self.gradx).flatten()
    842 # For transition states, the quality factor decreases in both directions
    843 if params.transition and Quality > 1.0:

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/geometric/internal.py:1991, in InternalCoordinates.calcGrad(self, xyz, gradx)
   1989 def calcGrad(self, xyz, gradx):
   1990     q0 = self.calculate(xyz)
-> 1991     Ginv = self.GInverse(xyz)
   1992     Bmat = self.wilsonB(xyz)
   1993     # Internal coordinate gradient
   1994     # Gq = np.matrix(Ginv)*np.matrix(Bmat)*np.matrix(gradx).T

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/geometric/internal.py:3674, in DelocalizedInternalCoordinates.GInverse(self, xyz)
   3673 def GInverse(self, xyz):
-> 3674     return self.GInverse_SVD(xyz)

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/geometric/internal.py:1850, in InternalCoordinates.GInverse_SVD(self, xyz, sqrt, invMW)
   1848     G = self.GMatrix(xyz, invMW)
   1849     time_G = click()
-> 1850     U, S, VT = np.linalg.svd(G)
   1851     time_svd = click()
   1852 except np.linalg.LinAlgError:

File ~/miniconda3/envs/echem/lib/python3.12/site-packages/numpy/linalg/linalg.py:1681, in svd(a, full_matrices, compute_uv, hermitian)
   1678         gufunc = _umath_linalg.svd_n_s
   1680 signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
-> 1681 u, s, vh = gufunc(a, signature=signature, extobj=extobj)
   1682 u = u.astype(result_t, copy=False)
   1683 s = s.astype(_realType(result_t), copy=False)

KeyboardInterrupt: 

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.85974413 1.85974413 0.         1.85974413 1.85974413] will be used as initial guess.        
* Info * Fitting the dihedral parameters...                                                                               
* Info * Optimizing dihedral via least squares fitting...                                                                 
* Info * New fitted barriers: [1.98197174 1.98197174 1.97103413 1.98197174 1.98197174]                                    
* Info * Validating the fitted force field...                                                                             
                                                                                                                          
* Info * Summary of validation                                                                                            
* Info * ---------------------                                                                                            
* Info * Maximum difference: 2.368 kJ/mol                                                                                 
* Info * Standard deviation: 0.969 kJ/mol                                                                                 
                                                                                                                          
../../_images/d40f171a5f3da714d5ed49a97be8021746399060c2eaf509bc0e570f27742e28.png
* 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.98197174 1.98197174 1.97103413 0.         1.98197174 1.98197174] will be used as initial guess.
* Info * Fitting the dihedral parameters...                                                                               
* Info * Optimizing dihedral via least squares fitting...                                                                 
* Info * New fitted barriers: [2.00263461 2.00263461 1.90578785 1.23978392 2.00263461 2.00263461]                         
* 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                                                                                 
                                                                                                                          
../../_images/17b4bfcb882b635126a252911d633dcf34cd1c9d16da60153478a3532a3ed0f3.png
* 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')