Source code for Scripts.DatabankLib.settings.elements

"""
Elements assigneer functions.

Helps Molecule to define brutto formula from existing mapping files and element guesser.
"""

import json
import os
import re

import MDAnalysis as mda
import periodictable

from DatabankLib import NMLDB_DATA_PATH
from DatabankLib.core import System


[docs] def uname2element(mapping_name: str) -> str: """Guess element name from universal atom name. :param mapping_name: name of the mapping, e.g. 'M_C1', 'M_G' :raises KeyError: if element cannot be determined from mapping name :return: element name, e.g. 'C', 'N', ... """ # removes numbers and '_M' from the end of string name1 = re.sub(r"[0-9]*_M$", "", mapping_name) # removes M_X12Y23... sequence from the beginning name2 = re.sub(r"^M_([A-Z]{1,2}[0-9]{1,4})*", "", name1) # name2 is the atom and electrons are assigned to this atom if name2 == "G": # G is a glycerol carbon so change G to C name2 = "C" elif name2 in ["X", "D"]: # X is dummy, D is a drude particle return "Dummy" # Dummy is used for unknown elements in MDAnalysis try: _ = getattr(periodictable, name2).number except AttributeError as e: msg = f"This mapping name cannot be read by our rules: {mapping_name}" raise KeyError(msg) from e else: return name2
[docs] def guess_elements(system: System, u: mda.Universe) -> None: """Assign elements to atoms in the MDAnalysis Universe based on the system's composition. :param system: Databank System object :type system: System :param u: MDAnalysis Universe object :type u: mda.Universe """ for _mol, comp in system["COMPOSITION"].items(): # get the name of molecule used in simulation files mol_simu_name = comp["NAME"] mol = system.content[_mol] # if lipid is split to multiple residues _moltype = f"moltype {comp['MOLTYPE']} and " if "MOLTYPE" in comp else "" try: ua_dict_f = os.path.join(NMLDB_DATA_PATH, "lipid_json_buildH", system["UNITEDATOM_DICT"][_mol] + ".json") with open(ua_dict_f) as f: ua_dict = json.load(f) except (KeyError, TypeError): ua_dict = False for uname, props in mol.mapping_dict.items(): res = str(props.get("RESIDUE", "")) _residue = f"resname {mol_simu_name if res == '' else res} and " _name = f"name {props['ATOMNAME']}" selstr = _moltype + _residue + _name selection = u.select_atoms(selstr) if ua_dict: if props["ATOMNAME"] in ua_dict: elname = ua_dict[props["ATOMNAME"]][0] if elname == "CH": elname = "CH1" if elname == "CHdoublebond": elname = "CH1" elif selection.n_atoms == 0: continue # this is an implicit atom else: elname = uname2element(uname) else: # if not united-atom, ALL atoms must be present if selection.n_atoms == 0: msg = f"Selection '{selstr}' did not match any atoms in the universe." raise KeyError(msg) elname = uname2element(uname) selection.atoms.elements = selection.n_atoms * [elname]
# end mapping loop # end molecules loop