API reference

Signatures, documented defaults, and expandable source below are generated from this checkout on every build. Start with the linked tutorials for worked examples, interpretation, and known limitations. Public methods lacking original docstrings are included so they are discoverable.

Private modules beginning with _ are implementation details. Names accidentally re-exported by wildcard imports (for example NumPy or imported helpers) are not a supported matmalib API. Import the classes and functions below explicitly.

Molecular objects

Tutorial: Molecules and geometry. Compatibility notes: Core behavior.

matmalib.core.Atom(coord: list = None, type: str = None, valence: int = None, index: int = None)

Source code in matmalib/core.py

def __init__(self, coord:list = None, type:str = None, valence:int = None, index:int = None):
    self.coord = np.array(coord)
    self.type = type
    self.id = index
    self.mass = None
    self.atoms = None

    if self.coord is not None:
        self.x = self.coord[0]
        self.y = self.coord[1]
        self.z = self.coord[2]

coord = np.array(coord) instance-attribute

type = type instance-attribute

id = index instance-attribute

mass = None instance-attribute

atoms = None instance-attribute

x = self.coord[0] instance-attribute

y = self.coord[1] instance-attribute

z = self.coord[2] instance-attribute

matmalib.core.Fragment()

Source code in matmalib/core.py

def __init__(self):
    self.atoms = []

atoms = [] instance-attribute

xyz property

Mass property

matmalib.core.State(**kwargs)

Source code in matmalib/core.py

def __init__(self, **kwargs):
    self.__dict__.update(kwargs)

xyz property

Mass property

mass property

from_file(path, software) classmethod

Source code in matmalib/core.py

@classmethod
def from_file(cls, path, software):

    builder = StateBuilder(path=path, software=software)
    system =  builder.build()
    system.id = os.path.basename(os.path.dirname(path)) or os.path.basename(path)
    return system

element_symbol(A) staticmethod

Return the symbol for atomic number 1–118 (integer or string).

Source code in matmalib/core.py

@staticmethod
def element_symbol(A):
    """Return the symbol for atomic number 1–118 (integer or string)."""
    return _ELEMENT_SYMBOLS[str(A)]

add_atoms(atom)

Source code in matmalib/core.py

def add_atoms(self, atom):
    self.atoms.append(atom)

remove_atoms(atom)

Source code in matmalib/core.py

def remove_atoms(self, atom):
    self.atoms.remove(atom)

measure_internal(atom_list)

Source code in matmalib/core.py

def measure_internal(self, atom_list):
    return StateEditor.measure(atom_list=atom_list, state=self)

set_internal(atom_list, val)

Source code in matmalib/core.py

def set_internal(self, atom_list, val):

    editor = StateEditor(self)

    if len(atom_list) == 2:
        editor.set_distance(atom_list=atom_list, new_distance=val)        
    elif len(atom_list) == 3:
        editor.set_angle(atom_list=atom_list, new_ang=val)
    elif len(atom_list) == 4:
        editor.set_dihedral(atom_list=atom_list, new_dih=val)

connectivity_matrix(distXX=1.65, distXH=1.15)

Creates a connectivity matrix of the molecule. A connectivity matrix holds the information of which atoms are bonded and to what.

Parameters:

Name	Type	Description	Default
distXX		The max distance between two atoms (not hydrogen) to be considered a bond	1.65
distXH		The max distance between any atom and a hydrogen atom to be considered a bond	1.15

Source code in matmalib/core.py

def connectivity_matrix(self, distXX=1.65, distXH=1.15):

    """ Creates a connectivity matrix of the molecule. A connectivity matrix holds the information of which atoms are bonded and to what.

    :param distXX: The max distance between two atoms (not hydrogen) to be considered a bond
    :param distXH: The max distance between any atom and a hydrogen atom to be considered a bond
    """

    Nat = self.NAtoms
    self.conn_mat = np.zeros((Nat, Nat))

    for at1 in self.atoms:
        for at2 in self.atoms:

            dist = StateEditor.get_distance([at1,at2])

            if at1 == at2:
                pass

            elif (at1.type == 'H' or at2.type == 'H') and dist < distXH:
                self.conn_mat[at1.id, at2.id] = 1;
                self.conn_mat[at2.id, at1.id] = 1
            elif (at1.type != 'H' and at2.type != 'H') and dist < distXX:
                self.conn_mat[at1.id, at2.id] = 1;
                self.conn_mat[at2.id, at1.id] = 1

    # Remove bifurcated Hs:
    for at1 in self.atoms:
        if at1.type == 'H' and np.sum(self.conn_mat[at1.id, :]) > 1:

            at2list = np.where(self.conn_mat[at1.id, :] == 1)
            at2list = at2list[0].tolist()

            at2dist = [round(StateEditor.get_distance([at1, self.atoms[at2x]]), 3) for at2x in at2list]
            for at, dist in zip(at2list, at2dist):
                if self.atoms[at].type == 'H':
                    at2list.remove(at)
                    at2dist.remove(dist)

            at2 = at2list[at2dist.id(min(at2dist))]
            for at2x in at2list:
                if at2x != at2:
                    print('remove', self.id, at2x, at1, at2)
                    self.conn_mat[at1.id, at2x] = 0;
                    self.conn_mat[at2x, at1.id] = 0

    graph = nx.from_numpy_array(self.conn_mat)

    self.Nmols = nx.number_connected_components(graph)
    self.graph = graph

system_split()

Source code in matmalib/core.py

def system_split(self):

    g = self.graph
    components = [g.subgraph(c).copy() for c in nx.connected_components(g)]

    frags = []
    for comp in components:
        frag = Fragment()
        frag.graph = comp

        atom_indices = comp.nodes()
        for index in atom_indices:
            atom = self.atoms[index]
            frag.atoms.append(atom)

        frags.append(frag)

    self.mols = frags

show(width=600, height=600)

Source code in matmalib/core.py

def show(self, width=600, height=600):

    import py3Dmol as p3D

    """ Displays a 3D rendering of the conformer using Py3Dmol

    :param width: the width of the display window
    :param height: the height of the display window
    """

    XYZ = "{0:3d}\n{1:s}\n".format(self.NAtoms, self.id)
    for at, xyz in zip(self.atoms, self.xyz):
        XYZ += "{0:3s}{1:10.3f}{2:10.3f}{3:10.3f}\n".format(at.type, xyz[0], xyz[1], xyz[2])
    xyzview = p3D.view(width=width, height=height)
    xyzview.addModel(XYZ, 'xyz')
    xyzview.setStyle({'stick': {}})
    xyzview.zoomTo()
    xyzview.show()

matmalib.core.Data(**kwargs)

Source code in matmalib/core.py

def __init__(self, **kwargs):

    self.__dict__.update(kwargs)

from_file(path, software, **kwargs) classmethod

Source code in matmalib/core.py

@classmethod
def from_file(cls, path, software, **kwargs):

    builder = DataBuilder(path=path, software=software, **kwargs)

    return builder.build()

matmalib.core.Space(states: list = None, theory: dict = None, input_path: str = None, output_path: str = None)

Source code in matmalib/core.py

def __init__(self, states:list = None, theory:dict = None, input_path:str = None, output_path:str = None):
    self.states = [] if states is None else states

    if theory is None:
        theory = {
            'disp': True,
            'nprocs': 24,
            'mem': '60B',
            'method': 'PBE1PBE',
            'basis_set': '6-311G**',
            'jobtype': 'opfreq',
            'charge': 0,
            'multiplicity': 1
        }
    self.theory = theory
    self.input_path = input_path
    self.output_path = output_path

states = [] if states is None else states instance-attribute

theory = theory instance-attribute

input_path = input_path instance-attribute

output_path = output_path instance-attribute

energies property

free_energies property

enthalpies property

from_xyz(path, *, glycans=False) classmethod

Read every XYZ frame into a State (or a Glycan when glycans=True).

Source code in matmalib/core.py

@classmethod
def from_xyz(cls, path, *, glycans=False):
    """Read every XYZ frame into a State (or a Glycan when glycans=True)."""
    from .glycan import Glycan, _xyz_frames
    states = list(_xyz_frames(path))
    if not states:
        raise ValueError("The trajectory contains no frames.")
    for i, state in enumerate(states):
        state.id = str(i)
        state.path = str(path)
    if glycans:
        states = [Glycan.from_state(state) for state in states]
    return cls(states=states)

rmsd_matrix(**kwargs)

Return aligned pairwise RMSDs; heavy atoms are used by default.

Source code in matmalib/core.py

def rmsd_matrix(self, **kwargs):
    """Return aligned pairwise RMSDs; heavy atoms are used by default."""
    from .conformers import rmsd_matrix
    return rmsd_matrix(self, **kwargs)

cluster(threshold=1.0, **kwargs)

Return a ClusterResult with memberships and representative conformers.

Source code in matmalib/core.py

def cluster(self, threshold=1.0, **kwargs):
    """Return a ClusterResult with memberships and representative conformers."""
    from .conformers import cluster_conformers
    return cluster_conformers(self, threshold, **kwargs)

remove_duplicates(threshold=0.1, **kwargs)

Return a new Space of complete-linkage representatives; keep self intact.

Every pair within a duplicate group is at most threshold apart. Other options match cluster(), except the linkage method is fixed.

Source code in matmalib/core.py

def remove_duplicates(self, threshold=0.1, **kwargs):
    """Return a new Space of complete-linkage representatives; keep self intact.

    Every pair within a duplicate group is at most threshold apart.
    Other options match cluster(), except the linkage method is fixed.
    """
    return self.cluster(threshold, method="complete", **kwargs).representatives

to_file(software: str = 'gaussian')

Source code in matmalib/core.py

def to_file(self, software:str = 'gaussian'):
    path = self.output_path
    theory = self.theory

    builder = InputBuilder(path=path, theory=theory, software=software, space=self)
    builder.build(path=path, theory=theory, software=software, space=self)

sum(states: list, quant: str = 'E')

Source code in matmalib/core.py

def sum(self, states:list, quant:str = 'E'):

    if isinstance(states, int):
        states_list = [state for state in self.states]
    else:
        states_list = states

    if quant not in ['E', 'F', 'H']:
        print('Error: Quantity not recognized')
        return None

    total = 0
    for state in states_list:
        if not hasattr(state, quant):
            print('Error: State does not have attribute', quant)
            return None

        val = getattr(state, quant)
        total += val

    return total

from_dir(path, software='gaussian') classmethod

Source code in matmalib/core.py

@classmethod
def from_dir(cls, path, software='gaussian'):

    dirs = [d for d in os.listdir(path) if os.path.isdir(f"{path}/{d}")]
    states = []

    for subdir in dirs:

        if software == 'gaussian':
            state = State.from_file(f"{path}/{subdir}/input.log", software='gaussian')
        # Include other software formats in the future.

        states.append(state)

    return Space(states)

diff(states=None, quant: str = 'E')

Source code in matmalib/core.py

def diff(self, states=None, quant:str = 'E'):

    if states:
        if isinstance(states, list):

            if isinstance(states[0], int):
                states_list = [self.states[n] for n in states]
            else:
                states_list = states

    else:
        states_list = self.states

    if quant not in ['E', 'F', 'H']:
        raise ValueError(f"Error: Quantity {quant} not recognized!")

    for state in states_list:
        if not hasattr(state, quant):
            raise AttributeError(f"Error: State does not have attribute {quant}")

    diff_list = []
    # val_0 = getattr(states_list[0], quant)

    vals = [getattr(state, quant) for state in states_list]
    val_0 = np.min(vals)

    for state in states_list:
        val = getattr(state, quant) - val_0
        diff_list.append(val * 627.5095)

    return diff_list

add(state)

Source code in matmalib/core.py

def add(self, state):

    if isinstance(state, list):
        for st in state:
            self.states.append(st)
    else:
        self.states.append(state)

remove(state)

Source code in matmalib/core.py

def remove(self, state):
    type = type(state)
    if type == list:
        for st in state:
            if type(st) != int:
                self.states.remove(st)
            else:
                self.states.pop(st)

    elif type == int:
        self.states.pop(state)
    else:
        self.states.remove(state)

sort(val: str = 'E')

Source code in matmalib/core.py

def sort(self, val:str = 'E'):
    try:
        self.states.sort(key=lambda x: getattr(x, val))
    except:
        print(f"Error Attribute {val} not recognized")

table(sort='E')

Source code in matmalib/core.py

def table(self, sort='E'):
    from tabulate import tabulate
    try:
        self.states.sort(key=lambda x: getattr(x, sort))
        state_array = []
        for i, state in enumerate(self.states):
            row = [i, state.id]
            if hasattr(state, 'E'):
                row.append(state.E)
            if hasattr(state, 'F'):
                row.append(state.F)
            if hasattr(state, 'H'):
                row.append(state.H)
            state_array.append(row)
        print(tabulate(state_array, headers=['Index', 'Name', 'Energy', 'Free Energy', 'Enthalpy'], tablefmt='fancy_grid'))
    except:
        print(f"Error Attribute {sort} not recognized")

Geometry and file builders

Tutorials: Geometry editing, File I/O, and Ensemble energies. InputBuilder may be imported explicitly from matmalib.calc.

matmalib.util.StateEditor(state)

Source code in matmalib/util.py

def __init__(self, state):

    self.state = state

state = state instance-attribute

get_distance(atom_list, state=None) staticmethod

Finds the distance between two atoms

Parameters:

Name	Type	Description	Default
atom_list		Two Atom objects or two zero-based atom positions.	required
state		State used to resolve integer atom positions.	None

Returns:

Type	Description
	(float) the distance between 2 atoms

Source code in matmalib/util.py

@staticmethod
def get_distance(atom_list, state=None):

    """ Finds the distance between two atoms
    :param atom_list: Two Atom objects or two zero-based atom positions.
    :param state: State used to resolve integer atom positions.
    :return: (float) the distance between 2 atoms
    """

    from .core import Atom

    # If at1 or at2 are integers, extract these values from the supplied state.
    if len(atom_list) != 2: print(
        f"Warning! {len(atom_list)} atoms specified for function get_distance, expected 2.")

    if all(isinstance(at, int) for at in atom_list) and state != None:
        at1 = state.atoms[atom_list[0]]
        at2 = state.atoms[atom_list[1]]

    elif all(isinstance(at, Atom) for at in atom_list):
        at1 = atom_list[0]
        at2 = atom_list[1]

    else:
        raise ValueError('Error: Specify all Atoms, or all integers for atom list!')

    dist = 0
    for dim in range(3):
        dist += (at1.coord[dim] - at2.coord[dim]) ** 2

    return math.sqrt(dist)

adjacent_atoms(conn_mat, at) staticmethod

returns all adjacent atoms to a specific atom in a conformation

Parameters:

Name	Type	Description	Default
conn_mat		the connectivity matrix	required
at		index of the selected atom	required

Returns:

Type	Description
	indices of all atoms adjacent to the selected atom

Source code in matmalib/util.py

@staticmethod
def adjacent_atoms(conn_mat, at):

    """returns all adjacent atoms to a specific atom in a conformation
    :param conn_mat: the connectivity matrix
    :param at: index of the selected atom
    :return: indices of all atoms adjacent to the selected atom
    """

    return np.nonzero(conn_mat[at, :])[0]

get_angle(atom_list, state=None) staticmethod

Calculate an angle between three atoms: angle = acos(dot(X,Y)/(norm(X)*norm(Y)))

Parameters:

Name	Type	Description	Default
atom_list		Three Atom objects or three zero-based atom positions.	required
state		State used to resolve integer atom positions.	None

Returns:

Type	Description
	A tuple of angle in degrees and rotation-axis vector.

Source code in matmalib/util.py

@staticmethod
def get_angle(atom_list, state=None):

    """Calculate an angle between three atoms:
    angle = acos(dot(X,Y)/(norm(X)*norm(Y)))
    :param atom_list: Three Atom objects or three zero-based atom positions.
    :param state: State used to resolve integer atom positions.
    :return: A tuple of angle in degrees and rotation-axis vector.
    """

    from .core import Atom

    if len(atom_list) != 3: print(
        f"Warning! {len(atom_list)} atoms specified for function get_distance, expected 3.")

    if all(isinstance(at, int) for at in atom_list) and state != None:
        at1 = state.atoms[atom_list[0]]
        at2 = state.atoms[atom_list[1]]
        at3 = state.atoms[atom_list[2]]

    elif all(isinstance(at, Atom) for at in atom_list):
        at1 = atom_list[0]
        at2 = atom_list[1]
        at3 = atom_list[2]

    else:
        raise ValueError('Error: Specify all Atoms, or all integers for atom list!')

    r1 = [x - y for x, y in zip(at1.coord, at2.coord)]
    r2 = [x - y for x, y in zip(at3.coord, at2.coord)]

    norm1 = StateEditor.get_distance([at1, at2])
    norm2 = StateEditor.get_distance([at2, at3])

    norm = norm1 * norm2
    dot = np.dot(r1, r2) / norm
    angle = math.acos(dot)
    axor = np.cross(r1, r2)

    return angle * 180 / np.pi, axor

get_dihedral(atom_list, state=None, ax: bool = False) staticmethod

Measure the signed torsion of four ordered atoms, in degrees.

Uses the same atan2 convention as glycan torsions: a positive right-handed rotation of the fourth atom about the second-to-third atom direction increases the angle. Values lie in [-180, 180]. This corrects the opposite sign and translation-dependent sign bug in the legacy implementation.

Parameters:

Name	Type	Description	Default
atom_list		Four Atom objects or four zero-based atom positions.	required
state		State used to resolve integer atom positions.	None
ax	bool	Also return the unnormalized cross product of the two plane normals, retained for compatibility. It is zero for exactly planar torsions; it is not a stable rotation axis.	False

Returns:

Type	Description
	Signed degrees, or (degrees, plane-normal cross product).

Raises:

Type	Description
ValueError	For invalid selections, nonfinite coordinates, coincident central atoms, or collinear triples.

Source code in matmalib/util.py

@staticmethod
def get_dihedral(atom_list, state=None, ax: bool = False):
    """Measure the signed torsion of four ordered atoms, in degrees.

    Uses the same atan2 convention as glycan torsions: a positive
    right-handed rotation of the fourth atom about the second-to-third
    atom direction increases the angle. Values lie in [-180, 180].
    This corrects the opposite sign and translation-dependent sign bug
    in the legacy implementation.

    :param atom_list: Four Atom objects or four zero-based atom positions.
    :param state: State used to resolve integer atom positions.
    :param ax: Also return the unnormalized cross product of the two
        plane normals, retained for compatibility. It is zero for
        exactly planar torsions; it is not a stable rotation axis.
    :returns: Signed degrees, or (degrees, plane-normal cross product).
    :raises ValueError: For invalid selections, nonfinite coordinates,
        coincident central atoms, or collinear triples.
    """

    from .core import Atom

    if len(atom_list) != 4:
        raise ValueError("A dihedral requires exactly four atoms.")
    if all(isinstance(at, _Integral) for at in atom_list) and state is not None:
        if any(at < 0 or at >= len(state.atoms) for at in atom_list):
            raise ValueError("Dihedral atom positions are out of range.")
        atoms = [state.atoms[at] for at in atom_list]
    elif all(isinstance(at, Atom) for at in atom_list):
        atoms = atom_list
    else:
        raise ValueError("Specify all Atom objects, or all integer positions with a State.")

    coordinates = np.asarray([at.coord for at in atoms], dtype=float)
    angle = float(_dihedral(coordinates))
    if ax:
        bonds = np.diff(coordinates, axis=0)
        normals = (np.cross(bonds[0], bonds[1]), np.cross(bonds[1], bonds[2]))
        return angle, np.cross(*normals)
    return angle

measure(atom_list, state=None) staticmethod

Source code in matmalib/util.py

@staticmethod
def measure(atom_list, state=None):

    if len(atom_list) == 2:
        return StateEditor.get_distance(atom_list, state=state)
    elif len(atom_list) == 3:
        return StateEditor.get_angle(atom_list, state=state)
    elif len(atom_list) == 4:
        return StateEditor.get_dihedral(atom_list, state=state)

determine_carried_atoms(state, atom_list) staticmethod

Source code in matmalib/util.py

@staticmethod
def determine_carried_atoms(state, atom_list):

    import networkx as nx
    from .core import Atom

    """Find all atoms necessary to be carried over during rotation
    of an atom
    :param at1: (list) the xyz coordinates of an atom
    :param at2: (list) the xyz coordinates of another atom
    :param conn_matt: the connectivity matrix of a conformer
    """
    #   1. Zero the connections in connectivity matrix
    # tmp_conn = np.copy(conf.conn_mat)

    if len(atom_list) != 2: print(
        f"Warning! {len(atom_list)} atoms specified for function get_distance, expected 2.")

    if all(isinstance(at, int) for at in atom_list) and state != None:
        at1 = state.atoms[atom_list[0]]
        at2 = state.atoms[atom_list[1]]

    elif all(isinstance(at, Atom) for at in atom_list):
        at1 = atom_list[0]
        at2 = atom_list[1]

    else:
        raise ValueError('Error: Specify all Atoms, or all integers for atom list!')

    coord1 = at1.coord
    coord2 = at2.coord

    broke_bond = False
    if state.conn_mat[at1.id, at2.id] == 1 or state.conn_mat[at2.id, at1.id] == 1:
        StateEditor.disconnect_atoms(state, at1.id, at2.id)
        broke_bond = True

    # cm = state.graph # This doesn't work.
    cm = nx.graph.Graph(state.conn_mat)

    if nx.is_connected(cm):
        print("Matrix still connected")
        return

    #   2. Determine the connected atoms:
    for subgraph in nx.connected_components(cm):
        if at2.id in subgraph:
            # print('subgraph:',subgraph)
            carried_atoms = [state.atoms[at] for at in subgraph]
    if broke_bond: StateEditor.connect_atoms(state, at1.id, at2.id)

    return carried_atoms

