openff/interchange/drivers/openmm.py
def _get_openmm_energies(
    omm_sys: openmm.System,
    box_vectors,
    positions,
    round_positions=None,
    hard_cutoff=False,
    electrostatics: bool = True,
) -> EnergyReport:
    """Given a prepared `openmm.System`, run a single-point energy calculation."""
    """\
    if hard_cutoff:
        omm_sys = _set_nonbonded_method(
            omm_sys, "cutoff", electrostatics=electrostatics
        )
    else:
        omm_sys = _set_nonbonded_method(omm_sys, "PME")
    """

    for idx, force in enumerate(omm_sys.getForces()):
        force.setForceGroup(idx)

    integrator = openmm.VerletIntegrator(1.0 * unit.femtoseconds)
    context = openmm.Context(omm_sys, integrator)

    if box_vectors is not None:
        if not isinstance(box_vectors, (unit.Quantity, list)):
            box_vectors = box_vectors.magnitude * unit.nanometer
        context.setPeriodicBoxVectors(*box_vectors)

    if isinstance(positions, unit.Quantity):
        # Convert list of Vec3 into a NumPy array
        positions = np.asarray(positions.value_in_unit(unit.nanometer)) * unit.nanometer
    else:
        positions = positions.magnitude * unit.nanometer

    if round_positions is not None:
        rounded = np.round(positions, round_positions)
        context.setPositions(rounded)
    else:
        context.setPositions(positions)

    raw_energies = dict()
    omm_energies = dict()

    for idx in range(omm_sys.getNumForces()):
        state = context.getState(getEnergy=True, groups={idx})
        raw_energies[idx] = state.getPotentialEnergy()
        del state

    # This assumes that only custom forces will have duplicate instances
    for key in raw_energies:
        force = omm_sys.getForce(key)
        if type(force) == openmm.HarmonicBondForce:
            omm_energies["HarmonicBondForce"] = raw_energies[key]
        elif type(force) == openmm.HarmonicAngleForce:
            omm_energies["HarmonicAngleForce"] = raw_energies[key]
        elif type(force) == openmm.PeriodicTorsionForce:
            omm_energies["PeriodicTorsionForce"] = raw_energies[key]
        elif type(force) in [
            openmm.NonbondedForce,
            openmm.CustomNonbondedForce,
            openmm.CustomBondForce,
        ]:
            if "Nonbonded" in omm_energies:
                omm_energies["Nonbonded"] += raw_energies[key]
            else:
                omm_energies["Nonbonded"] = raw_energies[key]

    # Fill in missing keys if interchange does not have all typical forces
    for required_key in [
        "HarmonicBondForce",
        "HarmonicAngleForce",
        "PeriodicTorsionForce",
        "NonbondedForce",
    ]:
        if not any(required_key in val for val in omm_energies):
            pass  # omm_energies[required_key] = 0.0 * kj_mol

    del context
    del integrator

    report = EnergyReport()

    report.update_energies(
        {
            "Bond": omm_energies.get("HarmonicBondForce", 0.0 * kj_mol),
            "Angle": omm_energies.get("HarmonicAngleForce", 0.0 * kj_mol),
            "Torsion": _canonicalize_torsion_energies(omm_energies),
            "Nonbonded": omm_energies.get(
                "Nonbonded", _canonicalize_nonbonded_energies(omm_energies)
            ),
        }
    )

    report.energies.pop("vdW")
    report.energies.pop("Electrostatics")

    return report
Given a prepared `openmm.System`, run a single-point energy calculation.
Given a prepared `openmm.System`, run a single-point energy calculation.
Given a prepared ` openmm . System ` run a single - point energy calculation .
def _get_openmm_energies(
    omm_sys: openmm.System,
    box_vectors,
    positions,
    round_positions=None,
    hard_cutoff=False,
    electrostatics: bool = True,
) -> EnergyReport:
    for idx, force in enumerate(omm_sys.getForces()):
        force.setForceGroup(idx)
    integrator = openmm.VerletIntegrator(1.0 * unit.femtoseconds)
    context = openmm.Context(omm_sys, integrator)
    if box_vectors is not None:
        if not isinstance(box_vectors, (unit.Quantity, list)):
            box_vectors = box_vectors.magnitude * unit.nanometer
        context.setPeriodicBoxVectors(*box_vectors)
    if isinstance(positions, unit.Quantity):
        positions = np.asarray(positions.value_in_unit(unit.nanometer)) * unit.nanometer
    else:
        positions = positions.magnitude * unit.nanometer
    if round_positions is not None:
        rounded = np.round(positions, round_positions)
        context.setPositions(rounded)
    else:
        context.setPositions(positions)
    raw_energies = dict()
    omm_energies = dict()
    for idx in range(omm_sys.getNumForces()):
        state = context.getState(getEnergy=True, groups={idx})
        raw_energies[idx] = state.getPotentialEnergy()
        del state
    for key in raw_energies:
        force = omm_sys.getForce(key)
        if type(force) == openmm.HarmonicBondForce:
            omm_energies["HarmonicBondForce"] = raw_energies[key]
        elif type(force) == openmm.HarmonicAngleForce:
            omm_energies["HarmonicAngleForce"] = raw_energies[key]
        elif type(force) == openmm.PeriodicTorsionForce:
            omm_energies["PeriodicTorsionForce"] = raw_energies[key]
        elif type(force) in [
            openmm.NonbondedForce,
            openmm.CustomNonbondedForce,
            openmm.CustomBondForce,
        ]:
            if "Nonbonded" in omm_energies:
                omm_energies["Nonbonded"] += raw_energies[key]
            else:
                omm_energies["Nonbonded"] = raw_energies[key]
    for required_key in [
        "HarmonicBondForce",
        "HarmonicAngleForce",
        "PeriodicTorsionForce",
        "NonbondedForce",
    ]:
        if not any(required_key in val for val in omm_energies):
            pass  
    del context
    del integrator
    report = EnergyReport()
    report.update_energies(
        {
            "Bond": omm_energies.get("HarmonicBondForce", 0.0 * kj_mol),
            "Angle": omm_energies.get("HarmonicAngleForce", 0.0 * kj_mol),
            "Torsion": _canonicalize_torsion_energies(omm_energies),
            "Nonbonded": omm_energies.get(
                "Nonbonded", _canonicalize_nonbonded_energies(omm_energies)
            ),
        }
    )
    report.energies.pop("vdW")
    report.energies.pop("Electrostatics")
    return report
