Dear Anthony,

I dont think the problem is related to the starting geometry. It should in general not matter where you put the molecules in the beginning, since you want to compare the energies of three well defined points. These are (1) the isolated species, (2) the fully optimized product and (3) the optimized transition state. In your case, it seem that 1 2 and 3 are connected by a barrierless proton-transfer and hence the optimization converges directly to the product minimum, no matter where you start your optimization.

The underlying problem is that you are trying to compute an ionic reaction that occurs in water without any solvent model. Hydronium ions only exist in bulk water where they are stabilized by a large number of additional water molecules that allows the positive charge to be distributed. In gas phase (your calculation), a sole hydronium ion is very high in energy (unstable) and hence will loose the additional proton as soon as you provide it with a base that allows for a better delocalization of the surplus charge, which is what you are observing.

Try using a solvent model like the PCM and some additional water molecules, which will presumably resovle the problem and provide at least a qualitatively correct description. Nevertheless, to get somewhat close to a quantitative picture, you will need a much better model for the environment that stabilizes the hydronium ion (e.g. QM/MM, preferrably polarizable) as well as a much bigger basis set (def2-TZVP, cc-pVTZ).

Regards,

Jan

First picture is the AFTER

Second picture is the BEFORE

Dear Anthony,

The hydronium ion (H3O+) should be placed at least 4 angstrom from the closest atom in the organic reactant. Otherwise, you will get the product upon optimization.

Hoping this will be helpful,

Rafik

Dear Rafik,

I'm currently running a modredundant scan geometry of bond distance between the organic carbonyl oxygen and the hydronium ion, displacing the ion greater than 4 angstrom. I hope to see the energy drop so I have an idea where I should be optimising as such a distance. 

If this work, I will be extremely interested to learn how a transition point can be found from the displacement of the hydrogen atom off the ion. I'll put my results here. 

Thanks again.

Anthony

Dear Rafik,

Do you have any suggestions or tips on how to optimise the organic with the ion more than 4 angstroms away? On each step of the optimisation algorithm, the ion is moved, seemingly at random, around the molecule until it final coordinates with a very polar atom and thus loses one of the hydrogens. I have tried putting the hydronium ion 5.5 angstroms away and directly opposite the carbonyl oxygen.

Thanks again, Anthony

Dear Anthony,

In the attached Figure 2, I noticed that the hydronium ion is placed in the right side of the organic compound which makes it easy for the proton to be transferred to the carbonyl oxygen. It is suggested that you place it at the left side of the organic compound in a distance of 4 A.

If this approach will succeed, then you can record the energy as the global minimum (GM) energy, then you can start with the reaction coordinate. 

Rafik

Dear Anthony,

I dont think the problem is related to the starting geometry. It should in general not matter where you put the molecules in the beginning, since you want to compare the energies of three well defined points. These are (1) the isolated species, (2) the fully optimized product and (3) the optimized transition state. In your case, it seem that 1 2 and 3 are connected by a barrierless proton-transfer and hence the optimization converges directly to the product minimum, no matter where you start your optimization.

The underlying problem is that you are trying to compute an ionic reaction that occurs in water without any solvent model. Hydronium ions only exist in bulk water where they are stabilized by a large number of additional water molecules that allows the positive charge to be distributed. In gas phase (your calculation), a sole hydronium ion is very high in energy (unstable) and hence will loose the additional proton as soon as you provide it with a base that allows for a better delocalization of the surplus charge, which is what you are observing.

Try using a solvent model like the PCM and some additional water molecules, which will presumably resovle the problem and provide at least a qualitatively correct description. Nevertheless, to get somewhat close to a quantitative picture, you will need a much better model for the environment that stabilizes the hydronium ion (e.g. QM/MM, preferrably polarizable) as well as a much bigger basis set (def2-TZVP, cc-pVTZ).

Regards,

Jan

Many thanks, Jan. That was a brilliant answer. 

I suspected that it was to do with the distribution of the charge on hydronium ion. I hadn't thought of it from the perspective, though, of the vacuum environment in which I have it in. 

Sadly, at this stage along the reaction coordinate, I have been using gas phase. I don't particularly want to go back and recalculate for a PCM model. I take the approach of building up an explicit solvent environment around the region of interest. 

Many thanks again

I have done this kind of calculations more than twenty years. You have to find the GM in a gas phase and NOT in the presence of a solvent. In the GM, the two reactants should be far away each from the other such that no interactions are exist between them.