connect_atoms(state, at1: int, at2: int) staticmethod

Connects two atoms in the conn_mat

Parameters:

Name	Type	Description	Default
state		State with an associated conn_mat.	required
at1	int	First zero-based matrix index.	required
at2	int	Second zero-based matrix index.	required

Returns:

Type	Description
	None

Source code in matmalib/util.py

@staticmethod
def connect_atoms(state, at1: int, at2: int):

    """Connects two atoms in the conn_mat
    :param state: State with an associated conn_mat.
    :param at1: First zero-based matrix index.
    :param at2: Second zero-based matrix index.
    :return: None
    """

    state.conn_mat[at1, at2] = 1
    state.conn_mat[at2, at1] = 1

disconnect_atoms(state, at1: int, at2: int) staticmethod

Disconnects two atoms in the conn_mat

Parameters:

Name	Type	Description	Default
state		State with an associated conn_mat.	required
at1	int	First zero-based matrix index.	required
at2	int	Second zero-based matrix index.	required

Returns:

Type	Description
	None

Source code in matmalib/util.py

@staticmethod
def disconnect_atoms(state, at1: int, at2: int):

    """Disconnects two atoms in the conn_mat
    :param state: State with an associated conn_mat.
    :param at1: First zero-based matrix index.
    :param at2: Second zero-based matrix index.
    :return: None
    """

    state.conn_mat[at1, at2] = 0
    state.conn_mat[at2, at1] = 0

set_distance(atom_list, new_distance)

Set a new distance between two atoms

Parameters:

Name	Type	Description	Default
atom_list		Two zero-based atom positions.	required
new_distance		value of bond distance to be set	required

Source code in matmalib/util.py

def set_distance(self, atom_list, new_distance):

    """Set a new distance between two atoms
    :param atom_list: Two zero-based atom positions.
    :param new_distance: value of bond distance to be set
    """

    if len(atom_list) != 2 or not all(isinstance(at, int) for at in atom_list):
        raise ValueError("The position needs to be defined by 2 integers")

    at1 = self.state.atoms[atom_list[0]]
    at2 = self.state.atoms[atom_list[1]]

    # Identify displacement vector:

    displacement_vector = at2.coord - at1.coord
    norm_vector = np.sqrt(np.sum(displacement_vector ** 2))
    normalized_vector = displacement_vector / norm_vector

    old_distance = StateEditor.get_distance([at1, at2], state=self.state)
    translation = normalized_vector * (new_distance - old_distance)

    carried_atoms = StateEditor.determine_carried_atoms(self.state, [at1.id, at2.id])

    for at in carried_atoms:
        self.state.atoms[at.id].coord += translation

set_angle(atom_list, new_ang)

Set a new angle between three atoms

Parameters:

Name	Type	Description	Default
atom_list		Three zero-based atom positions.	required
new_ang		Bond angle in degrees to be set.	required

Source code in matmalib/util.py

def set_angle(self, atom_list, new_ang):

    """Set a new angle between three atoms
    :param atom_list: Three zero-based atom positions.
    :param new_ang: Bond angle in degrees to be set.
    """

    if len(atom_list) != 3 or not all(isinstance(at, int) for at in atom_list):
        raise ValueError("The position needs to be defined by 3 integers")

    at1 = self.state.atoms[atom_list[0]]
    at2 = self.state.atoms[atom_list[1]]
    at3 = self.state.atoms[atom_list[2]]

    #   Determine the axis of rotation:

    old_ang, axor = StateEditor.get_angle([at1, at2, at3], state=self.state)
    norm_axor = np.sqrt(np.sum(axor ** 2))
    normalized_axor = axor / norm_axor

    #   Each carried_atom is rotated by euler-rodrigues formula:
    #   Also, I move the midpoint of the bond to the mid atom
    #   the rotation step and then move the atom back.

    rot_angle = np.pi * (new_ang - old_ang) / 180.
    translation = at2.coord

    # apply rotations to at3.
    rot = expm(np.cross(np.eye(3), normalized_axor * (rot_angle)))
    carried_atoms = StateEditor.determine_carried_atoms(self.state, [at2.id, at3.id])

    for at in carried_atoms:
        at.coord = np.dot(rot, at.coord - translation)
        at.coord = at.coord + translation

set_dihedral(atom_list, new_dih, incr=False, axis_pos='bond', threshold=0.1)

Rotate the moving fragment using the signed dihedral convention.

In bond mode, rotate the component on the third atom's side of the central bond about the second-to-third atom direction. The first atom must remain fixed and the fourth must belong to that component. Cutting the bond must separate its endpoints; rings are rejected. Connectivity is never modified. Exact planar targets are supported.

The legacy terminal mode rotates only atoms beyond the fourth atom about a parallel axis through that atom. It leaves the defining quartet fixed, so it does not set that quartet's measured torsion.

Parameters:

Name	Type	Description	Default
atom_list		Four zero-based atom positions.	required
new_dih		Finite target angle in degrees, or an increment.	required
incr		Treat new_dih as an increment to the current dihedral.	False
axis_pos		Use the central bond ('bond') or terminal pivot ('term').	'bond'
threshold		Finite nonnegative tolerance in degrees for a no-op.	0.1

Returns:

Type	Description
	None; modifies the State's atom coordinates in place.

Raises:

Type	Description
ValueError	For invalid inputs, degenerate geometry, or a bond that does not separate the required fixed and moving atoms.

Source code in matmalib/util.py

def set_dihedral(self, atom_list, new_dih, incr=False, axis_pos="bond", threshold=0.1):
    """Rotate the moving fragment using the signed dihedral convention.

    In bond mode, rotate the component on the third atom's side of the
    central bond about the second-to-third atom direction. The first
    atom must remain fixed and the fourth must belong to that component.
    Cutting the bond must separate its endpoints; rings are rejected.
    Connectivity is never modified. Exact planar targets are supported.

    The legacy terminal mode rotates only atoms beyond the fourth atom
    about a parallel axis through that atom. It leaves the defining
    quartet fixed, so it does not set that quartet's measured torsion.

    :param atom_list: Four zero-based atom positions.
    :param new_dih: Finite target angle in degrees, or an increment.
    :param incr: Treat new_dih as an increment to the current dihedral.
    :param axis_pos: Use the central bond ('bond') or terminal pivot ('term').
    :param threshold: Finite nonnegative tolerance in degrees for a no-op.
    :returns: None; modifies the State's atom coordinates in place.
    :raises ValueError: For invalid inputs, degenerate geometry, or a
        bond that does not separate the required fixed and moving atoms.
    """
    import networkx as nx

    if len(atom_list) != 4 or not all(isinstance(at, _Integral) for at in atom_list):
        raise ValueError("A dihedral must be defined by four integer atom positions.")
    if len(set(atom_list)) != 4 or any(at < 0 or at >= len(self.state.atoms) for at in atom_list):
        raise ValueError("Supply four distinct atom positions within the State.")
    if axis_pos not in ("bond", "term"):
        raise ValueError("axis_pos must be 'bond' or 'term'.")
    if not isinstance(new_dih, _Real) or not np.isfinite(new_dih):
        raise ValueError("The dihedral target must be a finite real number.")
    if not isinstance(threshold, _Real) or not np.isfinite(threshold) or threshold < 0:
        raise ValueError("threshold must be finite and nonnegative.")

    old_dih = self.get_dihedral(atom_list, state=self.state)
    delta = float(new_dih % 360.0)
    if not incr:
        delta -= old_dih
    delta = (delta + 180.0) % 360.0 - 180.0
    if abs(delta) <= threshold:
        return

    atoms = [self.state.atoms[at] for at in atom_list]
    axis = np.asarray(atoms[2].coord, dtype=float) - atoms[1].coord
    axis /= np.linalg.norm(axis)
    translation = ((atoms[1].coord + atoms[2].coord) / 2
                   if axis_pos == "bond" else atoms[3].coord.copy())

    matrix = np.asarray(getattr(self.state, "conn_mat", None), dtype=float)
    size = len(self.state.atoms)
    if (matrix.shape != (size, size) or not np.all(np.isfinite(matrix))
            or not np.array_equal(matrix, matrix.T)):
        raise ValueError("Supply a finite symmetric connectivity matrix indexed by atom position.")
    graph = nx.from_numpy_array(matrix)
    fixed, moving = atom_list[1:3] if axis_pos == "bond" else atom_list[2:4]
    if not graph.has_edge(fixed, moving):
        raise ValueError("The selected rotation bond is missing from the connectivity matrix.")
    graph.remove_edge(fixed, moving)
    carried = nx.node_connected_component(graph, moving)
    if fixed in carried:
        raise ValueError("Cannot rotate a ring bond: cutting it does not separate the atoms.")
    if axis_pos == "bond":
        if atom_list[0] in carried or atom_list[3] not in carried:
            raise ValueError("The first and fourth atoms must lie on opposite sides of the rotation bond.")
    else:
        carried.remove(atom_list[3])
        if carried.intersection(atom_list):
            raise ValueError("Terminal rotation must leave the defining quartet fixed.")
    if not carried:
        return

    carried = sorted(carried)
    coordinates = np.asarray([self.state.atoms[at].coord for at in carried], dtype=float)
    if coordinates.shape != (len(carried), 3) or not np.all(np.isfinite(coordinates)):
        raise ValueError("The moving atoms must have finite three-dimensional coordinates.")
    rotation = expm(np.cross(np.eye(3), axis * np.radians(delta)))
    rotated = (coordinates - translation) @ rotation.T + translation
    for index, coordinate in zip(carried, rotated):
        self.state.atoms[index].coord = coordinate

matmalib.util.StateBuilder(path, software)

Source code in matmalib/util.py

def __init__(self, path, software):
    self.path = path
    self.software = software

path = path instance-attribute

software = software instance-attribute

element_symbol(A) staticmethod

Return the symbol for an atomic number (1–118).

Source code in matmalib/util.py

@staticmethod
def element_symbol(A):
    """Return the symbol for an atomic number (1–118)."""
    return _ELEMENT_SYMBOLS[str(A)]

gaussian()

Parses log files from Gaussian 16 and appends information inside the log file to the mol object.

Source code in matmalib/util.py

def gaussian(self):

    """
    Parses log files from Gaussian 16 and appends information inside the log file to the mol object.

    """
    flags = {'freq_flag': False, 'nmr_flag': False, 'opt_flag': False, 'jcoup_flag': False, 'normal_mode': False,
             'read_geom': False}
    job_type = None

    # temprorary variables to hold the data
    freq = []
    ints = []
    vibs = []
    geom = []
    atoms = []
    nmr = []
    self.NAtoms = None
    mode_num = -1
    mode_lines = -1
    atom_index = 0
    final_geom = False

    for line in open(self.path, 'r').readlines():

        if not job_type and re.search('^ #', line):

            if "opt" in line:
                if "freq" in line:
                    job_type = 'optfreq'
                else:
                    job_type = 'opt'
            elif "freq" in line:
                if "opt" in line:
                    job_type = 'optfreq'

                else:
                    job_type = 'freq'
                    flags["freq_flag"] = True

            elif "nmr" in line:
                job_type = 'nmr'
            else:
                job_type = 'sp'

        if self.NAtoms is None and re.search('^ NAtoms=', line):
            self.NAtoms = int(line.split()[1])

        if job_type == 'optfreq' or job_type == "freq":

            if flags['freq_flag'] == False and re.search('Normal termination', line): flags['freq_flag'] = True
            # We skip the opt part of optfreq job, all info is in the freq part

            if flags['freq_flag'] == True:

                if re.search('SCF Done', line):
                    self.E = float(line.split()[4])
                elif re.search('Sum of electronic and zero-point Energies', line):
                    self.Ezpe = float(line.split()[6])
                elif re.search('Sum of electronic and thermal Enthalpies', line):
                    self.H = float(line.split()[6])
                elif re.search('Sum of electronic and thermal Free Energies', line):
                    self.F = float(line.split()[7])

                elif re.search('Standard orientation:', line):
                    final_geom = True

                elif re.search('Coordinates', line) and len(geom) == 0:
                    flags['read_geom'] = True

                elif flags['read_geom'] == True and re.search(
                        r'^\s*\d+\s+\d+\s+\d+\s+[+-]?\d+\.\d{6}\s+[+-]?\d+\.\d{6}\s+[+-]?\d+\.\d{6}\s*$',
                        line) and final_geom:

                    if len(line.split()) != 6:
                        continue

                    geom.append([float(x) for x in line.split()[3:6]])
                    coord = [float(x) for x in line.split()[3:6]]
                    type = StateBuilder.element_symbol(line.split()[1])
                    atom = matmacore.Atom(coord, type)

                    atoms.append(atom)
                    atom_index = atom_index + 1

                    if int(line.split()[0]) == self.NAtoms:
                        flags['read_geom'] = False

                elif re.search(' Deg. of freedom', line):
                    self.NVibs = int(line.split()[3])

                elif re.search('^ Frequencies', line):
                    freq_line = line.strip()
                    for f in freq_line.split()[2:5]: freq.append(float(f))
                    flags['normal_mode'] = False

                elif re.search('^ IR Inten', line):
                    ir_line = line.strip()
                    for i in ir_line.split()[3:6]: ints.append(float(i))

                if 'Atom  AN' in line:
                    flags['normal_mode'] = True  # locating normal modes of a frequency
                    mode_1 = []
                    mode_2 = []
                    mode_3 = []
                    # continue


                elif re.search(r'^\s*\d+\s+\d+(?:\s+[+-]?\d+\.\d+){3,}\s*$', line) and flags['normal_mode'] == True:

                    if len(line.split()) >= 5:
                        mode_1.append([float(x) for x in line.split()[2:5]])

                    if len(line.split()) >= 8:
                        mode_2.append([float(x) for x in line.split()[5:8]])

                    if len(line.split()) >= 11:
                        mode_3.append([float(x) for x in line.split()[8:11]])

                    mode_num += 1

                elif mode_num > mode_lines and flags['normal_mode'] == True:

                    flags['normald_mode'] = False
                    for m in [mode_1, mode_2, mode_3]: vibs.append(np.array(m))

        elif job_type == 'opt':

            if re.search('SCF Done', line): E = float(line.split()[4])
            if re.search('Optimization completed.', line):
                self.E = E;
                flags['opt_flag'] = True
            if flags['opt_flag'] == True:
                if re.search('Standard orientation:', line):
                    flags['read_geom'] = True
                    final_geom = True

                elif flags['read_geom'] == True and re.search(
                        r'^\s*\d+\s+\d+\s+\d+\s+[+-]?\d+\.\d{6}\s+[+-]?\d+\.\d{6}\s+[+-]?\d+\.\d{6}\s*$',
                        line) and final_geom:

                    if len(line.split()) != 6:
                        continue

                    geom.append([float(x) for x in line.split()[3:6]])
                    coord = [float(x) for x in line.split()[3:6]]
                    type = StateBuilder.element_symbol(line.split()[1])
                    atom = matmacore.Atom(coord, type)

                    atoms.append(atom)
                    atom_index = atom_index + 1

                    if int(line.split()[0]) == self.NAtoms:
                        flags['read_geom'] = False

        elif job_type == 'nmr':

            if re.search('SCF Done', line):
                self.E = float(line.split()[4])
            elif re.search('Coordinates', line) and len(geom) == 0:
                flags['read_geom'] = True

            elif flags['read_geom'] == True and re.search(
                    r'^\s*\d+\s+\d+\s+\d+\s+[+-]?\d+\.\d{6}\s+[+-]?\d+\.\d{6}\s+[+-]?\d+\.\d{6}\s*$',
                    line):

                if len(line.split()) != 6:
                    continue

                geom.append([float(x) for x in line.split()[3:6]])
                coord = [float(x) for x in line.split()[3:6]]
                type = StateBuilder.element_symbol(line.split()[1])
                atom = matmacore.Atom(coord, type)

                atoms.append(atom)
                atom_index = atom_index + 1

                if int(line.split()[0]) == self.NAtoms:
                    flags['read_geom'] = False

            elif re.search('Total nuclear spin-spin coupling J', line):
                spin = [[] for i in range(self.NAtoms)]
                flags['jcoup_flag'] = True

            elif flags['jcoup_flag'] == True and re.search(r'-?\d\.\d+[Dd][+\-]\d\d?', line):
                for x in line.split()[1:]:
                    spin[int(line.split()[0]) - 1].append(float(x.replace('D', 'E')))

            elif flags['jcoup_flag'] == True and re.search('End of Minotr F.D. properties file', line):
                flags['jcoup_flag'] = False

        elif job_type == 'sp':

            if re.search('SCF Done', line):
                self.E = float(line.split()[4])
            elif re.search('Standard orientation:', line):
                flags['read_geom'] = True

            elif flags['read_geom'] == True and re.search(
                    r'^\s*\d+\s+\d+\s+\d+\s+[+-]?\d+\.\d{6}\s+[+-]?\d+\.\d{6}\s+[+-]?\d+\.\d{6}\s*$',
                    line):

                if len(line.split()) != 6:
                    continue

                geom.append([float(x) for x in line.split()[3:6]])
                coord = [float(x) for x in line.split()[3:6]]
                type = StateBuilder.element_symbol(line.split()[1])
                atom = matmacore.Atom(coord, type)

                atoms.append(atom)
                atom_index = atom_index + 1

                if int(line.split()[0]) == self.NAtoms:
                    flags['read_geom'] = False

    # postprocessing:
    if job_type == 'freq' or job_type == 'optfreq':
        self.Freq = np.array(freq)
        self.Ints = np.array(ints)
        self.Vibs = np.zeros((self.NVibs, self.NAtoms, 3))

        for i in range(self.NVibs): self.Vibs[i, :, :] = vibs[i]

    if job_type == 'nmr':
        for at in spin:
            while len(at) < self.NAtoms: at.append(0)
        self.nmr = np.tril(spin)

    for i, atom in enumerate(atoms):
        atoms[i].id = i

    self.atoms = atoms

cp2k_geom()

Parses information from CP2K. Right now only works for AIMD Trajectories, plan to implement opt compatability later.

Source code in matmalib/util.py

def cp2k_geom(self):

    """ Parses information from CP2K. Right now only works for AIMD Trajectories, plan to implement opt compatability later.
    """

    with open(self.path, 'r') as geom:
        first_line = geom.readline()
        NAtoms = int(first_line.split()[0])

        raw_coords = geom.readlines()[-NAtoms:]

        coords = np.array([line.split()[1:] for line in raw_coords], dtype=float)
        types = np.array([line.split()[0] for line in raw_coords])

        geom.close()

    index = 0
    atoms = []
    for coord, type in zip(coords, types):
        atom = matmacore.Atom(coord, type, index)
        atoms.append(atom)
        index += 1

    self.NAtoms = NAtoms
    self.atoms = atoms

orca()

Parses information from ORCA input.log file

Source code in matmalib/util.py

def orca(self):

    """ Parses information from ORCA input.log file
    """

    with open(self.path, "r") as f:
        for line in f:
            if "Geometry Optimization Run" in line:
                output_file = "input.orca.xyz"

                with open("/".join([self.path, output_file]), "r") as f2:
                    first_line = f2.readline()
                    Natoms = int(first_line.split()[0])

                    second_line = f2.readlines()[0:1]
                    energy = np.array([i.split()[-1] for i in second_line], dtype=float)

                    raw_coords = f2.readlines()[-Natoms:]

                    coords = np.array([line.split()[1:] for line in raw_coords], dtype=float)
                    types = np.array([line.split()[0] for line in raw_coords])

                self.Natoms = Natoms
                self.xyz = coords
                self.type_list = atoms

                index = 0
                atoms = []
                for coord, type in zip(coords, types):
                    atom = matmacore.Atom(coord, type, index)
                    atoms.append(atom)
                    index += 1

                self.atoms = atoms

                self.E = energy

            elif "Single Point Calculation" and "FINAL SINGLE POINT ENERGY" in line:
                # print(line.split()[-1])
                job_type = "sp"
                energy = np.array(line.split()[-1], dtype=float)
                self.E = energy
                return self

            else:
                print("Unable to read ORCA log file")

build()

Source code in matmalib/util.py

def build(self):

    from .core import State

    if self.software == 'cp2k':
        self.cp2k_geom()
    elif self.software == 'gaussian':
        self.gaussian()
    elif self.software == 'orca':
        self.orca()
    else:
        raise ValueError(f"Unsupported software: {self.software}")
    state = State(**self.__dict__)
    state.connectivity_matrix()
    state.system_split()
    return state

matmalib.util.DataBuilder(path, software, **kwargs)

Source code in matmalib/util.py

