How to use the parmed.geometry.center_of_mass function in ParmEd

self,
        reference_pdb="reference.pdb",
        solvated_pdb="output.pdb",
        solvated_xml="system.xml",
        dummy_pdb="output.pdb",
        dummy_xml="output.xml",
    ):

        reference_structure = pmd.load_file(reference_pdb, structure=True)

        # Determine the offset coordinates for the new dummy atoms.
        if self.guest == "release":

            host_coordinates = reference_structure[f":{self.host_yaml['resname'].upper()}"].coordinates
            # Cheap way to get the center of geometry
            offset_coordinates = pmd.geometry.center_of_mass(host_coordinates,
                                                             masses=np.ones(len(host_coordinates)))

        else:
            guest_angle_restraint_mask = self.guest_yaml["restraints"]["guest"][-1]["restraint"][
                "atoms"
            ].split()

            offset_coordinates = reference_structure[f':{self.guest_yaml["name"].upper()} | :{self.host_yaml["resname"].upper()}']\
                [guest_angle_restraint_mask[1]].coordinates

        # First add dummy atoms to structure
        logger.debug(f"Adding dummy atoms to {solvated_pdb}")
        try:

            self._add_dummy_to_PDB(solvated_pdb, dummy_pdb, offset_coordinates,
                                   dummy_atom_tuples=[(0, 0, -6.0),

Parameters
    ----------
    structure : parmed.Structure
        Molecular structure containing coordinates
    mask : str
        Atom selection
    
    Returns
    -------
    np.array
        Coordinates of the selection center of mass
    """

    mask_coordinates = structure[mask].coordinates
    mask_masses = [atom.mass for atom in structure[mask].atoms]
    mask_com = pmd.geometry.center_of_mass(
        np.asarray(mask_coordinates), np.asarray(mask_masses)
    )
    return mask_com

Returns
    -------
    float
        The angle between the masks and the axis.
    """

    if "x" in axis.lower():
        axis = np.array([1.0, 0.0, 0.0])
    elif "y" in axis.lower():
        axis = np.array([0.0, 1.0, 0.0])
    elif "z" in axis.lower():
        axis = np.array([0.0, 0.0, 1.0])

    mask1_coordinates = structure[mask1].coordinates
    mask1_masses = [atom.mass for atom in structure[mask1].atoms]
    mask1_com = pmd.geometry.center_of_mass(
        np.asarray(mask1_coordinates), np.asarray(mask1_masses)
    )

    mask2_coordinates = structure[mask2].coordinates
    mask2_masses = [atom.mass for atom in structure[mask2].atoms]
    mask2_com = pmd.geometry.center_of_mass(
        np.asarray(mask2_coordinates), np.asarray(mask2_masses)
    )

    vector = mask2_com + -1.0 * mask1_com
    theta = np.arccos(
        np.dot(vector, axis) / (np.linalg.norm(vector) * np.linalg.norm(axis))
    )

    return theta