Skip to content

Sync with steven #27

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ftclark3 wants to merge 15 commits into
base: master
Choose a base branch
from
Open

Sync with steven #27

Show file tree
Hide file tree
Changes from 10 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
7b1ee50
allow user to specify debye length and nonbonded cutoff for dna.Elect…
ftclark3 May 27, 2026
ded7e97
allow user to specify nonbonded cutoff, debye length, and excluded vo…
ftclark3 May 27, 2026
06b0184
revert to former default nonbonded cutoff in dna.Electrostatics
ftclark3 May 27, 2026
f8d3170
clarified comments
ftclark3 May 28, 2026
834069e
allow user to specify debye length and nonbonded cutoff
ftclark3 May 28, 2026
7ea49ff
added steven's dna force classes
ftclark3 May 28, 2026
7e7b881
updated StringProteinDNA with steven's logical changes commented out.…
ftclark3 May 28, 2026
90c9a57
condense protein distance distance bias classes
ftclark3 May 28, 2026
5f1e763
added protein-dna biasing functions based on electrostatics and the p…
ftclark3 May 28, 2026
253318c
added steven's scripts
ftclark3 May 28, 2026
5c51e84
Update cutoff_distance to use nanometer units
ftclark3 May 28, 2026
697ae35
Update protein_dna.py
ftclark3 May 28, 2026
68b0c88
setting up new branch
ftclark3 Jun 1, 2026
7c8e235
Merge remote-tracking branch 'ftclark3/sync_with_steven' into sync_wi…
ftclark3 Jun 1, 2026
17468e7
updated imports
ftclark3 Jun 1, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
147 changes: 142 additions & 5 deletions open3SPN2/force/dna.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -487,10 +487,12 @@ def defineInteraction(self):


class Electrostatics(DNAForce, openmm.CustomNonbondedForce):
def __init__(self, dna, force_group=13, temperature=300*unit.kelvin, salt_concentration=100*unit.millimolar, OpenCLPatch=True):
def __init__(self, dna, force_group=13, temperature=300*unit.kelvin, salt_concentration=100*unit.millimolar, ldby = None, cutoff_distance = None, OpenCLPatch=True):
self.force_group = force_group
self.T = temperature
self.C = salt_concentration
self.ldby = ldby
self.cutoff_distance = cutoff_distance
super().__init__(dna, OpenCLPatch=OpenCLPatch)

def reset(self):
Expand All @@ -506,19 +508,33 @@ def reset(self):
ec = 1.60217653E-19 * unit.coulomb # proton charge
pv = 8.8541878176E-12 * unit.farad / unit.meter # dielectric permittivity of vacuum

ldby = np.sqrt(dielectric * pv * kb * T / (2.0 * Na * ec ** 2 * C))
if self.ldby is None:
ldby = np.sqrt(dielectric * pv * kb * T / (2.0 * Na * ec ** 2 * C))
else:
ldby = self.ldby

ldby = ldby.in_units_of(unit.nanometer)

if self.cutoff_distance == None:
cutoff_distance = 5 # for backward compatibility
else:
cutoff_distance = self.cutoff_distance

cutoff_nm = cutoff_distance.value_in_unit(unit.nanometer)
Comment thread
sourcery-ai[bot] marked this conversation as resolved.

denominator = 4 * np.pi * pv * dielectric / (Na * ec ** 2)
denominator = denominator.in_units_of(unit.kilocalorie_per_mole**-1 * unit.nanometer**-1)
#print(ldby, denominator)
print(ldby, denominator)

electrostaticForce = openmm.CustomNonbondedForce("""energy;
energy=q1*q2*exp(-r/dh_length)/denominator/r;""")
electrostaticForce.addPerParticleParameter('q')
electrostaticForce.addGlobalParameter('dh_length', ldby)
electrostaticForce.addGlobalParameter('denominator', denominator)

electrostaticForce.setCutoffDistance(5)
electrostaticForce.setCutoffDistance(cutoff_nm)
print(f"dna screening length {ldby} nm")
print(f"dna electrostatic cutoff {electrostaticForce.getCutoffDistance()} nm")
if self.periodic:
electrostaticForce.setNonbondedMethod(electrostaticForce.CutoffPeriodic)
else:
Expand Down Expand Up @@ -546,4 +562,125 @@ def defineInteraction(self):
self.force.addParticle(parameters)

