How to use the parmed.load_file function in ParmEd

import parmed as pmd
import numpy as np
import logging as log

from paprika import restraints
from paprika import analysis

logger = log.getLogger()
logger.setLevel(log.DEBUG)
log.basicConfig(format="%(asctime)s %(message)s", datefmt="%Y-%m-%d %I:%M:%S %p")

inputpdb = pmd.load_file("cb6_but_gb_apr_ref_data/vac.pdb")

# Distance restraint
rest1 = restraints.DAT_restraint()
rest1.continuous_apr = True
rest1.amber_index = True
rest1.topology = inputpdb
rest1.mask1 = ":CB6@O"
rest1.mask2 = ":BUT@C1"
rest1.attach["target"] = 4.5
rest1.attach["fraction_list"] = [0.00, 0.04, 0.181, 0.496, 1.000]
rest1.attach["fc_final"] = 5.0
rest1.pull["fc"] = rest1.attach["fc_final"]
rest1.pull["target_initial"] = rest1.attach["target"]
rest1.pull["target_final"] = 18.5
rest1.pull["num_windows"] = 19
rest1.initialize()

def runNCMC(platform_name, nstepsNC, nprop, outfname):

    logger = init_logger()

    #Generate the ParmEd Structure
    prmtop = utils.get_data_filename('blues', 'tests/data/vacDivaline.prmtop')
    inpcrd = utils.get_data_filename('blues', 'tests/data/vacDivaline.inpcrd')
    struct = parmed.load_file(prmtop, xyz=inpcrd)
    logger.info('Structure: %s' % struct.topology)

    #Define some options
    opt = { 'temperature' : 300.0, 'friction' : 1, 'dt' : 0.004,
            'hydrogenMass' : 3.024,
            'nIter' : 100, 'nstepsNC' : nstepsNC, 'nstepsMD' : 1000, 'nprop' : nprop,
            'nonbondedMethod' : 'PME', 'nonbondedCutoff': 10,
            'constraints': 'HBonds',
            'trajectory_interval' : 100, 'reporter_interval' : 250,
            'ncmc_traj' : None, 'write_move' : False,
            'platform' : platform_name,
            'outfname' : 'vacDivaline',
            'verbose' : False}

    for k,v in opt.items():
        logger.debug('Options: {} = {}'.format(k,v))

