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.874 kJ/mol 0.000 kJ/mol 0.874
* Info * 200.0 deg 3.298 kJ/mol 0.097 kJ/mol 3.201
* Info * 210.0 deg 7.154 kJ/mol 0.799 kJ/mol 6.355
* Info * 220.0 deg 12.108 kJ/mol 2.385 kJ/mol 9.723
* Info * 230.0 deg 17.595 kJ/mol 4.878 kJ/mol 12.717
* Info * 240.0 deg 22.921 kJ/mol 8.047 kJ/mol 14.873
* Info * 250.0 deg 27.395 kJ/mol 11.478 kJ/mol 15.917
* Info * 260.0 deg 30.446 kJ/mol 14.445 kJ/mol 16.001
* Info * 270.0 deg 31.709 kJ/mol 15.853 kJ/mol 15.856
* Info * 280.0 deg 31.072 kJ/mol 15.105 kJ/mol 15.967
* Info * 290.0 deg 28.685 kJ/mol 12.570 kJ/mol 16.115
* Info * 300.0 deg 24.938 kJ/mol 9.354 kJ/mol 15.584
* Info * 310.0 deg 20.400 kJ/mol 6.386 kJ/mol 14.013
* Info * 320.0 deg 15.722 kJ/mol 4.246 kJ/mol 11.476
* Info * 330.0 deg 11.535 kJ/mol 3.135 kJ/mol 8.399
* Info * 340.0 deg 8.332 kJ/mol 2.897 kJ/mol 5.435
* Info * 350.0 deg 6.382 kJ/mol 3.082 kJ/mol 3.300
* Info * 360.0 deg 5.711 kJ/mol 3.165 kJ/mol 2.546
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 16.115 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.22536935 2.22536935 2.22536935 2.22536935]
* 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.401 kJ/mol 0.000 kJ/mol 0.401
* Info * 200.0 deg 1.465 kJ/mol 0.097 kJ/mol 1.368
* Info * 210.0 deg 3.237 kJ/mol 0.799 kJ/mol 2.438
* Info * 220.0 deg 5.634 kJ/mol 2.385 kJ/mol 3.249
* Info * 230.0 deg 8.400 kJ/mol 4.878 kJ/mol 3.522
* Info * 240.0 deg 11.169 kJ/mol 8.047 kJ/mol 3.122
* Info * 250.0 deg 13.559 kJ/mol 11.478 kJ/mol 2.081
* Info * 260.0 deg 15.249 kJ/mol 14.445 kJ/mol 0.804
* Info * 270.0 deg 16.040 kJ/mol 15.853 kJ/mol 0.187
* Info * 280.0 deg 15.875 kJ/mol 15.105 kJ/mol 0.770
* Info * 290.0 deg 14.849 kJ/mol 12.570 kJ/mol 2.279
* Info * 300.0 deg 13.186 kJ/mol 9.354 kJ/mol 3.833
* Info * 310.0 deg 11.205 kJ/mol 6.386 kJ/mol 4.818
* Info * 320.0 deg 9.248 kJ/mol 4.246 kJ/mol 5.002
* Info * 330.0 deg 7.617 kJ/mol 3.135 kJ/mol 4.482
* Info * 340.0 deg 6.499 kJ/mol 2.897 kJ/mol 3.602
* Info * 350.0 deg 5.910 kJ/mol 3.082 kJ/mol 2.828
* Info * 360.0 deg 5.711 kJ/mol 3.165 kJ/mol 2.546
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 5.002 kJ/mol
* Info * Standard deviation: 1.527 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.22536935 2.22536935 0. 2.22536935 2.22536935] will be used as initial guess.
* Info * Fitting the dihedral parameters...
* Info * Optimizing dihedral via least squares fitting...
* Info * New fitted barriers: [2.30913649 2.30913649 1.16030449 2.30913649 2.30913649]
* 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.404 kJ/mol 0.000 kJ/mol 0.404
* Info * 200.0 deg 1.473 kJ/mol 0.097 kJ/mol 1.376
* Info * 210.0 deg 3.249 kJ/mol 0.799 kJ/mol 2.450
* Info * 220.0 deg 5.639 kJ/mol 2.385 kJ/mol 3.254
* Info * 230.0 deg 8.379 kJ/mol 4.878 kJ/mol 3.500
* Info * 240.0 deg 11.091 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.303
* Info * 280.0 deg 15.164 kJ/mol 15.105 kJ/mol 0.058
* Info * 290.0 deg 13.883 kJ/mol 12.570 kJ/mol 1.313
* Info * 300.0 deg 11.948 kJ/mol 9.354 kJ/mol 2.595
* Info * 310.0 deg 9.692 kJ/mol 6.386 kJ/mol 3.305
* Info * 320.0 deg 7.476 kJ/mol 4.246 kJ/mol 3.230
* Info * 330.0 deg 5.620 kJ/mol 3.135 kJ/mol 2.484
* Info * 340.0 deg 4.327 kJ/mol 2.897 kJ/mol 1.430
* Info * 350.0 deg 3.627 kJ/mol 3.082 kJ/mol 0.545
* Info * 360.0 deg 3.391 kJ/mol 3.165 kJ/mol 0.226
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 3.500 kJ/mol
* Info * Standard deviation: 1.279 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.30913649 2.30913649 1.16030449 0. 2.30913649 2.30913649] will be used as initial guess.
* Info * Fitting the dihedral parameters...
* Info * Optimizing dihedral via least squares fitting...
* Info * New fitted barriers: [2.30964534 2.30964534 1.11148563 1.79468455 2.30964534 2.30964534]
* Info * Validating the fitted force field...
* Info * Dihedral MM energy(rel) QM energy(rel) diff
* Info * ---------------------------------------------------------------
* Info * 180.0 deg 0.015 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.009 kJ/mol 0.097 kJ/mol -0.088
* Info * 210.0 deg 0.580 kJ/mol 0.799 kJ/mol -0.219
* Info * 220.0 deg 2.187 kJ/mol 2.385 kJ/mol -0.199
* Info * 230.0 deg 4.933 kJ/mol 4.878 kJ/mol 0.055
* Info * 240.0 deg 8.442 kJ/mol 8.047 kJ/mol 0.395
* Info * 250.0 deg 11.955 kJ/mol 11.478 kJ/mol 0.477
* Info * 260.0 deg 14.580 kJ/mol 14.445 kJ/mol 0.135
* Info * 270.0 deg 15.618 kJ/mol 15.853 kJ/mol -0.235
* Info * 280.0 deg 14.820 kJ/mol 15.105 kJ/mol -0.285
* Info * 290.0 deg 12.485 kJ/mol 12.570 kJ/mol -0.085
* Info * 300.0 deg 9.348 kJ/mol 9.354 kJ/mol -0.006
* Info * 310.0 deg 6.309 kJ/mol 6.386 kJ/mol -0.078
* Info * 320.0 deg 4.098 kJ/mol 4.246 kJ/mol -0.148
* Info * 330.0 deg 3.035 kJ/mol 3.135 kJ/mol -0.100
* Info * 340.0 deg 2.955 kJ/mol 2.897 kJ/mol 0.057
* Info * 350.0 deg 3.320 kJ/mol 3.082 kJ/mol 0.237
* Info * 360.0 deg 3.504 kJ/mol 3.165 kJ/mol 0.339
* Info * Summary of validation
* Info * ---------------------
* Info * Maximum difference: 0.477 kJ/mol
* Info * Standard deviation: 0.212 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')