# add neighbor exclusion
addNonBondedExclusions(self.dna, self.force, self.OpenCLPatch)
addNonBondedExclusions(self.dna, self.force, self.OpenCLPatch)


class group_constraint_by_position(DNAForce):
"""
Apply a harmonic restraint on the center of mass of selected DNA atoms to keep them near (x0, y0, z0).
"""

def __init__(self, dna, k=1*unit.kilocalorie_per_mole,
x0=10*unit.angstrom, y0=10*unit.angstrom, z0=10*unit.angstrom,
appliedToResidues=None, force_group=24, OpenCLPatch=True):

self.force_group = force_group
Comment on lines +568 to +577

@sourcery-ai sourcery-ai Bot May 28, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

suggestion (bug_risk): Constraint uses sum of mass-weighted coordinates, not center of mass as advertised

The current CV uses sum(m_i x_i) directly, so the effective reference position scales with the number and masses of atoms instead of being invariant like a true COM (which would be sum(m_i x_i) / sum(m_i)). To match the documented behavior, either:

  • divide each coordinate sum by total_mass in the CustomCVForce (e.g., (sum_x/total_mass - x0)^2), or
  • normalize the per-particle weights to mass/total_mass and keep the existing expression.

Aligning the implementation with the docstring will prevent downstream confusion and unintended parameter changes.

Suggested implementation:

class group_constraint_by_position(DNAForce):
    """
    Apply a harmonic restraint on the center of mass of selected DNA atoms to keep them near (x0, y0, z0).

    The actual CV is the center-of-mass (COM) position:
        COM = sum_i(m_i * r_i) / sum_i(m_i)
    not just the raw mass-weighted position sum_i(m_i * r_i).
    """

    def __init__(self, dna, k=1*unit.kilocalorie_per_mole, 
                 x0=10*unit.angstrom, y0=10*unit.angstrom, z0=10*unit.angstrom, 
                 appliedToResidues=None, force_group=24, OpenCLPatch=True):

        # Store configuration
        self.dna = dna
        self.force_group = force_group
        self.k = k
        self.x0 = x0
        self.y0 = y0
        self.z0 = z0
        self.appliedToResidues = appliedToResidues
        self.OpenCLPatch = OpenCLPatch

        # Build the OpenMM force object
        self._create_force()

    def _create_force(self):
        """
        Create a CustomCVForce that restrains the center of mass of the selected atoms.
        """
        import openmm as mm
        from openmm import unit as omm_unit

        system = self.dna.system
        topology = self.dna.topology

        # Collect the particle indices that belong to the selected residues
        selected_particle_indices = []
        if self.appliedToResidues is None:
            # Default: all DNA atoms; assumes dna.atom_list or equivalent is available.
            # Adapt this selection logic to your existing DNAForce conventions.
            for atom in topology.atoms():
                if atom.residue.name.startswith("D") or atom.residue.name.startswith("R"):
                    selected_particle_indices.append(atom.index)
        else:
            residues_set = set(self.appliedToResidues)
            for atom in topology.atoms():
                if atom.residue.index in residues_set:
                    selected_particle_indices.append(atom.index)

        if len(selected_particle_indices) == 0:
            raise ValueError("group_constraint_by_position: no atoms selected for COM restraint.")

        # Compute total mass of the group
        masses = [system.getParticleMass(i) for i in selected_particle_indices]
        total_mass = sum(masses)  # this is an OpenMM unit-bearing quantity

        # Create per-coordinate CVs: sum(m_i * x_i), sum(m_i * y_i), sum(m_i * z_i)
        # These are still mass-weighted sums; we convert to COM by dividing by total_mass in the main CV expression.
        sum_x_force = mm.CustomCVForce("sum(mass * x)")
        sum_y_force = mm.CustomCVForce("sum(mass * y)")
        sum_z_force = mm.CustomCVForce("sum(mass * z)")

        sum_x_force.addPerParticleParameter("mass")
        sum_y_force.addPerParticleParameter("mass")
        sum_z_force.addPerParticleParameter("mass")

        for idx, m in zip(selected_particle_indices, masses):
            sum_x_force.addParticle(idx, [m])
            sum_y_force.addParticle(idx, [m])
            sum_z_force.addParticle(idx, [m])

        # Main COM restraint CV:
        # Use COM = sum(m_i * coord_i) / total_mass for each coordinate
        com_restraint = mm.CustomCVForce(
            "0.5 * k * ( (sum_x/total_mass - x0)^2 + (sum_y/total_mass - y0)^2 + (sum_z/total_mass - z0)^2 )"
        )

        # Add the collective variables that provide sum_x, sum_y, sum_z
        com_restraint.addCollectiveVariable("sum_x", sum_x_force)
        com_restraint.addCollectiveVariable("sum_y", sum_y_force)
        com_restraint.addCollectiveVariable("sum_z", sum_z_force)

        # Global parameters: force constant, reference COM position, and total mass of the group
        # Note: total_mass is constant for this group; exposing it as a global parameter
        # makes the COM definition explicit and invariant to group size.
        com_restraint.addGlobalParameter("k", self.k)
        com_restraint.addGlobalParameter("x0", self.x0)
        com_restraint.addGlobalParameter("y0", self.y0)
        com_restraint.addGlobalParameter("z0", self.z0)
        com_restraint.addGlobalParameter("total_mass", total_mass)

        com_restraint.setForceGroup(self.force_group)

        # Register this force with the DNAForce infrastructure
        self.force = com_restraint
        self.addForce(self.force, self.OpenCLPatch)
  1. The selection logic for selected_particle_indices must be aligned with how other DNAForce subclasses in this file choose their atoms (e.g., there might already be helper methods or precomputed atom lists you should reuse instead of the placeholder topology-based loop).
  2. The current implementation of group_constraint_by_position in the PR likely already constructs CustomCVForce instances for sum(m_i * x_i) etc.; if so, you can instead only change:
    • the main CV expression from (sum_x - x0)^2 to (sum_x/total_mass - x0)^2 (and similarly for y, z), and
    • add a global parameter total_mass initialized to the Python-computed total mass for the group.
  3. Ensure that any existing calls to addForce or class initialization in DNAForce are preserved; if DNAForce expects a different initialization pattern (e.g., overriding setup() rather than calling _create_force() in __init__), adapt _create_force() to that pattern instead of calling addForce directly as shown.

