add atom names option for hexadecapole

ichor_core/ichor/core/multipoles/hexadecapole.py

@@ -1,5 +1,5 @@
 from collections import Counter
-from typing import Tuple
+from typing import List, Tuple, Union
 
 import numpy as np
 from ichor.core.common import constants
@@ -657,7 +657,7 @@ def displace_hexadecapole_cartesian(
 
 
 def recover_molecular_hexadecapole(
- atoms: "ichor.core.atoms.Atoms", # noqa F821
+ atoms: Union["ichor.core.Atoms", List["Atoms"]], # noqa F821
  ints_dir: "ichor.core.files.IntDirectory", # noqa F821
  convert_to_debye_angstrom_squared=True,
  convert_to_cartesian=True,
@@ -668,7 +668,7 @@ def recover_molecular_hexadecapole(
  from ichor.core.multipoles.octupole import octupole_spherical_to_cartesian
  from ichor.core.multipoles.quadrupole import quadrupole_spherical_to_cartesian
 
- atoms = atoms.to_bohr()
+ atoms = [a.to_bohr() for a in atoms]
 
  prefactor = 8 / 35
 
@@ -749,6 +749,7 @@ def recover_molecular_hexadecapole(
 def get_gaussian_and_aimall_molecular_hexadecapole(
  gaussian_output: "ichor.core.files.GaussianOutput", # noqa F821
  ints_directory: "ichor.core.files.IntsDir", # noqa: F821
+ atom_names: list = None,
 ):
  """Gets the Gaussian hexadecapole moment and converts it to traceless (still in Debye Angstrom^3)
  Also gets the AIMAll recovered molecule hexadecapole moment from atomic ones.
@@ -766,13 +767,25 @@ def get_gaussian_and_aimall_molecular_hexadecapole(
  must also be used in the AIMAll calculation.
  :param ints_directory: A IntsDirectory instance containing the
  AIMAll .int files for the same geometry that was used in Gaussian
+ :param atom_names: Optional list of atom names, which represent a subset of
+ the atoms. The atomic multipole moments for this subset of atoms will be summed
  :return: A tuple of 3x3x3x3 np.ndarrays, where the first is the Gaussian
  hexadecapole moment and the second is the AIMAll recovered hexadecapole moment.
  """
 
  # in angstroms, convert to bohr
  atoms = gaussian_output.atoms
  atoms = atoms.to_bohr()
+
+ if atom_names:
+ # ensure that the passed in atom names are a subset of the all of the atom names
+ if not set(atom_names).issubset(set(atoms.names)):
+ raise ValueError(
+ f"The passed atom names : {atom_names} must be a subset of all the atom names {atoms.names}."
+ )
+
+ atoms = [i for i in atoms if i.name in atom_names]
+
  # in debye angstrom squared
  raw_gaussian_octupole = np.array(gaussian_output.molecular_hexadecapole)
  # convert to xxx xxy xxz xyy xyz xzz yyy yyz yzz zzz