def __init__(self, path, software, **kwargs):

    self.path = path
    self.software = software

    for k, v in kwargs.items():
        setattr(self, k, v)

path = path instance-attribute

software = software instance-attribute

cp2k_md(colvar=None, timestep=None)

Source code in matmalib/util.py

def cp2k_md(self, colvar=None, timestep=None):

    from .core import State
    from .core import Atom

    self.time_unit = 'fs'  # Default unit in CP2K

    if getattr(self, 'observable', None) == 'cremer_pople':
        from .glycan import _xyz_frames, cremer_pople_trajectory
        timestep = getattr(self, 'timestep', None) if timestep is None else timestep
        if timestep is None:
            raise ValueError('Specify timestep for reading CP2K trajectories.')
        result = cremer_pople_trajectory(
            _xyz_frames(self.path), getattr(self, 'rings', None), timestep=timestep,
            degrees=getattr(self, 'degrees', True), path=self.path, software='cp2k',
            conn_mat=getattr(self, 'conn_mat', None), distXX=getattr(self, 'distXX', 1.65),
            distXH=getattr(self, 'distXH', 1.15), classify=getattr(self, 'classify', False),
        )
        self.__dict__.update(result.__dict__)
        return

    if re.search('.xyz', self.path):
        # Reads colvars straight from the xyz file. Requires specification of timestep and colvar coordinate.

        colvar = getattr(self, 'colvar', None) if colvar is None else colvar
        timestep = getattr(self, 'timestep', None) if timestep is None else timestep

        if colvar == None:
            raise ValueError('Error: Specify the atoms for your collective variables.')

        if timestep == None:
            raise ValueError('Error: Specify timestep for reading CP2K trajectories.')

        with open(self.path, 'r') as geom:
            line_1 = geom.readline()
            self.Natoms = int(line_1.split()[0])

            line_cnt = self.Natoms + 2
            line_num = 1
            conf_array = []
            conf_raw_coord = []

            for line in geom:

                for col in colvar:
                    if line_num % line_cnt == (col + 2) % line_cnt:
                        conf_raw_coord.append(line.split())

                if len(conf_raw_coord) == len(colvar):
                    conf_array.append(conf_raw_coord)
                    conf_raw_coord = []

                line_num += 1

        conf_array = np.array(conf_array)
        mol_space = []
        colvar_data = []

        for conf in conf_array:
            index = 0
            state = State()
            atoms = []

            for row in conf:
                type = row[0]
                xyz = [float(i) for i in row[1:]]

                atom = Atom(coord=xyz, type=type, index=index)
                atoms.append(atom)

                index += 1

            state.atoms = atoms
            mol_space.append(state)

            col_data = StateEditor.measure(atoms)
            colvar_data.append(col_data)

        time = [i / (1 / timestep) for i in range(len(colvar_data))]

        data = np.array(list(zip(time, colvar_data)))
        self.data = data

    elif re.search('.metadynLog', self.path):
        # Reads colvars from the metadynLog file.

        data = np.loadtxt(self.path)

        self.time_unit = 'fs'
        self.data = data

gromacs()

Parses information from gromacs *.xvg file

Source code in matmalib/util.py

def gromacs(self):
    """ Parses information from gromacs *.xvg file
    """

    # I don't use Gromacs enough to know how the output looks like. Right now, I am just reading the .xvg file:

    if re.search('.xvg', self.path):

        with open(self.path, 'r') as f:
            i = 0
            for line in f.readlines():
                if re.search('#', line) or re.search('@', line):
                    i = i + 1

                if re.search('xaxis', line):
                    time_unit = line.split()[-1].strip('()"')
                    self.time_unit = time_unit
            f.close()

        data = np.loadtxt(self.path, skiprows=i)

        self.data = data

    if re.search('.xpm', self.path):

        f = open(self.path, 'r')
        pattern = None;
        matrix_lett = [];
        matrix_dict = {}
        grid = None
        for line in f.readlines():

            if re.search(r'(\d\s\d\b)', line) and grid is None:
                grid = int(line.split()[1])
                pattern = grid * '[A-Z]'

            if pattern and re.search(pattern, line):
                replL1 = line.replace('",', '')
                replL2 = line.replace('"', '')
                matrix_lett.append('\n'.join(replL2.split()))

            if re.search('((".*")[0-9])', line):
                replL1 = line.replace('"', ' ')
                letters, energy = values = str(replL1.split()[0]), float(replL1.split()[4])
                matrix_dict.update(dict.fromkeys(letters, energy))

        matrix = np.zeros([grid, grid])
        for i in range(grid):
            for j in range(grid):
                matrix[grid - i - 1, j] = matrix_dict[matrix_lett[i][j]]

        self.data = matrix

csv(skiprows=0, dtype=float, delimiter: str = ',')

Parse the CSV file at this builder's path into its data attribute.

Parameters:

Name	Type	Description	Default
skiprows		(int) Number of rows to skip	0
dtype		(type) Data type	float
delimiter	str	(str) Delimiter	','

Source code in matmalib/util.py

def csv(self, skiprows=0, dtype=float, delimiter: str = ','):
    """
    Parse the CSV file at this builder's path into its data attribute.
    :param skiprows: (int) Number of rows to skip
    :param dtype: (type) Data type
    :param delimiter: (str) Delimiter
    """

    skiprows = getattr(self, 'skiprows', skiprows)
    dtype = getattr(self, 'dtype', dtype)
    delimiter = getattr(self, 'delimiter', delimiter)

    self.data = np.loadtxt(self.path, skiprows=skiprows, dtype=dtype, delimiter=delimiter)

build()

Source code in matmalib/util.py

def build(self):

    from .core import Data

    if self.software == 'cp2k':
        self.cp2k_md()
    elif self.software == 'gromacs':
        self.gromacs()
    elif self.software == 'csv':
        self.csv()
    else:
        raise ValueError(f"Unsupported software: {self.software}")

    return Data(**self.__dict__)

matmalib.calc.InputBuilder(path: str = None, theory: dict = None, software: str = 'gaussian', space=None)

Source code in matmalib/calc.py

def __init__(self, path:str = None, theory:dict = None, software:str = 'gaussian', space = None):
    self.path = path
    self.theory = theory
    self.software = software
    self.space = space

path = path instance-attribute

theory = theory instance-attribute

software = software instance-attribute

space = space instance-attribute

element_number(A) staticmethod

Return the atomic number for an element symbol from H through Og.

Source code in matmalib/calc.py

@staticmethod
def element_number(A):
    """Return the atomic number for an element symbol from H through Og."""
    return _ELEMENT_NUMBERS[A]

gaussian()

Source code in matmalib/calc.py

def gaussian(self):
    space = self.space
    theory = self.theory
    names = self.names

    if names is not None:
        for i in range(space.states):
            space.states[i].id = names[i]

    for state in space.states:
        if theory['disp'] in [True, 'EmpiricalDispersion=GD3BJ', 'GD3BJ']:
            theory['disp'] = 'EmpiricalDispersion=GD3BJ'
        elif theory['disp'] in ['EmpiricalDispersion=GD3', 'GD3']:
            theory['disp'] = 'EmpiricalDispersion=GD3'
        else:
            theory['disp'] = ' '


        input_dir = os.path.join(self.path,state.id)
        if not os.path.isdir(input_dir):
            os.mkdir(input_dir)
        input_file = os.path.join(input_dir,'input.com')

        f = open(input_file, 'w')
        f.write('%nproc=' + str(theory['nprocs']) + '\n')
        f.write('%mem=' + theory['mem'] + '\n')
        if 'other_options' in theory.keys():
            f.write(' '.join(['#P', theory['method'], theory['basis_set'], theory['jobtype'], theory['other_options'],
                          theory['disp'], '\n']))
        else:
            f.write(
                ' '.join(['#P', theory['method'], theory['basis_set'], theory['jobtype'],
                          theory['disp'], '\n']))
        f.write('\n')
        f.write(state.id + '\n')
        f.write('\n ')
        f.write(str(theory['charge']) + ' ' + str(theory['multiplicity']) + '\n')
        for at, xyz in zip(state.atoms, state.xyz):
            line = '{0:5s} {1:10.3f} {2:10.3f} {3:10.3f}\n'.format(at.type, xyz[0], xyz[1], xyz[2])
            f.write(line)
        f.write(' ')

        if 'extra' in theory.keys():
            if theory['extra'] == None:
                f.close()
            else:
                f.write('\n')
                f.write(theory['extra'] + '\n')
                f.write(' ')
            f.close()
        else:
            f.close()

fhiaims()

Source code in matmalib/calc.py

def fhiaims(self):
    space = self.space
    theory = self.theory

    for state in space:

        if not os.path.isdir(os.path.join(self.path, state.id)):
            os.mkdir(os.path.join(self.path, state.id))

        control_file = os.path.join(self.path, state.id, 'control.in')
        geom_file = os.path.join(self.path, state.id, '/geometry.in')

        c = open(control_file, 'w')
        c.write('xc ' + str(theory['xc']) + '\n')
        c.write(theory['disp'] + '\n')
        c.write('charge ' + str(theory['charge']) + '\n')
        c.write(theory['jobtype'] + '\n')
        c.write(theory['convergence_options'] + '\n')
        c.write('density_update_method ' + theory['density_update_method'] + '\n')
        c.write('check_cpu_consistency ' + theory['check_cpu_consistency'] + '\n')
        diff_atoms = set([atom.type for atom in state.atoms])

        for at in diff_atoms:
            EN = "{0:02d}".format(InputBuilder.element_number(at))
            with open('/exports/apps/fhi-aims.210226/species_defaults/' + theory[
                'basis_set'] + '/' + EN + '_' + at + '_default', 'r') as light:
                for line in light.readlines():
                    c.write(line)
        c.close()

        g = open(geom_file, 'w')
        for n, at, xyz in zip(range(state.NAtoms), state.atoms, state.xyz):
            if n in theory['extra']:
                freeze = 'constrain_relaxation .true.'
            else:
                freeze = ''
            line = 'atom      {0:10.3f} {1:10.3f} {2:10.3f} {3:3s}{4:s}\n'.format(xyz[0], xyz[1], xyz[2], at,
                                                                                  freeze)
            g.write(line)
        g.close()

xyz()

Source code in matmalib/calc.py

def xyz(self):
    space = self.space
    for state in space:

        input_dir = os.path.join(self.path,state.id)
        if not os.path.isdir(input_dir):
            os.mkdir(input_dir)
        input_file = os.path.join(input_dir,'geometry.xyz')

        f = open(input_file, 'w')
        f.write('{0:3d}\n'.format(state.NAtoms))
        f.write(f'{state.id}\n')
        for at, xyz in zip(state.atoms, state.xyz):
            line = '{0:5s} {1:10.3f} {2:10.3f} {3:10.3f}\n'.format(at.type, xyz[0], xyz[1], xyz[2])
            f.write(line)
        f.close()

build(space, theory, software, path, names: list = None)

Source code in matmalib/calc.py

def build(self, space, theory, software, path, names:list = None):
    self.space = space
    self.theory = theory
    self.software = software
    self.path = path
    self.names = names

    if software == 'gaussian':
        self.gaussian()
    elif software == 'fhiaims':
        self.fhiaims()
    elif software == 'xyz':
        self.xyz()

Carbohydrates and trajectories

Tutorials: Glycans, Ring puckering, and Trajectories. Glycan inherits the State methods above.

matmalib.glycan.Glycan(atoms, *, rings=None, monomers=None, distXX=1.65, distXH=1.15, **kwargs)

Bases: State

A carbohydrate State whose monomers are Fragment objects.

rings and monomers are parallel sequences of zero-based positions in atoms (not Atom.id values). Each ring must be in cyclic order; each monomer lists all atoms assigned to that sugar. Assign bridging atoms to only one monomer. Unassigned atoms remain in the State and are available through unassigned_atoms (for example solvent or counterions).

With no assignments, detect sugars using CarPpy's connectivity-based approach. Supplied conn_mat or graph uses atom-list positions; otherwise bonds are inferred with distXX=1.65 and distXH=1.15 (angstrom defaults). Automatic rings start at O, followed by the anomeric carbon and the remaining ring carbons. Ketoses start at C2, aldoses at C1. Supply both rings and monomers to bypass detection for unusual chemistry.

Source code in matmalib/glycan.py

def __init__(self, atoms, *, rings=None, monomers=None, distXX=1.65, distXH=1.15, **kwargs):
    super().__init__(**kwargs)
    self.atoms = list(atoms)
    self.NAtoms = len(self.atoms)
    assignments = None
    links = []
    if monomers is not None and rings is None:
        raise ValueError("Supply ordered rings when specifying monomers explicitly.")
    if rings is not None:
        rings = _rings(rings, self.NAtoms)
    if rings is None or monomers is None:
        rings, detected_members, assignments, graph, links = assign_sugars(
            self.atoms, rings=rings, conn_mat=getattr(self, "conn_mat", None),
            graph=getattr(self, "graph", None), distXX=distXX, distXH=distXH,
        )
        if monomers is None:
            monomers = detected_members
        self.graph = graph
        self.Nmols = nx.number_connected_components(graph)
    self.rings = _rings(rings, self.NAtoms)
    memberships = tuple(
        _indices(group, name="Monomer", size=self.NAtoms) for group in monomers
    )
    if len(memberships) != len(self.rings):
        raise ValueError("Supply one monomer atom assignment per ring.")
    assigned = set()
    fragments = []
    for index, (ring, members) in enumerate(zip(self.rings, memberships)):
        if not set(ring).issubset(members):
            raise ValueError(f"Monomer {index} does not contain its ring atoms.")
        if assigned.intersection(members):
            raise ValueError("Monomers must not share atoms.")
        assigned.update(members)
        fragment = Fragment()
        fragment.id = index
        fragment.atom_indices = members
        fragment.ring_indices = ring
        fragment.atoms = [self.atoms[i] for i in members]
        fragment.ring_atoms = [self.atoms[i] for i in ring]
        fragment.ring_type = "pyranose" if len(ring) == 6 else "furanose"
        fragment.atom_assignment = dict(assignments[index]) if assignments is not None else {}
        fragment.named_atoms = {name: self.atoms[i] for name, i in fragment.atom_assignment.items()}
        if hasattr(self, "graph"):
            fragment.graph = self.graph.subgraph(members).copy()
        fragments.append(fragment)
    self.monomers = fragments
    self.fragments = fragments
    self.unassigned_atoms = [atom for i, atom in enumerate(self.atoms) if i not in assigned]
    self.glycosidic_bonds = links

atoms = list(atoms) instance-attribute

NAtoms = len(self.atoms) instance-attribute

graph = graph instance-attribute

Nmols = nx.number_connected_components(graph) instance-attribute

rings = _rings(rings, self.NAtoms) instance-attribute

monomers = fragments instance-attribute

fragments = fragments instance-attribute

unassigned_atoms = [atom for i, atom in enumerate(self.atoms) if i not in assigned] instance-attribute

glycosidic_bonds = links instance-attribute

from_state(state, *, rings=None, monomers=None, **kwargs) classmethod

Wrap a State, retaining its metadata and sharing its Atom objects.

Source code in matmalib/glycan.py

@classmethod
def from_state(cls, state, *, rings=None, monomers=None, **kwargs):
    """Wrap a State, retaining its metadata and sharing its Atom objects."""
    metadata = dict(state.__dict__)
    metadata.pop("atoms", None)
    metadata.pop("rings", None)
    metadata.pop("monomers", None)
    metadata.update(kwargs)
    return cls(state.atoms, rings=rings, monomers=monomers, **metadata)

from_file(path, software='xyz', *, rings=None, monomers=None, **kwargs) classmethod

Read a system and detect sugars; XYZ/CP2K uses the last saved frame.

Source code in matmalib/glycan.py

@classmethod
def from_file(cls, path, software="xyz", *, rings=None, monomers=None, **kwargs):
    """Read a system and detect sugars; XYZ/CP2K uses the last saved frame."""
    if software in ("xyz", "cp2k"):
        state = None
        for state in _xyz_frames(path):
            pass
        if state is None:
            raise ValueError("The trajectory contains no frames.")
        state.path = str(path)
        state.software = software
        state.id = str(path)
    else:
        state = State.from_file(str(path), software)
    return cls.from_state(state, rings=rings, monomers=monomers, **kwargs)

cremer_pople(*, degrees=True)

Return one puckering dictionary per monomer for the current geometry.

Source code in matmalib/glycan.py

def cremer_pople(self, *, degrees=True):
    """Return one puckering dictionary per monomer for the current geometry."""
    xyz = self.xyz
    return [cremer_pople(xyz[list(ring)], degrees=degrees) for ring in self.rings]

identify_ring_puckers(**kwargs)

Return CP coordinates, canonical pucker labels, and distances per ring.

Six-membered labels follow CarPpy's positional convention: 1–5 are the ring carbons after O, and 6 is O. For ketoses these positions are not the chemical carbon numbers. Five-membered labels use 0–4 in ring order. See classify_ring_pucker for tolerances and conventions.

Source code in matmalib/glycan.py

def identify_ring_puckers(self, **kwargs):
    """Return CP coordinates, canonical pucker labels, and distances per ring.

    Six-membered labels follow CarPpy's positional convention: 1–5 are
    the ring carbons after O, and 6 is O. For ketoses these positions are
    not the chemical carbon numbers. Five-membered labels use 0–4 in
    ring order. See classify_ring_pucker for tolerances and conventions.
    """
    xyz = self.xyz
    return [classify_ring_pucker(xyz[list(ring)], **kwargs) for ring in self.rings]

identify_sugars()

Return base-residue names, D/L configuration and alpha/beta anomer.

Common aldohexopyranoses and their C2 amino/N-acetyl, C6 deoxy and uronic-acid variants are supported. Other chemistry remains Unknown.

Source code in matmalib/glycan.py

def identify_sugars(self):
    """Return base-residue names, D/L configuration and alpha/beta anomer.

    Common aldohexopyranoses and their C2 amino/N-acetyl, C6 deoxy and
    uronic-acid variants are supported. Other chemistry remains Unknown.
    """
    from ._glycan_analysis import identify_sugars
    return identify_sugars(self)

glycosidic_torsions(*, degrees=True)

Return phi/psi/omega and their atom quadruplets for each linkage.

Source code in matmalib/glycan.py

def glycosidic_torsions(self, *, degrees=True):
    """Return phi/psi/omega and their atom quadruplets for each linkage."""
    from ._glycan_analysis import glycosidic_torsions
    return glycosidic_torsions(self, degrees=degrees)

backbone_torsions(*, degrees=True)

Return hydroxymethyl omega and its atom quadruplet per monomer.

Source code in matmalib/glycan.py

def backbone_torsions(self, *, degrees=True):
    """Return hydroxymethyl omega and its atom quadruplet per monomer."""
    from ._glycan_analysis import backbone_torsions
    return backbone_torsions(self, degrees=degrees)

snfg(**kwargs)

Return an SNFG (figure, axes); accepts draw_snfg plotting options.

Source code in matmalib/glycan.py

def snfg(self, **kwargs):
    """Return an SNFG (figure, axes); accepts draw_snfg plotting options."""
    from .snfg import draw_snfg
    return draw_snfg(self, **kwargs)

cremer_pople_trajectory(frames, **kwargs)

Analyze frames using this Glycan's fixed atom assignments.

Accepts a Space, an iterable of States, or coordinate arrays with shape (n_frames, n_atoms, 3). Every frame must use this State's atom order. See :func:cremer_pople_trajectory for output and time conventions.

Source code in matmalib/glycan.py

def cremer_pople_trajectory(self, frames, **kwargs):
    """Analyze frames using this Glycan's fixed atom assignments.

    Accepts a Space, an iterable of States, or coordinate arrays with shape
    (n_frames, n_atoms, 3). Every frame must use this State's atom order.
    See :func:`cremer_pople_trajectory` for output and time conventions.
    """
    elements = tuple(atom.type for atom in self.atoms)

    def checked_frames():
        for index, frame in enumerate(getattr(frames, "states", frames)):
            if isinstance(frame, State):
                if tuple(atom.type for atom in frame.atoms) != elements:
                    raise ValueError(f"Frame {index} has different elements or atom order.")
                xyz = frame.xyz
            else:
                xyz = np.asarray(frame, dtype=float)
            if xyz.shape != (self.NAtoms, 3):
                raise ValueError(f"Frame {index} must have shape ({self.NAtoms}, 3).")
            yield xyz

    result = cremer_pople_trajectory(checked_frames(), self.rings, **kwargs)
    result.glycan = self
    return result

matmalib.glycan.cremer_pople_trajectory(frames, rings=None, *, timestep=1.0, times=None, time_unit='fs', degrees=True, path='', software=None, elements=None, conn_mat=None, distXX=1.65, distXH=1.15, classify=False)

Analyze ordered rings in a trajectory and return a plotting-ready Data.

frames is a Space, an iterable of States, or an iterable of full coordinate arrays. Without rings, detect sugars on the first State (or supply elements for coordinate arrays). conn_mat can provide reference connectivity; otherwise CarPpy distance cutoffs infer bonds. rings lists zero-based atom positions in cyclic order, held fixed for every frame. Frames must have identical atom order and whole, unwrapped rings; no periodic-boundary correction is performed.

