Minimizer

Monte Carlo energy minimization for DNA structures using the PMCpy engine.

mdna.minimizer.Minimizer

Monte Carlo energy minimization of DNA structures using the PMCpy engine.

The Minimizer relaxes a Nucleic object by running rigid base-pair Monte Carlo simulations with sequence-dependent elastic potentials and optional excluded-volume interactions. Three equilibration strategies are available:

• Full equilibration (default) — adaptive convergence monitoring.
• Simple equilibration — fixed-sweep protocol, faster but less thorough.
• Writhe equilibration — maintains linking number while equilibrating writhe for circular DNA (requires compiled PyLk Cython extensions).

After minimization the internal frames are updated in-place and the full MC trajectory can be retrieved via :meth:get_MC_traj.

Parameters:

Name	Type	Description	Default
nucleic	Nucleic	A :class:~mdna.nucleic.Nucleic instance whose attributes (frames, sequence, topology, …) are copied into the Minimizer.	required

Raises:

Type	Description
ImportError	If the PMCpy Run class cannot be imported.

Source code in mdna/minimizer.py

class Minimizer:
    """Monte Carlo energy minimization of DNA structures using the PMCpy engine.

    The Minimizer relaxes a ``Nucleic`` object by running rigid base-pair Monte
    Carlo simulations with sequence-dependent elastic potentials and optional
    excluded-volume interactions.  Three equilibration strategies are available:

    * **Full equilibration** (default) — adaptive convergence monitoring.
    * **Simple equilibration** — fixed-sweep protocol, faster but less thorough.
    * **Writhe equilibration** — maintains linking number while equilibrating
      writhe for circular DNA (requires compiled PyLk Cython extensions).

    After minimization the internal frames are updated in-place and the full MC
    trajectory can be retrieved via :meth:`get_MC_traj`.

    Args:
        nucleic: A :class:`~mdna.nucleic.Nucleic` instance whose attributes
            (frames, sequence, topology, …) are copied into the Minimizer.

    Raises:
        ImportError: If the PMCpy ``Run`` class cannot be imported.
    """

    def __init__(self, nucleic: 'Nucleic') -> None:
        # Dynamically set attributes from the nucleic instance
        self.__dict__.update(nucleic.__dict__)

        # Check if the required import is available
        if not self._check_import():
            raise ImportError("Run class from pmcpy.run.run is not available.")

    def _check_import(self) -> bool:
        """Check whether the PMCpy ``Run`` class is importable.

        Returns:
            bool: True if import succeeds, False otherwise.
        """
        try:
            from .PMCpy.pmcpy.run.run import Run
            self.Run = Run  # Store the imported class in the instance
            return True
        except ImportError as e:
            print(f"ImportError: {e}")
            return False

    def _initialize_mc_engine(self):
        """Create and return a PMCpy ``Run`` instance configured from the current state.

        Returns:
            Run: An initialized PMCpy Monte Carlo engine.
        """
        # Get the positions and triads of the current frame
        pos = self.frames[:,self.frame,0,:]
        triads = self.frames[:,self.frame,1:,:].transpose(0,2,1) # flip row vectors to column vectors
        print('Circular:',self.circular)
        # Initialize the Monte Carlo engine
        mc = self.Run(triads=triads,positions=pos,
                        sequence=self.sequence,
                        closed=self.circular,
                        endpoints_fixed=self.endpoints_fixed,
                        fixed=self.fixed,
                        temp=self.temperature,
                        exvol_rad=self.exvol_rad)
        return  mc

    def _update_frames(self) -> None:
        """Write optimized positions and triads back into the stored frames array."""
        # update the spline with new positions and triads
        self.frames[:,self.frame,0,:] = self.out['positions'] # set the origins of the frames
        self.frames[:,self.frame,1:,:] = self.out['triads'].transpose(0,2,1) # set the triads of the frames as row vectors

    def _get_positions_and_triads(self):
        """Extract final positions and triads from the MC output.

        Also stores the last-frame triads and positions in ``self.out`` for
        subsequent use by :meth:`_update_frames`.

        Returns:
            tuple[numpy.ndarray, numpy.ndarray]: ``(positions, triads)`` arrays
                with shapes ``(n_frames, n_bp, 3)`` and
                ``(n_frames, n_bp, 3, 3)`` respectively, where triads are
                returned as row-vector convention.
        """
        # get the positions and triads of the simulation
        positions = self.out['confs'][:,:,:3,3] 
        triads = self.out['confs'][:,:,:3,:3]

        # get the last frames of the simulation
        self.out['triads'] = triads[-1]
        self.out['positions'] = positions[-1]
        return positions, triads.transpose(0,1,3,2) # flip column vectors to row vectors

    def minimize(self,  frame: int = -1, exvol_rad : float = 2.0, temperature : int = 300,  simple : bool = False, equilibrate_writhe : bool = False, endpoints_fixed : bool = True, fixed : List[int] = [], dump_every : int = 20, plot : bool = False) -> None:
        """Run Monte Carlo energy minimization on the DNA structure.

        Initializes the PMCpy engine with the current base-pair frames and
        sequence, performs MC equilibration, and updates the internal frames
        with the relaxed configuration.

        Args:
            frame: Index of the stored frame to minimize (default ``-1``, the
                last frame).
            exvol_rad: Excluded-volume radius in nm.  Base pairs closer than
                this distance incur a repulsive energy penalty.
            temperature: Temperature in Kelvin for the Metropolis acceptance
                criterion.
            simple: If True, use simple (fixed-sweep) equilibration instead of
                the adaptive convergence protocol.
            equilibrate_writhe: If True, additionally equilibrate the writhe
                for circular DNA.  Requires ``simple=True``.
            endpoints_fixed: If True, the first and last base pairs are held
                fixed during the simulation.
            fixed: List of base-pair indices to keep fixed.
            dump_every: Save a trajectory snapshot every *n* MC sweeps.
            plot: If True, plot the energy trace during full equilibration.

        Raises:
            ValueError: If ``equilibrate_writhe=True`` with ``simple=False``.
        """
        # Set the parameters
        self.endpoints_fixed = endpoints_fixed
        self.fixed = fixed
        self.exvol_rad = exvol_rad
        self.temperature = temperature
        self.frame = frame
        print('Minimize the DNA structure:\nsimple equilibration =', simple, '\nequilibrate writhe =', equilibrate_writhe, '\nexcluded volume radius =', exvol_rad, '\ntemperature =', temperature)
        minimizer = self._initialize_mc_engine()    

        # Run the Monte Carlo simulation
        if equilibrate_writhe:
            self.out = minimizer.equilibrate_simple(equilibrate_writhe=equilibrate_writhe,dump_every=dump_every)
        elif not simple and equilibrate_writhe:
            raise ValueError("Minimization of writhe is only supported for simple equilibration.")
        else:
            self.out = minimizer.equilibrate(dump_every=dump_every,plot_equi=plot)

        # Update the reference frames
        positions, triads = self._get_positions_and_triads()
        self._update_frames()

    def get_MC_traj(self) -> 'md.Trajectory':
        """Build an MDTraj Trajectory from the MC sampling snapshots.

        Each frame contains Argon atoms at base-pair center positions and
        Helium (dummy) atoms offset along the major-groove vector, connected
        by bonds.  This allows visualization of the MC relaxation path.

        Returns:
            md.Trajectory: An MDTraj trajectory with ``2 * n_bp`` atoms per
                frame and ``n_snapshots`` frames.
        """
        # Get the xyz coordinates of the new spline
        xyz = self.out['confs'][:, :, :3, 3]

        # Get the triads and calculate the major vector positions
        triads = self.out['confs'][:, :, :3, :3].transpose(0, 1, 3, 2)
        major_vector = triads[:, :, 0, :]
        major_positions = xyz + major_vector*1.1  # Scale the major vector by 1.1

        # Concatenate the original xyz and the major_positions
        all_positions = np.concatenate((xyz, major_positions), axis=1)

        # Create a topology for the new spline
        topology = md.Topology()
        # Add a chain to the topology
        chain = topology.add_chain()
        # Add argon atoms to the topology
        num_atoms = xyz.shape[1]
        for i in range(num_atoms):
            residue = topology.add_residue(name='Ar', chain=chain)
            topology.add_atom('Ar', element=md.element.argon, residue=residue)

        # Add dummy atoms to the topology
        for i in range(num_atoms):
            residue = topology.add_residue(name='DUM', chain=chain)
            topology.add_atom('DUM', element=md.element.helium, residue=residue)  # Using helium as a placeholder element

        # Add bonds to the topology
        for i in range(num_atoms - 1):
            topology.add_bond(topology.atom(i), topology.atom(i + 1))

        # Add bonds between each atom and its corresponding dummy atom
        for i in range(num_atoms):
            topology.add_bond(topology.atom(i), topology.atom(num_atoms + i))

        if self.circular:
            # Add a bond between the first and last atom
            topology.add_bond(topology.atom(0), topology.atom(num_atoms - 1))

        # Create a trajectory from the all_positions coordinates and the topology
        traj = md.Trajectory(all_positions, topology=topology)

        return traj

    def run(self, cycles: int, dump_every: int = 20, start_id: int = 0) -> np.ndarray:
        """Run the Monte Carlo simulation"""
        raise NotImplementedError("This method is not implemented yet.")

