Molecular mechanics#
Force field generation#
import veloxchem as vlx
import numpy as np
# load B3LYP optimized geometry in the xyz format
molecule = vlx.Molecule.read_xyz_file("../input_files/hs276_optim_b3lyp_dev2-svp.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.
This part of the code load the MM Forfield Generator and create an initial topology for the molecule object. By default, RESP charges are computed at the B3LYP/6-31G as recommended.
ff_gen = vlx.MMForceFieldGenerator()
ff_gen.create_topology(molecule)
ff_gen.write_gromacs_files('HS-276_initial', 'MOL')
* Info * Using 6-31G* basis set for RESP charges...
* Info * Sum of partial charges is not a whole number.
* Info * Compensating by removing 1.000e-06 from the largest charge.
* 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.067 deg
* Info * Updated bond angle 1-14-15 (ca-na-cc) to 126.448 deg
* Info * Updated bond angle 2-1-14 (ca-ca-na) to 131.608 deg
* Info * Updated bond angle 7-9-10 (ca-ca-cc) to 135.484 deg
* Info * Updated bond angle 9-1-14 (ca-ca-na) to 107.725 deg
* Info * Updated bond angle 14-15-17 (na-cc-cd) to 128.299 deg
* Info * Updated bond angle 19-21-22 (cd-cc-cc) to 129.072 deg
* Info * Updated bond angle 21-22-23 (cc-cc-cd) to 128.378 deg
The Force-Field Generator can identify rotatable bonds.
rot_bonds = ff_gen.rotatable_bonds
print(rot_bonds)
[[14, 15], [21, 22], [28, 29], [29, 31], [31, 32]]
Force-Field Definition and Reparametrization#
We will here focus on the reparametrization of the barrier around the rotatable bond [21, 22]. We use the file 16-21-22-27.xyz from the data folder containing the results from a relaxed scan around rotatable bond [21, 22], i. e. optimized geometries and energies (in the comment line for each geometries).
ff_gen.reparameterize_dihedrals(rotatable_bond=(21,22),
scan_file="../input_files/16-21-22-27.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 * ../input_files/16-21-22-27.xyz
* Info * Performing dihedral scan for MM baseline by excluding the involved dihedral barriers...
* Info * Target dihedral angle: 16-21-22-27
* Info * 180.0 deg...
* Info * 190.0 deg...
* Info * 200.0 deg...
* Info * 210.0 deg...
* Info * 220.0 deg...
* Info * 230.0 deg...
* Info * 240.0 deg...
* Info * 250.0 deg...
* Info * 260.0 deg...
* Info * 270.0 deg...
* Info * 280.0 deg...
* Info * 290.0 deg...
* Info * 300.0 deg...
* Info * 310.0 deg...
* Info * 320.0 deg...
* Info * 330.0 deg...
* Info * 340.0 deg...
* Info * 350.0 deg...
* Info * 360.0 deg...
* 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 * 180.0 deg 0.000 kJ/mol 0.056 kJ/mol -0.056
* Info * 190.0 deg 0.872 kJ/mol 0.000 kJ/mol 0.872
* Info * 200.0 deg 3.294 kJ/mol 0.097 kJ/mol 3.197
* Info * 210.0 deg 7.151 kJ/mol 0.799 kJ/mol 6.352
* Info * 220.0 deg 12.105 kJ/mol 2.385 kJ/mol 9.720
* Info * 230.0 deg 17.593 kJ/mol 4.878 kJ/mol 12.715
* Info * 240.0 deg 22.920 kJ/mol 8.047 kJ/mol 14.873
* Info * 250.0 deg 27.394 kJ/mol 11.478 kJ/mol 15.916
* Info * 260.0 deg 30.445 kJ/mol 14.445 kJ/mol 16.000
* Info * 270.0 deg 31.708 kJ/mol 15.853 kJ/mol 15.855
* Info * 280.0 deg 31.070 kJ/mol 15.105 kJ/mol 15.965
* Info * 290.0 deg 28.682 kJ/mol 12.570 kJ/mol 16.112
* Info * 300.0 deg 24.933 kJ/mol 9.354 kJ/mol 15.580
* Info * 310.0 deg 20.393 kJ/mol 6.386 kJ/mol 14.007
* Info * 320.0 deg 15.714 kJ/mol 4.246 kJ/mol 11.469
* Info * 330.0 deg 11.526 kJ/mol 3.135 kJ/mol 8.390
* Info * 340.0 deg 8.324 kJ/mol 2.897 kJ/mol 5.426
* Info * 350.0 deg 6.379 kJ/mol 3.082 kJ/mol 3.297
* Info * 360.0 deg 5.710 kJ/mol 3.165 kJ/mol 2.546
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 16.112 kJ/mol
* Info * Standard deviation: 5.743 kJ/mol
* Info * Fitting the dihedral parameters...
* Info * Optimizing dihedral via least squares fitting...
* Info * New fitted barriers: [2.22525967 2.22525967 2.22525967 2.22525967]
* Info * Validating the fitted force field...
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 180.0 deg 0.000 kJ/mol 0.056 kJ/mol -0.056
* Info * 190.0 deg 0.400 kJ/mol 0.000 kJ/mol 0.400
* Info * 200.0 deg 1.461 kJ/mol 0.097 kJ/mol 1.364
* Info * 210.0 deg 3.233 kJ/mol 0.799 kJ/mol 2.434
* Info * 220.0 deg 5.631 kJ/mol 2.385 kJ/mol 3.246
* Info * 230.0 deg 8.398 kJ/mol 4.878 kJ/mol 3.520
* Info * 240.0 deg 11.168 kJ/mol 8.047 kJ/mol 3.121
* Info * 250.0 deg 13.557 kJ/mol 11.478 kJ/mol 2.079
* Info * 260.0 deg 15.247 kJ/mol 14.445 kJ/mol 0.802
* Info * 270.0 deg 16.038 kJ/mol 15.853 kJ/mol 0.185
* Info * 280.0 deg 15.873 kJ/mol 15.105 kJ/mol 0.768
* Info * 290.0 deg 14.845 kJ/mol 12.570 kJ/mol 2.275
* Info * 300.0 deg 13.181 kJ/mol 9.354 kJ/mol 3.827
* Info * 310.0 deg 11.198 kJ/mol 6.386 kJ/mol 4.811
* Info * 320.0 deg 9.240 kJ/mol 4.246 kJ/mol 4.994
* Info * 330.0 deg 7.608 kJ/mol 3.135 kJ/mol 4.473
* Info * 340.0 deg 6.491 kJ/mol 2.897 kJ/mol 3.593
* Info * 350.0 deg 5.906 kJ/mol 3.082 kJ/mol 2.824
* Info * 360.0 deg 5.710 kJ/mol 3.165 kJ/mol 2.546
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 4.994 kJ/mol
* Info * Standard deviation: 1.525 kJ/mol
* Info * Dihedral MM parameters have been reparameterized and updated in the topology.
The rotational barrier is now much improved but the relative energies between the two conformer minima are still not well reproduced, leading to significant errors in the statistical dihedral distributions. As a remedy to this situation, an additional dihedral potential can be added.
ff_gen.add_dihedral((16, 21, 22, 27),
barrier=0.0,
phase=180.0,
periodicity=1)
* Info * Added dihedral 16-21-22-27
We can now restart the reparametrization procedure with this added dihedral angle.
ff_gen.reparameterize_dihedrals(rotatable_bond=(21,22),
scan_file="../input_files/16-21-22-27.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 * ../input_files/16-21-22-27.xyz
* Info * Performing dihedral scan for MM baseline by excluding the involved dihedral barriers...
* Info * Target dihedral angle: 16-21-22-27
* Info * 180.0 deg...
* Info * 190.0 deg...
* Info * 200.0 deg...
* Info * 210.0 deg...
* Info * 220.0 deg...
* Info * 230.0 deg...
* Info * 240.0 deg...
* Info * 250.0 deg...
* Info * 260.0 deg...
* Info * 270.0 deg...
* Info * 280.0 deg...
* Info * 290.0 deg...
* Info * 300.0 deg...
* Info * 310.0 deg...
* Info * 320.0 deg...
* Info * 330.0 deg...
* Info * 340.0 deg...
* Info * 350.0 deg...
* Info * 360.0 deg...
* Info * Dihedral barriers [2.22525967 2.22525967 0. 2.22525967 2.22525967] will be used as initial guess.
* Info * Fitting the dihedral parameters...
* Info * Optimizing dihedral via least squares fitting...
* Info * New fitted barriers: [2.30904737 2.30904737 1.15807886 2.30904737 2.30904737]
* Info * Validating the fitted force field...
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 180.0 deg 0.000 kJ/mol 0.056 kJ/mol -0.056
* Info * 190.0 deg 0.402 kJ/mol 0.000 kJ/mol 0.402
* Info * 200.0 deg 1.470 kJ/mol 0.097 kJ/mol 1.373
* Info * 210.0 deg 3.246 kJ/mol 0.799 kJ/mol 2.447
* Info * 220.0 deg 5.637 kJ/mol 2.385 kJ/mol 3.252
* Info * 230.0 deg 8.378 kJ/mol 4.878 kJ/mol 3.499
* Info * 240.0 deg 11.092 kJ/mol 8.047 kJ/mol 3.044
* Info * 250.0 deg 13.387 kJ/mol 11.478 kJ/mol 1.909
* Info * 260.0 deg 14.940 kJ/mol 14.445 kJ/mol 0.495
* Info * 270.0 deg 15.550 kJ/mol 15.853 kJ/mol -0.302
* Info * 280.0 deg 15.164 kJ/mol 15.105 kJ/mol 0.059
* Info * 290.0 deg 13.883 kJ/mol 12.570 kJ/mol 1.313
* Info * 300.0 deg 11.946 kJ/mol 9.354 kJ/mol 2.593
* Info * 310.0 deg 9.689 kJ/mol 6.386 kJ/mol 3.302
* Info * 320.0 deg 7.472 kJ/mol 4.246 kJ/mol 3.226
* Info * 330.0 deg 5.615 kJ/mol 3.135 kJ/mol 2.479
* Info * 340.0 deg 4.323 kJ/mol 2.897 kJ/mol 1.425
* Info * 350.0 deg 3.628 kJ/mol 3.082 kJ/mol 0.546
* Info * 360.0 deg 3.394 kJ/mol 3.165 kJ/mol 0.229
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 3.499 kJ/mol
* Info * Standard deviation: 1.278 kJ/mol
* Info * Dihedral MM parameters have been reparameterized and updated in the topology.
Another dihedral angle can be added for a better agreement with the QM Potential Energy Surface.
ff_gen.add_dihedral((16, 21, 22, 27),
barrier=0.0,
phase=0.0,
periodicity=4)
* Info * Added dihedral 16-21-22-27
And start a new parametrization with this added dihedral angle.
ff_gen.reparameterize_dihedrals(rotatable_bond=(21,22),
scan_file="../input_files/16-21-22-27.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 * ../input_files/16-21-22-27.xyz
* Info * Performing dihedral scan for MM baseline by excluding the involved dihedral barriers...
* Info * Target dihedral angle: 16-21-22-27
* Info * 180.0 deg...
* Info * 190.0 deg...
* Info * 200.0 deg...
* Info * 210.0 deg...
* Info * 220.0 deg...
* Info * 230.0 deg...
* Info * 240.0 deg...
* Info * 250.0 deg...
* Info * 260.0 deg...
* Info * 270.0 deg...
* Info * 280.0 deg...
* Info * 290.0 deg...
* Info * 300.0 deg...
* Info * 310.0 deg...
* Info * 320.0 deg...
* Info * 330.0 deg...
* Info * 340.0 deg...
* Info * 350.0 deg...
* Info * 360.0 deg...
* Info * Dihedral barriers [2.30904737 2.30904737 1.15807886 0. 2.30904737 2.30904737] will be used as initial guess.
* Info * Fitting the dihedral parameters...
* Info * Optimizing dihedral via least squares fitting...
* Info * New fitted barriers: [2.30954291 2.30954291 1.10951333 1.79256159 2.30954291 2.30954291]
* Info * Validating the fitted force field...
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 180.0 deg 0.016 kJ/mol 0.056 kJ/mol -0.040
* Info * 190.0 deg 0.000 kJ/mol 0.000 kJ/mol 0.000
* Info * 200.0 deg 0.008 kJ/mol 0.097 kJ/mol -0.089
* Info * 210.0 deg 0.580 kJ/mol 0.799 kJ/mol -0.219
* Info * 220.0 deg 2.189 kJ/mol 2.385 kJ/mol -0.196
* Info * 230.0 deg 4.936 kJ/mol 4.878 kJ/mol 0.058
* Info * 240.0 deg 8.446 kJ/mol 8.047 kJ/mol 0.399
* Info * 250.0 deg 11.957 kJ/mol 11.478 kJ/mol 0.479
* Info * 260.0 deg 14.581 kJ/mol 14.445 kJ/mol 0.136
* Info * 270.0 deg 15.619 kJ/mol 15.853 kJ/mol -0.234
* Info * 280.0 deg 14.821 kJ/mol 15.105 kJ/mol -0.284
* Info * 290.0 deg 12.486 kJ/mol 12.570 kJ/mol -0.084
* Info * 300.0 deg 9.349 kJ/mol 9.354 kJ/mol -0.004
* Info * 310.0 deg 6.310 kJ/mol 6.386 kJ/mol -0.077
* Info * 320.0 deg 4.098 kJ/mol 4.246 kJ/mol -0.147
* Info * 330.0 deg 3.033 kJ/mol 3.135 kJ/mol -0.102
* Info * 340.0 deg 2.952 kJ/mol 2.897 kJ/mol 0.055
* Info * 350.0 deg 3.321 kJ/mol 3.082 kJ/mol 0.239
* Info * 360.0 deg 3.507 kJ/mol 3.165 kJ/mol 0.343
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 0.479 kJ/mol
* Info * Standard deviation: 0.213 kJ/mol
* Info * Dihedral MM parameters have been reparameterized and updated in the topology.
Rotatable Bond |
Max difference (kJ/mol) |
Std deviation (kJ/mol) |
|---|---|---|
[21, 22] |
0.477 |
0.212 |
With this excellent agreement between the QM and MM Potential Energy Surface, we can now save the topology as, i. e., Gromacs ouptut files.
ff_gen.write_gromacs_files('HS-276_final', 'MOL')