Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make functionalities insensitive to the capitalization of element names #21

Merged
merged 2 commits into from
Nov 3, 2024
Merged

Make functionalities insensitive to the capitalization of element names #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
36ea56e
Make fragment matching case-insensitive on element name
padix-key Nov 1, 2024
7f62318
Make parametrization of interactions case-insensitive
padix-key Nov 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/hydride/fragments.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -301,7 +301,8 @@ def _fragment(atoms, mask=None):
fragments = [None] * atoms.array_length()

all_bond_indices, all_bond_types = atoms.bonds.get_all_bonds()
elements = atoms.element
# Always convert to upper case to make the fragment matching case-insensitive
elements = np.char.upper(atoms.element)
charges = atoms.charge
coord = atoms.coord

Expand Down
99 changes: 50 additions & 49 deletions src/hydride/relax.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@ class MinimumFinder:
If this parameter is set, periodic boundary conditions are
taken into account (minimum-image convention), based on
the box vectors given with this parameter.

See also
--------
relax_hydrogen
Expand All @@ -195,13 +195,13 @@ class MinimumFinder:
float32 force_cutoff=10.0, box=None):
cdef int i, j, k, pair_i
cdef int32 atom_i, atom_j, bonded_atom_i

if atoms.array_length() != groups.shape[0]:
raise ValueError(
f"There are {atoms.array_length()} atoms, "
f"but {groups.shape[0]} group indicators"
)

self._prev_group_energies = None
self._prev_coord = atoms.coord.copy()
self._groups = np.asarray(groups)
Expand All @@ -228,7 +228,7 @@ class MinimumFinder:
radius=force_cutoff
)


cdef int32[:,:] bond_indices = atoms.bonds.get_all_bonds()[0]


Expand All @@ -245,7 +245,8 @@ class MinimumFinder:
(adj_indices.shape[0] * adj_indices.shape[1]),
dtype=np.int32
)
cdef uint8[:] mask = relevant_mask.astype(np.uint8)
# Always convert to upper case to make the parametrization case-insensitive
elements = np.char.upper(atoms.element)
pair_i = 0
for i in range(relevant_indices.shape[0]):
atom_i = relevant_indices[i]
Expand All @@ -265,8 +266,8 @@ class MinimumFinder:
continue
# Do not include interactions with atoms
# that have unknown, e.g. 'placeholder' elements
if atoms.element[atom_j] not in NB_VALUES.keys():
warnings.warn(f"Unknown element '{atoms.element[atom_j]}'")
if elements[atom_j] not in NB_VALUES.keys():
warnings.warn(f"Unknown element '{elements[atom_j]}'")
continue
interaction_pairs[pair_i, 0] = atom_i
interaction_pairs[pair_i, 1] = atom_j
Expand All @@ -277,7 +278,7 @@ class MinimumFinder:
self._interaction_groups = np.asarray(interaction_groups)[:pair_i]
# H-H-Interaction are included twice in the standard interaction
# pairs, which would give wrong global energy
# -> Deduplicate these pairs
# -> Deduplicate these pairs
self._dedup_interaction_mask = self._deduplicate_pairs()


Expand All @@ -291,7 +292,7 @@ class MinimumFinder:
cdef float32[:] elec_param = np.zeros(pair_i, dtype=np.float32)
for i in range(pair_i):
elec_param[i] = 332.0673 * (
charges[interaction_pairs[i, 0]] *
charges[interaction_pairs[i, 0]] *
charges[interaction_pairs[i, 1]]
)
self._elec_param = np.asarray(elec_param)
Expand All @@ -300,7 +301,7 @@ class MinimumFinder:
nb_values = np.array(
[
NB_VALUES.get(element, (np.nan, np.nan))
for element in atoms.element
for element in elements
],
dtype=np.float32
)
Expand All @@ -310,9 +311,9 @@ class MinimumFinder:
cdef float32[:] r_12 = np.zeros(pair_i, dtype=np.float32)
cdef float32[:] sq_eps = np.zeros(pair_i, dtype=np.float32)
# Special handlding for potential hydrogen bonds:
# If hydrogen in bound to a donor element the optimal distance
# If hydrogen in bound to a donor element the optimal distance
# to the possible acceptor is decreased
cdef uint8[:] hbond_mask = np.isin(atoms.element, HBOND_ELEMENTS) \
cdef uint8[:] hbond_mask = np.isin(elements, HBOND_ELEMENTS) \
.astype(np.uint8)
cdef float32 hbond_factor
# Calculate parameters for each H-X interaction
Expand All @@ -333,15 +334,15 @@ class MinimumFinder:
else:
hbond_factor = 1.0
r_6[i] = (hbond_factor * 0.5 * (
radii[atom_i] +
radii[atom_i] +
radii[atom_j]
))**6
r_12[i] = r_6[i]**2
sq_eps[i] = scales[atom_i] * scales[atom_j]
self._r_6 = np.asarray(r_6)
self._r_12 = np.asarray(r_12)
self._eps = np.sqrt(np.asarray(sq_eps))