data[:, 0] contains explicit times or arange(n_frames)*timestep. columns names every column: time, then monomer_0.Q, monomer_0.theta, monomer_0.phi, monomer_0.q2, monomer_0.q3, and so on. Five-membered rings omit theta and q3. Angles default to degrees; lengths retain input units. Frames are consumed once, without storing the full coordinate trajectory. With classify=True, also attach pucker_labels, pucker_families, and pucker_distances arrays of shape (n_frames, n_rings), following classify_ring_pucker's positional naming and degree-distance conventions. The numeric data columns are unchanged.

Source code in matmalib/glycan.py

def cremer_pople_trajectory(
    frames, rings=None, *, timestep=1.0, times=None, time_unit="fs", degrees=True,
    path="", software=None, elements=None, conn_mat=None, distXX=1.65, distXH=1.15,
    classify=False,
):
    """Analyze ordered rings in a trajectory and return a plotting-ready Data.

    ``frames`` is a Space, an iterable of States, or an iterable of full
    coordinate arrays. Without ``rings``, detect sugars on the first State
    (or supply ``elements`` for coordinate arrays). ``conn_mat`` can provide
    reference connectivity; otherwise CarPpy distance cutoffs infer bonds.
    ``rings`` lists zero-based atom positions in cyclic order, held fixed for
    every frame. Frames must have identical atom order
    and whole, unwrapped rings; no periodic-boundary correction is performed.

    ``data[:, 0]`` contains explicit ``times`` or ``arange(n_frames)*timestep``.
    ``columns`` names every column: time, then monomer_0.Q, monomer_0.theta,
    monomer_0.phi, monomer_0.q2, monomer_0.q3, and so on. Five-membered rings
    omit theta and q3. Angles default to degrees; lengths retain input units.
    Frames are consumed once, without storing the full coordinate trajectory.
    With ``classify=True``, also attach ``pucker_labels``, ``pucker_families``,
    and ``pucker_distances`` arrays of shape (n_frames, n_rings), following
    classify_ring_pucker's positional naming and degree-distance conventions.
    The numeric data columns are unchanged.
    """
    frames = iter(getattr(frames, "states", frames))
    glycan = None
    if rings is None:
        first = next(frames, None)
        if first is None:
            raise ValueError("The trajectory contains no frames.")
        if isinstance(first, State):
            reference = first
        else:
            xyz = np.asarray(first, dtype=float)
            if elements is None:
                raise ValueError("Supply elements or explicit rings for coordinate-array trajectories.")
            elements = tuple(elements)
            if xyz.shape != (len(elements), 3):
                raise ValueError("elements must match the coordinate array's atom count.")
            reference = State(atoms=[Atom(coord=coord, type=element, index=i)
                                    for i, (coord, element) in enumerate(zip(xyz, elements))])
        if isinstance(reference, Glycan) and conn_mat is None:
            glycan = reference
        else:
            options = dict(distXX=distXX, distXH=distXH)
            if conn_mat is not None:
                options["conn_mat"] = conn_mat
            glycan = Glycan.from_state(reference, **options)
        rings = glycan.rings
        frames = chain((first,), frames)
    rings = _rings(rings)
    if times is None and (not np.isfinite(timestep) or timestep <= 0):
        raise ValueError("timestep must be finite and positive.")
    names = [list(("Q", "theta", "phi", "q2", "q3") if len(ring) == 6
                  else ("Q", "phi", "q2")) for ring in rings]
    columns = ("time",) + tuple(
        f"monomer_{i}.{name}" for i, fields in enumerate(names) for name in fields
    )
    rows = []
    pucker_labels, pucker_families, pucker_distances = [], [], []
    shape = None
    elements = None
    for index, frame in enumerate(frames):
        if isinstance(frame, State):
            current_elements = tuple(atom.type for atom in frame.atoms)
            if elements is not None and elements != current_elements:
                raise ValueError(f"Frame {index} has different elements or atom order.")
            elements = current_elements
            xyz = frame.xyz
        else:
            xyz = np.asarray(frame, dtype=float)
        if xyz.ndim != 2 or xyz.shape[1] != 3 or not np.all(np.isfinite(xyz)):
            raise ValueError(f"Frame {index} must contain finite (n_atoms, 3) coordinates.")
        if shape is None:
            shape = xyz.shape
            _rings(rings, size=len(xyz))
        elif xyz.shape != shape:
            raise ValueError(f"Frame {index} has a different number of atoms.")
        row = []
        labels, families, distances = [], [], []
        for ring, fields in zip(rings, names):
            try:
                analyze = classify_ring_pucker if classify else cremer_pople
                result = analyze(xyz[list(ring)], degrees=degrees)
            except ValueError as exc:
                raise ValueError(f"Frame {index}, ring {ring}: {exc}") from exc
            row.extend(result[name] for name in fields)
            if classify:
                labels.append(result["label"])
                families.append(result["family"])
                distances.append(result["distance"])
        rows.append(row)
        if classify:
            pucker_labels.append(labels)
            pucker_families.append(families)
            pucker_distances.append(distances)
    if not rows:
        raise ValueError("The trajectory contains no frames.")
    if times is None:
        times = np.arange(len(rows)) * timestep
    times = np.asarray(times, dtype=float)
    if times.shape != (len(rows),) or not np.all(np.isfinite(times)):
        raise ValueError("times must contain one finite value per frame.")
    data = Data(
        data=np.column_stack((times, rows)), columns=columns, rings=rings,
        time_unit=time_unit, angle_unit="degrees" if degrees else "radians",
        path=str(path), software=software, glycan=glycan,
    )
    if classify:
        data.pucker_labels = np.asarray(pucker_labels)
        data.pucker_families = np.asarray(pucker_families)
        data.pucker_distances = np.asarray(pucker_distances)
    return data

matmalib.puckering.cremer_pople(coordinates, *, degrees=True)

Calculate ring puckering for one geometry or a batch of geometries.

coordinates has shape (..., ring_size, 3) with atoms in cyclic order. All rings in a batch must have the same size (five or six). Coordinates must describe whole, unwrapped rings. Amplitudes retain the input length unit; angles are degrees by default, or radians when degrees=False.

The returned dictionary contains Q, q2, and phi. Six-membered rings also have the signed amplitude q3 and polar angle theta. phi lies in [0, 360) and theta in [0, 180] (or [0, 2*pi) and [0, pi]). Undefined angles are NaN: phi for planar rings and pure chairs, and theta for planar rings. A degenerate mean plane raises ValueError.

