TS-search with explicit solvent

Question

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?

Answer 1

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():

    geom = manage_xyz.read_xyz("geom.xyz")
    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,
            ad_idx=0,
            multiplicity=1,  # your system multiplicity
            )

    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,
            addtr=True,    # Adding translations and rotations!
            topology=top,
            )

    nifty.printcool("Building Delocalized Internal Coordinates")
    coord_obj1 = DelocalizedInternalCoordinates.from_options(
            xyz=xyz,
            atoms=atoms,
            addtr = True,  # Adding translations and rotations!
            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))
    manage_xyz.write_xyzs_w_comments('output.xyz',geoms,energies,scale=1.)
    manage_xyz.write_xyz('minimized.xyz',geoms[-1],energies[-1],scale=1.)

if __name__=='__main__':
    main()