Guessing transition states with force fireld interpolation#
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 \(E_{1}\) and \(E_{2}\) 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 \(\lambda\), 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 vlx
import 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_TS([rea1, rea2], [pro1, pro2])
tsguesser.show_results(results)
Show code cell output
* Info * Building forcefields. Disable mute_ff_builder to see detailed output.
* Info * System has charge 0.0 and multiplicity 1. Provide correct values if this is wrong.
* Info * Building MM systems for the transition state guess. Disable mute_sys_builder to see detailed output.
* Info * Saving reactant and product forcefield as json to ts_data
* Info * Saving systems as xml to ts_data
Starting initial MM scan
MM parameters:
mm steps: 1000
mm temperature: 600 K
mm step size: 0.001 ps
folder name: ts_data
saving mm traj: 0
Lambda | E_int | E1 | E2 | Em
---------------------------------------------
0.00 -85.0 -91.0 540.5 -91.0
0.05 -77.1 -86.8 503.2 -57.3
0.10 -70.2 -75.5 468.9 -21.1
0.15 -63.9 -58.9 436.1 15.3
0.20 -58.3 -38.3 403.5 50.1
0.25 -52.9 -14.4 371.1 82.0
0.30 -47.9 12.6 338.0 110.2
0.35 -43.3 42.2 304.0 133.8
0.40 -39.7 75.1 268.0 152.2
0.45 -36.1 112.3 228.9 164.8
0.50 -33.6 166.6 176.0 171.3
0.55 -36.4 231.2 115.7 167.6
0.60 -42.2 280.7 72.7 155.9
0.65 -49.6 324.5 37.2 137.8
0.70 -58.8 365.1 6.5 114.1
0.75 -68.6 403.0 -19.9 85.8
0.80 -79.6 438.8 -42.6 53.7
0.85 -91.2 473.9 -61.5 18.8
0.90 -103.8 508.1 -76.3 -17.8
0.95 -117.1 542.8 -86.0 -54.5
1.00 -131.7 577.9 -89.4 -89.4
* Info * Found highest MM E: 171.292 at Lammba: 0.5.
* Info * Saving results to ts_results.h5
Starting SCF scan
Lambda | SCF Energy | Difference
----------------------------------------
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
Cell In[2], line 12
9 pro1.set_charge(-1)
11 tsguesser = vlx.TransitionStateGuesser()
---> 12 results = tsguesser.find_TS([rea1, rea2], [pro1, pro2])
13 tsguesser.show_results(results)
File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/tsguesser.py:150, in TransitionStateGuesser.find_TS(self, reactant, product, scf, **ff_kwargs)
147 # Scan SCF
148 if scf:
149 # assert False, 'Not implemented yet'
--> 150 self.scan_scf(self.results)
152 return self.molecule, self.results
File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/tsguesser.py:491, in TransitionStateGuesser.scan_scf(self, results)
489 for i, l in enumerate(self.lambda_vec):
490 geom = structures[i]
--> 491 scf_E = self._get_scf_energy(geom)
492 if i == 0:
493 ref = scf_E
File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/tsguesser.py:533, in TransitionStateGuesser._get_scf_energy(self, positions)
531 self.scf_drv.ostream.mute()
532 basis = MolecularBasis.read(self.molecule, self.scf_basis)
--> 533 scf_results = self.scf_drv.compute(self.molecule, basis)
535 if not self.scf_drv.is_converged:
536 self.ostream.print_warning(
537 "SCF did not converge, increasing convergence threshold to 1.0e-4 and maximum itterations to 200."
538 )
File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/scfdriver.py:575, in ScfDriver.compute(self, molecule, ao_basis, min_basis)
572 self.acc_type = 'DIIS'
574 # check molecule
--> 575 molecule_sanity_check(molecule)
577 # check dft setup
578 dft_sanity_check(self, 'compute')
File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/sanitychecks.py:52, in molecule_sanity_check(mol)
44 def molecule_sanity_check(mol):
45 """
46 Checks molecule for charge/multiplicity combination and geometry.
47
48 :param mol:
49 The molecule.
50 """
---> 52 assert_msg_critical(
53 mol.check_multiplicity(),
54 'Molecule: Incompatible multiplicity and number of electrons')
55 assert_msg_critical(mol.check_proximity(0.1), 'Molecule: Atoms too close')
File ~/miniconda3/envs/echem/lib/python3.12/site-packages/veloxchem/errorhandler.py:50, in assert_msg_critical(condition, msg)
40 """
41 Asserts that the condition is true. Otherwise terminates the program with
42 an error message.
(...) 47 The error message.
48 """
49 if __debug__ and MPI.COMM_WORLD.Get_size() == 1:
---> 50 assert condition, msg
51 else:
52 if not condition:
AssertionError: Molecule: Incompatible multiplicity and number of electrons
To better understand the workflow and it’s limitations, let’s go through it step by step.
Step by step calculation of \(S_N 2\) 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 * Building MM systems for the transition state guess. Disable mute_ff_build to see detailed output.
* Info * Saving reactant and product forcefield as json to ts_data_1758618280
* Info * Saving systems as xml to ts_data_1758618280/systems
The optimised structures can be inspected as follows:
tsguesser.reactant.molecule.show()
tsguesser.product.molecule.show()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
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 \(E1\) decreases (or \(E2\) 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.
mm_results = tsguesser.scan_mm()
Starting MM scan for lambda values [0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0]
MM parameters:
MD steps: 1000
MD temperature: 600 K
MD step size: 0.001 ps
folder name: ts_data_1758618280
saving MD traj: False
Lambda | E1 | E2 | V
---------------------------------------------
0.00 -90.2 544.5 -90.2
0.05 -86.0 506.4 -56.4
0.10 -74.7 471.9 -20.0
0.15 -58.2 438.8 16.4
0.20 -37.6 406.3 51.1
0.25 -13.7 374.0 83.2
0.30 13.0 340.8 111.4
0.35 42.8 306.4 135.1
0.40 75.7 270.3 153.5
0.45 112.6 231.6 166.2
0.50 176.5 169.8 173.1
0.55 232.1 117.2 168.9
0.60 281.0 74.8 157.3
0.65 325.1 39.0 139.1
0.70 365.2 8.4 115.5
0.75 403.0 -18.1 87.2
0.80 438.6 -40.7 55.1
0.85 473.7 -59.7 20.3
0.90 507.9 -74.4 -16.2
0.95 542.1 -84.2 -52.8
1.00 578.3 -87.7 -87.7
* Info * Found highest MM E: 173.114 at Lammba: 0.5.
* Info * Saving results to ts_results.h5
***********
* Warning * File ts_results.h5 already exists. Overwriting it.
***********
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_scf()
* Info * Disable mute_scf to see detailed output.
Starting SCF scan
SCF parameters:
basis: def2-svp
DFT xc fun: b3lyp
Lambda | SCF Energy | Difference
----------------------------------------
0.00 -8173664.35916 0.000
0.05 -8173662.82513 1.534
0.10 -8173654.75778 9.601
0.15 -8173643.81409 20.545
0.20 -8173631.75058 32.609
0.25 -8173619.87759 44.482
0.30 -8173608.84388 55.515
0.35 -8173598.97950 65.380
0.40 -8173590.71415 73.645
0.45 -8173584.49864 79.861
0.50 -8173582.89609 81.463
0.55 -8173582.69938 81.660
0.60 -8173583.52124 80.838
0.65 -8173585.82457 78.535
0.70 -8173589.56852 74.791
0.75 -8173594.39379 69.965
0.80 -8173600.14252 64.217
0.85 -8173605.86392 58.495
0.90 -8173611.16007 53.199
0.95 -8173614.80039 49.559
1.00 -8173614.02489 50.334
* Info * Found highest SCF E: -8173582.699 at Lambda: 0.55.
* Info * Saving results to ts_results.h5
***********
* Warning * File ts_results.h5 already exists. Overwriting it.
***********
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()
3Dmol.js failed to load for some reason. Please check your browser console for error messages.
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 \(\lambda\) 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 \(\lambda\) 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.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 * Building MM systems for the transition state guess. Disable mute_ff_build to see detailed output.
* Info * Saving reactant and product forcefield as json to ts_data_1758621028
* Info * Saving systems as xml to ts_data_1758621028/systems
Starting MM scan for lambda values [0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0]
MM parameters:
MD steps: 1000
MD temperature: 600 K
MD step size: 0.001 ps
folder name: ts_data_1758621028
saving MD traj: False
Lambda | E1 | E2 | V
---------------------------------------------
0.00 -127.7 632.1 -127.7
0.05 -122.1 581.1 -87.0
0.10 -108.7 536.8 -44.2
0.15 -90.5 491.1 -3.3
0.20 -66.6 450.3 36.8
0.25 -38.7 409.9 73.4
0.30 -7.4 369.5 105.6
0.35 27.4 328.3 132.7
0.40 66.1 285.5 153.8
0.45 109.8 239.4 168.1
0.50 159.6 190.3 175.0
0.55 221.8 134.3 173.7
0.60 289.7 77.7 162.5
0.65 340.7 36.8 143.2
0.70 388.6 1.0 117.3
0.75 434.7 -30.2 86.1
0.80 478.5 -56.7 50.3
0.85 523.0 -79.2 11.1
0.90 567.2 -96.6 -30.2
0.95 612.8 -108.5 -72.4
1.00 660.9 -109.1 -109.1
* Info * Found peak MM E: 174.951 at Lambda: 0.5.
* Info * Doing conformer search from Lambda: 0.35 to Lambda: 0.65.
Starting MM scan with conformer search for lambda values [0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65]
conf. steps: 10000
conf. snapshots: 10
MD temperature: 600 K
MD step size: 0.001 ps
folder name: ts_data_1758621028
saving MD traj: False
Lambda | E1 | E2 | V | E_int
---------------------------------------------
0.35 8.0 318.0 116.5 3
0.40 46.6 275.7 138.2 2
0.45 89.5 230.4 152.9 3
0.50 145.6 174.8 160.2 3
0.55 213.8 113.7 158.7 3
0.60 271.9 64.6 147.5 3
0.65 324.4 22.9 128.4 5
* Info * Found highest MM E: 160.176 at Lammba: 0.5.
* Info * Saving results to ts_results.h5
***********
* Warning * File ts_results.h5 already exists. Overwriting it.
***********
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],
breaking_bonds={(1, 3)},
reactant_partial_charges=rea_charges,
product_partial_charges=pro_charges,
)
* Info * Building forcefields. Disable mute_ff_builder to see detailed output.
Reaction summary
2 breaking bonds:
ReaType ProType ID - ReaType ProType ID
c3 c2 14 - hc ho 31
oh oh 1 - ho ho 3
2 forming bonds:
ReaType ProType ID - ReaType ProType ID
oh oh 1 - hc ho 31
ho ho 3 - o oh 24
* Info * System has charge 0.0 and multiplicity 1. Provide correct values if this is wrong.
* Info * Building MM systems for the transition state guess. Disable mute_ff_build to see detailed output.
* Info * Saving reactant and product forcefield as json to ts_data
* Info * Saving systems as xml to ts_data
The output shows that the program found the correct forming and breaking bonds. We can then proceed with the mm scan.
mm_results = tsguesser.scan_mm()
tsguesser.show_results(mm_results)
Starting initial MM scan
MM parameters:
mm steps: 1000
mm temperature: 600 K
mm step size: 0.001 ps
folder name: ts_data
saving mm traj: 0
Lambda | E1 | E2 | V
---------------------------------------------
0.00 -64.6 1092.9 -64.6
0.05 -53.6 845.3 -8.7
0.10 -42.8 785.3 40.0
0.15 -11.3 724.4 99.0
0.20 26.9 666.2 154.7
0.25 69.4 610.7 204.7
0.30 120.3 554.0 250.4
0.35 164.3 506.7 284.1
0.40 221.9 454.0 314.7
0.45 268.2 403.1 328.9
0.50 318.4 352.9 335.7
0.55 373.4 300.8 333.5
0.60 425.4 251.1 320.8
0.65 483.5 199.2 298.7
0.70 566.4 114.4 250.0
0.75 641.9 72.4 214.8
0.80 708.4 25.4 162.0
0.85 768.8 -27.0 92.4
0.90 848.3 -63.7 27.5
0.95 974.7 -90.3 -37.0
1.00 1176.5 -103.3 -103.3
* Info * Found highest MM E: 335.652 at Lammba: 0.5.
* Info * Saving results to ts_results.h5
***********
* Warning * File ts_results.h5 already exists. Overwriting it.
***********
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.
water = vlx.Molecule.read_smiles('O')
rea = vlx.Molecule.read_smiles('C(C(=O)O)([O-])')
pro = vlx.Molecule.read_smiles('C(=C(O)[O-])O')
tsguesser = vlx.TransitionStateGuesser()
tsguesser.mute_ff_build = False
tsguesser.ffbuilder.optimize_mol = False # If True, does an xtb optimization of every provided molecule object before reparameterisation.
tsguesser.ffbuilder.reparameterize_bonds = False # If True, reparameterizes unknown bond force constants with the Seminario method.
tsguesser.ffbuilder.reparameterize_angles = False # If True, reparameterizes unknown angle force constants with the Seminario method.
tsguesser.ffbuilder.optimize_ff = True # If True, does an mm optimization of the combined reactant and product after reparameterisation.
tsguesser.ffbuilder.mm_opt_constrain_bonds = True # If true, constrains nonbonded interactions of atoms involved in breaking/forming bonds during the mm optimization.
tsguesser.ffbuilder.water_model = 'spce' # Water model to use in case a water molecule is provided
tsguesser.ffbuilder.product_mapping = None # Optional mapping of product atoms to reactant atoms
tsguesser.ffbuilder.mute_scf = False # If True, mutes the scf output during xtb, resp and hessian calculations
tsguesser.ffbuilder.skip_reaction_matching = False # If True, skips the matching of product atoms to the reactant and assumes the order of the given structures to be correct
tsguesser.ffbuilder.hessian_basis = 'def2-SV_P_' # Basis set for hessian calculation. For now cannot contain f-functions
tsguesser.ffbuilder.hessian_xc_fun = 'B3LYP' # xc functional for hessian calculation
tsguesser.lambda_vec = [
0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65,
0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1
] # Lambda values to scan
tsguesser.scf_xcfun = "b3lyp" # xc functional for scf calculations
tsguesser.scf_basis = 'def2-svp' # Basis set for scf calculations
tsguesser.mute_scf = True # If True, mutes the scf output during scf calculations
tsguesser.mute_ff_build = True # If True, mutes the output from the force field builder
tsguesser.mm_temperature = 600 # Temperature for mm md runs in K
tsguesser.mm_steps = 1000 # Number of steps for mm md runs
tsguesser.mm_step_size = 0.001 # Step size for mm md runs in ps
tsguesser.save_mm_traj = False # If True, saves the mm md trajectory as an xyz file in the ts_data folder
tsguesser.scf_drv = None # Custom scf driver to be used for all scf calculations
tsguesser.folder_name = 'ts_data' # Folder name for all ts data
tsguesser.force_conformer_search = False # If True, forces a conformer search on all lambda values
tsguesser.discont_conformer_search = True # If False, performs a conformer search where E2 is increasing or E1 is decreasing, which indicates a discontinuity in the PES
tsguesser.peak_conformer_search = False # If True, performs a conformer search at the peak of the PES
tsguesser.peak_conformer_search_range = 1 # Range of lambda values around the peak to perform the conformer search
tsguesser.conformer_steps = 10000 # Number of steps for the conformer search md runs
tsguesser.conformer_snapshots = 10 # Number of snapshots to take during the conformer search
tsguesser.results_file = 'ts_results.h5' # File to save the results
tsguesser.find_TS(
[water, rea], # reactants
[water, pro], # products
# If an scf calculation should be done. Defaults to True
scf=True,
# Set of tuples with indices (1-based) of bonds that have to break and are not allowed to re-connect
breaking_bonds={(1, 3)},
# Optional partial charges for the reactant. Should be a list of list of floats matching the list of reactants
reactant_partial_charges=None,
# Optional partial charges for the product. Should be a list of list of floats matching the list of products
product_partial_charges=None,
# Optional hessians for the reactant. Should be a list of np.arrays matching the list of reactants
reactant_hessians=None,
# Optional hessians for the product. Should be a list of np.arrays matching the list of products
product_hessians=None,
# Total reactant multiplicity to override override the guess
reactant_total_multiplicity=-1,
# Total product multiplicity to override override the guess
product_total_multiplicity=-1,
)