What to do with (large) imaginary frequencies for constrained minimum structures?

Question

I am performing DFT calculations using ORCA 4.0.1 on an enzyme active site model. The model contains 89 atoms including the substrate (see Animation 1), five of which are fixed in space (the spherical atoms in Animation 1). I am using dispersion-corrected B3LYP with 6-31g(d,p) basis set, with the "RIJCOSX" approximation, and with a CPCM model for the solvent (dielectric constant $\epsilon = 4)$. From what I have seen in the literature, this method is often applied to enzyme models with success, and this is the reason for using it. Once I start computing reaction paths, I plan to try different functionals and basis sets to see if this give drastically different results. My geometry optimization tolerances are

Energy change: 5e-6;
Max gradient: 2e-3;
MRS gradient: 5e-4;
Max displacement: 4e-3;
RMS displacement: 2e-3

The two gradient tolerances have been loosened slightly in order to converge the calculations. By using the default gradient tolerances, the calculations simply do not, I have found (especially the maximum gradient). Additionally, it has been absolutely necessary to compute a better guess for the initial Hessian for the optimization, done at the B3LYP/3-21G-rijcosx level.

The fixed atoms are kept fixed to maintain the orientation of the active site residues. Only atoms where the side chain or backbone has been truncated are kept fixed. Without fixing these atoms, the residues will most likely deviate way too much from the crystal structure, and no longer represent the enzyme environment we try to simulate.

Now, due to the constrained atoms, I expected imaginary frequencies to occur. However, for one of my presumably converged structures, I have quite strong imaginary frequencies. The list of frequencies start as follows: -220.645934, -69.117105, -65.597557, -61.314719, -48.720497, -40.439949, -32.562415, -21.395791, -16.637113, 6.126145. The two strongest negative frequencies are mainly a sigma bond rotation of a methyl group (of which the carbon is fixed in space during the optimization).

It is difficult to find any good literature on what to do with and how to interpret such negative frequencies in constrained optimizations. Is it "safe" to disregard the frequencies that mainly involved fixed atoms? Should I converge my geometry tighter to come closer to the true minimum? That's not so easy, I find, as the molecule is quite flexible. It took me quite some time to come up with a method that converges, which involve computing a good guess for the initial Hessian used in the optimization (computed at the B3lYP/3-21G level). I hope the imaginary frequencies are artefacts of the constraints, and that I can safely disregard them. However, they seem to have quite higher frequencies than what other people report in the literature (many say that their imaginary frequencies are around 20 cm^-1).

This has to be a common challenge for people doing electronic structure calculations on large and constrained systems, but very few address this issue specifically in their publications.

I did read through this question, but it addresses not quite the same issue that I face.

Animation 1. The active site model. A total of 89 atoms.

Animation 2. Imaginary mode with frequency -220.645934 cm$^{-1}$.

What I am trying right now

I have three jobs waiting in the queue right now:

1. I have made a small rotation of the methyl groups strongly involved in the strong imaginary mode

2. I have started another geometry optimization of the structure given above (exhibiting the imaginary frequencies), but by reading the computed Hessian (the "exact" one, from same level of theory as optimization)

3. I have restarted an identical optimization, but by computing the "exact" Hessian from the start, giving a much better guess of which direction to go in.

I can add here later what these attempts result in.

Answer 1

The first I noticed is you said:

The two gradient tolerances have been loosened slightly in order to converge the calculations. By using the default gradient tolerances, the calculations simply do not, I have found (especially the maximum gradient).

At first sight, the structure of your molecule should not leads to serious problems, but the above is a bad signal, and turns the results unreliable (at least doubtful). The constrains must be considered in the calculation of vibrational frequencies for they make sense.

Maybe no one can give you a foolproof steps series to follow. My suggestions are below.

What to do/check

• If comments are right, and Orca does not care about fixing atoms for the calculation of the vibrational frequencies, then you need to change the software for this task (because in this case results are not meaningful) (This also answer the question "Is it "safe" to disregard the frequencies that mainly involved fixed atoms?"). The simpler way to check this is by looking at the number of freedom degrees reported (or deduced form the output). I was referring directly to vibrational f.d. because low frequency modes may be treated like rotations. I doubt that Orca lacks of this basic functionality. In such case I bet that this is described their User Manual.

• The step you are doing are reasonable. My guess is that maybe you are far from a true minimum, in some unstable geometry, and that is why you have problems with accuracy changes. You said: "Should I converge my geometry tighter to come closer to the true minimum? ". Yes, it is always advisable to get good frequencies. So I would use not only tighter geometry convergence criteria but the TIGHT convergence option for the SCF part.

• In the same lines, you can try using modern semi-empirical methods and then make the fine tuning with DFT. Even using a simpler DFT functional.

• I agree with comments that you are not using the best basis set out there, but I think that this is not a problem here.

• If you want to discard weird behaviors you can try with another functional, but I do not think that this is the problem (because it is known that B3LYP performs well in compounds with the elements involved here).

In brief. I would check if Orca compute correctly the frequencies, if not: change the software. If yes: discard those results, try to converge the structure using cheaper methods (if you are lack of computational resources) from different initial conditions, store those results and start new geometry optimizations and vibrational frequencies calculations starting from those previous converged states. The last step using the tight convergence criteria for SCF and normal convergence criteria for geometries, and using the target basis set (for energy and its derivatives).