get_MC_traj()

Build an MDTraj Trajectory from the MC sampling snapshots.

Each frame contains Argon atoms at base-pair center positions and Helium (dummy) atoms offset along the major-groove vector, connected by bonds. This allows visualization of the MC relaxation path.

Returns:

Type	Description
Trajectory	md.Trajectory: An MDTraj trajectory with 2 * n_bp atoms per frame and n_snapshots frames.

Source code in mdna/minimizer.py

def get_MC_traj(self) -> 'md.Trajectory':
    """Build an MDTraj Trajectory from the MC sampling snapshots.

    Each frame contains Argon atoms at base-pair center positions and
    Helium (dummy) atoms offset along the major-groove vector, connected
    by bonds.  This allows visualization of the MC relaxation path.

    Returns:
        md.Trajectory: An MDTraj trajectory with ``2 * n_bp`` atoms per
            frame and ``n_snapshots`` frames.
    """
    # Get the xyz coordinates of the new spline
    xyz = self.out['confs'][:, :, :3, 3]

    # Get the triads and calculate the major vector positions
    triads = self.out['confs'][:, :, :3, :3].transpose(0, 1, 3, 2)
    major_vector = triads[:, :, 0, :]
    major_positions = xyz + major_vector*1.1  # Scale the major vector by 1.1

    # Concatenate the original xyz and the major_positions
    all_positions = np.concatenate((xyz, major_positions), axis=1)

    # Create a topology for the new spline
    topology = md.Topology()
    # Add a chain to the topology
    chain = topology.add_chain()
    # Add argon atoms to the topology
    num_atoms = xyz.shape[1]
    for i in range(num_atoms):
        residue = topology.add_residue(name='Ar', chain=chain)
        topology.add_atom('Ar', element=md.element.argon, residue=residue)

    # Add dummy atoms to the topology
    for i in range(num_atoms):
        residue = topology.add_residue(name='DUM', chain=chain)
        topology.add_atom('DUM', element=md.element.helium, residue=residue)  # Using helium as a placeholder element

    # Add bonds to the topology
    for i in range(num_atoms - 1):
        topology.add_bond(topology.atom(i), topology.atom(i + 1))

    # Add bonds between each atom and its corresponding dummy atom
    for i in range(num_atoms):
        topology.add_bond(topology.atom(i), topology.atom(num_atoms + i))

    if self.circular:
        # Add a bond between the first and last atom
        topology.add_bond(topology.atom(0), topology.atom(num_atoms - 1))

    # Create a trajectory from the all_positions coordinates and the topology
    traj = md.Trajectory(all_positions, topology=topology)

    return traj

