Finding transition states is of crucial importance for computational chemistry, and there are many methods for doing so. However, most of these methods require some initial guess of the transition state structure. If the guess is not good enough, the method might not converge. Reaction Force Field interpolation allows for generating such guesses at a very cheap computational level by relying mostly on molecular mechanics.
Force field interpolation works by describing the reactant and product states with classical force fields and respectively. By taking a linear combination
one can crudely describe the Potential Energy Surface (PES) along the reaction. By doing an MM optimisation of the structure at various values of , the relative MM energy can be obtained. The optimised structure can also be used to perform a single-point quantum mechanical (QM) calculation, which will give an more accurate energy. The structure at which these energies are highest is a good guess for the transition state structure. This workflow can be executed with just a couple lines of code. Let’s consider the Bromine substitution of Ethyl chloride:
The following code will perform the entire workflow and visualise the results:
import veloxchem as vlximport veloxchem as vlx
rea1 = vlx.Molecule.read_smiles("[Br]")
rea2 = vlx.Molecule.read_smiles("CCCl")
rea1.set_charge(-1)
pro1 = vlx.Molecule.read_smiles("[Cl]")
pro2 = vlx.Molecule.read_smiles("CCBr")
pro1.set_charge(-1)
tsguesser = vlx.TransitionStateGuesser()
results = tsguesser.find_transition_state([rea1, rea2], [pro1, pro2])
tsguesser.show_results(results)Output
* Info * Building forcefields. Disable mute_ff_build to see detailed output.
Reaction summary
1 breaking bonds:
ReaType ProType ID - ReaType ProType ID
c3 c3 3 - cl cl 4
1 forming bonds:
ReaType ProType ID - ReaType ProType ID
br br 1 - c3 c3 3
* Info * System has charge -1.0 and multiplicity 1. Provide correct values if this is wrong.
* Info * Saving reactant and product forcefield as json to ts_data_1775638362
* Info * Building MM systems for the transition state guess. Disable mute_ff_build to see detailed output.
* Info * Saving systems as xml to ts_data_1775638362/systems
* Info * Saving systems to /home/david/KTH_Postdoc/echem/docs/mol_struct/ts_data_1775638362/systems
Starting MM scan
MD steps: 1000
conf. snapshots: 1
MD temperature: 600 K
MD step size: 0.001 ps
folder name: ts_data_1775638362
saving MD traj: False
Lambda | E1 (kj/mol) | E2 (kj/mol) | V (kj/mol) | n_conf
------------------------------------------------------------
0.00 -69.96 556.68 -69.96 1
0.05 -67.94 510.30 -39.03 1
0.10 -63.26 464.65 -10.47 1
0.15 -55.43 419.62 15.83 1
0.20 -44.34 376.34 39.80 1
0.25 -28.05 332.86 62.18 1
0.30 -9.48 291.35 80.77 1
0.35 13.48 250.38 96.39 1
0.40 39.31 211.46 108.17 1
0.45 72.57 171.08 116.90 1
0.50 126.39 118.60 122.50 1
0.55 169.10 81.20 120.75 1
0.60 217.97 42.76 112.84 1
0.65 268.02 10.61 100.70 1
0.70 319.27 -16.66 84.12 1
0.75 371.54 -39.20 63.49 1
0.80 423.68 -54.94 40.79 1
0.85 478.32 -70.39 11.92 1
0.90 532.96 -79.50 -18.25 1
0.95 586.94 -84.61 -51.03 1
1.00 640.11 -86.38 -86.38 1
* Info * Found highest MM E: 122.499 at Lammba: 0.5 and conformer index: 0.
* Info * Saving results to ts_results_1775638362.h5
To better understand the workflow and it’s limitations, let’s go through it step by step.
Step by step calculation of reaction¶
To create interpolated forcefields, we need to provide reactant and product structures as molecule objects. The program will then calculate partial charges, but it is also possible to provide those manually. The program figures out which bonds are broken and formed. It then does an MM optimisation of the combined reactand and product systems which will be used as starting points for the scan.
rea1 = vlx.Molecule.read_smiles("[Br]")
rea2 = vlx.Molecule.read_smiles("CCCl")
rea1.set_charge(-1)
pro1 = vlx.Molecule.read_smiles("[Cl]")
pro2 = vlx.Molecule.read_smiles("CCBr")
pro1.set_charge(-1)
tsguesser = vlx.TransitionStateGuesser()
tsguesser.build_forcefields([rea1, rea2], [pro1, pro2])* Info * Building forcefields. Disable mute_ff_build to see detailed output.
Reaction summary
1 breaking bonds:
ReaType ProType ID - ReaType ProType ID
c3 c3 3 - cl cl 4
1 forming bonds:
ReaType ProType ID - ReaType ProType ID
br br 1 - c3 c3 3
* Info * System has charge -1.0 and multiplicity 1. Provide correct values if this is wrong.
* Info * Saving reactant and product forcefield as json to ts_data_1775638402
{'breaking_bonds': {(2, 3)},
'forming_bonds': {(0, 2)},
'static_bonds': {(1, 2), (1, 4), (1, 5), (1, 6), (2, 7), (2, 8)},
'reactant': <veloxchem.mmforcefieldgenerator.MMForceFieldGenerator at 0x7fbb2f4aec90>,
'product': <veloxchem.mmforcefieldgenerator.MMForceFieldGenerator at 0x7fbb2f4358e0>}The optimised structures can be inspected as follows:
tsguesser.reactant.molecule.show()
tsguesser.product.molecule.show()If this looks reasonable, we can proceed to scan the reaction. By default, the reaction is split up in 21 steps. On every step, the structure is minimised to the new forcefield, then sampled for 1000 steps at 600K, and then optimised again. If the value of decreases (or increases), this indicates that the obtained structures are could be part of different valleys in the energy landscape. The program will then try to do a conformer search to find the lowest lying conformer.
tsguesser.build_systems()
mm_results = tsguesser.scan_mm()* Info * Building MM systems for the transition state guess. Disable mute_ff_build to see detailed output.
* Info * Saving systems as xml to ts_data_1775638402/systems
* Info * Saving systems to /home/david/KTH_Postdoc/echem/docs/mol_struct/ts_data_1775638402/systems
Starting MM scan
MD steps: 1000
conf. snapshots: 1
MD temperature: 600 K
MD step size: 0.001 ps
folder name: ts_data_1775638402
saving MD traj: False
Lambda | E1 (kj/mol) | E2 (kj/mol) | V (kj/mol) | n_conf
------------------------------------------------------------
0.00 -70.00 557.39 -70.00 1
0.05 -68.03 510.35 -39.11 1
0.10 -63.34 464.57 -10.55 1
0.15 -53.99 418.69 16.91 1
0.20 -42.78 374.92 40.76 1
0.25 -29.69 333.85 61.20 1
0.30 -11.23 292.63 79.93 1
0.35 11.72 251.57 95.67 1
0.40 39.58 211.09 108.19 1
0.45 74.68 169.32 117.27 1
0.50 124.24 121.24 122.74 1
0.55 170.20 79.06 120.07 1
0.60 217.95 42.68 112.79 1
0.65 267.82 10.63 100.65 1
0.70 319.07 -16.66 84.06 1
0.75 371.47 -39.24 63.44 1
0.80 424.68 -57.15 39.22 1
0.85 478.49 -70.47 11.87 1
0.90 532.90 -79.59 -18.34 1
0.95 586.91 -84.69 -51.11 1
1.00 638.24 -84.26 -84.26 1
* Info * Found highest MM E: 122.738 at Lammba: 0.5 and conformer index: 0.
* Info * Saving results to ts_results_1775638402.h5
The results can then be visualised.
tsguesser.show_results(mm_results)While the MM energies give a reasonable description of the energetic contributions from coloumbic and van der Waals interactions, they do not capture bond breaking and formation. To get a better estimate of the energy, we can perform single-point QM calculations on the MM optimised structures. Here, we use DFT with the B3LYP functional and the ‘def2-SVP’ basis set which are the default settings.
tsguesser.scf_xcfun = 'b3lyp'
tsguesser.scf_basis = 'def2-svp'
scf_results = tsguesser.scan_qm()* Info * Disable mute_scf to see detailed output.
Starting QM scan
QM parameters:
Basis: def2-svp
DFT xc fun: PBE0
Max conf.: 5
Lambda | Conf. i | Rel. E (kJ/mol) | MM V (kJ/mol)
------------------------------------------------------------
0.00 1 0.000 -69.997
0.05 1 2.242 -39.107
0.10 1 8.392 -10.553
0.15 1 18.074 16.910
0.20 1 30.481 40.761
0.25 1 44.823 61.197
0.30 1 60.546 79.926
0.35 1 76.708 95.667
0.40 1 92.064 108.187
0.45 1 105.133 117.269
0.50 1 114.952 122.738
0.55 1 116.638 120.073
0.60 1 112.774 112.787
0.65 1 103.911 100.649
0.70 1 91.821 84.062
0.75 1 78.780 63.440
0.80 1 66.107 39.219
0.85 1 56.152 11.871
0.90 1 49.591 -18.338
0.95 1 47.493 -51.105
1.00 1 50.402 -84.261
* Info * Found highest QM E: -8171996.406 at Lambda: 0.55 and conformer index: 0.
* Info * Saving results to ts_results_1775638402.h5
tsguesser.show_results(scf_results)And these results can be visualised as well. Note how in this (simple) example, the peak of the MM energies coincides with the SCF peak. In more complex reactions, this is not always the case, but it is usually close and the MM peak might already be a good enough guess.
Then we can proceed with a proper QM transition state optimisation.
basis = vlx.MolecularBasis.read(tsguesser.molecule, "def2-svp")
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()
scf_drv.xcfun = "b3lyp"
scf_results = scf_drv.compute(tsguesser.molecule, basis)
scf_opt_drv = vlx.OptimizationDriver(scf_drv)
scf_opt_drv.ostream.mute()
scf_opt_drv.transition = True #We set the option to find transition states
opt_results = scf_opt_drv.compute(tsguesser.molecule, basis, scf_results)
vlx.Molecule.read_xyz_string(opt_results['final_geometry']).show()Conformer searching¶
Some systems might have many conformers along the PES, and in order to find the correct transition state, it is important to find the lowest lying conformer. When peak_conformer_search is set to True, a conformer search will be performed at the peak of the PES over a range of values. It is also possible to automatically perform a conformer search when a discontinuity in the PES is detected by setting discont_conformer_search to True. Lastly, setting force_conformer_search to True will perform a conformer search at every value. All of these functions are turned off by default because they are more time consuming. Depending on the system and the use case, the user has to decide which options are appropriate. In the following example, the longer carbon backbone gives the system more conformational freedom. Rescanning the entire PES is not necessary, but performing a conformer search at the peak is useful to obtain the correct conformer for the transition state guess.
rea1 = vlx.Molecule.read_smiles("[Br]")
rea2 = vlx.Molecule.read_smiles("CC(C)CC(Cl)C")
rea1.set_charge(-1)
pro1 = vlx.Molecule.read_smiles("[Cl]")
pro2 = vlx.Molecule.read_smiles("CC(C)CC(Br)C")
pro1.set_charge(-1)
tsguesser = vlx.TransitionStateGuesser()
tsguesser.build_forcefields([rea1, rea2], [pro1, pro2])
tsguesser.build_systems()
tsguesser.peak_conformer_search = True
tsguesser.peak_conformer_search_range = 2
results = tsguesser.scan_mm()
tsguesser.show_results(results)* Info * Building forcefields. Disable mute_ff_build to see detailed output.
Reaction summary
1 breaking bonds:
ReaType ProType ID - ReaType ProType ID
c3 c3 6 - cl cl 7
1 forming bonds:
ReaType ProType ID - ReaType ProType ID
br br 1 - c3 c3 6
* Info * System has charge -1.0 and multiplicity 1. Provide correct values if this is wrong.
* Info * Saving reactant and product forcefield as json to ts_data_1775639120
* Info * Building MM systems for the transition state guess. Disable mute_ff_build to see detailed output.
* Info * Saving systems as xml to ts_data_1775639120/systems
* Info * Saving systems to /home/david/KTH_Postdoc/echem/docs/mol_struct/ts_data_1775639120/systems
Starting MM scan
MD steps: 1000
conf. snapshots: 1
MD temperature: 600 K
MD step size: 0.001 ps
folder name: ts_data_1775639120
saving MD traj: False
Lambda | E1 (kj/mol) | E2 (kj/mol) | V (kj/mol) | n_conf
------------------------------------------------------------
0.00 -96.91 644.98 -96.91 1
0.05 -92.72 597.28 -58.22 1
0.10 -83.41 547.77 -20.29 1
0.15 -69.29 500.58 16.19 1
0.20 -51.67 452.36 49.14 1
0.25 -29.58 404.09 78.84 1
0.30 -3.27 356.77 104.74 1
0.35 26.83 307.93 125.22 1
0.40 69.88 247.99 141.12 1
0.45 112.12 199.02 151.23 1
0.50 159.24 152.83 156.03 1
0.55 210.57 104.10 152.01 1
0.60 271.94 71.98 151.96 1
0.65 322.52 17.76 124.43 1
0.70 379.67 -15.98 102.72 1
0.75 436.60 -47.23 73.73 1
0.80 495.29 -78.08 36.60 1
0.85 552.96 -97.91 -0.28 1
0.90 611.51 -114.33 -41.75 1
0.95 667.30 -124.90 -85.29 1
1.00 719.40 -127.86 -127.86 1
* Info * Found peak MM E: 156.032 at Lambda: 0.5.
* Info * Doing conformer search from Lambda: 0.4 to Lambda: 0.6.
Starting MM scan with conformer search for lambda values [0.4, 0.45, 0.5, 0.55, 0.6]
MD steps: 1000
conf. snapshots: 10
MD temperature: 600 K
MD step size: 0.001 ps
folder name: ts_data_1775639120
saving MD traj: False
Lambda | E1 (kj/mol) | E2 (kj/mol) | V (kj/mol) | n_conf
------------------------------------------------------------
0.40 67.68 249.46 140.39 5
0.45 110.06 200.43 150.73 3
0.50 158.03 151.22 154.62 2
0.55 211.09 103.81 152.09 3
0.60 267.60 55.76 140.50 2
* Info * Found highest MM E: 154.623 at Lammba: 0.5 and conformer index: 1.
* Info * Saving results to ts_results_1775639120.h5
Forced bond breaking¶
To correctly interpolate between the reactant and the product, each atom in the reactant and product need to match. This mapping is detected automatically, but can be provided manually. From this mapping the breaking and forming bonds are derived. It is even possible to specify bonds that should be forced to break which are then not allowed to reconnect. This is useful for example in certain catalysis reactions, where both the reactant and the product contain the same catalyst molecule. The optimal mapping would not break or form any bonds with these molecules, but the reaction of interest should involve this catalyst. The following example demonstrates this with methanol catalysing the proton transfer of a larger complex. In order to run the example faster, the partial charges are pre-computed and provided manually.
evb = vlx.EvbDriver()
q_rea = [[
-0.63376,
0.155898,
0.40852,
0.023114,
0.023114,
0.023114,
],
[
-0.053121, -0.224161, -0.146708, -0.198581, -0.104851, 0.057048,
0.211108, -0.00552, 0.582421, -0.578475, -0.239007, -0.200459,
0.076475, 0.046494, -0.014331, -0.064532, -0.241586, -0.879622,
0.102603, 0.135941, 0.12827, 0.137736, 0.085068, 0.010794,
0.056581, 0.118934, 0.118934, 0.118934, 0.071642, 0.071642,
0.084769, 0.084769, 0.084769, 0.068894, 0.068894, 0.068894,
0.11978, 0.11978, 0.11978
]]
q_pro = [
[
-0.637505,
0.152377,
0.411487,
0.024547,
0.024547,
0.024547,
],
[
-0.088958, -0.19184, -0.147759, -0.194697, -0.060166, 0.060322,
0.057387, -0.279844, 0.484962, -0.274851, -0.152366,
-0.7389260000000002, 0.034466, -0.007965, -0.126105, -0.288341,
-0.009113, -0.684542, 0.112558, 0.144675, 0.133772, 0.134584, 0.1315,
0.137143, 0.100035, 0.100035, 0.100035, 0.071572, 0.071572, 0.100344,
0.100344, 0.100344, 0.139586, 0.139586, 0.139586, 0.066736, 0.066736,
0.066736, 0.450857
],
]
rea1 = vlx.Molecule.read_smiles("OC")
rea2 = vlx.Molecule.read_smiles("C1=CC=CC=C1C(C(C(=O)OC)C[N+](C)(C)C)([O-])")
pro1 = vlx.Molecule.read_smiles("OC")
pro2 = vlx.Molecule.read_smiles("C1=CC=CC=C1C(C(=C(OC)[O-])C[N+](C)(C)C)O")
rea_charges = q_rea
pro_charges = q_pro
breaking_bonds = {(0, 2)}
tsguesser = vlx.TransitionStateGuesser()
tsguesser.build_forcefields(
[rea1, rea2],
[pro1, pro2],
forced_breaking_bonds={(1, 3)},
reactant_partial_charges=rea_charges,
product_partial_charges=pro_charges,
)* Info * Building forcefields. Disable mute_ff_build to see detailed output.
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[15], line 50
46 rea_charges = q_rea
47 pro_charges = q_pro
48 breaking_bonds = {(0, 2)}
49 tsguesser = vlx.TransitionStateGuesser()
---> 50 tsguesser.build_forcefields(
51 [rea1, rea2],
52 [pro1, pro2],
53 forced_breaking_bonds={(1, 3)},
File ~/.conda/envs/echem/lib/python3.12/site-packages/veloxchem/tsguesser.py:170, in TransitionStateGuesser.build_forcefields(self, reactant, product, **build_forcefields_kwargs)
167 self.ostream.flush()
168 self.ostream.mute()
--> 170 self.reactant, self.product, self.forming_bonds, self.breaking_bonds, reactants, products, product_mapping = self.ffbuilder.build_forcefields(
171 reactant=reactant,
172 product=product,
173 **build_forcefields_kwargs,
174 )
176 if self.mute_ff_build:
177 self.ostream.unmute()
File ~/.conda/envs/echem/lib/python3.12/site-packages/veloxchem/reaffbuilder.py:180, in ReactionForceFieldBuilder.build_forcefields(self, reactant, product, reactant_partial_charges, product_partial_charges, reactant_hessians, product_hessians, reactant_total_multiplicity, product_total_multiplicity, forced_breaking_bonds, forced_forming_bonds, product_mapping)
175 forced_breaking_bonds = {(bond[0] - 1, bond[1] - 1)
176 for bond in forced_breaking_bonds}
177 forced_forming_bonds = {(bond[0] - 1, bond[1] - 1)
178 for bond in forced_forming_bonds}
--> 180 product_ff, product_mapping = self._match_reactant_and_product(
181 reactant_ff,
182 reactant_ff.molecule.get_element_ids(),
183 product_ff,
184 product_ff.molecule.get_element_ids(),
185 forced_breaking_bonds,
186 forced_forming_bonds,
187 )
188 elif product_mapping is not None:
189 self.ostream.print_info(
190 f"Skipping reaction matching because the mapping {product_mapping} is already provided"
191 )
File ~/.conda/envs/echem/lib/python3.12/site-packages/veloxchem/reaffbuilder.py:544, in ReactionForceFieldBuilder._match_reactant_and_product(self, reactant_ff, rea_elems, product_ff, pro_elems, breaking_bonds, forming_bonds)
535 total_mapping, breaking_bonds, forming_bonds = rm.get_mapping(
536 reactant_ff,
537 rea_elems,
(...) 541 forming_bonds,
542 ) # type: ignore
543 if total_mapping is None:
--> 544 raise ValueError(
545 "Could not find a mapping between the reactant and product force fields."
546 )
547 total_mapping = {v: k for k, v in total_mapping.items()}
548 print_mapping = {k + 1: v + 1 for k, v in total_mapping.items()}
ValueError: Could not find a mapping between the reactant and product force fields.The output shows that the program found the correct forming and breaking bonds. We can then proceed with the mm scan.
tsguesser.build_systems()
mm_results = tsguesser.scan_mm()
tsguesser.show_results(mm_results)---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[16], line 1
----> 1 tsguesser.build_systems()
2 mm_results = tsguesser.scan_mm()
3 tsguesser.show_results(mm_results)
File ~/.conda/envs/echem/lib/python3.12/site-packages/veloxchem/tsguesser.py:224, in TransitionStateGuesser.build_systems(self, constraints)
218 self.ostream.print_info(
219 "Building MM systems for the transition state guess. Disable mute_ff_build to see detailed output."
220 )
221 self.ostream.flush()
223 self.systems, self.topology, _ = sysbuilder.build_systems(
--> 224 self.reactant,
225 self.product,
226 list(self.lambda_vector),
227 self.sys_builder_configuration,
228 constraints,
229 )
230 if self.mute_ff_build:
231 self.ostream.print_blank()
AttributeError: 'TransitionStateGuesser' object has no attribute 'reactant'Available parameters and options¶
There are a bunch more options available to customise the TS-guesser. Below is a complete list of all options with their default values.
To be added.