How to use ParmEd - 10 common examples

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))

gaff_xml_filename = openmoltools.utils.get_data_filename("parameters/gaff.xml")
    forcefield = app.ForceField('amber99sbildn.xml', 'tip3p.xml', gaff_xml_filename)
    forcefield.registerTemplateGenerator(forcefield_generators.gaffTemplateGenerator)

    #Setup proteim
    modeller = app.Modeller(pdbfile.topology, pdbfile.positions)
    modeller.addHydrogens(forcefield, pH=7.0)

    ###Solvate and Ionize...?

    #Generate OpenMM system with FF parameters
    topology = modeller.getTopology()
    system = forcefield.createSystem( modeller.getTopology() )

    #Load System to ParmEd Structure
    structure = pmd.openmm.load_topology(topology, system)

    return structure

for (index, atom) in enumerate(mol.GetAtoms()):
        [x,y,z] = mol.GetCoords(atom)
        positions[index,0] = x*unit.angstroms
        positions[index,1] = y*unit.angstroms
        positions[index,2] = z*unit.angstroms

    #Initialize SMIRFF parameters
    # using smarty ForceField object
    smirnoff_xml_filename ="smirff99Frosst.ffxml"
    forcefield = ForceField(smirnoff_xml_filename)

    #Create OpenMM System
    system = forcefield.createSystem(topology, [mol])

    #Load System to ParmEd Structure
    structure = pmd.openmm.load_topology(topology, system)

    return structure

def __init__(self, constraints=app.HBonds, rigid_water=True, nonbondedCutoff=8.0 * unit.angstroms, use_dispersion_correction=True, nonbondedMethod=app.PME, hydrogenMass=None, switch_width=None, ewaldErrorTolerance=5E-4, **kwargs):

        TestSystem.__init__(self, **kwargs)

        try:
            from parmed.amber import AmberParm
        except ImportError as e:
            print("DHFR test system requires Parmed (`import parmed`).")
            raise(e)

        prmtop_filename = get_data_filename("data/dhfr/prmtop")
        crd_filename = get_data_filename("data/dhfr/inpcrd")

        # Initialize system.
        self.prmtop = AmberParm(prmtop_filename, crd_filename)
        system = self.prmtop.createSystem(constraints=constraints, nonbondedMethod=nonbondedMethod, rigidWater=rigid_water, nonbondedCutoff=nonbondedCutoff, hydrogenMass=hydrogenMass)

        # Extract topology
        self.topology = self.prmtop.topology

        # Set dispersion correction use.
        forces = {system.getForce(index).__class__.__name__: system.getForce(index) for index in range(system.getNumForces())}
        forces['NonbondedForce'].setUseDispersionCorrection(use_dispersion_correction)
        forces['NonbondedForce'].setEwaldErrorTolerance(ewaldErrorTolerance)

        if switch_width is not None:
            forces['NonbondedForce'].setUseSwitchingFunction(True)
            forces['NonbondedForce'].setSwitchingDistance(nonbondedCutoff - switch_width)


        positions = self.prmtop.positions

def test_atomtyping(self, mol_name):
        top_filename = '{}.top'.format(mol_name)
        gro_filename = '{}.gro'.format(mol_name)
        top_path = os.path.join(self.resource_dir, top_filename)
        gro_path = os.path.join(self.resource_dir, gro_filename)

        structure = pmd.gromacs.GromacsTopologyFile(top_path, xyz=gro_path,
                                                    parametrize=False)
        known_opls_types = [atom.type for atom in structure.atoms]

        print("Typing {}...".format(mol_name))
        typed_structure = OPLSAA.apply(structure)

        generated_opls_types = list()
        for i, atom in enumerate(typed_structure.atoms):
            message = ('Found multiple or no OPLS types for atom {} in {} ({}): {}\n'
                       'Should be atomtype: {}'.format(
                i, mol_name, top_filename, atom.type, known_opls_types[i]))
            assert atom.type, message
            generated_opls_types.append(atom.type)

        both = zip(generated_opls_types, known_opls_types)

def by_name(cls, basename, parameterized=False):

        mol2_path, top_path, gro_path = cls.available()[basename]
        if parameterized:
            # structure = pmd.gromacs.GromacsTopologyFile(top_path, xyz=gro_path,
            #                                         parametrize=parameterized)
            structure = pmd.gromacs.GromacsTopologyFile(top_path, xyz=gro_path,
                                                    parametrize=False)
            structure.title = structure.title.replace(' GAS', '')
        else:
            # structure = pmd.gromacs.GromacsTopologyFile(pdb_path, xyz=pdb_path,
            #                                         parametrize=False)
            structure = pmd.load_file(mol2_path, structure=True)
            structure.title = structure.title.replace(' GAS', '')

        return structure

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 )