remove unnecessary code

ichor_core/ichor/core/database/query_database.py

@@ -432,7 +432,7 @@ def write_processed_one_atom_data_to_csv(
  if global_forces_array is not None:
  b_matrix = form_b_matrix(atoms, alf, central_atom_index)
  negative_dE_df = convert_to_feature_forces(
- global_forces_array, b_matrix, alf, central_atom_index
+ global_forces_array, b_matrix
  )
  else:
  negative_dE_df = [None] * n_features
@@ -543,7 +543,7 @@ def write_processed_one_atom_data_to_csv(
  if global_forces_array is not None:
  b_matrix = form_b_matrix(atoms, alf, central_atom_index)
  negative_dE_df = convert_to_feature_forces(
- global_forces_array, b_matrix, alf, central_atom_index
+ global_forces_array, b_matrix
  )
  else:
  negative_dE_df = [None] * n_features

ichor_core/ichor/core/models/gaussian_energy_derivative_wrt_features.py

@@ -148,6 +148,11 @@ def b_matrix_true_finite_differences(atoms, system_alf, central_atom_idx=0, h=1e
  before calculating the features, defaults to 1e-6
 
  returns: A tuple of the analytical and finite differences B matrix
+
+ .. note::
+ The ordering of the columns (corresponding to 3N Cartesian coordinates) is always is the
+ original ordering corresponding to the Atoms instance that is passed in. This ensures that
+ the correct ordering is used when converting to feature forces.
  """
  from ichor.core.calculators import default_feature_calculator
 
@@ -214,9 +219,7 @@ def form_g_inverse(g_matrix: np.ndarray):
  return g_inv
 
 
-def convert_to_feature_forces(
- global_cartesian_forces: np.ndarray, b_matrix, system_alf, central_atom_idx
-):
+def convert_to_feature_forces(global_cartesian_forces: np.ndarray, b_matrix):
  """
  Compute -dE/df (since the global Cartesian forces are negative of the derivative of the potential).
  If making models with these values, need to take the NEGATIVE of -dE/df, as the
@@ -225,63 +228,39 @@ def convert_to_feature_forces(
  dE/df are the values that need to be used when adding derivatives to GP model.
 
  :param global_cartesian_forces: A 2D numpy array of shape (N_atoms, 3) containing the global Cartesian forces.
- The rows of this array are swapped internally to match the rows of the b-matrix.
  :param b_matrix: Wilson B matrix to be used to calculate, as well as dE/df. Should be of shape
  (3N-6) x 3N where N is the number of atoms.
  :param system_alf: A list of ALF instances containing the ALF of every atom.
  :param central_atom_idx: The index of the atom (0-indexed) which was the central ALF atom for which
  features (and feature derivatives) were calculated.
 
  .. note::
- The ordering of the forces should match up with the ordering of the b_matrix.
+ The ordering of the forces should match up with the ordering of the B matrix.
+ The 3N-6 rows of the B matrix correspond to the feature indices of the atom
+ The 3N columns of the B matrix correspond to the x,y,z coordinates of each atom,
+ where the ordering of the atoms is exactly the same as what is used to calculate
+ the B matrix.
 
  See Using Redundant Internal Coordinates to Optimize Equilibrium Geometries and Transition States
  https://doi.org/10.1002/(SICI)1096-987X(19960115)17:1<49::AID-JCC5>3.0.CO;2-0
  https://doi.org/10.1063/1.462844
  """
 
- # copy array so that original array is not altered unintentionally.
- copied_forces_array = global_cartesian_forces.copy()
- natoms = global_cartesian_forces.shape[0]
-
- # original row indices
- original_row_indices = [i for i in range(natoms)]
-
- # contains indices of central atom, x-axis atom, xy-plane atom
- alf_atom_indices = list(system_alf[central_atom_idx])
- if None in alf_atom_indices:
- alf_atom_indices.remove(None)
-
- # contains all other atom indices not in the alf
- non_local_atom_indices = [
- t for t in range(natoms) if t not in system_alf[central_atom_idx]
- ]
- # the correct order which the forces array should be in
- atom_indices_new_order = alf_atom_indices + non_local_atom_indices
-
- # swap rows of forces array
- copied_forces_array[
- [atom_indices_new_order, original_row_indices], :
- ] = copied_forces_array[[original_row_indices, atom_indices_new_order], :]
- copied_forces_array = copied_forces_array.flatten()
-
  g_matrix = form_g_matrix(b_matrix)
  # can use np.linalg.pinv here as well
  g_inv = form_g_inverse(g_matrix)
 
  # note that if the forces are passed in, will get the
  # feature forces, which are the -ve of the gradient
  # the gradient is what is used for models
- feature_forces = g_inv @ b_matrix @ copied_forces_array
+ feature_forces = g_inv @ b_matrix @ global_cartesian_forces
 
  return feature_forces
 
 
 def convert_to_cartesian_forces(
  dE_df_array: np.ndarray,
  b_matrix: np.ndarray,
- system_alf: List[ALF],
- central_atom_idx: int,
 ):
  """Converts from local 'feature' forces to global Cartesian forces, as given by Gaussian.
 
@@ -298,29 +277,7 @@ def convert_to_cartesian_forces(
  if dE_df_array.ndim == 0:
  dE_df_array = dE_df_array[..., np.newaxis]
 
- natoms = int(b_matrix.shape[-1] / 3)
-
- # original row indices
- original_row_indices = [i for i in range(natoms)]
-
- # contains indices of central atom, x-axis atom, xy-plane atom
- alf_atom_indices = list(system_alf[central_atom_idx])
- if None in alf_atom_indices:
- alf_atom_indices.remove(None)
-
- # contains all other atom indices not in the alf
- non_local_atom_indices = [
- t for t in range(natoms) if t not in system_alf[central_atom_idx]
- ]
- # the correct order which the forces array should be in
- atom_indices_new_order = alf_atom_indices + non_local_atom_indices
-
  # form a N_atoms x 3 Cartesian Forces array
  dE_dCart = np.matmul(b_matrix.T, dE_df_array).reshape(-1, 3)
 
- # swap rows of Cartesian array to match the original Atoms ordering
- dE_dCart[[original_row_indices, atom_indices_new_order], :] = dE_dCart[
- [atom_indices_new_order, original_row_indices], :
- ]
-
  return dE_dCart

ichor_core/ichor/core/models/predict_forces_for_all_atoms.py

@@ -13,7 +13,7 @@ def predict_fflux_forces_for_all_atoms(
  models: "ichor.core.models.Models", # noqa F821
  system_alf: List[ALF], # noqa F821
 ) -> np.ndarray:
- """Predicts the forces that FFLUX predicts (which are written to IQA_FORCES file).
+ """Predicts the Cartesian forces that FFLUX predicts (which are written to IQA_FORCES file).
 
  .. note::
  The atoms instance is converted to Bohr internally because the forces are per Bohr in FFLUX.