OpenChemistry / avogadrolibs

Avogadro libraries provide 3D rendering, visualization, analysis and data processing useful in computational chemistry, molecular modeling, bioinformatics, materials science, and related areas.
https://two.avogadro.cc/
BSD 3-Clause "New" or "Revised" License
444 stars 171 forks source link

molecular orbitals are incorrect for a Cartesian Pople basis set #602

Closed ajs99778 closed 3 years ago

ajs99778 commented 3 years ago

Avogadro version: (please complete the following information from the About box):

Desktop version: (please complete the following information):

Describe the bug Molecular orbitals are incorrect. I compared Avogadro2 to cubegen, and several orbitals are noticeably different.

To Reproduce I've attached an fchk file. MO 30 is one of the more noticeable ones.

Expected behavior I think the issue is with 10F shells. Looking at some code (I assume Avogadro and Avogadro2 use similar code to process orbitals, but don't want to look around the code too much), Avogadro seems to use this order for the 10F shell:

  1. x^3
  2. x^2 * y
  3. x^2 * z
  4. x * y^2
  5. x y z
  6. x * z^2
  7. y^3
  8. y^2 * z
  9. y * z^2
  10. z^3

Unless these are getting reordered somewhere, this could be where things are going wrong. The order in the fchk file is:

  1. x^3
  2. y^3
  3. z^3
  4. x * y^2
  5. x^2 * y
  6. x^2 * z
  7. x * z^2
  8. y * z^2
  9. y^2 * z
  10. x y z

Gaussian uses 7F by default, so this probably wouldn't impact many people.

Screenshots There are two images in the attached folder. One is an image from Avogadro2 for MO 30. The other is a cube file for MO 30 generated by cubegen opened in ChimeraX.

Additional context

welcome[bot] commented 3 years ago

Thanks for opening your first issue here! Please try to include example files and screenshots if possible. If you're looking for support, please post on our forum: https://discuss.avogadro.cc/

welcome[bot] commented 3 years ago

Thanks for opening your first issue here! Please try to include example files and screenshots if possible. If you're looking for support, please post on our forum: https://discuss.avogadro.cc/

ghutchis commented 3 years ago

Yes, there are two orders, and different software use both different ordering and different normalization practices. The gau2grid documentation has a good description: https://gau2grid.readthedocs.io/en/latest/order.html

ghutchis commented 3 years ago

(Incidentally, Avogadro will read cube files just fine.)

ghutchis commented 3 years ago

I'm about to release 1.94, so this won't make that version. However, I intend on a more complete revisit to the MO code in a few weeks (including support for higher angular momentum orbitals, e.g., g and h).

ajs99778 commented 3 years ago

There might also be something going on with ORCA files in Avogadro. Here's MO 28 for benzene in Avogadro vs. a cube file made with orca_plot: orbits-tz-28 orbits-tz-28-orca Every other orbital I've looked at for this appears the same as orca_plot, but I only checked a few of them. I've also had issues opening other ORCA output files with MO info. I've attached orbits-tz.out (and corresponding gbw file for orca_plot), which has the odd MO 28, and sto-3g.out, which causes Avogadro v1.2.0 to crash when I try to open it. orca_files.zip

ghutchis commented 3 years ago

But the top image is from Avogadro 1.x, correct? That looks like a problem finding the isosurface, not in the cube generation. It's having trouble with some of the nodes at that resolution -- which is why the other orbitals look fine.

IMHO, the best thing for ORCA support is to add MKL support (#327). It shouldn't be too hard to duplicate the Molden reader...

ajs99778 commented 3 years ago

Yes, the top image is from Avogadro 1.2.0. Bottom is from ChimeraX again. I tried each of the resolution options available in Avogadro, and it always looked basically the same. I also opened the cube file from orca_plot in Avogadro and Avogadro 2, and it looks the same as the bottom image. orca_plot defaults to a 40x40x40 point cube, which works out to around 0.29x0.30x0.27 Angstrom resolution for this structure. I've also tried 100x100x100 points, which works out to ~0.1 Angstrom resolution. It looks like the bottom image in Avogadro and Avogadro 2.

I have found a few other orbitals in orbits-tz.out that differ from expected. The issue might be with parsing the orca file, which could probably be solved with an MKL parser as you suggest. The blocks with MO's that appear funky in Avogadro have a coefficient that looks like "-##.#####", which leaves no white space between the columns. Here's a section for MO's 43-48 (orca 0-indexes the orbitals):

                     42        43        44        45        46        47   
                   0.32620   0.32956   0.35528   0.35707   0.37766   0.40425
                   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
                  --------  --------  --------  --------  --------  --------
  0C   1s        -0.000710  0.000151  0.006035  0.020287 -0.035811  0.000214
  0C   2s        -0.001695  0.000333  0.015166  0.051217 -0.079271  0.000472
...
  1C   4s         0.027384 -0.003733  1.352275 -1.397765 -0.568898  0.032879
  1C   5s        -0.948696  0.322365 -5.545041  5.614592-33.159026 -0.155188
...

orbital 47 (46 by ORCA's numbering) also looks different in Avogadro than the cube file from orca_plot. If this is an issue with finding the isosurface, it seems to happen pretty frequently alongside MO coefficients less than -10. The only cases I've seen where an MO with a coefficient < -10 looks the same as the cube file orca_plot generates is when that orbital is in the first column (e.g. MO 120 by orca's numbering).

ghutchis commented 3 years ago

Reading the orbitals from the ORCA output goes through Open Babel, so we'd need to fix up that parsing code, plus get the orbitals output to CJSON for Avogadro2. That's a good idea, but would take more work than adding an MKL parser. (Particularly since IIRC MKL is v. similar to Molden.)

Let's make sure Avo2 is working on the correct data. Then we can question whether the marching cubes is having an issue on the isosurface.

I don't have time today to mess with MKL parsing. Basically, you'd copy https://github.com/OpenChemistry/avogadrolibs/blob/master/avogadro/quantumio/molden.cpp

ghutchis commented 3 years ago

I've also updated the avogadro-cclib script, which reads orbital / basis sets from ORCA (and others).

I'd still like to get support for MKL and Orca 5.0 new json format though.

github-actions[bot] commented 3 years ago

Here are the build results macOS.dmg Ubuntu-2004.tar.gz Win64.exe Artifacts will only be retained for 90 days.

ajs99778 commented 3 years ago

Cartesian orbitals from FCHK files look good to me!

I saw that they added some json format stuff to ORCA 5.0, but I can't find any info in the manual about getting ORCA to write a json file.

ghutchis commented 3 years ago

It's a new program orca2_json instead of orca2_mkl to process the .gbw file. I haven't tested it enough, but the .mkl files are in progress now.