self.k = k
self.x0 = x0
self.y0 = y0
self.z0 = z0
self.appliedToResidues = appliedToResidues

# Call the superclass initializer
super().__init__(dna, OpenCLPatch=OpenCLPatch)

def reset(self):
# Convert parameters
k_constraint = self.k.value_in_unit(unit.kilojoule_per_mole)
x0 = self.x0.value_in_unit(unit.nanometer)
y0 = self.y0.value_in_unit(unit.nanometer)
z0 = self.z0.value_in_unit(unit.nanometer)

# Define mass-weighted sum forces
sum_x = openmm.CustomExternalForce("mass * x")
sum_y = openmm.CustomExternalForce("mass * y")
sum_z = openmm.CustomExternalForce("mass * z")

sum_x.addPerParticleParameter("mass")
sum_y.addPerParticleParameter("mass")
sum_z.addPerParticleParameter("mass")

# Define harmonic restraint on COM
harmonic = openmm.CustomCVForce(
f"{k_constraint} * ((sum_x - {x0})^2 + (sum_y - {y0})^2 + (sum_z - {z0})^2)"
)
harmonic.addCollectiveVariable("sum_x", sum_x)
harmonic.addCollectiveVariable("sum_y", sum_y)
harmonic.addCollectiveVariable("sum_z", sum_z)
harmonic.setForceGroup(self.force_group)

self.force = harmonic

def defineInteraction(self):
"""
Adds selected DNA atoms to the sum_x, sum_y, sum_z forces with appropriate mass weighting.
"""
total_mass = 0.0
for i in range(len(self.dna.atoms)):
#mass = self.dna.system.getParticleMass(i).value_in_unit(unit.dalton)
mass = 1

# Add particle to sum_x, sum_y, sum_z
Comment on lines +614 to +623

@sourcery-ai sourcery-ai Bot May 28, 2026

Copy link
Copy Markdown

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

suggestion (bug_risk): appliedToResidues and actual masses are ignored in group_constraint_by_position

In defineInteraction, all DNA atoms are added and mass = 1 is used for every particle (with the getParticleMass call commented out). If residue-based selection is meant to matter, the indices should be filtered using appliedToResidues before adding them to the CVs. Also, using unit mass for all atoms turns the CV into an unweighted average position rather than a true COM. Please either restore getParticleMass (or equivalent) or clearly document that this CV is intentionally unweighted, and ensure appliedToResidues is respected when building the index list.

