Polymer

mdinterface.core.polymer.Polymer

Bases: Specie

A polymer chain built from one or more repeating monomer Specie objects.

Inherits from :class:~mdinterface.core.specie.Specie. The monomers are assembled into a linear chain with :func:build_polymer, then the parent class handles topology and force-field setup. Optionally, LigParGen is called at each junction point to refine atom types and charges.

Parameters:

Name	Type	Description	Default
monomers	Specie or list of Specie	Monomer unit(s) to polymerize. A list defines a co-polymer sequence. Each monomer's ASE Atoms must carry a polymerize array marking the leaving atom (e.g. H or F) at each chain end: value 1 = head, 2 = tail. Those atoms are deleted during assembly and the bond forms between the heavy atoms they were attached to.	None
nrep	int	Number of monomer repetitions in the chain.	None
sequence	list of int	Explicit monomer sequence indices when monomers is a list (e.g. [0, 1, 0, 1] alternates two monomers).	None
refine_polymer	bool	If True, run LigParGen at every junction point to obtain accurate OPLS-AA parameters for the chain interior. Requires LigParGen.	False
offset	bool	If True, shift the total charge to match the nominal charge after topology refinement.	False
ending	str	Element symbol used to cap dangling bonds at chain termini during LigParGen snippet calculations.	"H"

Notes

Parameters not listed above (charges, atom_types, lj, cutoff, name, lammps_data, fix_missing, chg_scaling, pbc, ligpargen, tot_charge) are forwarded to :class:~mdinterface.core.specie.Specie.

Examples:

::

from mdinterface import Polymer
chain = Polymer(monomer, nrep=10)
chain = Polymer([monomer_A, monomer_B], sequence=[0,1,0,1,0,1])

Source code in mdinterface/core/polymer.py