I have published more than 40 papers in very good journals using this approach.

Your calculations should be in the gas phase and in NOT in a solvent using PCM.

H3O+ exists in the gas phase and if some body thinks that it is not he has to go to search the literature.

If you use PCM you will never get the energies in the gas phase.

It is worth noting, that the solvation models such as PCM are relatively new. Therefore, how our great scientists were able to handle this kind of calculations without using PCM in the years before the invention of such models?

Rafik

Hi Rafik,

I was able to optimise the ion + organic in gas phase at a distance of around 8 angstrom. The hydronium ion is meant to lose a hydrogen atom which protonates the oxygen of the carbonyl group as the nitrogen attack the carbonyl carbon. The nitrogen will not give up its hydrogen, making the nitrogen + charged. 

How am I to set up the transition state of this reaction if my hydronium ion is so far from the carbonyl. 

Many thanks

Hi Anthony,

You start an approach of the H towards the O of the CO from 8 A in increment of 1 A for 4-5 cycles (until reaching 3-4 A), then you run another reaction coordinate with a smaller increment 0.1 A.

Rafik

Hi Rafik,

Is there a way to automate the switch from 1 angstrom steps to 0.1 angstrom steps? Or will this require two separate calculations? 

Hi Anthony,

May be there is, but I do not know how. I suggest to do it in two reaction coordinate in order to make sure that at distance of 3-4 A you still have the two separated reactants.

Rafik

Anthony,

I suggest you to avoid using explicit environments, even if it seems like the lesser of two evils right now. In my experience, this becomes very messy due to a number of reasons and, at least in my experience, never led to any usable results. Firstly, the PES of explicitly solvated systems is usually flat and riddled with local minima, which makes it very difficult and time intensive to converge optimizations. Secondly, the results you obtain even if you manage to converge the calculations are of very limited meaning since what you find is only one of many solutions. It will strongly depend on the conformation of the solvent shell as well as on number of water molecules and unfortunately, the results hardly ever converge w.r.t. this number.

Thirdly, and this is usually the most problematic thing, the solvent shell will most certainly differ between isolated species, product and transition state. Lets say two water molecules turn around adding one or two hydrogen bonds and lets say this happens between educt and product. This will make up for tens of kJ/mol and you wont be able to distinguish between the actual reaction energy and the contribution of the solvent shell. Hence, you will need some kind of sampling, i.e.  averaging over tens or hundreds of calculations using different explicit solvation shells to account for this (you will see why if you compare the energetic outcome for two seemingly similar explicit model systems).

Even if it seems like a ton of work now, I promise you the explicit route is more work and its not only the easiest but also the most reliable solution to recalculation everything with a PCM right now. The gas-phase geometries you already have will most certainly provide a very reasonable starting point for the optimizations with a PCM. It might nevertheless be a good idea to add at least a few water molecules at the most polar hydrogen bond acceptors and donors since PCMs aren't very good at describing these. Its just that with a PCM that describes the bulk electrostatics, the results may actually converge w.r.t. to a number of explicit solvent molecules you can afford to calculate :)

Good luck!

Jan

Are you interested in this reaction in  solution or in the gas phase? If it's the latter, then don't use any continuum solvent model or explicit water molecules. If, as with most protein work, you're modeling something in solution, then you should take Jan's advice and include both. The dielectric constant is radically different in each case and will make a huge difference in the forces felt by partial charges such as oxygen lone pairs and double bonds and protons on hydronium. H3O+ is actually a pretty lousy model for a solvated proton (see the literature for a LOT of discussion of this fact) and so including explicit water molecules is important if you want to get the physics remotely right. Electronic structure calculations don't deal well with very localized positive charges, so including explicit water molecules that will allow the proton to delocalize is important. How do you decide how many water molecules to use? Well, we have information about the thermodynamics of a protonated carbonyl versus a non-protonated through experimental pKa values. Find the number of explicit water molecules that give you the appropriate difference in free energies between those states and that"s just about the best you can do (within a certain choice of functional and basis set) to ensure that you get the kinetics right as well. 

If you don't want to go back and redo everything with a continuum solvent model and with explicit waters delocalizing positive charges, then it's unlikely that you'll get reliable results for solution phase chemistry. That might be fine in some cases, but expect to hear about it from a reviewer if you send it off to publication. 

Hi,

It is simply use order cleanup + at build the molecules and sign track self-consistent field and increase of the cyclic number of energy optimization.