Uses the original Cremer–Pople convention (JACS 97, 1354–1358, 1975, https://doi.org/10.1021/ja00839a011): the mean-plane normal is R' x R'', with sine-weighted R' and cosine-weighted R''. The phase is defined by q2cos(phi) = sqrt(2/N)sum(zcos(4pij/N)) and q2sin(phi) = -sqrt(2/N)sum(zsin(4pij/N)), for j starting at zero. Changing the starting atom or traversal direction changes the angles.

Source code in matmalib/puckering.py

def cremer_pople(coordinates, *, degrees=True):
    """Calculate ring puckering for one geometry or a batch of geometries.

    ``coordinates`` has shape ``(..., ring_size, 3)`` with atoms in cyclic
    order. All rings in a batch must have the same size (five or six).
    Coordinates must describe whole, unwrapped rings. Amplitudes retain the
    input length unit; angles are degrees by default, or radians when
    ``degrees=False``.

    The returned dictionary contains ``Q``, ``q2``, and ``phi``. Six-membered
    rings also have the signed amplitude ``q3`` and polar angle ``theta``.
    ``phi`` lies in [0, 360) and ``theta`` in [0, 180] (or [0, 2*pi) and
    [0, pi]). Undefined angles are NaN: phi for planar rings and pure chairs,
    and theta for planar rings. A degenerate mean plane raises ValueError.

    Uses the original Cremer–Pople convention (JACS 97, 1354–1358, 1975,
    https://doi.org/10.1021/ja00839a011): the mean-plane normal is R' x R'',
    with sine-weighted R' and cosine-weighted R''. The phase is defined by
    q2*cos(phi) = sqrt(2/N)*sum(z*cos(4*pi*j/N)) and
    q2*sin(phi) = -sqrt(2/N)*sum(z*sin(4*pi*j/N)), for j starting at zero.
    Changing the starting atom or traversal direction changes the angles.
    """
    xyz = np.asarray(coordinates, dtype=float)
    if xyz.ndim < 2 or xyz.shape[-1] != 3 or xyz.shape[-2] not in (5, 6):
        raise ValueError("Expected coordinates with shape (..., 5 or 6, 3).")
    if not np.all(np.isfinite(xyz)):
        raise ValueError("Ring coordinates must be finite.")

    size = xyz.shape[-2]
    centered = xyz - xyz.mean(axis=-2, keepdims=True)
    phase = 2 * np.pi * np.arange(size) / size
    sine_vector = np.einsum("...ij,i->...j", centered, np.sin(phase))
    cosine_vector = np.einsum("...ij,i->...j", centered, np.cos(phase))
    normal = np.cross(sine_vector, cosine_vector)
    normal_length = np.linalg.norm(normal, axis=-1)
    scale = np.sum(centered**2, axis=(-2, -1))
    tolerance = 64 * np.finfo(float).eps
    if np.any(normal_length <= tolerance * scale):
        raise ValueError("Ring geometry has a degenerate Cremer–Pople mean plane.")
    normal = normal / normal_length[..., None]
    heights = np.einsum("...ij,...j->...i", centered, normal)

    cosine_mode = np.sqrt(2 / size) * np.einsum(
        "...i,i->...", heights, np.cos(2 * phase)
    )
    sine_mode = -np.sqrt(2 / size) * np.einsum(
        "...i,i->...", heights, np.sin(2 * phase)
    )
    q2 = np.hypot(cosine_mode, sine_mode)
    amplitude_tolerance = tolerance * np.sqrt(scale)
    phi = np.mod(np.arctan2(sine_mode, cosine_mode), 2 * np.pi)
    # Modulo can round a small negative angle up to exactly 2*pi.
    phi = np.where(phi >= 2 * np.pi, 0.0, phi)
    phi = np.where(q2 <= amplitude_tolerance, np.nan, phi)
    result = {"Q": q2, "q2": q2, "phi": np.degrees(phi) if degrees else phi}

    if size == 6:
        q3 = np.einsum("...i,i->...", heights, (-1.0) ** np.arange(size)) / np.sqrt(size)
        amplitude = np.hypot(q2, q3)
        theta = np.where(
            amplitude <= amplitude_tolerance, np.nan, np.arctan2(q2, q3)
        )
        result.update(Q=amplitude, q3=q3, theta=np.degrees(theta) if degrees else theta)

    return result

matmalib.puckering.ring_dihedrals(coordinates)

CarPpy's three ring-puckering angles, in degrees (not CP theta).

Input shape is (..., 6, 3), in O, C1, C2, C3, C4, C5 order. These are the transformed diagonal dihedrals used by CarPpy's ring_canon routine.

Source code in matmalib/puckering.py

def ring_dihedrals(coordinates):
    """CarPpy's three ring-puckering angles, in degrees (not CP theta).

    Input shape is (..., 6, 3), in O, C1, C2, C3, C4, C5 order. These are
    the transformed diagonal dihedrals used by CarPpy's ring_canon routine.
    """
    xyz = np.asarray(coordinates, dtype=float)
    if xyz.shape[-2:] != (6, 3):
        raise ValueError("CarPpy ring dihedrals require ordered six-membered rings.")
    indices = np.array([[5, 3, 1, 2], [1, 5, 3, 4], [3, 1, 5, 0]])
    angles = dihedral(xyz[..., indices, :])
    return np.sign(angles) * (180 - np.abs(angles))

matmalib.puckering.classify_ring_pucker(coordinates, *, degrees=True, planar_tolerance=1e-06, max_distance=None)

Return CP coordinates plus a nearest canonical pucker and its distance.

Six-membered rings use CarPpy's 38 canonical three-angle references; names use C1–C5 and 6 for O, with O first in the input. carppy_theta is the measured three-angle vector and distance is its Euclidean distance from the reference, in degrees. It is distinct from the CP polar theta.

Five-membered rings use 20 ideal envelope/twist phases with indices 0–4 in input order. 0E raises atom 0; 0T1 raises 0 and lowers 1 relative to the CP normal. Distance is the wrapped phase difference in degrees. Angles in the returned CP coordinates follow degrees.

planar_tolerance uses the coordinate length unit. Planar rings are labelled planar. max_distance can reject a remote nearest reference with label unassigned; nearest_label still records the nearest candidate. These assignments describe proximity to ideal conformations, not exact conformational identities.

Source code in matmalib/puckering.py

def classify_ring_pucker(coordinates, *, degrees=True, planar_tolerance=1e-6, max_distance=None):
    """Return CP coordinates plus a nearest canonical pucker and its distance.

    Six-membered rings use CarPpy's 38 canonical three-angle references; names
    use C1–C5 and 6 for O, with O first in the input. ``carppy_theta`` is the
    measured three-angle vector and ``distance`` is its Euclidean distance
    from the reference, in degrees. It is distinct from the CP polar theta.

    Five-membered rings use 20 ideal envelope/twist phases with indices 0–4
    in input order. ``0E`` raises atom 0; ``0T1`` raises 0 and lowers 1 relative
    to the CP normal. Distance is the wrapped phase difference in degrees.
    Angles in the returned CP coordinates follow ``degrees``.

    ``planar_tolerance`` uses the coordinate length unit. Planar rings are
    labelled planar. ``max_distance`` can reject a remote nearest reference
    with label unassigned; nearest_label still records the nearest candidate.
    These assignments describe proximity to ideal conformations, not exact
    conformational identities.
    """
    xyz = np.asarray(coordinates, dtype=float)
    if xyz.shape not in ((5, 3), (6, 3)):
        raise ValueError("Pucker classification requires one five- or six-membered ring.")
    if not np.isfinite(planar_tolerance) or planar_tolerance < 0:
        raise ValueError("planar_tolerance must be finite and nonnegative.")
    if max_distance is not None and (not np.isfinite(max_distance) or max_distance < 0):
        raise ValueError("max_distance must be finite and nonnegative.")
    result = cremer_pople(xyz, degrees=degrees)
    if result["Q"] <= planar_tolerance:
        return dict(result, label="planar", nearest_label=None, family="planar", distance=np.nan)
    if len(xyz) == 6:
        measured = ring_dihedrals(xyz)
        labels = list(CARPPY_PUCKERS)
        distances = np.linalg.norm(np.asarray(list(CARPPY_PUCKERS.values())) - measured, axis=1)
        best = int(np.argmin(distances))
        label = labels[best]
        family = next(name for code, name in (("C", "chair"), ("B", "boat"),
                      ("S", "skew"), ("H", "half-chair"), ("E", "envelope")) if code in label)
        result["carppy_theta"] = measured
        result["method"] = "carppy-three-angle"
    else:
        phase = float(result["phi"] if degrees else np.degrees(result["phi"]))
        angular = 4 * np.pi * np.arange(5) / 5
        labels, phases, families = [], [], []
        for index in range(5):
            for sign in (1, -1):
                height = np.zeros(5)
                height[index] = sign
                labels.append(f"{index}E" if sign == 1 else f"E{index}")
                phases.append(np.degrees(np.arctan2(-height @ np.sin(angular), height @ np.cos(angular))))
                families.append("envelope")
                other = (index + 1) % 5
                height[other] = -sign
                labels.append(f"{index}T{other}" if sign == 1 else f"{other}T{index}")
                phases.append(np.degrees(np.arctan2(-height @ np.sin(angular), height @ np.cos(angular))))
                families.append("twist")
        distances = np.abs((phase - np.asarray(phases) + 180) % 360 - 180)
        best = int(np.argmin(distances))
        label, family = labels[best], families[best]
        result["method"] = "cremer-pople-phase"
    distance = float(distances[best])
    result.update(label=label if max_distance is None or distance <= max_distance else "unassigned",
                  nearest_label=label, family=family, distance=distance)
    return result

matmalib.snfg.draw_snfg(glycan, *, labels=None, ax=None, positions=None, show_linkages=True, show_names=False, symbol_size=0.32, symbol_linewidth=2.0, linkage_linewidth=2.0, bond_length=1.2, diagonal_bond_scale=1.2)

Draw a Glycan as SNFG symbols; return (figure, axes).

Common aldohexopyranoses are named from connectivity and stereochemistry. Unsupported chemistry uses the white Unknown hexagon. labels may be a sequence (one standard SNFG abbreviation per monomer) or an index/name mapping overriding selected names, e.g. {0: 'Fru'}. Explicit unknown abbreviations raise ValueError. Names describe base residues; protecting groups and all modifications are not inventoried by automatic naming.

Reducing ends appear at the right. Automatic linkages extend from each acceptor towards its donor according to the acceptor carbon position: 2 down (270 degrees), 3 lower left (225), 4 left (180), 6 upper left (135), and 8 up (90), independent of branching. bond_length sets the horizontal and vertical center-to-center distance in drawing coordinates (default 1.2). diagonal_bond_scale multiplies this length for diagonal links (default 1.2, or 20 percent longer); use 1 for equal lengths. Alpha linkages are dashed; beta and unknown linkages are solid. symbol_linewidth and linkage_linewidth set stroke widths in points. positions optionally maps every monomer index to (x, y); it is required for non-tree topologies. The lines between explicit positions are drawn straight without snapping or rescaling. Unknown or other acceptor positions and overlapping fixed-length layouts also require explicit positions. Linkage annotations give alpha/beta and donor–acceptor carbon positions. Anomers that cannot be determined are shown as ?. Nondefault D/L configurations and furanose rings are annotated. Figures support normal Matplotlib export.

Source code in matmalib/snfg.py

def draw_snfg(glycan, *, labels=None, ax=None, positions=None, show_linkages=True,
              show_names=False, symbol_size=.32, symbol_linewidth=2.0, linkage_linewidth=2.0,
              bond_length=1.2, diagonal_bond_scale=1.2):
    """Draw a Glycan as SNFG symbols; return ``(figure, axes)``.

    Common aldohexopyranoses are named from connectivity and stereochemistry.
    Unsupported chemistry uses the white Unknown hexagon. ``labels`` may be
    a sequence (one standard SNFG abbreviation per monomer) or an index/name
    mapping overriding selected names, e.g. ``{0: 'Fru'}``. Explicit unknown
    abbreviations raise ValueError. Names describe base residues; protecting
    groups and all modifications are not inventoried by automatic naming.

    Reducing ends appear at the right. Automatic linkages extend from each
    acceptor towards its donor according to the acceptor carbon position:
    2 down (270 degrees), 3 lower left (225), 4 left (180), 6 upper left
    (135), and 8 up (90), independent of branching. ``bond_length`` sets the
    horizontal and vertical center-to-center distance in drawing coordinates
    (default 1.2). ``diagonal_bond_scale`` multiplies this length for diagonal
    links (default 1.2, or 20 percent longer); use 1 for equal lengths.
    Alpha linkages are dashed;
    beta and unknown linkages are solid. ``symbol_linewidth`` and ``linkage_linewidth`` set
    stroke widths in points. ``positions`` optionally maps every
    monomer index to (x, y); it is required for non-tree topologies. The
    lines between explicit positions are drawn straight without snapping or
    rescaling. Unknown or other acceptor positions and overlapping fixed-length
    layouts also require explicit positions.
    Linkage annotations give alpha/beta and donor–acceptor carbon positions. Anomers
    that cannot be determined are shown as ?. Nondefault D/L configurations
    and furanose rings are annotated. Figures support normal Matplotlib export.
    """
    if not np.isfinite(symbol_size) or symbol_size <= 0:
        raise ValueError("symbol_size must be finite and positive.")
    for name, value in (("symbol_linewidth", symbol_linewidth), ("linkage_linewidth", linkage_linewidth),
                        ("bond_length", bond_length), ("diagonal_bond_scale", diagonal_bond_scale)):
        if not np.isfinite(value) or value <= 0:
            raise ValueError(f"{name} must be finite and positive.")
    count = len(glycan.monomers)
    identities = identify_sugars(glycan)
    names = [info["name"] if info["name"] in _SYMBOLS else "Unknown" for info in identities]
    if labels is not None:
        if isinstance(labels, Mapping):
            overrides = labels.items()
        else:
            if isinstance(labels, str) or len(labels) != count:
                raise ValueError("labels must contain one SNFG abbreviation per monomer.")
            overrides = enumerate(labels)
        for index, name in overrides:
            if not isinstance(index, (int, np.integer)) or isinstance(index, bool) or not 0 <= index < count:
                raise ValueError("Label keys must be valid monomer indices.")
            if not isinstance(name, str) or name.replace("Nac", "NAc") not in _SYMBOLS:
                raise ValueError(f"Unsupported SNFG abbreviation: {name!r}.")
            names[index] = name.replace("Nac", "NAc")
    if positions is None:
        positions = _layout(glycan, bond_length, diagonal_bond_scale)
    if not isinstance(positions, Mapping) or set(positions) != set(range(count)):
        raise ValueError("positions must map every monomer index to one (x, y) coordinate.")
    xy = np.asarray([positions[i] for i in range(count)], dtype=float)
    if xy.shape != (count, 2) or not np.all(np.isfinite(xy)) or len(np.unique(xy, axis=0)) != count:
        raise ValueError("Symbol positions must be finite, distinct (x, y) coordinates.")
    span = np.ptp(xy, axis=0)
    if ax is None:
        fig, ax = plt.subplots(figsize=(max(2.5, span[0] + 1.5), max(2, span[1] + 1.5)))
    else:
        fig = ax.figure
    for link in glycan.glycosidic_bonds:
        start, end = xy[[link["donor"], link["acceptor"]]]
        anomer = identities[link["donor"]]["anomer"]
        ax.plot([start[0], end[0]], [start[1], end[1]], color="black", linewidth=linkage_linewidth,
                linestyle="--" if anomer == "alpha" else "-", zorder=1)
        if show_linkages:
            anomer = {"alpha": "α", "beta": "β"}.get(anomer, "?")
            donor, acceptor = linkage_positions(glycan, link)
            midpoint = (start + end) / 2
            delta = end - start
            normal = np.array([-delta[1], delta[0]]) / np.linalg.norm(delta)
            if normal[1] < 0:
                normal *= -1
            ax.annotate(f"{anomer}{donor}-{acceptor}", midpoint, xytext=12 * normal,
                        textcoords="offset points", ha="center", va="center", fontsize=10,
                        bbox=dict(facecolor="white", edgecolor="none", pad=.5), zorder=2)
    for index, (name, info) in enumerate(zip(names, identities)):
        _glyph(ax, xy[index], name, symbol_size, symbol_linewidth)
        marks = []
        if info["configuration"] and info["configuration"] != _default_configuration(name):
            marks.append(info["configuration"])
        if info["ring_type"] == "furanose":
            marks.append("f")
        if marks:
            ax.text(*xy[index], "".join(marks), ha="center", va="center", fontsize=9,
                    color="black", zorder=5, bbox=dict(facecolor="white", edgecolor="none", pad=.2))
        if show_names:
            ax.annotate(name, xy[index], xytext=(0, -22), textcoords="offset points",
                        ha="center", va="top", fontsize=10)
    margin = max(.8, symbol_size * 2)
    ax.set(xlim=(xy[:, 0].min() - margin, xy[:, 0].max() + margin),
           ylim=(xy[:, 1].min() - margin, xy[:, 1].max() + margin), aspect="equal")
    ax.set_axis_off()
    return fig, ax

Conformer analysis

Tutorial: Conformers and clustering.

matmalib.conformers.rmsd(first, second, *, atom_indices=None, heavy_only=False, align=True)

Return RMSD in input length units, using a proper rotation when aligned.

Inputs are States or coordinate arrays. Atom correspondence is positional; this does not perform symmetry-based atom permutation or mirror reflection. Neither input is modified.

Source code in matmalib/conformers.py

def rmsd(first, second, *, atom_indices=None, heavy_only=False, align=True):
    """Return RMSD in input length units, using a proper rotation when aligned.

    Inputs are States or coordinate arrays. Atom correspondence is positional;
    this does not perform symmetry-based atom permutation or mirror reflection.
    Neither input is modified.
    """
    indices = _selection(first, atom_indices, heavy_only)
    _compatible(first, second)
    return _aligned_rmsd(_coordinates(first)[indices], _coordinates(second)[indices], align)

matmalib.conformers.rmsd_matrix(states, *, atom_indices=None, heavy_only=True, align=True)

Symmetric pairwise RMSD matrix for a Space or iterable of States.

Source code in matmalib/conformers.py

def rmsd_matrix(states, *, atom_indices=None, heavy_only=True, align=True):
    """Symmetric pairwise RMSD matrix for a Space or iterable of States."""
    states = list(getattr(states, "states", states))
    if not states:
        raise ValueError("At least one conformer is required.")
    indices = _selection(states[0], atom_indices, heavy_only)
    xyz = []
    for state in states:
        _compatible(states[0], state)
        xyz.append(_coordinates(state)[indices])
    distances = np.zeros((len(states), len(states)))
    for first in range(len(states)):
        for second in range(first):
            distances[first, second] = distances[second, first] = _aligned_rmsd(xyz[first], xyz[second], align)
    return distances

matmalib.conformers.cluster_conformers(states, threshold=1.0, *, method='single', heavy_only=True, atom_indices=None, energy='E', energy_window=None, energy_unit='hartree')

Cluster aligned RMSDs with single, complete, average, or weighted linkage.

threshold uses coordinate length units. An optional energy_window is in kcal/mol above the minimum; energy_unit describes stored energies (hartree, kcal/mol, or kJ/mol). Complete finite energies are required for filtering. Representatives are lowest-energy members when all energies are available, otherwise medoids (minimum sum of within-cluster RMSDs). Set energy=None to select medoids even when energies are present.

Source code in matmalib/conformers.py

def cluster_conformers(states, threshold=1.0, *, method="single", heavy_only=True,
                       atom_indices=None, energy="E", energy_window=None,
                       energy_unit="hartree"):
    """Cluster aligned RMSDs with single, complete, average, or weighted linkage.

    ``threshold`` uses coordinate length units. An optional ``energy_window``
    is in kcal/mol above the minimum; ``energy_unit`` describes stored energies
    (hartree, kcal/mol, or kJ/mol). Complete finite energies are required for
    filtering. Representatives are lowest-energy members when all energies
    are available, otherwise medoids (minimum sum of within-cluster RMSDs).
    Set ``energy=None`` to select medoids even when energies are present.
    """
    from .core import Space

    source = states
    states = list(getattr(states, "states", states))
    if not np.isfinite(threshold) or threshold < 0:
        raise ValueError("threshold must be finite and nonnegative.")
    if method not in ("single", "complete", "average", "weighted"):
        raise ValueError("Use single, complete, average, or weighted linkage for RMSD distances.")
    distances = rmsd_matrix(states, atom_indices=atom_indices, heavy_only=heavy_only)
    energies = None
    if energy is not None:
        values = [getattr(state, energy, None) for state in states]
        if all(value is not None and np.ndim(value) == 0 for value in values):
            try:
                numeric = np.asarray(values, dtype=float)
            except (TypeError, ValueError):
                numeric = np.full(len(values), np.nan)
            if np.all(np.isfinite(numeric)):
                energies = numeric
    included = np.arange(len(states))
    if energy_window is not None:
        if not np.isfinite(energy_window) or energy_window < 0:
            raise ValueError("energy_window must be finite and nonnegative.")
        if energies is None:
            raise ValueError("Energy filtering requires finite energies for every conformer.")
        factors = {"hartree": 627.509474, "kcal/mol": 1.0, "kJ/mol": 1 / 4.184}
        if energy_unit not in factors:
            raise ValueError("energy_unit must be hartree, kcal/mol, or kJ/mol.")
        included = included[(energies - energies.min()) * factors[energy_unit] <= energy_window]
    if len(included) == 1:
        tree, raw = np.empty((0, 4)), np.ones(1, dtype=int)
    else:
        condensed = squareform(distances[np.ix_(included, included)], checks=False)
        tree = linkage(condensed, method=method, optimal_ordering=True)
        raw = fcluster(tree, threshold, criterion="distance")
    groups = sorted((tuple(included[raw == label]) for label in np.unique(raw)), key=min)
    labels = np.full(len(states), -1, dtype=int)
    representatives = []
    for label, group in enumerate(groups):
        labels[list(group)] = label
        scores = energies[list(group)] if energies is not None else distances[np.ix_(group, group)].sum(axis=1)
        representatives.append(int(group[int(np.argmin(scores))]))
    selected = Space(states=[states[i] for i in representatives], theory=getattr(source, "theory", None))
    return ClusterResult(labels, tuple(groups), tuple(representatives), selected, distances,
                         tree, included, float(threshold), method)

matmalib.conformers.ClusterResult(labels: np.ndarray, members: tuple, representative_indices: tuple, representatives: object, distances: np.ndarray, linkage: np.ndarray, included_indices: np.ndarray, threshold: float, method: str) dataclass

Cluster labels/members use original indices; excluded labels are -1.

distances covers all input states. Rows of linkage refer to included_indices. representatives is a new Space sharing the selected State objects; no input conformer is reordered or modified.

labels: np.ndarray instance-attribute

members: tuple instance-attribute

representative_indices: tuple instance-attribute

representatives: object instance-attribute

distances: np.ndarray instance-attribute

linkage: np.ndarray instance-attribute

included_indices: np.ndarray instance-attribute

threshold: float instance-attribute

method: str instance-attribute

Plotting

Tutorials: Plotting and styling, Free-energy surfaces, and IR spectra.

matmalib.viz.Plot()

A class that creates plots from a numpy array. It can make reaction coordinate diagrams, RMSD and RMSF trajectory plots, Free energy plots, Scatter Plots, and SNFG Figures.

Construct a plot object with default configuration and colors.

Source code in matmalib/viz.py

def __init__(self):

    """
    Construct a plot object with default configuration and colors.
    """

    config_dict = {
        'xrange': None,
        'yrange': None,
        'xtick': None,
        'ytick': None,
        'xlabel': None,
        'ylabel': None,
        'title': None,
        'font': None,
        'xextend': None,
        'yextend': None,
        'title_fontsize': None,
        'axis_fontsize': None,
        'tick_fontsize': None,
        'headers': None,
    }

    config = Plot.Config(**config_dict)
    self.config = config

    self.config_dict = config_dict

    # RK: Just set some default colors for ease of use, adjust as needed.
    self.colors = ["#4856d1", "#f3a667", '#e31a1c', '#1f78b4', '#53546c']

config = config instance-attribute

config_dict = config_dict instance-attribute

colors = ['#4856d1', '#f3a667', '#e31a1c', '#1f78b4', '#53546c'] instance-attribute

Config(xrange: list = None, yrange: list = None, xtick: list = None, ytick: list = None, xlabel: str = None, ylabel: str = None, title: str = None, font: str = None, xextend: float = None, yextend: float = None, title_fontsize: float = None, axis_fontsize: float = None, tick_fontsize: float = None, headers: list = None) dataclass

xrange: list = None class-attribute instance-attribute

yrange: list = None class-attribute instance-attribute

xtick: list = None class-attribute instance-attribute

ytick: list = None class-attribute instance-attribute

xlabel: str = None class-attribute instance-attribute

ylabel: str = None class-attribute instance-attribute

title: str = None class-attribute instance-attribute

font: str = None class-attribute instance-attribute

xextend: float = None class-attribute instance-attribute

yextend: float = None class-attribute instance-attribute

title_fontsize: float = None class-attribute instance-attribute

axis_fontsize: float = None class-attribute instance-attribute

tick_fontsize: float = None class-attribute instance-attribute

headers: list = None class-attribute instance-attribute

snfg(glycan, **kwargs)

Draw SNFG symbols; return (figure, axes) and retain them on this Plot.

Source code in matmalib/viz.py

def snfg(self, glycan, **kwargs):
    """Draw SNFG symbols; return (figure, axes) and retain them on this Plot."""
    from .snfg import draw_snfg
    self.fig, self.ax = draw_snfg(glycan, **kwargs)
    return self.fig, self.ax

cluster_dendrogram(clusters, *, labels=None, ax=None, **kwargs)

Plot a ClusterResult's linkage tree in RMSD units.

Optional labels cover all original conformers; energy-excluded states are omitted from the tree. Remaining options go to SciPy dendrogram.

Source code in matmalib/viz.py

def cluster_dendrogram(self, clusters, *, labels=None, ax=None, **kwargs):
    """Plot a ClusterResult's linkage tree in RMSD units.

    Optional labels cover all original conformers; energy-excluded states
    are omitted from the tree. Remaining options go to SciPy dendrogram.
    """
    from scipy.cluster.hierarchy import dendrogram
    if labels is None:
        labels = [str(i) for i in range(len(clusters.labels))]
    if len(labels) != len(clusters.labels):
        raise ValueError("Supply one dendrogram label per original conformer.")
    selected = [labels[i] for i in clusters.included_indices]
    if ax is None:
        self.fig, self.ax = plt.subplots()
    else:
        self.fig, self.ax = ax.figure, ax
    if len(selected) == 1:
        if kwargs.get("orientation", "top") in ("left", "right"):
            self.ax.set_yticks([0], selected)
            self.ax.set_xlim(0, max(1, clusters.threshold * 1.1))
        else:
            self.ax.set_xticks([0], selected)
            self.ax.set_ylim(0, max(1, clusters.threshold * 1.1))
    else:
        kwargs.setdefault("color_threshold", clusters.threshold)
        dendrogram(clusters.linkage, labels=selected, ax=self.ax, **kwargs)
    if kwargs.get("orientation", "top") in ("left", "right"):
        self.ax.set_xlabel("RMSD (coordinate units)")
        self.ax.axvline(clusters.threshold, linestyle="--", color="gray", linewidth=1)
    else:
        self.ax.set_ylabel("RMSD (coordinate units)")
        self.ax.axhline(clusters.threshold, linestyle="--", color="gray", linewidth=1)
    return self.fig, self.ax

cmap(color_num: int = None, offset: float = 0, map: str = 'ice', reverse: bool = False)

Generates and processes a colormap with optional offsetting logic.

Parameters:

Name	Type	Description	Default
color_num	int	(int) Number of discrete colors.	None
offset	float	(float) Fractional offset to shift the colormap.	0
map	str	(str) Name of the colormap from the colormaps library. Also configures the palette for cremer_pople_sphere and updates an already-drawn sphere and its colorbar. An explicit cmap argument when drawing a new sphere overrides this palette for that call only.	'ice'

Source code in matmalib/viz.py

def cmap(self, color_num: int = None, offset: float = 0, map: str = 'ice', reverse: bool = False):
    """
    Generates and processes a colormap with optional offsetting logic.
    :param color_num: (int) Number of discrete colors.
    :param offset: (float) Fractional offset to shift the colormap.
    :param map: (str) Name of the colormap from the colormaps library.

    Also configures the palette for cremer_pople_sphere and updates
    an already-drawn sphere and its colorbar. An explicit cmap argument
    when drawing a new sphere overrides this palette for that call only.
    """
    # Check if the colormap exists in colormaps
    if not hasattr(colormaps, map):
        raise ValueError(f"Colormap '{map}' not found in colormaps library!")

    # Fetch colormap
    colors_obj = getattr(colormaps, map)

    if color_num is not None:
        color_num += 1
    # Ensure the colormap has an array of colors
    if not hasattr(colors_obj, 'colors'):
        raise ValueError(f"The selected colormap '{map}' does not have a valid 'colors' attribute!")

    colormap_colors = colors_obj.colors

    # Validating the shape of colormap_colors
    if len(colormap_colors[0]) != 3:
        raise ValueError(f"Expected RGB colors in the colormap, but got shape {np.array(colormap_colors).shape}.")

    # Applying offset manually

    if offset != 0:
        new_colors = []

        for color in colormap_colors:

            new_color = []
            for color_elm in color:
                color_elm -= offset

                if color_elm > 1:
                    color_elm = 1

                if color_elm < 0:
                    color_elm = 0

                new_color.append(color_elm)
            new_colors.append(new_color)

        colormap_colors = new_colors

    if reverse:
        colormap_colors = list(reversed(colormap_colors))

    # Discretize the colormap to the required number of colors
    if color_num is not None:
        discrete_colors = np.linspace(0, len(colormap_colors) - 1, color_num, dtype=int)
        self.colors = [colormap_colors[i] for i in discrete_colors]
    else:
        self.colors = np.asarray(colormap_colors).tolist()

    self._pucker_cmap = ListedColormap(self.colors, name=map)
    heatmap = getattr(self, "heatmap", None)
    if heatmap is not None:
        heatmap.set_cmap(self._pucker_cmap)
        heatmap.figure.canvas.draw_idle()

savefig(filename='fig', format: str = 'png')

Source code in matmalib/viz.py

def savefig(self, filename='fig', format: str = 'png'):
    self.fig.savefig(f"{self.path}/{filename}.{format}", dpi=300, bbox_inches='tight')

set_colors(colors: list = None)

Source code in matmalib/viz.py

def set_colors(self, colors: list = None):
    self.colors = colors
    if colors is not None:
        self._pucker_cmap = ListedColormap(colors, name="custom")
    elif hasattr(self, "_pucker_cmap"):
        del self._pucker_cmap

set_config(conf: dict)

Source code in matmalib/viz.py

def set_config(self, conf: dict):
    old_conf = self.config_dict

    for key, value in old_conf.items():
        if key not in conf.keys():
            conf[key] = value

    self.config = Plot.Config(**conf)

set_axes(ax: matplotlib.pyplot.axes)

Source code in matmalib/viz.py

def set_axes(self, ax: matplotlib.pyplot.axes):
    from matplotlib import rc

    config = self.config
    config_dict = asdict(config)

    for key, value in config_dict.items():

        if key == 'xrange' and value is not None:
            ax.set_xlim(value[0], value[1])
        if key == 'yrange' and value is not None:
            ax.set_ylim(value[0], value[1])
        if key == 'xticks' and value is not None:
            ax.set_xticks(value)
        if key == 'yticks' and value is not None:
            ax.set_yticks(value)
        if key == 'xlabel' and value is not None:
            ax.set_xlabel(value)
        if key == 'ylabel' and value is not None:
            ax.set_ylabel(value)
        if key == 'title' and value is not None:
            ax.set_title(value)
        if key == 'font' and value is not None:
            mpl.rcParams['font.sans-serif'] = value
            mpl.rcParams['font.family'] = "sans-serif"
        if key == 'axis fontsize' and value is not None:
            mpl.rcParams['axes.labelsize'] = value
        if key == 'title fontsize' and value is not None:
            ax.title.set_size(value)
        if key == 'tick fontsize' and value is not None:
            ax.tick_params(labelsize=value, axis='both')

trajectory(data_list, var_name='colvar', col=1, average=None, title=None, hist=True, alpha=None, calc_qa=False, overlap=False, headers=None, *, projection=None, monomer=0, pucker_labels=False, cmap=None, bins=(72, 36), num_levels=8, quantity='free_energy', temp=300, energy_unit='kcal/mol')

Plots MD trajectory with histogram. Takes in data for CP2K or Gromacs via Mol.

Parameters:

Name	Type	Description	Default
data_list		(Mol / List) Either a Data object, or a list of Data objs if you want to overlay data.	required
var_name		(list) Name of the collective variable you are plotting on your y-axis.	'colvar'
col		(int) Index of the column containing your colvar data, in the case that you have multiple. For a Cremer–Pople free-energy map, use cremer_pople_sphere. projection="mollweide" forwards to that method, including monomer, pucker_labels, cmap, bins, num_levels, quantity, temp and energy_unit. In that view col, var_name and hist do not apply; smoothing, overlap and calc_qa are unsupported. With projection=None or "rectilinear", the time-series view is used.	1

Source code in matmalib/viz.py

def trajectory(self, data_list, var_name='colvar', col=1, average=None, title=None, hist=True, alpha=None,
               calc_qa=False, overlap=False, headers=None, *, projection=None, monomer=0,
               pucker_labels=False, cmap=None, bins=(72, 36), num_levels=8,
               quantity="free_energy", temp=300, energy_unit="kcal/mol"):

    """ Plots MD trajectory with histogram. Takes in data for CP2K or Gromacs via Mol.
    :param data_list: (Mol / List) Either a Data object, or a list of Data objs if you want to overlay data.
    :param var_name: (list) Name of the collective variable you are plotting on your y-axis.
    :param col: (int) Index of the column containing your colvar data, in the case that you have multiple.

    For a Cremer–Pople free-energy map, use cremer_pople_sphere.
    ``projection="mollweide"`` forwards to that method, including monomer,
    pucker_labels, cmap, bins, num_levels, quantity, temp and energy_unit.
    In that view col, var_name and hist do not
    apply; smoothing, overlap and calc_qa are unsupported. With
    projection=None or "rectilinear", the time-series view is used.
    """
    if projection not in (None, "rectilinear", "mollweide"):
        raise ValueError("projection must be None, 'rectilinear', or 'mollweide'.")
    if not isinstance(data_list, list):
        data_list = [data_list]

    if headers is None and self.config.headers is not None:
        headers = self.config.headers

    if projection == "mollweide":
        if (average is not None and np.any(np.asarray(average) > 1)) or overlap or calc_qa:
            raise ValueError("Mollweide heatmaps do not support average, overlap, or calc_qa.")
        return self.cremer_pople_sphere(data_list, monomer=monomer, title=title, alpha=alpha,
                                       headers=headers, pucker_labels=pucker_labels, cmap=cmap,
                                       bins=bins, num_levels=num_levels, quantity=quantity, temp=temp,
                                       energy_unit=energy_unit)

    self.path = data_list[0].path
    # CP2K default timestep unit is in fs, Gromacs is in ps:
    # We convert these to ps and nm respectively:

    #         if mol_list[0].software == 'cp2`k':
    #             time_unit = 'ps'`

    #         elif mol_list[0].software == 'gromacs':
    #             time_unit = self.time_unit

    fig, ax = plt.subplots(1, 2, figsize=(11, 3), gridspec_kw={'width_ratios': [3.5, 1]})

    i = 0

    if alpha == None:
        alpha = [0.8] * len(data_list)

    if average == None:
        average = [0] * len(data_list)

    elif not isinstance(average, list):
        average = [average] * len(data_list)

    if headers == None:
        headers = [None for i in range(len(data_list))]

    for header, datum in zip(headers, data_list):

        residue = False

        if datum.time_unit == 'fs':
            time = (datum.data[:, 0] / 1000).tolist()  # fs -> ps for CP2K
            time_label = 'ps'

        elif datum.time_unit == 'ps':
            time = (datum.data[:, 0] / 1000).tolist()  # ps -> ns for GROMACS
            time_label = 'ns'

        elif datum.time_unit == 'ns':
            time = (datum.data[:, 0]).tolist()  # if GROMACS already in ns don't convert
            time_label = 'ns'

        elif datum.time_unit == 'Residue':
            time = (datum.data[:, 0]).tolist()  # if GROMACS in residue, also don't convert
            time_label = 'residue'
            residue = True

        colvar = datum.data[:, col].tolist()

        timestep = np.abs(time[0] - time[1])

        if not hasattr(self, 'colors'):
            self.cmap(color_num=len(data_list), offset=0, map='ice')

        color = self.colors

        if average[i] > 1:
            array_len = len(colvar)
            # conv_kernel = np.ones(average[i])/array_len
            conv_kernel = np.ones(average[i]) / average[i]
            colvar_conv = np.convolve(colvar, conv_kernel, mode='valid').tolist()
            time = time[:-1 * average[i] + 1]

        if overlap == True:

            ax[0].plot(time, colvar_conv, linewidth=2, color=color[i], alpha=alpha[i], label=header)
            ax[0].plot(time, colvar[:(len(colvar_conv))], linewidth=0.8, color=color[i], alpha=alpha[i] * .3)
            ax[1].hist(colvar, bins='rice', orientation="horizontal", color=color[i], alpha=alpha[i])

        elif overlap == False and average[i] > 1:
            ax[0].plot(time, colvar_conv, linewidth=0.8, color=color[i], alpha=alpha[i], label=header)

        else:
            ax[0].plot(time, colvar, linewidth=0.8, color=color[i], alpha=alpha[i], label=header)
            ax[1].hist(colvar, bins='rice', orientation="horizontal", color=color[i], alpha=alpha[i],
                       label=np.round(np.average(colvar)))

        if calc_qa == True:

            nbins = 50
            hist = np.histogram(colvar[500:], nbins, range=(min(colvar), max(colvar)))

            dmin = np.argmin(hist[0][15:23]) + 15
            bs = np.sum(hist[0][:dmin + 1]);
            us = np.sum(hist[0][dmin + 1:])

            if us == 0:
                Qa = 1000.0; boundary = 0.0

            else:
                Qa = float(bs) / float(us);
                boundary = float(dmin) / 10

            # ax[1].fill_between([0, ax[1].get_xlim()[1]], boundary, boundary+0.1, color='0.8')

            # Only annotate Qas if one trajectory is entered, otherwise print them.
            if len(data_list) == 1:
                ax[1].axhline(y=boundary, color='gray', linestyle='-', alpha=0.5, linewidth=5)

                textstr = r'$Q_a$={0:3.2f}'.format(Qa)
                ax[1].text(0.55 * ax[1].get_xlim()[1], 0.95 * ax[1].get_ylim()[1], textstr, fontsize=14,
                           verticalalignment='top')

            else:
                print(f"mol{i + 1} Qa = {np.round(Qa, 3)}")

        # if len(data_list) == 1:
        #    ax[1].set_title(f"average = {np.round(np.average(colvar), 3)}", fontsize = 10)

        # else:

        # x= (np.round(np.average(colvar), 3))
        # ax[1].legend()
        # print(x)

        # print(f"average = {np.round(np.average(colvar), 2)}")

        i = i + 1

    if headers[0] is not None:
        ax[0].legend(loc='upper right')

    if not residue:
        ax[0].set_xlabel(f"time ({time_label}); stepsize = {timestep}{time_label}")

    else:
        ax[0].set_xlabel(f"{time_label}")

    ax[0].set_ylabel(var_name)

    if title != None:
        ax[0].set_title(f"{title}", fontsize=10)

    if hist == False:
        fig.delaxes(ax[1])

    if not residue:
        xmax = ax[0].get_xlim()[1]
        xmax = xmax + 1
        ax[0].set_xlim(0, xmax)

    ax[1].set_xlabel('structures')

    self.set_axes(ax[0])

    if ax[1].get_ylim() != ax[0].get_ylim():
        ax[1].set_ylim(ax[0].get_ylim())

    # Hard code the y axis of the histogram to align with the trajectory:
    for key, value in self.config_dict.items():
        if key == 'yrange' and value is not None:
            ax[1].set_ylim(value[0], value[1])

    plt.tight_layout()
    self.fig = fig
    self.ax = ax

cremer_pople_sphere(data_list, *, monomer=0, cmap=None, title=None, alpha=None, headers=None, pucker_labels=False, figsize=(8, 4), bins=(72, 36), num_levels=8, quantity='free_energy', temp=300, energy_unit='kcal/mol')

Plot relative free energy on a Mollweide Cremer–Pople sphere.

Accept a Data or list of Data returned by cremer_pople_trajectory. Select a zero-based six-membered monomer (default 0). Named theta and phi columns are required; angle_unit is degrees or radians. Each valid row contributes one conformer, irrespective of time or order. Multiple datasets are pooled into one histogram. bins is a positive integer or a (phi, theta) pair of bin counts; the default (72, 36) gives 5-degree cells. By default, each occupied bin has free energy -R * temp * log(N / Nmax), using R = 8.314 J/(mol K) and temp=300 kelvin, matching Plot.fes. energy_unit defaults to "kcal/mol"; "kJ/mol" is also supported (1 kcal = 4.184 kJ). The most populated bin is zero; empty bins are masked. This uses counts per angular cell, with no spherical-area correction or reweighting. quantity="counts" restores the population map.

Filled contours interpolate between sampled bin centers using num_levels color bands (default 8). In free-energy mode, a finer drawing grid interpolates energies using occupied neighbors and rounds the boundary at half occupancy. Isolated samples remain visible and empty bin centers stay white. This changes only the rendering, not the stored populations or energies. A light dashed grid marks CP angles. Use the palette configured by Plot.cmap() or Plot.set_colors(), or a pastel green palette by default. An explicit cmap Matplotlib name or object overrides it for this plot. The vertical colorbar on the right shows free energy in the selected unit, or Number of conformers in count mode.

Longitude is 180-phi degrees and latitude is 90-theta; tick labels show CP angles in degrees. Phi wraps at 0/360 degrees. Pure chairs are counted at the poles using phi=180, where phase is arbitrary; planar frames and other undefined angles are omitted. Set pucker_labels=True to annotate the 38 ideal puckers. Labels assume O, C1, ..., C5 ring order, with O denoting position 6 in classification results.

figsize defaults to (8, 4) inches and alpha to 1. Alpha must be common to all pooled datasets. Headers are joined into a default title. Supports title and label/font configuration; configured axis limits are unsupported. Return (figure, axes), also stored as fig and ax on this Plot; heatmap and colorbar hold the filled contours and Colorbar. conformer_counts holds the exact integer histogram in displayed (latitude, longitude) order, with latitude_edges and longitude_edges in radians. In free-energy mode, free_energy holds the masked energy grid in the same order and in the selected energy_unit.

Source code in matmalib/viz.py

def cremer_pople_sphere(self, data_list, *, monomer=0, cmap=None, title=None,
                       alpha=None, headers=None, pucker_labels=False, figsize=(8, 4),
                       bins=(72, 36), num_levels=8, quantity="free_energy", temp=300,
                       energy_unit="kcal/mol"):
    """Plot relative free energy on a Mollweide Cremer–Pople sphere.

    Accept a Data or list of Data returned by cremer_pople_trajectory.
    Select a zero-based six-membered ``monomer`` (default 0). Named theta
    and phi columns are required; angle_unit is degrees or radians. Each
    valid row contributes one conformer, irrespective of time or order.
    Multiple datasets are pooled into one histogram. ``bins`` is a positive
    integer or a (phi, theta) pair of bin counts; the default (72, 36) gives
    5-degree cells. By default, each occupied bin has free energy
    ``-R * temp * log(N / Nmax)``, using R = 8.314 J/(mol K) and
    ``temp=300`` kelvin, matching Plot.fes. ``energy_unit`` defaults to
    "kcal/mol"; "kJ/mol" is also supported (1 kcal = 4.184 kJ).
    The most populated bin is
    zero; empty bins are masked. This uses counts per angular cell, with
    no spherical-area correction or reweighting. ``quantity="counts"``
    restores the population map.

    Filled contours interpolate between sampled bin centers using
    ``num_levels`` color bands (default 8). In free-energy mode, a finer
    drawing grid interpolates energies using occupied neighbors and rounds
    the boundary at half occupancy. Isolated samples remain visible and
    empty bin centers stay white. This changes only the rendering, not the
    stored populations or energies. A light dashed grid marks CP angles.
    Use the palette configured by Plot.cmap() or Plot.set_colors(), or a
    pastel green palette by default. An explicit ``cmap`` Matplotlib name
    or object overrides it for this plot. The vertical colorbar on the
    right shows free energy in the selected unit, or Number of conformers
    in count mode.

    Longitude is 180-phi degrees and latitude is 90-theta; tick labels show
    CP angles in degrees. Phi wraps at 0/360 degrees. Pure chairs are
    counted at the poles using phi=180, where phase is arbitrary; planar
    frames and other undefined angles are omitted. Set pucker_labels=True
    to annotate the 38 ideal puckers. Labels assume
    O, C1, ..., C5 ring order, with O denoting position 6 in classification
    results.

    ``figsize`` defaults to (8, 4) inches and ``alpha`` to 1. Alpha must be
    common to all pooled datasets. Headers are joined into a default title.
    Supports title and label/font configuration; configured axis limits
    are unsupported. Return (figure, axes), also stored as fig and ax on
    this Plot; heatmap and colorbar hold the filled contours and Colorbar.
    conformer_counts holds the exact integer histogram in displayed
    (latitude, longitude) order, with latitude_edges and longitude_edges
    in radians. In free-energy mode, free_energy holds the masked energy
    grid in the same order and in the selected energy_unit.
    """
    from ._pucker_reference import CREMER_POPLE_PUCKERS

    if not isinstance(data_list, list):
        data_list = [data_list]
    if not data_list:
        raise ValueError("Supply at least one Cremer–Pople trajectory.")
    if isinstance(monomer, bool) or not isinstance(monomer, Integral) or monomer < 0:
        raise ValueError("monomer must be a nonnegative integer.")
    if self.config.xrange is not None or self.config.yrange is not None:
        raise ValueError("Mollweide axes have fixed limits; remove xrange and yrange from the plot config.")
    try:
        bin_counts = (bins, bins) if isinstance(bins, Integral) else tuple(bins)
    except TypeError:
        raise ValueError("bins must be a positive integer or a (phi, theta) pair of positive integers.") from None
    if len(bin_counts) != 2 or any(isinstance(n, bool) or not isinstance(n, Integral) or n <= 0
                                   for n in bin_counts):
        raise ValueError("bins must be a positive integer or a (phi, theta) pair of positive integers.")
    if isinstance(num_levels, bool) or not isinstance(num_levels, Integral) or num_levels <= 0:
        raise ValueError("num_levels must be a positive integer.")
    if quantity not in ("free_energy", "counts"):
        raise ValueError("quantity must be 'free_energy' or 'counts'.")
    if not isinstance(temp, Real) or isinstance(temp, bool) or not np.isfinite(temp) or temp <= 0:
        raise ValueError("temp must be finite and positive, in kelvin.")
    if energy_unit not in ("kcal/mol", "kJ/mol"):
        raise ValueError("energy_unit must be 'kcal/mol' or 'kJ/mol'.")
    if headers is None and self.config.headers is not None:
        headers = self.config.headers
    if headers is None:
        headers = [None] * len(data_list)
    if len(headers) != len(data_list):
        raise ValueError("Supply one header per dataset.")
    if alpha is None or np.isscalar(alpha):
        alpha = [1.0 if alpha is None else alpha] * len(data_list)
    if len(alpha) != len(data_list):
        raise ValueError("Supply one alpha value per dataset, or a single shared value.")
    alpha = np.asarray(alpha, dtype=float)
    if not np.all(np.isfinite(alpha) & (alpha >= 0) & (alpha <= 1)):
        raise ValueError("alpha must lie between 0 and 1.")
    if not np.all(alpha == alpha[0]):
        raise ValueError("Pooled heatmaps require the same alpha for all datasets.")

    if cmap is None:
        colormap = getattr(self, "_pucker_cmap", None)
        if colormap is None:
            colormap = LinearSegmentedColormap.from_list(
                "pucker_green", ["#b7efbe", "#88e2b2", "#58cfa4"]
            )
    else:
        colormap = mpl.colormaps.get_cmap(cmap)
    phi_edges = np.linspace(0, 2 * np.pi, bin_counts[0] + 1)
    theta_edges = np.linspace(0, np.pi, bin_counts[1] + 1)
    counts = np.zeros(bin_counts, dtype=np.int64)
    for datum in data_list:
        columns = tuple(getattr(datum, "columns", ()))
        theta_name, phi_name = f"monomer_{monomer}.theta", f"monomer_{monomer}.phi"
        if theta_name not in columns or phi_name not in columns:
            raise ValueError(
                f"Mollweide projection requires six-membered Cremer–Pople data with "
                f"{theta_name} and {phi_name} columns."
            )
        values = np.asarray(datum.data, dtype=float)
        if values.ndim != 2 or values.shape[0] == 0 or values.shape[1] != len(columns):
            raise ValueError("Trajectory data must be a nonempty 2D array matching its columns.")
        theta, phi = values[:, [columns.index(theta_name), columns.index(phi_name)]].T
        unit = getattr(datum, "angle_unit", None)
        if unit == "degrees":
            theta, phi = np.radians(theta), np.radians(phi)
        elif unit != "radians":
            raise ValueError("Cremer–Pople angle_unit must be 'degrees' or 'radians'.")
        if np.any(np.isfinite(theta) & ((theta < 0) | (theta > np.pi))):
            raise ValueError("Cremer–Pople theta must lie between 0 and 180 degrees (0 and pi radians).")

        # Phase is arbitrary at a pole: count all pure chairs together.
        poles = np.isclose(theta, 0, atol=1e-12, rtol=0) | np.isclose(theta, np.pi, atol=1e-12, rtol=0)
        phi = np.where(poles, np.pi, phi)
        valid = np.isfinite(theta) & np.isfinite(phi)
        histogram, _, _ = np.histogram2d(phi[valid] % (2 * np.pi), theta[valid],
                                        bins=(phi_edges, theta_edges))
        counts += histogram.astype(np.int64)

    maximum = int(counts.max())
    fig, ax = plt.subplots(figsize=figsize, subplot_kw={"projection": "mollweide"},
                           layout="constrained")
    fig.set_facecolor("white")
    ax.set_facecolor("white")
    # Both CP angles decrease along the displayed axes.
    self.conformer_counts = counts.T[::-1, ::-1].copy()
    self.longitude_edges = np.pi - phi_edges[::-1]
    self.latitude_edges = np.pi / 2 - theta_edges[::-1]
    longitude = (self.longitude_edges[:-1] + self.longitude_edges[1:]) / 2
    latitude = (self.latitude_edges[:-1] + self.latitude_edges[1:]) / 2

    self.free_energy = None
    if quantity == "free_energy":
        occupied = self.conformer_counts > 0
        self.free_energy = np.ma.masked_all(self.conformer_counts.shape, dtype=float)
        if maximum:
            gas_constant = .008314 / (4.184 if energy_unit == "kcal/mol" else 1)
            self.free_energy[occupied] = gas_constant * temp * (
                np.log(maximum) - np.log(self.conformer_counts[occupied])
            )
        upper = max(1., float(np.ceil(self.free_energy.max()))) if maximum else 1.
        norm = Normalize(0, upper)
        levels = np.linspace(0, upper, num_levels + 1)
        longitude, latitude, field = _pucker_energy_contours(
            self.free_energy, longitude, latitude,
        )
    else:
        norm = Normalize(0, max(1, maximum))
        # A half-count threshold retains isolated samples in count mode.
        levels = np.linspace(.5, maximum + .5, min(num_levels, maximum) + 1) if maximum else [0, 1]
        values = self.conformer_counts if maximum else np.ma.masked_all(self.conformer_counts.shape)

        # Wrap the seam and extend edge rows to the poles in count mode.
        seam = (values[:, :1] + values[:, -1:]) / 2
        field = np.ma.concatenate((seam, values, seam), axis=1)
        field = np.ma.concatenate((field[:1], field, field[-1:]))
        longitude = np.r_[-np.pi, longitude, np.pi]
        latitude = np.r_[-np.pi / 2, latitude, np.pi / 2]
    heatmap = ax.contourf(longitude, latitude, field, levels=levels,
                          cmap=colormap, norm=norm, alpha=alpha[0], zorder=1,
                          corner_mask=True)
    heatmap.set_transform(_PuckerContourTransform(ax.transProjection) + ax.transAffine + ax.transAxes)

    ax.set_longitude_grid(30)
    ax.set_latitude_grid(30)
    # Geographic axes place the secondary x labels below the equator,
    # leaving space above it for boat and skew-boat annotations.
    ax.tick_params(axis="x", labeltop=True, labelbottom=False)
    ax.tick_params(axis="both", labelsize=8)
    ax.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: f"{180 - np.degrees(x):.0f}°"))
    ax.yaxis.set_major_formatter(FuncFormatter(lambda y, pos: f"{90 - np.degrees(y):.0f}°"))
    ax.set_xlabel(r"$\phi$ (degrees)", fontsize=9)
    ax.set_ylabel(r"$\theta$ (degrees)", fontsize=9)
    ax.set_axisbelow(False)
    ax.grid(True, color="0.75", linewidth=.65, linestyle="--", zorder=3)
    ax.spines["geo"].set_color("0.3")
    ax.spines["geo"].set_linewidth(.8)

    if pucker_labels:
        reference_points = []
        for label, (theta, phi) in CREMER_POPLE_PUCKERS.items():
            family = next(code for code in "CBSHE" if code in label)
            above, below = label.replace("6", "O").split(family)
            text = rf"$^{{{above}}}\mathrm{{{family}}}_{{{below}}}$"
            # Show seam references at both edges, where they identify the
            # same point on the sphere. Offset labels into the map.
            longitudes = (-np.pi, np.pi) if phi == 0 else (np.radians(180 - phi),)
            for longitude in longitudes:
                latitude = np.radians(90 - theta)
                reference_points.append((longitude, latitude))
                dx = 4 if longitude == -np.pi else -4 if longitude == np.pi else 0
                dy = -7 if theta == 0 else -10 if family == "E" else 6
                ax.annotate(text, (longitude, latitude), xytext=(dx, dy),
                            textcoords="offset points", fontsize=8, color="0.25",
                            ha="left" if dx > 0 else "right" if dx < 0 else "center",
                            va="top" if theta == 0 else "bottom", zorder=4,
                            bbox=dict(facecolor="white", edgecolor="none", alpha=0.7, pad=0.5))
        ax.scatter(*np.asarray(reference_points).T, s=6, color="0.45", zorder=2)

    if title is None and any(header is not None for header in headers):
        title = ", ".join(str(header) for header in headers if header is not None)
    if title is not None:
        ax.set_title(title, fontsize=11, pad=6)
    self.set_axes(ax)
    self.heatmap = heatmap
    if quantity == "free_energy":
        ticks = np.linspace(0, upper, min(num_levels, 5) + 1) if maximum else [0]
        colorbar_label = f"Free energy ({energy_unit})"
    else:
        ticks = MaxNLocator(nbins=5, integer=True).tick_values(0, max(1, maximum))
        ticks = ticks[(ticks >= 1) & (ticks <= maximum)] if maximum else [0]
        colorbar_label = "Number of conformers"
    self.colorbar = fig.colorbar(heatmap, ax=ax,
                                 orientation="vertical", fraction=0.035, pad=0.035,
                                 shrink=.8, ticks=ticks)
    self.colorbar.set_label(colorbar_label, fontsize=9)
    self.colorbar.ax.tick_params(labelsize=8)
    self.colorbar.outline.set_linewidth(.7)
    self.path = getattr(data_list[0], "path", "")
    self.fig, self.ax = fig, ax
    return fig, ax

