potential workflow from start to end

README.md

@@ -100,6 +100,23 @@ Values are accessed directly as `phys.value` or `phys.function()`.
 | [phys.functions](https://pablogila.github.io/ATON/aton/phys/functions.html) | Functions to sort and analyse element data, and to update the atoms dictionary |  
 
 
+## Quantum rotations
+
+The **QRotor** module is used to study the energy levels of quantum rotations,
+such as methyl and amine groups.
+
+### [aton.qrotor](https://pablogila.github.io/ATON/aton/qrotor.html)
+
+| | |
+| --- | --- |
+| [qrotor.classes](https://pablogila.github.io/ATON/aton/qrotor/classes.html)     | Definition of the `QSys` and `QExp` classes |
+| [qrotor.constants](https://pablogila.github.io/ATON/aton/qrotor/constants.html) | Bond lengths and inertias |
+| [qrotor.rotate](https://pablogila.github.io/ATON/aton/qrotor/rotate.html)       | Rotate specific atoms from structural files |
+| [qrotor.potential](https://pablogila.github.io/ATON/aton/qrotor/potential.html) | Potential definitions and loading functions |
+| [qrotor.solve](https://pablogila.github.io/ATON/aton/qrotor/solve.html)         | Solve rotation eigenvalues and eigenvectors |
+| [qrotor.plot](https://pablogila.github.io/ATON/aton/qrotor/plot.html)           | Plotting functions |
+
+
 ## Spectra analysis
 
 The **spx** module includes tools for spectral analysis from

aton/_version.py

@@ -10,5 +10,5 @@
 
 """
 
-__version__ = 'v0.0.1rc32'
+__version__ = 'v0.0.1rc34'
 

aton/qrotor/__init__.py

@@ -1,11 +1,7 @@
 """
-# ⚠️ UNDER DEVELOPMENT
-
-THIS PACKAGE IS STILL UNDER HEAVY DEVELOPMENT, DON'T USE IT
-
 # QRotor
  
-The QRotor subpackage is used to study the energy levels of quantum rotations, such as methyl and amine groups.
+The QRotor module is used to study the energy levels of quantum rotations, such as methyl and amine groups.
 
 # Index
 

aton/qrotor/classes.py

@@ -17,7 +17,6 @@
 
 
 import numpy as np
-from copy import deepcopy
 from .constants import *
 from aton.st import alias
 from aton.spx import Plotting
@@ -36,7 +35,8 @@ def __init__(
             save_eigenvectors: bool = False,
             gridsize: int = None,
             grid = None,
-            potential_name: str = None,
+            B: float = None,
+            potential_name: str = '',
             potential_constants: list = None,
             potential_values = None,
             ):
@@ -58,7 +58,11 @@ def __init__(
         self.gridsize: int = gridsize
         """Number of points in the grid."""
         self.grid = grid
-        """The grid with the points to be used in the calculation. Can be set automatically over $2 \\Pi$ with `QSys.set_grid()`."""
+        """The grid with the points to be used in the calculation.
+
+        Can be set automatically over $2 \\Pi$ with `QSys.set_grid()`.
+        Units should be in radians.
+        """
         if not B:
             B = self.B
         self.B: float = B
@@ -71,8 +75,10 @@ def __init__(
         """List of constants to be used in the calculation of the potential energy, in the `aton.qrotor.potential` module."""
         self.potential_values = potential_values
         """Numpy array with the potential values for each point in the grid.
+
         Can be calculated with a function available in the `qrotor.potential` module,
         or loaded externally with the `qrotor.potential.load()` function.
+        Units should be in eV.
         """
         self.potential_offset: float = None
         """`min(V)` before offset correction when `QSys.correct_potential_offset = True`"""

aton/qrotor/plot.py

@@ -6,19 +6,22 @@
 
 # Index
 
-`reduced_energies()`  
-`potential()`
-`energies_DEV()`      NOT IMPLEMENTED  
-`energy_DEV()`        NOT IMPLEMENTED  
-`convergence_DEV()`   NOT IMPLEMENTED  
-`eigenvectors_DEV()`  NOT IMPLEMENTED  
+| | |
+| --- | --- |
+| `reduced_energies()` | Reduced energies E/B as a function of the reduced potential V/B |
+| `potential()`        | Potential values as a function of the angle |
+| `energies_DEV()`     | NOT IMPLEMENTED |
+| `energy_DEV()`       | NOT IMPLEMENTED |
+| `convergence_DEV()`  | NOT IMPLEMENTED |
+| `eigenvectors_DEV()` | NOT IMPLEMENTED |
 
 ---
 """
 
 
 from .classes import *
 import matplotlib.pyplot as plt
+from copy import deepcopy
 
 
 def reduced_energies(data:QExp):
@@ -38,9 +41,23 @@ def reduced_energies(data:QExp):
     plt.show()
 
 
-def potential(data:QExp):                     #############   TODO IMPLEMENT
-    """Plot the potential values of the system"""
-    pass
+def potential(data):
+    """Plot the potential values of the system.
+
+    Input `data` can be a `QSys` or `QExp` object.
+    """
+    dataset = deepcopy(data)
+    if isinstance(dataset, QSys):
+        dataset = QExp(systems=[dataset], comment=data.comment)
+    if not isinstance(dataset, QExp):
+        raise TypeError(f'aton.qrotor.plot.potential requires a QExp or QSys object as input!')
+    for system in dataset.systems:
+        plt.plot(system.grid, system.potential_values, marker='', linestyle='-')
+    plt.xlabel('Angle / rad')
+    plt.ylabel('Potential energy / eV')
+    plt.xticks([0, np.pi/2, np.pi, 3*np.pi/2, 2*np.pi], ['0', r'$\frac{\pi}{2}$', r'$\pi$', r'$\frac{3\pi}{2}$', r'$2\pi$'])
+    plt.title(dataset.comment)
+    plt.show()
 
 
 def energies_DEV(data:QExp):
@@ -87,8 +104,7 @@ def energy_DEV(data:QExp):
     if not data.comment or (len(data.variables) == 1 and data.variables[0].comment):
         plt.title(f'{data.variables[0].comment}')
 
-    plt.xticks([0, np.pi/2, np.pi, 3*np.pi/2, 2*np.pi],
-                ['0', r'$\frac{\pi}{2}$', r'$\pi$', r'$\frac{3\pi}{2}$', r'$2\pi$'])
+    plt.xticks([0, np.pi/2, np.pi, 3*np.pi/2, 2*np.pi], ['0', r'$\frac{\pi}{2}$', r'$\pi$', r'$\frac{3\pi}{2}$', r'$2\pi$'])
 
     unique_potentials = []
     unique_elements = []

aton/qrotor/potential.py

@@ -24,6 +24,7 @@
 from . import constants
 import numpy as np
 import os
+from copy import deepcopy
 from scipy.interpolate import CubicSpline
 import aton.st.alias as alias
 import aton.st.file as file
@@ -32,15 +33,16 @@
 
 
 def load(
-        filepath,
+        filepath:str='potential.dat',
         system:QSys=None,
         angle:str='deg',
         energy:str='eV',
         ) -> QSys:
     """Read a potential rotational energy dataset.
 
-    The file should contain two columns with `angle` and `potential` values.
+    The file in `filepath` should contain two columns with `angle` and `potential` values.
     Degrees and eV are assumed as default units unless stated in `angle` and `energy`.
+    Units will be converted automatically to radians and eV.
     """
     file_path = file.get(filepath)
     system = QSys() if system is None else system
@@ -52,8 +54,8 @@ def load(
         if line.startswith('#'):
             continue
         position, potential = line.split()
-        positions.append(float(position))
-        potentials.append(float(potential))
+        positions.append(float(position.strip()))
+        potentials.append(float(potential.strip()))
     # Save angles to numpy arrays
     if angle.lower() in alias.units['deg']:
         positions = np.radians(positions)
@@ -64,8 +66,10 @@ def load(
     # Save energies to numpy arrays
     if energy.lower() in alias.units['meV']:
         potentials = np.array(potentials) * 1000
-    elif energy.lower() in alias.units['ev']:
+    elif energy.lower() in alias.units['eV']:
         potentials = np.array(potentials)
+    elif energy.lower() in alias.units['Ry']:
+        potentials = np.array(potentials) * phys.Ry_to_eV
     else:
         raise ValueError(f"Energy unit '{energy}' not recognized.")
     # Set the system
@@ -91,26 +95,29 @@ def from_qe(
     """
     folder = file.get_dir(folder)
     files = file.get_list(folder=folder, filters=filters, abspath=True)
-    potential_data = '# Angle / deg, Potential / eV\n'
+    potential_data = '# Angle/deg    Potential/eV\n'
     potential_data_list = []
+    print('Extracting the potential as a function of the angle...')
     for filepath in files:
+        filename = os.path.basename(filepath)
         filepath = file.get(filepath=filepath, filters='.out', return_anyway=True)
         if not filepath:  # Not an output file, skip it
             continue
         content = qe.read_out(filepath)
         if not content['Success']:  # Ignore unsuccessful calculations
+            print(f'{filename} was unsuccessful :(')
             continue
         energy = content['Energy'] * phys.Ry_to_eV
-        filename = os.path.basename(filepath)
         splits = filename.split('_')
         angle = splits[-1].replace('.out', '')
         angle = float(angle)
         potential_data_list.append((angle, energy))
+        print(f'{filename} added')
     # Sort by angle
     potential_data_list_sorted = sorted(potential_data_list, key=lambda x: x[0])
     # Append the sorted values as a string
     for angle, energy in potential_data_list_sorted:
-        potential_data += f'{angle}, {energy}\n'
+        potential_data += f'{angle}    {energy}\n'
     with open(output, 'w') as f:
         f.write(potential_data)
     print(f'Saved angles and potential values at {output}')
@@ -127,22 +134,37 @@ def interpolate(system:QSys) -> QSys:
     new_V = cubic_spline(new_grid)
     system.grid = new_grid
     system.potential_values = new_V
+    print(f"Interpolated potential to a grid of size {gridsize}")
     return system
 
 
 # Redirect to the desired potential energy function
 def solve(system:QSys):
-    """Solves `aton.qrotor.classes.QSys.potential_values`, according to the `aton.qrotor.classes.QSys.potential_name`."""
-    if system.potential_name.lower() == 'titov2023':
-        return titov2023(system)
-    elif system.potential_name.lower() == 'zero':
-        return zero(system)
-    elif system.potential_name.lower() == 'sine':
-        return sine(system)
-    elif system.potential_values:  # Re
-        return system.potential_values
-    else:
-        raise ValueError(f'Unrecognised potential_name ({system.potential_name}) and no potential_values found in system!')
+    """Solves `aton.qrotor.classes.QSys.potential_values`,
+    according to the `aton.qrotor.classes.QSys.potential_name`.
+    Returns the new `potential_values`.
+
+    If `QSys.potential_name` is not present or not recognised,
+    the current `QSys.potential_values` are used.
+
+    If a bigger `QSys.gridsize` is provided,
+    the potential is also interpolated to the new gridsize.
+    """
+    data = deepcopy(system)
+    # Is there a potential_name?
+    if not data.potential_name:
+        if not data.potential_values.any():
+            raise ValueError(f'No potential_name and no potential_values found in the system!')
+    elif data.potential_name.lower() == 'titov2023':
+        data.potential_values = titov2023(data)
+    elif data.potential_name.lower() == 'zero':
+        data.potential_values = zero(data)
+    elif data.potential_name.lower() == 'sine':
+        data.potential_values = sine(data)
+    # At least there should be potential_values
+    elif not data.potential_values.any():
+        raise ValueError("Unrecognised potential_name '{data.potential_name}' and no potential_values found")
+    return data.potential_values
 
 
 def zero(system:QSys):

aton/qrotor/solve.py

@@ -21,6 +21,7 @@
 
 from .classes import *
 from .potential import solve as solve_potential
+from .potential import interpolate
 from copy import deepcopy
 import time
 import numpy as np
@@ -58,6 +59,8 @@ def potential(system:QSys) -> QSys:
     Then it applies extra operations, such as removing the potential offset
     if `aton.qrotor.classes.QSys.correct_potential_offset = True`.
     """
+    if system.gridsize > len(system.grid):
+        system = interpolate(system)
     V = solve_potential(system)
     if system.correct_potential_offset is True:
         offset = min(V)
@@ -102,7 +105,7 @@ def energies(var, filename:str=None) -> QExp:
     elif isinstance(var, QExp):
         data = deepcopy(var)
     else:
-        raise ValueError('Input must be a QSys or QExp object.')
+        raise TypeError('Input must be a QSys or QExp object.')
     for system in data.systems:
         system = potential(system)
         system = schrodinger(system)

docs/aton.html

@@ -71,6 +71,7 @@ <h2>Contents</h2>
   <ul>
     <li><a href="#interfaces-for-ab-initio-codes">Interfaces for <em>ab-initio</em> codes</a></li>
     <li><a href="#physico-chemical-constants">Physico-chemical constants</a></li>
+    <li><a href="#quantum-rotations">Quantum rotations</a></li>
     <li><a href="#spectra-analysis">Spectra analysis</a></li>
     <li><a href="#general-text-edition">General text edition</a></li>
     <li><a href="#system-tools">System tools</a></li>
@@ -97,7 +98,7 @@ <h2>Submodules</h2>
             </ul>
 
 
-            <footer>ATON v0.0.1rc32 documentation</footer>
+            <footer>ATON v0.0.1rc34 documentation</footer>
 
         <a class="attribution" title="pdoc: Python API documentation generator" href="https://pdoc.dev" target="_blank">
             built with <span class="visually-hidden">pdoc</span><img
@@ -249,6 +250,48 @@ <h3 id="atonphys"><code><a href="aton/phys.html">aton.phys</a></code></h3>
 </tbody>
 </table>
 
+<h2 id="quantum-rotations">Quantum rotations</h2>
+
+<p>The <strong>QRotor</strong> module is used to study the energy levels of quantum rotations,
+such as methyl and amine groups.</p>
+
+<h3 id="atonqrotor"><code><a href="aton/qrotor.html">aton.qrotor</a></code></h3>
+
+<table>
+<thead>
+<tr>
+  <th></th>
+  <th></th>
+</tr>
+</thead>
+<tbody>
+<tr>
+  <td><code><a href="aton/qrotor/classes.html">aton.qrotor.classes</a></code></td>
+  <td>Definition of the <code>QSys</code> and <code>QExp</code> classes</td>
+</tr>
+<tr>
+  <td><code><a href="aton/qrotor/constants.html">aton.qrotor.constants</a></code></td>
+  <td>Bond lengths and inertias</td>
+</tr>
+<tr>
+  <td><code><a href="aton/qrotor/rotate.html">aton.qrotor.rotate</a></code></td>
+  <td>Rotate specific atoms from structural files</td>
+</tr>
+<tr>
+  <td><code><a href="aton/qrotor/potential.html">aton.qrotor.potential</a></code></td>
+  <td>Potential definitions and loading functions</td>
+</tr>
+<tr>
+  <td><code><a href="aton/qrotor/solve.html">aton.qrotor.solve</a></code></td>
+  <td>Solve rotation eigenvalues and eigenvectors</td>
+</tr>
+<tr>
+  <td><code><a href="aton/qrotor/plot.html">aton.qrotor.plot</a></code></td>
+  <td>Plotting functions</td>
+</tr>
+</tbody>
+</table>
+
 <h2 id="spectra-analysis">Spectra analysis</h2>
 
 <p>The <strong>spx</strong> module includes tools for spectral analysis from

docs/aton/_version.html

@@ -76,7 +76,7 @@ <h2>API Documentation</h2>
     </ul>
 
 
-            <footer>ATON v0.0.1rc32 documentation</footer>
+            <footer>ATON v0.0.1rc34 documentation</footer>
 
         <a class="attribution" title="pdoc: Python API documentation generator" href="https://pdoc.dev" target="_blank">
             built with <span class="visually-hidden">pdoc</span><img
@@ -116,7 +116,7 @@ <h1 class="modulename">
 </span><span id="L-10"><a href="#L-10"><span class="linenos">10</span></a>
 </span><span id="L-11"><a href="#L-11"><span class="linenos">11</span></a><span class="sd">&quot;&quot;&quot;</span>
 </span><span id="L-12"><a href="#L-12"><span class="linenos">12</span></a>
-</span><span id="L-13"><a href="#L-13"><span class="linenos">13</span></a><span class="n">__version__</span> <span class="o">=</span> <span class="s1">&#39;v0.0.1rc32&#39;</span>
+</span><span id="L-13"><a href="#L-13"><span class="linenos">13</span></a><span class="n">__version__</span> <span class="o">=</span> <span class="s1">&#39;v0.0.1rc34&#39;</span>
 </span></pre></div>