What happened when you used the default value? What software are you using?

As both Eugene and McKay have hinted at, you are going to need to do *a lot* of experimentation and as McKay already pointed out, it does sound like you have hit a problem that tightening the convergence isn't going to help you anymore...

There are a few things besides this you might want to consider which might help your calculation. In your case, the first thing I would recommend is using another co-ordinate system (e.g.,  internal co-ordinates) rather than inputting the atoms as Cartesian co-ords (which I assume you are doing?). This often helps avoid these looping convergence problems because the software is then minimizing bond angles and distances rather than moving atoms X, Y and Z. This is also more efficient and will optimize in fewer steps... some software will automatically convert your cartesians into internals, or ...if you don't want to do all of that, sometimes just changing the position of an atom or two by a small amount and re-running the code can completely remove such an issue... 

However, (pretty) please bear in mind that - despite its wide use - the level of theory you have chosen more than likely isn't going to give you results that accurately portray that of your molecule nor your band-gap very well... in such large molecules I would expect dispersion forces would make a large difference (which can nowadays be included easily, but again take longer), and also there are no polarization functions so probably the hydrogen-bonding will not be properly accounted for, so the structure may be completely off, etc, etc... now, you might be able to get away with this if you are looking for a general trend, but I would highly recommend a slightly higher level (longer) calculation might be worth doing if you are able to...

Only method of tests and mistakes can help You, Dumitru.

It depends strongly on your molecule.

The fact that the there is a jump between 0.002 and 0.005 that repeats suggests you may be in a loop. Look at the intermediate geometries during the optimization. If there is an infinite loop, you will see the atoms hopping back and forth to the same positions. If this happens, you can tweak the geometry slightly and run it again, or you can try using a smaller stepsize for atomic movements during the optimization.

That would indeed explain why the tighter convergence criteria is taking so long. In general, tighter convergence will be more accurate, but I must agree with Dr. Kislyakov that you should test various methods to see which best duplicates what you see in the lab or in literature. I would start with a smaller representative system when testing different method parameters.

What happened when you used the default value? What software are you using?

As both Eugene and McKay have hinted at, you are going to need to do *a lot* of experimentation and as McKay already pointed out, it does sound like you have hit a problem that tightening the convergence isn't going to help you anymore...

There are a few things besides this you might want to consider which might help your calculation. In your case, the first thing I would recommend is using another co-ordinate system (e.g.,  internal co-ordinates) rather than inputting the atoms as Cartesian co-ords (which I assume you are doing?). This often helps avoid these looping convergence problems because the software is then minimizing bond angles and distances rather than moving atoms X, Y and Z. This is also more efficient and will optimize in fewer steps... some software will automatically convert your cartesians into internals, or ...if you don't want to do all of that, sometimes just changing the position of an atom or two by a small amount and re-running the code can completely remove such an issue... 

However, (pretty) please bear in mind that - despite its wide use - the level of theory you have chosen more than likely isn't going to give you results that accurately portray that of your molecule nor your band-gap very well... in such large molecules I would expect dispersion forces would make a large difference (which can nowadays be included easily, but again take longer), and also there are no polarization functions so probably the hydrogen-bonding will not be properly accounted for, so the structure may be completely off, etc, etc... now, you might be able to get away with this if you are looking for a general trend, but I would highly recommend a slightly higher level (longer) calculation might be worth doing if you are able to...

Dear all, thank you very much for suggestions.

The equilibrium geometry was located, but with a memory leak. And the results are quit strange in the sense of molecular orbital picture. I am working with Gamess.

Will re-run with a higher basis set with default convergence tolerance. I suppose I need both d and p polarisation functions for heavy atoms and hydrogen? And what about diffuse functions? I think using both polarisation and diffuse functions will be too much for my machine.

The think is that I am not really a theoretician, but rather a synthetic and physical chemist. These calculations I need for rationalizing the design of molecular systems to further synthesize. So Yes, I am mostly looking at the general trend by introducing some small changes in the system.