Suggested implementation:

        self.force = harmonic

    def defineInteraction(self):
        """
        Adds selected DNA atoms to the sum_x, sum_y, sum_z collective variables with mass weighting
        to compute a true center-of-mass (COM). If ``appliedToResidues`` is provided, only atoms
        belonging to those residues are included.
        """
        total_mass = 0.0

        # Build the list of atom indices to include, honoring appliedToResidues if provided.
        if self.appliedToResidues is not None:
            atom_indices = []
            for i, atom in enumerate(self.dna.atoms):
                # Assumes atoms expose a residue identifier as `resid`.
                # If your topology uses a different attribute, update this accordingly.
                resid = getattr(atom, "resid", None)
                if resid in self.appliedToResidues:
                    atom_indices.append(i)
        else:
            atom_indices = range(len(self.dna.atoms))

        for i in atom_indices:
            # Use the actual particle mass so the CV represents a true COM, not an unweighted average.
            mass = self.dna.system.getParticleMass(i).value_in_unit(unit.dalton)

            # Add particle to sum_x, sum_y, sum_z with its mass as the weight
            self.force.getCollectiveVariable(0).addParticle(i, [mass])
            self.force.getCollectiveVariable(1).addParticle(i, [mass])
            self.force.getCollectiveVariable(2).addParticle(i, [mass])

            total_mass += mass

        if total_mass == 0:
            raise ValueError("No atoms were selected for the group constraint; check appliedToResidues.")
  1. Ensure that the atom objects in self.dna.atoms expose the residue identifier as atom.resid. If a different attribute is used (e.g. atom.residue.id or atom.residue.index), adjust the resid = getattr(atom, "resid", None) line and the subsequent membership test accordingly.
  2. Confirm that unit is imported in this module (e.g. from openmm import unit or import openmm.unit as unit) so that value_in_unit(unit.dalton) works as expected.
  3. If the codebase already has a helper for residue-based atom selection (e.g. on DNAForce or dna), you may prefer to replace the explicit loop that builds atom_indices with that helper to keep residue-selection logic consistent across forces.

self.force.getCollectiveVariable(0).addParticle(i, [mass])
self.force.getCollectiveVariable(1).addParticle(i, [mass])
self.force.getCollectiveVariable(2).addParticle(i, [mass])

total_mass += mass

if total_mass == 0:
raise ValueError("No atoms were selected for the group constraint; check appliedToResidues.")


class measure_from_position(DNAForce):
def __init__(self, dna, x0=10*unit.angstrom, y0=10*unit.angstrom, z0=10*unit.angstrom, appliedToResidues=None, force_group=4):
self.force_group = force_group
self.x0 = x0
self.y0 = y0
self.z0 = z0
self.appliedToResidues = appliedToResidues
super().__init__(dna)

def reset(self):
# Convert parameters
x0 = self.x0.value_in_unit(unit.nanometer)
y0 = self.y0.value_in_unit(unit.nanometer)
z0 = self.z0.value_in_unit(unit.nanometer)

# Define mass-weighted sum forces
sum_x = openmm.CustomExternalForce("mass * x")
sum_y = openmm.CustomExternalForce("mass * y")
sum_z = openmm.CustomExternalForce("mass * z")

sum_x.addPerParticleParameter("mass")
sum_y.addPerParticleParameter("mass")
sum_z.addPerParticleParameter("mass")

# Define harmonic restraint on COM
harmonic = openmm.CustomCVForce(
f"(sum_x - {x0})^2 + (sum_y - {y0})^2 + (sum_z - {z0})^2"
)
harmonic.addCollectiveVariable("sum_x", sum_x)
harmonic.addCollectiveVariable("sum_y", sum_y)
harmonic.addCollectiveVariable("sum_z", sum_z)
harmonic.setForceGroup(self.force_group)

self.force = harmonic

def defineInteraction(self):
"""
Adds selected DNA atoms to the sum_x, sum_y, sum_z forces with appropriate mass weighting.
"""
total_mass = 0.0
for i in range(len(self.dna.atoms)):
#mass = self.dna.system.getParticleMass(i).value_in_unit(unit.dalton)\
mass = 1

# Add particle to sum_x, sum_y, sum_z
self.force.getCollectiveVariable(0).addParticle(i, [mass])
self.force.getCollectiveVariable(1).addParticle(i, [mass])
self.force.getCollectiveVariable(2).addParticle(i, [mass])

total_mass += mass

if total_mass == 0:
raise ValueError("No atoms were selected for the group constraint; check appliedToResidues.")
Loading