Add read_gtensor and read_hyperfine_coupling to txtSileORCA

[tfrederiksen]

This PR enables easy reading of the electronic g-tensor and the hyperfine tensors from ORCA.

• Added tests for new/changed functions?
• Ran isort . and black . [24.2.0] at top-level
• Documentation for functionality in docs/
• Changes documented in CHANGELOG.md

[tfrederiksen]

I'm unsure about how we should handle units. ORCA reports the hfi-coupling in units of MHz. We could convert with a factor eV2MHz = 2.417989242e8 or, perhaps better, utilize the unit functionality of sisl (where MHz is currently missing).

[codecov]

Codecov Report

Attention: Patch coverage is 99.10714% with 1 lines in your changes are missing coverage. Please review.

Project coverage is 86.85%. Comparing base (a04a6ff) to head (e989f31).
Report is 5 commits behind head on main.

❗ Current head e989f31 differs from pull request most recent head b7ea782. Consider uploading reports for the commit b7ea782 to get more accurate results

Files	Patch %	Lines
src/sisl/io/orca/txt.py	98.61%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #722      +/-   ##
==========================================
+ Coverage   86.82%   86.85%   +0.03%     
==========================================
  Files         399      399              
  Lines       50835    50936     +101     
==========================================
+ Hits        44137    44240     +103     
+ Misses       6698     6696       -2     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[zerothi]

I just had a glance at the molecule3_property.txt output, some comments and suggestions:

• the GTensor seems to be something that can be calculated for several property indices (geom/prop):

  $ EPRNMR_GTensor
  description: Electronic g tensor
  geom. index: 1
  prop. index: 1
  Source density: 1 SCF
  Spin multiplicity: 2

  Would it make sense to extend this to allow reading more? Can ORCA do this?

  Great that this is an SCF quantity extracted via SileBinder :)
  How does this look in an MD + SCF? Our main usage with these indices is to refer to MD steps, and then scf steps are handled in the parsing method, so what does ORCA do? SCF-step output, or?

• the same thing can be said about the ATensor? Also, the ATensor is an atomic quantity, so basically it is indexed by the atom.
  I think the return here should be a complete array containing the data for all atoms.
  Perhaps using xarray would be useful here? It can be set-up in a tabulated fashion, which seems useful in an xarray context (also makes post-processing data much simpler!)

Thanks for this @tfrederiksen !

[zerothi]

I'm unsure about how we should handle units. ORCA reports the hfi-coupling in units of MHz. We could convert with a factor eV2MHz = 2.417989242e8 or, perhaps better, utilize the unit functionality of sisl (where MHz is currently missing).

Definitely we should add these things to the base units in sisl!

[zerothi]

And, should hfi be named hyperfine_atensor? Would this be more used in the litt.? The name hfi seems very short (read non-descriptive)

[tfrederiksen]

And, should hfi be named hyperfine_atensor? Would this be more used in the litt.? The name hfi seems very short (read non-descriptive)

I think people working with hyperfine interactions are used to the abbreviation "HFI", but I'm also fine with choosing a more descriptive name, perhaps read_hyperfine_coupling.

[zerothi]

I think people working with hyperfine interactions are used to the abbreviation "HFI", but I'm also fine with choosing a more descriptive name, perhaps read_hyperfine_coupling.

I would be in favor of that 👍

[tfrederiksen]

I just had a glance at the molecule3_property.txt output, some comments and suggestions:

the GTensor seems to be something that can be calculated for several property indices (geom/prop):
> $ EPRNMR_GTensor
> description: Electronic g tensor
> geom. index: 1
> prop. index: 1
> Source density: 1 SCF
> Spin multiplicity: 2

Would it make sense to extend this to allow reading more? Can ORCA do this?

I don't quite see the relevance of that. Within the EPRNMR module there is (as far as I understand) just a single g-tensor corresponding to the given (global) electron spin state computed. In case of molecule3 this is a spin doublet.

Great that this is an SCF quantity extracted via SileBinder :)

Do you mean that the use of SileBinder for stepping through the A-tensors is not the intended usage (only SCF stuff)?

How does this look in an MD + SCF? Our main usage with these indices is to refer to MD steps, and then scf steps are handled in the parsing method, so what does ORCA do? SCF-step output, or?

I am not sufficiently familiar with ORCA to answer this. I think in practice, the way we have used it, is to first determine a molecular geometry, then to switch to a specialized basis and to compute the g-tensor and HFIs.

the same thing can be said about the ATensor? Also, the ATensor is an atomic quantity, so basically it is indexed by the atom.

In the case of molecule3 indeed we've asked to compute the hyperfine couplings for all atoms through this specification in the input:

%EPRNMR
GTENSOR TRUE
NUCLEI = ALL C {AISO, ADIP}
NUCLEI = ALL H {AISO, ADIP}
END

However, the user may choose to compute for just a subset of nuclei. According to the ORCA manual:

% For example:
% calculates the hyperfine coupling for all nitrogen atoms
Nuclei = all N { aiso, adip, fgrad, rho};
% calculates the spin spin coupling constants for all carbon atoms
% assuming all carbon atoms are 13C
Nuclei = all C { ssall, ist = 13};
% etc