The "desert island" basis set is 6-31+G(d,p), in the sense that this is a very good basis set for optimizations. Of course, you need to check if your system is not too big for this basis set. For inorganic and bioinorganic systems, I stick with with 6-31G(d,p), and many in the literature do as well. If hydrogens do not matter that much (no medium range interactions), you can go down to 6-31G(d). For orbitals and energies, you should go to a larger triple zeta basis set at least, but only for single point energy calculations. Now, related to convergence criteria, I think you should just stick to the standard given ones. Tighter convergence is usually used to get a better value for entropy, and I don't think you are calculating that, are you?

Regarding the dispersion forces that Dr. Bennett mentioned, you can use Grimme's D3 dispersion scheme. I'm not sure it is implemented in GAMESS, but you can download Grimme's standalone software and add the resulting D3 corrections on the final energy. It won't affect orbital energies this way, though.

http://www.thch.uni-bonn.de/tc/index.php?section=downloads&subsection=DFT-D3&lang=english

Barring any methodological issues in the functional and basis set, increasing convergence criteria without increasing the grid is likely leading to the problem observed. A large grid is needed when there are very soft modes in large molecules, such as methyl rotations or other very soft torsions.

The "desert island" and 0.01 convergence tolerance worked well, and not so long as I expected.

Thank you everybody for suggestions.

I think everyone else already answered everything - but just my 2 cents...

The take-home message : "Use the highest level of theory, the largest basis set, and the largest grid size that you can afford to run on your system"

The polarization functions (d,p) are 'probably' more important than adding the diffuse functions if you absolutely have to choose, and will likely be less problematic in terms of additional convergence time.

I 'personally' don't have a lot of confidence in calculations that are performed using basis sets below triple-zeta quality, which in your case would mean going a step further still, to 6-311G(d,p). BUT, this may be impractical for your size molecule and CPU power. Also, bear in mind perhaps that's just my own prejudice from my own experiences... similarly, for what you want, the 0.01 convergence tolerance is 'probably' fine, if you ever want to predict the vibrational frequencies you are going to need to go a few orders of magnitude lower. But for UV-VIS this is probably okay. However, I will warn you that for such a calculation be aware that you may only be at a local minimum, not necessarily the global minimum... ... in order to check you are at a local minimum, you should run a frequency (Hessian) on your converged structure and make sure there are no imaginary frequencies. If there are, you will need to be more strict with your convergence criteria (or increase your grid size... or use a better level of theory).

The D3 dispersion scheme is implemented in GAMESS, although I think CAMB3LYP was designed to add dispersion to B3LYP and from what I call is a better choice. Try it. See how much it slows things down. Is it worth it?

btw, GAMESS can also convert your coordinates into internals, there is a keyword in the manual how and its as simple as adding that into your input file. 

As Bradley wrote, increasing the grid size is always recommended if you have any convergence issues but will increase your run-time. In my experience this is extremely helpful getting tricky molecules to converge (wavefunction), and since the energies are more accurate, also somewhat helpful for geometry optimization - again, as he said, particularly for soft modes - BUT at the 0.01 convergence criterion, I'm not sure it would be worth it here? Your molecule may well have a methyl group hanging out in slightly the wrong place, but would that greatly affect the UV-VIS? Try it. See how much is slows things down. Is it worth it?

I understand you said you are primarily not a computational chemist (neither am I), but it would still be a good exercise (or good practice) to notice that there aren't large deviations in your structures calculated at different theories, or if they are to try and understand which may be the more reliable.

Also, another suggestion is that you can use a lower level of theory to approximately get the geometry of your molecule, and then use a higher level of theory using this geometry as a starting point for your higher level calculations. This can save a lot of time. However, you still may have the local vs. global minimum problem, but to solve that problem for such a molecule, this would become your day job...

