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`).
remove pre-computation of fixed set of sequences in DomainPool
Currently, we aren't getting the full benefit from the option `DomainPool.replace_with_close_sequences`. The reason is that it does not pick a sequence uniformly at random from *all* those of the correct length, and crucially also satisfying NumpyConstraints, but instead from those (perhaps much smaller number) that were last generated when calling `_reset_precomputed_sequences`.

One could imagine there is a "easy to find through local search" downhill path from a random sequence to a very particular sequence (i.e., one very unlikely to be included in any given call to `_reset_precomputed_sequences`).

Here's a better idea for how to implement `DomainPool.replace_with_close_sequences`, which allows any sequence with the correct length and satisfying the NumpyConstraints to be picked, yet that samples uniformly at random from all such sequences that are a given Hamming distance from the previous sequence.

First, generate 100,000 sequences when the `DomainPool` is first created, not to store, but merely to estimate a success probability of passing the `NumpyConstraint`'s. Sometimes the constraints are very light, and this probably will be large, like 1/2 or 1/3. In such a case, there's little point in running large numbers of sequences through the `NumpyConstraint`'s. Generating sequences one at a time, and then checking against all `NumpyConstraint`'s, will be fast enough, because we expected on average only to generate a few sequences before one passes. That said, perhaps this needs to be more dynamic. For example, we might have a sequence where most of its neighbors would not pass, so the success probably of its neighborhood is much lower than average. But perhaps this can be handled by dynamically adjusting. For example, we always start with generating just a small number of sequences, but increase their number by a constant factor each time we fail to find any satisfying the `NumpyConstraint`'s. Hopefully this leads to fairly quickly finding at least one.

Hmm, now I'm not sure where to go with this. Leaving my thoughts for now, will come back and edit later.
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.
86 remove pre computation of fixed set of sequences in domain pool

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`).


Discussion / hints:

Fix patch:
diff --git a/dsd/constraints.py b/dsd/constraints.py
index 373f2b00..b03b09ac 100644
--- a/dsd/constraints.py
+++ b/dsd/constraints.py
@@ -63,9 +63,6 @@
 length_key = 'length'
 strand_name_in_strand_pool_key = 'strand_name'
 sequences_key = 'sequences'
-num_sampled_key = 'num_sampled'
-max_samples_key = 'max_samples'
-num_times_sequences_reset_key = 'num_times_sequences_reset'
 replace_with_close_sequences_key = 'replace_with_close_sequences'
 hamming_probability_key = 'hamming_probability'
 num_sequences_to_generate_upper_limit_key = 'num_sequences_to_generate_upper_limit'
@@ -84,9 +81,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 +118,6 @@ def _configure_logger() -> None:
 
 _configure_logger()
 