def _get_openmm_energies ( omm_sys : openmm . System , box_vectors , positions , round_positions = None , hard_cutoff = False , electrostatics : bool = True , ) -> EnergyReport : """\
    if hard_cutoff:
        omm_sys = _set_nonbonded_method(
            omm_sys, "cutoff", electrostatics=electrostatics
        )
    else:
        omm_sys = _set_nonbonded_method(omm_sys, "PME")
    """ for idx , force in enumerate ( omm_sys . getForces ( ) ) : force . setForceGroup ( idx ) integrator = openmm . VerletIntegrator ( 1.0 * unit . femtoseconds ) context = openmm . Context ( omm_sys , integrator ) if box_vectors is not None : if not isinstance ( box_vectors , ( unit . Quantity , list ) ) : box_vectors = box_vectors . magnitude * unit . nanometer context . setPeriodicBoxVectors ( * box_vectors ) if isinstance ( positions , unit . Quantity ) : positions = np . asarray ( positions . value_in_unit ( unit . nanometer ) ) * unit . nanometer else : positions = positions . magnitude * unit . nanometer if round_positions is not None : rounded = np . round ( positions , round_positions ) context . setPositions ( rounded ) else : context . setPositions ( positions ) raw_energies = dict ( ) omm_energies = dict ( ) for idx in range ( omm_sys . getNumForces ( ) ) : state = context . getState ( getEnergy = True , groups = { idx } ) raw_energies [ idx ] = state . getPotentialEnergy ( ) del state for key in raw_energies : force = omm_sys . getForce ( key ) if type ( force ) == openmm . HarmonicBondForce : omm_energies [ "HarmonicBondForce" ] = raw_energies [ key ] elif type ( force ) == openmm . HarmonicAngleForce : omm_energies [ "HarmonicAngleForce" ] = raw_energies [ key ] elif type ( force ) == openmm . PeriodicTorsionForce : omm_energies [ "PeriodicTorsionForce" ] = raw_energies [ key ] elif type ( force ) in [ openmm . NonbondedForce , openmm . CustomNonbondedForce , openmm . CustomBondForce , ] : if "Nonbonded" in omm_energies : omm_energies [ "Nonbonded" ] += raw_energies [ key ] else : omm_energies [ "Nonbonded" ] = raw_energies [ key ] for required_key in [ "HarmonicBondForce" , "HarmonicAngleForce" , "PeriodicTorsionForce" , "NonbondedForce" , ] : if not any ( required_key in val for val in omm_energies ) : pass del context del integrator report = EnergyReport ( ) report . update_energies ( { "Bond" : omm_energies . get ( "HarmonicBondForce" , 0.0 * kj_mol ) , "Angle" : omm_energies . get ( "HarmonicAngleForce" , 0.0 * kj_mol ) , "Torsion" : _canonicalize_torsion_energies ( omm_energies ) , "Nonbonded" : omm_energies . get ( "Nonbonded" , _canonicalize_nonbonded_energies ( omm_energies ) ) , } ) report . energies . pop ( "vdW" ) report . energies . pop ( "Electrostatics" ) return report
Given a prepared `openmm.System`, run a single-point energy calculation.
Given a prepared ` openmm . System ` run a single - point energy calculation .
"""Given a prepared `openmm.System`, run a single-point energy calculation.""" """\
    if hard_cutoff:
        omm_sys = _set_nonbonded_method(
            omm_sys, "cutoff", electrostatics=electrostatics
        )
    else:
        omm_sys = _set_nonbonded_method(omm_sys, "PME")
    """ # Convert list of Vec3 into a NumPy array # This assumes that only custom forces will have duplicate instances # Fill in missing keys if interchange does not have all typical forces # omm_energies[required_key] = 0.0 * kj_mol
omm_sys
openmm.System
box_vectors
positions
round_positions
hard_cutoff
electrostatics
bool