def calculate_global_energy(self, coord):
"""
Expand Down Expand Up @@ -420,7 +421,7 @@ class MinimumFinder:
prev_coord[i, 0] = next_coord[i, 0]
prev_coord[i, 1] = next_coord[i, 1]
prev_coord[i, 2] = next_coord[i, 2]

# Prepare the reference energies for the next call of
# 'select_minimum()'
# Do this after this call of select_minimum(), instead of the beginning
Expand All @@ -431,10 +432,10 @@ class MinimumFinder:
)
self._prev_group_energies = self._sum_for_groups(prev_energies)
global_energy = np.sum(prev_energies[self._dedup_interaction_mask])

return self._prev_coord, global_energy, \
np.frombuffer(accept_next, dtype=bool)


def _calculate_energies(self, prev_coord, next_coord):
"""
Expand All @@ -457,15 +458,15 @@ class MinimumFinder:
for the corresponding atom group.
"""
cdef int i

distances = struc.distance(
prev_coord[self._interaction_pairs[:,1]],
next_coord[self._interaction_pairs[:,0]],
self._box
)
)
return self._pairwise_energy_function(distances)


def _sum_for_groups(self, float32[:] values):
"""
Sum up values (e.g. energies) for interaction pairs into
Expand All @@ -475,21 +476,21 @@ class MinimumFinder:
----------
values : ndarray, shape=(p,), dtype=np.float32
The input values.

Returns
-------
values : ndarray, shape=(g,), dtype=np.float32
The summed values.
"""
cdef int i

cdef float32[:] group_values = np.zeros(
self._n_groups, dtype=np.float32
)
cdef int32[:] groups = self._interaction_groups
for i in range(values.shape[0]):
group_values[groups[i]] += values[i]

return np.asarray(group_values)


Expand All @@ -502,7 +503,7 @@ class MinimumFinder:
----------
distances : ndarray, shape(p,), dtype=np.float32
The distances for each pair of interacting atoms.

Returns
-------
energy : ndarray, shape(p,), dtype=np.float32
Expand All @@ -516,7 +517,7 @@ class MinimumFinder:
-2 * self._r_6 / distances**6 + self._r_12 / distances**12
)
).astype(np.float32, copy=False)


def _deduplicate_pairs(self):
"""
Expand Down Expand Up @@ -601,7 +602,7 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
the box vectors given with this parameter.
If `box` is set to true, the box vectors are taken from the
``box`` attribute of `atoms` instead.

Returns
-------
relaxed_coord : ndarray, shape=(n,) or shape=(m,n), dtype=np.float32
Expand All @@ -617,11 +618,11 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
Notes
-----
The potential consists of the following terms:

.. math::

V = V_\text{el} + V_\text{nb}

V_\text{el} = 332.067
\sum_i^\text{H} \sum_j^\text{All}
\frac{q_i q_j}{D_{ij}}
Expand All @@ -631,7 +632,7 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
\left(
\frac{r_{ij}^{12}}{D_{ij}^{12}} - 2\frac{r_{ij}^6}{D_{ij}^6}
\right)

where :math:`D_{ij}` is the distance between the atoms :math:`i`
and :math:`j`. :math:`\epsilon_{ij}` and :math:`r_{ij}` are the
well depth and optimal distance between these atoms, respectively,
Expand All @@ -640,20 +641,20 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
.. math::

\epsilon_{ij} = \sqrt{ \epsilon_i \epsilon_j},

r_{ij} = \frac{r_i + r_j}{2}.

:math:`\epsilon_{i/j}` and :math:`r_{i/j}` are taken from the
*Universal Force Field* [1]_ [2]_.

References
----------