fes_1d(data, *, smooth_window=None, show_raw=True, show_basins=True, energy_unit='kJ/mol', ax=None, color='#235789', label=None)

Plot the 1D FES returned by mlmd.analyze_metadynamics.

Uses the analysis script's local cubic smoothing, basin shading and projected-barrier marker. Smoothing preserves unsampled gaps and does not modify the raw profile or its thermodynamics. Set smooth_window=0 to disable it. Return (figure, axes), also retained on this Plot. Plot.fes continues to handle two-dimensional population plots.

Source code in matmalib/viz.py

def fes_1d(self, data, *, smooth_window=None, show_raw=True, show_basins=True,
           energy_unit="kJ/mol", ax=None, color="#235789", label=None):
    """Plot the 1D FES returned by ``mlmd.analyze_metadynamics``.

    Uses the analysis script's local cubic smoothing, basin shading and
    projected-barrier marker. Smoothing preserves unsampled gaps and does
    not modify the raw profile or its thermodynamics. Set smooth_window=0
    to disable it. Return (figure, axes), also retained on this Plot.
    ``Plot.fes`` continues to handle two-dimensional population plots.
    """
    from ._mlmd_analysis_core import Profile, smooth_profile
    factors = {"kJ/mol": 1.0, "kcal/mol": 1 / 4.184, "eV": 1 / 96.48533212331002}
    if energy_unit not in factors:
        raise ValueError("energy_unit must be kJ/mol, kcal/mol, or eV.")
    if not hasattr(data, "free_energy") or not hasattr(data, "x"):
        raise TypeError("fes_1d expects Data from mlmd.analyze_metadynamics().")
    x = np.asarray(data.x, dtype=float)
    raw = np.asarray(data.free_energy, dtype=float)
    if (x.ndim != 1 or raw.shape != x.shape or len(x) < 2 or not np.isfinite(x).all()
            or np.any(np.diff(x) <= 0) or not np.isfinite(raw).any()):
        raise ValueError("A 1D FES needs an increasing CV grid and sampled free energies.")
    if getattr(data, "energy_unit", "kJ/mol") != "kJ/mol":
        raise ValueError("Analysis input energies must be in kJ/mol.")
    window = getattr(data, "smooth_window", 11) if smooth_window is None else smooth_window
    if isinstance(window, bool) or not isinstance(window, Integral):
        raise ValueError("smooth_window must be an integer.")
    smooth = smooth_profile(Profile(x, raw), window)
    if ax is None:
        self.fig, self.ax = plt.subplots(figsize=(8, 4.5), constrained_layout=True)
    else:
        self.fig, self.ax = ax.figure, ax
    factor = factors[energy_unit]
    if window and show_raw:
        self.ax.plot(x, np.where(np.isfinite(raw), raw * factor, np.nan),
                     color="#9aa5b1", linewidth=1, alpha=0.45, label="Raw FES")
    self.ax.plot(x, np.where(np.isfinite(smooth), smooth * factor, np.nan),
                 color=color, linewidth=2,
                 label=label or ("Smoothed FES" if window else "FES"))
    summary = getattr(data, "summary", {})
    if show_basins:
        for name, shade in (("reactant", "#38995b"), ("product", "#d48122")):
            if f"{name}_range" in summary:
                text = name.capitalize()
                if summary.get("basin_assignment") == "inferred":
                    text += " (inferred)"
                self.ax.axvspan(*summary[f"{name}_range"], color=shade, alpha=0.15, label=text)
        if "transition_state_cv" in summary:
            self.ax.axvline(summary["transition_state_cv"], color="#a13c4b",
                            linestyle="--", label="Projected barrier")
    self.ax.set(xlabel=f"{getattr(data, 'cv', 'CV')} [{getattr(data, 'cv_unit', 'input units')}]",
                ylabel=f"Relative free energy [{energy_unit}]")
    if window or label or (show_basins and "reactant_range" in summary):
        self.ax.legend(frameon=False)
    self.set_axes(self.ax)
    self.path = getattr(data, "path", None)
    self.free_energy = raw.copy()
    return self.fig, self.ax