It is therefore not clear which nuclei to be found in the output.

I think the return here should be a complete array containing the data for all atoms.
Perhaps using xarray would be useful here? It can be set-up in a tabulated fashion, which seems useful in an xarray context (also makes post-processing data much simpler!)

I tend to think xarray is unnecessary here because at the end of the day, the tensors are just simple 3x3 matrices.

[zerothi]

I don't quite see the relevance of that. Within the EPRNMR module there is (as far as I understand) just a single g-tensor corresponding to the given (global) electron spin state computed. In case of molecule3 this is a spin doublet.

It was more if ORCA enabled one to split stuff into molecules, and then this would return for geom.index 2 or something similar.
But lets keep it. Thanks.

Great that this is an SCF quantity extracted via SileBinder :)

Do you mean that the use of SileBinder for stepping through the A-tensors is not the intended usage (only SCF stuff)?
No, just a comment, but, see below.

How does this look in an MD + SCF? Our main usage with these indices is to refer to MD steps, and then scf steps are handled in the parsing method, so what does ORCA do? SCF-step output, or?

I am not sufficiently familiar with ORCA to answer this. I think in practice, the way we have used it, is to first determine a molecular geometry, then to switch to a specialized basis and to compute the g-tensor and HFIs.

Ok, the point of the silebinder would be to do MD values. So read_gtensor[1] would yield the gtensor at MD step 2 (starting at 1).
So in this case if gtensor will always only be written once, then it shouldn't have the SileBinder() location.

the same thing can be said about the ATensor? Also, the ATensor is an atomic quantity, so basically it is indexed by the atom.

In the case of molecule3 indeed we've asked to compute the hyperfine couplings for all atoms through this specification in the input:

%EPRNMR
GTENSOR TRUE
NUCLEI = ALL C {AISO, ADIP}
NUCLEI = ALL H {AISO, ADIP}
END

However, the user may choose to compute for just a subset of nuclei. According to the ORCA manual:

% For example:
% calculates the hyperfine coupling for all nitrogen atoms
Nuclei = all N { aiso, adip, fgrad, rho};
% calculates the spin spin coupling constants for all carbon atoms
% assuming all carbon atoms are 13C
Nuclei = all C { ssall, ist = 13};
% etc

It is therefore not clear which nuclei to be found in the output.

I think the return here should be a complete array containing the data for all atoms.
Perhaps using xarray would be useful here? It can be set-up in a tabulated fashion, which seems useful in an xarray context (also makes post-processing data much simpler!)

I tend to think xarray is unnecessary here because at the end of the day, the tensors are just simple 3x3 matrices.

Ok, so disregard my all atoms return value. It should then return for all available atoms.

So here, a few comments. AFAIK, the current implementation would do: read_hyperfine_interaction[1] would yield the 2nd entry for whatever atom the user have requested.
This goes against the main utilization of the read_*[:] method, which should be MD indexable things.

So I would suspect the read_hyperfine_interaction to return a list of elements for all atoms requested.
So line 205 should read: self.step_to("EPRNMR_ATensor", ...) and then loop over Number of stored nuclei.

I would still advocate on moving this to xarray. For instance to easily search for the atom with the lowest eigenvalue. Or some other query, but we could easily start with a list of dicts with values if you wish.

A question, is the pre-factor something that should be factored on the Raw A tensor? If so shouldn't this just be done and then let it be discarded?

[tfrederiksen]

I am not sufficiently familiar with ORCA to answer this. I think in practice, the way we have used it, is to first determine a molecular geometry, then to switch to a specialized basis and to compute the g-tensor and HFIs.

Ok, the point of the silebinder would be to do MD values. So read_gtensor[1] would yield the gtensor at MD step 2 (starting at 1). So in this case if gtensor will always only be written once, then it shouldn't have the SileBinder() location.

OK I see. Will change this.

Ok, so disregard my all atoms return value. It should then return for all available atoms.

So here, a few comments. AFAIK, the current implementation would do: read_hyperfine_interaction[1] would yield the 2nd entry for whatever atom the user have requested. This goes against the main utilization of the read_*[:] method, which should be MD indexable things.

So I would suspect the read_hyperfine_interaction to return a list of elements for all atoms requested. So line 205 should read: self.step_to("EPRNMR_ATensor", ...) and then loop over Number of stored nuclei.

OK, will do.

I would still advocate on moving this to xarray. For instance to easily search for the atom with the lowest eigenvalue. Or some other query, but we could easily start with a list of dicts with values if you wish.

A list of dicts would be easier for me at this stage (how we currently have it in our scripts).

A question, is the pre-factor something that should be factored on the Raw A tensor? If so shouldn't this just be done and then let it be discarded?

To my understanding the raw tensor already contains "everything" (dipolar + Fermi contact contributions). I am not sure why the output also provides a prefactor value (varies per nuclear species) but it is a factor that the code has used to get the raw tensor. Maybe some users will want it, I don't know.

[zerothi]

Great! Thanks!

[tfrederiksen]

Maybe ready now?

[zerothi]

Many thanks @tfrederiksen !