# TS-search with explicit solvent

I'm looking for an advice how the best to set up the system and optimize a transition state with explicit solvent.

I believe that solvent molecules are participating in the reaction by specific solvation/coordination and assist in one of the reaction steps. Geometry optimization of the pre-reaction complex converged to a local minima of the system featuring solvent coordination on a different spot (hydrogen bonding).

I do not want to introduce too many solvent molecules, as solvent behavior is extremely dynamic.

Could you, please, recommend possible ways to optimize TS with explicit solvent, or some literature as a reference?

• I edited the grammar to make sense. -- think about this as if you are sending a professor a question, you don't want any grammatical errors and you want to be clear as possible. May 19, 2020 at 6:56
• What method are you using to find the TS? May 20, 2020 at 12:50
• I start with single ended Growing String Method (pyGSM) by @CodyAldaz, and then run TS-search. May 20, 2020 at 19:09
• I'll answer this soon. One last question about editing, what do you mean you can't find a "local minima of the solvent assisted transformation"? -- it seems you can find many minima but they just aren't what you want. Could you post a picture? May 20, 2020 at 21:47
• @CodyAldaz, You are right, my wording is still not accurate. I cannot find a Desired local minima, while all calculations fall to any other geometry, always driven by hydrogen bond formation. When I tried to constrain distance from the solvent to three nearest atoms in the molecule (to fix initial structure), SE_GSM started right, but then switched to a wrong direction. May 21, 2020 at 3:18

Explicit solvent can be of importance, like you mention, whenever strong coordination or hydrogen bonding is involved. For example, transition metals often coordinate to solvent and not including them in the calculation would give very poor results. Similarly hydrogen bonding can affect the electrostatics and result in poor properties like reaction paths, pKa, etc.

There are generally two ways to introduce explicit solvent into a simulation:

1. Introduce one or two explicit solvent model in the areas of strong coordination, and use an implicit solvent model for the remaining interactions.
2. Treat all of the solvent explitly, for example with QM/MM. Normally, this involves running some molecular dynamics with periodic boundary conditions, and then optimizing transition states for several snapshots of the dynamics. However, this method introduces so many degrees of freedom that it's no longer possible to apply the harmonic approximation to calculate entropy, and it's difficult to gather enough statistics. An example of this can be found in Houk's Solvent perturbed transition state calculations [1]

Since you don't want to do the latter (I also would not recommend this for a first pass either). Your best bet is to try to coordinate a few solvent and optimize.

Finally, any time you are dealing with intermolecular systems (e.g. reactant and one or more solvent molecules) I recommend the translation and rotation internal coordinate system (TRIC)[2]. This coordinate system is similar to the very popular Delocalized Internal Coordinate (DLC) system but is superior for intermolecular optimization.

This can be accomplished via the package Geometric [3], but since you are already using pyGSM, this can also be accomplished with a pyGSM script[4]. I've attached the optimization script below.

References:

1. Yang, Z., Doubleday, C. & Houk, K. N. QM/MM Protocol for Direct Molecular Dynamics of Chemical Reactions in Solution: The Water-Accelerated Diels-Alder Reaction. J. Chem. Theory Comput. 11, 5606–5612 (2015).
2. Wang, L. P. & Song, C. Geometry optimization made simple with translation and rotation coordinates. J. Chem. Phys. 144, 214108 (2016)
3. https://github.com/leeping/geomeTRIC
4. https://github.com/ZimmermanGroup/pyGSM
from pygsm.level_of_theories.qchem import QChem
from pygsm.potential_energy_surfaces import PES
from pygsm.optimizers import *
from pygsm.wrappers import Molecule
from pygsm.utilities import *
from coordinate_systems import Topology,PrimitiveInternalCoordinates,DelocalizedInternalCoordinates

def main():

xyz = manage_xyz.xyz_to_np(geom)

nifty.printcool(" Building the LOT")
lot = QChem.from_options(
lot_inp_file="qstart", # a Q-Chem input file with only the $$rem$$ arguments
geom=geom,
)

nifty.printcool(" Building the PES")
pes = PES.from_options(
lot=lot,
)

nifty.printcool("Building the topology")
atom_symbols  = manage_xyz.get_atoms(geom)
ELEMENT_TABLE = elements.ElementData()
atoms = [ELEMENT_TABLE.from_symbol(atom) for atom in atom_symbols]
top = Topology.build_topology(xyz,atoms)

nifty.printcool("Building Primitive Internal Coordinates")
p1 = PrimitiveInternalCoordinates.from_options(
xyz=xyz,
atoms=atoms,
topology=top,
)

nifty.printcool("Building Delocalized Internal Coordinates")
coord_obj1 = DelocalizedInternalCoordinates.from_options(
xyz=xyz,
atoms=atoms,
primitives=p1,
)

nifty.printcool("Building Molecule")
reactant = Molecule.from_options(
geom=geom,
PES=pes,
coord_obj = coord_obj1,
Form_Hessian=True,
)

print(" Done creating molecule")
optimizer =  eigenvector_follow.from_options(Linesearch='backtrack',OPTTHRESH=0.0005,DMAX=0.5,abs_max_step=0.5,conv_Ediff=0.1)

print("initial energy is {:5.4f}".format(reactant.energy))
geoms,energies = optimizer.optimize(
molecule=reactant,
refE=reactant.energy,
opt_steps=500,
verbose=True,
)

print("Final energy is {:5.4f}".format(reactant.energy))