minimize(frame=-1, exvol_rad=2.0, temperature=300, simple=False, equilibrate_writhe=False, endpoints_fixed=True, fixed=[], dump_every=20, plot=False)

Run Monte Carlo energy minimization on the DNA structure.

Initializes the PMCpy engine with the current base-pair frames and sequence, performs MC equilibration, and updates the internal frames with the relaxed configuration.

Parameters:

Name	Type	Description	Default
frame	int	Index of the stored frame to minimize (default -1, the last frame).	-1
exvol_rad	float	Excluded-volume radius in nm. Base pairs closer than this distance incur a repulsive energy penalty.	2.0
temperature	int	Temperature in Kelvin for the Metropolis acceptance criterion.	300
simple	bool	If True, use simple (fixed-sweep) equilibration instead of the adaptive convergence protocol.	False
equilibrate_writhe	bool	If True, additionally equilibrate the writhe for circular DNA. Requires simple=True.	False
endpoints_fixed	bool	If True, the first and last base pairs are held fixed during the simulation.	True
fixed	List[int]	List of base-pair indices to keep fixed.	[]
dump_every	int	Save a trajectory snapshot every n MC sweeps.	20
plot	bool	If True, plot the energy trace during full equilibration.	False

Raises:

Type	Description
ValueError	If equilibrate_writhe=True with simple=False.

