[BUG] Using ORCA output files

[RaphaelRobidas]

Describe the bug
I am trying to use chemtools with ORCA output files. Looking at the code, I've seen that both .wfn and .mkl files are parsed.

After using orca_2aim to generate the .wfn file, I try to load it using CondensedConceptualDFT.from_file(file_path, model="quadratic"), but get this error:

Traceback (most recent call last):
  File "cdft.py", line 11, in <module>
    tool = CondensedConceptualDFT.from_file(file_path, model="quadratic")
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/toolbox/conceptual.py", line 452, in from_file
    molecules = cls.load_file(fname)
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/toolbox/conceptual.py", line 142, in load_file
    molecule = Molecule.from_file(fnames)
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/wrappers/molecule.py", line 111, in from_file
    iodata = IOData.from_file(os.path.join(os.getcwd(), str(fname)))
  File "/home/raphael/1XX_General/109_Software/109.05_Python/mambaforce/envs/chemtools37/lib/python3.7/site-packages/horton/io/iodata.py", line 308, in from_file
    result.update(load_wfn(filename))
  File "/home/raphael/1XX_General/109_Software/109.05_Python/mambaforce/envs/chemtools37/lib/python3.7/site-packages/horton/io/wfn.py", line 206, in load_wfn
    mo_count, mo_occ, mo_energy, coefficients, energy = load_wfn_low(filename)
  File "/home/raphael/1XX_General/109_Software/109.05_Python/mambaforce/envs/chemtools37/lib/python3.7/site-packages/horton/io/wfn.py", line 94, in load_wfn_low
    num_mo, num_primitives, num_atoms = helper_num(f)
  File "/home/raphael/1XX_General/109_Software/109.05_Python/mambaforce/envs/chemtools37/lib/python3.7/site-packages/horton/io/wfn.py", line 47, in helper_num
    assert line.startswith('GAUSSIAN')
AssertionError

The line begins with GTO instead of GAUSSIAN. If I remove this assertion, I get this error instead:

Traceback (most recent call last):
  File "cdft.py", line 11, in <module>
    tool = CondensedConceptualDFT.from_file(file_path, model="quadratic")
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/toolbox/conceptual.py", line 452, in from_file
    molecules = cls.load_file(fname)
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/toolbox/conceptual.py", line 142, in load_file
    molecule = Molecule.from_file(fnames)
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/wrappers/molecule.py", line 111, in from_file
    iodata = IOData.from_file(os.path.join(os.getcwd(), str(fname)))
  File "/home/raphael/1XX_General/109_Software/109.05_Python/mambaforce/envs/chemtools37/lib/python3.7/site-packages/horton/io/iodata.py", line 308, in from_file
    result.update(load_wfn(filename))
  File "/home/raphael/1XX_General/109_Software/109.05_Python/mambaforce/envs/chemtools37/lib/python3.7/site-packages/horton/io/wfn.py", line 208, in load_wfn
    permutation = get_permutation_basis(type_assignment)
  File "/home/raphael/1XX_General/109_Software/109.05_Python/mambaforce/envs/chemtools37/lib/python3.7/site-packages/horton/io/wfn.py", line 163, in get_permutation_basis
    permutation = get_permutation_orbital(type_assignment)
  File "/home/raphael/1XX_General/109_Software/109.05_Python/mambaforce/envs/chemtools37/lib/python3.7/site-packages/horton/io/wfn.py", line 121, in get_permutation_orbital
    length = degeneracy[value]
KeyError: 4

If I use the .mkl file generated by orca_2mkl instead, I get this output instead:

User:           raphael
Machine info:   Linux Hydre 6.1.0-9-amd64 #1 SMP PREEMPT_DYNAMIC Debian 6.1.27-1 (2023-05-08) x86_64
Time:           2023-06-08T08:54:15.508964
Python version: 3.7.12 | packaged by conda-forge | (default, Oct 26 2021, 06:08:21) [GCC 9.4.0]
HORTON version: 2.1.0
Current Dir:    /home/raphael/tmp/orca
Command line:   cdft.py
HORTON module:  /home/raphael/1XX_General/109_Software/109.05_Python/mambaforce/envs/chemtools37/lib
                /python3.7/site-packages/horton/log.py