class Polymer(Specie):
    """
    A polymer chain built from one or more repeating monomer Specie objects.

    Inherits from :class:`~mdinterface.core.specie.Specie`.  The monomers are
    assembled into a linear chain with :func:`build_polymer`, then the parent
    class handles topology and force-field setup.  Optionally, LigParGen is
    called at each junction point to refine atom types and charges.

    Parameters
    ----------
    monomers : Specie or list of Specie
        Monomer unit(s) to polymerize.  A list defines a co-polymer sequence.
        Each monomer's ASE ``Atoms`` must carry a ``polymerize`` array marking
        the leaving atom (e.g. H or F) at each chain end: value ``1`` = head,
        ``2`` = tail.  Those atoms are deleted during assembly and the bond
        forms between the heavy atoms they were attached to.
    nrep : int, optional
        Number of monomer repetitions in the chain.
    sequence : list of int, optional
        Explicit monomer sequence indices when *monomers* is a list (e.g.
        ``[0, 1, 0, 1]`` alternates two monomers).
    refine_polymer : bool, default False
        If True, run LigParGen at every junction point to obtain accurate
        OPLS-AA parameters for the chain interior.  Requires LigParGen.
    offset : bool, default False
        If True, shift the total charge to match the nominal charge after
        topology refinement.
    ending : str, default "H"
        Element symbol used to cap dangling bonds at chain termini during
        LigParGen snippet calculations.

    Notes
    -----
    Parameters not listed above (``charges``, ``atom_types``, ``lj``,
    ``cutoff``, ``name``, ``lammps_data``, ``fix_missing``, ``chg_scaling``,
    ``pbc``, ``ligpargen``, ``tot_charge``) are forwarded to
    :class:`~mdinterface.core.specie.Specie`.

    Examples
    --------
    ::

        from mdinterface import Polymer
        chain = Polymer(monomer, nrep=10)
        chain = Polymer([monomer_A, monomer_B], sequence=[0,1,0,1,0,1])
    """

    def __init__(self, monomers=None, charges=None, atom_types=None, bonds=None,
                 angles=None, dihedrals=None, impropers=None, lj={}, cutoff=1.0,
                 name=None, lammps_data=None, fix_missing=False, chg_scaling=1.0,
                 pbc=False, ligpargen=False, tot_charge=None, nrep=None,
                 sequence=None, refine_polymer=False, offset=False, ending="H"):

        # initialize polymer stuff
        self._snippet_cache = {}
        self._sequence = sequence

        # polymerize
        polymer = build_polymer(monomers, sequence=sequence, nrep=nrep)

        # Initialize the parent class with polymerized monomers
        super().__init__(atoms=polymer, charges=charges, atom_types=atom_types,
                         bonds=bonds, angles=angles, dihedrals=dihedrals,
                         impropers=impropers, lj=lj, cutoff=cutoff, name=name,
                         lammps_data=lammps_data, fix_missing=fix_missing,
                         chg_scaling=chg_scaling, pbc=pbc, ligpargen=ligpargen,
                         tot_charge=tot_charge)

        if refine_polymer:
            self.refine_polymer_topology(Nmax=12, offset=offset, ending=ending)

        return

    # return list of elements adjacent to a connection point
    def _get_connection_elements(self):

        # Get elements where there is a connection
        centers = np.argwhere((self.atoms.arrays["is_connected"] ==1) |
                              (self.atoms.arrays["is_connected"] ==2)).flatten()

        pairs = []
        for center in centers:
            if any(center in pp for pp in pairs):
                continue
            distances = self.atoms.get_distances(center, centers)
            idx = centers[np.argsort(distances)[1]]
            pairs.append([center, idx])

        return pairs

    def _get_start_end(self):

        str_idx = np.argwhere(self.atoms.arrays["is_connected"] == -1).flatten()
        end_idx = np.argwhere(self.atoms.arrays["is_connected"] == -2).flatten()

        return np.concatenate([str_idx, end_idx])


    def _update_connection(self, center, partner, Nmax, charges, ending="H"):

        # make a lil snippet
        snippet, snippet_idxs = make_snippet(self, center, Nmax, ending=ending)

        # get local indexes (within dihedral from center)
        ldxs = list(set(np.concatenate(self.find_relevant_distances(4, centers=center))))

        # find mapping between indexes
        mapping = [np.argwhere(snippet_idxs == ll)[0][0] for ll in ldxs]

        # Check if snippet already exists in cache
        snippet_hash = ''.join(snippet.get_chemical_symbols())
        if snippet_hash in self._snippet_cache:

            cached_snippet = copy.deepcopy(self._snippet_cache[snippet_hash])

            sn_atoms     = cached_snippet["sn_atoms"]
            sn_atypes    = cached_snippet["sn_atypes"]
            sn_bonds     = cached_snippet["sn_bonds"]
            sn_angles    = cached_snippet["sn_angles"]
            sn_dihedrals = cached_snippet["sn_dihedrals"]
            sn_impropers = cached_snippet["sn_impropers"]
            new_charges  = cached_snippet["new_charges"]

        # if not, ligpargen it
        else:

            # snippet charge
            if "nominal_charge" in snippet.arrays:
                sn_charge = snippet.arrays["nominal_charge"].sum()
            else:
                raise ValueError("No 'nominal_charge' found in snippet!")
                # sn_charge = None

            # run ligpargen
            sn_atoms, sn_atypes, sn_bonds, sn_angles, sn_dihedrals, sn_impropers =\
                run_ligpargen(snippet, charge=sn_charge, is_snippet=True)

            # get charges
            new_charges = sn_atoms.get_initial_charges()

            # Store the new snippet and its charges in cache
            self._snippet_cache[snippet_hash] = {
                "new_charges"  : new_charges,
                "sn_atoms"     : sn_atoms,
                "sn_atypes"    : sn_atypes,
                "sn_bonds"     : sn_bonds,
                "sn_angles"    : sn_angles,
                "sn_dihedrals" : sn_dihedrals,
                "sn_impropers" : sn_impropers,
                "snippet_idxs" : snippet_idxs,
                "local_idxs"   : ldxs,
                "snippet"      : snippet,
                "sn_charge"    : sn_charge
                }

        # update topology of section
        original_idxs = self._sids[snippet_idxs]
        local_idxs = self._sids[ldxs]

        # remap topology of local indexes
        new_bonds, new_angles, new_dihedrals, new_impropers = remap_snippet_topology(
            original_idxs, sn_atoms, sn_atypes, sn_bonds, sn_angles, sn_dihedrals,
            sn_impropers, local_idxs)

        charges[ldxs] = new_charges[mapping]

        # Correct OPLS types for the two junction atoms (they had H caps
        # substituting their bonded partner during monomer LigParGen runs).
        self._update_junction_lj_types([center, partner], snippet_idxs, sn_atypes)

        self._add_to_topology(bonds=new_bonds, angles=new_angles,
                              dihedrals=new_dihedrals, impropers=new_impropers)

        return

    # main driver that refines charges across the whole polymer
    def refine_polymer_topology(self, Nmax=12, offset=False, ending="H"):
        """
        Refines the charges for the specified species.

        Parameters:
        specie (ase.Atoms): The species object containing the atoms.
        Nmax (int): The maximum number of neighbors to consider.
        offset (bool): Whether to apply an offset to the charges.

        Returns:
        np.ndarray: The refined charges.
        """

        # Clean topology first to remove any invalid interactions from polymerization
        self._cleanup_topology()

        # get charges and connection elements
        charges = self.charges
        pairs = self._get_connection_elements()

        str_end_idxs = self._get_start_end()

        # get charges at every point
        for pair in pairs:
            ci, cj = int(pair[0]), int(pair[1])
            self._update_connection(cj, ci, Nmax, charges, ending=ending)
        # for center in str_end_idxs:
            # self._update_connection(center, Nmax, charges, ending=ending)

        # bring back to zero
        if offset:

            if "nominal_charge" in self.atoms.arrays:
                target = self.atoms.arrays["nominal_charge"].sum()
            else:
                raise Warning("No nominal charge found, assume it is 0.")

            charges -= ((charges.sum() - target) / len(charges))

        self.atoms.set_initial_charges(charges)
        return