Good luck

Note:

Once I start computing reaction paths, I plan to try different functionals and basis sets to see if this give drastically different results.

Why do you expect that something like this would happen? You can save time by using a little large basis set and a better functional once you get the B3LYP/6-31g(d,p) results. It should not take too much time.

---

EDIT DUE TO COMMENTS

This answers generated some discussion about B3LYP or not B3LYP?

In that respect, my viewpoint is as follows:

• B3LYP is not the most accurate functional (there should be no doubt about that).
• It is reasonably good for geometries in systems like yours. Errors in geometries can be larger due to small basis set than due to B3LYP, or at least of the same order.
• While there are functional with better results for geometries, you also must take care of how large basis sets those functional needs to be a better option. I mean, maybe using 5$\zeta$ basis sets they are better, but not using a small basis set. (and maybe not)
• The above discussion can be very long, but in any case, you are using approximations that (in my opinion) are much stronger than the errors due to B3LYP. For example, the permittivity value, or the geometry of the constained parts, etc..
• Also, even if you could perform a FullCI/CBS you shoud notice that the simulated system is different to the real system, so extremely exacts results in the simulated system may no improve your knowledge about reality.
• So, if the functionals are accurate enough, what it is important here is obtaining the corresponding results for the functional/basis set. The functionales can be more or less problematic. B3LYP should not be problematic for your system. My guess is what I said in the answer.
• Because of the above, I think that B3LYP is reasonable and enough for your case.

Answer 2

Large imaginary frequencies, while you are looking for a local minimum, are usually an indicator that your guessed starting geometry is not good. Values larger than $10i~\pu{cm-1}$ are only observed for transition states, and only if you have one of these modes, it is actually meaningful.

When you are fixing atom coordinates, then you are (in a practical sense) employing an external potential. If the program you are using does not account for that constraint, it is not suited for these kinds of optimisations.^note
At the very least, your thermochemistry will be wrong. If your aiming for anything else than a rough guesstimate, then this is not good enough (if you are contempt with that, why bother with frequencies).

In general there are a few things you can do to improve convergence. One of the most important ones and often underrated is the SCF convergence. If you have a sloppy converged SCF, then also gradients will be sloppy, which leads to poor performance of the geometry optimiser. Tighten it to 1E-7 to achieve better results.
Often numerical noise is a problem in calculations, make sure that you are using a tight (tighter, the tightest) grid for DFT calculations. In principle the same as above applies; if it is to coarse, your SCF will converge sloppy, etc.
This should at least let you converge the geometry on default cutoffs (gradient, displacement). If you then have imaginary frequencies, you should use path following algorithms to get rid of these; and if they are not numerical noise, this should not be a problem.

Use dispersion corrections. There is no excuse not using them (except missing parameters, but then use a different method).

From my experience, solvent models can drastically influence the convergence, in many cases it may lead to oscillation. Try optimising without it first, turn it back on later.

You have stated that this is a part of an enzyme. My first guess seeing this large imaginary modes, was that you chose to cut the moieties at the wrong position, substituting more demanding groups with hydrogen atoms. While this is common practise for many computational models, if you have problems with the molecular structure, this is a sign that you cut away something important. Increase your cut-out; don't replace with hydrogen, but with methyl instead (in doubt fix those, too).

Leaving something out of the calculation can lead to all sorts of artefacts. You should have a look at a methodology where you can describe the system in a whole. There a various applications developed just for this case. Something along the lines of DFTB+ might do wonders for your system. Another possibility is using QM/MM, and ORCA is capable of that.
If you want to describe reactions, in solution, sooner, rather than later, you need to look at these options. Aqueous media are especially challenging.

B3LYP is a terrible functional. This is probably more often true than not. It is not only one of the slowest performers (for hybrid functionals), it also has been shown to produce unreliable results (just search B3LYP failure). Starting with a hybrid functional is overkill for this system size in any case. Always use a non-hybrid functional (BP86-D3(BJ), B97D3, M06-L, etc.) first, it saves a lot of time, resources, brain cells. Work your way up Jacob's ladder, not down.
Use def2- basis sets. The close cooperation with Ahlrich's group has many of the great features developed there (re: turbomole) readily implemented in ORCA. This includes the extremely well performing basis set, which do already come with especially prepared auxiliary basis sets for density fitting (RI).

TL;DR: Tighten SCF; tighten grid; use dispersion; use a cheaper functional; use better basis; increase system size (you can fix adjacent moieties).

---

Note: If the program does not correct for these constraints, then the derived gradients won't be good enough to achieve convergence.

Answer 3

I do many calculations on the active site model, but using Gaussian instead. Like what you have also mentioned in the question, the constraint on the truncated atoms caused some imaginary frequencies. However, I was greatly surprised by the -200 cm-1. Try to optimize the structure and calculate the freq at the level B3LYP/6-31**. I would expect much much smaller imaginary there.

update: I thought the calculations were done at the level of B3lYP/3-21G. Sorry about that. Have you found the solution now?