Initialized: <horton.gbasis.cext.GOBasis object at 0x14b38ee92910>
  Number of basis functions         : 138
  Number of normalization constants : 183
  Maximum shell type                : 3
  Center     0 :  S8 S8 S1 S1 S1 P3 P1 P1 P1 Dp1 Dp1 Dp1 Fp1 Fp1
  Center     1 :  S3 S1 S1 S1 P1 P1 P1 Dp1 Dp1
  Center     2 :  S3 S1 S1 S1 P1 P1 P1 Dp1 Dp1
  Center     3 :  S3 S1 S1 S1 P1 P1 P1 Dp1 Dp1
  Center     4 :  S3 S1 S1 S1 P1 P1 P1 Dp1 Dp1

Detected incorrect normalization of orbitals loaded from a file.
Trying to fix it as if it was a file generated by ORCA.
Initialized: <horton.gbasis.cext.GOBasis object at 0x14b38ee92280>
  Number of basis functions         : 138
  Number of normalization constants : 183
  Maximum shell type                : 3
  Center     0 :  S8 S8 S1 S1 S1 P3 P1 P1 P1 Dp1 Dp1 Dp1 Fp1 Fp1
  Center     1 :  S3 S1 S1 S1 P1 P1 P1 Dp1 Dp1
  Center     2 :  S3 S1 S1 S1 P1 P1 P1 Dp1 Dp1
  Center     3 :  S3 S1 S1 S1 P1 P1 P1 Dp1 Dp1
  Center     4 :  S3 S1 S1 S1 P1 P1 P1 Dp1 Dp1

Detected typical ORCA errors in file. Fixing them...
Traceback (most recent call last):
  File "cdft.py", line 11, in <module>
    tool = CondensedConceptualDFT.from_file(file_path, model="quadratic")
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/toolbox/conceptual.py", line 453, in from_file
    return cls.from_molecule(molecules, model, approach, scheme, **kwargs)
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/toolbox/conceptual.py", line 481, in from_molecule
    gcdft = GlobalConceptualDFT.from_molecule(molecule, model)
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/toolbox/conceptual.py", line 249, in from_molecule
    dict_energy = get_dict_energy(molecule)
  File "/home/raphael/1XX_General/109_Software/109.02_Chemistry/chemtools/chemtools/toolbox/utils.py", line 190, in get_dict_energy
    nelec + 1: molecule.energy + lumo_e,
TypeError: unsupported operand type(s) for +: 'NoneType' and 'float'

To Reproduce

from chemtools import CondensedConceptualDFT

file_path = "sp.mkl" # or "sp.wfn"
tool = CondensedConceptualDFT.from_file(file_path, model="quadratic")

.wfn and .mkl files

Expected behavior
I would expect chemtools to parse the files correctly like it does for .fchk files and calculate the relevant properties.

System Information:

• OS: Debian 11, Linux 6.1.0-9-amd64
• Python version: 3.7.12 (conda)
• Numpy version: 1.21.6
• Scipy version: 1.7.3

[m-rauen]

Hello Raphael! I've been having the same problem as you reported above, I think I found a solution to it but without explanation. It came to my attention from other friends' comments that orca_2aim can't always convert and/or provide a good wave function file. In that sense, I decided to try other alternatives and I stumble upon this post on Research Gate.

There are 2 valid alternatives in the comments section of the RG post mentioned before. 1) the approach mentioned by Tian Lu which uses Multiwfn software (a.k.a. his software); 2) the approach mentioned by Barend Greyling, which uses the Molden2AIM utility program.

Maybe this can be useful for your problem. Try it out and let me know!

[PaulWAyers]

I think the preferred method would be to make Orca generate a molden file, then generated a *.wfn with iodata. It may be that iodata handles molden files better than chemtools, because that is a more modern implementation.

In general, you should avoid using *.wfn files whenever you can. Because they are (very) restricted in terms of the basis sets they can handle, and force you to uncontracted contracted basis functions, post-processing of *.wfn files is unnecessarily inefficient.

@FarnazH or @leila-pujal would have more insight into this.