Skip to content

Commit

Permalink
orca2easyspin update (#229)
Browse files Browse the repository at this point in the history 
* orca2easyspin: restructure to allow version-dependent reading from main or property output files

* orca2easyspin: text-based prop file import S, g

* orca2easyspin: add support for A and Q

* orca2easyspin: small edits

* orac2easyspin: implement reading main output file

* orca2easyspin: improve file selection logic

* orca2easyspin: support multiple output structures

* orca2easyspin: calculate Q tensors from efg

* orca2easyspin: fix some tests

* orca2easyspin: some improvements

* orca2easyspin: impose axis convention on Q tensor

* orca2easyspin: centralize hf cutoff code, etc

* orca2easyspin: many improvements, tests

* orca2easyspin: sync main and prop data structures

* orca2easyspin: update documentation

* orca2easyspin etc: various small edits

* orca2easyspin: add direction test

* orca2easyspin: add direction test

* orca2easyspin: small edits, add error message for missing S

* orca2easyspin: update tests, better treatment of eigenvectors
  • Loading branch information
@stestoll
stestoll authored Feb 24, 2022
1 parent b8e92de commit e2e45f1
Show file tree
Hide file tree
Showing 70 changed files with 35,036 additions and 1,305 deletions.
63 changes: 38 additions & 25 deletions docsrc/orca.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,19 +26,21 @@

<section>

<h1>Import ORCA data</h1>
<h1>Import ORCA calculation results</h1>

<p>
<a href="http://orcaforum.kofo.mpg.de/" target="_blank">ORCA</a> is a quantum chemistry program that can calculate molecular EPR properties such as the g tensor, A tensors, electric field gradients, and D tensors.
<a href="http://orcaforum.kofo.mpg.de/" target="_blank">ORCA</a> is a quantum chemistry program that can calculate molecular EPR properties such as the g tensor, the zero-field splitting tensor, hyperfine coupling tensors, and electric field gradients for nuclear quadrupole coupling tensors.
</p>

<p>
Here, we show how you can import these ORCA-calculated tensors from ORCA output files into EasySpin, using the function <a class="esf" href="orca2easyspin.html">orca2easyspin</a>.
</p>

<!-- --------------------------------------------------------------------- -->
<div class="subtitle">Calculating EPR properties using ORCA</div>

<p>
In order to get EPR parameters from an ORCA calculation, you have to tell ORCA to calculate these parameters. Here is a very simple example of an ORCA input file for the hydroxyl radical. Let's name the file <code>hydroxyl.oif</code>:
In order to get EPR parameters from an ORCA calculation, you have to tell ORCA to calculate these parameters. Here is a simple example of an ORCA input file for the hydroxyl radical. Let's name the file <code>hydroxyl.oif</code>:

<pre class="nohighlight">
! UKS B3LYP 6-31G
Expand All @@ -54,47 +56,52 @@ <h1>Import ORCA data</h1>
</pre>

<p>
The <code>%eprnmr...end</code> block specifies which EPR parameters you want ORCA to calculate. <code>gtensor 1</code> instructs ORCA to calculate the g tensor. The <code>Nuclei</code> line tells ORCA to calculate for all hydrogens the following properties: Fermi contact hyperfine coupling (<code>aiso</code>), dipolar hyperfine coupling (<code>adip</code>), orbital contribution to hyperfine coupling (<code>aorb</code>), electric field gradient at nucleus (<code>fgrad</code>), and spin density at nucleus (<code>rho</code>).
The <code>%eprnmr...end</code> block specifies which EPR parameters you want ORCA to calculate. <code>gtensor 1</code> instructs ORCA to calculate the g tensor. The <code>Nuclei</code> line tells ORCA to calculate for all hydrogens the following properties: Fermi contact hyperfine coupling (<code>aiso</code>), dipolar hyperfine coupling (<code>adip</code>), orbital contribution to hyperfine coupling (<code>aorb</code>), the electric field gradient tensor at the nucleus (<code>fgrad</code>), and the spin density at the nucleus (<code>rho</code>).
</p>

<p>
For details about ORCA calculations, see the ORCA manual.
</p>

<p>
Next, you need to run the ORCA calculation. On the Windows/Linux/MacOS command/shell prompt, type
Next, run the ORCA calculation. On the Windows/Linux/MacOS command/shell prompt, type

<pre class="matlab">
orca hydroxyl.oif > hydroxyl.oof
</pre>

<p>
This will generate a long output text file <code>hydroxyl.oof</code>. Additionally, it generates a smaller binary file <code>hydroxyl.prop</code> that contains all the calculated properties.
This generate a long output text file <code>hydroxyl.oof</code>. Additionally, it generates a smaller file <code>hydroxyl_property.txt</code> that contains all the calculated properties. (Prior to version 5, ORCA generated a binary property file instead <code>hydroxyl.prop</code>.)
</p>

<!-- --------------------------------------------------------------------- -->
<div class="subtitle">Importing the results of the ORCA calculation</div>

<p>
You can now use EasySpin's function <a class="esf" href="orca2easyspin.html">orca2easyspin</a> to read the ORCA output and generate a spin system structure for you. EasySpin will not parse through the long text output file, but will directly read the binary propery file <code>hydroxyl.prop</code>. Here is how it works:
Next, use EasySpin's function <a class="esf" href="orca2easyspin.html">orca2easyspin</a> to read the ORCA output and generate a spin system structure for you. Here is how it works:

<pre class="matlab">
Sys = orca2easyspin('hydroxyl.prop')
Sys = orca2easyspin('hydroxyl.oof')
</pre>
<pre class="mloutput">
Spin not provided in the ORCA hydroxyl.prop file. Assuming S = 1/2.
Sys =
S: 0.5000
xyz: [2x3 double]
g: [2.0023 2.0065 2.0774]
gFrame: [0.1744 0 0]
Nucs: 'H'
A: [-9.7731 -116.1546 -140.4956]
AFrame: [2.9791e-12 1.5708 1.3965]
Q: [-0.0866 -0.1350 0.2216]
QFrame: [-2.9672 0 0]
Sys =
struct with fields:

S: 0.5000
g: [2.0021 2.0071 2.0768]
gFrame: [0.7254 0 0]
Nucs: 'H'
A: [21.9560 -102.4794 -140.8393]
AFrame: [0 1.5708 5.4378]
Q: [-0.0644 -0.0897 0.1541]
QFrame: [0.7254 0 0]
data: [1×1 struct]
NucsIdx: 1
</pre>

<p>
<code>Sys</code> is an EasySpin <a href="spinsystem.html">spin system structure</a>. All the fields are in the required units (MHz for tensors, radians for Euler angles). <code>Sys</code> is almost ready for use in EasySpin. The only additional information you need to supply is some <a href="broadenings.html">line broadening</a>, e.g. in <code>Sys.lwpp</code>.
</p>

<pre class="matlab">
Sys.lwpp = 1; % mT
Expand All @@ -107,21 +114,27 @@ <h1>Import ORCA data</h1>
<div class="subtitle">Hyperfine and quadrupole data for isotope mixtures</div>

<p>
If you look at the <code>Sys</code> structure above, you will see that the <code>Sys.Nucs</code> field contains an element (hydrogen), without specifying a specific isotope (1H or 2H). In this case, EasySpin simulates all the spectra of all isotopologues with significant natural abundance and combines the results.
In the imported <code>Sys</code> structure above, the <code>Sys.Nucs</code> field contains an element (hydrogen), without specifying a specific isotope (1H or 2H). In this case, EasySpin simulates all the spectra of all isotopologues with significant natural abundance and combines the results. You can adjust the isotopologue cutoff using <code>Opt.IsoCutoff</code>.
</p>

<p>
Of course, the hyperfine values are isotope-specific. For example, in the absence of any isotope effects, the hyperfine values of 2H are about 6.5 times smaller than those of 1H. The same holds for nuclear quadrupole coupling constants: They also depend on the isotope.
Of course, the hyperfine values are isotope-specific. For example, in the absence of any isotope effects, the hyperfine values of <sup>2</sup>H are about 6.5 times smaller than those of <sup>1</sup>>H. The same holds for nuclear quadrupole coupling constants: They also depend on the isotope.
</p>

<p>
EasySpin has simple rules to decide which isotope of a given element the values in <code>Sys.A</code> and in <code>Sys.Q</code> refer to:
EasySpin uses the following convention to decide isotope of a given element the values in <code>Sys.A</code> and in <code>Sys.Q</code> refer to:

<ol>
<li><code>Sys.A</code> refers to the most abundant among the isotopes with spin 1/2 or larger.
<li><code>Sys.Q</code> refers to the most abundant among the isotopes with spin 1 or larger.
<li><code>Sys.A</code> refers to the most abundant isotope among the isotopes with spin 1/2 or larger.
<li><code>Sys.Q</code> refers to the most abundant isotope among the isotopes with spin 1 or larger.
</ol>

<p>
All conversions from these leading isotopes to the others are automatically performed by EasySpin internally.
All conversions from these reference isotopes to the others are automatically performed by EasySpin internally.
</p>

<p>
For almost all elements that have an isotope with spin 1 or larger, the reference isotopes for A and Q happen to be the same. For example, for Cu, both are <sup>63</sup>Cu. The only exceptions are H (<sup>1</sup>H for A and <sup>2</sup>H for Q), Xe (<sup>129</sup>Xe for A and <sup>131</sup>Xe for Q), and Hg (<sup>199</sup>Hg for A and <sup>201</sup>Hg for Q). Several elements (He, C, F, Si, P, Fe, Se, etc.) do not have isotopes with spin 1 or larger.
</p>

<hr>
Expand Down
38 changes: 25 additions & 13 deletions docsrc/orca2easyspin.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,15 +29,15 @@
<div class="functitle">orca2easyspin</div>

<p>
Loads EPR properties calculated by ORCA.
Import EPR properties from ORCA output file.
</p>

<!-- ================================================================ -->
<div class="subtitle">Syntax</div>

<pre class="matlab">
Sys = orca2easyspin(PropFileName)
... = orca2easyspin(PropFileName,hfCutoff)
Sys = orca2easyspin(OrcaFileName)
... = orca2easyspin(OrcaFileName,hfCutoff)
</pre>

<p>
Expand All @@ -47,46 +47,58 @@
<div class="subtitle">Description</div>

<p>
This function loads EPR properties calculated by <a href="http://cec.mpg.de/forum" target="_blank">ORCA</a>. ORCA is a quantum chemistry program with extensive support for the calculation of EPR properties such as the g tensor, A tensors, electric field gradients and D tensors.
This function loads EPR properties calculated by <a href="http://cec.mpg.de/forum" target="_blank">ORCA</a>. ORCA is a quantum chemistry program with extensive support for the calculation of EPR properties such as the g tensor, the zero-field splitting tensor, hyperfine tensors, and electric field gradients for nuclear quadrupole tensors.

<p>
<code>PropFileName</code> is the filename of the binary property file generated by an ORCA calculation. This file has <code>*.prop</code> extension.
<code>OrcaFileName</code> is the filename, including extension, of the main output file generated by ORCA.
</p>

<p>
<code>Sys</code> is an EasySpin <a href="spinsystem.html">spin system structure</a> containing all the EPR parameters calculated by ORCA. Assuming all EPR properties of the spin system have been calculated by ORCA, only <a href="broadenings.html"></a>line broadenings</a> need to be added.
Alternative, the name of the property output file can be used. This property output file is text-based, has the extension <code>_property.txt</code>, and contains calculated properties in a condensed format compared to the main outut file. For older ORCA version prior to 5, the property output file is binary and has the extension <code>.prop</code>. There is no significant difference, or advantage, of using the property file compared to using the main output file.
</p>

<p>
<code>Sys</code> is an EasySpin <a href="spinsystem.html">spin system structure</a> containing all the spin Hamiltonian EPR parameters calculated by ORCA. Assuming all EPR properties of the spin system have been calculated by ORCA, only <a href="broadenings.html"></a>line broadenings</a> need to be added to be ready for an EPR spectral simulation.
</p>

<p>
In addition to the spin Hamiltonian parameters <code>Sys</code> also contains a list of atom coordinates in <code>Sys.xyz</code>, a list of element numbers in <code>Sys.Atoms</code>. The field <code>Sys.NucsIdx</code> list the atom index for every nucleus listed in <code>Sys.Nucs</code>. For example, if <code>Sys.Nucs='H'</code> and <code>Sys.NucsIdx=7</code>, this means that the hydrogen is atom number 7, and its coordinates are <code>Sys.xyz(7,:)</code>.
In addition to the spin Hamiltonian parameters <code>Sys</code> also contains additional information imported from ORCA, such as a list of atom coordinates, a list of element numbers.
</p>

<p>
<code>hfCutoff</code>, if given, specifies a cutoff level for hyperfine couplings. Any atoms with hyperfine couplings equal to or smaller than <code>hfCutoff</code> (in MHz) will be omitted from the spin system. This is useful for large molecules, where many of the ORCA-calculated hyperfine couplings would not be resolved in EPR.
The field <code>Sys.NucsIdx</code> list the atom index for every nucleus listed in <code>Sys.Nucs</code>. For example, if <code>Sys.Nucs='H'</code> and <code>Sys.NucsIdx=7</code>, this indicates that the hydrogen is atom number 7.
</p>

<p>
<code>hfCutoff</code>, if given, specifies a cutoff level for hyperfine couplings. Any atoms with hyperfine couplings equal to or smaller than <code>hfCutoff</code> (in MHz) will be omitted from the spin system. This is useful for large molecules, where many of the ORCA-calculated hyperfine couplings are small and would not be resolvable in an EPR spectrum.
</p>

<!-- ================================================================ -->
<div class="subtitle">Examples</div>

<p>
To load a spin system from an ORCA calculation with results in <code>mymol.oof</code> (ORCA output text file) and <code>mymol.prop</code> (associated binary ORCA property file), type
To load a spin system from an ORCA calculation with results in <code>molecule.out</code> (main ORCA output text file) or <code>molecule_property.txt</code> (associated ORCA property output file), use one of the following:
</p>

<pre class="matlab">
Sys = orca2easyspin('mymol'); % loads spin system from mymol.prop
Sys = orca2easyspin('molecule.out'); % loads spin system from main ORCA output file
Sys = orca2easyspin('molecule_property.txt'); % loads spin system from property output file
</pre>

<p>
To skip nuclei with hyperfine couplings less than 0.8 MHz, specify a second parameter:
</p>

<pre class="matlab">
Sys = orca2easyspin('mymol',0.8); % 0.8 MHz hyperfine cut-off
Sys = orca2easyspin('molecule.out',0.8); % exclude nuclei with <=0.8 MHz hyperfine cut-off
</pre>

<p>
In order to remove all nuclei except one or two of interest, use the function <a class="esf" href="nucspinkeep.html">nucspinkeep</a>:

<pre class="matlab">
Sys = orca2easyspin('mymol'); % load entire spin system from ORCA results
Sys = nucspinkeep(Sys,[1 2]); % remove all nuclei except the first and the second
Sys = orca2easyspin('molecule.out'); % load entire spin system from ORCA output
Sys = nucspinkeep(Sys,[1 2]); % remove all nuclei except the first and the second
</pre>

<!-- ====================================================== -->
Expand Down
1 change: 1 addition & 0 deletions docsrc/releases.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,7 @@ <h1>Changes from release to release</h1>
<li>The field <code>Sys.aF</code> in <a class="esf" href="spinsystem.html">spin systems</a> is no longer supported. Use <code>Sys.B4</code> and <code>Sys.B4Frame</code> instead.</li>
<li>Interconvert between spherical and cartesian tensors with the new functions <a class="esf" href="tensor_cart2sph.html">tensor_cart2sph</a> and <a class="esf" href="tensor_sph2cart.html">tensor_sph2cart</a>.</li>
<li><a class="esf" href="rotview.html">rotview</a> provides a graphical interface to visualize and explore Euler angles and their associated rotations.</li>
<li><a class="esf" href="orca2easyspin.html">orca2easyspin</a> supports ORCA5 and can do read in the main output file.
</ul>

<p>Major bug fixes</p>
Expand Down
Loading

0 comments on commit e2e45f1

Please sign in to comment.