refine_polymer_topology(Nmax=12, offset=False, ending='H')

Refines the charges for the specified species.

Parameters: specie (ase.Atoms): The species object containing the atoms. Nmax (int): The maximum number of neighbors to consider. offset (bool): Whether to apply an offset to the charges.

Returns: np.ndarray: The refined charges.

Source code in mdinterface/core/polymer.py

def refine_polymer_topology(self, Nmax=12, offset=False, ending="H"):
    """
    Refines the charges for the specified species.

    Parameters:
    specie (ase.Atoms): The species object containing the atoms.
    Nmax (int): The maximum number of neighbors to consider.
    offset (bool): Whether to apply an offset to the charges.

    Returns:
    np.ndarray: The refined charges.
    """

    # Clean topology first to remove any invalid interactions from polymerization
    self._cleanup_topology()

    # get charges and connection elements
    charges = self.charges
    pairs = self._get_connection_elements()

    str_end_idxs = self._get_start_end()

    # get charges at every point
    for pair in pairs:
        ci, cj = int(pair[0]), int(pair[1])
        self._update_connection(cj, ci, Nmax, charges, ending=ending)
    # for center in str_end_idxs:
        # self._update_connection(center, Nmax, charges, ending=ending)

    # bring back to zero
    if offset:

        if "nominal_charge" in self.atoms.arrays:
            target = self.atoms.arrays["nominal_charge"].sum()
        else:
            raise Warning("No nominal charge found, assume it is 0.")

        charges -= ((charges.sum() - target) / len(charges))

    self.atoms.set_initial_charges(charges)
    return