ghutchis / conformer-benchmark

Data and scripts for comprehensive benchmark of conformer relative energies
https://doi.org/10.1002/qua.26381
MIT License
23 stars 4 forks source link

Geometries are not fully optimized #1

Open diogomart opened 3 years ago

diogomart commented 3 years ago

Hello,

The maximum force on any atom of each conformer is, on average, 0.03139 Hartree/Bohr, which is about 70-fold larger than Gaussian's convergence criterion at 0.00045 Hartree/Bohr. This was calculated at the B3LYP/6-311G** level of theory for 5748 of the neutral conformers. Among these conformers, the lowest value of the maximum force on any atom of each conformer was 0.012 Hartree/Bohr.

An optimization of omegacsd_CNBPCT/rmsd92, a conformer depicted in Fig. 8 of the paper associated with this repository, using ORCA 4.2.1 and B3LYP-D3BJ def2-SVP (which is the level of theory used for optimization according to the paper), shows that the geometry of that conformer is not fully optimized. The energy decreases by about 7 kcal/mol after three geometry optimization steps, and more steps are needed for a full optimization. The calculated energy of this conformer as available in this repository is -4130.633327915281 Hartree, which is very close (~0.2 kcal/mol) from the energy tabulated in file energies/b3lyp-svp.txt.

ghutchis commented 3 years ago

The maximum force on any atom of each conformer is, on average, 0.03139 Hartree/Bohr, which is about 70-fold larger than Gaussian's convergence criterion at 0.00045 Hartree/Bohr. This was calculated at the B3LYP/6-311G** level of theory for 5748 of the neutral conformers.

If anything, the paper shows pretty clearly that there is a difference between B3LYP and B3LYP-D3BJ, as well as between basis sets. Moreover, when you switch programs, you switch grids, integral cutoffs and a host of other settings (RI, RICOS, etc.) -- especially with regards to geometry optimization criteria. Almost every QM program has its own criteria for declaring "done" in a geometry optimization, and that can also depend on the program version or optimization method (GDIIS, etc.).

So no, I would not expect geometries optimized w.r.t. Orca v3 B3LYP-D3BJ def2-SVP (Grid4 IIRC) to be optimized w.r.t Gaussian B3LYP/6-311G**

An optimization of omegacsd_CNBPCT/rmsd92, a conformer depicted in Fig. 8 of the paper associated with this repository, using ORCA 4.2.1 and B3LYP-D3BJ def2-SVP (which is the level of theory used for optimization according to the paper), shows that the geometry of that conformer is not fully optimized.

The geometries are from https://doi.org/10.1002/qua.25512 - which used Orca 3. While ~7 kcal/mol seems quite high, I also wouldn't be surprised if some of the geometry optimizations are not completely converged - there were over 7,000 optimizations. We did check for completion and re-submit as needed, but when performing batch calculations, there are undoubtedly issues.

In any case, the main point of the paper (and the repository) is to provide a set of single-point conformer geometries and energies on meaningful molecules.

In this particular set, we're certainly not going to re-optimize anything, since all corresponding single-point energies for that geometry would need to be re-computed.

We can certainly investigate tighter geometry optimization criteria for the next benchmark.

diogomart commented 3 years ago

I'll re-optimize with Gaussian and maybe run the DLPNO single point. I'll make that data available. This is a really nice dataset and the fully optimized geometries may be useful to others.

ghutchis commented 3 years ago

If you want to re-optimize, I'd highly suggest using B3LYP-D3BJ or -D4. There's minimal overhead and the geometries are significantly better. I'd also definitely go with def2-SVP or def2-TZVP family instead of the Pople basis sets.

diogomart commented 3 years ago

I'll use D3BJ. I never used the def2- basis sets but will definitely consider them. Thanks!