fes(data, cols=[1, 2], temp=300, num_levels=8)

Plots MD FES. Takes in data for CP2K or Gromacs via Mol.

Parameters:

Name	Type	Description	Default
data		(Data) Class Data.	required
cols		(int) Index of the 2 columns containing your colvar data, in the case that you have more than 2.	[1, 2]

Source code in matmalib/viz.py

def fes(self, data, cols=[1, 2], temp=300, num_levels=8):

    """ Plots MD FES. Takes in data for CP2K or Gromacs via Mol.
    :param data: (Data) Class Data.
    :param cols: (int) Index of the 2 columns containing your colvar data, in the case that you have more than 2.
    """

    self.path = data.path
    Temp = temp;
    R = 8.314  # J/K mol

    colvar1 = data.data[:, cols[0]]
    colvar2 = data.data[:, cols[1]]

    Hall, x_edges, y_edges = np.histogram2d(colvar1, colvar2, bins=72)

    Hall = - R * Temp * np.log(Hall)
    hmin = np.min(Hall)

    Hall_rel = 0.001 * (Hall.T - hmin)

    vmin, vmax = 0, np.ceil(np.nanmax(Hall_rel[~np.isinf(Hall_rel)]))
    MHall = np.ma.masked_greater(Hall_rel, vmax)

    fig, ax = plt.subplots(figsize=(6, 6))

    # colors = self.colors
    # cmap = ListedColormap(colors)

    self.set_axes(ax)

    num_levels = num_levels
    plot = ax.contourf(x_edges[:-1], y_edges[:-1], MHall, colors=self.colors, zorder=1, levels=num_levels)
    # plot = ax.contourf(x_edges[:-1], y_edges[:-1], Hall.T, cmap=cmap, zorder=1, levels=num_levels)

    cb_ticks = np.linspace(vmin, vmax, 6)
    cb = fig.colorbar(plot, ax=ax, ticks=cb_ticks, pad=0.05, shrink=0.6)
    cb.ax.set_yticklabels([f"{tick:.1f}" for tick in cb_ticks], fontsize=12)
    cb.set_label("\n Free energy [kJ]", fontsize=14)

    # Enable grid that aligns with ticks
    ax.grid(True, ls='--', zorder=10.0)

    x_range = ax.get_xlim()
    y_range = ax.get_ylim()

    dx = x_range[1] - x_range[0]
    dy = y_range[1] - y_range[0]

    ax.set_aspect(dx/dy, adjustable='box')

    fig.tight_layout()
    self.fig = fig
    self.ax = ax

contour(xpm_data_list, limit=16)

Plots contour maps for a provided list of moles from xpm files.

Parameters:

Name	Type	Description	Default
xpm_data_list		Data object or list of Data objects read from XPM files.	required
limit		(Int) The upper limit on your energy scale	16

Source code in matmalib/viz.py

def contour(self, xpm_data_list, limit=16):

    """ Plots contour maps for a provided list of moles from xpm files.
    :param xpm_data_list: Data object or list of Data objects read from XPM files.
    :param limit: (Int) The upper limit on your energy scale
    """

    if not isinstance(xpm_data_list, list):
        xpm_data_list = [xpm_data_list]

    if len(xpm_data_list) == 1:
        self.path = xpm_data_list[0].path
    else:
        self.path = xpm_data_list[0].path

    Mat = []

    for xpm_data in xpm_data_list:
        M = xpm_data.data
        MM = np.ma.masked_greater(M, limit - 1)
        Mat.append(MM)

    #         if xpm_mol1 != None and xpm_mol2 != None and xpm_mol3 != None:
    #             M1 = xpm_mol1.data
    #             M2 = xpm_mol2.data
    #             M3 = xpm_mol3.data
    #             MM1 = np.ma.masked_greater(M1, limit-1)
    #             MM2 = np.ma.masked_greater(M2, limit-1)
    #             MM3 = np.ma.masked_greater(M3, limit-1)
    #             Mat = [MM1, MM2, MM3]

    #         elif xpm_mol1 != None and xpm_mol2 != None:
    #             M1 = xpm_mol1.data
    #             M2 = xpm_mol2.data
    #              #DiffM = M1-M2
    #             MM1 = np.ma.masked_greater(M1, limit-1)
    #             MM2 = np.ma.masked_greater(M2, limit-1)
    #             Mat = [MM1, MM2]

    #         else:
    #             raise ValueError("Give me a correct number of some chunky matrices")
    fig, axes = plt.subplots(1, len(Mat), figsize=(4 * len(Mat) + (len(Mat) - 1) * 1, 4), sharex=True, sharey=True)
    # fig, axes = plt.subplots(1,len(Mat), figsize=(4*len(Mat) + 1.5, 4), sharex=True, sharey=True)

    levels = np.linspace(0, limit, 9)  # 8 levels between 0 and limit
    if len(Mat) > 1:
        for n, ax in enumerate(axes):

            self.set_axes(ax)
            dx = ax.get_xlim()[1] - ax.get_xlim()[0]
            dy = ax.get_ylim()[1] - ax.get_ylim()[0]

            ax.set_aspect(dx/dy, adjustable='box')
            ax.grid(True, ls='--', zorder=10.0)

            # Create the contourf plot with consistent levels
            plot = ax.contourf(Mat[n], levels=levels, colors=self.colors, zorder=1)

            # Add a color bar with consistent boundaries and ticks
            cb = fig.colorbar(plot, ax=ax, pad=0.025, aspect=20, ticks=levels, shrink=0.7)
            cb.set_ticklabels(["{0:3.1f}".format(x) for x in levels])
    else:
        self.set_axes(axes)
        dx = axes.get_xlim()[1] - axes.get_xlim()[0]
        dy = axes.get_ylim()[1] - axes.get_ylim()[0]

        axes.set_aspect(dx / dy, adjustable='box')
        axes.grid(True, ls='--', zorder=10.0)

        # Create the contourf plot with consistent levels
        plot = axes.contourf(Mat[0], levels=levels, colors=self.colors, zorder=1)

        # Add a color bar with consistent boundaries and ticks
        cb = fig.colorbar(plot, ax=axes, pad=0.025, aspect=20, ticks=levels, shrink=0.7)
        cb.set_ticklabels(["{0:3.1f}".format(x) for x in levels])

    fig.tight_layout()

    self.fig = fig
    self.ax = axes

foo_plot(data, SCR=None, w=0.5)

Not a bar plot. Mostly used for SCR stuff

Parameters:

Name	Type	Description	Default
data		(Data) Data objects generated from csv files.	required
SCR		(list) Specifies specific SCRs to plot from your data. Otherwise, default is to plot them all.	None
w		(float) Width of the lines in the foo plot.	0.5

Source code in matmalib/viz.py

def foo_plot(self, data, SCR=None, w=0.5):

    """ Not a bar plot. Mostly used for SCR stuff
    :param data: (Data) Data objects generated from csv files.
    :param SCR: (list) Specifies specific SCRs to plot from your data. Otherwise, default is to plot them all.
    :param w: (float) Width of the lines in the foo plot.
    """

    color = self.colors
    self.path = data.path

    if SCR == None:

        SCR = []
        for line in data.data[1:]:
            SCR.append(line.split()[0])

    SUG = data.data[0].split();
    data_dict = {}

    for line in data.data[1:]:
        # Doctor Founder says not to transform the data, only plot the data:
        # data[line.split()[0]] = [log(float(i)*(float(i)+1)/C0,10) for i in line.split()[1:]]

        data_dict[line.split()[0]] = [float(i) for i in line.split()[1:]]

    fig, ax = plt.subplots(figsize=(8.0, 4.0))

    ax.tick_params(axis='both', which='both', bottom=True, top=False, labelbottom=True, right=False, left=False,
                   labelleft=True, labelright=False)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['left'].set_visible(False)
    ax.xaxis.set_tick_params(direction='out')

    x_pos = np.arange(len(SUG))
    xmax = len(SUG)

    ax.set_xlim(-1, xmax)

    ymin = 0;
    ymax = np.round(max(val for sublist in data_dict.values() for val in sublist))
    ax.set_ylim(ymin, ymax + 0.1)
    yticks = np.linspace(ymin, ymax, 7)
    ax.set_yticks(yticks)
    # ax.set_yticklabels(yticks)

    ax.set_xticks(x_pos)
    ax.set_xticklabels(SUG)

    # ax.set_xlabel(r'time [ns]')
    for i in yticks:  ax.plot([-1, xmax], [i, i], '0.75', lw=0.5)

    ax.grid(axis='y', color='grey', linestyle='-', linewidth=0.5)

    for n, scr in enumerate(SCR):
        for i in range(len(SUG)):
            ax.plot([x_pos[i] - w / 2, x_pos[i] + w / 2], [data_dict[scr][i], data_dict[scr][i]], color=color[n],
                    lw=2)

            # ax.bar ( x_pos[i], data[scr][i], align='center', color=color[n])

    plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")

    self.set_axes(ax)

    fig.tight_layout()
    self.fig = fig
    self.ax = ax

scatter(data=None, headers=None, format: str = '.', figsize=(6, 6))

Generates a scatter plot from data

Source code in matmalib/viz.py

def scatter(self, data=None, headers=None, format: str = '.', figsize=(6, 6)):

    """
    Generates a scatter plot from data
    """

    if headers is not None and self.config.headers is not None:
        headers = self.config.headers

    data_list = []
    if isinstance(data, list):
        for set in data:
            data_list.append(set.data)
    else:
        data = data.data
        data_list.append(data)

    if headers != None:
        desc = headers

    # x_extend = 0
    # y_extend = 0

    colors = self.colors if self.colors is not None else ['b', 'r', 'g', 'c', 'm', 'y', 'k']
    fig, ax = plt.subplots(1, 1, figsize=figsize)

    for i, data in enumerate(data_list):

        if data.ndim == 2:
            data_x = data[:, 0]
            data_ys = []
            row_len = len(data[0, :])
        else:
            data_x = data[0]
            data_ys = []
            row_len = len(data)

        nullfmt = NullFormatter()  # no labels

        # fig, ax = plt.subplots(1, figsize=figsize)

        for col in range(1, row_len):
            if data.ndim == 2:
                data_y = data[:, col]
                data_ys.append(data_y)
            else:
                data_y = data[col]
                data_ys.append(data_y)

            fit = np.polyfit(data_x, data_y, 1)
            val = np.polyval(fit, data_x)

            if headers is not None:
                if len(data_list) == 1:
                    ax.scatter(data_x, data_y, marker=format, label=desc[col - 1], color=colors[col - 1])
                    ax.legend()
                else:
                    ax.scatter(data_x, data_y, marker=format, label=desc[i], color=colors[i])
                    ax.legend()

            else:
                if len(data_list) == 1:
                    ax.scatter(data_x, data_y, marker=format, color=colors[col - 1])
                else:
                    ax.scatter(data_x, data_y, marker=format, color=colors[i])

    ax.tick_params(axis='both', which='both', bottom=True, top=False, labelbottom=True, right=False, left=True,
                   labelleft=True)
    for s in ['top', 'right']:
        ax.spines[s].set_visible(False)
    ax.xaxis.set_tick_params(direction='out')
    ax.yaxis.set_tick_params(direction='out')

    self.set_axes(ax)

    xmin, xmax = ax.get_xlim()[0], ax.get_xlim()[1]
    ymin, ymax = ax.get_ylim()[0], ax.get_ylim()[1]

    # if x_extend != 0.0 :
    #    xmin -= x_extend ; xmax += x_extend
    #    ax.set_xlim(xmin, xmax)
    #    print(xmin, xmax)

    # if y_extend != 0.0:
    #    ax.set_ylim(ymin-y_extend, ymax)

    # ax.plot([xmin, xmax], [ymin-y_extend, ymin-y_extend], 'k', lw=1)
    # ax.plot([xmin+ax_fit, xmin+ax_fit], [ymin ,ymax], 'k', lw=1)

    plt.tight_layout()
    self.fig = fig
    self.ax = ax

gaussian_broadening(state, broaden=1, resolution=1)

Performs gaussian broadening on IR spectrum generates attribute self.IR - np.array with dimmension 4000/resolution consisting gaussian-boraden spectrum

Parameters:

Name	Type	Description	Default
broaden		(float) gaussian broadening in wn-1	1
resolution		(float) resolution of the spectrum (number of points for 1 wn) defaults is 1, needs to be fixed in plotting	1

Source code in matmalib/viz.py

def gaussian_broadening(self, state, broaden=1, resolution=1):

    """ Performs gaussian broadening on IR spectrum
    generates attribute self.IR - np.array with dimmension 4000/resolution consisting gaussian-boraden spectrum

    :param broaden: (float) gaussian broadening in wn-1
    :param resolution: (float) resolution of the spectrum (number of points for 1 wn) defaults is 1, needs to be fixed in plotting
    """

    IR = np.zeros((int(4000 / resolution) + 1,))
    X = np.linspace(0, 4000, int(4000 / resolution) + 1)
    for f, i in zip(state.Freq, state.Ints):  IR += i * np.exp(-0.5 * ((X - f) / int(broaden)) ** 2)
    state.IR = np.vstack((X, IR)).T  # tspec

ir(state, xmin=900, xmax=1700, normal_modes: bool = False)

Plots the IR spectrum in xmin -- xmax range, x-axis is multiplied by scaling factor, everything is normalized to 1. If exp_data is specified, then the top panel is getting plotted too. Need to add output directory. Default name is self._id

Source code in matmalib/viz.py

def ir(self, state, xmin=900, xmax=1700, normal_modes: bool = False):

    """ Plots the IR spectrum in xmin -- xmax range,
    x-axis is multiplied by scaling factor, everything
    is normalized to 1. If exp_data is specified,
    then the top panel is getting plotted too.
    Need to add output directory. Default name is self._id
    """

    import matplotlib.pyplot as plt
    from matplotlib.ticker import NullFormatter
    import math

    colors = self.colors
    background_colors = []

    for color in colors:
        color[0] += 0.25
        color[1] += 0.25
        color[2] += 0.25
        background_colors.append(color)

    fig, ax = plt.subplots(1, figsize=(math.ceil(10 * (xmax - xmin) / (1500)), 3))

    self.gaussian_broadening(state, resolution=0.001, broaden=2)

    # left, width = 0.02, 0.98 ; bottom, height = 0.15, 0.8
    # ax  = [left, bottom, width, height ]
    # ax  = plt.axes(ax)
    exten = 20

    ax.tick_params(axis='both', which='both', bottom=True, top=False, labelbottom=True, right=False, left=False,
                   labelleft=False)
    ax.spines['top'].set_visible(False);
    ax.spines['right'].set_visible(False);
    ax.spines['left'].set_visible(False)
    ax.xaxis.set_tick_params(direction='out')
    ax.yaxis.set_major_formatter(NullFormatter())
    ax.set_ylim(0, 1.1)

    xticks = np.linspace(xmin, xmax, int((xmax - xmin + 2 * exten) / 100) + 1)
    ax.set_xticks(xticks)
    ax.set_xticklabels([int(x) for x in xticks], fontsize=10)
    ax.set_xlim(xmin - exten, xmax + exten + 10)

    shift = 0.05
    incr = (state.IR[-1, 0] - state.IR[0, 0]) / (len(state.IR) - 1)
    scale_t = 1 / np.amax(state.IR[int(xmin / incr):int(xmax / incr) + 100, 1])

    Xsc = state.IR[:, 0];
    IRsc = state.IR[:, 1] * scale_t
    ir_theo = ax.plot(Xsc, IRsc + shift, color='k', linewidth=0.5)
    ax.fill_between(Xsc, np.linspace(shift, shift, len(IRsc)), IRsc + shift, color=background_colors[0], alpha=0.5)

    if normal_modes == True:
        for l in range(len(state.Freq)):
            ax.plot([state.Freq[l], state.Freq[l]], [shift, state.Ints[l] * scale_t + shift], linewidth=2,
                    color=colors[0])

    fig.tight_layout()
    current_path = os.getcwd()

space(list_of_spaces, names=None, quant: str = 'E')

Source code in matmalib/viz.py

def space(self, list_of_spaces, names=None, quant: str = 'E'):

    fig, ax = plt.subplots()

    n = 0

    for space in list_of_spaces:

        # if quant == 'E': data = space.energies
        # elif quant == 'F': data = space.free_energies
        # elif quant == 'H': data = space.enthalpies

        data = space.diff(quant=quant)

        if names == None:
            name = [str(n)] * len(data)
        else:
            name = names[n]

        ax.scatter(name, data, marker='_', linewidths=3, s=700)

        n = n + 1

    if n != 1:
        pad = 1 / (n - 1)
        ax.set_xmargin(pad)

    plt.ylabel(f"{quant} (Ha)")
    plt.xlabel('conformer id')

    self.set_axes(ax)

Machine-learning molecular dynamics

The mlmd interface provides geometry minimization, equilibration, production dynamics, checkpoint restart, periodic solvation, and one-dimensional metadynamics analysis. Its configuration classes are nested under mlmd; simulation dependencies load only when required.

Tutorials: run a simulation, solvate a system, metadynamics, restart, and analyze free energy. The configuration reference explains every configuration class with source-derived defaults; the script feature map records the integration's scope.

matmalib.mlmd.mlmd(model_path=None, *, device=None, compile_model=False, calculator=None, model_kwargs=None)

Simulation interface accepting XYZ/PDB paths or matmalib State objects.

model_path accepts one AIMNet inference file, or a sequence of files whose predictions are averaged. With no path, all four AIMNet2-2025 B97-3c/D3 members are used. Registry downloads are managed by AIMNet. Other MLIPs can be used by passing an ASE calculator instance. Input systems are copied; no trajectory is accumulated in memory.

Source code in matmalib/mlmd.py

def __init__(self, model_path=None, *, device=None, compile_model=False,
             calculator=None, model_kwargs=None):
    if calculator is not None and model_path is not None:
        raise ValueError("Choose model_path or calculator, not both.")
    self.model_path = model_path
    self.device = device
    self.compile_model = compile_model
    self.calculator = calculator
    self.model_kwargs = dict(model_kwargs or {})
    self._calculators = {}
    from .core import Data
    self.data = Data()  # structured diagnostics from the most recent operation

