Repository: UC-Davis-molecular-computing/nuad (language: python)

Issue:
refactor constraints to reuse code
Currently there is nearly identical code for finding violations of `DomainConstraint`'s, `StrandConstraint`'s, `DomainPairConstraint`'s, and `StrandPairConstraint`'s.

For implementing other issues such as #92, it will be easier to modify if we refactor to share code between these.

One way is for the constraints to have an identical function signature for `evaluate`, e.g., each takes a tuple of strings for sequence(s) and a `DesignPart`, which is a union type of `Domain | Strand | Tuple[Domain, Domain] | Tuple[Strand, Strand]`, for example (or perhaps the sequence can be a union of `str` and `Tuple[str, str]` to be consistent between them, or conversely the second argument could be a `Tuple[DesignPart]`, where `DesignPart` is a union of `Domain` and `Strand`; this latter would also capture `ComplexConstraint`).
move constraint evaluate and summary into one function, so that creating reports doesn't require calling constraints again
Currently, we separate the evaluate of a constraint, and the summarizing of it in a report, into two steps. This is inefficient (though mitigated somewhat by lru_cached calls to the underlying NUPACK functions) since creating a report requires re-evaluating every constraint, in essence. (Most constraint `summary` calls do essentially the same thing as `evaluate`.)

By combining `evaluate` and `summary` into a single function, the `Violation` created when a constraint is violated can store the summary, and then creating a report will be much faster, just asking the `Violation` for the summary.

This might lock us into only reporting violations. Perhaps we can also keep track of non-violations in a different class.
move constraint evaluate and summary into one function, so that creating reports doesn't require calling constraints again
Currently, we separate the evaluate of a constraint, and the summarizing of it in a report, into two steps. This is inefficient (though mitigated somewhat by lru_cached calls to the underlying NUPACK functions) since creating a report requires re-evaluating every constraint, in essence. (Most constraint `summary` calls do essentially the same thing as `evaluate`.)

By combining `evaluate` and `summary` into a single function, the `Violation` created when a constraint is violated can store the summary, and then creating a report will be much faster, just asking the `Violation` for the summary.

This might lock us into only reporting violations. Perhaps we can also keep track of non-violations in a different class.


Discussion / hints:

Fix patch:
diff --git a/dsd/constraints.py b/dsd/constraints.py
index 373f2b00..2ae625ad 100644
--- a/dsd/constraints.py
+++ b/dsd/constraints.py
@@ -84,9 +84,9 @@
 T = TypeVar('T')
 KeyFunction = Callable[[T], Any]
 
-Complex = Tuple['Strand', ...]
-"""A Complex is a group of :any:`Strand`'s, in general that we expect to be bound by complementary 
-:any:`Domain`'s."""
+# Complex = Tuple['Strand', ...]
+# """A Complex is a group of :any:`Strand`'s, in general that we expect to be bound by complementary
+# :any:`Domain`'s."""
 
 all_dna_bases: Set[str] = {'A', 'C', 'G', 'T'}
 """
@@ -121,8 +121,6 @@ def _configure_logger() -> None:
 
 _configure_logger()
 
-T = TypeVar('T')
-
 
 def all_pairs(values: Iterable[T],
               with_replacement: bool = True,
@@ -969,11 +967,34 @@ def mandatory_field(ret_type: Type, json_map: Dict, main_key: str, *legacy_keys:
     raise ValueError(msg)
 
 
+class Part(ABC):
+
+    @abstractmethod
+    def name_of_part_type(self) -> str:
+        pass
+
+
 DomainLabel = TypeVar('DomainLabel')
 
 
 @dataclass
-class Domain(JSONSerializable, Generic[DomainLabel]):
+class DomainPair(Part, Generic[DomainLabel]):
+    domain1: Domain
+    domain2: Domain
+
+    @property
+    def name(self) -> str:
+        return f'{self.domain1.name}, {self.domain2.name}'
+
+    def name_of_part_type(self) -> str:
+        return 'domain pair'
+
+    def individual_parts(self) -> Tuple[Domain, Domain]:
+        return self.domain1, self.domain2
+
+
+@dataclass
+class Domain(JSONSerializable, Part, Generic[DomainLabel]):
     """
     Represents a contiguous substring of the DNA sequence of a :any:`Strand`, which is intended
     to be either single-stranded, or to bind fully to the Watson-Crick complement of the :any:`Domain`.
@@ -1001,10 +1022,6 @@ class Domain(JSONSerializable, Generic[DomainLabel]):
     length and individual :any:`DomainConstraint`'s.
     """
 
-    # TODO: Change name to "_sequence", add init=False
-    # TODO: `set_sequence`
-    #        - if fixed, error
-    #        - if dependent, error
     # TODO: `set_sequence_recursive_up`
 
     #        - if parent is not none, make recursive call to set_sequence_recursive_up
@@ -1077,9 +1094,6 @@ def __init__(self, name: str, pool: Optional[DomainPool] = None, sequence: Optio
         self.dependent = dependent
         self._subdomains = subdomains
 
-        self.__post_init__()
-
-    def __post_init__(self) -> None:
         if self.name.endswith('*'):
             raise ValueError('Domain name cannot end with *\n'
                              f'domain name = {self.name}')
@@ -1101,6 +1115,9 @@ def __post_init__(self) -> None:
         for subdomain in self._subdomains:
             subdomain.parent = self
 
+    def name_of_part_type(self) -> str:
+        return 'domain'
+
     def to_json_serializable(self, suppress_indent: bool = True) -> Union[NoIndent, Dict[str, Any]]:
         """
         :return:
@@ -1237,7 +1254,6 @@ def length(self) -> int:
                              'Assign a DomainPool (which has a length field) to give this Domain a length.')
         return self._pool.length
 
-    @property
     def sequence(self) -> str:
         """
         :return: DNA sequence of this domain (unstarred version)
@@ -1248,8 +1264,7 @@ def sequence(self) -> str:
                              f'sequence: {self._sequence}')
         return self._sequence
 
-    @sequence.setter
-    def sequence(self, new_sequence: str) -> None:
+    def set_sequence(self, new_sequence: str) -> None:
         """
         :param new_sequence: new DNA sequence to set
         """
@@ -1269,7 +1284,7 @@ def sequence(self, new_sequence: str) -> None:
                     f'Domain {self} is length {self.length} but subdomains {self._subdomains} has total '
                     f'length of {sd_total_length}')
         self._sequence = new_sequence
-        self._starred_sequence = dv.wc(self.sequence)
+        self._starred_sequence = dv.wc(new_sequence)
         self._set_subdomain_sequences(new_sequence)
         self._set_parent_sequence(new_sequence)
 
@@ -1321,7 +1336,7 @@ def set_fixed_sequence(self, fixed_sequence: str) -> None:
         :param fixed_sequence: new fixed DNA sequence to set
         """
         self._sequence = fixed_sequence
-        self._starred_sequence = dv.wc(self.sequence)
+        self._starred_sequence = dv.wc(fixed_sequence)
         self._set_subdomain_sequences(fixed_sequence)
         self._set_parent_sequence(fixed_sequence)
         self.fixed = True
