Skip to content

Commit

Permalink
first stable release!
Browse files Browse the repository at this point in the history
  • Loading branch information
pablogila committed Jan 22, 2025
1 parent 98668a4 commit e79aae6
Show file tree
Hide file tree
Showing 41 changed files with 1,579 additions and 1,476 deletions.
2 changes: 1 addition & 1 deletion aton/_version.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,5 +10,5 @@
"""

__version__ = 'v0.0.1rc34'
__version__ = 'v0.0.1'

8 changes: 6 additions & 2 deletions aton/qrotor/classes.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ def __init__(
comment: str = None,
group: str = 'CH3',
E_levels: int = 5,
units = None, ################ TODO CHECK THAT THIS WORKS. previously = []
units = None, ################ TODO CHECK THAT THIS WORKS. previously = []
correct_potential_offset: bool = True,
save_eigenvectors: bool = False,
gridsize: int = None,
Expand Down Expand Up @@ -122,7 +122,11 @@ def summary(self):
}

def set_grid(self, gridsize:int=None):
"""Sets the `QSys.grid` to the specified `gridsize`."""
"""Sets the `QSys.grid` to the specified `gridsize`.
It removes the previous grid,
do **not** use this method if the potential was loaded externally!
"""
if gridsize is not None:
self.gridsize = gridsize
if self.gridsize is None:
Expand Down
10 changes: 5 additions & 5 deletions aton/qrotor/constants.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,17 +37,17 @@
"""Rotation radius of the amine group, in meters."""

# Inertia, SI units
I_CH = 3 * (phys.atoms['H'].mass * r_CH**2)
I_CH = 3 * (phys.atoms['H'].mass * phys.amu_to_kg * r_CH**2)
"""Inertia of CH3, in uma·m^2."""
I_CD = 3 * (phys.atoms['H'].isotope[2].mass * r_CH**2)
I_CD = 3 * (phys.atoms['H'].isotope[2].mass * phys.amu_to_kg * r_CH**2)
"""Inertia of CD3, in uma·m^2."""
I_NH = 3 * (phys.atoms['H'].mass * r_NH**2)
I_NH = 3 * (phys.atoms['H'].mass * phys.amu_to_kg * r_NH**2)
"""Inertia of NH3, in uma·m^2."""
I_ND = 3 * (phys.atoms['H'].isotope[2].mass * r_NH**2)
I_ND = 3 * (phys.atoms['H'].isotope[2].mass * phys.amu_to_kg * r_NH**2)
"""Inertia of ND3, in uma·m^2."""

# Rotational energy.
B_CH = ((phys.hbar_eV**2) / (2 * I_CH)) * phys.J_to_eV
B_CH = ((phys.hbar**2) / (2 * I_CH)) * phys.J_to_eV
"""Rotational energy of CH3, in eV·s/uma·m^2."""
B_CD = ((phys.hbar**2) / (2 * I_CD)) * phys.J_to_eV
"""Rotational energy of CD3, in eV·s/uma·m^2."""
Expand Down
49 changes: 30 additions & 19 deletions aton/qrotor/potential.py
Original file line number Diff line number Diff line change
Expand Up @@ -35,13 +35,13 @@
def load(
filepath:str='potential.dat',
system:QSys=None,
angle:str='deg',
energy:str='eV',
angle_unit:str='deg',
energy_unit:str='eV',
) -> QSys:
"""Read a potential rotational energy dataset.
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`.
The file in `filepath` should contain two columns with angle and potential energy values.
Degrees and eV are assumed as default units unless stated in `angle_unit` and `energy_unit`.
Units will be converted automatically to radians and eV.
"""
file_path = file.get(filepath)
Expand All @@ -57,21 +57,21 @@ def load(
positions.append(float(position.strip()))
potentials.append(float(potential.strip()))
# Save angles to numpy arrays
if angle.lower() in alias.units['deg']:
if angle_unit.lower() in alias.units['deg']:
positions = np.radians(positions)
elif angle.lower() in alias.units['rad']:
elif angle_unit.lower() in alias.units['rad']:
positions = np.array(positions)
else:
raise ValueError(f"Angle unit '{angle}' not recognized.")
raise ValueError(f"Angle unit '{angle_unit}' not recognized.")
# Save energies to numpy arrays
if energy.lower() in alias.units['meV']:
if energy_unit.lower() in alias.units['meV']:
potentials = np.array(potentials) * 1000
elif energy.lower() in alias.units['eV']:
elif energy_unit.lower() in alias.units['eV']:
potentials = np.array(potentials)
elif energy.lower() in alias.units['Ry']:
elif energy_unit.lower() in alias.units['Ry']:
potentials = np.array(potentials) * phys.Ry_to_eV
else:
raise ValueError(f"Energy unit '{energy}' not recognized.")
raise ValueError(f"Energy unit '{energy_unit}' not recognized.")
# Set the system
system.grid = np.array(positions)
system.gridsize = len(positions)
Expand All @@ -81,45 +81,56 @@ def load(

def from_qe(
folder=None,
filters:str=None,
output:str='potential.dat',
include:list=['.out'],
ignore:list=['slurm-'],
) -> None:
"""Creates a potential data file from Quantum ESPRESSO outputs.
The angle in degrees is extracted from the output filenames,
which must follow `whatever_ANGLE.out`.
Outputs from SCF calculations must be located in the provided `folder` (CWD if None),
and can be filtered with `filters`.
Outputs from SCF calculations must be located in the provided `folder` (CWD if None).
Files can be filtered by those containing the specified `filters`,
excluding those containing any string from the `ignore` list.
The `output` name is `potential.dat` by default.
"""
folder = file.get_dir(folder)
files = file.get_list(folder=folder, filters=filters, abspath=True)
files = file.get_list(folder=folder, include=include, ignore=ignore, abspath=True)
potential_data = '# Angle/deg Potential/eV\n'
potential_data_list = []
print('Extracting the potential as a function of the angle...')
print('----------------------------------')
counter_success = 0
counter_errors = 0
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 :(')
print(f'x {filename}')
counter_errors += 1
continue
energy = content['Energy'] * phys.Ry_to_eV
splits = filename.split('_')
angle = splits[-1].replace('.out', '')
angle = float(angle)
potential_data_list.append((angle, energy))
print(f'{filename} added')
print(f'OK {filename}')
counter_success += 1
# 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'
with open(output, 'w') as f:
f.write(potential_data)
print('----------------------------------')
print(f'Succesful calculations (OK): {counter_success}')
print(f'Faulty calculations (x): {counter_errors}')
print('----------------------------------')
print(f'Saved angles and potential values at {output}')
return None

Expand Down Expand Up @@ -153,7 +164,7 @@ def solve(system:QSys):
data = deepcopy(system)
# Is there a potential_name?
if not data.potential_name:
if not data.potential_values.any():
if not any(data.potential_values):
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)
Expand All @@ -162,7 +173,7 @@ def solve(system:QSys):
elif data.potential_name.lower() == 'sine':
data.potential_values = sine(data)
# At least there should be potential_values
elif not data.potential_values.any():
elif not any(data.potential_values):
raise ValueError("Unrecognised potential_name '{data.potential_name}' and no potential_values found")
return data.potential_values

Expand Down
78 changes: 39 additions & 39 deletions aton/qrotor/solve.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,11 @@
| | |
| --- | --- |
| `laplacian_matrix()` | Calculate the second derivative matrix for a given grid |
| `hamiltonian_matrix()` | Calculate the hamiltonian matrix of the system |
| `energies()` | Solve the system(s) for the `QSys` or `QExp` object |
| `potential()` | Solve the potential values of the system |
| `schrodinger()` | Solve the Schrödiger equation for the system |
| `energies()` | Solve the system(s) for the `QSys` or `QExp` object |
| `hamiltonian_matrix()` | Calculate the hamiltonian matrix of the system |
| `laplacian_matrix()` | Calculate the second derivative matrix for a given grid |
---
"""
Expand All @@ -29,27 +29,21 @@
import aton


def laplacian_matrix(grid):
"""Calculates the Laplacian (second derivative) matrix for a given `grid`."""
x = grid
diagonals = [-2*np.ones(len(x)), np.ones(len(x)), np.ones(len(x))]
laplacian_matrix = sparse.spdiags(diagonals, [0, -1, 1], format='lil')
# Periodic boundary conditions
laplacian_matrix[0, -1] = 1
laplacian_matrix[-1, 0] = 1
dx = x[1] - x[0]
laplacian_matrix /= dx**2
return laplacian_matrix


def hamiltonian_matrix(system:QSys):
"""Calculates the Hamiltonian matrix for a given `aton.qrotor.classes.QSys` object."""
V = system.potential_values.tolist()
potential = sparse.diags(V, format='lil')
B = system.B
x = system.grid
H = -B * laplacian_matrix(x) + potential
return H
def energies(var, filename:str=None) -> QExp:
"""Solves the Schrödinger equation for a given `aton.qrotor.classes.QSys` or `aton.qrotor.classes.QExp` object."""
if isinstance(var, QSys):
data = QExp()
data.systems = [deepcopy(var)]
elif isinstance(var, QExp):
data = deepcopy(var)
else:
raise TypeError('Input must be a QSys or QExp object.')
for system in data.systems:
system = potential(system)
system = schrodinger(system)
if filename:
aton.st.file.save(data, filename)
return data


def potential(system:QSys) -> QSys:
Expand Down Expand Up @@ -97,19 +91,25 @@ def schrodinger(system:QSys) -> QSys:
return system


def energies(var, filename:str=None) -> QExp:
"""Solves the Schrödinger equation for a given `aton.qrotor.classes.QSys` or `aton.qrotor.classes.QExp` object."""
if isinstance(var, QSys):
data = QExp()
data.systems = [deepcopy(var)]
elif isinstance(var, QExp):
data = deepcopy(var)
else:
raise TypeError('Input must be a QSys or QExp object.')
for system in data.systems:
system = potential(system)
system = schrodinger(system)
if filename:
aton.st.file.save(data, filename)
return data
def hamiltonian_matrix(system:QSys):
"""Calculates the Hamiltonian matrix for a given `aton.qrotor.classes.QSys` object."""
V = system.potential_values.tolist()
potential = sparse.diags(V, format='lil')
B = system.B
x = system.grid
H = -B * laplacian_matrix(x) + potential
return H


def laplacian_matrix(grid):
"""Calculates the Laplacian (second derivative) matrix for a given `grid`."""
x = grid
diagonals = [-2*np.ones(len(x)), np.ones(len(x)), np.ones(len(x))]
laplacian_matrix = sparse.spdiags(diagonals, [0, -1, 1], format='lil')
# Periodic boundary conditions
laplacian_matrix[0, -1] = 1
laplacian_matrix[-1, 0] = 1
dx = x[1] - x[0]
laplacian_matrix /= dx**2
return laplacian_matrix

61 changes: 40 additions & 21 deletions aton/st/file.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ def get(
if os.path.isfile(filepath):
return os.path.abspath(filepath)
elif os.path.isdir(filepath):
files = get_list(filepath, filters, abspath=True)
files = get_list(folder=filepath, include=filters, abspath=True)
elif return_anyway:
return None
else:
Expand All @@ -66,36 +66,55 @@ def get(


def get_list(
folder:str,
filters=None,
abspath:bool=True
folder:str=None,
include=None,
ignore=None,
abspath:bool=True,
also_folders:bool=False,
) -> list:
"""Return the files inside a `folder`, applying optional `filters`.
"""Return the files inside a `folder`, applying optional filters.
Only filenames containing all strings in the `include` list will be returned.
Filenames containing any string from the `ignore` list will be ignored.
The full paths are returned by default; to get only the base names, set `abspath = False`.
The CWD folder is used by default if no `folder` is provided.
It also returns folders if `also_folders = True`.
"""
if not folder:
folder = os.getcwd()
if os.path.isfile(folder):
folder = os.path.dirname(folder)
if not os.path.isdir(folder):
raise FileNotFoundError('Directory not found: ' + folder)
folder = os.path.abspath(folder)
files = os.listdir(folder)
# Apply filters or not
if filters is not None:
target_files = []
if not isinstance(filters, list):
filters = [str(filters)]
for filter_i in filters:
filter_i = os.path.basename(filter_i)
for f in files:
if filter_i in f:
target_files.append(f)
files = target_files
if not files:
return []
# Absolute paths?
if abspath:
filepaths = []
for f in files:
filepaths.append(os.path.join(folder, f))
files = filepaths
files = [os.path.join(folder, f) for f in files]
# Should we keep only files?
if not also_folders:
files = [f for f in files if not os.path.isdir(f if abspath else os.path.join(folder, f))]
# Apply filters if provided
if include is not None:
# Ensure include filters is always a list
if not isinstance(include, list):
include = [str(include)]
# Normalize filter names
include = [os.path.basename(i) for i in include]
# Only keep files that contain all filters
files = [f for f in files if all(filter_str in os.path.basename(f) for filter_str in include)]
# Remove files containing any string from the ignore list
if ignore is not None:
# Ensure ignore filters is always a list
if not isinstance(ignore, list):
ignore = [str(ignore)]
# Normalize ignoring filter names
ignore = [os.path.basename(i) for i in ignore]
# Exclude the corresponding files
files = [f for f in files if not any(filter_str in os.path.basename(f) for filter_str in ignore)]
return files


Expand Down Expand Up @@ -223,7 +242,7 @@ def copy_to_folders(
"""
if folder is None:
folder = os.getcwd()
old_files = get_list(folder, extension)
old_files = get_list(folder=folder, include=extension)
if old_files is None:
raise ValueError('No ' + extension + ' files found in path!')
for old_file in old_files:
Expand Down
2 changes: 1 addition & 1 deletion docs/aton.html
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ <h2>Submodules</h2>
</ul>


<footer>ATON v0.0.1rc34 documentation</footer>
<footer>ATON v0.0.1 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
Expand Down
4 changes: 2 additions & 2 deletions docs/aton/_version.html
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ <h2>API Documentation</h2>
</ul>


<footer>ATON v0.0.1rc34 documentation</footer>
<footer>ATON v0.0.1 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
Expand Down Expand Up @@ -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.1rc34&#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.1&#39;</span>
</span></pre></div>


Expand Down
Loading

0 comments on commit e79aae6

Please sign in to comment.