Lastly, you can also just use that geometry and run MUCH higher level single-point calculations to more accurately obtain the energies of the HOMO/LUMO needed. You start seeing the notation similar to this (note the double //):

B3LYP/6-311+G(d,p)//B3LYP/6-31(d,p) which states the level of theory used for the energy and geometry, respectively. Since this is a single point energy you can often just run the highest level calculation possible on your system... again, if the lower level of theory is too low, the geometry and hence calculation may be meaningless, but it can serve its purpose of getting more accurate energies, or energy levels.

Dear Chris,

thank you a lot for your "2 cents". Very helpful actually.

When I'm doing an optimization of a molecule I'm following the next steps:

1. Running conformer distribution search with Spartan (AM1 level) and plotting the structure vs energy. From this I choose few "best" conformers and run a short HF/3-21G geometry optimization to be sure I take the global minimum. Sometimes if the energies are really close the second conformer may become first (I mean a lower energy), usually less than 5 kcal/mol.

2. Running B3LYP/3-21G geometry optimization with Gamess, starting from AM1 - Spartan result. If there are two (or more) very similar (as energy) results from previous step, calculating them both.

3. Running a B3LYP/6-31G (now 6-31+G(d,p)) geometry optimization.

4. Here I suppose I could do a B3LYP/6-311+(d,p) single point calculation to get more accurate energy levels.

Does this sound reasonable?

Hi Dumitru,

Short answer: Yes.

    Unfortunately I have no experience with Spartan, nor the accuracy of the AM1 level of theory but from the sound of it you are trying to do what you can...

     Nor am I  the best person to comment on exhaustive searches for the global minimum on the potential energy surface - as you are probably aware, the lowest energy conformer at the HF/3-21G level may not be the lowest energy at higher levels of theory... however, your overall approach seems fairly sound to me but perhaps someone else can comment...

I think the B3LYP/6-31+G(d,p) is a sensible choice given the size of your molecule.

The energy levels will be only very slightly more accurate with the larger basis set but better is better... perhaps what might be more interesting - if you have the time - is trying a few different functionals (e.g., PBE or something) at that minimized geometry to see how differently they calculate the energy levels you are interested in.

You can probably speed up your whole progress by switching to a better performing and more balanced basis set, like def2-SVP or CC-PVDZ. Pople basis sets are antiques and should be used for reproduction only.

Also note that adding polarization and diffuse functions does not necessarily create a better basis set. It merely compensates (or introduces) some errors.

If you want faster results drop the HF step and use semi empirical methods instead. Or drop them all together:

Use pure DFT (BP86, M06-L) instead of the worst hybrid functional introduced. After that you can use a good hybrid like pbe0. 

Increase grid size for faster geometry convergence. Better wave function means better gradient, better gradient is a better step size, hence faster convergence. If you are analyzing the orbitals, the wave function should be well converged anyway.

You would want to check if you have a minimum, the best analysis will be worth nothing at a transition state, so run frequencies, but they take time.

There is no strict answer to this question - it depends on many factors, but generally it is recommended to start from default value, perform the optimization, and then enhance the perematers for reoptimization.

in reply to Martin's advice, which was generally good, but two words of caution:

i) certainly cc-pVDZ is not superior to 6-31+G(d,p) - there is a ton of literature on this fact (although cc-pVTZ is regarded as better than 6-311+G(d,p)). def2-SVP may not be a bad choice but is not always implemented in modern software as a default option, whereas pople basis sets are always available.

ii) I don't think B3LYP is 'the worst' functional - far from it, and it really depends on what application you are using it for. I agree it is overused and often in applications where it performs horribly. However, it still reproduces vibrational frequencies far better than many modern functionals (admittedly, by cancellation of errors) . It's also 'not bad' at geometry optimization either... However, in terms of reproducing energy levels yes, try a few different functionals - these are pretty good suggestions. Bear in mind that many modern functionals (e.g., M06, etc.) do perform a lot better when you to use larger grids, and thus they do take a bit longer to run.