A third question would be:

Whilst forming a reaction mechanism using DFT/QM methods, does the protonation/deprotonation warrant a reactant - transition_state - product sub mechanism? Could one assume that the nitrogens in the imidazole structure (part of a greater struture) are in a flux of protonated/deprotonated due to the solvent environment until something happens to the structure, such as the removal of a neighbouring hydroxyl group?

Dear Anthony,

I'm not sure if I understood well your question but I will try to answer. You can always calculate a deprotonation transition state (although most of the times, this barrier is very small and difficult to find). If you know your reaction conditions, you can take the best brönsted base of your system and compute the deprotonation (it can be the solvent). If the product is much more unstable than the protonated species, you would have recombination, but you can calculate the species separately and sum the energy to compare the stability of your protonated-deprotonated species.

Another option to see if your system is more or less deprotonated is to compute the pKa value, and compare it with the pH used in your reaction, assuming that your species are in equilibrium. To compute it, you need the energy value for the proton in water (your reaction should be in water, but according to your description, that is the case). You can see an example here: ACS Catal., 2017, 7 (3), pp 1712–1719.

I would agree with both of the previous answers: 1. it is often acceptable to assume that non-carbon proton transfers will not be rate-limiting, so knowing the equilibrium constant (easily obtainable from the relative free energies of the reactant and product) is kinetically sufficient. 2. If you want to find it, you need to give the proton to something in your system, and solvent is always a valid choice (unless perhaps you are in an enzyme pocket).  I have found that the choice between manually stepping the N-H bond length while the H is coordinated to a solvent LP is often the easiest way to locate the TS, which can then be optimized by the standard opt=TS method; you can use an IRC, but it is only faster in some cases (in my experience, those involving a reaction coordinate with contributions from several bond stretches/bends).

Calculating pKas can be a tricky business, and if it were me I wouldn't mix solvated-proton energies calculated by others with energies I calculated; I would calculate all the energies myself to take advantage of what error cancellation I could in my computational system.  Also, knowing the position of the equilibrium tells you nothing about the speed of the equilibrated process; that is why it may be sufficient to know that, unless you have specific reasons to suspect otherwise, the proton transfer to/from N will be faster than other processes in your chain.

Thank you, both. There is a lot of food for thought there. 

Dear Ted, you spoke of manually stepping the N-H bond length whilst the H is coordinated to solvent. How did you manage to optimise the proton in the solvent environment to begin this process? Consider, that my system is positively charged (from the proton in the environment as the product or from the proton of -N+(H)=C as the reactant structure).

I am finding that the -N=C nitrogen of my compound is such a strong attractor for the proton positioned within the explicit water solvent region of the PCM environment, that eventually the proton works its way out during the optimisation steps and back onto the nitrogen, resulting in the -N+(H)=C structure. 

I assumed as I wasn't able to optimise a product (protonated water environment with -N=C structure to the compound) that there would be no point doing a modredundant bond scan. However, I am giving it a try, but starting with the proton on the nitrogen and working onto a neighbouring water molecule. I hope I find something from that scan.

It is quite possible that the proton is much less stable on water than on your nitrogen; however, many reactions proceed through steps with unfavorable equilibria, so that is not a problem in and of itself.  It will tend, however, to make the TS late, meaning you may have to step it quite far – I have often experienced having to step into a region where my calculation had two non-overlapping solvent cavities, although the TS has so far always been in the region before cavity separation.  The validity of energies in a two-cavity region is questionable; the correct product energies would be the sum of the solvated species (neutral N and hydronium) calculated separately to avoid BSSE and possible weirdness from two separate cavities that are nevertheless to close to one another that no solvent molecule could be between them.  This sort of computational problem is related to the classical issue of ion separation in a medium, where it was found that there is in fact an energy barrier between the contact ion pair and the solvent-separated ion pair.

I just want to verify that you are conducting a relaxed scan; you might get an approximate answer from a rigid scan before doing a relaxed scan to save computational time (if your system is large and/or floppy), but to find the TS, I strongly recommend a relaxed scan.  I have been quite surprised by the exact geometric course of a reaction on occasion, but a relaxed scan always tells the tale.

I've just started two relaxed bond scans between:

1) increasing distance between the N-H bond directed into the solvent

2) decreasing distance between the N-H...OH2 

Given your very useful description of proton distance into the solvent cavity, I suspect that the increasing distance between N-H (1) will hopefully reveal at least the product I am looking for.  

I understood from your earlier post(s) that you have an explicit water molecule nearby; are you setting it up with a guess at a coordinated geometry (coordinated to the departing proton, that is)?

Yes, it is a guess, but to be very certain that I am being clear the description is as following. I have a charged molecule (part of the molecular is a charged imidazole; -N+(H)=) surrounded by a small solvent layer (eight water molecules) and a PCM environment. I've *always* been able to optimise to this part of the reaction mechanism. The next step, deprotonation of the imidazole region, I've never been able to optimise. So my approaches to the scan bond, I'm scanning by taking the charge (proton) off the imidazole:

1) by increasing the N-H bond distance and just seeing where the H goes (hopefully in an energy well in the solvent).

2) by decreasing the distance between the hydrogen of the N-H of the imidazole with a neighbouring water molecule. I suspect, though, that with strength in the affinity for N...H will simply cause the neighbouring water molecule to "pull" into the N-H bond with no N-H covalent disassociation.

Happy to hear your thoughts on this setup. 

Thanks

My systems tend to be too computationally demanding for any more explicit solvent molecules than are absolutely necessary, so I've never tried just stretching the X-H bond to see what grabs it.  You might very well still arrive at the TS, but with higher intermediate energies than starting with a water molecule hydrogen-bonded to the proton to be removed.  If the water molecules won't stick to the imidazolium H, I guess you have no choice (I would be curious why they don't, though).

If you are using an unrestricted QM method or one that can become unrestricted if it results in a lower electronic energy, I would watch for an increase in higher-spin-order contributions to the solution, as that would be a warning that, in the geometric configuration at that point, the electronic solution is headed toward a homolytic dissociation.  I would be surprised if you had that problem, though.

Thanks for your input. Thankfully, my reaction is closed-shell/restricted, quite simple reactions but just a long cascade of them and a few stages in which I need to juggle with protonation states. 

I'll be sure to tell you how the bond stretch scans go.