model_path = model_path instance-attribute

device = device instance-attribute

compile_model = compile_model instance-attribute

calculator = calculator instance-attribute

model_kwargs = dict(model_kwargs or {}) instance-attribute

data = Data() instance-attribute

SystemConfig(cell: Any = None, pbc: bool | tuple[bool, bool, bool] | None = None, charge: float | None = None, multiplicity: int | None = None, min_distance: float = 0.5, fixed_atoms: tuple[int, ...] = (), fix_com: bool = False, shake_hydrogens: bool = False, shake_cutoff: float = 1.3, shake_tolerance: float = 1e-06) dataclass

cell: Any = None class-attribute instance-attribute

pbc: bool | tuple[bool, bool, bool] | None = None class-attribute instance-attribute

charge: float | None = None class-attribute instance-attribute

multiplicity: int | None = None class-attribute instance-attribute

min_distance: float = 0.5 class-attribute instance-attribute

fixed_atoms: tuple[int, ...] = () class-attribute instance-attribute

fix_com: bool = False class-attribute instance-attribute

shake_hydrogens: bool = False class-attribute instance-attribute

shake_cutoff: float = 1.3 class-attribute instance-attribute

shake_tolerance: float = 1e-06 class-attribute instance-attribute

MinimizationConfig(cell: Any = None, pbc: bool | tuple[bool, bool, bool] | None = None, charge: float | None = None, multiplicity: int | None = None, min_distance: float = 0.5, fixed_atoms: tuple[int, ...] = (), fix_com: bool = False, shake_hydrogens: bool = False, shake_cutoff: float = 1.3, shake_tolerance: float = 1e-06, optimizer: str = 'FIRE', fmax: float = 0.05, steps: int = 10000, maxstep: float = 0.1, optimizer_kwargs: dict[str, Any] = dict(), require_convergence: bool = True, trajectory: str | PathLike | None = None, trajectory_interval: int = 1, log_file: str | PathLike | None = None, overwrite: bool = False) dataclass

Bases: SystemConfig

optimizer: str = 'FIRE' class-attribute instance-attribute

fmax: float = 0.05 class-attribute instance-attribute

steps: int = 10000 class-attribute instance-attribute

maxstep: float = 0.1 class-attribute instance-attribute

optimizer_kwargs: dict[str, Any] = field(default_factory=dict) class-attribute instance-attribute

require_convergence: bool = True class-attribute instance-attribute

trajectory: str | PathLike | None = None class-attribute instance-attribute

trajectory_interval: int = 1 class-attribute instance-attribute

log_file: str | PathLike | None = None class-attribute instance-attribute

overwrite: bool = False class-attribute instance-attribute

MetadynamicsConfig(plumed_file: str | PathLike | None = None, input_lines: tuple[str, ...] | None = None, save_input: str | PathLike | None = None, cv: str = 'DISTANCE', cv_label: str = 'cv', atoms: tuple[int, ...] = (1, 2), metad_label: str = 'metad', metad_arg: str | None = None, sigma: float = 0.02, height_kj_mol: float = 2.0, pace: int = 300, biasfactor: float = 10.0, temperature_K: float | None = None, hills_file: str | PathLike = 'HILLS', colvar_file: str | PathLike = 'COLVAR', colvar_stride: int = 10, colvar_arg: str | None = None, wall_type: str | None = None, wall_label: str = 'wall', wall_arg: str | None = None, wall_at: float = 0.35, wall_kappa: float = 150000.0, wall_exp: int = 2, wall_eps: float = 1.0, wall_offset: float = 0.0, calc_rct: bool | None = None, grid_auto: bool | None = None, grid_min: float = 0.0, grid_max: float | None = None, grid_bins: int | None = None, grid_spacing: float = 0.002, log_file: str | PathLike | None = None) dataclass

plumed_file: str | PathLike | None = None class-attribute instance-attribute

input_lines: tuple[str, ...] | None = None class-attribute instance-attribute

save_input: str | PathLike | None = None class-attribute instance-attribute

cv: str = 'DISTANCE' class-attribute instance-attribute

cv_label: str = 'cv' class-attribute instance-attribute

atoms: tuple[int, ...] = (1, 2) class-attribute instance-attribute

metad_label: str = 'metad' class-attribute instance-attribute

metad_arg: str | None = None class-attribute instance-attribute

sigma: float = 0.02 class-attribute instance-attribute

height_kj_mol: float = 2.0 class-attribute instance-attribute

pace: int = 300 class-attribute instance-attribute

biasfactor: float = 10.0 class-attribute instance-attribute

temperature_K: float | None = None class-attribute instance-attribute

hills_file: str | PathLike = 'HILLS' class-attribute instance-attribute

colvar_file: str | PathLike = 'COLVAR' class-attribute instance-attribute

colvar_stride: int = 10 class-attribute instance-attribute

colvar_arg: str | None = None class-attribute instance-attribute

wall_type: str | None = None class-attribute instance-attribute

wall_label: str = 'wall' class-attribute instance-attribute

wall_arg: str | None = None class-attribute instance-attribute

wall_at: float = 0.35 class-attribute instance-attribute

wall_kappa: float = 150000.0 class-attribute instance-attribute

wall_exp: int = 2 class-attribute instance-attribute

wall_eps: float = 1.0 class-attribute instance-attribute

wall_offset: float = 0.0 class-attribute instance-attribute

calc_rct: bool | None = None class-attribute instance-attribute

grid_auto: bool | None = None class-attribute instance-attribute

grid_min: float = 0.0 class-attribute instance-attribute

grid_max: float | None = None class-attribute instance-attribute

grid_bins: int | None = None class-attribute instance-attribute

grid_spacing: float = 0.002 class-attribute instance-attribute

log_file: str | PathLike | None = None class-attribute instance-attribute

MDConfig(cell: Any = None, pbc: bool | tuple[bool, bool, bool] | None = None, charge: float | None = None, multiplicity: int | None = None, min_distance: float = 0.5, fixed_atoms: tuple[int, ...] = (), fix_com: bool = False, shake_hydrogens: bool = False, shake_cutoff: float = 1.3, shake_tolerance: float = 1e-06, steps: int = 10000, timestep_fs: float = 0.5, temperature_K: float = 300.0, ensemble: str = 'NVT', thermostat: str = 'langevin', barostat: str = 'mtk', thermostat_time_fs: float = 100.0, barostat_time_fs: float = 1000.0, pressure_bar: float = 1.0, compressibility_per_bar: float = 4.57e-05, thermostat_chain: int = 3, barostat_chain: int = 3, seed: int | None = None, initialize_velocities: bool | None = None, remove_com_motion: bool = True, hydrogen_mass: float | None = None, heavy_mass_multiplier: float = 1.0, max_velocity_angstrom_fs: float | None = None, max_force_eV_A: float | None = 150.0, trajectory: str | PathLike = 'trajectory.xyz', trajectory_interval: int = 10, wrap_output: bool = False, log_file: str | PathLike | None = 'md.log', log_interval: int = 10, checkpoint_interval: int = 1000, checkpoint_file: str | PathLike | None = 'auto', mode: str = 'md', metadynamics: mlmd.MetadynamicsConfig | None = None, minimization: mlmd.MinimizationConfig | None = None, equilibration: mlmd.EquilibrationConfig | None = None, continue_equilibration_bias: bool = False, overwrite: bool = False) dataclass

Bases: SystemConfig

steps: int = 10000 class-attribute instance-attribute

timestep_fs: float = 0.5 class-attribute instance-attribute

temperature_K: float = 300.0 class-attribute instance-attribute

ensemble: str = 'NVT' class-attribute instance-attribute

thermostat: str = 'langevin' class-attribute instance-attribute

barostat: str = 'mtk' class-attribute instance-attribute

thermostat_time_fs: float = 100.0 class-attribute instance-attribute

barostat_time_fs: float = 1000.0 class-attribute instance-attribute

pressure_bar: float = 1.0 class-attribute instance-attribute

compressibility_per_bar: float = 4.57e-05 class-attribute instance-attribute

thermostat_chain: int = 3 class-attribute instance-attribute

barostat_chain: int = 3 class-attribute instance-attribute

seed: int | None = None class-attribute instance-attribute

initialize_velocities: bool | None = None class-attribute instance-attribute

remove_com_motion: bool = True class-attribute instance-attribute

hydrogen_mass: float | None = None class-attribute instance-attribute

heavy_mass_multiplier: float = 1.0 class-attribute instance-attribute

max_velocity_angstrom_fs: float | None = None class-attribute instance-attribute

max_force_eV_A: float | None = 150.0 class-attribute instance-attribute

trajectory: str | PathLike = 'trajectory.xyz' class-attribute instance-attribute

trajectory_interval: int = 10 class-attribute instance-attribute

wrap_output: bool = False class-attribute instance-attribute

log_file: str | PathLike | None = 'md.log' class-attribute instance-attribute

log_interval: int = 10 class-attribute instance-attribute

checkpoint_interval: int = 1000 class-attribute instance-attribute

checkpoint_file: str | PathLike | None = 'auto' class-attribute instance-attribute

mode: str = 'md' class-attribute instance-attribute

metadynamics: mlmd.MetadynamicsConfig | None = None class-attribute instance-attribute

minimization: mlmd.MinimizationConfig | None = None class-attribute instance-attribute

equilibration: mlmd.EquilibrationConfig | None = None class-attribute instance-attribute

continue_equilibration_bias: bool = False class-attribute instance-attribute

overwrite: bool = False class-attribute instance-attribute

EquilibrationConfig(cell: Any = None, pbc: bool | tuple[bool, bool, bool] | None = None, charge: float | None = None, multiplicity: int | None = None, min_distance: float = 0.5, fixed_atoms: tuple[int, ...] = (), fix_com: bool = False, shake_hydrogens: bool = False, shake_cutoff: float = 1.3, shake_tolerance: float = 1e-06, steps: int = 2000, timestep_fs: float = 0.5, temperature_K: float = 300.0, ensemble: str = 'NVT', thermostat: str = 'langevin', barostat: str = 'mtk', thermostat_time_fs: float = 10.0, barostat_time_fs: float = 1000.0, pressure_bar: float = 1.0, compressibility_per_bar: float = 4.57e-05, thermostat_chain: int = 3, barostat_chain: int = 3, seed: int | None = None, initialize_velocities: bool | None = None, remove_com_motion: bool = True, hydrogen_mass: float | None = None, heavy_mass_multiplier: float = 1.0, max_velocity_angstrom_fs: float | None = None, max_force_eV_A: float | None = None, trajectory: str | PathLike | None = None, trajectory_interval: int = 10, wrap_output: bool = False, log_file: str | PathLike | None = 'equilibration.log', log_interval: int = 10, checkpoint_interval: int = 1000, checkpoint_file: str | PathLike | None = 'auto', mode: str = 'md', metadynamics: mlmd.MetadynamicsConfig | None = None, minimization: mlmd.MinimizationConfig | None = None, equilibration: mlmd.EquilibrationConfig | None = None, continue_equilibration_bias: bool = False, overwrite: bool = False) dataclass

Bases: MDConfig

steps: int = 2000 class-attribute instance-attribute

thermostat_time_fs: float = 10.0 class-attribute instance-attribute

max_force_eV_A: float | None = None class-attribute instance-attribute

trajectory: str | PathLike | None = None class-attribute instance-attribute

log_file: str | PathLike | None = 'equilibration.log' class-attribute instance-attribute

SolvationConfig(cell: float | tuple[float, float, float] = 30.0, density_g_ml: float = 1.0, counterion: Any = None, counterion_gap: float = 0.5, packmol: str | PathLike = 'packmol', tolerance: float = 2.0, nloop: int = 1000, timeout: float = 600.0, seed: int = 12345, volume_samples: int = 200000, vdw_radii: dict[str, float] = dict(), work_dir: str | PathLike | None = None, dry_run: bool = False, charge: float | None = None, multiplicity: int = 1, overwrite: bool = False) dataclass

cell: float | tuple[float, float, float] = 30.0 class-attribute instance-attribute

density_g_ml: float = 1.0 class-attribute instance-attribute

counterion: Any = None class-attribute instance-attribute

counterion_gap: float = 0.5 class-attribute instance-attribute

packmol: str | PathLike = 'packmol' class-attribute instance-attribute

tolerance: float = 2.0 class-attribute instance-attribute

nloop: int = 1000 class-attribute instance-attribute

timeout: float = 600.0 class-attribute instance-attribute

seed: int = 12345 class-attribute instance-attribute

volume_samples: int = 200000 class-attribute instance-attribute

vdw_radii: dict[str, float] = field(default_factory=dict) class-attribute instance-attribute

work_dir: str | PathLike | None = None class-attribute instance-attribute

dry_run: bool = False class-attribute instance-attribute

charge: float | None = None class-attribute instance-attribute

multiplicity: int = 1 class-attribute instance-attribute

overwrite: bool = False class-attribute instance-attribute

AnalysisConfig(colvar: Any = None, hills: str | PathLike | None = None, fes_file: str | PathLike | None = None, plumed_input: str | PathLike | None = None, plumed: str | PathLike | None = None, cv: str | None = None, temperature: float | None = None, energy_unit: str | None = None, time_unit: str | None = None, cv_unit: str | None = None, bins: int = 100, grid_range: tuple[float, float] | None = None, energy_column: str = 'file.free', rbias_column: str | None = None, bias_column: str | None = None, rct_column: str | None = None, additional_bias: tuple[str, ...] | None = None, start_time: float | None = None, end_time: float | None = None, upper_wall: tuple[float, float] | None = None, keep_wall: bool = False, reactant: tuple[float, float] | None = None, product: tuple[float, float] | None = None, kappa: float = 1.0, infrequent: bool = False, acceleration_column: str | None = None, smooth_window: int = 11, plot: bool = False, output_dir: str | PathLike | None = None, overwrite: bool = False) dataclass

colvar: Any = None class-attribute instance-attribute

hills: str | PathLike | None = None class-attribute instance-attribute

fes_file: str | PathLike | None = None class-attribute instance-attribute

plumed_input: str | PathLike | None = None class-attribute instance-attribute

plumed: str | PathLike | None = None class-attribute instance-attribute

cv: str | None = None class-attribute instance-attribute

temperature: float | None = None class-attribute instance-attribute

energy_unit: str | None = None class-attribute instance-attribute

time_unit: str | None = None class-attribute instance-attribute

cv_unit: str | None = None class-attribute instance-attribute

bins: int = 100 class-attribute instance-attribute

grid_range: tuple[float, float] | None = None class-attribute instance-attribute

energy_column: str = 'file.free' class-attribute instance-attribute

rbias_column: str | None = None class-attribute instance-attribute

bias_column: str | None = None class-attribute instance-attribute

rct_column: str | None = None class-attribute instance-attribute

additional_bias: tuple[str, ...] | None = None class-attribute instance-attribute

start_time: float | None = None class-attribute instance-attribute

end_time: float | None = None class-attribute instance-attribute

upper_wall: tuple[float, float] | None = None class-attribute instance-attribute

keep_wall: bool = False class-attribute instance-attribute

reactant: tuple[float, float] | None = None class-attribute instance-attribute

product: tuple[float, float] | None = None class-attribute instance-attribute

kappa: float = 1.0 class-attribute instance-attribute

infrequent: bool = False class-attribute instance-attribute

acceleration_column: str | None = None class-attribute instance-attribute

smooth_window: int = 11 class-attribute instance-attribute

plot: bool = False class-attribute instance-attribute

output_dir: str | PathLike | None = None class-attribute instance-attribute

overwrite: bool = False class-attribute instance-attribute

minimize(system, config=None, *, output=None, output_format=None)

Relax geometry; return a State, or a Path when output is supplied.

output_format='state' returns a State even when writing an XYZ. output_format='xyz' requires an output path. Nonconvergence raises unless require_convergence=False; the result records convergence.

Source code in matmalib/mlmd.py

def minimize(self, system, config=None, *, output=None, output_format=None):
    """Relax geometry; return a State, or a Path when ``output`` is supplied.

    ``output_format='state'`` returns a State even when writing an XYZ.
    ``output_format='xyz'`` requires an output path. Nonconvergence raises
    unless ``require_convergence=False``; the result records convergence.
    """
    self._require_ase()
    from ._mlmd_dynamics import minimize
    return minimize(self, system, config or self.MinimizationConfig(),
                    output, output_format)

equilibrate(system, config=None, *, output=None, output_format=None)

Run equilibration and return its final State or XYZ path.

Final masses and momenta are retained for a subsequent production run. An optional equilibration trajectory is streamed to its configured path. Checkpoints save positions and velocities every 1,000 steps and at completion by default, with or without PLUMED.

Source code in matmalib/mlmd.py

def equilibrate(self, system, config=None, *, output=None, output_format=None):
    """Run equilibration and return its final State or XYZ path.

    Final masses and momenta are retained for a subsequent production run.
    An optional equilibration trajectory is streamed to its configured path.
    Checkpoints save positions and velocities every 1,000 steps and at
    completion by default, with or without PLUMED.
    """
    self._require_ase()
    from ._mlmd_dynamics import equilibrate
    return equilibrate(self, system, config or self.EquilibrationConfig(),
                       output, output_format)

run_md(system, config=None, *, trajectory=None) -> Path

Run normal or biased dynamics; return the streamed XYZ trajectory path.

All ensembles save positions/velocities checkpoints every 1,000 steps and at completion by default. This includes ordinary MD and generated or custom PLUMED input; PLUMED is not needed for ordinary MD checkpointing.

Source code in matmalib/mlmd.py

def run_md(self, system, config=None, *, trajectory=None) -> Path:
    """Run normal or biased dynamics; return the streamed XYZ trajectory path.

    All ensembles save positions/velocities checkpoints every 1,000 steps and
    at completion by default. This includes ordinary MD and generated or
    custom PLUMED input; PLUMED is not needed for ordinary MD checkpointing.
    """
    self._require_ase()
    from dataclasses import replace
    from ._mlmd_dynamics import run_md
    config = config or self.MDConfig()
    if trajectory is not None:
        config = replace(config, trajectory=trajectory)
    return run_md(self, system, config)

restart_md(checkpoint, *, steps, trajectory=None, log_file=None, checkpoint_file='auto') -> Path

Continue a checkpoint for steps additional integration steps.

Restore velocities, RNG and thermostat/barostat state. Output is a new XYZ segment with continuous step/time labels. Use the same calculator and model files as the original run. Generated metadynamics input also restores its saved HILLS/COLVAR history into a new directory.

Source code in matmalib/mlmd.py

def restart_md(self, checkpoint, *, steps, trajectory=None, log_file=None,
               checkpoint_file="auto") -> Path:
    """Continue a checkpoint for ``steps`` additional integration steps.

    Restore velocities, RNG and thermostat/barostat state. Output is a new
    XYZ segment with continuous step/time labels. Use the same calculator
    and model files as the original run. Generated metadynamics input also
    restores its saved HILLS/COLVAR history into a new directory.
    """
    self._require_ase()
    from ._mlmd_checkpoint import restart
    return restart(self, checkpoint, steps, trajectory, log_file, checkpoint_file)

solvate(solute, solvent, config=None, *, output=None, output_format=None)

Pack a periodic solvent box; inputs accept XYZ, PDB, or State.

Density is total mass / box volume. The final State contains a solvation_report with counts, actual density and packing diagnostics. dry_run=True requires work_dir and returns Data with reusable packing inputs and the report; omit output/output_format for a dry run.

Source code in matmalib/mlmd.py

def solvate(self, solute, solvent, config=None, *, output=None, output_format=None):
    """Pack a periodic solvent box; inputs accept XYZ, PDB, or State.

    Density is total mass / box volume. The final State contains a
    ``solvation_report`` with counts, actual density and packing diagnostics.
    ``dry_run=True`` requires ``work_dir`` and returns Data with reusable
    packing inputs and the report; omit output/output_format for a dry run.
    """
    self._require_ase()
    from ._mlmd_solvation import solvate
    return solvate(self, solute, solvent, config or self.SolvationConfig(),
                   output, output_format)

analyze_metadynamics(config=None, *, output_format='data')

Analyze a 1D FES and optional infrequent kinetics from PLUMED files.

Return a matmalib Data object accepted by Plot.fes_1d. It is also stored as self.data. With output_format='files', require config.output_dir and return that directory after writing CSV/JSON.

Source code in matmalib/mlmd.py

def analyze_metadynamics(self, config=None, *, output_format="data"):
    """Analyze a 1D FES and optional infrequent kinetics from PLUMED files.

    Return a matmalib Data object accepted by ``Plot.fes_1d``. It is also
    stored as ``self.data``. With ``output_format='files'``, require
    ``config.output_dir`` and return that directory after writing CSV/JSON.
    """
    from ._mlmd_analysis import analyze
    return analyze(self, config or self.AnalysisConfig(), output_format)

Element mappings

matmalib.elements.ELEMENT_SYMBOLS maps string atomic numbers "1"–"118" to element symbols. ELEMENT_NUMBERS maps those symbols back to integer atomic numbers. State.element_symbol, StateBuilder.element_symbol, and InputBuilder.element_number wrap these mappings; they do not supply atomic masses. See element utilities.