Source code in mdna/minimizer.py

def minimize(self,  frame: int = -1, exvol_rad : float = 2.0, temperature : int = 300,  simple : bool = False, equilibrate_writhe : bool = False, endpoints_fixed : bool = True, fixed : List[int] = [], dump_every : int = 20, plot : bool = False) -> None:
    """Run Monte Carlo energy minimization on the DNA structure.

    Initializes the PMCpy engine with the current base-pair frames and
    sequence, performs MC equilibration, and updates the internal frames
    with the relaxed configuration.

    Args:
        frame: Index of the stored frame to minimize (default ``-1``, the
            last frame).
        exvol_rad: Excluded-volume radius in nm.  Base pairs closer than
            this distance incur a repulsive energy penalty.
        temperature: Temperature in Kelvin for the Metropolis acceptance
            criterion.
        simple: If True, use simple (fixed-sweep) equilibration instead of
            the adaptive convergence protocol.
        equilibrate_writhe: If True, additionally equilibrate the writhe
            for circular DNA.  Requires ``simple=True``.
        endpoints_fixed: If True, the first and last base pairs are held
            fixed during the simulation.
        fixed: List of base-pair indices to keep fixed.
        dump_every: Save a trajectory snapshot every *n* MC sweeps.
        plot: If True, plot the energy trace during full equilibration.

    Raises:
        ValueError: If ``equilibrate_writhe=True`` with ``simple=False``.
    """
    # Set the parameters
    self.endpoints_fixed = endpoints_fixed
    self.fixed = fixed
    self.exvol_rad = exvol_rad
    self.temperature = temperature
    self.frame = frame
    print('Minimize the DNA structure:\nsimple equilibration =', simple, '\nequilibrate writhe =', equilibrate_writhe, '\nexcluded volume radius =', exvol_rad, '\ntemperature =', temperature)
    minimizer = self._initialize_mc_engine()    

    # Run the Monte Carlo simulation
    if equilibrate_writhe:
        self.out = minimizer.equilibrate_simple(equilibrate_writhe=equilibrate_writhe,dump_every=dump_every)
    elif not simple and equilibrate_writhe:
        raise ValueError("Minimization of writhe is only supported for simple equilibration.")
    else:
        self.out = minimizer.equilibrate(dump_every=dump_every,plot_equi=plot)

    # Update the reference frames
    positions, triads = self._get_positions_and_triads()
    self._update_frames()

run(cycles, dump_every=20, start_id=0)

Run the Monte Carlo simulation

Source code in mdna/minimizer.py

def run(self, cycles: int, dump_every: int = 20, start_id: int = 0) -> np.ndarray:
    """Run the Monte Carlo simulation"""
    raise NotImplementedError("This method is not implemented yet.")