-T = TypeVar('T')
-
 
 def all_pairs(values: Iterable[T],
               with_replacement: bool = True,
@@ -546,6 +541,11 @@ def remove_violating_sequences(self, seqs: dn.DNASeqList) -> dn.DNASeqList:
         return constraint.remove_violating_sequences(seqs)
 
 
+log_numpy_generation = True
+
+
+# log_numpy_generation = False
+
 @dataclass
 class DomainPool:
     """
@@ -563,6 +563,23 @@ class DomainPool:
     length: int
     """Length of DNA sequences generated by this :any:`DomainPool`."""
 
+    replace_with_close_sequences: bool = False
+    """
+    If True, instead of picking a sequence uniformly at random from all those satisfying the constraints
+    when returning a sequence from :meth:`DomainPool.generate_sequence`,
+    one is picked "close" in Hamming distance to the previous sequence of the :any:`Domain`.
+    The field :data:`DomainPool.hamming_probability` is used to pick a distance at random, after which
+    a sequence that distance from the previous sequence is selected to return from .
+    """
+
+    hamming_probability: Dict[int, float] = field(default_factory=dict)
+    """
+    Dictionary that specifies probability of taking a new sequence from the pool that is some integer 
+    number of bases different from the previous sequence (Hamming distance). 
+    """
+
+
+
     numpy_constraints: List[NumpyConstraint] = field(
         compare=False, hash=False, default_factory=list, repr=False)
     """
@@ -594,54 +611,22 @@ class DomainPool:
     Optional; default is empty.
     """
 
-    sequences: dn.DNASeqList = field(compare=False, hash=False,
-                                     default_factory=lambda: dn.DNASeqList(length=0), repr=False)
-    # list of available sequences.
-    # We sample with replacement (uniform if :data:`DomainPool.replace_with_close_sequences` = False,
-    # according to :data:`DomainPool.hamming_probability` otherwise)
-    # This is represented as a DNASeqList for efficiency when calculating Hamming distances
-    # when the option `replace_with_close_sequences` is True.
-
-    num_sampled: int = 0
-
-    max_samples: Optional[int] = None
-
-    num_times_sequences_reset: int = 0
-
-    replace_with_close_sequences: bool = False
-    """
-    If True, instead of picking a sequence uniformly at random from all those satisfying the constraints
-    when returning a sequence from :meth:`DomainPool.generate_sequence`,
-    one is picked "close" in Hamming distance to the previous sequence of the :any:`Domain`.
-    The field :data:`DomainPool.hamming_probability` is used to pick a distance at random, after which
-    a sequence that distance from the previous sequence is selected to return from .
-    """
-
-    hamming_probability: Dict[int, float] = field(default_factory=dict)
-    """
-    Dictionary that specifies probability of taking a new sequence from the pool that is some integer 
-    number of bases different from the previous sequence (Hamming distance). 
-    """
-
-    num_sequences_to_generate_upper_limit: int = 4 ** 12
-    """
-    Maximum number of sequences to generate at once. If :py:data:`DomainPool.num_sequences_to_generate` 
-    exceeds this number, the user is asked to re-adjust constraint strictness. The number may be changed to 
-    satisfy individual device capabilities.  
-    """
-
-    num_sequences_to_generate: int = 10 ** 5
-    """
-    Number of sequences to be generated at once. Will be increased if no sequences satisfying constraints 
-    are found after generating the current specified amount of sequences. 
-    """
-
     def __post_init__(self) -> None:
         if len(self.hamming_probability) == 0:  # sets default probability distribution if the user does not
+            # exponentially decreasing probability of making i+1 (since i starts at 0) base changes
+            # for i in range(self.length):
+            #     self.hamming_probability[i + 1] = 1 / 2 ** (i + 1)
+            # self.hamming_probability[self.length] *= 2
+
+            # linearly decreasing probability of making i+1 (since i starts at 0) base changes
+            total = 0.0
             for i in range(self.length):
-                # exponentially decreasing probability of making i+1 (since i starts at 0) base changes
-                self.hamming_probability[i + 1] = 1 / 2 ** (i + 1)
-            self.hamming_probability[self.length] *= 2
+                prob = 1 / (i + 1)
+                self.hamming_probability[i + 1] = prob
+                total += prob
+            # normalize to be a probability measure
+            for length in self.hamming_probability:
+                self.hamming_probability[length] /= total
 
         idx = 0
         for numpy_constraint in self.numpy_constraints:
@@ -662,83 +647,49 @@ def __post_init__(self) -> None:
                                  f'but the element at index {idx} is of type {type(seq_constraint)}')
             idx += 1
 
-    def to_json(self, include_sequences: bool) -> str:
-        json_map = self.to_json_serializable(include_sequences)
+    def __hash__(self) -> int:
+        return hash((self.name, self.length))
+
+    def __eq__(self, other: Any) -> bool:
+        if not isinstance(other, DomainPool):
+            return False
+        return self.name == other.name and self.length == other.length
+
+    def to_json(self) -> str:
+        json_map = self.to_json_serializable()
         json_str = json.dumps(json_map, indent=2)
         return json_str
 
-    def to_json_serializable(self, include_sequences: bool) -> Dict[str, Any]:
+    def to_json_serializable(self, suppress_indent: bool) -> Dict[str, Any]:
         dct = {
             name_key: self.name,
             length_key: self.length,
-            num_sampled_key: self.num_sampled,
-            max_samples_key: self.max_samples,
-            num_times_sequences_reset_key: self.num_times_sequences_reset,
             replace_with_close_sequences_key: self.replace_with_close_sequences,
             hamming_probability_key: self.hamming_probability,
-            num_sequences_to_generate_upper_limit_key: self.num_sequences_to_generate_upper_limit,
-            num_sequences_to_generate_key: self.num_sequences_to_generate
         }
-        if include_sequences:
-            dct[sequences_key] = self.sequences.to_list()
         return dct
 
     @staticmethod
     def from_json_serializable(json_map: Dict[str, Any]) -> DomainPool:
         name = json_map[name_key]
         length = json_map[length_key]
-        num_sampled = json_map[num_sampled_key]
-        max_samples = json_map[max_samples_key]
-        num_times_sequences_reset = json_map[num_times_sequences_reset_key]
         replace_with_close_sequences = json_map[replace_with_close_sequences_key]
         hamming_probability_str_keys = json_map[hamming_probability_key]
         hamming_probability = {int(key): val for key, val in hamming_probability_str_keys.items()}
-        num_sequences_to_generate_upper_limit = json_map[num_sequences_to_generate_upper_limit_key]
-        num_sequences_to_generate = json_map[num_sequences_to_generate_key]
-        sequences_list = json_map[sequences_key]
-        sequences = dn.DNASeqList(seqs=sequences_list)
-        return DomainPool(name=name, length=length, sequences=sequences,
-                          num_sampled=num_sampled, max_samples=max_samples,
-                          num_times_sequences_reset=num_times_sequences_reset,
+        return DomainPool(name=name, length=length,
                           replace_with_close_sequences=replace_with_close_sequences,
                           hamming_probability=hamming_probability,
-                          num_sequences_to_generate_upper_limit=num_sequences_to_generate_upper_limit,
-                          num_sequences_to_generate=num_sequences_to_generate
                           )
 
-    def _reset_precomputed_sequences(self, rng: np.random.Generator) -> None:
-        # precomputes a new list of sequences satisfying numpy constraints and sequence constraints
-        self.sequences, use_random_subset = self._generate_sequences_satisfying_numpy_constraints(rng)
-        self._filter_sequence_constraints()
-        self.num_sampled = 0
-        self.num_times_sequences_reset += 1
-
-        if use_random_subset:
-            # sampling 2 * self.sequences.numseqs uniformly at random will leave 1/e^2 fraction unsampled
-            # in expectation: https://www.cs.purdue.edu/homes/hmaji/teaching/Spring%202017/lectures/03.pdf
-            self.max_samples = 2 * self.sequences.numseqs
-        else:
-            # if we didn't pick a random subset, we have all the sequences of this length satisfying
-            # the constraints, so never reset and just keep sampling them with replacement forever
-            self.max_samples = None
-
-    def _filter_sequence_constraints(self) -> None:
+    def _filter_sequence_constraints(self, seqs: dn.DNASeqList) -> None:
         if len(self.sequence_constraints) == 0:
             return
         idxs_to_keep = []
-        for idx in range(self.sequences.numseqs):
-            seq = self.sequences.get_seq_str(idx)
+        for idx in range(seqs.numseqs):
+            seq = seqs.get_seq_str(idx)
             if self.satisfies_sequence_constraints(seq):
                 idxs_to_keep.append(idx)
-        self.sequences.keep_seqs_at_indices(idxs_to_keep)
-
-    def __hash__(self) -> int:
-        return hash((self.name, self.length))
-
-    def __eq__(self, other: Any) -> bool:
-        if not isinstance(other, DomainPool):
-            return False
-        return self.name == other.name and self.length == other.length
+        seqs.keep_seqs_at_indices(idxs_to_keep)
 
     def satisfies_sequence_constraints(self, sequence: str) -> bool:
         """
@@ -747,34 +698,6 @@ def satisfies_sequence_constraints(self, sequence: str) -> bool:
         """
         return all(constraint(sequence) for constraint in self.sequence_constraints)
 
-    def find_hamming_distances(self, previous_sequence: str) -> Dict[int, dn.DNASeqList]:
-        """
-        Makes a dictionary mapping a Hamming distance to a :any:`DNASeqList`, corresponding to how many bases
-        different the sequence is from `previous_sequence` (i.e., the Hamming distance).
-
-        :param previous_sequence:
-            DNA sequence to be find Hamming distances from
-        :return:
-            dictionary mapping Hamming distances to sequences (represented as :any:`DNASeqList`)
-            that are that Hamming distance from `previous_sequence`
-        """
-        return self.sequences.hamming_map(previous_sequence)
-
-    def sequences_at_hamming_distance(self, previous_sequence: str,
-                                      distance: int) -> dn.DNASeqList:
-        """
-        Returns a :any:`DNASeqList`, which are `distance` bases
-        different from `previous_sequence` (i.e., the Hamming distance).
-
-        :param previous_sequence:
-            DNA sequence to be find Hamming distances from
-        :param distance:
-            target Hamming distance
-        :return:
-            :any:`DNASeqList` of sequences that are that Hamming distance from `previous_sequence`
-        """
-        return self.sequences.sequences_at_hamming_distance(previous_sequence, distance)
-
     def generate_sequence(self, rng: np.random.Generator, previous_sequence: Optional[str] = None) -> str:
         """
         Returns a DNA sequence of given length satisfying :py:data:`DomainPool.numpy_constraints` and
@@ -800,135 +723,148 @@ def generate_sequence(self, rng: np.random.Generator, previous_sequence: Optiona
             DNA sequence of given length satisfying :py:data:`DomainPool.numpy_constraints` and
             :py:data:`DomainPool.sequence_constraints`
         """
+        # import time
+        # b = time.perf_counter_ns()
         if not self.replace_with_close_sequences or previous_sequence is None:
-            # takes a completely random sequence from domain pool
             sequence = self._get_next_sequence_satisfying_numpy_and_sequence_constraints(rng)
         else:
-            # import time
-            # before = time.perf_counter_ns()
-            sequence = self._sample_hamming_distance_from_sequence(previous_sequence, rng, True)
-            # after = time.perf_counter_ns()
-            # print(f'time spent sampling Hamming neighbors: {(after - before) / 1e6:.1f} ms')
-
-        self.num_sampled += 1
+            sequence = self._sample_hamming_distance_from_sequence(previous_sequence, rng)
+        # a = time.perf_counter_ns()
+        # ms = (a-b)/10e6
+        # print(f'sample one seq time: {ms} ms')
 
         return sequence
 
-    def _sample_hamming_distance_from_sequence(self, previous_sequence: str, rng: np.random.Generator,
-                                               pick_distance_first: bool) -> str:
+    def _sample_hamming_distance_from_sequence(self, previous_sequence: str, rng: np.random.Generator) -> str:
         # all possible distances from 1 to len(previous_sequence) are calculated.
-        # if pick_distance_first is true, we save some time by not gathering sequences at all distances
-        # into a dict. Instead we first sample a distance and then look up only those sequences at that
-        # distance to sample. The tradeoff is that if there are no sequences at that distance, we have
-        # to re-sample. Given that we now keep every sequence in the DomainPool available (i.e.,
-        # sampling with replacement) until we regenerate them all, hopefully this is fairly efficient.
-
-        if self.sequences.numseqs == 0 or (
-                self.max_samples is not None and self.num_sampled >= self.max_samples):
-            if self.sequences.numseqs > 0:
-                logger.info('Twice as many pool sequences have been sampled with replacement from a '
-                            'randomly chosen subset. Regenerating fresh sequences.')
-            self._reset_precomputed_sequences(rng)
-        # takes neighbor to previous sequence; difference in bases randomly chosen
+
         hamming_probabilities = np.array(list(self.hamming_probability.values()))
 
-        if pick_distance_first:
-            # pick a distance at random, then re-pick if no sequences are at that distance
-            sequence = None
-            previous_sequence_1d_array = dn.seq2arr(previous_sequence)
-            available_distances_list = list(range(1, len(previous_sequence) + 1))
-
-            # For efficiency we inline the logic of self.sequences_at_hamming_distance()
-            # The next line is the most expensive part of the calculation, so we only do it once,
-            # rather than calling self.sequences_at_hamming_distance() repeatedly, which would recalculate
-            # the next line for each sampled distance even though it doesn't change.
-            computed_distances = np.sum(
-                np.bitwise_xor(self.sequences.seqarr, previous_sequence_1d_array) != 0, axis=1)
-
-            while sequence is None:
-                available_distances_arr = np.array(available_distances_list)
-                existing_hamming_probabilities = hamming_probabilities[available_distances_arr - 1]
-                prob_sum = existing_hamming_probabilities.sum()
-                existing_hamming_probabilities /= prob_sum
-                sampled_distance = rng.choice(available_distances_arr, p=existing_hamming_probabilities)
-
-                # sequences = self.sequences_at_hamming_distance(previous_sequence, distance)
-                indices_at_distance = computed_distances == sampled_distance
-                arr = self.sequences.seqarr[indices_at_distance]
-                sequences = dn.DNASeqList(seqarr=arr)
-
-                if sequences.numseqs == 0:
-                    logger.debug(f'found no sequences Hamming distance {sampled_distance} '
-                                 f'from {previous_sequence}; sampling a new distance')
-                    available_distances_list.remove(sampled_distance)
-                else:
-                    sequence = rng.choice(sequences)
+        # pick a distance at random, then re-pick if no sequences are at that distance
+        available_distances_list = list(range(1, len(previous_sequence) + 1))
+        num_to_generate = 1000
+        seqs: Optional[dn.DNASeqList] = None
+        sequence: Optional[str] = None
 
-        else:
-            # first determine all sequences at all distances, then pick distance at random from list of
-            # distances that have at least one sequence at that distance (takes more time usually)
-            neighbors = self.find_hamming_distances(previous_sequence)
+        while sequence is None:
+            if len(available_distances_list) == 0:
+                raise ValueError('out of Hamming distances to try, quitting')
 
-            # Some distances may not exist, so scale the probabilities of the remaining so they sum to one
-            available_distances_arr = np.array(list(neighbors.keys()))
+            # sample a Hamming distance that we haven't tried yet
+            available_distances_arr = np.array(available_distances_list)
             existing_hamming_probabilities = hamming_probabilities[available_distances_arr - 1]
             prob_sum = existing_hamming_probabilities.sum()
             existing_hamming_probabilities /= prob_sum
+            sampled_distance: int = rng.choice(available_distances_arr, p=existing_hamming_probabilities)
+
+            while seqs is None or len(seqs) == 0:
+                bases = self._bases_to_use()
+                length = self.length
+
+                num_ways_to_choose_subsequence_indices = dn.comb(length, sampled_distance)
+                num_different_bases = len(bases) - 1
+                num_subsequences = num_different_bases ** sampled_distance
+                num_sequences_at_sampled_distance = num_ways_to_choose_subsequence_indices * num_subsequences
+
+                if num_to_generate > num_sequences_at_sampled_distance:
+                    num_to_generate = num_sequences_at_sampled_distance
+
+                if num_to_generate >= num_sequences_at_sampled_distance // 4:
+                    # if we want sufficiently many random sequences, just generate all and sample a subset
+                    all_seqs = dn.DNASeqList(
+                        hamming_distance_from_sequence=(sampled_distance, previous_sequence), alphabet=bases,
+                        shuffle=True, rng=rng)
+                    generated_all_seqs = True
+                else:
+                    # otherwise sample num_to_generate with replacement
+                    all_seqs = dn.DNASeqList(
+                        hamming_distance_from_sequence=(sampled_distance, previous_sequence), alphabet=bases,
+                        shuffle=True, num_random_seqs=num_to_generate, rng=rng)
+                    generated_all_seqs = False
+
+                seqs = self._filter_numpy_constraints(all_seqs)
+                self._log_numpy_generation(length, num_to_generate, len(seqs))
+                self._filter_sequence_constraints(seqs)
+
+                if num_to_generate >= 10 ** 9 and len(seqs) == 0:
+                    logger.info("We've generated over 1 million random DNA sequences and not "
+                                "found one that passed your NumpyConstraints and "
+                                f"SequenceConstraints at Hamming distance {sampled_distance} from "
+                                f"the previous sequence {previous_sequence}. Trying another distance")
+                    available_distances_list.remove(sampled_distance)
+                num_to_generate *= 2
 
-            # randomly chooses a Hamming distance to change previous_sequence by
-            sampled_distance = rng.choice(available_distances_arr, p=existing_hamming_probabilities)
-            assert len(neighbors[sampled_distance]) > 0
-            sequence = rng.choice(neighbors[sampled_distance])
+                if len(seqs) == 0 and generated_all_seqs:
+                    # found no sequences passing constraints at distance `sampled_distance`,
+                    # need to try a new Hamming distance
+                    break
 
-        return sequence
+            if len(seqs) > 0:
+                sequence = seqs[0]
 
-    def _get_next_sequence_satisfying_numpy_and_sequence_constraints(self, rng: np.random.Generator) -> str:
-        if self.sequences.numseqs == 0 or self.num_sampled >= len(self.sequences):
-            self._reset_precomputed_sequences(rng)
-        idx = int(rng.integers(self.sequences.numseqs))
-        sequence = self.sequences[idx]
-        self.num_sampled += 1
         return sequence
 
-    def _generate_sequences_satisfying_numpy_constraints(self, rng: np.random.Generator) \
-            -> Tuple[dn.DNASeqList, bool]:
+        # For efficiency we inline the logic of self.sequences_at_hamming_distance()
+        # The next line is the most expensive part of the calculation, so we only do it once,
+        # rather than calling self.sequences_at_hamming_distance() repeatedly, which would recalculate
+        # the next line for each sampled distance even though it doesn't change.
+        # computed_distances = np.sum(
+        #     np.bitwise_xor(self.sequences.seqarr, previous_sequence_1d_array) != 0, axis=1)
+
+        # while sequence is None:
+        #     available_distances_arr = np.array(available_distances_list)
+        #     existing_hamming_probabilities = hamming_probabilities[available_distances_arr - 1]
+        #     prob_sum = existing_hamming_probabilities.sum()
+        #     existing_hamming_probabilities /= prob_sum
+        #     sampled_distance = rng.choice(available_distances_arr, p=existing_hamming_probabilities)
+        #
+        #     indices_at_distance = computed_distances == sampled_distance
+        #     arr = self.sequences.seqarr[indices_at_distance]
+        #     sequences = dn.DNASeqList(seqarr=arr)
+        #
+        #     if sequences.numseqs == 0:
+        #         logger.debug(f'found no sequences Hamming distance {sampled_distance} '
+        #                      f'from {previous_sequence}; sampling a new distance')
+        #         available_distances_list.remove(sampled_distance)
+        #     else:
+        #         sequence = rng.choice(sequences)
+        #
+        # return sequence
+
+    def _get_next_sequence_satisfying_numpy_and_sequence_constraints(self, rng: np.random.Generator) -> str:
+        num_to_generate = 1000
+        seqs: Optional[dn.DNASeqList] = None
+
+        while seqs is None or seqs.numseqs == 0:
+            seqs = self._generate_sequences_satisfying_numpy_constraints(rng, num_to_generate)
+            self._filter_sequence_constraints(seqs)
+            if num_to_generate > 10 ** 10:
+                raise NotImplementedError("We've generated over 10 billion random DNA sequences and not "
+                                          "found one that passed your NumpyConstraints and "
+                                          "SequenceConstraints. Try relaxing the constraints so that it's "
+                                          "more likely a random sequence satisfies the constraints.")
+            num_to_generate *= 2
+
+        return seqs[0]
+
+    def _generate_sequences_satisfying_numpy_constraints(self, rng: np.random.Generator,
+                                                         num_to_generate: int) -> dn.DNASeqList:
         bases = self._bases_to_use()
         length = self.length
-        # For lengths at most _length_threshold_numpy, we generate all DNA sequences in advance.
-        # Above this value, a random subset of DNA sequences will be generated.
-        _length_threshold_numpy = math.floor(math.log(self.num_sequences_to_generate, 4))
-        # _length_threshold_numpy = 10
-        use_random_subset = length > _length_threshold_numpy
-        if not use_random_subset:
-            seqs = dn.DNASeqList(length=length, alphabet=bases, shuffle=True, rng=rng)
-            num_starting_seqs = seqs.numseqs
-        else:
-            num_starting_seqs = self.num_sequences_to_generate
-            seqs = dn.DNASeqList(length=length, alphabet=bases, shuffle=True,
-                                 num_random_seqs=num_starting_seqs, rng=rng)
+        _length_threshold_numpy = math.floor(math.log(num_to_generate, 4))
+        seqs = dn.DNASeqList(length=length, alphabet=bases, shuffle=True,
+                             num_random_seqs=num_to_generate, rng=rng)
         seqs_satisfying_numpy_constraints = self._filter_numpy_constraints(seqs)
-        while seqs_satisfying_numpy_constraints.numseqs == 0 and \
-                self.num_sequences_to_generate < self.num_sequences_to_generate_upper_limit:
-            # 4 ** 13 sequences or more takes over a minute to generate and an excessive amount of memory
-            print(f'No valid sequences found in {self.num_sequences_to_generate} sequences examined. '
-                  f'Increasing number of sequences generated by a factor of 10.')
-            self.num_sequences_to_generate *= 10
-            num_starting_seqs = self.num_sequences_to_generate
-            seqs = dn.DNASeqList(length=length, alphabet=bases, shuffle=True,
-                                 num_random_seqs=num_starting_seqs, rng=rng)
-            seqs_satisfying_numpy_constraints = self._filter_numpy_constraints(seqs)
-        else:
-            if self.num_sequences_to_generate >= self.num_sequences_to_generate_upper_limit:
-                raise ValueError("Too many sequences need to be generated to find a satisfactory sequence. "
-                                 "Please adjust constraints. ")
-        num_decimals = len(str(self.num_sequences_to_generate))
-
-        logger.info(f'generated {num_starting_seqs:{num_decimals}} sequences '
-                    f'of length {length:2}, '
-                    f'of which {len(seqs_satisfying_numpy_constraints):{num_decimals}} '
-                    f'passed the numpy sequence constraints'
-                    f'{" (generated at random)" if use_random_subset else ""}')
-        return seqs_satisfying_numpy_constraints, use_random_subset
+        self._log_numpy_generation(length, num_to_generate, len(seqs_satisfying_numpy_constraints))
+        return seqs_satisfying_numpy_constraints
+
+    def _log_numpy_generation(self, length: int, num_to_generate: int, num_passed: int):
+        if log_numpy_generation:
+            num_decimals = len(str(num_to_generate))
+            logger.debug(f'generated {num_to_generate:{num_decimals}} sequences '
+                         f'of length {length:2}, '
+                         f'of which {num_passed:{num_decimals}} '
+                         f'passed the numpy sequence constraints')
 
     def _bases_to_use(self) -> Collection[str]:
         # checks explicitly for NumpyRestrictBasesConstraint
@@ -969,11 +905,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 +960,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 +1032,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 +1053,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 +1192,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 +1202,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 +1222,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 +1274,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 +1470,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 +1580,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 +1609,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 +1775,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.
@@ -1864,7 +1869,7 @@ def to_json_serializable(self, suppress_indent: bool = True) -> Union[NoIndent,
     def from_json_serializable(json_map: Dict[str, Any],
                                domain_with_name: Dict[str, Domain[DomainLabel]],
                                label_decoder: Callable[[Any], StrandLabel] = (lambda label: label),
-                               ) -> 'Strand[StrandLabel, DomainLabel]':
+                               ) -> Strand[StrandLabel, DomainLabel]:
         """
         :return:
             :any:`Strand` represented by dict `json_map`, assuming it was created by
@@ -1945,7 +1950,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
@@ -2589,15 +2594,12 @@ def to_json(self) -> str:
 
     @staticmethod
     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.
-        :param pool_with_name:
-            If specified should map a name to the :any:`DomainPool` with that name.
         :param domain_label_decoder:
             Function that transforms JSON representation of :py:data:`Domain.label` into the proper type.
         :param strand_label_decoder:
@@ -2608,8 +2610,7 @@ def from_json(json_str: str,
         """
         json_map = json.loads(json_str)
         design: Design[StrandLabel, DomainLabel] = Design.from_json_serializable(
-            json_map, pool_with_name=pool_with_name,
-            domain_label_decoder=domain_label_decoder, strand_label_decoder=strand_label_decoder)
+            json_map, domain_label_decoder=domain_label_decoder, strand_label_decoder=strand_label_decoder)
         return design
 
     def to_json_serializable(self, suppress_indent: bool = True) -> Dict[str, Any]:
@@ -2623,12 +2624,11 @@ def to_json_serializable(self, suppress_indent: bool = True) -> Dict[str, Any]:
         return {
             strands_key: [strand.to_json_serializable(suppress_indent) for strand in self.strands],
             domains_key: [domain.to_json_serializable(suppress_indent) for domain in self.domains],
-            domain_pools_num_sampled_key: {pool.name: pool.num_sampled for pool in self.domain_pools()}
+            domain_pools_key: [pool.to_json_serializable(suppress_indent) for pool in self.domain_pools()]
         }
 
     @staticmethod
     def from_json_serializable(json_map: Dict[str, Any],
-                               pool_with_name: Optional[Dict[str, DomainPool]] = None,
                                domain_label_decoder: Callable[[Any], DomainLabel] = lambda label: label,
                                strand_label_decoder: Callable[[Any], StrandLabel] = lambda label: label,
                                ) -> 'Design[StrandLabel, DomainLabel]':
@@ -2636,10 +2636,6 @@ def from_json_serializable(json_map: Dict[str, Any],
         :param json_map:
             JSON serializable object encoding this :any:`Design`, as returned by
             :py:meth:`Design.to_json_serializable`.
-        :param pool_with_name:
-            dict mapping name to :any:`DomainPool` with that name; required to rehydrate :any:`Domain`'s.
-            If None, then a DomainPool with no constraints is created with the name and domain length
-            found in the JSON.
         :param domain_label_decoder:
             Function that transforms JSON representation of :py:data:`Domain.label` into the proper type.
         :param strand_label_decoder:
@@ -2648,6 +2644,10 @@ def from_json_serializable(json_map: Dict[str, Any],
             :any:`Design` represented by dict `json_map`, assuming it was created by
             :py:meth:`Design.to_json_serializable`. No constraints are populated.
         """
+        pools_json = mandatory_field(Design, json_map, domain_pools_key)
+        pools: List[DomainPool] = [DomainPool.from_json_serializable(pool_json) for pool_json in pools_json]
+        pool_with_name: Dict[str, DomainPool] = {pool.name: pool for pool in pools}
+
         domains_json = mandatory_field(Design, json_map, domains_key)
         domains: List[Domain] = [
             Domain.from_json_serializable(domain_json, pool_with_name=pool_with_name,
@@ -2720,345 +2720,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 +2917,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 +3088,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 +3124,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 +3144,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
 
-    def __init__(self, constraint: Constraint, domains: Iterable[Domain], score: float):
+    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 +3164,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 +3193,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.
 
-    # Set of all :any:`Violations` that are associated to non-fixed :any:`constraint.Domain`'s.
+    num_checked: Dict[Constraint, int]
+
+    # number of instances of each :any:`Constraint` that were checked
+
+    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 +3231,13 @@ 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 +3253,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 +3281,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 +3325,56 @@ 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
+    def has_nonfixed_violations(self) -> bool:
+        # :return: whether there are any nonfixed Violations in this ViolationSet
+        return self.num_violations_nonfixed() > 0
+
+
+@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 +3399,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 +3407,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 +3441,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
 
-_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.'
+    @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 = "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."""
+    """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.
-    """
+    Specify :data:`Constraint._evaluate` in the constructor."""
 
-    summary: Callable[[Domain, Domain],
-                      str] = lambda _, __: _no_summary_string
-
-    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.
-    """
+    """Constraint that applies to a pair of :any:`Strand`'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.
-    """
+    Specify :data:`Constraint._evaluate` in the constructor."""
 
-    summary: Callable[[Iterable[Tuple[Domain, Domain]], bool],
-                      ConstraintReport] = lambda _: _no_summary_string
-
-    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
+    Specify :data:`Constraint._evaluate_bulk` in the constructor.
 
-    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
-
-    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 +3627,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 +3654,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`.
+    """
+
+    def __post_init__(self) -> None:
+        if self._evaluate_bulk is None:
+            raise ValueError('_evaluate_bulk must be specified for a StrandsConstraint')
+        super().__post_init__()
 
-    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.
+    @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.
+
+    Specify :data:`Constraint._evaluate_bulk` in the constructor.
+    """
+
+    def __post_init__(self) -> None:
+        if self._evaluate_bulk is None:
+            raise ValueError('_evaluate_bulk must be specified for a DomainPairsConstraint')
+        super().__post_init__()
+
+    @staticmethod
+    def part_name() -> str:
+        return 'domain pair'
+
 
-    evaluate: Callable[[Iterable[Strand]],
-                       List[Tuple[OrderedSet[Domain], float]]] = lambda _: []
+@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.
 
-    summary: Callable[[Iterable[Strand], bool],
-                      ConstraintReport] = lambda _: _no_summary_string
+    Specify :data:`Constraint._evaluate_bulk` in the constructor.
+    """
 
-    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
+    def __post_init__(self) -> None:
+        if self._evaluate_bulk is None:
+            raise ValueError('_evaluate_bulk must be specified for a StrandPairsConstraint')
+        super().__post_init__()
 
-    def generate_summary(self, strands: Iterable[Strand], report_only_violations: bool) -> ConstraintReport:
-        return (self.summary)(strands, report_only_violations)  # noqa
+    @staticmethod
+    def part_name() -> str:
+        return 'strand pair'
 
 
 @dataclass(frozen=True, eq=False)  # type: ignore
@@ -3929,33 +3708,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 +3820,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 +3920,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 +3933,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 +3958,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 +3982,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 +4016,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 +4052,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 +4069,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 +4083,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 +4105,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 +4121,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 +4148,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 +4221,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 +4293,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 +4313,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 +4327,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 +4366,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 +4394,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 +4414,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 +4428,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 +4467,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 +4487,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 +4510,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 +4559,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 +5994,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 +6145,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 +6162,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 +6263,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 +6301,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 +6338,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 +6366,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 +6487,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 +6508,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/np.py b/dsd/np.py
index ca80c3ec..52ee772d 100644
--- a/dsd/np.py
+++ b/dsd/np.py
@@ -36,9 +36,11 @@ def idx2seq(idx: int, length: int) -> str:
     return ''.join(seq)
 
 
-def seq2arr(seq: str) -> np.ndarray:
+def seq2arr(seq: str, base2bits_local: Optional[Dict[str, int]] = None) -> np.ndarray:
     """Convert seq (string with DNA alphabet) to numpy array with integers 0,1,2,3."""
-    return np.array([base2bits[base] for base in seq], dtype=np.ubyte)
+    if base2bits_local is None:
+        base2bits_local = base2bits
+    return np.array([base2bits_local[base] for base in seq], dtype=np.ubyte)
 
 
 def seqs2arr(seqs: Sequence[str]) -> np.ndarray:
@@ -61,6 +63,19 @@ def arr2seq(arr: np.ndarray) -> str:
     return ''.join(bases_ch)
 
 
+def make_array_with_all_sequences(length: int, digits: Sequence[int]) -> np.ndarray:
+    num_digits = len(digits)
+    num_seqs = num_digits ** length
+    powers_num_digits = [num_digits ** k for k in range(length)]
+    digits = np.array(digits, dtype=np.ubyte)
+
+    arr = np.zeros((num_seqs, length), dtype=np.ubyte)
+    for i, j, c in zip(reversed(powers_num_digits), powers_num_digits, list(range(length))):
+        arr[:, c] = np.tile(np.repeat(digits, i), j)
+
+    return arr
+
+
 def make_array_with_all_dna_seqs(length: int, bases: Collection[str] = ('A', 'C', 'G', 'T')) -> np.ndarray:
     """Return 2D numpy array with all DNA sequences of given length in
     lexicographic order. Bases contains bases to be used: ('A','C','G','T') by
@@ -70,50 +85,44 @@ def make_array_with_all_dna_seqs(length: int, bases: Collection[str] = ('A', 'C'
     is a 2D array, where each row represents a DNA sequence, and that row
     has one byte per base."""
 
-    if not set(bases) <= {'A', 'C', 'G', 'T'}:
-        raise ValueError(f"bases must be a subset of {'A', 'C', 'G', 'T'}; cannot be {bases}")
     if len(bases) == 0:
         raise ValueError('bases cannot be empty')
+    if not set(bases) <= {'A', 'C', 'G', 'T'}:
+        raise ValueError(f"bases must be a subset of {'A', 'C', 'G', 'T'}; cannot be {bases}")
 
-    num_bases = len(bases)
-    num_seqs = num_bases ** length
-
-    #     shift = np.arange(2*(length-1), -1, -2)
-    #     nums = np.repeat(np.arange(numseqs), length)
-    #     nums2D = nums.reshape([numseqs, length])
-    #     shifts = np.tile(0b11 << shift, numseqs)
-    #     shifts2D = shifts.reshape([numseqs, length])
-    #     arr = (shifts2D & nums2D) >> shift
-
-    # the former code took up too much memory (using int32 or int64)
-    # the following code makes sure it's 1 byte per base
-    powers_num_bases = [num_bases ** k for k in range(length)]
-    #         bases = np.array([_base2bits['A'], _base2bits['C'], _base2bits['G'], _base2bits['T']], dtype=np.ubyte)
     base_bits = [base2bits[base] for base in bases]
-    bases = np.array(base_bits, dtype=np.ubyte)
-
-    list_of_arrays = False
-    if list_of_arrays:
-        # this one seems to be faster but takes more memory, probably because just before the last command
-        # there are two copies of the array in memory at once
-        columns = []
-        for i, j, c in zip(reversed(powers_num_bases), powers_num_bases, list(range(length))):
-            columns.append(np.tile(np.repeat(bases, i), j))
-        arr = np.vstack(columns).transpose()
-    else:
-        # this seems to be slightly slower but takes less memory, since it
-        # allocates only the final array, plus one extra column of that
-        # array at a time
-        arr = np.empty((num_seqs, length), dtype=np.ubyte)
-        for i, j, c in zip(reversed(powers_num_bases), powers_num_bases, list(range(length))):
-            arr[:, c] = np.tile(np.repeat(bases, i), j)
+    digits = np.array(base_bits, dtype=np.ubyte)
 
-    return arr
+    return make_array_with_all_sequences(length, digits)
+    # num_bases = len(bases)
+    # num_seqs = num_bases ** length
+    #
+    # # the former code took up too much memory (using int32 or int64)
+    # # the following code makes sure it's 1 byte per base
+    # powers_num_bases = [num_bases ** k for k in range(length)]
+    #
+    # list_of_arrays = False
+    # if list_of_arrays:
+    #     # this one seems to be faster but takes more memory, probably because just before the last command
+    #     # there are two copies of the array in memory at once
+    #     columns = []
+    #     for i, j, c in zip(reversed(powers_num_bases), powers_num_bases, list(range(length))):
+    #         columns.append(np.tile(np.repeat(bases, i), j))
+    #     arr = np.vstack(columns).transpose()
+    # else:
+    #     # this seems to be slightly slower but takes less memory, since it
+    #     # allocates only the final array, plus one extra column of that
+    #     # array at a time
+    #     arr = np.empty((num_seqs, length), dtype=np.ubyte)
+    #     for i, j, c in zip(reversed(powers_num_bases), powers_num_bases, list(range(length))):
+    #         arr[:, c] = np.tile(np.repeat(bases, i), j)
+    #
+    # return arr
 
 
-def make_array_with_random_subset_of_dna_seqs(length: int, num_seqs: int,
-                                              rng: np.random.Generator,
-                                              bases: Collection[str] = ('A', 'C', 'G', 'T')) -> np.ndarray:
+def make_array_with_random_subset_of_dna_seqs(
+        length: int, num_seqs: int, rng: np.random.Generator = default_rng,
+        bases: Collection[str] = ('A', 'C', 'G', 'T')) -> np.ndarray:
     """
     Return 2D numpy array with random subset of size `num_seqs` of DNA sequences of given length.
     Bases contains bases to be used: ('A','C','G','T') by default, but can be set to a subset of these.
@@ -121,12 +130,25 @@ def make_array_with_random_subset_of_dna_seqs(length: int, num_seqs: int,
     Uses the encoding described in the documentation for DNASeqList. The result is a 2D array,
     where each row represents a DNA sequence, and that row has one byte per base.
 
-    :param length: length of each row
-    :param num_seqs: number of rows
-    :param bases: DNA bases to use
-    :param rng: numpy random number generator (type returned by numpy.random.default_rng())
-    :return: 2D numpy array with random subset of size `num_seqs` of DNA sequences of given length
+    Sequences returned will be unique (i.e., sampled without replacement) and in a random order
+
+    :param length:
+        length of each row
+    :param num_seqs:
+        number of rows
+    :param bases:
+        DNA bases to use
+    :param rng:
+        numpy random number generator (type returned by numpy.random.default_rng())
+    :return:
+        2D numpy array with random subset of size `num_seqs` of DNA sequences of given length
     """
+    if length < 0:
+        raise ValueError(f'length = {num_seqs} must be nonnegative')
+    elif length == 0:
+        return np.array([[]], dtype=np.ubyte)
+    if num_seqs <= 0:
+        raise ValueError(f'num_seqs = {num_seqs} must be positive')
     if not set(bases) <= {'A', 'C', 'G', 'T'}:
         raise ValueError(f"bases must be a subset of {'A', 'C', 'G', 'T'}; cannot be {bases}")
     if len(bases) == 0:
@@ -134,12 +156,250 @@ def make_array_with_random_subset_of_dna_seqs(length: int, num_seqs: int,
     elif len(bases) == 1:
         raise ValueError('bases must have at least two elements')
 
+    max_possible = len(bases) ** length
+    if num_seqs > max_possible:
+        raise ValueError(f'num_seqs = {num_seqs} is greater than the total number {max_possible} '
+                         f'of sequences of length {length} using alphabet {bases}, so we cannot guarantee '
+                         f'that many unique sequences. Please set num_seqs <= {max_possible}.')
+
+    # If we want sufficiently many sequences, then it's simpler to just generate all sequences
+    # of that length and choose a random subset of size num_seqs.
+    if num_seqs >= max_possible / 4:
+        all_seqs = make_array_with_all_dna_seqs(length=length, bases=bases)
+        # https://stackoverflow.com/a/27815343/5339430
+        idxs = rng.choice(all_seqs.shape[0], num_seqs, replace=False)
+        sampled_seqs = all_seqs[idxs]
+        return sampled_seqs
+
+    # This comment justifies why we sample 2*num_seqs sequences randomly (with replacement) in order
+    # to get our goal of at least num_seqs *unique* sequences. Define
+    #
+    #   m = number of sequences we sample with replacement
+    #   n = len(bases)^length = total number of sequences of that length
+    #   k = num_seqs = desired number of unique sequences
+    #
+    # If we sample m sequences with replacement, out of n total, this is tossing m balls into n bins.
+    # We want m sufficiently large that at least k bins are non-empty.
+    # We use the Poisson approximation to balls-in-bins
+    #   (Probability and Computing, Mitzenbacher and Upfal, 2nd edition, Section 5.4),
+    # where each bin is modeled as getting P(m/n) balls, where P(m/n) is a Poisson r.v. with rate m/n.
+    # The sum of Poisson r.v.'s is Poisson, so the number of non-empty bins is Poisson with rate m.
+    # We use this bound:
+    #   https://doi.org/10.1109/TIT.2006.890791, Appendix II
+    # for Poisson r.v.'s P(y)  with rate m to be less than a constant k, if k < m:
+    #
+    #   Pr[P(m) <= k] <= e^{-m} * (e*m / k)^k
+    #
+    # They state only the bound for the upper tail, but the proof works for the lower tail as well.
+    #
+    # Setting m = c*k for c > 1, we have
+    #   Pr[P(m) <= k]
+    #     <= e^{-c*k} * (e*c*k / k)^k
+    #      = (e^{-c})^k * (e*c)^k
+    #      = (e^{1-c}*c)^k
+    # Setting c = 2 gives e^{1-c}*c ~ 0.736, so  Pr[P(m) <= k] <= 0.75^k.
+
     base_bits = np.array([base2bits[base] for base in bases], dtype=np.ubyte)
+    num_seqs_to_sample = 2 * num_seqs
+    unique_sorted_arr = None
+
+    # odds are low to have a collision, so for simplicity we just repeat the whole process if needed
+    while unique_sorted_arr is None or len(unique_sorted_arr) < num_seqs:
+        arr = rng.choice(a=base_bits, size=(num_seqs_to_sample, length))
+        unique_sorted_arr = np.unique(arr, axis=0)
+        if len(unique_sorted_arr) < num_seqs:
+            print(f'WARNING: did not find {num_seqs} unique sequences. If you are seeing this warning '
+                  f'repeatedly, check the parameters to make_array_with_random_subset_of_dna_seqs.')
+
+    # we probably have too many, so pick a random subset of sequences to return
+    idxs = rng.choice(unique_sorted_arr.shape[0], num_seqs, replace=False)
+    sampled_seqs = unique_sorted_arr[idxs]
+    return sampled_seqs
+
+
+# https://stackoverflow.com/questions/4941753/is-there-a-math-ncr-function-in-python
+# In Python 3.8 there's math.comb, but this is about 3x faster somehow.
+import operator as op
+from functools import reduce
+
+
+def comb(n: int, k: int) -> int:
+    # n choose k = n! / (k! * (n-k)!)
+    k = min(k, n - k)
+    numer = reduce(op.mul, range(n, n - k, -1), 1)
+    denom = reduce(op.mul, range(1, k + 1), 1)
+    return numer // denom
+
+
+
+def make_array_with_all_dna_seqs_hamming_distance(
+        dist: int, seq: str, bases: Collection[str] = ('A', 'C', 'G', 'T')) -> np.ndarray:
+    """
+    Return 2D numpy array with all DNA sequences of given length in lexicographic order. Bases contains
+    bases to be used: ('A','C','G','T') by default, but can be set to a subset of these.
+
+    Uses the encoding described in the documentation for DNASeqList. The result is a 2D array, where each
+    row represents a DNA sequence, and that row has one byte per base.
+    """
+    length = len(seq)
+    assert 1 <= dist <= length
+
+    num_bases = len(bases)
+
+    if num_bases == 0:
+        raise ValueError('bases cannot be empty')
+    if not set(bases) <= {'A', 'C', 'G', 'T'}:
+        raise ValueError(f"bases must be a subset of {'A', 'C', 'G', 'T'}; cannot be {bases}")
+
+    num_ways_to_choose_subsequence_indices = comb(length, dist)
+    num_different_bases = len(bases) - 1
+    num_subsequences = num_different_bases ** dist
+    num_seqs = num_ways_to_choose_subsequence_indices * num_subsequences
+
+    # for simplicity of modular arithmetic, we use integers 0,...,len(bases)-1 to represent the bases,
+    # then map these back to the correct subset of 0,1,2,3 when we are done
+    offsets = range(1, num_different_bases + 1)
+    subseq_offsets = make_array_with_all_sequences(length=dist, digits=offsets)
+    assert len(subseq_offsets) == num_subsequences
+    subseq_offsets_repeats = \
+        np.tile(subseq_offsets.flatten(), num_ways_to_choose_subsequence_indices).reshape(num_seqs, dist)
+
+    # all (length choose dist) indices where we could change the bases
+    idxs = combnr_idxs(length, dist)
+    assert len(idxs) == num_ways_to_choose_subsequence_indices
+    idxs_repeat = np.tile(idxs, num_subsequences).reshape(num_seqs, length)
+    assert len(idxs_repeat) == num_seqs
+
+    # map subset of bases used to *prefix* of 0,1,2,3
+    base2bits_local = {base: digit for base, digit in zip(bases, range(4))}
+    seq_as_arr = seq2arr(seq, base2bits_local=base2bits_local)
+    new_arr = np.tile(seq_as_arr, num_seqs).reshape(num_seqs, length)
+
+    new_arr[idxs_repeat] += subseq_offsets_repeats.flatten()
+    new_arr %= num_bases
+
+    # now map back to correct subset of 0,1,2,3 to represent bases
+    for base, digit in zip(['A', 'C', 'G', 'T'], range(4)):
+        if base not in bases:
+            idxs_to_inc = new_arr >= digit
+            new_arr[idxs_to_inc] += 1
+
+    return new_arr
+
+
+def combnr_idxs(length: int, number: int) -> np.ndarray:
+    # Gives 2D Boolean numpy array, with `length` columns and (`length` choose `number`) rows,
+    # representing all ways to set exactly `number` elements of the row True and the others to False.
+    # Useful for indexing into a same-shape numpy array, changing exactly `number` elements in each row.
+    #
+    # :param length:
+    #     number of columns
+    # :param number:
+    #     number of True values in each row
+    # :return:
+    #     numpy array, with `length` columns and (`length` choose `number`) rows,
+    #     representing all ways to set exactly `number` elements of the row True and the others to False.
+    x = np.array(np.meshgrid(*([np.arange(0, length)] * number))).T.reshape(-1, number)
+    z = np.sum(np.identity(length)[x], 1, dtype=bool).astype(int)
+    return np.unique(z[np.sum(z, axis=1) == number], axis=0).astype(bool)
+
+
+def make_array_with_random_subset_of_dna_seqs_hamming_distance(
+        num_seqs: int, dist: int, seq: str, rng: np.random.Generator = default_rng,
+        bases: Collection[str] = ('A', 'C', 'G', 'T')) -> np.ndarray:
+    """
+    Return 2D numpy array with random subset of size `num_seqs` of DNA sequences of given length.
+    Bases contains bases to be used: ('A','C','G','T') by default, but can be set to a subset of these.
+
+    Uses the encoding described in the documentation for DNASeqList. The result is a 2D array,
+    where each row represents a DNA sequence, and that row has one byte per base.
+
+    Sampled *with* replacement, so the same row may appear twice in the returned array
+
+    :param num_seqs:
+        number of sequences to generate
+    :param dist:
+        Hamming distance to be from `seq`
+    :param seq:
+        sequence to generate other sequences close to
+    :param bases:
+        DNA bases to use
+    :param rng:
+        numpy random number generator (type returned by numpy.random.default_rng())
+    :return:
+        2D numpy array with random subset of size `num_seqs` of DNA sequences of given length
+    """
+    if not set(bases) <= {'A', 'C', 'G', 'T'}:
+        raise ValueError(f"bases must be a subset of {'A', 'C', 'G', 'T'}; cannot be {bases}")
+    if len(bases) == 0:
+        raise ValueError('bases cannot be empty')
+    elif len(bases) == 1:
+        raise ValueError('bases must have at least two elements')
+
+    length = len(seq)
+
+    if dist < 1:
+        raise ValueError(f'dist must be positive, but dist = {dist}')
+    if dist > length:
+        raise ValueError(f'should have dist <= len("{seq}") = {length}, but dist = {dist}')
 
-    arr = rng.choice(a=base_bits, size=(num_seqs, length))
-    unique_sorted_arr = np.unique(arr, axis=0)
+    num_different_bases = len(bases) - 1
 
-    return unique_sorted_arr
+    # map subset of bases used to *prefix* of 0,1,2,3
+    base2bits_local = {base: digit for base, digit in zip(bases, range(4))}
+    seq_as_arr = seq2arr(seq, base2bits_local=base2bits_local)
+    seqs = np.tile(seq_as_arr, num_seqs).reshape(num_seqs, length)
+
+    # for simplicity of modular arithmetic, we use integers 0,...,len(bases)-1 to represent the bases,
+    # then map these back to the correct subset of 0,1,2,3 when we are done
+    subseq_offsets = rng.integers(low=1, high=num_different_bases + 1, size=(num_seqs, dist))
+    assert len(subseq_offsets) == num_seqs
+
+    # all (length choose dist) indices where we could change the bases
+    idxs = random_choice_noreplace(np.arange(length), dist, num_seqs, rng)
+    assert len(idxs) == num_seqs
+    changes = rng.integers(1, num_different_bases + 1, size=idxs.shape, dtype=np.uint8)
+    seqs[np.arange(num_seqs)[:, None], idxs] += changes
+    seqs = np.mod(seqs, num_different_bases + 1)
+
+    # now map back to correct subset of 0,1,2,3 to represent bases
+    for base, digit in zip(['A', 'C', 'G', 'T'], range(4)):
+        if base not in bases:
+            idxs_to_inc = seqs >= digit
+            seqs[idxs_to_inc] += 1
+
+    # use the next two lines to return only unique rows
+    # seqs = np.unique(seqs, axis=0)
+    # rng.shuffle(seqs)
+
+    return seqs
+
+
+# def random_hamming(sequence: Union[List[int], np.ndarray], distance: int, number: int,
+#                    rng: np.random.Generator) -> np.ndarray:
+#     sequence = np.array(sequence, dtype=np.uint8)
+#     length = len(sequence)
+#     seqrepeats = np.tile(sequence, number).reshape((number, length))
+#     places = random_choice_noreplace(np.arange(length), distance, number, rng)
+#     changes = rng.integers(1, 4, size=places.shape, dtype=np.uint8)
+#     seqrepeats[np.arange(number)[:, None], places] += changes
+#     seqrepeats = np.mod(seqrepeats, 4)
+#     return seqrepeats
+
+
+# https://stackoverflow.com/a/59328647/5339430
+def random_choice_noreplace(l: np.ndarray, n_sample: int, num_draw: int,
+                            rng: np.random.Generator) -> np.ndarray:
+    '''
+    l: 1-D array or list
+    n_sample: sample size for each draw
+    num_draw: number of draws
+
+    Intuition: Randomly generate numbers, get the index of the smallest n_sample number for each row.
+    '''
+    l = np.array(l)
+    random_array_floats = rng.random((num_draw, len(l)))
+    return l[np.argpartition(random_array_floats, n_sample - 1, axis=-1)[:, :n_sample]]
 
 
 # @lru_cache(maxsize=10000000)
@@ -413,7 +673,8 @@ def __init__(self,
                  seqs: Optional[Sequence[str]] = None,
                  seqarr: np.ndarray = None,
                  filename: Optional[str] = None,
-                 rng: np.random.Generator = default_rng):
+                 rng: np.random.Generator = default_rng,
+                 hamming_distance_from_sequence: Optional[Tuple[int, str]] = None):
         """
         Creates a set of DNA sequences, all of the same length.
 
@@ -431,25 +692,42 @@ def __init__(self,
 
         - `filename`
 
-        :param length: length of sequences; `num_seqs` and `alphabet` can also be specified along with it
-        :param num_random_seqs: number of sequences to generate; if not specified, then all sequences
-                                of length `length` using bases from `alphabet` are generated
-        :param shuffle: whether to shuffle sequences
-        :param alphabet: alphabet; must be a subset of {'A', 'C', 'G', 'T'}
-        :param seqs: sequence (e.g., list or tuple) of strings, all of the same length
-        :param seqarr: 2D NumPy array, with axis 0 moving between sequences,
-                       and axis 1 moving between consecutive DNA bases in a sequence
-        :param filename: name of file containing a :any:`DNASeqList`
-                         as written by :py:meth:`DNASeqList.write_to_file`
-        :param rng: numpy random number generator (type returned by numpy.random.default_rng())
+        - `hamming_distance_from_sequence`
+
+        :param length:
+            length of sequences; `num_seqs` and `alphabet` can also be specified along with it
+        :param hamming_distance_from_sequence:
+            if specified and equal to `(dist, seq)`,
+            then only sequences at Hamming distance `dist` from `seq` will be generated.
+            Raises error if `length`, `seqs`, `seqarr`, or `filename` is specified.
+        :param num_random_seqs:
+            number of sequences to generate; if not specified, then all sequences
+            of length `length` using bases from `alphabet` are generated.
+            Sequences are sampled *with* replacement, so the same sequence may appear twice.
+        :param shuffle:
+            whether to shuffle sequences
+        :param alphabet:
+            a subset of {'A', 'C', 'G', 'T'}
+        :param seqs:
+            sequence (e.g., list or tuple) of strings, all of the same length
+        :param seqarr:
+            2D NumPy array, with axis 0 moving between sequences,
+            and axis 1 moving between consecutive DNA bases in a sequence
+        :param filename:
+            name of file containing a :any:`DNASeqList`
+            as written by :py:meth:`DNASeqList.write_to_file`
+        :param rng:
+            numpy random number generator (type returned by numpy.random.default_rng())
         """
-        for v1, v2 in it.combinations([length, seqs, seqarr, filename], 2):
+        for v1, v2 in it.combinations([length, seqs, seqarr, filename, hamming_distance_from_sequence], 2):
             if v1 is not None and v2 is not None:
-                raise ValueError('exactly one of length, seqs, seqarra, or filename can be non-None')
+                raise ValueError('exactly one of length, seqs, seqarr, filename, or '
+                                 'hamming_distance_from_sequence must be non-None')
         self.rng = rng
         if seqarr is not None:
             self.seqarr = seqarr
             self._update_size()
+
         elif seqs is not None:
             if len(seqs) == 0:
                 raise ValueError('seqs must have positive length')
@@ -459,18 +737,33 @@ def __init__(self,
                     raise ValueError('All sequences in seqs must be equal length')
             self.numseqs = len(seqs)
             self.seqarr = seqs2arr(seqs)
+
         elif filename is not None:
             self._read_from_file(filename)
+
         elif length is not None:
-            self.seqlen = length
             if num_random_seqs is None:
-                self.seqarr = make_array_with_all_dna_seqs(self.seqlen, alphabet)
+                self.seqarr = make_array_with_all_dna_seqs(length=length, bases=alphabet)
             else:
                 self.seqarr = make_array_with_random_subset_of_dna_seqs(
-                    self.seqlen, num_random_seqs, self.rng, alphabet)
-            self.numseqs = len(self.seqarr) if length > 0 else 0
+                    length=length, num_seqs=num_random_seqs, rng=self.rng, bases=alphabet)
+            self.seqlen = length
+            self.numseqs = len(self.seqarr) if self.seqlen > 0 else 1
+
+        elif hamming_distance_from_sequence is not None:
+            dist, seq = hamming_distance_from_sequence
+            if num_random_seqs is None:
+                self.seqarr = make_array_with_all_dna_seqs_hamming_distance(dist=dist, seq=seq,
+                                                                            bases=alphabet)
+            else:
+                self.seqarr = make_array_with_random_subset_of_dna_seqs_hamming_distance(
+                    num_seqs=num_random_seqs, dist=dist, seq=seq, rng=self.rng, bases=alphabet)
+            self.seqlen = len(seq)
+            self.numseqs = len(self.seqarr) if self.seqlen > 0 else 1
+
         else:
-            raise ValueError('at least one of length, seqs, seqarr, or filename must be specified')
+            raise ValueError('at least one of length, seqs, seqarr, filename, or '
+                             'hamming_distance_from_sequence must be specified')
 
         self.shift = np.arange(2 * (self.seqlen - 1), -1, -2)
 
@@ -528,7 +821,7 @@ def __repr__(self) -> str:
         return 'DNASeqSet(seqs={})'.format(str([self[i] for i in range(self.numseqs)]))
 
     def __str__(self) -> str:
-        if self.numseqs <= 64:
+        if self.numseqs <= 256:
             ret = [self.get_seq_str(i) for i in range(self.numseqs)]
             return ','.join(ret)
         else:
diff --git a/dsd/search.py b/dsd/search.py
index 9cfe1c1c..8b0213a8 100644
--- a/dsd/search.py
+++ b/dsd/search.py
@@ -54,18 +54,15 @@
 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
+from typing import List, Tuple, Sequence, FrozenSet, Optional, Dict, Callable, Iterable, Any, \
+    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 +87,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
@@ -109,7 +108,7 @@ def new_process_pool(cpu_count: int):
 
 
 def default_output_directory() -> str:
-    return os.path.join('output', f'{script_name_no_ext()}--{timestamp()}')
+    return os.path.join('out', f'{script_name_no_ext()}--{timestamp()}')
 
 
 def _violations_of_constraints(design: Design,
@@ -158,169 +157,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)
+        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)
 
-        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)
+        # 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
 
-    # 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
-
-    # 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)
-
-        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 +232,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 +254,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 +313,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 +366,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 +385,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
+        raise NotImplementedError('TODO: implement parallelization')
 
-        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
-
-    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
-
-
-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)
+        raise AssertionError()
+    return score, summary
 
-    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}')
+        raise NotImplementedError()
 
-        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
-
-    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
-
-
-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
+T = TypeVar('T')
 
-    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,
@@ -927,53 +501,88 @@ def _sequences_fragile_format_output_to_file(design: Design,
         f'{strand.sequence(delimiter="-")}' for strand in design.strands)
 
 
-def _write_sequences(design: Design, directory_intermediate: str, directory_final: str,
-                     filename_with_iteration_no_ext: str, filename_final_no_ext: str,
+def _write_intermediate_files(*, design: dc.Design, params: SearchParameters, rng: numpy.random.Generator,
+                              num_new_optimal: int, directories: _Directories,
+                              violation_set: dc.ViolationSet) -> None:
+    num_new_optimal_padded = f'{num_new_optimal}' if params.num_digits_update is None \
+        else f'{num_new_optimal:0{params.num_digits_update}d}'
+
+    _write_design_json(design, params=params, directories=directories,
+                       num_new_optimal_padded=num_new_optimal_padded)
+
+    _write_rng_state(rng, params=params, directories=directories,
+                     num_new_optimal_padded=num_new_optimal_padded)
+
+    _write_sequences(design, params=params, directories=directories,
+                     num_new_optimal_padded=num_new_optimal_padded)
+
+    _write_report(design, params=params, directories=directories,
+                  num_new_optimal_padded=num_new_optimal_padded, violation_set=violation_set)
+
+
+def _write_sequences(design: Design, params: SearchParameters, directories: _Directories,
+                     num_new_optimal_padded: str,
                      include_group: bool = True) -> None:
+    directory_intermediate = directories.sequence
+    directory_final = directories.out
+    # filename_with_iteration_no_ext = directories.indexed_sequences_full_filename_noext(num_new_optimal_padded)
+    # filename_final_no_ext = directories.best_sequences_full_filename_noext()
+    filename_with_iteration_no_ext = f'{directories.sequences_filename_no_ext}' + \
+                                     f'-{num_new_optimal_padded}'
+    filename_final_no_ext = f'current-best-{directories.sequences_filename_no_ext}'
     sequences_content = _sequences_fragile_format_output_to_file(design, include_group)
+    if not params.save_sequences_for_all_updates:
+        directory_intermediate = filename_with_iteration_no_ext = None
     _write_text_intermediate_and_final_files(directory_final, directory_intermediate,
                                              filename_final_no_ext, filename_with_iteration_no_ext,
                                              sequences_content, '.txt')
 
 
-def _write_dsd_design_json(design: Design, directories: _Directories, directory_final: str,
-                           filename_final_no_ext: str, num_new_optimal_padded: str) -> None:
-    directory_intermediate = directories.dsd_design
-    filename_with_iteration_no_ext = f'{directories.dsd_design_filename_no_ext}-{num_new_optimal_padded}'
+def _write_design_json(design: Design, params: SearchParameters, directories: _Directories,
+                       num_new_optimal_padded: str) -> None:
+    filename_final_no_ext = f'current-best-{directories.design_filename_no_ext}'
+    directory_intermediate = directories.design
+    filename_with_iteration_no_ext = f'{directories.design_filename_no_ext}-{num_new_optimal_padded}'
     json_str = design.to_json()
-    _write_text_intermediate_and_final_files(directory_final, directory_intermediate,
+    if not params.save_designs_for_all_updates:
+        directory_intermediate = filename_with_iteration_no_ext = None
+    _write_text_intermediate_and_final_files(directories.out, directory_intermediate,
                                              filename_final_no_ext, filename_with_iteration_no_ext,
                                              json_str, '.json')
 
 
-def _write_rng_state(rng: numpy.random.Generator, directories: _Directories,
+def _write_rng_state(rng: numpy.random.Generator, params: SearchParameters, directories: _Directories,
                      num_new_optimal_padded: str) -> None:
+    filename_final_no_ext = f'current-best-{directories.rng_state_filename_no_ext}'
     directory_intermediate = directories.rng_state
     filename_with_iteration_no_ext = f'{directories.rng_state_filename_no_ext}-{num_new_optimal_padded}'
     state = rng.bit_generator.state
     json_str = json.dumps(state, indent=2)
-    _write_text_intermediate_and_final_files(None, directory_intermediate,
-                                             None, filename_with_iteration_no_ext,
+    if not params.save_designs_for_all_updates:
+        directory_intermediate = filename_with_iteration_no_ext = None
+    _write_text_intermediate_and_final_files(directories.out, directory_intermediate,
+                                             filename_final_no_ext, filename_with_iteration_no_ext,
                                              json_str, '.json')
 
 
-def _write_domain_pools(domain_pools: Iterable[dc.DomainPool], directories: _Directories) -> None:
-    directory_intermediate = directories.domain_pools
+def _write_report(design: Design, params: SearchParameters, directories: _Directories,
+                  num_new_optimal_padded: str, violation_set: dc.ViolationSet) -> None:
+    directory_intermediate = directories.report
+    directory_final = directories.out
+    filename_with_iteration_no_ext = f'{directories.report_filename_no_ext}' + \
+                                     f'-{num_new_optimal_padded}'
+    filename_final_no_ext = f'current-best-{directories.report_filename_no_ext}'
 
-    for pool in domain_pools:
-        # to avoid filling drive with large files of sequences, only write domain pool when it resets
-        try:
-            highest_idx = _find_highest_index_in_directory(directories.domain_pools, pool.name, 'json')
-        except ValueError:
-            highest_idx = -1
+    report = f'''\
+Report on constraints
+=====================
+''' + design.summary_of_constraints(params.report_only_violations, violation_set=violation_set)
 
-        if highest_idx < pool.num_times_sequences_reset:
-            filename_no_ext = f'{pool.name}-{pool.num_times_sequences_reset}'
-            domain_pools_json_map = pool.to_json_serializable(include_sequences=True)
-            domain_pools_json_str = json.dumps(domain_pools_json_map, indent=2)
-            _write_text_intermediate_and_final_files(None, directory_intermediate,
-                                                     None, filename_no_ext,
-                                                     domain_pools_json_str, '.json')
+    if not params.save_designs_for_all_updates:
+        directory_intermediate = filename_with_iteration_no_ext = None
+    _write_text_intermediate_and_final_files(directory_final, directory_intermediate,
+                                             filename_final_no_ext, filename_with_iteration_no_ext,
+                                             report, '.txt')
 
 
 def _write_text_intermediate_and_final_files(directory_final: Optional[str], directory_intermediate: str,
@@ -989,24 +598,6 @@ def _write_text_intermediate_and_final_files(directory_final: Optional[str], dir
             file.write(content)
 
 
-def _write_report(design: Design, directory_intermediate: str, directory_final: str,
-                  filename_with_iteration: str, filename_final: str,
-                  report_only_violations: bool, violation_set: dc.ViolationSet) -> None:
-    report_str = design.summary_of_constraints(report_only_violations, violation_set=violation_set)
-    report = f'''\
-Report on constraints
-=====================
-{report_str}
-'''
-
-    for directory, filename in zip([directory_intermediate, directory_final],
-                                   [filename_with_iteration, filename_final]):
-        path = os.path.join(directory, filename)
-        with open(path, 'w') as file:
-            # file.write(sequences_content)
-            file.write(report)
-
-
 def _clear_directory(directory: str, force_overwrite: bool) -> None:
     files_relative = os.listdir(directory)
     files_and_directories = [os.path.join(directory, file) for file in files_relative]
@@ -1053,21 +644,19 @@ class _Directories:
     out: str
 
     # directories "fully qualified relative to project root": out joined with "subdirectory" strings below
-    dsd_design: str = field(init=False)
+    design: str = field(init=False)
     rng_state: str = field(init=False)
-    domain_pools: str = field(init=False)
     report: str = field(init=False)
     sequence: str = field(init=False)
 
     # relative to out directory
-    dsd_design_subdirectory: str = field(init=False, default='dsd_designs')
+    dsd_design_subdirectory: str = field(init=False, default='designs')
     rng_state_subdirectory: str = field(init=False, default='rng_state')
-    domain_pools_subdirectory: str = field(init=False, default='domain_pools')
     report_subdirectory: str = field(init=False, default='reports')
     sequence_subdirectory: str = field(init=False, default='sequences')
 
     # names of files to write (in subdirectories, and also "current-best" versions in out
-    dsd_design_filename_no_ext: str = field(init=False, default='dsd-design')
+    design_filename_no_ext: str = field(init=False, default='design')
     rng_state_filename_no_ext: str = field(init=False, default='rng-state')
     sequences_filename_no_ext: str = field(init=False, default='sequences')
     report_filename_no_ext: str = field(init=False, default='report')
@@ -1075,14 +664,20 @@ class _Directories:
     debug_file_handler: Optional[logging.FileHandler] = field(init=False, default=None)
     info_file_handler: Optional[logging.FileHandler] = field(init=False, default=None)
 
-    def all_subdirectories(self) -> List[str]:
-        return [self.dsd_design, self.rng_state, self.domain_pools, self.report, self.sequence]
+    def all_subdirectories(self, params: SearchParameters) -> List[str]:
+        result = []
+        if params.save_designs_for_all_updates:
+            result.extend([self.design, self.rng_state])
+        if params.save_sequences_for_all_updates:
+            result.append(self.report)
+        if params.save_reports_for_all_updates:
+            result.append(self.sequence)
+        return result
 
     def __init__(self, out: str, debug: bool, info: bool) -> None:
         self.out = out
-        self.dsd_design = os.path.join(self.out, self.dsd_design_subdirectory)
+        self.design = os.path.join(self.out, self.dsd_design_subdirectory)
         self.rng_state = os.path.join(self.out, self.rng_state_subdirectory)
-        self.domain_pools = os.path.join(self.out, self.domain_pools_subdirectory)
         self.report = os.path.join(self.out, self.report_subdirectory)
         self.sequence = os.path.join(self.out, self.sequence_subdirectory)
 
@@ -1096,6 +691,35 @@ def __init__(self, out: str, debug: bool, info: bool) -> None:
             self.info_file_handler.setLevel(logging.INFO)
             dc.logger.addHandler(self.info_file_handler)
 
+    @staticmethod
+    def indexed_full_filename_noext(filename_no_ext: str, directory: str, idx: int) -> str:
+        relative_filename = f'{filename_no_ext}-{idx}.json'
+        full_filename = os.path.join(directory, relative_filename)
+        return full_filename
+
+    def best_full_filename_noext(self, filename_no_ext: str) -> str:
+        relative_filename = f'current-best-{filename_no_ext}.json'
+        full_filename = os.path.join(self.out, relative_filename)
+        return full_filename
+
+    def indexed_design_full_filename_noext(self, idx: int) -> str:
+        return self.indexed_full_filename_noext(self.design_filename_no_ext, self.design, idx)
+
+    def indexed_rng_full_filename_noext(self, idx: int) -> str:
+        return self.indexed_full_filename_noext(self.rng_state_filename_no_ext, self.rng_state, idx)
+
+    def indexed_sequences_full_filename_noext(self, idx: int) -> str:
+        return self.indexed_full_filename_noext(self.sequences_filename_no_ext, self.sequence, idx)
+
+    def best_design_full_filename_noext(self) -> str:
+        return self.best_full_filename_noext(self.design_filename_no_ext)
+
+    def best_rng_full_filename_noext(self) -> str:
+        return self.best_full_filename_noext(self.rng_state_filename_no_ext)
+
+    def best_sequences_full_filename_noext(self) -> str:
+        return self.best_full_filename_noext(self.sequences_filename_no_ext)
+
 
 def _check_design(design: dc.Design) -> None:
     # verify design is legal
@@ -1157,17 +781,6 @@ class SearchParameters:
     whenever a new optimal sequence assignment is found.
     """
 
-    report_delay: float = 60.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
-    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
-    search finds many new improved designs frequently at the start of the search.
-    By setting this to 0, a new report will be written every time the design improves.
-    """
-
     on_improved_design: Callable[[int], None] = lambda _: None
     """
     Function to call whenever the design improves. Takes an integer as input indicating the number
@@ -1201,6 +814,9 @@ class SearchParameters:
 
     report_only_violations: bool = True
     """
+    NOTE: due to a recent change in how the search is conducted, setting this parameter to False is not
+    currently supported. 
+    
     If True, does not give report on each constraint that was satisfied; only reports violations
     and summary of all constraint checks of a certain type (e.g., how many constraint checks there were).
     """
@@ -1225,7 +841,7 @@ class SearchParameters:
     (for example, for sequences)
     sequences-0.txt, sequences-1.txt, ...,
     sequences-9.txt, sequences-10.txt, ...
-    If -nd 3 is specified, for instance, they will be written
+    If num_digits_update=3 is specified, for instance, they will be written
     sequences-000.txt, sequences-001.txt, ...,
     sequences-009.txt, sequences-010.txt, ...,
     sequences-099.txt, sequences-100.txt, ...,
@@ -1239,6 +855,27 @@ class SearchParameters:
     Log warning about sequences that are fixed, indicating they will not be re-assigned during the search.
     """
 
+    save_reports_for_all_updates: bool = True
+    """
+    A report on the most recently updated :any:`Design` is always written to a file 
+    `current-best-report.txt`. If this is True, then in the folder `reports`, a file unique to that update
+    is also written. Set to False to use less space on disk. 
+    """
+
+    save_designs_for_all_updates: bool = True
+    """
+    A serialized (JSON) description of the most recently updated :any:`Design` is always written to 
+    a file `current-best-design.json`. If this is True, then in the folder `dsd_designs`, a file unique to 
+    that update is also written. Set to False to use less space on disk. 
+    """
+
+    save_sequences_for_all_updates: bool = True
+    """
+    A list of sequences for each :any:`Strand` of most recently updated :any:`Design` is always written to 
+    a file `current-best-sequences.txt`. If this is True, then in the folder `sequences`, a file unique to 
+    that update is also written. Set to False to use less space on disk. 
+    """
+
 
 def search_for_dna_sequences(design: dc.Design, params: SearchParameters) -> None:
     """
@@ -1304,9 +941,7 @@ def search_for_dna_sequences(design: dc.Design, params: SearchParameters) -> Non
     design.compute_derived_fields()
     _check_design(design)
 
-    directories = _setup_directories(
-        debug=params.debug_log_file, info=params.info_log_file, force_overwrite=params.force_overwrite,
-        restart=params.restart, out_directory=params.out_directory)
+    directories = _setup_directories(params)
 
     if params.random_seed is not None:
         rng = np.random.default_rng(params.random_seed)
@@ -1346,19 +981,15 @@ def search_for_dna_sequences(design: dc.Design, params: SearchParameters) -> Non
 
         if not params.restart:
             # write initial sequences and report
-            _write_intermediate_files(design=design, rng=rng, num_new_optimal=num_new_optimal,
-                                      write_report=True, directories=directories,
-                                      report_only_violations=params.report_only_violations,
-                                      num_digits_update=params.num_digits_update,
-                                      violation_set=violation_set_opt)
+            _write_intermediate_files(design=design, params=params, rng=rng, num_new_optimal=num_new_optimal,
+                                      directories=directories, violation_set=violation_set_opt)
 
         # this helps with logging if we execute no iterations
         violation_set_new = violation_set_opt
 
         iteration = 0
-        time_of_last_improvement: float = -1.0
 
-        while len(violation_set_opt.non_fixed_violations) > 0 and \
+        while violation_set_opt.has_nonfixed_violations() and \
                 (params.max_iterations is None or iteration < params.max_iterations):
             _check_cpu_count(cpu_count)
 
@@ -1380,7 +1011,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
 
@@ -1394,23 +1025,8 @@ def search_for_dna_sequences(design: dc.Design, params: SearchParameters) -> Non
                 if score_delta < 0:  # increment whenever we actually improve the design
                     num_new_optimal += 1
                     on_improved_design(num_new_optimal)  # type: ignore
-
-                    current_time: float = time.time()
-                    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 current_time - time_of_last_improvement >= params.report_delay):  # > report_delay seconds since last report
-                        time_of_last_improvement = current_time
-                        write_report = True
-                        logger.debug('writing report')
-                    else:
-                        logger.debug('skipping report')
-
-                    _write_intermediate_files(design=design, rng=rng, num_new_optimal=num_new_optimal,
-                                              write_report=write_report, directories=directories,
-                                              report_only_violations=params.report_only_violations,
-                                              num_digits_update=params.num_digits_update,
+                    _write_intermediate_files(design=design, params=params, rng=rng,
+                                              num_new_optimal=num_new_optimal, directories=directories,
                                               violation_set=violation_set_opt)
 
             iteration += 1
@@ -1445,16 +1061,17 @@ def _check_cpu_count(cpu_count: int) -> None:
         _process_pool = new_process_pool(cpu_count)
 
 
-def _setup_directories(*, debug: bool, info: bool, force_overwrite: bool, restart: bool,
-                       out_directory: str) -> _Directories:
+def _setup_directories(params: SearchParameters) -> _Directories:
+    out_directory = params.out_directory
     if out_directory is None:
         out_directory = default_output_directory()
-    directories = _Directories(out=out_directory, debug=debug, info=info)
+    directories = _Directories(out=out_directory, debug=params.debug_log_file,
+                               info=params.info_log_file)
     if not os.path.exists(directories.out):
         os.makedirs(directories.out)
-    if not restart:
-        _clear_directory(directories.out, force_overwrite)
-    for subdir in directories.all_subdirectories():
+    if not params.restart:
+        _clear_directory(directories.out, params.force_overwrite)
+    for subdir in directories.all_subdirectories(params):
         if not os.path.exists(subdir):
             os.makedirs(subdir)
     return directories
@@ -1478,20 +1095,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
@@ -1547,22 +1165,26 @@ def _restart_from_directory(directories: _Directories, design: dc.Design) \
     # NOTE: restarts from highest index found in dsd_design subdirectory, NOT from "current-best" files,
     # which are ignored. This applies to both the design and the RNG state
 
-    # get DomainPools
-    pool_with_name = read_domain_pools(directories)
+    highest_idx_design = 0
 
-    # returns highest index found, as well as name of corresponding design file
-    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,
-                                   f'{directories.dsd_design_filename_no_ext}-{highest_idx_design}.json')
+    if os.path.isdir(directories.design):
+        # returns highest index found in design subdirectory
+        highest_idx_design = _find_highest_index_in_directory(directories.design,
+                                                              directories.design_filename_no_ext, 'json')
+        design_filename = directories.indexed_design_full_filename_noext(highest_idx_design)
+        rng_filename = directories.indexed_rng_full_filename_noext(highest_idx_design)
+    else:
+        # otherwise we go with contents of "current-best-*.json"
+        design_filename = directories.best_design_full_filename_noext()
+        rng_filename = directories.best_rng_full_filename_noext()
+
+    # read design
     with open(design_filename, 'r') as file:
         design_json_str = file.read()
-    design_stored = dc.Design.from_json(design_json_str, pool_with_name=pool_with_name)
+    design_stored = dc.Design.from_json(design_json_str)
     dc.verify_designs_match(design_stored, design, check_fixed=False)
 
-    # get RNG state
-    rng_filename = os.path.join(directories.rng_state,
-                                f'{directories.rng_state_filename_no_ext}-{highest_idx_design}.json')
+    # read RNG state
     with open(rng_filename, 'r') as file:
         rng_state_json = file.read()
     rng_state = json.loads(rng_state_json)
@@ -1575,54 +1197,23 @@ def _restart_from_directory(directories: _Directories, design: dc.Design) \
     # design_stored.copy_constraints_from(design)
     design.copy_sequences_from(design_stored)
 
-    design_json = json.loads(design_json_str)
-    stored_pool_idxs = design_json[dc.domain_pools_num_sampled_key]
-    for pool in design_stored.domain_pools():
-        idx = stored_pool_idxs[pool.name]
-        pool.num_sampled = idx
-
     return highest_idx_design, rng
 
 
-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
-
-    # first find pool names from dsd_design files
-    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,
-                                   f'{directories.dsd_design_filename_no_ext}-{highest_idx_design}.json')
-    with open(design_filename, 'r') as file:
-        design_json_str = file.read()
-    design_json = json.loads(design_json_str)
-    domains_json = design_json[dc.domains_key]
-    pool_names: Set[str] = set()
-    for domain_json in domains_json:
-        if dc.domain_pool_name_key in domain_json:
-            pool_name = domain_json[dc.domain_pool_name_key]
-            pool_names.add(pool_name)
-
-    pool_with_name = {}
-    for pool_name in pool_names:
-        highest_idx_pool = _find_highest_index_in_directory(directories.domain_pools,
-                                                            pool_name, 'json')
-        pool_filename = os.path.join(directories.domain_pools, f'{pool_name}-{highest_idx_pool}.json')
-        with open(pool_filename, 'r') as file:
-            pool_json_str = file.read()
-        pool_json = json.loads(pool_json_str)
-        pool = dc.DomainPool.from_json_serializable(pool_json)
-        assert pool.name == pool_name
-        pool_with_name[pool_name] = pool
-
-    return pool_with_name
-
-
 def _find_highest_index_in_directory(directory: str, filename_start: str, ext: str) -> int:
     # return highest index of filename (name matches "<filename_start>-<index>.<ext>"
-    filenames = [filename
-                 for filename in os.listdir(directory)
-                 if os.path.isfile(os.path.join(directory, filename))]
+    # raises ValueError if none exists
+    try:
+        list_dir = os.listdir(directory)
+    except FileNotFoundError:
+        list_dir = None
+    if list_dir is not None and len(list_dir) > 0:
+        filenames = [filename
+                     for filename in list_dir
+                     if os.path.isfile(os.path.join(directory, filename))]
+    else:
+        raise ValueError(f'no files in directory "{directory}" '
+                         f'match the pattern "{filename_start}-<index>.{ext}";\n')
 
     pattern = re.compile(filename_start + r'-(\d+)\.' + ext)
     filenames_matching = [filename for filename in filenames if pattern.search(filename)]
@@ -1644,38 +1235,6 @@ def _find_highest_index_in_directory(directory: str, filename_start: str, ext: s
     return max_index
 
 
-def _write_intermediate_files(*, design: dc.Design, rng: numpy.random.Generator,
-                              num_new_optimal: int, write_report: bool,
-                              directories: _Directories, report_only_violations: bool,
-                              num_digits_update: Optional[int],
-                              violation_set: dc.ViolationSet) -> None:
-    num_new_optimal_padded = f'{num_new_optimal}' if num_digits_update is None \
-        else f'{num_new_optimal:0{num_digits_update}d}'
-
-    _write_dsd_design_json(design, directories=directories, num_new_optimal_padded=num_new_optimal_padded,
-                           directory_final=directories.out,
-                           filename_final_no_ext=f'current-best-{directories.dsd_design_filename_no_ext}')
-    _write_rng_state(rng, directories, num_new_optimal_padded)
-    _write_domain_pools(design.domain_pools_to_domain_map.keys(), directories)
-
-    _write_sequences(design,
-                     directory_intermediate=directories.sequence,
-                     directory_final=directories.out,
-                     filename_with_iteration_no_ext=f'{directories.sequences_filename_no_ext}'
-                                                    f'-{num_new_optimal_padded}',
-                     filename_final_no_ext=f'current-best-{directories.sequences_filename_no_ext}')
-
-    if write_report:
-        _write_report(design,
-                      directory_intermediate=directories.report,
-                      directory_final=directories.out,
-                      filename_with_iteration=f'{directories.report_filename_no_ext}'
-                                              f'-{num_new_optimal_padded}.txt',
-                      filename_final=f'current-best-{directories.report_filename_no_ext}.txt',
-                      report_only_violations=report_only_violations,
-                      violation_set=violation_set)
-
-
 def _pfunc_killall() -> None:
     import subprocess
     delim = '#'
@@ -1843,8 +1402,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 +1413,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..254c1fe2 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,53 +69,56 @@ 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')
+        threshold=-2.0, temperature=52, short_description='DomainPairRNA')
 
     domain_pair_nupack_constraint = dc.nupack_domain_pair_constraint(
-        threshold=-0.5, temperature=52, short_description='DomainPairNoCompl',
-        threaded=threaded)
+        threshold=-0.5, temperature=52, short_description='DomainPairNUPACK',
+        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='StrandPairRNA', 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='StrandPairNUPACK', parallel=parallel, weight=0.1)
 
     design = dc.Design(strands,
                        constraints=[
-                           # strand_individual_ss_constraint,
+                           # domain_nupack_ss_constraint,
+                           strand_individual_ss_constraint,
                            # strand_pair_nupack_constraint,
                            # domain_pair_nupack_constraint,
-                           domain_pairs_rna_duplex_constraint,
-                           # dc.domains_not_substrings_of_each_other_domain_pair_constraint(),
+                           # domain_pairs_rna_duplex_constraint,
                            # strand_pairs_rna_duplex_constraint,
+                           # dc.domains_not_substrings_of_each_other_domain_pair_constraint(),
                        ])
 
     numpy_constraints: List[NumpyConstraint] = [
-        # dc.NearestNeighborEnergyConstraint(-9.5, -9.0, 52.0),
+        dc.NearestNeighborEnergyConstraint(-9.3, -9.0, 52.0),
         # dc.BaseCountConstraint(base='G', high_count=1),
         # dc.BaseEndConstraint(bases=('C', 'G')),
-        # dc.RunsOfBasesConstraint(['C', 'G'], 4),
+        dc.RunsOfBasesConstraint(['C', 'G'], 4),
         # dc.RunsOfBasesConstraint(['A', 'T'], 4),
         # dc.BaseEndConstraint(bases=('A', 'T')),
         # dc.BaseEndConstraint(bases=('C', 'G'), distance_from_end=1),
         # dc.BaseAtPositionConstraint(bases='T', position=3),
-        dc.ForbiddenSubstringConstraint(['GGGG', 'CCCC']),
+        # dc.ForbiddenSubstringConstraint(['GGGG', 'CCCC']),
         # dc.RestrictBasesConstraint(bases=['A', 'T', 'C']),
     ]
 
@@ -131,38 +134,46 @@ def main() -> None:
 
     replace_with_close_sequences = True
     # replace_with_close_sequences = False
-    domain_pool_10 = dc.DomainPool(f'length-{10}_domains', 10,
+    domain_pool_10 = dc.DomainPool(f'length-10_domains', 10,
                                    numpy_constraints=numpy_constraints,
                                    replace_with_close_sequences=replace_with_close_sequences,
                                    )
-    domain_pool_11 = dc.DomainPool(f'length-{11}_domains', 11,
+    domain_pool_11 = dc.DomainPool(f'length-11_domains', 11,
                                    numpy_constraints=numpy_constraints,
                                    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,
-                                 report_only_violations=True,
                                  random_seed=random_seed,
                                  max_iterations=None,
-                                 # force_overwrite=True,
-                                 # report_delay=0.0,
+                                 save_sequences_for_all_updates=False,
+                                 save_reports_for_all_updates=False,
+                                 save_designs_for_all_updates=False,
+                                 force_overwrite=True,
                                  )
-    params.force_overwrite = True
     ds.search_for_dna_sequences(design, params)
 
 
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)