@@ -1517,48 +1532,57 @@ def domains_not_substrings_of_each_other_domain_pair_constraint(
              (in particular, if they are equal length, then they are not the same domain)
     """
 
-    def domains_not_substrings_of_each_other(s1: str, s2: str,
-                                             domain1: Optional[Domain], domain2: Optional[Domain]) -> float:
-        if len(s1) > len(s2):
-            s1, s2 = s2, s1
-        if domain1 is not None:
-            assert s1 == domain1.sequence
-        if domain2 is not None:
-            assert s2 == domain2.sequence
-        if check_complements:
-            c1 = dv.wc(s1)
-            # by symmetry, only need to check c1 versus s2 for WC complement, since
-            # (s1 not in s2 <==> c1 in c2) and (c1 in s2 <==> s1 in c2)
-            return 1.0 if s1 in s2 or c1 in s2 else 0.0
-        else:
-            return 1.0 if s1 in s2 else 0.0
-
-    def summary(domain1: Domain, domain2: Domain) -> str:
-        s1 = domain1.sequence
-        s2 = domain2.sequence
+    # def evaluate(s1: str, s2: str, domain1: Optional[Domain], domain2: Optional[Domain]) -> float:
+    def evaluate(seqs: Tuple[str, ...],
+                 domains: Optional[Tuple[Domain, Domain]]) -> Tuple[float, str]:  # noqa
+        s1, s2 = seqs
         if len(s1) > len(s2):
             s1, s2 = s2, s1
+        summary = ''
+        score = 0.0
         passed = True
-        result = 'nothing to report'
         if s1 in s2:
-            result = f'{s1} is a substring of {s2}'
+            score = 1.0
+            summary = f'{s1} is a substring of {s2}'
             passed = False
         if check_complements:
+            # by symmetry, only need to check c1 versus s2 for WC complement, since
+            # (s1 not in s2 <==> c1 in c2) and (c1 in s2 <==> s1 in c2)
             c1 = dv.wc(s1)
             if c1 in s2:
                 msg = f'{c1} is a substring of {s2}'
                 if not passed:
-                    result += f'; {msg}'
+                    summary += f'; {msg}'
                 else:
-                    result = msg
-                # passed = False
-        return result
+                    summary = msg
+                score += 1.0
+        return score, summary
+
+    # def summary(domain1: Domain, domain2: Domain) -> str:
+    #     s1 = domain1.sequence()
+    #     s2 = domain2.sequence()
+    #     if len(s1) > len(s2):
+    #         s1, s2 = s2, s1
+    #     passed = True
+    #     result = 'nothing to report'
+    #     if s1 in s2:
+    #         result = f'{s1} is a substring of {s2}'
+    #         passed = False
+    #     if check_complements:
+    #         c1 = dv.wc(s1)
+    #         if c1 in s2:
+    #             msg = f'{c1} is a substring of {s2}'
+    #             if not passed:
+    #                 result += f'; {msg}'
+    #             else:
+    #                 result = msg
+    #             # passed = False
+    #     return result
 
     return DomainPairConstraint(description='domains not substrings of each other',
                                 short_description=short_description,
                                 weight=weight,
-                                evaluate=domains_not_substrings_of_each_other,
-                                summary=summary)
+                                _evaluate=evaluate)
 
 
 @dataclass
@@ -1618,7 +1642,7 @@ def to_json_serializable(self, suppress_indent: bool = True,
         return NoIndent(dct)
 
     @staticmethod
-    def from_json(json_map: Dict[str, Any]) -> IDTFields:
+    def from_json_serializable(json_map: Dict[str, Any]) -> IDTFields:
         scale = mandatory_field(IDTFields, json_map, idt_scale_key)
         purification = mandatory_field(IDTFields, json_map, idt_purification_key)
         plate = json_map.get(idt_plate_key)
@@ -1647,7 +1671,47 @@ def _check_idt_string_not_none_or_empty(value: str, field_name: str) -> None:
 
 
 @dataclass
-class Strand(JSONSerializable, Generic[StrandLabel, DomainLabel]):
+class StrandPair(Part, Generic[StrandLabel, DomainLabel]):
+    strand1: Strand
+    strand2: Strand
+
+    @property
+    def name(self) -> str:
+        return f'{self.strand1.name}, {self.strand2.name}'
+
+    def name_of_part_type(self) -> str:
+        return 'strand pair'
+
+    def individual_parts(self) -> Tuple[Strand, Strand]:
+        return self.strand1, self.strand2
+
+
+@dataclass
+class Complex(Part, Generic[StrandLabel, DomainLabel]):
+    strands: Tuple[Strand, ...]
+
+    @property
+    def name(self) -> str:
+        return ', '.join(strand.name for strand in self.strands)
+
+    def name_of_part_type(self) -> str:
+        return 'complex'
+
+    def individual_parts(self) -> Tuple[Strand, ...]:
+        return self.strands
+
+    def __iter__(self) -> Iterator[Strand]:
+        return iter(self.strands)
+
+    def __len__(self) -> int:
+        return len(self.strands)
+
+    def __getitem__(self, i: int) -> Strand:
+        return self.strands[i]
+
+
+@dataclass
+class Strand(JSONSerializable, Generic[StrandLabel, DomainLabel], Part):
     """Represents a DNA strand, made of several :any:`Domain`'s. """
 
     domains: List[Domain[DomainLabel]]
@@ -1773,6 +1837,9 @@ def __init__(self,
 
         self.compute_derived_fields()
 
+    def name_of_part_type(self) -> str:
+        return 'strand'
+
     def clone(self, name: Optional[str]) -> Strand:
         """
         Returns a copy of this :any:`Strand`. The copy is "shallow" in that the :any:`Domain`'s are shared.
@@ -1945,7 +2012,7 @@ def assign_dna(self, sequence: str) -> None:
         for domain in self.domains:
             end = start + domain.length
             domain_sequence = sequence[start:end]
-            domain.sequence = domain_sequence
+            domain.set_sequence(domain_sequence)
             start = end
 
     @property
@@ -2592,7 +2659,7 @@ def from_json(json_str: str,
                   pool_with_name: Optional[Dict[str, DomainPool]] = None,
                   strand_label_decoder: Callable[[Any], StrandLabel] = lambda label: label,
                   domain_label_decoder: Callable[[Any], DomainLabel] = lambda label: label,
-                  ) -> 'Design[StrandLabel, DomainLabel]':
+                  ) -> Design[StrandLabel, DomainLabel]:
         """
         :param json_str:
             The string representing the :any:`Design` as a JSON object.
@@ -2720,345 +2787,63 @@ def summary_of_constraints(self, report_only_violations: bool, violation_set: Vi
             summaries.append(summary)
 
         score = violation_set.total_score()
-        score_unfixed = violation_set.total_nonfixed_score()
-        score_total_summary =   f'total score of constraint violations: {score:.1f}'
+        score_unfixed = violation_set.total_score_nonfixed()
+        score_total_summary = f'total score of constraint violations: {score:.1f}'
         score_unfixed_summary = f'total score of unfixed constraint violations: {score_unfixed:.1f}'
 
-        summary = score_total_summary + '\n'
-        if score_unfixed != score:
-            summary += score_unfixed_summary + '\n\n'
-        else:
-            summary += '\n'
+        summary = (score_total_summary + '\n'
+                   + (score_unfixed_summary + '\n\n' if score_unfixed != score else '\n')
+                   + '\n'.join(summaries))
 
-        return summary + '\n'.join(summaries)
+        return summary
 
     def summary_of_constraint(self, constraint: Constraint, report_only_violations: bool,
                               violation_set: ViolationSet) -> str:
-        # summary of constraint only if not a DomainConstraint in a DomainPool
-        # or a StrandConstraint in a StrandGroup
-        report: ConstraintReport
-        content: str
-        num_violations: int
-        num_checks: int
-        if isinstance(constraint, DomainConstraint):
-            report = self.summary_of_domain_constraint(constraint, report_only_violations)
-        elif isinstance(constraint, StrandConstraint):
-            report = self.summary_of_strand_constraint(constraint, report_only_violations)
-        elif isinstance(constraint, DomainPairConstraint):
-            report = self.summary_of_domain_pair_constraint(constraint, report_only_violations)
-        elif isinstance(constraint, StrandPairConstraint):
-            report = self.summary_of_strand_pair_constraint(constraint, report_only_violations)
-        elif isinstance(constraint, DomainsConstraint):
-            report = self.summary_of_domains_constraint(constraint, report_only_violations)
-        elif isinstance(constraint, StrandsConstraint):
-            report = self.summary_of_strands_constraint(constraint, report_only_violations)
-        elif isinstance(constraint, DomainPairsConstraint):
-            report = self.summary_of_domain_pairs_constraint(constraint, report_only_violations)
-        elif isinstance(constraint, StrandPairsConstraint):
-            report = self.summary_of_strand_pairs_constraint(constraint, report_only_violations)
-        elif isinstance(constraint, DesignConstraint):
-            report = self.summary_of_design_constraint(constraint, report_only_violations)
-        elif isinstance(constraint, ComplexConstraint):
-            report = self.summary_of_complex_constraint(constraint, report_only_violations)
-        else:
-            content = f'skipping summary of constraint {constraint.description}; ' \
-                      f'unrecognized type {type(constraint)}'
-            report = ConstraintReport(constraint=constraint, content=content, num_violations=0, num_checks=0)
-
-        report.constraint = constraint
-
-        if _no_summary_string in report.content:
-            report = ConstraintReport(constraint=constraint, content=_no_summary_string,
-                                      num_violations=0, num_checks=0)
-
-        summary = add_header_to_content_of_summary(report, violation_set)
-        return summary
-
-    def summary_of_domain_constraint(self, constraint: DomainConstraint, report_only_violations: bool) \
-            -> ConstraintReport:
-        num_violations = 0
-        num_checks = 0
-        domains_to_check = self.domains if constraint.domains is None else constraint.domains
-        max_domain_name_length = max(len(domain.name) for domain in domains_to_check)
-
-        summaries = []
-
-        fixed_domains = [domain for domain in domains_to_check if domain.fixed]
-        unfixed_domains = [domain for domain in domains_to_check if not domain.fixed]
-        for domains_to_check, header_name in [(unfixed_domains, 'unfixed domains'),
-                                              (fixed_domains, 'fixed domains')]:
-            if len(domains_to_check) > 0:
-                report = self._summary_of_domains_in_domain_constraint(
-                    constraint, report_only_violations, domains_to_check, max_domain_name_length)
-                summary = _small_header(header_name, "=") + f'\n{report.content}\n'
-                num_violations += report.num_violations
-                num_checks += report.num_checks
-                summaries.append(summary)
+        if isinstance(constraint, (DomainConstraint, StrandConstraint,
+                                   DomainPairConstraint, StrandPairConstraint, ComplexConstraint,
+                                   DomainsConstraint, StrandsConstraint,
+                                   DomainPairsConstraint, StrandPairsConstraint, ComplexesConstraint)):
+            summaries = []
+            num_violations = 0
+            num_checks = violation_set.num_checked[constraint]
+            part_type_name = constraint.part_name()
+
+            violations_nonfixed = violation_set.violations_nonfixed[constraint]
+            violations_fixed = violation_set.violations_fixed[constraint]
+            for violations, header_name in [(violations_nonfixed, f"unfixed {part_type_name}"),
+                                            (violations_fixed, f"fixed {part_type_name}")]:
+                if len(violations) == 0:
+                    continue
 
-        content = ''.join(summaries)
-        report = ConstraintReport(constraint=constraint, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
+                max_part_name_length = max(len(violation.part.name) for violation in violations)
+                num_violations += len(violations)
 
-    # this function reuses code between summarizing fixed and unfixed domains
-    @staticmethod
-    def _summary_of_domains_in_domain_constraint(constraint: DomainConstraint,
-                                                 report_only_violations: bool,
-                                                 domains: Iterable[Domain[DomainLabel]],
-                                                 max_domain_name_length: int) -> ConstraintReport:
-        num_violations = 0
-        num_checks = 0
-        lines_and_excesses: List[Tuple[str, float]] = []
-        for fixed_domain in domains:
-            num_checks += 1
-            seq = fixed_domain.sequence
-            excess = constraint(seq, fixed_domain)
-            passed = excess <= 0.0
-            if not passed:
-                num_violations += 1
-            if not report_only_violations or (report_only_violations and not passed):
-                summary = constraint.generate_summary(fixed_domain, False)
-                line = f'domain {fixed_domain.name:{max_domain_name_length}}: ' \
-                       f'{summary} ' \
-                       f'{"" if passed else " **violation**"}'
-                lines_and_excesses.append((line, excess))
-
-        # put in descending order of excess
-        lines_and_excesses.sort(key=lambda line_: ' **violation**' not in line_)  # put violations first
-
-        lines = (line for line, _ in lines_and_excesses)
-        content = '\n'.join(lines)
-        report = ConstraintReport(constraint=constraint, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
-
-    def summary_of_strand_constraint(self, constraint: StrandConstraint, report_only_violations: bool) \
-            -> ConstraintReport:
-        num_violations = 0
-        num_checks = 0
-        all_strands_to_check = self.strands if constraint.strands is None else constraint.strands
-        max_strand_name_length = max(len(strand.name) for strand in all_strands_to_check)
-
-        summaries = []
-
-        fixed_strands = [strand for strand in all_strands_to_check if strand.fixed]
-        unfixed_strands = [strand for strand in all_strands_to_check if not strand.fixed]
-        for strands_to_check, header_name in [(unfixed_strands, 'unfixed strands'),
-                                              (fixed_strands, 'fixed strands')]:
-            if len(strands_to_check) > 0:
-                report = self._summary_of_strands_in_strand_constraint(
-                    constraint, report_only_violations, strands_to_check, max_strand_name_length)
-                summary = _small_header(header_name, "=") + f'\n{report.content}\n'
-                num_violations += report.num_violations
-                num_checks += report.num_checks
-                summaries.append(summary)
+                lines_and_scores: List[Tuple[str, float]] = []
+                for violation in violations:
+                    line = f'{part_type_name} {violation.part.name:{max_part_name_length}}: ' \
+                           f'{violation.summary} '
+                    lines_and_scores.append((line, violation.score))
 
-        content = ''.join(summaries)
-        report = ConstraintReport(constraint=constraint, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
+                lines_and_scores.sort(key=lambda line_and_score: line_and_score[1], reverse=True)
 
-    # this function reuses code between summarizing fixed and unfixed strands
-    @staticmethod
-    def _summary_of_strands_in_strand_constraint(constraint: StrandConstraint,
-                                                 report_only_violations: bool,
-                                                 strands: Iterable[Strand[StrandLabel, DomainLabel]],
-                                                 max_strand_name_length: int) -> ConstraintReport:
-        num_violations = 0
-        num_checks = 0
-        lines_and_excesses: List[Tuple[str, float]] = []
-        for strand in strands:
-            num_checks += 1
-            excess = constraint(strand.sequence(), strand)
-            passed = excess <= 0.0
-            if not passed:
-                num_violations += 1
-            if not report_only_violations or (report_only_violations and not passed):
-                summary = constraint.generate_summary(strand, False)
-                line = f'strand {strand.name:{max_strand_name_length}}: ' \
-                       f'{summary} ' \
-                       f'{"" if passed else " **violation**"}'
-                lines_and_excesses.append((line, excess))
-
-        # put in descending order of excess
-        lines_and_excesses.sort(key=lambda line_and_excess: line_and_excess[1], reverse=True)
-
-        lines = (line for line, _ in lines_and_excesses)
-        content = '\n'.join(lines)
-        report = ConstraintReport(constraint=constraint, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
-
-    def summary_of_domain_pair_constraint(self, constraint: DomainPairConstraint,
-                                          report_only_violations: bool) -> ConstraintReport:
-        pairs_to_check = constraint.pairs if constraint.pairs is not None else all_pairs(self.domains)
-
-        max_domain_name_length = max(len(domain.name) for domain in _flatten(pairs_to_check))
-
-        num_violations = 0
-        num_checks = 0
-        lines_and_excesses: List[Tuple[str, float]] = []
-        for domain1, domain2 in pairs_to_check:
-            num_checks += 1
-            seq1, seq2 = domain1.sequence, domain2.sequence
-            excess = constraint(seq1, seq2, domain1, domain2)
-            passed = excess <= 0.0
-            if not passed:
-                num_violations += 1
-            if not report_only_violations or (report_only_violations and not passed):
-                summary = constraint.generate_summary((domain1, domain2), False)
-                line = (f'domains '
-                        f'{domain1.name:{max_domain_name_length}}, '
-                        f'{domain2.name:{max_domain_name_length}}: '
-                        f'{summary}'
-                        f'{"" if passed else "  **violation**"}')
-                lines_and_excesses.append((line, excess))
-
-        # put in descending order of excess
-        lines_and_excesses.sort(key=lambda line_and_excess: line_and_excess[1], reverse=True)
-
-        lines = (line for line, _ in lines_and_excesses)
-        content = '\n'.join(lines)
-        report = ConstraintReport(constraint=constraint, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
-
-    def summary_of_strand_pair_constraint(self, constraint: StrandPairConstraint,
-                                          report_only_violations: bool) -> ConstraintReport:
-        num_violations = 0
-        num_checks = 0
-        all_pairs_to_check: List[
-            Tuple[Strand, Strand]] = constraint.pairs if constraint.pairs is not None else all_pairs(
-            self.strands)
-
-        # distinguish between pairs in which both strands are fixed (so cannot remove violation)
-        # versus those pairs in which at least one element of the pair is unfixed
-        both_fixed_pairs = [(s1, s2) for s1, s2 in all_pairs_to_check if s1.fixed and s2.fixed]
-        one_unfixed_pairs = [(s1, s2) for s1, s2 in all_pairs_to_check if not (s1.fixed and s2.fixed)]
-
-        max_strand_name_length = max(len(strand.name) for strand in _flatten(all_pairs_to_check))
-
-        summaries = []
-
-        for pairs_to_check, header_name in [(one_unfixed_pairs, 'pairs with at least one unfixed'),
-                                            (both_fixed_pairs, 'pairs with both fixed')]:
-            if len(pairs_to_check) > 0:
-                report = self._summary_of_strand_pairs_in_strand_pair_constraint(
-                    constraint, report_only_violations, pairs_to_check, max_strand_name_length)
-                summary = _small_header(header_name, "=") + f'\n{report.content}\n'
-                num_violations += report.num_violations
-                num_checks += report.num_checks
+                lines = (line for line, _ in lines_and_scores)
+                content = '\n'.join(lines)
+                summary = _small_header(header_name, "=") + f'\n{content}\n'
                 summaries.append(summary)
 
-        content = ''.join(summaries)
-        report = ConstraintReport(constraint=constraint, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
+            content = ''.join(summaries)
+            report = ConstraintReport(constraint=constraint, content=content,
+                                      num_violations=num_violations, num_checks=num_checks)
 
-    # this function reuses code between summarizing fixed and unfixed strands
-    @staticmethod
-    def _summary_of_strand_pairs_in_strand_pair_constraint(constraint: StrandPairConstraint,
-                                                           report_only_violations: bool,
-                                                           pairs: Iterable[
-                                                               Tuple[Strand[StrandLabel, DomainLabel],
-                                                                     Strand[StrandLabel, DomainLabel]]],
-                                                           max_strand_name_length: int) -> ConstraintReport:
-        num_violations = 0
-        num_checks = 0
-        lines_and_excesses: List[Tuple[str, float]] = []
-        for strand1, strand2 in pairs:
-            num_checks += 1
-            excess = constraint(strand1.sequence(), strand2.sequence(), strand1, strand2)
-            passed = excess <= 0.0
-            if not passed:
-                num_violations += 1
-            if not report_only_violations or (report_only_violations and not passed):
-                summary = constraint.generate_summary((strand1, strand2), False)
-                line = (f'strands '
-                        f'{strand1.name:{max_strand_name_length}}, '
-                        f'{strand2.name:{max_strand_name_length}}: '
-                        f'{summary}'
-                        f'{"" if passed else "  **violation**"}')
-                lines_and_excesses.append((line, excess))
-
-        # put in descending order of excess
-        lines_and_excesses.sort(key=lambda line_and_excess: line_and_excess[1], reverse=True)
-
-        lines = (line for line, _ in lines_and_excesses)
-        content = '\n'.join(lines)
-        report = ConstraintReport(constraint=constraint, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
+        elif isinstance(constraint, DesignConstraint):
+            raise NotImplementedError()
+        else:
+            content = f'skipping summary of constraint {constraint.description}; ' \
+                      f'unrecognized type {type(constraint)}'
+            report = ConstraintReport(constraint=constraint, content=content, num_violations=0, num_checks=0)
 
-    @staticmethod
-    def summary_of_complex_constraint(constraint: 'ComplexConstraint',
-                                      report_only_violations: bool) -> ConstraintReport:
-        num_violations = 0
-        num_checks = 0
-
-        lines_and_excesses: List[Tuple[str, float]] = []
-        for strand_complex in constraint.complexes:
-            num_checks += 1
-            excess = constraint(strand_complex)
-            passed = excess <= 0.0
-            if not passed:
-                num_violations += 1
-            if not report_only_violations or (report_only_violations and not passed):
-                summary = constraint.generate_summary(strand_complex, False)
-                strand_names = ', '.join([f'{strand.name}' for strand in strand_complex])
-                line = (f'strand complex: '
-                        f'{strand_names}'
-                        f'{"" if passed else "  **violation**"}'
-                        f'\n{summary}')
-                lines_and_excesses.append((line, excess))
-
-        # put in descending order of excess
-        lines_and_excesses.sort(key=lambda line_and_excess: line_and_excess[1], reverse=True)
-
-        lines = (line for line, _ in lines_and_excesses)
-        content = '\n'.join(lines)
-        report = ConstraintReport(constraint=constraint, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
-
-    def summary_of_domains_constraint(self, constraint: DomainsConstraint,
-                                      report_only_violations: bool) -> ConstraintReport:
-        # summary = f'domains\n{constraint.generate_summary(self.domains)}'
-        report = constraint.generate_summary(self.domains, report_only_violations)
-        return report
-
-    def summary_of_strands_constraint(self, constraint: StrandsConstraint,
-                                      report_only_violations: bool) -> ConstraintReport:
-        # summary = f'strands\n{constraint.generate_summary(self.strands)}'
-        report = constraint.generate_summary(self.strands, report_only_violations)
-        return report
-
-    def summary_of_domain_pairs_constraint(self, constraint: DomainPairsConstraint,
-                                           report_only_violations: bool) -> ConstraintReport:
-        pairs_to_check = constraint.pairs if constraint.pairs is not None else all_pairs(self.domains)
-        # summary = f'domain pairs\n{constraint.generate_summary(pairs_to_check)}'
-        report = constraint.generate_summary(pairs_to_check, report_only_violations) \
-            if len(pairs_to_check) > 0 \
-            else ConstraintReport(constraint=constraint,
-                                  content='constraint.pairs is empty; nothing to report',
-                                  num_violations=0, num_checks=0)
-        return report
-
-    def summary_of_strand_pairs_constraint(self, constraint: StrandPairsConstraint,
-                                           report_only_violations: bool) -> ConstraintReport:
-        pairs_to_check = constraint.pairs if constraint.pairs is not None else all_pairs(self.strands)
-        report = constraint.generate_summary(pairs_to_check, report_only_violations) \
-            if len(pairs_to_check) > 0 \
-            else ConstraintReport(constraint=constraint,
-                                  content='constraint.pairs is empty; nothing to report',
-                                  num_violations=0, num_checks=0)
-        return report
-
-    def summary_of_design_constraint(self, constraint: DesignConstraint,
-                                     report_only_violations: bool) -> ConstraintReport:
-        # summary = f'design\n{constraint.generate_summary(self)}'
-        report = constraint.generate_summary(self, report_only_violations)
-        return report
+        summary = add_header_to_content_of_summary(report, violation_set)
+        return summary
 
     @staticmethod
     def from_scadnano_file(sc_filename: str,
@@ -3199,7 +2984,7 @@ def from_scadnano_design(sc_design: sc.Design[StrandLabel, DomainLabel],
                             if fix_assigned_sequences:
                                 dsd_domain.set_fixed_sequence(domain_sequence)
                             else:
-                                dsd_domain.sequence = domain_sequence
+                                dsd_domain.set_sequence(domain_sequence)
 
                 # set domain labels
                 for dsd_domain, sc_domain in zip(dsd_strand.domains, sc_strand.domains):
@@ -3370,18 +3155,19 @@ def copy_sequences_from(self, other: Design) -> None:
         for domain in self.domains:
             other_domain = other.domains_by_name[domain.name]
             if other_domain.fixed:
-                domain.set_fixed_sequence(other_domain.sequence)
+                domain.set_fixed_sequence(other_domain.sequence())
             elif other_domain.has_sequence():
-                domain.sequence = other_domain.sequence
+                domain.set_sequence(other_domain.sequence())
 
         # no need to compute_derived_fields if we already called it above,
         # since new sequences won't change derived fields
         if not computed_derived_fields:
             self.compute_derived_fields()
 
+
 def add_header_to_content_of_summary(report: ConstraintReport, violation_set: ViolationSet) -> str:
     score = violation_set.score_of_constraint(report.constraint)
-    score_unfixed = violation_set.nonfixed_score_of_constraint(report.constraint)
+    score_unfixed = violation_set.score_of_constraint_nonfixed(report.constraint)
 
     if score != score_unfixed:
         summary_score_unfixed = f'\n* unfixed score of violations: {score_unfixed:.2f}'
@@ -3405,14 +3191,14 @@ def add_header_to_content_of_summary(report: ConstraintReport, violation_set: Vi
 DesignPart = TypeVar('DesignPart',
                      Domain,
                      Strand,
-                     Tuple[Domain, Domain],  # noqa
-                     Tuple[Strand, Strand],  # noqa
-                     Complex,  # noqa
-                     Iterable[Domain],  # noqa
-                     Iterable[Strand],  # noqa
-                     Iterable[Tuple[Domain, Domain]],  # noqa
-                     Iterable[Tuple[Strand, Strand]],  # noqa
-                     Iterable[Complex],  # noqa
+                     DomainPair,
+                     StrandPair,
+                     Complex,
+                     # Iterable[Domain],
+                     # Iterable[Strand],
+                     # Iterable[Tuple[Domain, Domain]],
+                     # Iterable[Tuple[Strand, Strand]],
+                     # Iterable[Complex],
                      Design)
 
 
@@ -3425,12 +3211,18 @@ class Violation(Generic[DesignPart]):
     constraint: Constraint
     # :any:`Constraint` that was violated to result in this :any:`Violation`.
 
+    part: DesignPart
+    # DesignPart that caused this violation
+
     domains: FrozenSet[Domain]  # = field(init=False, hash=False, compare=False, default=None)
     # :any:`Domain`'s that were involved in violating :py:data:`Violation.constraint`
 
-    _unweighted_score: float
+    summary: str
+
+    score: float
 
-    def __init__(self, constraint: Constraint, domains: Iterable[Domain], score: float):
+    def __init__(self, constraint: Constraint, part: DesignPart, domains: Iterable[Domain],
+                 score: float, summary: str) -> None:
         # :param constraint:
         #     :any:`Constraint` that was violated to result in this
         # :param domains:
@@ -3439,16 +3231,15 @@ def __init__(self, constraint: Constraint, domains: Iterable[Domain], score: flo
         #     total "score" of this violation, typically something like an excess energy over a
         #     threshold, squared, multiplied by the :data:`Constraint.weight`
         object.__setattr__(self, 'constraint', constraint)
+        object.__setattr__(self, 'part', part)
         domains_frozen = frozenset(domains)
         object.__setattr__(self, 'domains', domains_frozen)
-        object.__setattr__(self, '_unweighted_score', score)
-
-    @property
-    def score(self) -> float:
-        return self.constraint.weight * self._unweighted_score
+        object.__setattr__(self, 'score', score)
+        object.__setattr__(self, 'summary', summary)
 
     def __repr__(self) -> str:
-        return f'Violation({self.constraint.short_description}, score={self._unweighted_score:.2f})'
+        return f'Violation({self.constraint.short_description}, score={self.score:.2f}, ' \
+               f'summary={self.summary})'
 
     def __str__(self) -> str:
         return repr(self)
@@ -3469,19 +3260,33 @@ class ViolationSet:
     # It is designed to be efficiently updateable when a single :any:`Domain` changes, to efficiently update
     # only those violations of :any:`Constraint`'s that could have been affected by the changed :any:`Domain`.
 
-    all_violations: OrderedSet[Violation] = field(default_factory=OrderedSet)
-    # Set of all :any:`Violation`'s.
+    violations_all: Dict[Constraint, OrderedSet[Violation]]
+    # Dict mapping each :any:`Constraint` to the set of all :any:`Violation`'s of it.
 
-    domain_to_violations: Dict[Domain, OrderedSet[Violation]] = field(
-        default_factory=lambda: defaultdict(OrderedSet))
+    domain_to_violations: Dict[Domain, OrderedSet[Violation]]
     # Dict mapping each :any:`constraint.Domain` to the set of all :any:`Violation`'s for which it is blamed
 
-    non_fixed_violations: OrderedSet[Violation] = field(default_factory=OrderedSet)
+    violations_nonfixed: Dict[Constraint, OrderedSet[Violation]]
+    # Dict mapping each :any:`Constraint` to the set of all :any:`Violations`
+    # that are associated to non-fixed :any:`constraint.Domain`'s.
+
+    violations_fixed: Dict[Constraint, OrderedSet[Violation]]
+    # Dict mapping each :any:`Constraint` to the set of all :any:`Violations`
+    # that are associated to fixed :any:`constraint.Domain`'s.
+
+    num_checked: Dict[Constraint, int]
+
+    # number of instances of each :any:`Constraint` that were checked
 
-    # Set of all :any:`Violations` that are associated to non-fixed :any:`constraint.Domain`'s.
+    def __init__(self) -> None:
+        self.violations_all = defaultdict(OrderedSet)
+        self.domain_to_violations = defaultdict(OrderedSet)
+        self.violations_nonfixed = defaultdict(OrderedSet)
+        self.violations_fixed = defaultdict(OrderedSet)
+        self.num_checked = defaultdict(int)
 
     def __repr__(self):
-        lines = "\n  ".join(map(str, self.all_violations))
+        lines = "\n  ".join(map(str, self.violations_all.values()))
         return f'ViolationSet(\n  {lines})'
 
     def __str__(self):
@@ -3493,10 +3298,14 @@ def update(self, new_violations: Dict[Domain, OrderedSet[Violation]]) -> None:
         # :param new_violations: dict mapping each :any:`Domain` to the set of :any:`Violation`'s
         #                        for which it is blamed
         for domain, domain_violations in new_violations.items():
-            self.all_violations.update(domain_violations)
+            for violation in domain_violations:
+                self.violations_all[violation.constraint].add(violation)
+                if not domain.fixed:
+                    self.violations_nonfixed[violation.constraint].add(violation)
+                else:
+                    self.violations_fixed[violation.constraint].add(violation)
             self.domain_to_violations[domain].update(domain_violations)
-            if not domain.fixed:
-                self.non_fixed_violations.update(domain_violations)
+
 
     def clone(self) -> ViolationSet:
         # Returns a deep-ish copy of this :any:`ViolationSet`.
@@ -3512,8 +3321,26 @@ def clone(self) -> ViolationSet:
         domain_to_violations_deep_copy = defaultdict(OrderedSet, self.domain_to_violations)
         for domain, violations in domain_to_violations_deep_copy.items():
             domain_to_violations_deep_copy[domain] = OrderedSet(violations)
-        return ViolationSet(OrderedSet(self.all_violations), domain_to_violations_deep_copy,
-                            OrderedSet(self.non_fixed_violations))
+
+        violations_all_deep_copy = defaultdict(OrderedSet, self.violations_all)
+        for constraint, violations in violations_all_deep_copy.items():
+            violations_all_deep_copy[constraint] = OrderedSet(violations)
+
+        violations_nonfixed_deep_copy = defaultdict(OrderedSet, self.violations_nonfixed)
+        for constraint, violations in violations_nonfixed_deep_copy.items():
+            violations_nonfixed_deep_copy[constraint] = OrderedSet(violations)
+
+        violations_fixed_deep_copy = defaultdict(OrderedSet, self.violations_fixed)
+        for constraint, violations in violations_fixed_deep_copy.items():
+            violations_fixed_deep_copy[constraint] = OrderedSet(violations)
+
+        result = ViolationSet()
+        result.violations_all = violations_all_deep_copy
+        result.domain_to_violations = domain_to_violations_deep_copy
+        result.violations_nonfixed = violations_nonfixed_deep_copy
+        result.violations_fixed = violations_fixed_deep_copy
+
+        return result
 
     def remove_violations_of_domain(self, domain: Domain) -> None:
         # Removes any :any:`Violation`'s blamed on `domain`.
@@ -3522,23 +3349,42 @@ def remove_violations_of_domain(self, domain: Domain) -> None:
         # XXX: need to make a copy of this set, since we are modifying the sets in place
         # (values in self.domain_to_violations)
         violations_of_domain = set(self.domain_to_violations[domain])
-        self.all_violations -= violations_of_domain
-        self.non_fixed_violations -= violations_of_domain
+
+        for violations in self.violations_all.values():
+            violations -= violations_of_domain
+        for violations in self.violations_nonfixed.values():
+            violations -= violations_of_domain
+        for violations in self.violations_fixed.values():
+            violations -= violations_of_domain
+
         for violations_of_other_domain in self.domain_to_violations.values():
             violations_of_other_domain -= violations_of_domain
+
         assert len(self.domain_to_violations[domain]) == 0
 
     def total_score(self) -> float:
         """
         :return: Total score of all violations.
         """
-        return sum(violation.score for violation in self.all_violations)
+        return sum(violation.score
+                   for violations in self.violations_all.values()
+                   for violation in violations)
+
+    def total_score_nonfixed(self) -> float:
+        # :return:
+        #     Total score of all violations attributed to :any:`constraint.Domain`'s with
+        #     :any:`constraint.Domain.fixed` = False.
+        return sum(violation.score
+                   for violations in self.violations_nonfixed.values()
+                   for violation in violations)
 
-    def total_nonfixed_score(self) -> float:
+    def total_score_fixed(self) -> float:
         # :return:
         #     Total score of all violations attributed to :any:`constraint.Domain`'s with
         #     :any:`constraint.Domain.fixed` = False.
-        return sum(violation.score for violation in self.non_fixed_violations)
+        return sum(violation.score
+                   for violations in self.violations_fixed.values()
+                   for violation in violations)
 
     def score_of_constraint(self, constraint: Constraint) -> float:
         """
@@ -3547,31 +3393,52 @@ def score_of_constraint(self, constraint: Constraint) -> float:
         :return:
             Total score of all violations due to `constraint`.
         """
-        return sum(violation.score for violation in self.all_violations if violation.constraint == constraint)
+        return sum(violation.score
+                   for violations in self.violations_all.values()
+                   for violation in violations
+                   if violation.constraint == constraint)
 
-    def nonfixed_score_of_constraint(self, constraint: Constraint) -> float:
+    def score_of_constraint_nonfixed(self, constraint: Constraint) -> float:
         # :param constraint:
         #     constraint to filter scores on
         # :return:
         #     Total score of all nonfixed violations due to `constraint`.
-        return sum(
-            violation.score for violation in self.non_fixed_violations if violation.constraint == constraint)
+        return sum(violation.score
+                   for violations in self.violations_nonfixed.values()
+                   for violation in violations
+                   if violation.constraint == constraint)
+
+    def score_of_constraint_fixed(self, constraint: Constraint) -> float:
+        # :param constraint:
+        #     constraint to filter scores on
+        # :return:
+        #     Total score of all fixed violations due to `constraint`.
+        return sum(violation.score
+                   for violations in self.violations_fixed.values()
+                   for violation in violations
+                   if violation.constraint == constraint)
 
     def num_violations(self) -> float:
         # :return: Total number of violations.
-        return len(self.all_violations)
+        return sum(len(violations) for violations in self.violations_all.values())
 
-    def num_nonfixed_violations(self) -> float:
+    def num_violations_nonfixed(self) -> float:
         # :return: Total number of nonfixed violations.
-        return len(self.non_fixed_violations)
+        return sum(len(violations) for violations in self.violations_nonfixed.values())
 
+    def num_violations_fixed(self) -> float:
+        # :return: Total number of fixed violations.
+        return sum(len(violations) for violations in self.violations_fixed.values())
 
-@dataclass(frozen=True, eq=False)  # type: ignore
+
+@dataclass(frozen=True, eq=False)
 class Constraint(ABC, Generic[DesignPart]):
     description: str
-    """Description of the constraint, e.g., 'strand has secondary structure exceeding -2.0 kcal/mol'."""
+    """
+    Description of the constraint, e.g., 'strand has secondary structure exceeding -2.0 kcal/mol'.
+    """
 
-    short_description: str = field(default='')
+    short_description: str = ''
     """
     Very short description of the constraint suitable for compactly logging to the screen, e.g., 'strand_ss'
     """
@@ -3596,7 +3463,7 @@ class Constraint(ABC, Generic[DesignPart]):
     threshold than to reduce (by the same amount) a violation only 1 kcal/mol in excess of its threshold.
     """
 
-    threaded: bool = False
+    parallel: bool = False
     """
     Whether or not to use parallelization across multiple processes to take advantage of multiple
     processors/cores. (Not applicable to some types of constraints. 
@@ -3604,12 +3471,32 @@ class Constraint(ABC, Generic[DesignPart]):
 
     sequence_only: bool = True
     """
-    If :py:data:`Constraint.threaded` is True, then this should be set to True so that only the 
+    If :py:data:`Constraint.parallel` is True, then this should be set to True so that only the 
     sequence is serialized when passing data to other processes for parallel computation.
     Otherwise, significant time is spent serializing objects such as :any:`Strand` or :any:`Domain`,
     which is slower than not using parallelization in the first place.
     """
 
+    _evaluate: Optional[Callable[[Tuple[str, ...], Optional[DesignPart]],
+                                 Tuple[float, str]]] = None
+    """
+    Function that evaluates the :any:`Constraint`. It takes as input a tuple of DNA sequences 
+    (Python strings) and an optional :any:`Part`, where :any:`Part` is one of 
+    :any:`Domain`, :any:`Strand`, :any:`DomainPair`, :any:`StrandPair`, or :any:`Complex`
+    (the latter being an alias for arbitrary-length tuple of :any:`Strand`'s).
+    
+    The second argument will be None if :data:`Constraint.parallel` is True (since it's more expensive
+    to serialize the :any:`Domain` and :any:`Strand` objects than strings for passing data to processes
+    executing in parallel). Thus, if the :any:`Constraint` needs to use more data about the :any:`Part`
+    than just its DNA sequence, by accessing the second argument, :data:`Constraint.parallel` should be set
+    to False.    
+    
+    Exactly one of _evaluate and _evaluate_bulk should be specified.
+    """
+
+    _evaluate_bulk: Optional[Callable[[Iterable[DesignPart]],
+                                      List[Tuple[DesignPart, float, str]]]] = None
+
     def __post_init__(self) -> None:
         if len(self.short_description) == 0:
             # self.short_description = self.description
@@ -3618,247 +3505,185 @@ def __post_init__(self) -> None:
         if self.weight <= 0:
             raise ValueError(f'weight must be positive but it is {self.weight}')
 
-        if self.threaded and not self.sequence_only:
-            raise ValueError('cannot have both threaded=True and sequence_only=False;\n'
+        if self.parallel and not self.sequence_only:
+            raise ValueError('cannot have both parallel=True and sequence_only=False;\n'
                              'if you want to use threading, the constraint must be sequence_only '
                              'for efficiency, since the pathos library is used for parallel processing, '
                              'and it is too inefficient to serialize anything more than the DNA sequence')
 
-    @abstractmethod
-    def generate_summary(self, design_part: DesignPart, report_only_violations: bool) -> ConstraintReport:
-        """
-        Method that helps to give a summary of how well parts of the :any:`Design` are performing in
-        satisfying this constraint. For example, useful for generating a report after assigning sequences
-        to a :any:`Design`.
+        if not isinstance(self, DesignConstraint):
+            if self._evaluate is None and self._evaluate_bulk is None:
+                raise ValueError('Exactly one of _evaluate or _evaluate_bulk must be specified, '
+                                 'but both are None.')
 
-        :param design_part:
-            part of :any:`Design` that this :any:`Constraint` references
-        :param report_only_violations:
-            Whether to report only violations of constraints, or all evaluations of constraints, including
-            those that passed.
-        :return:
-            :any:`ConstraintReport` summarizing of how well `design_part` "performs" for :any:`Constraint`.
-            For example, a :any:`StrandConstraint` checking the partition function energy of a single
-            :any:`Strand` may return a string such as "strand secondary structure: -2.3 kcal/mol".
-            A :any:`StrandPairsConstraint` checking all pairs of :any:`Strand`'s may return a longer
-            string reporting on every pair.
+            if self._evaluate is not None and self._evaluate_bulk is not None:
+                raise ValueError('Exactly one of _evaluate or _evaluate_bulk must be specified, '
+                                 'but both have been specified.')
+
+    def evaluate(self, seqs: Tuple[str, ...], part: Optional[DesignPart]) -> Tuple[float, str]:
         """
-        raise NotImplementedError('subclasses of Constraint must implement generate_summary')
+        Evaluates this :any:`Constraint` using function `_evaluate` supplied in constructor.
 
-    def _check_sequence_only(self, part: Any):
+        :param seqs:
+            sequence(s) of relevant :any:`Part`, e.g., if `part` is a pair of :any:`Strand`'s,
+            then `seqs` is a pair of strings
+        :param part:
+            the :any:`Part` to be evaluated. Might be None if parallelization is being used,
+            since it is cheaper to serialize only the sequence(s) than the entire :any:`Part`
+            for passing to other processes to evaluate in parallel.
+        :return:
+            pair (`excess`, `summary`), where `excess` is a float indicating how much the constraint
+            was violated (0.0 if satisfied) and `summary` is a string summarizing the violation (or lack
+            thereof), suitable for printing into a line of a report. For example, if measuring
+            complex free energy -2.5 kcal/mol of a strand and comparing against a threshold -1.0 kcal/mol,
+            `excess` might be the difference 1.5 between the energy and the threshold,
+            and `summary` might be the string "-2.5 kcal/mol".
+        """
+        if self._evaluate is None:
+            raise ValueError('Cannot call evaluate on a Constraint unless _evaluate is specified')
         if not self.sequence_only and part is None:
             raise AssertionError('if sequence_only is False, '
                                  'then the design part cannot be None, but it is')
 
+        excess, summary = (self._evaluate)(seqs, part)  # noqa
+        if excess < 0.0:
+            excess = 0.0
+        score = self.weight * self.score_transfer_function(excess)
+        return score, summary
+
+    def evaluate_bulk(self, parts: Iterable[DesignPart]) -> List[Tuple[DesignPart, float, str]]:
+        if self._evaluate_bulk is None:
+            raise ValueError('Cannot call evaluate_bulk on a Constraint unless _evaluate_bulk is specified')
+
+        results: List[Tuple[DesignPart, float, str]] = (self._evaluate_bulk)(parts)  # noqa
+        results_with_scores_transferred_and_weighted: List[Tuple[DesignPart, float, str]] = []
+        for part, excess, summary in results:
+            if excess < 0.0:
+                excess = 0.0
+            score = self.weight * self.score_transfer_function(excess)
+            results_with_scores_transferred_and_weighted.append((part, score, summary))
+        return results_with_scores_transferred_and_weighted
+
+    @staticmethod
+    @abstractmethod
+    def part_name() -> str:
+        """
+        :return:
+            name of the :any:`Part` that this :any:`Constraint` tests
+            (e.g., "domain", "strand pair")
+        """
+        raise NotImplementedError()
 
-_no_summary_string = f"No summary for this constraint. " \
-                     f"To generate one, pass a function as the parameter named " \
-                     f'"summary" when creating the Constraint.'
+
+_no_summary_string = "No summary for this constraint. " \
+                     "To generate one, pass a function as the parameter named " \
+                     '"summary" when creating the Constraint.'
 
 
 @dataclass(frozen=True, eq=False)
-class ConstraintWithDomains(Constraint[DesignPart], Generic[DesignPart]):
+class ConstraintWithDomains(Constraint[DesignPart], Generic[DesignPart]):  # noqa
     domains: Optional[Tuple[Domain, ...]] = None
     """
     Tuple of :any:`Domain`'s to check; if not specified, all :any:`Domain`'s in :any:`Design` are checked.
     """
 
-    def generate_summary(self, design_part: DesignPart, report_only_violations: bool) -> str:
-        raise NotImplementedError('subclasses of ConstraintWithDomains must implement generate_summary')
-
 
 @dataclass(frozen=True, eq=False)
-class ConstraintWithStrands(Constraint[DesignPart], Generic[DesignPart]):
+class ConstraintWithStrands(Constraint[DesignPart], Generic[DesignPart]):  # noqa
     strands: Optional[Tuple[Strand, ...]] = None
     """
     Tuple of :any:`Strand`'s to check; if not specified, all :any:`Strand`'s in :any:`Design` are checked.
     """
 
-    def generate_summary(self, design_part: DesignPart, report_only_violations: bool) -> str:
-        raise NotImplementedError('subclasses of ConstraintWithStrands must implement generate_summary')
-
 
 @dataclass(frozen=True, eq=False)  # type: ignore
 class DomainConstraint(ConstraintWithDomains[Domain]):
-    """Constraint that applies to a single :any:`Domain`."""
+    """Constraint that applies to a single :any:`Domain`.
 
-    evaluate: Callable[[str, Optional[Domain]],
-                       float] = lambda _: 0.0
+    Specify :data:`Constraint._evaluate` in the constructor."""
 
-    summary: Callable[[Domain],
-                      str] = lambda _: _no_summary_string
-
-    def __call__(self, sequence: str, domain: Optional[Domain]) -> float:
-        self._check_sequence_only(domain)
-        excess = (self.evaluate)(sequence, domain)  # noqa
-        if excess < 0:
-            return 0.0
-        score = (self.score_transfer_function)(excess)  # noqa
-        return score
+    def __post_init__(self) -> None:
+        if self._evaluate is None:
+            raise ValueError('_evaluate must be specified for a DomainConstraint')
+        super().__post_init__()
 
-    def generate_summary(self, domain: Domain, report_only_violations: bool) -> str:
-        return (self.summary)(domain)  # noqa
+    @staticmethod
+    def part_name() -> str:
+        return 'domain'
 
 
 @dataclass(frozen=True, eq=False)  # type: ignore
 class StrandConstraint(ConstraintWithStrands[Strand]):
-    """Constraint that applies to a single :any:`Strand`."""
+    """Constraint that applies to a single :any:`Strand`.
 
-    evaluate: Callable[[str, Optional[Strand]],
-                       float] = lambda _: 0.0
+    Specify :data:`Constraint._evaluate` in the constructor."""
 
-    summary: Callable[[Strand],
-                      str] = lambda _: _no_summary_string
-
-    def __call__(self, sequence: str, strand: Optional[Strand]) -> float:
-        self._check_sequence_only(strand)
-        excess = (self.evaluate)(sequence, strand)  # noqa
-        if excess < 0:
-            return 0.0
-        score = (self.score_transfer_function)(excess)  # noqa
-        return score
+    def __post_init__(self) -> None:
+        if self._evaluate is None:
+            raise ValueError('_evaluate must be specified for a StrandConstraint')
+        super().__post_init__()
 
-    def generate_summary(self, strand: Strand, report_only_violations: bool) -> str:
-        return self.summary(strand)  # noqa
+    @staticmethod
+    def part_name() -> str:
+        return 'strand'
 
 
 @dataclass(frozen=True, eq=False)
-class ConstraintWithDomainPairs(Constraint[DesignPart], Generic[DesignPart]):
+class ConstraintWithDomainPairs(Constraint[DesignPart], Generic[DesignPart]):  # noqa
     pairs: Optional[Tuple[Tuple[Domain, Domain], ...]] = None
     """
     List of pairs of :any:`Domain`'s to check; if not specified, all pairs in :any:`Design` are checked.
     """
 
-    def generate_summary(self, design_part: DesignPart, report_only_violations: bool) -> str:
-        raise NotImplementedError('subclasses of ConstraintWithStrandPairs must implement generate_summary')
-
 
 @dataclass(frozen=True, eq=False)
-class ConstraintWithStrandPairs(Constraint[DesignPart], Generic[DesignPart]):
+class ConstraintWithStrandPairs(Constraint[DesignPart], Generic[DesignPart]):  # noqa
     pairs: Optional[Tuple[Tuple[Strand, Strand], ...]] = None
     """
     List of pairs of :any:`Strand`'s to check; if not specified, all pairs in :any:`Design` are checked.
     """
 
-    def generate_summary(self, design_part: DesignPart, report_only_violations: bool) -> str:
-        raise NotImplementedError('subclasses of ConstraintWithStrandPairs must implement generate_summary')
-
 
 @dataclass(frozen=True, eq=False)  # type: ignore
 class DomainPairConstraint(ConstraintWithDomainPairs[Tuple[Domain, Domain]]):
-    """Constraint that applies to a pair of :any:`Domain`'s."""
-
-    evaluate: Callable[[str, str, Optional[Domain], Optional[Domain]],
-                       float] = lambda _, __: 0.0
-    """
-    Pairwise check to perform on :any:`Domain`'s.
-    Returns True if and only if the pair satisfies the constraint.
-    """
+    """Constraint that applies to a pair of :any:`Domain`'s.
 
-    summary: Callable[[Domain, Domain],
-                      str] = lambda _, __: _no_summary_string
+    Specify :data:`Constraint._evaluate` in the constructor."""
 
-    def __call__(self, sequence1: str, sequence2: str,
-                 domain1: Optional[Domain], domain2: Optional[Domain]) -> float:
-        excess = (self.evaluate)(sequence1, sequence2, domain1, domain2)  # noqa
-        if excess < 0:
-            return 0.0
-        score = (self.score_transfer_function)(excess)  # noqa
-        return score
+    def __post_init__(self) -> None:
+        if self._evaluate is None:
+            raise ValueError('_evaluate must be specified for a DomainPairConstraint')
+        super().__post_init__()
 
-    def generate_summary(self, domain_pair: Tuple[Domain, Domain], report_only_violations: bool) -> str:
-        domain1, domain2 = domain_pair
-        return (self.summary)(domain1, domain2)  # noqa
+    @staticmethod
+    def part_name() -> str:
+        return 'domain pair'
 
 
 @dataclass(frozen=True, eq=False)  # type: ignore
 class StrandPairConstraint(ConstraintWithStrandPairs[Tuple[Strand, Strand]]):
-    """Constraint that applies to a pair of :any:`Strand`'s."""
-
-    evaluate: Callable[[str, str, Optional[Strand], Optional[Strand]],
-                       float] = lambda _, __: 0.0
-    """
-    Pairwise check to perform on :any:`Strand`'s.
-    Returns float indicating how much the constraint is violated,
-    or 0.0 if the constraint is satisfied.
-    """
-
-    summary: Callable[[Strand, Strand],
-                      str] = lambda _, __: _no_summary_string
-
-    def __call__(self, sequence1: str, sequence2: str,
-                 strand1: Optional[Strand], strand2: Optional[Strand]) -> float:
-        self._check_sequence_only(strand1)
-        self._check_sequence_only(strand2)
-        excess = (self.evaluate)(sequence1, sequence2, strand1, strand2)  # noqa
-        if excess < 0:
-            return 0.0
-        score = (self.score_transfer_function)(excess)  # noqa
-        return score
-
-    def generate_summary(self, strand_pair: Tuple[Strand, Strand], report_only_violations: bool) -> str:
-        strand1, strand2 = strand_pair
-        return (self.summary)(strand1, strand2)  # noqa
-
-
-@dataclass(frozen=True, eq=False)  # type: ignore
-class DomainPairsConstraint(ConstraintWithDomainPairs[Iterable[Tuple[Domain, Domain]]]):
-    """
-    Similar to :any:`DomainsConstraint` but operates on a specified list of pairs of :any:`Domain`'s.
-    """
-
-    evaluate: Callable[[Iterable[Tuple[Domain, Domain]]],
-                       List[Tuple[OrderedSet[Domain], float]]] = lambda _: []
-    """
-    Pairwise check to perform on :any:`Domain`'s.
-    Returns True if and only if the all pairs in the input iterable satisfy the constraint.
-    """
+    """Constraint that applies to a pair of :any:`Strand`'s.
 
-    summary: Callable[[Iterable[Tuple[Domain, Domain]], bool],
-                      ConstraintReport] = lambda _: _no_summary_string
+    Specify :data:`Constraint._evaluate` in the constructor."""
 
-    def __call__(self, domain_pairs: Iterable[Tuple[Domain, Domain]]) \
-            -> List[Tuple[OrderedSet[Domain], float]]:
-        sets_excesses = (self.evaluate)(domain_pairs)  # noqa
-        sets_scores = _alter_scores_by_transfer(sets_excesses, self.score_transfer_function)
-        return sets_scores
+    def __post_init__(self) -> None:
+        if self._evaluate is None:
+            raise ValueError('_evaluate must be specified for a StrandPairConstraint')
+        super().__post_init__()
 
-    def generate_summary(self, domain_pairs: Iterable[Tuple[Domain, Domain]],
-                         report_only_violations: bool) -> ConstraintReport:
-        return (self.summary)(domain_pairs, report_only_violations)  # noqa
+    @staticmethod
+    def part_name() -> str:
+        return 'strand pair'
 
 
 @dataclass(frozen=True, eq=False)  # type: ignore
-class StrandPairsConstraint(ConstraintWithStrandPairs[Iterable[Tuple[Strand, Strand]]]):
-    """
-    Similar to :any:`StrandsConstraint` but operates on a specified list of pairs of :any:`Strand`'s.
-    """
-
-    evaluate: Callable[[Iterable[Tuple[Strand, Strand]]],
-                       List[Tuple[OrderedSet[Domain], float]]] = lambda _: []
-    """
-    Pairwise check to perform on :any:`Strand`'s.
-    Returns True if and only if the all pairs in the input iterable satisfy the constraint.
+class DomainsConstraint(ConstraintWithDomains[Iterable[Domain]]):
     """
+    Constraint that applies to a several :any:`Domain`'s.
 
-    summary: Callable[[Iterable[Tuple[Strand, Strand]], bool],
-                      ConstraintReport] = lambda _: _no_summary_string
-
-    def __call__(self, strand_pairs: Iterable[Tuple[Strand, Strand]]) \
-            -> List[Tuple[OrderedSet[Domain], float]]:
-        sets_excesses = (self.evaluate)(strand_pairs)  # noqa
-        sets_scores = _alter_scores_by_transfer(sets_excesses, self.score_transfer_function)
-        return sets_scores
+    Specify :data:`Constraint._evaluate_bulk` in the constructor.
 
-    def generate_summary(self, strand_pairs: Iterable[Tuple[Strand, Strand]],
-                         report_only_violations: bool) -> ConstraintReport:
-        return (self.summary)(strand_pairs, report_only_violations)  # noqa
-
-
-@dataclass(frozen=True, eq=False)  # type: ignore
-class DomainsConstraint(ConstraintWithDomains[Iterable[Domain]]):
-    """
-    Constraint that applies to a several :any:`Domain`'s. The difference with :any:`DomainConstraint` is that
+    The difference with :any:`DomainConstraint` is that
     the caller may want to process all :any:`Domain`'s at once, e.g., by giving many of them to a third-party
     program such as ViennaRNA, which may be more efficient than repeatedly calling a Python function.
 
@@ -3866,33 +3691,26 @@ class DomainsConstraint(ConstraintWithDomains[Iterable[Domain]]):
     initial violations of constraints, subsequent calls to this constraint only give the domain that was
     mutated, not the entire of :any:`Domain`'s in the whole :any:`Design`.
     Use :any:`DesignConstraint` for constraints that require every :any:`Domain` in the :any:`Design`.
-
-    Return value is a list of sets of :any:`Domain`'s. Each element of the list corresponds to one violation
-    of the :any:`DomainsConstraint`. The search will assign to a :any:`Domain` `d` a weight of
-    :py:data:`Constraint.weight` once for each set in this list that contains `d`. For example, if
-    :py:data:`Constraint.weight` is 1.0, the the return value is ``[{d1, d2}, {d2, d3}]``, then
-    ``d1`` and ``d3`` are assigned weight 1.0, and ``d2`` is assigned weight 2.0.
     """
 
-    evaluate: Callable[[Iterable[Domain]],
-                       List[Tuple[OrderedSet[Domain], float]]] = lambda _: []
-
-    summary: Callable[[Iterable[Domain], bool],
-                      ConstraintReport] = lambda _: _no_summary_string
-
-    def __call__(self, domains: Iterable[Domain]) -> List[Tuple[OrderedSet[Domain], float]]:
-        sets_excesses = (self.evaluate)(domains)  # noqa
-        sets_scores = _alter_scores_by_transfer(sets_excesses, self.score_transfer_function)
-        return sets_scores
+    def __post_init__(self) -> None:
+        if self._evaluate_bulk is None:
+            raise ValueError('_evaluate_bulk must be specified for a DomainsConstraint')
+        super().__post_init__()
 
-    def generate_summary(self, domains: Iterable[Domain], report_only_violations: bool) -> ConstraintReport:
-        return (self.summary)(domains, report_only_violations)  # noqa
+    @staticmethod
+    def part_name() -> str:
+        return 'domain'
 
 
 @dataclass(frozen=True, eq=False)  # type: ignore
 class StrandsConstraint(ConstraintWithStrands[Iterable[Strand]]):
     """
-    Constraint that applies to a several :any:`Strand`'s. The difference with :any:`StrandConstraint` is that
+    Constraint that applies to a several :any:`Strand`'s.
+
+    Specify :data:`Constraint._evaluate_bulk` in the constructor.
+
+    The difference with :any:`StrandConstraint` is that
     the caller may want to process all :any:`Strand`'s at once, e.g., by giving many of them to a third-party
     program such as ViennaRNA.
 
@@ -3900,27 +3718,52 @@ class StrandsConstraint(ConstraintWithStrands[Iterable[Strand]]):
     initial violations of constraints, subsequent calls to this constraint only give strands containing
     the domain that was mutated, not the entire of :any:`Strand`'s in the whole :any:`Design`.
     Use :any:`DesignConstraint` for constraints that require every :any:`Strand` in the :any:`Design`.
+    """
 
-    Return value is a list of sets of :any:`Domain`'s. Each element of the list corresponds to one violation
-    of the :any:`StrandsConstraint`. The search will assign to a :any:`Domain` `d` a weight of
-    :py:data:`Constraint.weight` once for each set in this list that contains `d`. For example, if
-    :py:data:`Constraint.weight` is 1.0, the the return value is ``[{d1, d2}, {d2, d3}]``, then
-    ``d1`` and ``d3`` are assigned weight 1.0, and ``d2`` is assigned weight 2.0.
+    def __post_init__(self) -> None:
+        if self._evaluate_bulk is None:
+            raise ValueError('_evaluate_bulk must be specified for a StrandsConstraint')
+        super().__post_init__()
+
+    @staticmethod
+    def part_name() -> str:
+        return 'strand'
+
+
+@dataclass(frozen=True, eq=False)  # type: ignore
+class DomainPairsConstraint(ConstraintWithDomainPairs[Tuple[Domain, Domain]]):
     """
+    Similar to :any:`DomainsConstraint` but operates on a specified list of pairs of :any:`Domain`'s.
 
-    evaluate: Callable[[Iterable[Strand]],
-                       List[Tuple[OrderedSet[Domain], float]]] = lambda _: []
+    Specify :data:`Constraint._evaluate_bulk` in the constructor.
+    """
 
-    summary: Callable[[Iterable[Strand], bool],
-                      ConstraintReport] = lambda _: _no_summary_string
+    def __post_init__(self) -> None:
+        if self._evaluate_bulk is None:
+            raise ValueError('_evaluate_bulk must be specified for a DomainPairsConstraint')
+        super().__post_init__()
 
-    def __call__(self, strands: Iterable[Strand]) -> List[Tuple[OrderedSet[Domain], float]]:
-        sets_excesses = (self.evaluate)(strands)  # noqa
-        sets_scores = _alter_scores_by_transfer(sets_excesses, self.score_transfer_function)
-        return sets_scores
+    @staticmethod
+    def part_name() -> str:
+        return 'domain pair'
+
+
+@dataclass(frozen=True, eq=False)  # type: ignore
+class StrandPairsConstraint(ConstraintWithStrandPairs[Iterable[Tuple[Strand, Strand]]]):
+    """
+    Similar to :any:`StrandsConstraint` but operates on a specified list of pairs of :any:`Strand`'s.
+
+    Specify :data:`Constraint._evaluate_bulk` in the constructor.
+    """
 
-    def generate_summary(self, strands: Iterable[Strand], report_only_violations: bool) -> ConstraintReport:
-        return (self.summary)(strands, report_only_violations)  # noqa
+    def __post_init__(self) -> None:
+        if self._evaluate_bulk is None:
+            raise ValueError('_evaluate_bulk must be specified for a StrandPairsConstraint')
+        super().__post_init__()
+
+    @staticmethod
+    def part_name() -> str:
+        return 'strand pair'
 
 
 @dataclass(frozen=True, eq=False)  # type: ignore
@@ -3929,33 +3772,53 @@ class DesignConstraint(Constraint[Design]):
     Constraint that applies to the entire :any:`Design`. This is used for any :any:`Constraint` that
     does not naturally fit the structure of the other types of constraints.
 
-    There is an optional parameter `domain_changed`, defaulting to None; if specified, the constraint
-    can restrict its check on the assumption that only `domain_changed` has changed since the last
-    time the violations were collected. In other words, it need only return sets of :any:`Domain`'s involving
-    violations of this :any:`DesignConstraint` that involve `domain_changed`. For example, if it checks
-    all triples of :any:`Domain`'s, there is no need to check any triple not containing `domain_changed`.
-
-    Return value is a list of sets of :any:`Domain`'s. Each element of the list corresponds to one violation
-    of the :any:`DomainsConstraint`. The search will assign to a :any:`Domain` `d` a weight of
-    :py:data:`Constraint.weight` once for each set in this list that contains `d`. For example, if
-    :py:data:`Constraint.weight` is 1.0, the the return value is ``[{d1, d2}, {d2, d3}]``, then
-    ``d1`` and ``d3`` are assigned weight 1.0, and ``d2`` is assigned weight 2.0.
+    Unlike other constraints, which specify either :data:`Constraint._evaluate` or
+    :data:`Constraint._evaluate_bulk`, a :any:`DesignConstraint` leaves both of these unspecified and
+    specifies :data:`DesignConstraint._evaluate_design` instead.
     """
 
-    evaluate: Callable[[Design, Optional[Domain]],
-                       List[Tuple[OrderedSet[Domain], float]]] = lambda _, __: []
-
-    summary: Callable[[Design, bool],
-                      ConstraintReport] = lambda _: _no_summary_string
+    _evaluate_design: Callable[[Design, Iterable[Domain]],
+                               List[Tuple[DesignPart, float, str]]] = lambda _: []
+    """
+    Evaluates the :any:`Design` (first argument), possibly taking into account which :any:`Domain`(s) have
+    changed in the last iteration (second argument).
+    
+    Returns a list of tuples (`part`, `score`, `summary`), 
+    one tuple per violation of the :any:`DesignConstraint`.
+    
+    `part` is the part of the :any:`Design` that caused the violation.
+    It must be one of :any:`Domain`, :any:`Strand`, pair of `Domain`'s, or tuple of :any:`Strand`'s.
+    
+    `score` is the score of the violation.
+    
+    `summary` is a 1-line summary of the violation to put into the generated reports.
+    """
 
-    def __call__(self, design: Design, domains_changed: Optional[Iterable[Domain]]) \
-            -> List[Tuple[OrderedSet[Domain], float]]:
-        sets_excesses = (self.evaluate)(design, domains_changed)  # noqa
-        sets_scores = _alter_scores_by_transfer(sets_excesses, self.score_transfer_function)
-        return sets_scores
+    def __post_init__(self) -> None:
+        if self._evaluate is not None:
+            raise ValueError('_evaluate should be None in a DesignConstraint')
+        if self._evaluate_bulk is not None:
+            raise ValueError('_evaluate_bulk should be None in a DesignConstraint')
+        if self._evaluate_design is None:
+            raise ValueError('_evaluate_design should be specified in a DesignConstraint')
+
+    def evaluate_design(self, design: Design, domains_changed: Iterable[Domain]) \
+            -> List[Tuple[DesignPart, float, str]]:
+        if self._evaluate_bulk is None:
+            raise ValueError('Cannot call evaluate_bulk on a Constraint unless _evaluate_bulk is specified')
+
+        results = (self._evaluate_bulk)(design, domains_changed)  # noqa
+        results_with_scores_transferred_and_weighted = []
+        for part, excess, summary in results:
+            if excess < 0.0:
+                excess = 0.0
+            weighted_score = self.weight * self.score_transfer_function(excess)
+            results_with_scores_transferred_and_weighted.append((part, weighted_score, summary))
+        return results_with_scores_transferred_and_weighted
 
-    def generate_summary(self, design: Design, report_only_violations: bool) -> ConstraintReport:
-        return (self.summary)(design, report_only_violations)  # noqa
+    @staticmethod
+    def part_name() -> str:
+        return 'whole design'
 
 
 def verify_designs_match(design1: Design, design2: Design, check_fixed: bool = True) -> None:
@@ -4021,14 +3884,83 @@ def convert_threshold(threshold: Union[float, Dict[T, float]], key: T) -> float:
     return threshold_value
 
 
-def nupack_strand_secondary_structure_constraint(
+def nupack_domain_complex_free_energy_constraint(
+        threshold: float,
+        temperature: float = dv.default_temperature,
+        sodium: float = dv.default_sodium,
+        magnesium: float = dv.default_magnesium,
+        weight: float = 1.0,
+        score_transfer_function: Callable[[float], float] = default_score_transfer_function,
+        parallel: bool = False,
+        description: Optional[str] = None,
+        short_description: str = 'strand_ss_nupack',
+        domains: Optional[Iterable[Domain]] = None) -> DomainConstraint:
+    """
+    Returns constraint that checks individual :any:`Domain`'s for excessive interaction using
+    NUPACK's pfunc.
+
+    NUPACK 4 must be installed. Installation instructions can be found at
+    https://piercelab-caltech.github.io/nupack-docs/start/.
+
+    :param threshold:
+        energy threshold in kcal/mol
+    :param temperature:
+        temperature in Celsius
+    :param sodium:
+        molarity of sodium (more generally, monovalent ions such as Na+, K+, NH4+)
+        in moles per liter
+    :param magnesium:
+        molarity of magnesium (Mg++) in moles per liter
+    :param weight:
+        how much to weigh this :any:`Constraint`
+    :param score_transfer_function:
+        See :py:data:`Constraint.score_transfer_function`.
+    :param parallel:
+        Whether to use threadds to parallelize.
+    :param domains:
+        :any:`Domain`'s to check; if not specified, all domains are checked.
+    :param description:
+        detailed description of constraint suitable for putting in report; if not specified
+        a reasonable default is chosen
+    :param short_description:
+        short description of constraint suitable for logging to stdout
+    :return:
+        the constraint
+    """
+
+    def evaluate(seqs: Tuple[str], domain: Optional[Domain]) -> Tuple[float, str]:
+        sequence = seqs[0]
+        energy = dv.complex_free_energy_single_strand(sequence, temperature, sodium, magnesium)
+        if domain is not None:
+            logger.debug(
+                f'domain ss threshold: {threshold:6.2f} '
+                f'secondary_structure_single_strand({domain.name, temperature}) = {energy:6.2f} ')
+        excess = threshold - energy
+        return max(0.0, excess), f'{energy:6.2f} kcal/mol'
+
+    if description is None:
+        description = f'NUPACK secondary structure of domain exceeds {threshold} kcal/mol'
+
+    if domains is not None:
+        domains = tuple(domains)
+
+    return DomainConstraint(description=description,
+                            short_description=short_description,
+                            weight=weight,
+                            score_transfer_function=score_transfer_function,
+                            _evaluate=evaluate,
+                            parallel=parallel,
+                            domains=domains)
+
+
+def nupack_strand_complex_free_energy_constraint(
         threshold: float,
         temperature: float = dv.default_temperature,
         sodium: float = dv.default_sodium,
         magnesium: float = dv.default_magnesium,
         weight: float = 1.0,
         score_transfer_function: Callable[[float], float] = default_score_transfer_function,
-        threaded: bool = False,
+        parallel: bool = False,
         description: Optional[str] = None,
         short_description: str = 'strand_ss_nupack',
         strands: Optional[Iterable[Strand]] = None) -> StrandConstraint:
@@ -4052,7 +3984,7 @@ def nupack_strand_secondary_structure_constraint(
         how much to weigh this :any:`Constraint`
     :param score_transfer_function:
         See :py:data:`Constraint.score_transfer_function`.
-    :param threaded:
+    :param parallel:
         Whether to use threadds to parallelize.
     :param strands:
         Strands to check; if not specified, all strands are checked.
@@ -4065,19 +3997,20 @@ def nupack_strand_secondary_structure_constraint(
         the constraint
     """
 
-    def evaluate(sequence: str, strand: Optional[Strand]) -> float:
-        energy = dv.secondary_structure_single_strand(sequence, temperature, sodium, magnesium)
+    def evaluate(seqs: Tuple[str], strand: Optional[Strand]) -> Tuple[float, str]:
+        sequence = seqs[0]
+        energy = dv.complex_free_energy_single_strand(sequence, temperature, sodium, magnesium)
         if strand is not None:
             logger.debug(
                 f'strand ss threshold: {threshold:6.2f} '
                 f'secondary_structure_single_strand({strand.name, temperature}) = {energy:6.2f} ')
         excess = threshold - energy
-        return max(0.0, excess)
+        return max(0.0, excess), f'{energy:6.2f} kcal/mol'
 
-    def summary(strand: Strand) -> str:
-        sequence = strand.sequence()
-        energy = dv.secondary_structure_single_strand(sequence, temperature, sodium, magnesium)
-        return f'{energy:6.2f} kcal/mol'
+    # def summary(strand: Strand) -> str:
+    #     sequence = strand.sequence()
+    #     energy = dv.complex_free_energy_single_strand(sequence, temperature, sodium, magnesium)
+    #     return f'{energy:6.2f} kcal/mol'
 
     if description is None:
         description = f'NUPACK secondary structure of strand exceeds {threshold} kcal/mol'
@@ -4089,16 +4022,15 @@ def summary(strand: Strand) -> str:
                             short_description=short_description,
                             weight=weight,
                             score_transfer_function=score_transfer_function,
-                            evaluate=evaluate,
-                            threaded=threaded,
-                            strands=strands,
-                            summary=summary)
+                            _evaluate=evaluate,
+                            parallel=parallel,
+                            strands=strands)
 
 
 def nupack_domain_pair_constraint(
         threshold: float,
         temperature: float = dv.default_temperature,
-        threaded: bool = False,
+        parallel: bool = False,
         weight: float = 1.0,
         score_transfer_function: Callable[[float], float] = default_score_transfer_function,
         description: Optional[str] = None,
@@ -4114,9 +4046,9 @@ def nupack_domain_pair_constraint(
         Energy threshold in kcal/mol.
     :param temperature:
         Temperature in Celsius
-    :param threaded:
+    :param parallel:
         Whether to test the each pair of :any:`Domain`'s in parallel (i.e., sets field
-        :py:data:`Constraint.threaded`)
+        :py:data:`Constraint.parallel`)
     :param weight:
         How much to weigh this :any:`Constraint`.
     :param score_transfer_function:
@@ -4148,11 +4080,13 @@ def nupack_domain_pair_constraint(
     def binding_closure(seq_pair: Tuple[str, str]) -> float:
         return dv.binding(seq_pair[0], seq_pair[1], temperature=temperature)
 
-    def evaluate(seq1: str, seq2: str, domain1: Optional[Domain], domain2: Optional[Domain]) -> float:
+    # def evaluate(seq1: str, seq2: str, domain1: Optional[Domain], domain2: Optional[Domain]) -> float:
+    def evaluate(seqs: Tuple[str, ...], domain_pair: Optional[DomainPair]) -> Tuple[float, str]:
+        seq1, seq2 = seqs
         name_pairs = [(None, None)] * 4
-        if domain1 is not None and domain2 is not None:
+        if domain_pair is not None:
             seq_pairs, name_pairs, _ = _all_pairs_domain_sequences_complements_names_from_domains(
-                [(domain1, domain2)])
+                [domain_pair])
         else:
             # If seq1==seq2, don't check d-d* or d*-d in this case, but do check d-d and d*-d*
             seq_pairs = [
@@ -4182,21 +4116,15 @@ def evaluate(seq1: str, seq2: str, domain1: Optional[Domain], domain2: Optional[
             excesses.append(excess)
 
         max_excess = max(excesses)
-        return max(0.0, max_excess)
 
-    def summary(domain1: Domain, domain2: Domain) -> str:
-        seq_pairs, domain_name_pairs, _ = _all_pairs_domain_sequences_complements_names_from_domains(
-            [(domain1, domain2)])
-        energies = []
-        for seq1, seq2 in seq_pairs:
-            energy = dv.binding(seq1, seq2, temperature=temperature)
-            energies.append(energy)
-        max_name_length = max(len(name) for name in _flatten(domain_name_pairs))
+        max_name_length = max(len(name) for name in _flatten(name_pairs))
         lines = [f'{name1:{max_name_length}}, '
                  f'{name2:{max_name_length}}: '
                  f' {energy:6.2f} kcal/mol'
-                 for (name1, name2), energy in zip(domain_name_pairs, energies)]
-        return '\n  ' + '\n  '.join(lines)
+                 for (name1, name2), energy in zip(name_pairs, energies)]
+        msg = '\n  ' + '\n  '.join(lines)
+
+        return max(0.0, max_excess), msg
 
     if pairs is not None:
         pairs = tuple(pairs)
@@ -4205,9 +4133,8 @@ def summary(domain1: Domain, domain2: Domain) -> str:
                                 short_description=short_description,
                                 weight=weight,
                                 score_transfer_function=score_transfer_function,
-                                evaluate=evaluate,
-                                summary=summary,
-                                threaded=threaded,
+                                _evaluate=evaluate,
+                                parallel=parallel,
                                 pairs=pairs)
 
 
@@ -4220,7 +4147,7 @@ def nupack_strand_pair_constraint(
         score_transfer_function: Callable[[float], float] = default_score_transfer_function,
         description: Optional[str] = None,
         short_description: str = 'strand_pair_nupack',
-        threaded: bool = False,
+        parallel: bool = False,
         pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
 ) -> StrandPairConstraint:
     """
@@ -4242,7 +4169,7 @@ def nupack_strand_pair_constraint(
         How much to weigh this :any:`Constraint`.
     :param score_transfer_function:
         See :py:data:`Constraint.score_transfer_function`.
-    :param threaded:
+    :param parallel:
         Whether to use threading to parallelize evaluating this constraint.
     :param description:
         Detailed description of constraint suitable for report.
@@ -4258,21 +4185,25 @@ def nupack_strand_pair_constraint(
     if description is None:
         description = f'NUPACK binding energy of strand pair exceeds {threshold} kcal/mol'
 
-    def evaluate(sequence1: str, sequence2: str,
-                 strand1: Optional[Strand], strand2: Optional[Strand]) -> float:
+    # def evaluate(sequence1: str, sequence2: str,
+    #              strand1: Optional[Strand], strand2: Optional[Strand]) -> float:
+    def evaluate(seqs: Tuple[str, ...],
+                 strand_pair: Optional[StrandPair]) -> Tuple[float, str]:
+        sequence1, sequence2 = seqs
         energy = dv.binding(sequence1, sequence2, temperature=temperature,
                             sodium=sodium, magnesium=magnesium)
-        if strand1 is not None and strand2 is not None:
+        if strand_pair is not None:
+            strand1, strand2 = strand_pair.individual_parts()
             logger.debug(
                 f'strand pair threshold: {threshold:6.2f} '
                 f'binding({strand1.name, strand2.name, temperature}) = {energy:6.2f} ')
         excess = threshold - energy
-        return max(0.0, excess)
+        return max(0.0, excess), f'{energy:6.2f} kcal/mol'
 
-    def summary(strand1: Strand, strand2: Strand) -> str:
-        energy = dv.binding(strand1.sequence(), strand2.sequence(), temperature=temperature,
-                            sodium=sodium, magnesium=magnesium)
-        return f'{energy:6.2f} kcal/mol'
+    # def summary(strand1: Strand, strand2: Strand) -> str:
+    #     energy = dv.binding(strand1.sequence(), strand2.sequence(), temperature=temperature,
+    #                         sodium=sodium, magnesium=magnesium)
+    #     return f'{energy:6.2f} kcal/mol'
 
     if pairs is not None:
         pairs = tuple(pairs)
@@ -4281,10 +4212,11 @@ def summary(strand1: Strand, strand2: Strand) -> str:
                                 short_description=short_description,
                                 weight=weight,
                                 score_transfer_function=score_transfer_function,
-                                threaded=threaded,
+                                parallel=parallel,
                                 pairs=pairs,
-                                evaluate=evaluate,
-                                summary=summary)
+                                _evaluate=evaluate,
+                                # summary=summary
+                                )
 
 
 def chunker(sequence: Sequence[T], chunk_length: Optional[int] = None, num_chunks: Optional[int] = None) -> \
@@ -4353,6 +4285,71 @@ def cpu_count(logical: bool = False) -> int:
     return count
 
 
+def rna_duplex_domain_pairs_constraint(
+        threshold: float,
+        temperature: float = dv.default_temperature,
+        weight: float = 1.0,
+        score_transfer_function: Callable[[float], float] = lambda x: x,
+        description: Optional[str] = None,
+        short_description: str = 'rna_dup_dom_pairs',
+        pairs: Optional[Iterable[Tuple[Domain, Domain]]] = None,
+        parameters_filename: str = dv.default_vienna_rna_parameter_filename) \
+        -> DomainPairsConstraint:
+    """
+    Returns constraint that checks given pairs of :any:`Domain`'s for excessive interaction using
+    Vienna RNA's RNAduplex executable.
+
+    :param threshold:
+        energy threshold
+    :param temperature:
+        temperature in Celsius
+    :param weight:
+        how much to weigh this :any:`Constraint`
+    :param score_transfer_function:
+        See :py:data:`Constraint.score_transfer_function`.
+    :param description:
+        long description of constraint suitable for printing in report file
+    :param short_description:
+        short description of constraint suitable for logging to stdout
+    :param pairs:
+        pairs of :any:`Domain`'s to compare; if not specified, checks all pairs
+    :param parameters_filename:
+        name of parameters file for ViennaRNA; default is
+        same as :py:meth:`vienna_nupack.rna_duplex_multiple`
+    :return:
+        constraint
+    """
+
+    if description is None:
+        description = f'RNAduplex energy for some domain pairs exceeds {threshold} kcal/mol'
+
+    def evaluate(domain_pairs: Iterable[DomainPair]) -> List[Tuple[DomainPair, float, str]]:
+        # if any(pair.domain1.sequence is None or pair.domain2.sequence is None for pair in domain_pairs):
+        #     raise ValueError('cannot evaluate domains unless they have sequences assigned')
+        sequence_pairs, _, _ = _all_pairs_domain_sequences_complements_names_from_domains(domain_pairs)
+        pairs_scores_summaries: List[Tuple[DomainPair, float, str]] = []
+        energies = dv.rna_duplex_multiple(sequence_pairs, logger, temperature, parameters_filename)
+
+        for pair, energy in zip(domain_pairs, energies):
+            excess = threshold - energy
+            if excess > 0.0:
+                summary = f'{energy:6.2f} kcal/mol'
+                pair_score_summary = (pair, excess, summary)
+                pairs_scores_summaries.append(pair_score_summary)
+        return pairs_scores_summaries
+
+    pairs_tuple = None
+    if pairs is not None:
+        pairs_tuple = tuple(pairs)
+
+    return DomainPairsConstraint(description=description,
+                                 short_description=short_description,
+                                 weight=weight,
+                                 score_transfer_function=score_transfer_function,
+                                 _evaluate_bulk=evaluate,
+                                 pairs=pairs_tuple)
+
+
 def rna_duplex_strand_pairs_constraint(
         threshold: float,
         temperature: float = dv.default_temperature,
@@ -4360,7 +4357,7 @@ def rna_duplex_strand_pairs_constraint(
         score_transfer_function: Callable[[float], float] = default_score_transfer_function,
         description: Optional[str] = None,
         short_description: str = 'rna_dup_strand_pairs',
-        threaded: bool = False,
+        parallel: bool = False,
         pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
         parameters_filename: str = dv.default_vienna_rna_parameter_filename) \
         -> StrandPairsConstraint:
@@ -4380,7 +4377,7 @@ def rna_duplex_strand_pairs_constraint(
         Long description of constraint suitable for putting into constraint report.
     :param short_description:
         Short description of constraint suitable for logging to stdout.
-    :param threaded:
+    :param parallel:
         Whether to test the each pair of :any:`Strand`'s in parallel.
     :param pairs:
         Pairs of :any:`Strand`'s to compare; if not specified, checks all pairs.
@@ -4394,109 +4391,36 @@ def rna_duplex_strand_pairs_constraint(
     if description is None:
         description = f'RNAduplex energy for some strand pairs exceeds {threshold} kcal/mol'
 
-    from dsd.stopwatch import Stopwatch
-
     num_threads = cpu_count() - 1  # this seems to be slightly faster than using all cores
+
     # we use ThreadPool instead of pathos because we're farming this out to processes through
     # subprocess module anyway, no need for pathos to boot up separate processes or serialize through dill
     thread_pool = ThreadPool(processes=num_threads)
 
-    def calculate_energies_unthreaded(sequence_pairs: Sequence[Tuple[str, str]]) -> List[float]:
+    def calculate_energies_unparallel(sequence_pairs: Sequence[Tuple[str, str]]) -> List[float]:
         return dv.rna_duplex_multiple(sequence_pairs, logger, temperature, parameters_filename)
 
     def calculate_energies(sequence_pairs: Sequence[Tuple[str, str]]) -> List[float]:
-        if threaded and len(sequence_pairs) > 1:
+        if parallel and len(sequence_pairs) > 1:
             lists_of_sequence_pairs = chunker(sequence_pairs, num_chunks=num_threads)
-            lists_of_energies = thread_pool.map(calculate_energies_unthreaded, lists_of_sequence_pairs)
+            lists_of_energies = thread_pool.map(calculate_energies_unparallel, lists_of_sequence_pairs)
             energies = _flatten(lists_of_energies)
         else:
-            energies = calculate_energies_unthreaded(sequence_pairs)
+            energies = calculate_energies_unparallel(sequence_pairs)
         return energies
 
-    def evaluate(strand_pairs: Iterable[Tuple[Strand, Strand]]) -> List[Tuple[OrderedSet[Domain], float]]:
-        stopwatch: Optional[Stopwatch] = Stopwatch()  # noqa
-        # stopwatch = None  # uncomment to not log time
-
-        sequence_pairs = [(s1.sequence(), s2.sequence()) for s1, s2 in strand_pairs]
+    def evaluate(strand_pairs: Iterable[StrandPair]) -> List[Tuple[StrandPair, float, str]]:
+        sequence_pairs = [(pair.strand1.sequence(), pair.strand2.sequence()) for pair in strand_pairs]
         energies = calculate_energies(sequence_pairs)
+        pairs_scores_summaries: List[Tuple[StrandPair, float, str]] = []
 
-        domain_sets_weights: List[Tuple[OrderedSet[Domain], float]] = []
-        for (strand1, strand2), energy in zip(strand_pairs, energies):
-            excess = energy_excess(energy, threshold)
+        for pair, energy in zip(strand_pairs, energies):
+            excess = threshold - energy
             if excess > 0.0:
-                domain_set_weights = (
-                    OrderedSet(strand1.unfixed_domains() + strand2.unfixed_domains()), excess)
-                domain_sets_weights.append(domain_set_weights)
-
-        if stopwatch is not None:
-            stopwatch.stop()
-            logger.debug(f'*** rna_duplex_strand_pairs_constraint ***')
-            logger.debug(f'*   description: {description}')
-            logger.debug(f'*   evaluated {len(sequence_pairs)} pairs of strands')
-            logger.debug(f'*   total time to evaluate: {stopwatch}')
-            logger.debug(f'*   energies: {sorted(energies)}')
-
-        return domain_sets_weights
-
-    def summary(strand_pairs: Iterable[Tuple[Strand, Strand]],
-                report_only_violations: bool) -> ConstraintReport:
-        num_violations = 0
-        num_checks = 0
-
-        max_name_length = max(len(strand.name) for strand in _flatten(strand_pairs))
-
-        # distinguish between pairs in which both strands are fixed (so cannot remove violation)
-        # versus those pairs in which at least one element of the pair is unfixed
-        both_fixed_pairs = [(s1, s2) for s1, s2 in strand_pairs if s1.fixed and s2.fixed]
-        one_unfixed_pairs = [(s1, s2) for s1, s2 in strand_pairs if not (s1.fixed and s2.fixed)]
-
-        summaries = []
-
-        for pairs_to_check, header_name in [(one_unfixed_pairs, 'pairs with at least one unfixed'),
-                                            (both_fixed_pairs, 'pairs with both fixed')]:
-            if len(pairs_to_check) > 0:
-                report = _summary_of_pairs(pairs_to_check, report_only_violations, max_name_length)
-                summary_ = _small_header(header_name, "=") + f'\n{report.content}\n'
-                num_violations += report.num_violations
-                num_checks += report.num_checks
-                summaries.append(summary_)
-
-        content = ''.join(summaries)
-        report = ConstraintReport(constraint=None, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
-
-    def _summary_of_pairs(strand_pairs: Iterable[Tuple[Strand, Strand]], report_only_violations: bool,
-                          max_name_length: int) -> ConstraintReport:
-        sequence_pairs = [(s1.sequence(), s2.sequence()) for s1, s2 in strand_pairs]
-        energies = calculate_energies(sequence_pairs)
-
-        strand_pairs_energies = zip(strand_pairs, energies)
-
-        num_checks = len(energies)
-        num_violations = 0
-        lines_and_excesses: List[Tuple[str, float]] = []
-        for (strand1, strand2), energy in strand_pairs_energies:
-            excess = energy_excess(energy, threshold)
-            passed = excess <= 0.0
-            if not passed:
-                num_violations += 1
-            if not report_only_violations or (report_only_violations and not passed):
-                line = (f'strands '
-                        f'{strand1.name:{max_name_length}}, '
-                        f'{strand2.name:{max_name_length}}: '
-                        f'{energy:6.2f} kcal/mol'
-                        f'{"" if passed else "  **violation**"}')
-                lines_and_excesses.append((line, excess))
-
-        # put in descending order of excess
-        lines_and_excesses.sort(key=lambda line_and_excess: line_and_excess[1], reverse=True)
-
-        lines = (line for line, _ in lines_and_excesses)
-        content = '\n'.join(lines)
-        report = ConstraintReport(constraint=None, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
+                summary = f'{energy:6.2f} kcal/mol'
+                pair_score_summary = (pair, excess, summary)
+                pairs_scores_summaries.append(pair_score_summary)
+        return pairs_scores_summaries
 
     pairs_tuple = None
     if pairs is not None:
@@ -4506,8 +4430,7 @@ def _summary_of_pairs(strand_pairs: Iterable[Tuple[Strand, Strand]], report_only
                                  short_description=short_description,
                                  weight=weight,
                                  score_transfer_function=score_transfer_function,
-                                 evaluate=evaluate,
-                                 summary=summary,
+                                 _evaluate_bulk=evaluate,
                                  pairs=pairs_tuple)
 
 
@@ -4535,7 +4458,7 @@ def rna_cofold_strand_pairs_constraint(
         score_transfer_function: Callable[[float], float] = default_score_transfer_function,
         description: Optional[str] = None,
         short_description: str = 'rna_dup_strand_pairs',
-        threaded: bool = False,
+        parallel: bool = False,
         pairs: Optional[Iterable[Tuple[Strand, Strand]]] = None,
         parameters_filename: str = dv.default_vienna_rna_parameter_filename) \
         -> StrandPairsConstraint:
@@ -4555,7 +4478,7 @@ def rna_cofold_strand_pairs_constraint(
         Long description of constraint suitable for putting into constraint report.
     :param short_description:
         Short description of constraint suitable for logging to stdout.
-    :param threaded:
+    :param parallel:
         Whether to test the each pair of :any:`Strand`'s in parallel.
     :param pairs:
         Pairs of :any:`Strand`'s to compare; if not specified, checks all pairs.
@@ -4569,77 +4492,36 @@ def rna_cofold_strand_pairs_constraint(
     if description is None:
         description = f'RNAcofold energy for some strand pairs exceeds {threshold} kcal/mol'
 
-    from dsd.stopwatch import Stopwatch
-
     num_threads = cpu_count() - 1  # this seems to be slightly faster than using all cores
+
+    # we use ThreadPool instead of pathos because we're farming this out to processes through
+    # subprocess module anyway, no need for pathos to boot up separate processes or serialize through dill
     thread_pool = ThreadPool(processes=num_threads)
 
-    def calculate_energies_unthreaded(sequence_pairs: Sequence[Tuple[str, str]]) -> List[float]:
+    def calculate_energies_unparallel(sequence_pairs: Sequence[Tuple[str, str]]) -> List[float]:
         return dv.rna_cofold_multiple(sequence_pairs, logger, temperature, parameters_filename)
 
     def calculate_energies(sequence_pairs: Sequence[Tuple[str, str]]) -> List[float]:
-        if threaded and len(sequence_pairs) > 1:
+        if parallel and len(sequence_pairs) > 1:
             lists_of_sequence_pairs = chunker(sequence_pairs, num_chunks=num_threads)
-            lists_of_energies = thread_pool.map(calculate_energies_unthreaded, lists_of_sequence_pairs)
+            lists_of_energies = thread_pool.map(calculate_energies_unparallel, lists_of_sequence_pairs)
             energies = _flatten(lists_of_energies)
         else:
-            energies = calculate_energies_unthreaded(sequence_pairs)
+            energies = calculate_energies_unparallel(sequence_pairs)
         return energies
 
-    def evaluate(strand_pairs: Iterable[Tuple[Strand, Strand]]) -> List[Tuple[OrderedSet[Domain], float]]:
-        stopwatch: Optional[Stopwatch] = Stopwatch()  # noqa
-        # stopwatch = None  # uncomment to not log time
-
-        sequence_pairs = [(s1.sequence(), s2.sequence()) for s1, s2 in strand_pairs]
-        domain_sets_weights: List[Tuple[OrderedSet[Domain], float]] = []
-        # domain_sets_weights: List[Tuple[OrderedSet[Domain], float]] = []
-
+    def evaluate(strand_pairs: Iterable[StrandPair]) -> List[Tuple[StrandPair, float, str]]:
+        sequence_pairs = [(pair.strand1.sequence(), pair.strand2.sequence()) for pair in strand_pairs]
         energies = calculate_energies(sequence_pairs)
+        pairs_scores_summaries: List[Tuple[StrandPair, float, str]] = []
 
-        for (strand1, strand2), energy in zip(strand_pairs, energies):
-            excess = energy_excess(energy, threshold)
+        for pair, energy in zip(strand_pairs, energies):
+            excess = threshold - energy
             if excess > 0.0:
-                domain_set_weights = (OrderedSet(strand1.unfixed_domains() + strand2.unfixed_domains()),
-                                      excess)
-                domain_sets_weights.append(domain_set_weights)
-
-        if stopwatch is not None:
-            stopwatch.stop()
-            logger.debug(f'*** rna_cofold_strand_pairs_constraint ***')
-            logger.debug(f'*   description: {description}')
-            logger.debug(f'*   evaluated {len(sequence_pairs)} pairs of strands')
-            logger.debug(f'*   total time to evaluate: {stopwatch}')
-            logger.debug(f'*   energies: {sorted(energies)}')
-
-        return domain_sets_weights
-
-    def summary(strand_pairs: Iterable[Tuple[Strand, Strand]],
-                report_only_violations: bool) -> ConstraintReport:
-        sequence_pairs = [(s1.sequence(), s2.sequence()) for s1, s2 in strand_pairs]
-        energies = calculate_energies(sequence_pairs)
-        max_name_length = max(len(strand.name) for strand in _flatten(strand_pairs))
-        strand_pairs_energies = zip(strand_pairs, energies)
-
-        num_checks = len(energies)
-        num_violations = 0
-        lines: List[str] = []
-        for (strand1, strand2), energy in strand_pairs_energies:
-            passed = energy_excess(energy, threshold) <= 0.0
-            if not passed:
-                num_violations += 1
-            if not report_only_violations or (report_only_violations and not passed):
-                line = (f'strands '
-                        f'{strand1.name:{max_name_length}}, '
-                        f'{strand2.name:{max_name_length}}: '
-                        f'{energy:6.2f} kcal/mol'
-                        f'{"" if passed else "  **violation**"}')
-                lines.append(line)
-
-        if not report_only_violations:
-            lines.sort(key=lambda line_: ' **violation**' not in line_)  # put violations first
-
-        return ConstraintReport(constraint=None, content='\n'.join(lines),
-                                num_violations=num_violations, num_checks=num_checks)
+                summary = f'{energy:6.2f} kcal/mol'
+                pair_score_summary = (pair, excess, summary)
+                pairs_scores_summaries.append(pair_score_summary)
+        return pairs_scores_summaries
 
     pairs_tuple = None
     if pairs is not None:
@@ -4649,13 +4531,12 @@ def summary(strand_pairs: Iterable[Tuple[Strand, Strand]],
                                  short_description=short_description,
                                  weight=weight,
                                  score_transfer_function=score_transfer_function,
-                                 evaluate=evaluate,
-                                 summary=summary,
+                                 _evaluate_bulk=evaluate,
                                  pairs=pairs_tuple)
 
 
 def _all_pairs_domain_sequences_complements_names_from_domains(
-        domain_pairs: Iterable[Tuple[Domain, Domain]]) \
+        domain_pairs: Iterable[DomainPair]) \
         -> Tuple[List[Tuple[str, str]], List[Tuple[str, str]], List[Tuple[Domain, Domain]]]:
     """
     :param domain_pairs:
@@ -4670,7 +4551,8 @@ def _all_pairs_domain_sequences_complements_names_from_domains(
     sequence_pairs: List[Tuple[str, str]] = []
     names: List[Tuple[str, str]] = []
     domains: List[Tuple[Domain, Domain]] = []
-    for d1, d2 in domain_pairs:
+    for pair in domain_pairs:
+        d1, d2 = pair.individual_parts()
         if d1 == d2:
             # don't check d-d* or d*-d in this case, but do check d-d and d*-d*
             starred_each = [(False, False), (True, True)]
@@ -4692,144 +4574,33 @@ def _flatten(list_of_lists: Iterable[Iterable[T]]) -> List[T]:
     return list(itertools.chain.from_iterable(list_of_lists))
 
 
-def rna_duplex_domain_pairs_constraint(
-        threshold: float,
-        temperature: float = dv.default_temperature,
-        weight: float = 1.0,
-        score_transfer_function: Callable[[float], float] = lambda x: x,
-        description: Optional[str] = None,
-        short_description: str = 'rna_dup_dom_pairs',
-        pairs: Optional[Iterable[Tuple[Domain, Domain]]] = None,
-        parameters_filename: str = dv.default_vienna_rna_parameter_filename) \
-        -> DomainPairsConstraint:
-    """
-    Returns constraint that checks given pairs of :any:`Domain`'s for excessive interaction using
-    Vienna RNA's RNAduplex executable.
-
-    :param threshold:
-        energy threshold
-    :param temperature:
-        temperature in Celsius
-    :param weight:
-        how much to weigh this :any:`Constraint`
-    :param score_transfer_function:
-        See :py:data:`Constraint.score_transfer_function`.
-    :param description:
-        long description of constraint suitable for printing in report file
-    :param short_description:
-        short description of constraint suitable for logging to stdout
-    :param pairs:
-        pairs of :any:`Domain`'s to compare; if not specified, checks all pairs
-    :param parameters_filename:
-        name of parameters file for ViennaRNA; default is
-        same as :py:meth:`vienna_nupack.rna_duplex_multiple`
-    :return:
-        constraint
-    """
-
-    if description is None:
-        description = f'RNAduplex energy for some strand pairs exceeds {threshold} kcal/mol'
-
-    def evaluate(domain_pairs: Iterable[Tuple[Domain, Domain]]) -> List[Tuple[OrderedSet[Domain], float]]:
-        if any(d1.sequence is None or d2.sequence is None for d1, d2 in domain_pairs):
-            raise ValueError('cannot evaluate domains unless they have sequences assigned')
-        sequence_pairs, _, _ = _all_pairs_domain_sequences_complements_names_from_domains(domain_pairs)
-        domain_sets_weights: List[Tuple[OrderedSet[Domain], float]] = []
-        energies = dv.rna_duplex_multiple(sequence_pairs, logger, temperature, parameters_filename)
-        for (domain1, domain2), energy in zip(domain_pairs, energies):
-            excess = threshold - energy
-            if excess > 0.0:
-                domain_set_weights = (OrderedSet([domain1, domain2]), excess)
-                domain_sets_weights.append(domain_set_weights)
-        return domain_sets_weights
-
-    def summary(domain_pairs: Iterable[Tuple[Domain, Domain]],
-                report_only_violations: bool) -> ConstraintReport:
-        sequence_pairs, domain_name_pairs, domains = \
-            _all_pairs_domain_sequences_complements_names_from_domains(domain_pairs)
-        energies = dv.rna_duplex_multiple(sequence_pairs, logger, temperature, parameters_filename)
-        max_name_length = max(len(name) for name in _flatten(domain_name_pairs))
-
-        num_checks = len(energies)
-        num_violations = 0
-        lines: List[str] = []
-        for (domain1, domain2), (name1, name2), energy in zip(domains, domain_name_pairs, energies):
-            passed = energy_excess_domains(energy, threshold, domain1, domain2) <= 0.0
-            if not passed:
-                num_violations += 1
-            if not report_only_violations or (report_only_violations and not passed):
-                line = (f'domains '
-                        f'{name1:{max_name_length}}, '
-                        f'{name2:{max_name_length}}: '
-                        f'{energy:6.2f} kcal/mol'
-                        f'{"" if passed else "  **violation**"}')
-                lines.append(line)
-
-        if not report_only_violations:
-            lines.sort(key=lambda line_: ' **violation**' not in line_)  # put violations first
-
-        content = '\n'.join(lines)
-        report = ConstraintReport(constraint=None, content=content,
-                                  num_violations=num_violations, num_checks=num_checks)
-        return report
-
-    pairs_tuple = None
-    if pairs is not None:
-        pairs_tuple = tuple(pairs)
-
-    return DomainPairsConstraint(description=description,
-                                 short_description=short_description,
-                                 weight=weight,
-                                 score_transfer_function=score_transfer_function,
-                                 evaluate=evaluate,
-                                 summary=summary,
-                                 pairs=pairs_tuple)
-
-
 #########################################################################################
 # ComplexConstraints defined below here
 #########################################################################################
 
 
 @dataclass(frozen=True, eq=False)
-class ConstraintWithComplexes(Constraint[DesignPart], Generic[DesignPart]):
+class ConstraintWithComplexes(Constraint[DesignPart], Generic[DesignPart]):  # noqa
     complexes: Tuple[Complex, ...] = ()
     """
     List of complexes (tuples of :any:`Strand`'s) to check.
     """
 
-    def generate_summary(self, design_part: DesignPart, report_only_violations: bool) -> str:
-        raise NotImplementedError('subclasses of ConstraintWithStrandPairs must implement generate_summary')
-
 
 @dataclass(frozen=True, eq=False)  # type: ignore
 class ComplexConstraint(ConstraintWithComplexes[Complex]):
-    """Constraint that applies to a complex (tuple of :any:`Strand`'s).
+    """
+    Constraint that applies to a complex (tuple of :any:`Strand`'s).
+
+    Specify :data:`Constraint._evaluate` in the constructor.
 
     Unlike other types of :any:`Constraint`'s such as :any:`StrandConstraint` or :any:`StrandPairConstraint`,
     there is no default list of :any:`Complex`'s that a :any:`ComplexConstraint` is applied to. The list of
-    :any:`Complex`'s must be specified manually in the constructor."""
-
-    evaluate: Callable[[Complex],
-                       float] = lambda _, __: 0.0
+    :any:`Complex`'s must be specified manually in the constructor.
     """
-    Evaluation to perform on complex (tuple of :any:`Strand`'s).
-    Returns float indicating how much the constraint is violated,
-    or 0.0 if the constraint is satisfied.
-    """
-
-    summary: Callable[[Complex],
-                      str] = lambda _, __: _no_summary_string
 
-    def __call__(self, strand_complex: Complex) -> float:
-        excess = (self.evaluate)(strand_complex)  # noqa
-        if excess < 0:
-            return 0.0
-        score = (self.score_transfer_function)(excess)  # noqa
-        return score
-
-    def generate_summary(self, strand_complex: Complex, report_only_violations: bool) -> str:
-        return (self.summary)(strand_complex)  # noqa
+    def part_name(self) -> str:
+        return 'complex'
 
 
 def _alter_scores_by_transfer(sets_excesses: List[Tuple[OrderedSet[Domain], float]],
@@ -4852,25 +4623,8 @@ class ComplexesConstraint(ConstraintWithComplexes[Iterable[Complex]]):
     (tuples of :any:`Strand`'s).
     """
 
-    evaluate: Callable[[Iterable[Complex]],
-                       List[Tuple[OrderedSet[Domain], float]]] = lambda _: []
-    """
-    Check to perform on an iterable of complexes (tuples of :any:`Strand`'s).
-    Returns True if and only if the all complexes in the input iterable satisfy the constraint.
-    """
-
-    summary: Callable[[Iterable[Complex], bool],
-                      ConstraintReport] = lambda _: _no_summary_string
-
-    def __call__(self, complexes: Iterable[Complex]) \
-            -> List[Tuple[OrderedSet[Domain], float]]:
-        sets_excesses = (self.evaluate)(strands)  # noqa
-        sets_scores = _alter_scores_by_transfer(sets_excesses, self.score_transfer_function)
-        return sets_scores
-
-    def generate_summary(self, complexes: Iterable[Complex],
-                         report_only_violations: bool) -> ConstraintReport:
-        return (self.summary)(strands, report_only_violations)  # noqa
+    def part_name(self) -> str:
+        return 'complex'
 
 
 class _AdjacentDuplexType(Enum):
@@ -6304,8 +6058,8 @@ def _get_base_pair_domain_endpoints_to_check(
     addr_translation_table: Dict[StrandDomainAddress, List[StrandDomainAddress]] = {}
 
     # Need to convert strands into strands lowest level subdomains
-    leafify_strand_complex = tuple(
-        [_leafify_strand(strand, addr_translation_table) for strand in strand_complex])
+    leafify_strand_complex = Complex(tuple(
+        [_leafify_strand(strand, addr_translation_table) for strand in strand_complex]))
 
     new_nonimplicit_base_pairs = []
     if nonimplicit_base_pairs:
@@ -6455,7 +6209,7 @@ def __get_base_pair_domain_endpoints_to_check(
     return base_pair_domain_endpoints_to_check
 
 
-def nupack_complex_secondary_structure_constraint(
+def nupack_complex_base_pair_probability_constraint(
         strand_complexes: List[Complex],
         nonimplicit_base_pairs: Optional[Iterable[BoundDomains]] = None,
         all_base_pairs: Optional[Iterable[BoundDomains]] = None,
@@ -6472,7 +6226,7 @@ def nupack_complex_secondary_structure_constraint(
         score_transfer_function: Callable[[float], float] = default_score_transfer_function,
         description: Optional[str] = None,
         short_description: str = 'complex_secondary_structure_nupack',
-        threaded: bool = False,
+        parallel: bool = False,
 ) -> ComplexConstraint:
     """Returns constraint that checks given base pairs probabilities in tuples of :any:`Strand`'s
 
@@ -6573,9 +6327,9 @@ def nupack_complex_secondary_structure_constraint(
         See :py:data:`Constraint.short_description` defaults to 'complex_secondary_structure_nupack'
     :type short_description:
         str, optional
-    :param threaded:
+    :param parallel:
         **TODO**: Implement this
-    :type threaded:
+    :type parallel:
         bool, optional
     :raises ImportError:
         If NUPACK 4 is not installed.
@@ -6611,23 +6365,35 @@ def nupack_complex_secondary_structure_constraint(
 
     strand_complex_template = strand_complexes[0]
 
-    if type(strand_complex_template) is not tuple:
+    if not isinstance(strand_complex_template, Complex):
         raise ValueError(
-            f"First element in strand_complexes was not a tuple of Strands. Please provide a tuple of Strands.")
+            f"First element in strand_complexes was not a Complex of Strands. "
+            f"Please provide a Complex of Strands.")
 
     for strand in strand_complex_template:
         if type(strand) is not Strand:
             raise ValueError(f"Complex at index 0 contained non-Strand object: {type(strand)}")
 
+    for strand_complex in strand_complexes:
+        for strand in strand_complex:
+            for domain in strand.domains:
+                if not domain.has_length():
+                    raise ValueError(f'''\
+Domain {domain.name} has no length yet. To use 
+nupack_complex_secondary_structure_constraint, each Domain must have a length 
+assigned, either by assigning it a DomainPool first, or by setting the Domain 
+to have a fixed DNA sequence by calling domain.set_fixed_sequence.''')
+
     for idx in range(1, len(strand_complexes)):
         strand_complex = strand_complexes[idx]
-        if type(strand_complex) is not tuple:
+        if not isinstance(strand_complex, Complex):
             raise ValueError(
-                f"Element at index {idx} was not a tuple of Strands. Please provide a tuple of Strands.")
+                f"Element {strand_complex} at index {idx} is not a Complex of Strands. "
+                f"Please provide a Complex of Strands.")
         if len(strand_complex) != len(strand_complex_template):
             raise ValueError(
-                f"Inconsistent complex structures: Complex at index {idx} contained {len(strand_complex)} strands, "
-                f"but complex at index 0 contained {len(strand_complex_template)} strands.")
+                f"Inconsistent complex structures: Complex at index {idx} contained {len(strand_complex)} "
+                f"strands, but complex at index 0 contained {len(strand_complex_template)} strands.")
         for s in range(len(strand_complex)):
             other_strand: Strand = strand_complex[s]
             template_strand: Strand = strand_complex_template[s]
@@ -6636,16 +6402,17 @@ def nupack_complex_secondary_structure_constraint(
                     f"Complex at index {idx} contained non-Strand object at index {s}: {type(other_strand)}")
             if len(other_strand.domains) != len(template_strand.domains):
                 raise ValueError(
-                    f"Strand {other_strand} (index {s} of strand_complexes at index {idx}) does not match the provided template"
-                    f"({template_strand}). "
-                    f"Strand {other_strand} contains {len(other_strand.domains)} domains but template strand {template_strand} contains "
-                    f"{len(template_strand.domains)} domains.")
+                    f"Strand {other_strand} (index {s} of strand_complexes at index {idx}) does not match "
+                    f"the provided template ({template_strand}). "
+                    f"Strand {other_strand} contains {len(other_strand.domains)} domains but template "
+                    f"strand {template_strand} contains {len(template_strand.domains)} domains.")
             for d in range(1, len(other_strand.domains)):
                 domain_length: int = other_strand.domains[d].length
                 template_domain_length: int = template_strand.domains[d].length
                 if domain_length != template_domain_length:
                     raise ValueError(
-                        f"Strand {other_strand} (the strand at index {s} of the complex located at index {idx} of strand_complexes) does not match the "
+                        f"Strand {other_strand} (the strand at index {s} of the complex located at index "
+                        f"{idx} of strand_complexes) does not match the "
                         f"provided template ({template_strand}): domain at index {d} is length "
                         f"{domain_length}, but expected {template_domain_length}.")
 
@@ -6663,31 +6430,47 @@ def nupack_complex_secondary_structure_constraint(
     if description is None:
         description = 'Base pair probability of complex'
 
-    def evaluate(strand_complex_: Complex) -> float:
+    def evaluate(seqs: Tuple[str, ...], strand_complex_: Complex) -> Tuple[float, str]:
+        assert len(seqs) == len(strand_complex)
         bps = _violation_base_pairs(strand_complex_)
         err_sq = 0.0
+        # eval
         for bp in bps:
             e = base_type_probability_threshold[bp.base_pair_type] - bp.base_pairing_probability
             assert e > 0
             err_sq += e ** 2
-        return err_sq
+        # summary
+        if len(bps) == 0:
+            summary = "\tAll base pairs satisfy thresholds."
+        else:
+            summary_list = []
+            for bp in bps:
+                i = bp.base_index1
+                j = bp.base_index2
+                p = bp.base_pairing_probability
+                t = bp.base_pair_type
+                summary_list.append(
+                    f'\t{i},{j}: {math.floor(100 * p)}% '
+                    f'(<{round(100 * base_type_probability_threshold[t])}% [{t}])')
+            summary = '\n'.join(summary_list)
+        return err_sq, summary
 
     # summary would print all the base pairs
     # * indices of the bases e.g 2,7: 97.3% (<99%);  9,13: 75% (<80%); 1,7: 2.1% (>1%)
     # * maybe consider puting second and after base pairs on new line with indent
-    def summary(strand_complex_: Complex) -> str:
-        bps = _violation_base_pairs(strand_complex_)
-        if len(bps) == 0:
-            return "\tAll base pairs satisfy thresholds."
-        summary_list = []
-        for bp in bps:
-            i = bp.base_index1
-            j = bp.base_index2
-            p = bp.base_pairing_probability
-            t = bp.base_pair_type
-            summary_list.append(
-                f'\t{i},{j}: {math.floor(100 * p)}% (<{round(100 * base_type_probability_threshold[t])}% [{t}])')
-        return '\n'.join(summary_list)
+    # def summary(strand_complex_: Complex) -> str:
+    #     bps = _violation_base_pairs(strand_complex_)
+    #     if len(bps) == 0:
+    #         return "\tAll base pairs satisfy thresholds."
+    #     summary_list = []
+    #     for bp in bps:
+    #         i = bp.base_index1
+    #         j = bp.base_index2
+    #         p = bp.base_pairing_probability
+    #         t = bp.base_pair_type
+    #         summary_list.append(
+    #             f'\t{i},{j}: {math.floor(100 * p)}% (<{round(100 * base_type_probability_threshold[t])}% [{t}])')
+    #     return '\n'.join(summary_list)
 
     def _violation_base_pairs(strand_complex_: Complex) -> List[_BasePair]:
         nupack_complex_result = dv.nupack_complex_base_pair_probabilities(strand_complex_,
@@ -6768,7 +6551,8 @@ def _violation_base_pairs(strand_complex_: Complex) -> List[_BasePair]:
                 # Determine if base pair is adjacent to exterior base pair
                 prob_thres = internal_base_pair_prob
                 bp_type = BasePairType.INTERIOR_TO_STRAND
-                if i == 1 and d1_5p_d2_3p_ext_bp_type is not BasePairType.INTERIOR_TO_STRAND or i == domain_length_ - 2 and d1_3p_d2_5p_ext_bp_prob_thres is not BasePairType.INTERIOR_TO_STRAND:
+                if i == 1 and d1_5p_d2_3p_ext_bp_type is not BasePairType.INTERIOR_TO_STRAND \
+                        or i == domain_length_ - 2 and d1_3p_d2_5p_ext_bp_prob_thres is not BasePairType.INTERIOR_TO_STRAND:
                     prob_thres = border_internal_base_pair_prob
                     bp_type = BasePairType.ADJACENT_TO_EXTERIOR_BASE_PAIR
 
@@ -6788,7 +6572,6 @@ def _violation_base_pairs(strand_complex_: Complex) -> List[_BasePair]:
                              short_description=short_description,
                              weight=weight,
                              score_transfer_function=score_transfer_function,
-                             threaded=threaded,
+                             parallel=parallel,
                              complexes=tuple(strand_complexes),
-                             evaluate=evaluate,
-                             summary=summary)
+                             _evaluate=evaluate)
diff --git a/dsd/search.py b/dsd/search.py
index 9cfe1c1c..a444cf09 100644
--- a/dsd/search.py
+++ b/dsd/search.py
@@ -54,18 +54,16 @@
 import shutil
 import sys
 import logging
-import pprint
 from collections import defaultdict, deque
 import collections.abc as abc
 from dataclasses import dataclass, field
 from typing import List, Tuple, Sequence, Set, FrozenSet, Optional, Dict, Callable, Iterable, Any, \
-    Deque, TypeVar
+    Deque, TypeVar, cast
 import statistics
 import textwrap
 import time
 import re
 import datetime
-from pprint import pprint
 
 import numpy.random
 from ordered_set import OrderedSet
@@ -90,8 +88,10 @@
 
 from dsd.constraints import Domain, Strand, Design, Constraint, DomainConstraint, StrandConstraint, \
     DomainPairConstraint, StrandPairConstraint, ConstraintWithDomainPairs, ConstraintWithStrandPairs, \
-    logger, all_pairs, all_pairs_iterator, ConstraintWithDomains, ConstraintWithStrands, \
-    ComplexConstraint, ConstraintWithComplexes, Complex
+    logger, all_pairs, ConstraintWithDomains, ConstraintWithStrands, \
+    ComplexConstraint, ConstraintWithComplexes, Complex, DomainsConstraint, StrandsConstraint, \
+    DomainPairsConstraint, StrandPairsConstraint, ComplexesConstraint, DesignPart, DesignConstraint, \
+    DomainPair, StrandPair
 import dsd.constraints as dc
 
 from dsd.stopwatch import Stopwatch
@@ -158,169 +158,58 @@ def _violations_of_constraints(design: Design,
             assert not domain_changed.fixed
             violation_set.remove_violations_of_domain(domain_changed)
 
-    # individual domain constraints across all domains in Design
-    for domain_constraint in design.domain_constraints:
-        domains_to_check = _determine_domains_to_check(design.domains, domains_changed, domain_constraint)
+    constraints: List[Constraint] = cast(List[Constraint], design.domain_constraints) + \
+                                    design.strand_constraints + \
+                                    design.domain_pair_constraints + \
+                                    design.strand_pair_constraints + \
+                                    design.domains_constraints + \
+                                    design.strands_constraints + \
+                                    design.domain_pairs_constraints + \
+                                    design.strand_pairs_constraints + \
+                                    design.design_constraints + \
+                                    design.complex_constraints  # noqa
 
-        current_score_gap = violation_set_old.total_score() - violation_set.total_score() \
-            if never_increase_score and violation_set_old is not None else None
-
-        domain_violations, quit_early_in_func = _violations_of_domain_constraint(
-            domains=domains_to_check, constraint=domain_constraint, current_score_gap=current_score_gap)
-        violation_set.update(domain_violations)
-
-        quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
-        assert quit_early == quit_early_in_func
-        if quit_early:
-            return violation_set
-
-    # individual strand constraints across all strands in Design
-    for strand_constraint in design.strand_constraints:
-        strands_to_check = _determine_strands_to_check(design.strands, domains_changed, strand_constraint)
+    for constraint in constraints:
+        parts_to_check = find_parts_to_check(constraint, design, domains_changed)
 
         current_score_gap = violation_set_old.total_score() - violation_set.total_score() \
             if never_increase_score and violation_set_old is not None else None
 
-        strand_violations, quit_early_in_func = _violations_of_strand_constraint(
-            strands=strands_to_check, constraint=strand_constraint, current_score_gap=current_score_gap)
-        violation_set.update(strand_violations)
-
-        quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
-        assert quit_early == quit_early_in_func
-        if quit_early:
-            return violation_set
-
-    # pairs of domains
-    for domain_pair_constraint in design.domain_pair_constraints:
-        current_score_gap = violation_set_old.total_score() - violation_set.total_score() \
-            if never_increase_score and violation_set_old is not None else None
-        domain_pair_violations, quit_early_in_func = _violations_of_domain_pair_constraint(
-            domains=design.domains, constraint=domain_pair_constraint, domains_changed=domains_changed,
-            current_score_gap=current_score_gap)
-        violation_set.update(domain_pair_violations)
-
-        quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
-        assert quit_early == quit_early_in_func
-        if quit_early:
-            return violation_set
-
-    # pairs of strands
-    for strand_pair_constraint in design.strand_pair_constraints:
-        current_score_gap = violation_set_old.total_score() - violation_set.total_score() \
-            if never_increase_score and violation_set_old is not None else None
-        strand_pair_violations, quit_early_in_func = _violations_of_strand_pair_constraint(
-            strands=design.strands, constraint=strand_pair_constraint, domains_changed=domains_changed,
-            current_score_gap=current_score_gap)
-        violation_set.update(strand_pair_violations)
-
-        quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
-        assert quit_early == quit_early_in_func
-        if quit_early:
-            return violation_set
+        violations, quit_early_in_func = _violations_of_constraint(
+            parts=parts_to_check, constraint=constraint, current_score_gap=current_score_gap,
+            domains_changed=domains_changed, design=design)
+        violation_set.update(violations)
 
-    # complexes
-    for complex_constraint in design.complex_constraints:
-        current_score_gap = violation_set_old.total_score() - violation_set.total_score() \
-            if never_increase_score and violation_set_old is not None else None
-        complex_violations, quit_early_in_func = _violations_of_complex_constraint(
-            constraint=complex_constraint, domains_changed=domains_changed,
-            current_score_gap=current_score_gap)
-        violation_set.update(complex_violations)
+        # TODO: wasteful; get this info as a second return from the first call
+        parts_to_check_total = find_parts_to_check(constraint, design, None)
+        violation_set.num_checked[constraint] = len(parts_to_check_total)
 
         quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
         assert quit_early == quit_early_in_func
         if quit_early:
             return violation_set
 
-    # constraints that process each domain, but all at once (e.g., to hand off in batch to RNAduplex)
-    for domains_constraint in design.domains_constraints:
-        domains_to_check = _determine_domains_to_check(design.domains, domains_changed, domains_constraint)
-
-        if log_names_of_domains_and_strands_checked:
-            logger.debug(f'$ for domains constraint {domains_constraint.description}, '
-                         f'checking these domains:\n'
-                         f'${pprint.pformat(domains_to_check, indent=pprint_indent)}')
-
-        sets_of_violating_domains_weights = domains_constraint(domains_to_check)
-        domains_violations = _convert_sets_of_violating_domains_to_violations(
-            domains_constraint, sets_of_violating_domains_weights)
-        violation_set.update(domains_violations)
-
-        quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
-        if quit_early:
-            return violation_set
+    return violation_set
 
-    # constraints that process each strand, but all at once (e.g., to hand off in batch to RNAduplex)
-    for strands_constraint in design.strands_constraints:
-        strands_to_check = _determine_strands_to_check(design.strands, domains_changed, strands_constraint)
-
-        if log_names_of_domains_and_strands_checked:
-            logger.debug(f'$ for strands constraint {strands_constraint.description}, '
-                         f'checking these strands:\n'
-                         f'${pprint.pformat(strands_to_check, indent=pprint_indent)}')
-
-        if len(strands_to_check) > 0:
-            sets_of_violating_domains_weights = strands_constraint(strands_to_check)
-            domains_violations = _convert_sets_of_violating_domains_to_violations(
-                strands_constraint, sets_of_violating_domains_weights)
-            violation_set.update(domains_violations)
-
-            quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
-            if quit_early:
-                return violation_set
-
-    # constraints that process all pairs of domains at once (e.g., to hand off in batch to RNAduplex)
-    for domain_pairs_constraint in design.domain_pairs_constraints:
-        domain_pairs_to_check = _determine_domain_pairs_to_check(design.domains, domains_changed,
-                                                                 domain_pairs_constraint)
-
-        if log_names_of_domains_and_strands_checked:
-            logger.debug(f'$ for domain pairs constraint {domain_pairs_constraint.description}, '
-                         f'checking these strand pairs:\n'
-                         f'${pprint.pformat(domain_pairs_to_check, indent=pprint_indent)}')
-
-        if len(domain_pairs_to_check) > 0:
-            sets_of_violating_domains_weights = domain_pairs_constraint(domain_pairs_to_check)
-            domains_violations = _convert_sets_of_violating_domains_to_violations(
-                domain_pairs_constraint, sets_of_violating_domains_weights)
-            violation_set.update(domains_violations)
-
-            quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
-            if quit_early:
-                return violation_set
-
-    # constraints that process all pairs of strands at once (e.g., to hand off in batch to RNAduplex)
-    for strand_pairs_constraint in design.strand_pairs_constraints:
-        strand_pairs_to_check = _determine_strand_pairs_to_check(design.strands, domains_changed,
-                                                                 strand_pairs_constraint)
-        if log_names_of_domains_and_strands_checked:
-            logger.debug(f'$ for strand pairs constraint {strand_pairs_constraint.description}, '
-                         f'checking these strand pairs:\n'
-                         f'${pprint.pformat(strand_pairs_to_check, indent=pprint_indent)}')
-
-        if len(strand_pairs_to_check) > 0:
-            sets_of_violating_domains_weights = strand_pairs_constraint(strand_pairs_to_check)
-            domains_violations = _convert_sets_of_violating_domains_to_violations(
-                strand_pairs_constraint, sets_of_violating_domains_weights)
-            violation_set.update(domains_violations)
-
-            quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
-            if quit_early:
-                return violation_set
-
-    # constraints that processes whole design at once (for anything not captured by the above, e.g.,
-    # processing all triples of strands)
-    for design_constraint in design.design_constraints:
-        sets_of_violating_domains_weights = design_constraint(design, domains_changed)
-        domains_violations = _convert_sets_of_violating_domains_to_violations(
-            design_constraint, sets_of_violating_domains_weights)
-        violation_set.update(domains_violations)
 
-        quit_early = _quit_early(never_increase_score, violation_set, violation_set_old)
-        if quit_early:
-            return violation_set
-
-    return violation_set
+def find_parts_to_check(constraint: dc.Constraint, design: dc.Design,
+                        domains_changed: Optional[Iterable[Domain]]) -> Sequence[dc.DesignPart]:
+    parts_to_check: Sequence[dc.DesignPart]
+    if isinstance(constraint, ConstraintWithDomains):
+        parts_to_check = _determine_domains_to_check(design.domains, domains_changed, constraint)
+    elif isinstance(constraint, ConstraintWithStrands):
+        parts_to_check = _determine_strands_to_check(design.strands, domains_changed, constraint)
+    elif isinstance(constraint, ConstraintWithDomainPairs):
+        parts_to_check = _determine_domain_pairs_to_check(design.domains, domains_changed, constraint)
+    elif isinstance(constraint, ConstraintWithStrandPairs):
+        parts_to_check = _determine_strand_pairs_to_check(design.strands, domains_changed, constraint)
+    elif isinstance(constraint, ConstraintWithComplexes):
+        parts_to_check = _determine_complexes_to_check(domains_changed, constraint)
+    elif isinstance(constraint, dc.DesignConstraint):
+        parts_to_check = []  # not used when checking DesignConstraint
+    else:
+        raise NotImplementedError()
+    return parts_to_check
 
 
 def _is_significantly_greater(x: float, y: float) -> bool:
@@ -344,29 +233,6 @@ def _at_least_one_domain_unfixed(pair: Tuple[Domain, Domain]) -> bool:
     return not (pair[0].fixed and pair[1].fixed)
 
 
-def _determine_domain_pairs_to_check(all_domains: Iterable[Domain],
-                                     domains_changed: Optional[Iterable[Domain]],
-                                     constraint: ConstraintWithDomainPairs) \
-        -> Sequence[Tuple[Domain, Domain]]:
-    """
-    Determines domain pairs to check between domains in `all_domains`.
-    If `domain_changed` is None, then this is all pairs where they are not both fixed if constraint.pairs
-    is None, otherwise it is constraint.pairs.
-    If `domain_changed` is not None, then among those pairs specified above,
-    it is all pairs where one of the two is `domain_changed`.
-    """
-    # either all pairs, or just constraint.pairs if specified
-    domain_pairs_to_check_if_domain_changed_none = constraint.pairs if constraint.pairs is not None \
-        else all_pairs_iterator(all_domains, with_replacement=True, where=_at_least_one_domain_unfixed)
-
-    # filter out those not containing domain_change if specified
-    domain_pairs_to_check = list(domain_pairs_to_check_if_domain_changed_none) if domains_changed is None \
-        else [(domain1, domain2) for domain1, domain2 in domain_pairs_to_check_if_domain_changed_none
-              if domain1 in domains_changed or domain2 in domains_changed]
-
-    return domain_pairs_to_check
-
-
 def _determine_domains_to_check(all_domains: Iterable[Domain],
                                 domains_changed: Optional[Iterable[Domain]],
                                 constraint: ConstraintWithDomains) -> Sequence[Domain]:
@@ -389,6 +255,58 @@ def _determine_domains_to_check(all_domains: Iterable[Domain],
     return domains_to_check
 
 
+def _determine_strands_to_check(all_strands: Iterable[Strand],
+                                domains_changed: Optional[Iterable[Domain]],
+                                constraint: ConstraintWithStrands) -> Sequence[Strand]:
+    """
+    Similar to _determine_domains_to_check but for strands.
+    """
+    # either all pairs, or just constraint.pairs if specified
+    strands_to_check_if_domain_changed_none = constraint.strands if constraint.strands is not None \
+        else [strand for strand in all_strands if not strand.fixed]
+
+    # filter out those not containing domain_change if specified
+    strands_to_check: List[Strand] = []
+    if domains_changed is None:
+        strands_to_check = strands_to_check_if_domain_changed_none
+    else:
+        for strand in strands_to_check_if_domain_changed_none:
+            for domain_changed in domains_changed:
+                if domain_changed in strand.domains:
+                    strands_to_check.append(strand)
+                    break
+
+    return strands_to_check
+
+
+def _determine_domain_pairs_to_check(all_domains: Iterable[Domain],
+                                     domains_changed: Optional[Iterable[Domain]],
+                                     constraint: ConstraintWithDomainPairs) \
+        -> Sequence[DomainPair]:
+    """
+    Determines domain pairs to check between domains in `all_domains`.
+    If `domain_changed` is None, then this is all pairs where they are not both fixed if constraint.pairs
+    is None, otherwise it is constraint.pairs.
+    If `domain_changed` is not None, then among those pairs specified above,
+    it is all pairs where one of the two is `domain_changed`.
+    """
+    # either all pairs, or just constraint.pairs if specified
+    if constraint.pairs is not None:
+        domain_pairs_to_check_if_domain_changed_none: List[DomainPair] = \
+            [DomainPair(pair[0], pair[1]) for pair in constraint.pairs]
+    else:
+        pairs = all_pairs(all_domains, with_replacement=True, where=_at_least_one_domain_unfixed)
+        domain_pairs_to_check_if_domain_changed_none = [DomainPair(pair[0], pair[1]) for pair in pairs]
+
+    # filter out those not containing domain_change if specified
+    domain_pairs_to_check = list(domain_pairs_to_check_if_domain_changed_none) if domains_changed is None \
+        else [domain_pair for domain_pair in
+              domain_pairs_to_check_if_domain_changed_none
+              if domain_pair.domain1 in domains_changed or domain_pair.domain2 in domains_changed]
+
+    return domain_pairs_to_check
+
+
 def _at_least_one_strand_unfixed(pair: Tuple[Strand, Strand]) -> bool:
     return not (pair[0].fixed and pair[1].fixed)
 
@@ -396,31 +314,35 @@ def _at_least_one_strand_unfixed(pair: Tuple[Strand, Strand]) -> bool:
 def _determine_strand_pairs_to_check(all_strands: Iterable[Strand],
                                      domains_changed: Optional[Iterable[Domain]],
                                      constraint: ConstraintWithStrandPairs) -> \
-        Sequence[Tuple[Strand, Strand]]:
+        Sequence[StrandPair]:
     """
     Similar to _determine_domain_pairs_to_check but for strands.
     """
     # either all pairs, or just constraint.pairs if specified
-    strand_pairs_to_check_if_domain_changed_none = constraint.pairs if constraint.pairs is not None \
-        else all_pairs(all_strands, where=_at_least_one_strand_unfixed)
+    if constraint.pairs is not None:
+        strand_pairs_to_check_if_domain_changed_none: List[StrandPair] = \
+            [StrandPair(pair[0], pair[1]) for pair in constraint.pairs]
+    else:
+        pairs = all_pairs(all_strands, where=_at_least_one_strand_unfixed)
+        strand_pairs_to_check_if_domain_changed_none = [StrandPair(pair[0], pair[1]) for pair in pairs]
 
     # filter out those not containing domain_change if specified
-    strand_pairs_to_check: List[Tuple[Strand, Strand]] = []
+    strand_pairs_to_check: List[StrandPair] = []
     if domains_changed is None:
         strand_pairs_to_check = strand_pairs_to_check_if_domain_changed_none
     else:
-        for strand1, strand2 in strand_pairs_to_check_if_domain_changed_none:
+        for strand_pair in strand_pairs_to_check_if_domain_changed_none:
             for domain_changed in domains_changed:
-                if domain_changed in strand1.domains or domain_changed in strand2.domains:
-                    strand_pairs_to_check.append((strand1, strand2))
+                if domain_changed in strand_pair.strand1.domains or \
+                        domain_changed in strand_pair.strand2.domains:
+                    strand_pairs_to_check.append(strand_pair)
                     break
 
     return strand_pairs_to_check
 
 
 def _determine_complexes_to_check(domains_changed: Optional[Iterable[Domain]],
-                                  constraint: ConstraintWithComplexes) -> \
-        Tuple[Complex]:
+                                  constraint: ConstraintWithComplexes) -> Tuple[Complex]:
     """
     Similar to _determine_domain_pairs_to_check but for complexes.
     """
@@ -445,30 +367,6 @@ def _determine_complexes_to_check(domains_changed: Optional[Iterable[Domain]],
         return tuple(complexes_to_check)
 
 
-def _determine_strands_to_check(all_strands: Iterable[Strand],
-                                domains_changed: Optional[Iterable[Domain]],
-                                constraint: ConstraintWithStrands) -> Sequence[Strand]:
-    """
-    Similar to _determine_domain_pairs_to_check but for strands.
-    """
-    # either all pairs, or just constraint.pairs if specified
-    strands_to_check_if_domain_changed_none = constraint.strands if constraint.strands is not None \
-        else [strand for strand in all_strands if not strand.fixed]
-
-    # filter out those not containing domain_change if specified
-    strands_to_check: List[Strand] = []
-    if domains_changed is None:
-        strands_to_check = strands_to_check_if_domain_changed_none
-    else:
-        for strand in strands_to_check_if_domain_changed_none:
-            for domain_changed in domains_changed:
-                if domain_changed in strand.domains:
-                    strands_to_check.append(strand)
-                    break
-
-    return strands_to_check
-
-
 def _strands_containing_domains(domains: Optional[Iterable[Domain]], strands: List[Strand]) -> List[Strand]:
     """
     :param domains:
@@ -488,435 +386,112 @@ def _strands_containing_domains(domains: Optional[Iterable[Domain]], strands: Li
         return list(strands_set)
 
 
-def _convert_sets_of_violating_domains_to_violations(
-        constraint: Constraint, sets_of_violating_domains: Iterable[Tuple[OrderedSet[Domain], float]]) \
-        -> Dict[Domain, OrderedSet[dc.Violation]]:
-    domains_violations: Dict[Domain, OrderedSet[dc.Violation]] = defaultdict(OrderedSet)
-    for domain_set, score in sets_of_violating_domains:
-        violation = dc.Violation(constraint, domain_set, score)
-        for domain in domain_set:
-            domain_violations = domains_violations[domain]
-            domain_violations.add(violation)
-    return domains_violations
-
-
 _empty_frozen_set: FrozenSet = frozenset()
 
 
-# XXX: Although this is written very generally for multiple domains; for most iterationg only one domain
-# changes and we are only checking that one domain in domains. So there's no optimization here for
-# quitting early since we are usually only checking a single constraint.
-def _violations_of_domain_constraint(domains: Iterable[Domain],
-                                     constraint: DomainConstraint,
-                                     current_score_gap: Optional[float],
-                                     ) -> Tuple[Dict[Domain, OrderedSet[dc.Violation]], bool]:
+def _violations_of_constraint(parts: Sequence[DesignPart],
+                              constraint: Constraint[DesignPart],
+                              current_score_gap: Optional[float],
+                              domains_changed: Optional[Iterable[Domain]] = None,
+                              design: Optional[Design] = None,  # only used with DesignConstraint
+                              ) -> Tuple[Dict[Domain, OrderedSet[dc.Violation]], bool]:
     violations: Dict[Domain, OrderedSet[dc.Violation]] = defaultdict(OrderedSet)
-    unfixed_domains = [domain for domain in domains if not domain.fixed]
-    violating_domains_scores: List[Optional[Tuple[Domain, float]]] = []
+    violating_parts_scores_summaries: List[Tuple[DesignPart, float, str]] = []
 
     score_discovered_here: float = 0.0
     quit_early = False
     num_threads = dc.cpu_count()
-    chunk_size = num_threads
-
-    if log_names_of_domains_and_strands_checked:
-        logger.debug(f'$ for domain constraint {constraint.description}, checking these domains:')
-        logger.debug(f'{pprint.pformat(unfixed_domains, indent=pprint_indent)}')
-
-    if (not constraint.threaded
-            or num_threads == 1
-            or len(unfixed_domains) == 1):
-        logger.debug(f'NOT using threading for domain constraint {constraint.description}')
-        for domain in unfixed_domains:
-            score: float = constraint(domain.sequence, domain)
-            if score > 0.0:
-                violating_domains_scores.append((domain, score))
-                if current_score_gap is not None:
-                    score_discovered_here += constraint.weight * score
-                    if _is_significantly_greater(score_discovered_here, current_score_gap):
-                        quit_early = True
-                        break
-    else:
-        logger.debug(f'using threading for domain constraint {constraint.description}')
-
-        domains_to_check = unfixed_domains
 
-        def sequence_to_score(sequence: str) -> float:
-            return constraint(sequence, None)
-
-        if current_score_gap is None:
-            sequences_to_check = (domain.sequence for domain in domains_to_check)
-            scores = list(_process_pool.map(sequence_to_score, sequences_to_check))
-            violating_domains_scores = [(domain, score) for domain, score in zip(domains_to_check, scores)
-                                        if score > 0]
-        else:
-            chunks = dc.chunker(domains_to_check, chunk_size)
-            for domain_chunk in chunks:
-                sequence_chunk = [domain.sequence for domain in domain_chunk]
-                scores = list(_process_pool.map(sequence_to_score, sequence_chunk))
-                violating_domains_chunk = [(strand, score) for strand, score in zip(domain_chunk, scores)
-                                           if score > 0]
-                violating_domains_scores.extend(violating_domains_chunk)
-
-                # quit early if possible
-                total_score_chunk = sum(score for _, score in violating_domains_chunk)
-                score_discovered_here += constraint.weight * total_score_chunk
+    if not constraint.parallel or num_threads == 1 or len(parts) == 1:
+        if isinstance(constraint,
+                      (DomainConstraint, StrandConstraint,
+                       DomainPairConstraint, StrandPairConstraint, ComplexConstraint)):
+            for part in parts:
+                score, summary = _evaluate_individual_part_constraint(constraint, part)
+                if score > 0.0:
+                    violating_parts_scores_summaries.append((part, score, summary))
+                    if current_score_gap is not None:
+                        score_discovered_here += score
+                        if _is_significantly_greater(score_discovered_here, current_score_gap):
+                            quit_early = True
+                            break
+
+        elif isinstance(constraint,
+                        (DomainsConstraint, StrandsConstraint,
+                         DomainPairsConstraint, StrandPairsConstraint,
+                         ComplexesConstraint, DesignConstraint)):
+            if isinstance(constraint, DesignConstraint):
+                violating_parts_scores_summaries = constraint.evaluate_design(design, domains_changed)
+            else:
+                # XXX: I don't understand the mypy error on the next line
+                violating_parts_scores_summaries = constraint.evaluate_bulk(parts)
+
+            # we can't quite this function early,
+            # but we can let the caller know to stop evaluating constraints
+            total_score = sum(score for _, score, _ in violating_parts_scores_summaries)
+            if current_score_gap is not None:
+                score_discovered_here += total_score
                 if _is_significantly_greater(score_discovered_here, current_score_gap):
                     quit_early = True
-                    break
-
-    for violating_domain_score in violating_domains_scores:
-        if violating_domain_score is not None:
-            violating_domain, score = violating_domain_score
-            violation = dc.Violation(constraint, [violating_domain], score)
-            violations[violating_domain].add(violation)
-
-    return violations, quit_early
-
-
-def _violations_of_strand_constraint(strands: Iterable[Strand],
-                                     constraint: StrandConstraint,
-                                     current_score_gap: Optional[float],
-                                     ) -> Tuple[Dict[Domain, OrderedSet[dc.Violation]], bool]:
-    """
-    :param strands:
-        Strands to check for violations
-    :param constraint:
-        Constraint to check.
-    :param current_score_gap:
-        Current gap between total score of constraint violations found so far and total score of
-        optimal design. Used for quitting early.
-    :return:
-        1. dict mapping each domain to the set of violations that blame it
-        2. bool indicating whether we quit early
-    """
-    strands_to_check = [strand for strand in strands if not strand.fixed]
-
-    violating_strands_scores: List[Tuple[Strand, float]] = []
+        else:
+            raise AssertionError(f'constraint {constraint} of unrecognized type {type(constraint)}')
 
-    score_discovered_here: float = 0.0
-    quit_early = False
-    num_threads = dc.cpu_count()
-    chunk_size = num_threads
-
-    if log_names_of_domains_and_strands_checked:
-        logger.debug(f'$ for strand constraint {constraint.description}, checking these strands:')
-        logger.debug(f'$ {pprint.pformat(strands_to_check, indent=pprint_indent)}')
-
-    if (not constraint.threaded
-            or num_threads == 1
-            or len(strands_to_check) == 1
-            or (current_score_gap is not None and chunk_size == 1)):
-        logger.debug(f'NOT using threading for strand constraint {constraint.description}')
-        for strand in strands_to_check:
-            score: float = constraint(strand.sequence(), strand)
-            if score > 0.0:
-                violating_strands_scores.append((strand, score))
-                if current_score_gap is not None:
-                    score_discovered_here += constraint.weight * score
-                    if _is_significantly_greater(score_discovered_here, current_score_gap):
-                        quit_early = True
-                        break
     else:
-        logger.debug(f'using threading for strand constraint {constraint.description}')
-        assert constraint.sequence_only  # should have been checked in constraint post_init
-
-        def sequence_to_score(sequence: str) -> float:
-            return constraint(sequence, None)
-
-        if current_score_gap is None:
-            sequences_to_check = (strand.sequence() for strand in strands_to_check)
-            scores = list(_process_pool.map(sequence_to_score, sequences_to_check))
-            violating_strands_scores = [(strand, score) for strand, score in zip(strands_to_check, scores)
-                                        if score > 0]
-        else:
-            chunks = dc.chunker(strands_to_check, chunk_size)
-            for strand_chunk in chunks:
-                sequence_chunk = [strand.sequence() for strand in strand_chunk]
-                scores = list(_process_pool.map(sequence_to_score, sequence_chunk))
-                violating_strands_chunk = [(strand, score) for strand, score in zip(strand_chunk, scores)
-                                           if score > 0]
-                violating_strands_scores.extend(violating_strands_chunk)
-
-                # quit early if possible
-                total_score_chunk = sum(score for _, score in violating_strands_chunk)
-                score_discovered_here += constraint.weight * total_score_chunk
-                if _is_significantly_greater(score_discovered_here, current_score_gap):
-                    quit_early = True
-                    break
+        raise NotImplementedError('TODO: implement parallelization')
 
-    violations: Dict[Domain, OrderedSet[dc.Violation]] = defaultdict(OrderedSet)
-    for strand, score in violating_strands_scores:
-        unfixed_domains_set = OrderedSet(strand.unfixed_domains())
-        violation = dc.Violation(constraint, unfixed_domains_set, score)
-        for domain in unfixed_domains_set:
+    for part, score, summary in violating_parts_scores_summaries:
+        domains = _domains_in_part(part, exclude_fixed=True)
+        violation = dc.Violation(constraint=constraint, part=part, domains=domains,
+                                 score=score, summary=summary)
+        for domain in domains:
             violations[domain].add(violation)
 
     return violations, quit_early
 
 
-T = TypeVar('T')
-
-
-def remove_none_from_list(lst: Iterable[Optional[T]]) -> List[T]:
-    return [elt for elt in lst if elt is not None]
-
-
-def _violations_of_domain_pair_constraint(domains: Iterable[Domain],
-                                          constraint: DomainPairConstraint,
-                                          domains_changed: Optional[Iterable[Domain]],
-                                          current_score_gap: Optional[float],
-                                          ) -> Tuple[Dict[Domain, OrderedSet[dc.Violation]], bool]:
-    # If specified, current_score_gap is the current difference between the score of violated constraints
-    # that have been found so far in the current iteration, compared to the total score of violated
-    # constraints in the optimal solution so far. It is positive
-    # (i.e., total_score_opt - total_score_cur_so_far)
-    # If specified and it is discovered while looping in this function that total_score_cur_so_far plus
-    # the score of violated constraints discovered in this function exceeds total_score_opt, quit early.
-    domain_pairs_to_check: Sequence[Tuple[Domain, Domain]] = \
-        _determine_domain_pairs_to_check(domains, domains_changed, constraint)
-
-    if log_names_of_domains_and_strands_checked:
-        logger.debug(f'$ for domain pair constraint {constraint.description}, checking these domain pairs:')
-        logger.debug(f'$ {pprint.pformat(domain_pairs_to_check, indent=pprint_indent)}')
-
-    violating_domain_pairs_scores: List[Optional[Tuple[Domain, Domain, float]]] = []
-
-    score_discovered_here: float = 0.0
-    quit_early = False
-
-    cpu_count = dc.cpu_count()
-
-    # since each domain pair check is already parallelized for the four domain pairs
-    # (d1,d2), (d1,w2), (w1,d2), (w1,w2), we take smaller chunks
-    chunk_size = cpu_count // 4
-
-    if (not constraint.threaded
-            or cpu_count == 1
-            or (current_score_gap is not None and chunk_size == 1)):
-        logger.debug(f'NOT using threading for domain pair constraint {constraint.description}')
-        for domain1, domain2 in domain_pairs_to_check:
-            assert not domain1.fixed or not domain2.fixed
-            score: float = constraint(domain1.sequence, domain2.sequence, domain1, domain2)
-            if score > 0.0:
-                violating_domain_pairs_scores.append((domain1, domain2, score))
-                if current_score_gap is not None:
-                    score_discovered_here += constraint.weight * score
-                    if _is_significantly_greater(score_discovered_here, current_score_gap):
-                        quit_early = True
-                        break
+def _evaluate_individual_part_constraint(constraint: dc.Constraint, part: dc.DesignPart) -> Tuple[float, str]:
+    if isinstance(constraint, (DomainConstraint, StrandConstraint)):
+        assert isinstance(part, (Domain, Strand))
+        score, summary = constraint.evaluate((part.sequence(),), part)
+    elif isinstance(constraint, (DomainPairConstraint, StrandPairConstraint, ComplexConstraint)):
+        assert isinstance(part, (DomainPair, StrandPair, Complex))
+        seqs = tuple(indv_part.sequence() for indv_part in part.individual_parts())
+        # XXX: not sure why mypy error on the next line
+        score, summary = constraint.evaluate(seqs, part)
     else:
-        logger.debug(f'using threading for domain pair constraint {constraint.description}')
-
-        def sequence_pair_to_score(seq_pair: Tuple[str, str]) -> float:
-            seq1, seq2 = seq_pair
-            return constraint(seq1, seq2, None, None)
-
-        if current_score_gap is None:
-            sequence_pairs_to_check = [(domain1.sequence, domain2.sequence)
-                                       for domain1, domain2 in domain_pairs_to_check]
-            scores = list(_process_pool.map(sequence_pair_to_score, sequence_pairs_to_check))
-            violating_domain_pairs_scores = [(domain1, domain2, score) for (domain1, domain2), score in
-                                             zip(domain_pairs_to_check, scores) if score > 0]
-
-        else:
-            chunks = dc.chunker(domain_pairs_to_check, chunk_size)
-            for domain_pair_chunk in chunks:
-                sequence_chunk = [(domain1.sequence, domain2.sequence)
-                                  for domain1, domain2 in domain_pair_chunk]
-                scores = list(_process_pool.map(sequence_pair_to_score, sequence_chunk))
-                violating_domain_pairs_chunk = [(domain1, domain2, score) for (domain1, domain2), score in
-                                                zip(domain_pair_chunk, scores) if score > 0]
-                violating_domain_pairs_scores.extend(violating_domain_pairs_chunk)
-
-                # quit early if possible
-                total_score_chunk = sum(score for _, _, score in violating_domain_pairs_chunk)
-                score_discovered_here += constraint.weight * total_score_chunk
-                if _is_significantly_greater(score_discovered_here, current_score_gap):
-                    quit_early = True
-                    break
-
-    violations: Dict[Domain, OrderedSet[dc.Violation]] = defaultdict(OrderedSet)
-    violating_domain_pair_score: Optional[Tuple[Domain, Domain, float]]
-    for violating_domain_pair_score in violating_domain_pairs_scores:
-        if violating_domain_pair_score is not None:
-            domain1, domain2, score = violating_domain_pair_score
-            unfixed_domains_set: Set[Domain] = set()
-            if not domain1.fixed:
-                unfixed_domains_set.add(domain1)
-            if not domain2.fixed:
-                unfixed_domains_set.add(domain2)
-            violation = dc.Violation(constraint, frozenset(unfixed_domains_set), score)
-            if not domain1.fixed:
-                violations[domain1].add(violation)
-            if not domain2.fixed:
-                violations[domain2].add(violation)
-
-    return violations, quit_early
-
+        raise AssertionError()
+    return score, summary
 
-def _violations_of_strand_pair_constraint(strands: Iterable[Strand],
-                                          constraint: StrandPairConstraint,
-                                          domains_changed: Optional[Iterable[Domain]],
-                                          current_score_gap: Optional[float],
-                                          ) -> Tuple[Dict[Domain, OrderedSet[dc.Violation]], bool]:
-    strand_pairs_to_check: Sequence[Tuple[Strand, Strand]] = \
-        _determine_strand_pairs_to_check(strands, domains_changed, constraint)
 
-    if log_names_of_domains_and_strands_checked:
-        logger.debug(f'$ for strand pair constraint {constraint.description}, checking these strand pairs:')
-        logger.debug(f'$ {pprint.pformat(strand_pairs_to_check, indent=pprint_indent)}')
-
-    violating_strand_pairs_scores: List[Tuple[Strand, Strand, float]] = []
-
-    score_discovered_here: float = 0.0
-    quit_early = False
-    cpu_count = dc.cpu_count()
-    chunk_size = cpu_count
-
-    if (not constraint.threaded
-            or cpu_count == 1
-            or (current_score_gap is not None and chunk_size == 1)):
-        logger.debug(f'NOT using threading for strand pair constraint {constraint.description}')
-        for strand1, strand2 in strand_pairs_to_check:
-            assert not strand1.fixed or not strand2.fixed
-            if constraint.sequence_only:
-                score = constraint(strand1.sequence(), strand2.sequence(), None, None)
-            else:
-                score = constraint(strand1.sequence(), strand2.sequence(), strand1, strand2)
-            if score > 0.0:
-                violating_strand_pairs_scores.append((strand1, strand2, score))
-                if current_score_gap is not None:
-                    score_discovered_here += constraint.weight * score
-                    if _is_significantly_greater(score_discovered_here, current_score_gap):
-                        quit_early = True
-                        break
+def _domains_in_part(part: dc.DesignPart, exclude_fixed: bool) -> List[Domain]:
+    """
+    :param part:
+        :any:`DesignPart` (e.g., :any:`Strand`, :any:`Domani`, Tuple[:any:`Strand`, :any:`Strand`])
+    :param exclude_fixed:
+        whether to exclude :any:`Domain`'s with :data:`Domain.fixed` == True
+    :return:
+        domains associated with part (e.g., all domains in :any:`Strand`)
+    """
+    if isinstance(part, Domain):
+        return [part] if not (exclude_fixed and part.fixed) else []
+    elif isinstance(part, Strand):
+        return part.domains if not exclude_fixed else part.unfixed_domains()
+    elif isinstance(part, DomainPair):
+        return list(domain for domain in part.individual_parts() if not (exclude_fixed and domain.fixed))
+    elif isinstance(part, (StrandPair, Complex)):
+        domains_per_strand = [strand.domains if exclude_fixed else strand.unfixed_domains()
+                              for strand in part.individual_parts()]
+        domain_iterable: Iterable[Domain] = _flatten(domains_per_strand)
+        return list(domain_iterable)
     else:
-        logger.debug(f'using threading for strand pair constraint {constraint.description}')
-
-        assert constraint.sequence_only  # should have been checked in StrandPairConstraint post_init
-
-        def sequence_pair_to_score(seq_pair: Tuple[str, str]) -> float:
-            seq1, seq2 = seq_pair
-            return constraint(seq1, seq2, None, None)
-
-        if current_score_gap is None:
-            sequence_pairs_to_check = [(strand1.sequence(), strand2.sequence())
-                                       for strand1, strand2 in strand_pairs_to_check]
-            scores = list(_process_pool.map(sequence_pair_to_score, sequence_pairs_to_check))
-            violating_strand_pairs_scores = [(strand1, strand2, score) for (strand1, strand2), score in
-                                             zip(strand_pairs_to_check, scores) if score > 0]
-
-        else:
-            chunks = dc.chunker(strand_pairs_to_check, chunk_size)
-            for strand_pair_chunk in chunks:
-                sequence_chunk = [(strand1.sequence(), strand2.sequence())
-                                  for strand1, strand2 in strand_pair_chunk]
-                scores = list(_process_pool.map(sequence_pair_to_score, sequence_chunk))
-                violating_strand_pairs_chunk = [(strand1, strand2, score) for (strand1, strand2), score in
-                                                zip(strand_pair_chunk, scores) if score > 0]
-                violating_strand_pairs_scores.extend(violating_strand_pairs_chunk)
-
-                # quit early if possible
-                total_score_chunk = sum(score for _, _, score in violating_strand_pairs_chunk)
-                score_discovered_here += constraint.weight * total_score_chunk
-                if _is_significantly_greater(score_discovered_here, current_score_gap):
-                    quit_early = True
-                    break
+        raise NotImplementedError()
 
-    violations: Dict[Domain, OrderedSet[dc.Violation]] = defaultdict(OrderedSet)
-    for strand1, strand2, score in violating_strand_pairs_scores:
-        unfixed_domains_set = OrderedSet(strand1.unfixed_domains() + strand2.unfixed_domains())
-        violation = dc.Violation(constraint, unfixed_domains_set, score)
-        for domain in unfixed_domains_set:
-            violations[domain].add(violation)
-
-    return violations, quit_early
 
+T = TypeVar('T')
 
-def _violations_of_complex_constraint(constraint: ComplexConstraint,
-                                      domains_changed: Optional[Iterable[Domain]],
-                                      current_score_gap: Optional[float],
-                                      ) -> Tuple[Dict[Domain, OrderedSet[dc.Violation]], bool]:
-    complexes_to_check: Tuple[Complex] = \
-        _determine_complexes_to_check(domains_changed, constraint)
-
-    if log_names_of_domains_and_strands_checked:
-        logger.debug(f'$ for complex constraint {constraint.description}, checking these complexes:')
-        logger.debug(f'$ {pprint.pformat(complexes_to_check, indent=pprint_indent)}')
-
-    violating_complexes_scores: List[Optional[Tuple[Complex, float]]] = []
-
-    score_discovered_here: float = 0.0
-    quit_early = False
-    cpu_count = dc.cpu_count()
-    chunk_size = cpu_count
-
-    if (not constraint.threaded
-            or cpu_count == 1
-            or (current_score_gap is not None and chunk_size == 1)):
-        logger.debug(f'NOT using threading for strand pair constraint {constraint.description}')
-        for strand_complex in complexes_to_check:
-            score = constraint(strand_complex)
-            if score > 0.0:
-                violating_complexes_scores.append((strand_complex, score))
-                if current_score_gap is not None:
-                    score_discovered_here += constraint.weight * score
-                    if _is_significantly_greater(score_discovered_here, current_score_gap):
-                        quit_early = True
-                        break
-            if quit_early:
-                # Need to break out of checking each strand in complex since we added complex already
-                break
-    else:
-        logger.debug(f'NOT using threading for strand pair constraint {constraint.description}')
-
-        def complex_score_if_violates(strand_complex_: Complex) \
-                -> Optional[Tuple[Complex, float]]:
-            # return strand pair if it violates the constraint, else None
-            score_ = constraint(strand_complex_)
-            if score_ > 0.0:
-                return strand_complex_, score_
-            else:
-                return None
-
-        if current_score_gap is None:
-            violating_complexes_scores = list(
-                _process_pool.map(complex_score_if_violates, complexes_to_check))
-        else:
-            chunks = dc.chunker(complexes_to_check, chunk_size)
-            for complex_chunk in chunks:
-                violating_complexes_chunk_with_none: List[Optional[Tuple[Complex, float]]] = \
-                    _process_pool.map(complex_score_if_violates, complex_chunk)
-                violating_complexes_chunk: List[Tuple[Complex, float]] = \
-                    remove_none_from_list(violating_complexes_chunk_with_none)
-                violating_complexes_scores.extend(violating_complexes_chunk)
-
-                # quit early if possible
-                total_score_chunk = sum(
-                    complex_pair_score[1]
-                    for complex_pair_score in violating_complexes_chunk
-                    if complex_pair_score is not None)
-                score_discovered_here += constraint.weight * total_score_chunk
-                if _is_significantly_greater(score_discovered_here, current_score_gap):
-                    quit_early = True
-                    break
-
-    violations: Dict[Domain, OrderedSet[dc.Violation]] = defaultdict(OrderedSet)
-    violating_complex_score: Optional[Tuple[Complex, float]]
-    for violating_complex_score in violating_complexes_scores:
-        if violating_complex_score is not None:
-            strand_complex, score = violating_complex_score
-            unfixed_domains_set_builder = set()
-            strand: Strand
-            for strand in strand_complex:
-                unfixed_domains_set_builder.update(strand.unfixed_domains())
-            unfixed_domains_set = frozenset(unfixed_domains_set_builder)
-            violation = dc.Violation(constraint, unfixed_domains_set, score)
-            for domain in unfixed_domains_set:
-                violations[domain].add(violation)
 
-    return violations, quit_early
+def remove_none_from_list(lst: Iterable[Optional[T]]) -> List[T]:
+    return [elt for elt in lst if elt is not None]
 
 
 def _sequences_fragile_format_output_to_file(design: Design,
@@ -1157,10 +732,18 @@ class SearchParameters:
     whenever a new optimal sequence assignment is found.
     """
 
-    report_delay: float = 60.0
+    report_delay: float = 0.0
     """
     Every time the design improves, a report on the constraints is written, as long as it has been as
-    `report_delay` seconds since the last report was written. Since writing a report requires evaluating
+    `report_delay` seconds since the last report was written. 
+    
+    UPDATE: the following is no longer true, since dsd has been re-written to not re-do the work of
+    re-evaluating constraints when writing reports, so it is now quite fast to write reports.
+    The default has been changed to 0.0, i.e., reports will be written on every design update,
+    but the below explaination is left here to understand why this was ever an option. It may be
+    removed in the future.
+    
+    Since writing a report requires evaluating
     all constraints, it requires more time than a single iteration, which requires evaluating only those
     constraints involving the :any:`constraints.Domain` whose DNA sequence was changed.
     Thus the default value of 60 seconds avoids spending too much time writing reports, since the
@@ -1358,7 +941,7 @@ def search_for_dna_sequences(design: dc.Design, params: SearchParameters) -> Non
         iteration = 0
         time_of_last_improvement: float = -1.0
 
-        while len(violation_set_opt.non_fixed_violations) > 0 and \
+        while len(violation_set_opt.violations_nonfixed) > 0 and \
                 (params.max_iterations is None or iteration < params.max_iterations):
             _check_cpu_count(cpu_count)
 
@@ -1380,7 +963,7 @@ def search_for_dna_sequences(design: dc.Design, params: SearchParameters) -> Non
             # based on total score of new constraint violations compared to optimal assignment so far,
             # decide whether to keep the change
             # score_delta = violation_set_new.total_score() - violation_set_opt.total_score()
-            score_delta = violation_set_new.total_nonfixed_score() - violation_set_opt.total_nonfixed_score()
+            score_delta = violation_set_new.total_score_nonfixed() - violation_set_opt.total_score_nonfixed()
             prob_keep_change = params.probability_of_keeping_change(score_delta)
             keep_change = rng.random() < prob_keep_change if prob_keep_change < 1 else True
 
@@ -1399,7 +982,7 @@ def search_for_dna_sequences(design: dc.Design, params: SearchParameters) -> Non
                     write_report = False
                     # don't write report unless it is
                     if (time_of_last_improvement < 0  # first iteration
-                            or len(violation_set_opt.all_violations) == 0  # last iteration (search is over)
+                            or len(violation_set_opt.violations_all) == 0  # last iteration (search is over)
                             or current_time - time_of_last_improvement >= params.report_delay):  # > report_delay seconds since last report
                         time_of_last_improvement = current_time
                         write_report = True
@@ -1478,20 +1061,21 @@ def _reassign_domains(domains_opt: List[Domain], scores_opt: List[float], max_do
     for domain in independent_domains:
         # set sequence of domain_changed to random new sequence from its DomainPool
         assert domain not in original_sequences
-        previous_sequence = domain.sequence
+        previous_sequence = domain.sequence()
         original_sequences[domain] = previous_sequence
-        domain.sequence = domain.pool.generate_sequence(rng, previous_sequence)
+        new_sequence = domain.pool.generate_sequence(rng, previous_sequence)
+        domain.set_sequence(new_sequence)
 
     dependent_domains = [domain for domain in domains_changed if domain.dependent]
     for domain in dependent_domains:
-        original_sequences[domain] = domain.sequence
+        original_sequences[domain] = domain.sequence()
 
     return domains_changed, original_sequences
 
 
 def _unassign_domains(domains_changed: Iterable[Domain], original_sequences: Dict[Domain, str]) -> None:
     for domain_changed in domains_changed:
-        domain_changed.sequence = original_sequences[domain_changed]
+        domain_changed.set_sequence(original_sequences[domain_changed])
 
 
 # used for debugging; early on, the algorithm for quitting early had a bug and was causing the search
@@ -1548,9 +1132,9 @@ def _restart_from_directory(directories: _Directories, design: dc.Design) \
     # which are ignored. This applies to both the design and the RNG state
 
     # get DomainPools
-    pool_with_name = read_domain_pools(directories)
+    pool_with_name = _read_domain_pools(directories)
 
-    # returns highest index found, as well as name of corresponding design file
+    # returns highest index found in dsd_design subdirectory
     highest_idx_design = _find_highest_index_in_directory(directories.dsd_design,
                                                           directories.dsd_design_filename_no_ext, 'json')
     design_filename = os.path.join(directories.dsd_design,
@@ -1584,7 +1168,7 @@ def _restart_from_directory(directories: _Directories, design: dc.Design) \
     return highest_idx_design, rng
 
 
-def read_domain_pools(directories: _Directories) -> Dict[str, dc.DomainPool]:
+def _read_domain_pools(directories: _Directories) -> Dict[str, dc.DomainPool]:
     # return dict mapping name of DomainPool to DomainPool, read from highest index of file with certain
     # name in directories.domain_pool subdirectory
 
@@ -1843,8 +1427,9 @@ def assign_sequences_to_domains_randomly_from_pools(design: Design,
         if overwrite_existing_sequences:
             if not domain.fixed:
                 at_least_one_domain_unfixed = True
-                domain.sequence = domain.pool.generate_sequence(rng, domain.sequence)
-                assert len(domain.sequence) == domain.pool.length
+                new_sequence = domain.pool.generate_sequence(rng, domain.sequence())
+                domain.set_sequence(new_sequence)
+                assert len(domain.sequence()) == domain.pool.length
             else:
                 logger.info(skip_nonfixed_msg)
         else:
@@ -1853,8 +1438,9 @@ def assign_sequences_to_domains_randomly_from_pools(design: Design,
                 # domain is not fixed so that we know it is eligible to be overwritten during the search
                 at_least_one_domain_unfixed = True
             if not domain.fixed and not domain.has_sequence():
-                domain.sequence = domain.pool.generate_sequence(rng)
-                assert len(domain.sequence) == domain.pool.length
+                new_sequence = domain.pool.generate_sequence(rng)
+                domain.set_sequence(new_sequence)
+                assert len(domain.sequence()) == domain.pool.length
             elif warn_fixed_sequences:
                 if domain.fixed:
                     logger.info(skip_fixed_msg)
diff --git a/dsd/vienna_nupack.py b/dsd/vienna_nupack.py
index b467b35a..caabab9e 100644
--- a/dsd/vienna_nupack.py
+++ b/dsd/vienna_nupack.py
@@ -23,6 +23,8 @@
 
 import numpy as np
 
+import dsd.constraints as dc
+
 os_is_windows = sys.platform == 'win32'
 
 parameter_set_directory = 'nupack_viennaRNA_parameter_files'
@@ -129,10 +131,7 @@ def pfunc(seqs: Union[str, Tuple[str, ...]],
     return dg
 
 
-Complex = Tuple['Strand', ...]
-
-
-def nupack_complex_base_pair_probabilities(strand_complex: Complex,
+def nupack_complex_base_pair_probabilities(strand_complex: 'dc.Complex',
                                            temperature: float = default_temperature,
                                            sodium: float = default_sodium,
                                            magnesium: float = default_magnesium) -> np.ndarray:
@@ -408,8 +407,7 @@ def rna_cofold_multiple(seq_pairs: Sequence[Tuple[str, str]],
     for line in lines[:-1]:
         dg_list.append(-float(line.split(':')[1].split('(')[1].split(')')[0]))
     if len(lines) - 1 != len(seq_pairs):
-        raise AssertionError(
-            'lengths do not match: #lines:{} #seqpairs:{}'.format(len(lines) - 1, len(seq_pairs)))
+        raise AssertionError(f'lengths do not match: #lines:{len(lines) - 1} #seqpairs:{len(seq_pairs)}')
 
     return dg_list
 
@@ -422,7 +420,7 @@ def wc(seq: str) -> str:
     return seq.translate(_wctable)[::-1]
 
 
-def secondary_structure_single_strand(
+def complex_free_energy_single_strand(
         seq: str, temperature: float = default_temperature, sodium: float = default_sodium,
         magnesium: float = default_magnesium) -> float:
     """Computes the complex free energy of a single strand.
@@ -435,7 +433,7 @@ def secondary_structure_single_strand(
 
 def binding_complement(seq: str, temperature: float = default_temperature, sodium: float = default_sodium,
                        magnesium: float = default_magnesium, subtract_indv: bool = True) -> float:
-    """Computes the complex free energy of a strand with its perfect WC complement.
+    """Computes the complex free energy of a strand with its perfect Watson-Crick complement.
 
     NUPACK 4 must be installed. Installation instructions can be found at
     https://piercelab-caltech.github.io/nupack-docs/start/.
@@ -489,19 +487,19 @@ def log_energy(energy: float) -> None:
 
 def domain_orthogonal(seq: str, seqs: Sequence[str], temperature: float, sodium: float,
                       magnesium: float, orthogonality: float,
-                      orthogonality_ave: float = -1, threaded: bool = False) -> bool:
+                      orthogonality_ave: float = -1, parallel: bool = False) -> bool:
     """test orthogonality of domain with all others and their wc complements
 
     NUPACK 4 must be installed. Installation instructions can be found at
     https://piercelab-caltech.github.io/nupack-docs/start/.
     """
-    if threaded:
+    if parallel:
         raise NotImplementedError()
 
     def binding_callback(s1: str, s2: str) -> float:
         return binding(s1, s2, temperature=temperature, sodium=sodium, magnesium=magnesium)
 
-    if threaded:
+    if parallel:
         results = [
             global_thread_pool.apply_async(binding_callback, args=(s, s))
             for s in (seq, wc(seq))]
@@ -519,7 +517,7 @@ def binding_callback(s1: str, s2: str) -> float:
             return False
     energy_sum = 0.0
     for altseq in seqs:
-        if threaded:
+        if parallel:
             results = [
                 global_thread_pool.apply_async(binding_callback,
                                                args=(seq1, seq2, temperature, sodium, magnesium))
@@ -556,22 +554,22 @@ def binding_callback(s1: str, s2: str) -> float:
 def domain_pairwise_concatenated_no_sec_struct(seq: str, seqs: Sequence[str], temperature: float,
                                                sodium: float, magnesium: float,
                                                concat: float, concat_ave: float = -1,
-                                               threaded: bool = False) -> bool:
+                                               parallel: bool = False) -> bool:
     """test lack of secondary structure in concatenated domains
 
     NUPACK 4 must be installed. Installation instructions can be found at
     https://piercelab-caltech.github.io/nupack-docs/start/.
     """
 
-    if threaded:
+    if parallel:
         raise NotImplementedError()
 
     energy_sum = 0.0
     for altseq in seqs:
         wc_seq = wc(seq)
         wc_altseq = wc(altseq)
-        if threaded:
-            results = [global_thread_pool.apply_async(secondary_structure_single_strand,
+        if parallel:
+            results = [global_thread_pool.apply_async(complex_free_energy_single_strand,
                                                       args=(seq1 + seq2, temperature, sodium, magnesium)) for
                        (seq1, seq2) in
                        [(seq, altseq),
@@ -589,29 +587,29 @@ def domain_pairwise_concatenated_no_sec_struct(seq: str, seqs: Sequence[str], te
                 return False
             energy_sum += sum(energies)
         else:
-            seq_alt = secondary_structure_single_strand(seq + altseq, temperature, sodium, magnesium)
+            seq_alt = complex_free_energy_single_strand(seq + altseq, temperature, sodium, magnesium)
             if seq_alt > concat:
                 return False
-            seq_wcalt = secondary_structure_single_strand(seq + wc_altseq, temperature, sodium, magnesium)
+            seq_wcalt = complex_free_energy_single_strand(seq + wc_altseq, temperature, sodium, magnesium)
             if seq_wcalt > concat:
                 return False
-            wcseq_alt = secondary_structure_single_strand(wc_seq + altseq, temperature, sodium, magnesium)
+            wcseq_alt = complex_free_energy_single_strand(wc_seq + altseq, temperature, sodium, magnesium)
             if wcseq_alt > concat:
                 return False
-            wcseq_wcalt = secondary_structure_single_strand(wc_seq + wc_altseq, temperature, sodium,
+            wcseq_wcalt = complex_free_energy_single_strand(wc_seq + wc_altseq, temperature, sodium,
                                                             magnesium)
             if wcseq_wcalt > concat:
                 return False
-            alt_seq = secondary_structure_single_strand(altseq + seq, temperature, sodium, magnesium)
+            alt_seq = complex_free_energy_single_strand(altseq + seq, temperature, sodium, magnesium)
             if alt_seq > concat:
                 return False
-            alt_wcseq = secondary_structure_single_strand(altseq + wc_seq, temperature, sodium, magnesium)
+            alt_wcseq = complex_free_energy_single_strand(altseq + wc_seq, temperature, sodium, magnesium)
             if alt_wcseq > concat:
                 return False
-            wcalt_seq = secondary_structure_single_strand(wc_altseq + seq, temperature, sodium, magnesium)
+            wcalt_seq = complex_free_energy_single_strand(wc_altseq + seq, temperature, sodium, magnesium)
             if wcalt_seq > concat:
                 return False
-            wcalt_wcseq = secondary_structure_single_strand(wc_altseq + wc_seq, temperature, sodium,
+            wcalt_wcseq = complex_free_energy_single_strand(wc_altseq + wc_seq, temperature, sodium,
                                                             magnesium)
             if wcalt_wcseq > concat:
                 return False
diff --git a/examples/many_strands_no_common_domains.py b/examples/many_strands_no_common_domains.py
index 23867fa9..6d133e87 100644
--- a/examples/many_strands_no_common_domains.py
+++ b/examples/many_strands_no_common_domains.py
@@ -7,7 +7,7 @@
 import dsd.constraints as dc  # type: ignore
 import dsd.vienna_nupack as dv  # type: ignore
 import dsd.search as ds  # type: ignore
-from dsd.constraints import NumpyConstraint, SequenceConstraint
+from dsd.constraints import NumpyConstraint
 
 
 # command-line arguments
@@ -69,41 +69,44 @@ def main() -> None:
     # strand i is    [----------|----------|----------|---------->
     strands = [dc.Strand([f's{i}', f'w{i}', f'n{i}', f'e{i}']) for i in range(num_strands)]
 
-    # fix all domains of strand 0 and one domain of strand 1
-    for domain in strands[0].domains:
-        domain.set_fixed_sequence('ACGTACGTAC')
-    strands[1].domains[0].set_fixed_sequence('ACGTACGTAC')
+    some_fixed = False
+    if some_fixed:
+        # fix all domains of strand 0 and one domain of strand 1
+        for domain in strands[0].domains:
+            domain.set_fixed_sequence('ACGTACGTAC')
+        strands[1].domains[0].set_fixed_sequence('ACGTACGTAC')
 
-    threaded = False
-    # threaded = True
+    parallel = False
+    # parallel = True
+
+    domain_nupack_ss_constraint = dc.nupack_domain_complex_free_energy_constraint(
+        threshold=-0.0, temperature=52, short_description='DomainSS')
 
     domain_pairs_rna_duplex_constraint = dc.rna_duplex_domain_pairs_constraint(
         threshold=-2.0, temperature=52, short_description='DomainPairNoCompl')
 
     domain_pair_nupack_constraint = dc.nupack_domain_pair_constraint(
         threshold=-0.5, temperature=52, short_description='DomainPairNoCompl',
-        threaded=threaded)
+        parallel=parallel)
 
     strand_pairs_rna_duplex_constraint = dc.rna_duplex_strand_pairs_constraint(
-        threshold=-1.0, temperature=52, short_description='StrandPairNoCompl', threaded=threaded)
-
-    # strand_individual_ss_constraint = dc.nupack_strand_secondary_structure_constraint(
-    #     threshold=-1.5, temperature=52, short_description='StrandSS', threaded=threaded)
+        threshold=-1.0, temperature=52, short_description='StrandPairNoCompl', parallel=parallel)
 
-    strand_individual_ss_constraint = dc.nupack_strand_secondary_structure_constraint(
-        threshold=-1.0, temperature=52, short_description='StrandSS', threaded=threaded)
+    strand_individual_ss_constraint = dc.nupack_strand_complex_free_energy_constraint(
+        threshold=-1.0, temperature=52, short_description='StrandSS', parallel=parallel)
 
     strand_pair_nupack_constraint = dc.nupack_strand_pair_constraint(
-        threshold=-3.5, temperature=52, short_description='StrandPairNoCompl', threaded=threaded, weight=0.1)
+        threshold=-3.5, temperature=52, short_description='StrandPairNoCompl', parallel=parallel, weight=0.1)
 
     design = dc.Design(strands,
                        constraints=[
-                           # strand_individual_ss_constraint,
-                           # strand_pair_nupack_constraint,
-                           # domain_pair_nupack_constraint,
+                           domain_nupack_ss_constraint,
+                           strand_individual_ss_constraint,
+                           strand_pair_nupack_constraint,
+                           domain_pair_nupack_constraint,
                            domain_pairs_rna_duplex_constraint,
+                           strand_pairs_rna_duplex_constraint,
                            # dc.domains_not_substrings_of_each_other_domain_pair_constraint(),
-                           # strand_pairs_rna_duplex_constraint,
                        ])
 
     numpy_constraints: List[NumpyConstraint] = [
@@ -140,19 +143,27 @@ def main() -> None:
                                    replace_with_close_sequences=replace_with_close_sequences,
                                    )
 
-    for strand in strands[1:]:  # skip all domains on strand 0 since all its domains are fixed
-        for domain in strand.domains[:2]:
-            if domain.name != 's1':  # skip for s1 since that domain is fixed
+    if some_fixed:
+        for strand in strands[1:]:  # skip all domains on strand 0 since all its domains are fixed
+            for domain in strand.domains[:2]:
+                if domain.name != 's1':  # skip for s1 since that domain is fixed
+                    domain.pool = domain_pool_10
+            for domain in strand.domains[2:]:
+                domain.pool = domain_pool_11
+    else:
+        for strand in strands:
+            for domain in strand.domains[:2]:
                 domain.pool = domain_pool_10
-        for domain in strand.domains[2:]:
-            domain.pool = domain_pool_11
+            for domain in strand.domains[2:]:
+                domain.pool = domain_pool_11
 
-    # have to set nupack_complex_secondary_structure_constraint after DomainPools are set
-    strand_complexes = [(strand,) for i,strand in enumerate(strands[2:])]
-    strand_base_pair_prob_constraint = dc.nupack_complex_secondary_structure_constraint(
+    # have to set nupack_complex_secondary_structure_constraint after DomainPools are set,
+    # so that we know the domain lengths
+    strand_complexes = [dc.Complex((strand,)) for i, strand in enumerate(strands[2:])]
+    strand_base_pair_prob_constraint = dc.nupack_complex_base_pair_probability_constraint(
         strand_complexes=strand_complexes)
 
-    design.add_constraints([strand_base_pair_prob_constraint])
+    # design.add_constraints([strand_base_pair_prob_constraint])
 
     params = ds.SearchParameters(out_directory=args.directory,
                                  restart=args.restart,
diff --git a/examples/sample_designer.py b/examples/sample_designer.py
index cd8d75b1..5f7a1de5 100644
--- a/examples/sample_designer.py
+++ b/examples/sample_designer.py
@@ -88,7 +88,7 @@ def main() -> None:
         threshold=-7.0, temperature=52, short_description='StrandPairCompl',
         pairs=[(strand0, strand2), (strand0, strand3), (strand1, strand2), (strand1, strand3)]
     )
-    strand_individual_ss_constraint = dc.nupack_strand_secondary_structure_constraint(
+    strand_individual_ss_constraint = dc.nupack_strand_complex_free_energy_constraint(
         threshold=-0.0, temperature=52, short_description='StrandSS')
 
     initial_design = dc.Design(strands,
diff --git a/examples/seesaw_gate.py b/examples/seesaw_gate.py
index 8ce34982..a61ea072 100644
--- a/examples/seesaw_gate.py
+++ b/examples/seesaw_gate.py
@@ -223,7 +223,7 @@ def reporter_base_strand(gate) -> dc.Strand:
 g_5_s_5_7_complex = (signal_5_7_strand, gate_5_base_strand)
 
 g_5_s_5_6_nonimplicit_base_pairs = [(signal_5_6_toehold_addr, gate_5_bound_toehold_3p_addr)]
-g_5_s_5_6_complex_constraint = dc.nupack_complex_secondary_structure_constraint(
+g_5_s_5_6_complex_constraint = dc.nupack_complex_base_pair_probability_constraint(
     strand_complexes=[g_5_s_5_6_complex, g_5_s_5_7_complex],
     nonimplicit_base_pairs=g_5_s_5_6_nonimplicit_base_pairs,
     # base_pair_probs={
@@ -257,7 +257,7 @@ def reporter_base_strand(gate) -> dc.Strand:
 #                      T*         S5*      s5*   T*
 g_5_s_2_5_complex = (signal_2_5_strand, gate_5_base_strand)
 g_5_s_2_5_nonimplicit_base_pairs = [(signal_2_5_toehold_addr, gate_5_bound_toehold_5p_addr)]
-g_5_s_2_5_complex_constraint = dc.nupack_complex_secondary_structure_constraint(
+g_5_s_2_5_complex_constraint = dc.nupack_complex_base_pair_probability_constraint(
     strand_complexes=[g_5_s_2_5_complex],
     nonimplicit_base_pairs=g_5_s_2_5_nonimplicit_base_pairs,
     # base_pair_probs={
@@ -293,7 +293,7 @@ def reporter_base_strand(gate) -> dc.Strand:
 #        INTERIOR_TO_STRAND   BLUNT_END
 #  s2*   T*        S5*       s5*
 t_2_5_w_5_complex = (waste_5_strand, threshold_2_5_base_strand)
-t_2_5_w_5_complex_constraint = dc.nupack_complex_secondary_structure_constraint(
+t_2_5_w_5_complex_constraint = dc.nupack_complex_base_pair_probability_constraint(
     strand_complexes=[t_2_5_w_5_complex])
 
 #
@@ -324,7 +324,7 @@ def reporter_base_strand(gate) -> dc.Strand:
 #                s2*   T*        S5*       s5*
 
 waste_2_5_complex = (signal_2_5_strand, threshold_2_5_base_strand)
-waste_2_5_complex_constraint = dc.nupack_complex_secondary_structure_constraint(
+waste_2_5_complex_constraint = dc.nupack_complex_base_pair_probability_constraint(
     strand_complexes=[waste_2_5_complex])
 
 #               S6       s6
@@ -350,7 +350,7 @@ def reporter_base_strand(gate) -> dc.Strand:
 #    INTERIOR_TO_STRAND   BLUNT_END
 #    T*        S6*       s6*
 reporter_6_complex = (waste_6_strand, reporter_6_base_strand)
-reporter_6_complex_constraint = dc.nupack_complex_secondary_structure_constraint(
+reporter_6_complex_constraint = dc.nupack_complex_base_pair_probability_constraint(
     strand_complexes=[reporter_6_complex])
 
 #       S5        s5    T          S6       s6
@@ -377,7 +377,7 @@ def reporter_base_strand(gate) -> dc.Strand:
 #                        INTERIOR_TO_STRAND  BLUNT_END
 #                       T*        S6*       s6*
 f_waste_6_complex = (signal_5_6_strand, reporter_6_base_strand)
-f_waste_6_complex_constraint = dc.nupack_complex_secondary_structure_constraint(
+f_waste_6_complex_constraint = dc.nupack_complex_base_pair_probability_constraint(
     strand_complexes=[f_waste_6_complex])
 
 
diff --git a/examples/square_root_circuit.py b/examples/square_root_circuit.py
index 8a34f60f..480dbbc9 100644
--- a/examples/square_root_circuit.py
+++ b/examples/square_root_circuit.py
@@ -43,7 +43,7 @@
     [three_letter_code_constraint, c_content_constraint, no_aaaaa_constraint, no_gggg_constraint])
 
 # Alias
-dc_complex_constraint = dc.nupack_complex_secondary_structure_constraint
+dc_complex_constraint = dc.nupack_complex_base_pair_probability_constraint
 
 # Stores all domains used in design
 TOEHOLD_DOMAIN: dc.Domain = dc.Domain('T', pool=TOEHOLD_DOMAIN_POOL)


Test patch:
diff --git a/examples/hamming_dist_test.py b/examples/hamming_dist_test.py
index 66d5c408..d04b2143 100644
--- a/examples/hamming_dist_test.py
+++ b/examples/hamming_dist_test.py
@@ -25,8 +25,8 @@ def main():
                                  max_iterations=None)
     params.force_overwrite = True  # directly deletes output folder contents w/o user input
     # params.report_delay = 0.0
-    design = dc.Design(strands, constraints=[dc.nupack_strand_secondary_structure_constraint(
-        threshold=-1.5, temperature=52, short_description='StrandSS', threaded=False)])
+    design = dc.Design(strands, constraints=[dc.nupack_strand_complex_free_energy_constraint(
+        threshold=-1.5, temperature=52, short_description='StrandSS', parallel=False)])
     ds.search_for_dna_sequences(design, params)
 
 
diff --git a/examples/limit_test.py b/examples/limit_test.py
index fc0c98b7..02a36d91 100644
--- a/examples/limit_test.py
+++ b/examples/limit_test.py
@@ -27,8 +27,8 @@ def main():
                                  max_iterations=None)
     params.force_overwrite = True  # directly deletes output folder contents w/o user input
     # params.report_delay = 0.0
-    design = dc.Design(strands, constraints=[dc.nupack_strand_secondary_structure_constraint(
-        threshold=-1.5, temperature=52, short_description='StrandSS', threaded=False)])
+    design = dc.Design(strands, constraints=[dc.nupack_strand_complex_free_energy_constraint(
+        threshold=-1.5, temperature=52, short_description='StrandSS', parallel=False)])
     ds.search_for_dna_sequences(design, params)
 
 
diff --git a/tests/test.py b/tests/test.py
index 110067c4..db8e7018 100644
--- a/tests/test.py
+++ b/tests/test.py
@@ -930,14 +930,14 @@ def test_assign_dna_sequence_to_parent(self):
         domains = self.sample_nested_domains()
         sequence = 'CATAGCTTTCTTGTTCTGATCGGAAC'
         a = domains['a']
-        a.sequence = sequence
-        self.assertEqual(sequence, a.sequence)
-        self.assertEqual(sequence[0: 11], domains['b'].sequence)
-        self.assertEqual(sequence[11:], domains['C'].sequence)
-        self.assertEqual(sequence[0:5], domains['E'].sequence)
-        self.assertEqual(sequence[5:11], domains['F'].sequence)
-        self.assertEqual(sequence[11:18], domains['g'].sequence)
-        self.assertEqual(sequence[18:], domains['h'].sequence)
+        a.set_sequence(sequence)
+        self.assertEqual(sequence, a.sequence())
+        self.assertEqual(sequence[0: 11], domains['b'].sequence())
+        self.assertEqual(sequence[11:], domains['C'].sequence())
+        self.assertEqual(sequence[0:5], domains['E'].sequence())
+        self.assertEqual(sequence[5:11], domains['F'].sequence())
+        self.assertEqual(sequence[11:18], domains['g'].sequence())
+        self.assertEqual(sequence[18:], domains['h'].sequence())
 
     def test_assign_dna_sequence_to_leaf(self):
         """
@@ -954,11 +954,11 @@ def test_assign_dna_sequence_to_leaf(self):
         domains = self.sample_nested_domains()
         E = domains['E']
         F = domains['F']
-        E.sequence = 'CATAG'
-        F.sequence = 'CTTTCC'
-        self.assertEqual('CATAG', E.sequence)
-        self.assertEqual('CTTTCC', F.sequence)
-        self.assertEqual('CATAGCTTTCC', domains['b'].sequence)
+        E.set_sequence('CATAG')
+        F.set_sequence('CTTTCC')
+        self.assertEqual('CATAG', E.sequence())
+        self.assertEqual('CTTTCC', F.sequence())
+        self.assertEqual('CATAGCTTTCC', domains['b'].sequence())
 
     def test_assign_dna_sequence_mixed(self):
         """
@@ -976,38 +976,38 @@ def test_assign_dna_sequence_mixed(self):
         E = domains['E']
         F = domains['F']
         C = domains['C']
-        E.sequence = 'CATAG'
-        F.sequence = 'CTTTCT'
-        C.sequence = 'TGTTCTGATCGGAAC'
+        E.set_sequence('CATAG')
+        F.set_sequence('CTTTCT')
+        C.set_sequence('TGTTCTGATCGGAAC')
 
         # Assert initial assignment is correct
-        self.assertEqual('CATAG''CTTTCT''TGTTCTGATCGGAAC', domains['a'].sequence)
-        self.assertEqual('CATAG''CTTTCT', domains['b'].sequence)
-        self.assertEqual('TGTTCTGATCGGAAC', domains['C'].sequence)
-        self.assertEqual('CATAG', domains['E'].sequence)
-        self.assertEqual('CTTTCT', domains['F'].sequence)
-        self.assertEqual('TGTTCTG', domains['g'].sequence)
-        self.assertEqual('ATCGGAAC', domains['h'].sequence)
+        self.assertEqual('CATAG''CTTTCT''TGTTCTGATCGGAAC', domains['a'].sequence())
+        self.assertEqual('CATAG''CTTTCT', domains['b'].sequence())
+        self.assertEqual('TGTTCTGATCGGAAC', domains['C'].sequence())
+        self.assertEqual('CATAG', domains['E'].sequence())
+        self.assertEqual('CTTTCT', domains['F'].sequence())
+        self.assertEqual('TGTTCTG', domains['g'].sequence())
+        self.assertEqual('ATCGGAAC', domains['h'].sequence())
 
         # Assert subsequent reassignment to leaf is correct
-        F.sequence = 'ATGTTT'
-        self.assertEqual('CATAG''ATGTTT''TGTTCTGATCGGAAC', domains['a'].sequence)
-        self.assertEqual('CATAG''ATGTTT', domains['b'].sequence)
-        self.assertEqual('TGTTCTGATCGGAAC', domains['C'].sequence)
-        self.assertEqual('CATAG', domains['E'].sequence)
-        self.assertEqual('ATGTTT', domains['F'].sequence)
-        self.assertEqual('TGTTCTG', domains['g'].sequence)
-        self.assertEqual('ATCGGAAC', domains['h'].sequence)
+        F.set_sequence('ATGTTT')
+        self.assertEqual('CATAG''ATGTTT''TGTTCTGATCGGAAC', domains['a'].sequence())
+        self.assertEqual('CATAG''ATGTTT', domains['b'].sequence())
+        self.assertEqual('TGTTCTGATCGGAAC', domains['C'].sequence())
+        self.assertEqual('CATAG', domains['E'].sequence())
+        self.assertEqual('ATGTTT', domains['F'].sequence())
+        self.assertEqual('TGTTCTG', domains['g'].sequence())
+        self.assertEqual('ATCGGAAC', domains['h'].sequence())
 
         # Assert subsequent reassignment to internal node is correct
-        C.sequence = 'GGGGGGGGGGGGGGG'
-        self.assertEqual('CATAG''ATGTTT''GGGGGGGGGGGGGGG', domains['a'].sequence)
-        self.assertEqual('CATAG''ATGTTT', domains['b'].sequence)
-        self.assertEqual('GGGGGGGGGGGGGGG', domains['C'].sequence)
-        self.assertEqual('CATAG', domains['E'].sequence)
-        self.assertEqual('ATGTTT', domains['F'].sequence)
-        self.assertEqual('GGGGGGG', domains['g'].sequence)
-        self.assertEqual('GGGGGGGG', domains['h'].sequence)
+        C.set_sequence('GGGGGGGGGGGGGGG')
+        self.assertEqual('CATAG''ATGTTT''GGGGGGGGGGGGGGG', domains['a'].sequence())
+        self.assertEqual('CATAG''ATGTTT', domains['b'].sequence())
+        self.assertEqual('GGGGGGGGGGGGGGG', domains['C'].sequence())
+        self.assertEqual('CATAG', domains['E'].sequence())
+        self.assertEqual('ATGTTT', domains['F'].sequence())
+        self.assertEqual('GGGGGGG', domains['g'].sequence())
+        self.assertEqual('GGGGGGGG', domains['h'].sequence())
 
     def test_error_assign_dna_sequence_to_parent_with_incorrect_size_subdomain(self):
         """
@@ -1025,7 +1025,7 @@ def test_error_assign_dna_sequence_to_parent_with_incorrect_size_subdomain(self)
 
         a: Domain = Domain('a', assign_domain_pool_of_size(15), dependent=True, subdomains=[B, C])
         with self.assertRaises(ValueError):
-            a.sequence = 'A' * 15
+            a.set_sequence('A' * 15)
 
     def test_construct_strand_using_dependent_subdomain(self) -> None:
         """Test constructing a strand using a dependent subdomain (not parent)