them to the end of the data file. If so, update data structures accordingly
and handle conversion appropriately.

    Notes
    -----
    Currently, this function ensures that - after AMBER conversion reorders
    water molecules with residue names 'WAT' to occur last in the resulting
    parameter/coordinate files - the internal data structures are updated to
    have the correct order in the relevant lists (labels, smiles_strings,
    n_monomers). If for some reason GROMACS conversion were removed, these
    would need to be updated elsewhere. (Probably this should be done anyway,
    as this is not really a GROMACS issue.)

            """
            # Read in AMBER format parameter/coordinate file and convert in gromacs
            gromacs_topology = parmed.load_file(self.prmtop_filename, self.inpcrd_filename )

            # Split the topology into components and check that we have the right number of components
            components = gromacs_topology.split()
            assert len(components)==len(self.n_monomers), "Number of monomers and number of components in the combined topology do not match."


            ####    HANDLE ORDERING OF WATER   ####
            # Check if any of the residues is named "WAT". If it is, antechamber will potentially have re-ordered it from where it was (it places residues named "WAT" at the end) so it may no longer appear in the order in which we expect.
            resnames = [ components[i][0].residues[0].name for i in range(len(components)) ]
            wat_present = False


            # Manage presence of WAT residues and possible re-ordering
            if 'WAT' in resnames:

                # If there is a water present, then we MIGHT have re-ordering.

if not self.guest == "release":
            # Align the host-guest complex so the first guest atom is at (0, 0, 0) and the second guest atom lies
            # along the positive z-axis.
            guest_angle_restraint_mask = self.guest_yaml["restraints"]["guest"][-1]["restraint"][
                "atoms"
            ].split()
            aligned_structure = align.zalign(
                structure, guest_angle_restraint_mask[1], guest_angle_restraint_mask[2]
            )
            aligned_structure.save(str(intermediate_pdb), overwrite=True)

        else:
            # Create a PDB file just for the host.

            host = pmd.load_file(str(input_pdb), structure=True)
            host_coordinates = host[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)))

            # Find the principal components, take the two largest, and find the vector orthogonal to that
            # (should be cross-product right hand rule, I think). Use that vector to align with the z-axis.
            # This may not generalize to non-radially-symmetric host molecules.

            aligned_coords = np.empty_like(structure.coordinates)
            for atom in range(len(structure.atoms)):
                aligned_coords[atom] = structure.coordinates[atom] - offset_coordinates
            structure.coordinates = aligned_coords

            inertia_tensor = np.dot(structure.coordinates.transpose(), structure.coordinates)
            eigenvalues, eigenvectors = np.linalg.eig(inertia_tensor)

def build_monomers(self):
        """Generate GAFF mol2 and frcmod files for each chemical."""
        for k, smiles_string in enumerate(self.smiles_strings):
            mol2_filename = self.gaff_mol2_filenames[k]
            frcmod_filename = self.frcmod_filenames[k]
            inpcrd_filename = self.inpcrd_filenames[k]
            prmtop_filename = self.prmtop_filenames[k]
            sdf_filename = self.sdf_filenames[k]
            if not (os.path.exists(mol2_filename) and os.path.exists(frcmod_filename)):
                #Convert SMILES strings to mol2 and frcmod files for antechamber
                openmoltools.openeye.smiles_to_antechamber(smiles_string, mol2_filename, frcmod_filename)
                #Correct the mol2 file partial atom charges to have a total net integer molecule charge  
                mol2f = parmed.formats.Mol2File
                mol2f.write(parmed.load_file(mol2_filename).fix_charges(),mol2_filename)
            #Generate amber coordinate and topology files for the unsolvated molecules
            mol_name = os.path.basename(mol2_filename).split('.')[0]
            openmoltools.amber.run_tleap(mol_name, mol2_filename,frcmod_filename, prmtop_filename, inpcrd_filename)
    
            #Read Mol2 File and write SDF file
            mol2tosdf.writeSDF(mol2_filename, sdf_filename, mol_name)

        #Generate unique residue names for molecules in mol2 files
        openmoltools.utils.randomize_mol2_residue_names( self.gaff_mol2_filenames )

top = struct.topology
    pos = np.array([list(i) for i in struct.positions.value_in_unit(unit.nanometers)])*unit.nanometers

    lig_atoms = list(range(2635,2649))
    struct.positions = pos
    sys = struct.createSystem(nonbondedMethod=openmm.app.PME)

    thermo = ThermodynamicState(sys, temperature=300*unit.kelvin)

    sampler = SamplerState(pos, box_vectors=struct.box_vectors)
    topography = Topography(topology=top, ligand_atoms=list(range(2635,2649)))
    a = BoreschBLUES(thermo, sampler, topography)
    a.determine_missing_parameters(thermo, sampler, topography)
    new_sys = a.restrain_state(thermo)
    other_sys = add_restraints(sys, struct, lig_atoms)
    struct2 = parmed.load_file('eqToluene.prmtop', xyz='posB.pdb')
    other_sys = add_restraints(other_sys, struct2, lig_atoms)
    print('new_sys', new_sys.getNumForces())
    print('new_sys', new_sys.getForces())
    print('other_sys', other_sys.getNumForces())
    print('other_sys', other_sys.getForces())

# Add the aligning dummy atoms to the solvated pdb files.
            self._add_dummy_atoms(
                index, solvate_complex.coordinate_file_path, reference_structure_path
            )

            # Extra step to create GAFF1/2 structures properly
            if isinstance(self._force_field_source, TLeapForceFieldSource):

                # Fix atom names of guest molecule for Tleap processing
                if self._paprika_setup.guest is self.taproom_guest_name:

                    gaff_version = self._force_field_source.leap_source.replace(
                        "leaprc.", ""
                    )

                    structure_mol = pmd.load_file(
                        os.path.join(
                            window_directory[: len(window_directory) - 13],
                            f"{self._paprika_setup.guest}.{gaff_version}.mol2",
                        )
                    )
                    structure_pdb = pmd.load_file(
                        self._solvated_coordinate_paths[index]
                    )

                    # Get atom names of guest molecule from restrained.pdb and *.gaff.mol2
                    mol_name = []
                    pdb_name = []
                    for original, guest in zip(
                        structure_mol,
                        structure_pdb[f":{self._paprika_setup.guest.upper()}"],
                    ):

def to_parmed(self):
        import parmed
        prmtoppath = os.path.join(tempfile.mkdtemp(), 'prmtop')
        self.prmtop.put(prmtoppath)
        pmd = parmed.load_file(prmtoppath)
        return pmd

"""
    import parmed

    # We require at least ParmEd 2.5.1 because of issues with the .mol2 writer (issue #691 on ParmEd) prior to that.
    try: #Try to get version tag
        ver = parmed.version
    except: #If too old for version tag, it is too old
        oldParmEd = Exception('ERROR: ParmEd is too old, please upgrade to 2.0.4 or later')
        raise oldParmEd
    if ver < (2,5,1):
        raise RuntimeError("ParmEd is too old, please upgrade to 2.0.4 or later")

    names = get_unique_names(len(mol2_filenames))

    for k, filename in enumerate(mol2_filenames):
        struct = parmed.load_file(filename)
        struct.name = names[k]
        mol2file = parmed.formats.Mol2File
        mol2file.write(struct, filename)

-----
        Reference for parmed.load_Structure *args and **kwargs
        https://parmed.github.io/ParmEd/html/structobj/parmed.formats.registry.load_file.html#parmed.formats.registry.load_file
        """
        if 'restart' in config['structure'].keys():
            rst7 = config['structure']['restart']
            config['Logger'].info('Restarting simulation from {}'.format(rst7))
            restart = parmed.amber.Rst7(rst7)
            config['structure'].pop('restart')

            structure = parmed.load_file(**config['structure'])
            structure.positions = restart.positions
            structure.velocities = restart.velocities
            structure.box = restart.box
        else:
            structure = parmed.load_file(**config['structure'])

        config['Structure'] = structure
        return config