.. [1] AK Rappé, CJ Casewit, KS Colwell, WA Goddard III and WM Skiff,
"UFF, a full periodic table force field for molecular mechanics
and molecular dynamics simulations."
J Am Chem Soc, 114, 10024-10035 (1992).

.. [2] T Ogawa and T Nakano,
"The Extended Universal Force Field (XUFF): Theory and applications."
CBIJ, 10, 111-133 (2010)
Expand All @@ -678,7 +679,7 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
return return_coord, np.zeros(0)
else:
return return_coord

if box is not None:
if box is True:
if atoms.box is None:
Expand Down Expand Up @@ -710,7 +711,7 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
rotation_axes[i, 1] = coord[bonded_atom_index]
matrix_indices[hydrogen_indices] = i
rotation_freedom[i] = is_free

cdef float32[:,:,:] rot_mat_v = np.zeros(
(len(rotatable_bonds), 3, 3), dtype=np.float32
)
Expand Down Expand Up @@ -747,13 +748,13 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
cdef float32 sin_a, cos_a, icos_a
cdef float32 x, y, z
n = 0

# Loop terminates via break if result converges
# or iteration number is exceeded
while True:
if iterations is not None and n >= iterations:
break

# Generate next hydrogen conformation
# Calculate rotation matrices for these angles
for mat_i in range(angles_v.shape[0]):
Expand Down Expand Up @@ -802,7 +803,7 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
next_coord_v[i, 1] += support_v[mat_i, 1]
next_coord_v[i, 2] += support_v[mat_i, 2]


# Calculate next conformation based on energy
curr_coord, curr_energy, accepted_angles = minimum_finder.select_minimum(
next_coord.astype(np.float32, copy=False)
Expand All @@ -823,7 +824,7 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
break
if not np.isnan(prev_energy) and curr_energy > prev_energy:
# The relaxation algorithm allows the case, that the energy
# oscillates between multiple almost-minimum energies due to its
# oscillates between multiple almost-minimum energies due to its
# discrete nature and so convergence is never reached
# To prevent this, the relaxation terminates, if the energy
# of the accepted is higher than the one before
Expand All @@ -834,9 +835,9 @@ def relax_hydrogen(atoms, iterations=None, mask= None,
if return_trajectory:
trajectory.append(prev_coord.copy())
energies.append(curr_energy)

n += 1

if return_trajectory:
return_coord = np.stack(trajectory)
else:
Expand All @@ -863,7 +864,7 @@ def _find_rotatable_bonds(atoms, mask=None):
Ignore bonds, where the index of the heavy atom in the mask is
False.
By default no bonds are ignored.

Returns
-------
rotatable_bonds : list of tuple(int, int, bool, ndarray)
Expand All @@ -887,7 +888,7 @@ def _find_rotatable_bonds(atoms, mask=None):
f"but there are {atoms.array_length()} atoms"
)
atom_mask = mask.astype(np.uint8)

if atoms.bonds is None:
raise struc.BadStructureError(
"The input structure must have an associated BondList"
Expand All @@ -898,7 +899,7 @@ def _find_rotatable_bonds(atoms, mask=None):

cdef uint8[:] is_hydrogen = (atoms.element == "H").astype(np.uint8)
cdef uint8[:] is_nitrogen = (atoms.element == "N").astype(np.uint8)

cdef list rotatable_bonds = []

cdef int32[:] hydrogen_indices
Expand All @@ -910,15 +911,15 @@ def _find_rotatable_bonds(atoms, mask=None):
cdef bint is_rotatable
for i in range(all_bond_indices.shape[0]):
hydrogen_indices = np.zeros(4, np.int32)

if is_hydrogen[i] or not atom_mask[i]:
continue

bonded_heavy_index = -1
bonded_heavy_btype = -1
is_rotatable = True
h_i = 0

# Check for number of bonded heavy atoms
# and store bonded hydrogen atoms
for j in range(all_bond_indices.shape[1]):
Expand All @@ -939,7 +940,7 @@ def _find_rotatable_bonds(atoms, mask=None):
# -> no rotational freedom
is_rotatable = False
break

# There must be at least one bonded hydrogen atom
if h_i == 0:
is_rotatable = False
Expand Down Expand Up @@ -978,7 +979,7 @@ def _find_rotatable_bonds(atoms, mask=None):
is_free = False
else:
is_free = False

# 180 degrees rotation makes only sense if there is only one
# hydrogen atom, as two hydrogen atoms would simply replace each
# other after rotation
Expand Down
Loading
Loading