From 72393e7c243dd431f511ad12be0c8b6945d5482e Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Mon, 8 Jun 2026 16:32:19 -0500
Subject: [PATCH 01/13] make vmd call optional

---
 frustratometer/classes/Frustratometer.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/frustratometer/classes/Frustratometer.py b/frustratometer/classes/Frustratometer.py
index 11c22c5..f4c32ba 100644
--- a/frustratometer/classes/Frustratometer.py
+++ b/frustratometer/classes/Frustratometer.py
@@ -421,6 +421,8 @@ def vmd(self, sequence: str = None, single:Union[str,np.array] = 'singleresidue'
                                       -self.frustration(kind=pair, sequence=sequence, aa_freq=aa_freq),
                                       max_connections=max_connections, movie_name=movie_name, still_image_name=still_image_name)
         frustration.call_vmd(self.pdb_file, tcl_script)
+        if call_vmd:
+            frustration.call_vmd(self.pdb_file, tcl_script)
 
     def view_pair_frustration(self, sequence:str = None, pair:str = 'mutational', aa_freq:np.array = None):
         """

From a82b3e2e7a92661c36cfd0fd3a167a5a209e0631 Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Mon, 8 Jun 2026 16:33:29 -0500
Subject: [PATCH 02/13] make vmd call optional

---
 frustratometer/classes/Frustratometer.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/frustratometer/classes/Frustratometer.py b/frustratometer/classes/Frustratometer.py
index f4c32ba..4288120 100644
--- a/frustratometer/classes/Frustratometer.py
+++ b/frustratometer/classes/Frustratometer.py
@@ -390,7 +390,8 @@ def auc(self):
         return frustration.compute_auc(self.roc())
 
     def vmd(self, sequence: str = None, single:Union[str,np.array] = 'singleresidue', pair:Union[str,np.array] = 'mutational',
-             aa_freq:np.array = None, correction:int = 0, max_connections:Union[int,None] = None, movie_name=None, still_image_name=None):
+             aa_freq:np.array = None, correction:int = 0, max_connections:Union[int,None] = None, movie_name=None, still_image_name=None,
+             call_vmd: bool = True):
         """
         Calculates frustration indices and superimposes frustration patterns onto PDB structure using the VMD software.
 

From a8c513f78d6054df693178b1a0918d7f06434438 Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Mon, 8 Jun 2026 16:34:01 -0500
Subject: [PATCH 03/13] make vmd call optional

---
 frustratometer/classes/Frustratometer.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/frustratometer/classes/Frustratometer.py b/frustratometer/classes/Frustratometer.py
index 4288120..32d66b4 100644
--- a/frustratometer/classes/Frustratometer.py
+++ b/frustratometer/classes/Frustratometer.py
@@ -421,7 +421,6 @@ def vmd(self, sequence: str = None, single:Union[str,np.array] = 'singleresidue'
                                       -self.frustration(kind=single, sequence=sequence, aa_freq=aa_freq),
                                       -self.frustration(kind=pair, sequence=sequence, aa_freq=aa_freq),
                                       max_connections=max_connections, movie_name=movie_name, still_image_name=still_image_name)
-        frustration.call_vmd(self.pdb_file, tcl_script)
         if call_vmd:
             frustration.call_vmd(self.pdb_file, tcl_script)
 

From 7bb6fb4fd1f691f1cd30527e25a3009e9475587d Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Thu, 11 Jun 2026 16:56:45 -0500
Subject: [PATCH 04/13] add __iter__ to SparseMatrix

---
 frustratometer/classes/Structure.py | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/frustratometer/classes/Structure.py b/frustratometer/classes/Structure.py
index 64bd390..1cd1f85 100644
--- a/frustratometer/classes/Structure.py
+++ b/frustratometer/classes/Structure.py
@@ -101,6 +101,22 @@ def __repr__(self):
         kind = 'mask' if self.data is None else 'distance'
         return f"SparseMatrix({kind}, nnz={len(self)}, L={self.shape})"
 
+    def __iter__(self):
+        # check that our attributes all make sense
+        assert self.row.shape == self.col.shape == (self.row.shape[0],) # tuple of length 1
+        if type(self.data) == np.ndarray:
+            assert self.data.shape == self.row.shape
+        elif not (self.data is None):
+            raise AssertionError(f'self.data was {self.data}')
+        N = self.row.shape[0]
+        counter = 0
+        while counter < N:
+            if self.data is None:
+                yield (self.row[counter], self.col[counter], self.data) 
+            else: # numpy.ndarray
+                yield (self.row[counter], self.col[counter], self.data[counter])
+            counter += 1
+
 class Structure:
 
     def __init__(self, pdb_file: Union[Path,str], chain: Union[str,None]=None, seq_selection: str = None, aligned_sequence: str = None, filtered_aligned_sequence: str = None,

From c081c010538f69e349eb050000cbbd896d0db9c6 Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Fri, 12 Jun 2026 17:22:09 -0500
Subject: [PATCH 05/13] added __and__ to SparseMatrix

---
 frustratometer/classes/Structure.py | 57 +++++++++++++++++++++++++----
 1 file changed, 50 insertions(+), 7 deletions(-)

diff --git a/frustratometer/classes/Structure.py b/frustratometer/classes/Structure.py
index 1cd1f85..b87df15 100644
--- a/frustratometer/classes/Structure.py
+++ b/frustratometer/classes/Structure.py
@@ -109,13 +109,56 @@ def __iter__(self):
         elif not (self.data is None):
             raise AssertionError(f'self.data was {self.data}')
         N = self.row.shape[0]
-        counter = 0
-        while counter < N:
-            if self.data is None:
-                yield (self.row[counter], self.col[counter], self.data) 
-            else: # numpy.ndarray
-                yield (self.row[counter], self.col[counter], self.data[counter])
-            counter += 1
+        n = 0
+        while n < N:
+            yield (self.row[n], self.col[n], None if self.data is None else self.data[n]) 
+            n += 1
+
+    def __and__(self, other):
+        # create a new SparseMatrix that acts as a boolean mask
+        # for the indices present in self and other
+        if not isinstance(other, SparseMatrix):
+            return NotImplemented
+        else:
+            if self.shape != other.shape:
+                raise ValueError("intersection of two boolean arrays of different shapes is ill-defined")
+            new_row, new_col = self._and_helper(other)
+            return SparseMatrix(new_row, new_col)
+
+    def __rand__(self, other):
+        if not isinstance(other, SparseMatrix):
+            return NotImplemented
+        else:
+            return self.__and__(other)
+
+    def __iand__(self, other):
+        # in-place modification of self
+        if not isinstance(other, SparseMatrix):
+            return NotImplemented
+        if self.shape != other.shape:
+            raise ValueError("intersection of two boolean arrays of different shapes is ill-defined")
+        self.row, self.col = self._and_helper(other) # modify self.row and self.col in place
+        self.data = None
+        return self
+
+    def _and_helper(self, other):
+        #new_row = []
+        #new_col = []
+        #other_set = set(zip(other.row, other.col))
+        #for i, j in zip(self.row, self.col):
+        #    if (i,j) in other_set:
+        #        new_row.append(i)
+        #        new_col.append(j)
+        #return np.array(new_row), np.array(new_col)
+        
+        # Vectorized membership using 64-bit integer keys; preserves order of self.
+        # Similar speed to the above option, may be more memory-efficient
+        # if numpy significantly outperforms native python
+        L = self.shape 
+        keys_self = (self.row.astype(np.int64) * L) + self.col.astype(np.int64)
+        keys_other = (other.row.astype(np.int64) * L) + other.col.astype(np.int64)
+        mask = np.in1d(keys_self, keys_other, assume_unique=False)
+        return self.row[mask].astype(np.intp), self.col[mask].astype(np.intp)
 
 class Structure:
 

From c7f3a1e5a8f5349afeb2a0b97716429346ce6d1d Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Fri, 12 Jun 2026 17:23:02 -0500
Subject: [PATCH 06/13] updated SparseMatrix __iter__

---
 frustratometer/classes/Structure.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/frustratometer/classes/Structure.py b/frustratometer/classes/Structure.py
index b87df15..ad83746 100644
--- a/frustratometer/classes/Structure.py
+++ b/frustratometer/classes/Structure.py
@@ -104,7 +104,7 @@ def __repr__(self):
     def __iter__(self):
         # check that our attributes all make sense
         assert self.row.shape == self.col.shape == (self.row.shape[0],) # tuple of length 1
-        if type(self.data) == np.ndarray:
+        if isinstance(self.data, np.ndarray):
             assert self.data.shape == self.row.shape
         elif not (self.data is None):
             raise AssertionError(f'self.data was {self.data}')

From c6e60e4c99cd16be9533d8b7c0515a884511ba7a Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Fri, 12 Jun 2026 17:50:09 -0500
Subject: [PATCH 07/13] update frustration.write_tcl_script to handle sparse
 matrices

---
 frustratometer/frustration/frustration.py | 56 +++++++++++++++++++----
 1 file changed, 47 insertions(+), 9 deletions(-)

diff --git a/frustratometer/frustration/frustration.py b/frustratometer/frustration/frustration.py
index 08fd1d8..433b35d 100644
--- a/frustratometer/frustration/frustration.py
+++ b/frustratometer/frustration/frustration.py
@@ -894,6 +894,10 @@ def write_tcl_script(pdb_file: Union[Path,str], chain: str, mask: np.array, dist
                     pair_frustration: np.array, tcl_script: Union[Path, str] ='frustration.tcl',max_connections: int =None, movie_name: Union[Path, str] =None, still_image_name: Union[Path, str] =None) -> Union[Path, str]:
     """
     Writes a tcl script that can be run with VMD to superimpose the frustration patterns onto the corresponding PDB structure. 
+    
+    This function is fastest when distance_matrix and mask inputs are dense.
+    Converting the matrices from sparse to dense or alternatively following the sparse matrix-
+    based logic adds about 0.4 seconds to this function call for NFKB (615 residues).
 
     Parameters
     ----------
@@ -902,9 +906,9 @@ def write_tcl_script(pdb_file: Union[Path,str], chain: str, mask: np.array, dist
     chain : str
         Select chain from pdb
     mask : np.array
-        A 2D Boolean array that determines which residue pairs should be considered in the energy computation. The mask should have dimensions (L, L), where L is the length of the sequence.
+        SparseMatrix or 2D Boolean array that determines which residue pairs should be considered in the energy computation. The mask should have dimensions (L, L), where L is the length of the sequence.
     distance_matrix : np.array
-        LxL array for sequence of length L, describing distances between contacts
+        SparseMatrix or LxL array for sequence of length L, describing distances between contacts
     distance_cutoff : float
         Maximum distance at which a contact occurs
     single_frustration : np.array
@@ -950,13 +954,47 @@ def write_tcl_script(pdb_file: Union[Path,str], chain: str, mask: np.array, dist
     for r, f in zip(residues, single_frustration):
         fo.write(f'[atomselect top "residue {int(r)}"] set beta {f}\n')
 
-    ii, jj = np.triu_indices(len(residues), k=1)
-    keep = mask[ii, jj] > 0
-    if distance_cutoff:
-        keep &= distance_matrix[ii, jj] <= distance_cutoff
-    ii, jj = ii[keep], jj[keep]
-    frust = pair_frustration[ii, jj]
-    dist = distance_matrix[ii, jj]
+    # if we're really worried about memory,
+    # we should loop over elements of mask
+    # and distance_matrix and check their indices
+    # and values, but distance_matrix shouldn't
+    # be very big (10^4 residues -> 10^8
+    # elements in distance matrix -> 800 MB,
+    # assuming 8-byte floats)
+    try:
+        mask = mask.to_dense(fill=False)
+    except AttributeError:
+        pass # it's already a numpy array
+    try:
+        distance_matrix = distance_matrix.to_dense(fill=np.inf)
+    except AttributeError:
+        pass # it's already a numpy array
+    def process_dense_arrays():
+        ii, jj = np.triu_indices(len(residues), k=1)
+        keep = mask[ii, jj] > 0
+        if distance_cutoff:
+            keep &= distance_matrix[ii, jj] <= distance_cutoff
+        ii, jj = ii[keep], jj[keep]
+        frust = pair_frustration[ii, jj]
+        dist = distance_matrix[ii, jj]
+        return frust, dist, ii, jj
+    if isinstance(distance_matrix,np.ndarray) and isinstance(mask,np.ndarray):
+        frust, dist, ii, jj = process_dense_arrays()
+    elif isinstance(distance_matrix,np.ndarray) and not isinstance(mask,np.ndarray):
+        mask = mask.to_dense(fill=False)
+        frust, dist, ii, jj = process_dense_arrays()
+    elif (not isinstance(distance_matrix,np.ndarray)) and isinstance(mask,np.ndarray):
+        distance_matrix = distance_matrix.to_dense(fill=np.inf)
+        frust, dist, ii, jj = process_dense_arrays()
+    else: # This is about the same speed as converting to dense matrix for NFKB (615 residues)
+          # note that this code is currently unreachable
+        if distance_cutoff:
+            mask &= distance_matrix.compute_mask(maximum_contact_distance=distance_cutoff)
+        upper_triangular_part = np.unique(np.sort(np.array([mask.row, mask.col]),axis=0),axis=1)
+        ii = upper_triangular_part[0,:]
+        jj = upper_triangular_part[1,:]
+        frust = pair_frustration[ii,jj] # pair_frustration is always dense
+        dist = distance_matrix.lookup(ii,jj)
 
     # Green = minimally frustrated (most negative first); red = frustrated (most positive
     # first). For each: sort, cap at max_connections, then draw the in-range contacts.

From fc5b2023a3ac8f742c7cf8fa94a7da18b35174ee Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Fri, 12 Jun 2026 18:23:32 -0500
Subject: [PATCH 08/13] reduce number of write() calls in write_tcl_script (may
 be marginally faster)

---
 frustratometer/frustration/frustration.py | 26 +++++++++++++++--------
 1 file changed, 17 insertions(+), 9 deletions(-)

diff --git a/frustratometer/frustration/frustration.py b/frustratometer/frustration/frustration.py
index 433b35d..044e081 100644
--- a/frustratometer/frustration/frustration.py
+++ b/frustratometer/frustration/frustration.py
@@ -930,7 +930,8 @@ def write_tcl_script(pdb_file: Union[Path,str], chain: str, mask: np.array, dist
     tcl_script : Path or str
         tcl script file
     """
-    fo = open(tcl_script, 'w+')
+    to_write = []
+    #fo = open(tcl_script, 'w+')
     single_frustration = np.nan_to_num(single_frustration,nan=0,posinf=0,neginf=0)
     pair_frustration = np.nan_to_num(pair_frustration,nan=0,posinf=0,neginf=0)
     
@@ -949,10 +950,10 @@ def write_tcl_script(pdb_file: Union[Path,str], chain: str, mask: np.array, dist
     ca_xyz[np.searchsorted(residues, ca_sel.getResindices())] = ca_sel.getCoords()
     cb_xyz[np.searchsorted(residues, cb_sel.getResindices())] = cb_sel.getCoords()
 
-    fo.write(f'[atomselect top all] set beta 0\n')
+    to_write.append(f'[atomselect top all] set beta 0\n')
     # Single residue frustration
     for r, f in zip(residues, single_frustration):
-        fo.write(f'[atomselect top "residue {int(r)}"] set beta {f}\n')
+        to_write.append(f'[atomselect top "residue {int(r)}"] set beta {f}\n')
 
     # if we're really worried about memory,
     # we should loop over elements of mask
@@ -1019,13 +1020,13 @@ def process_dense_arrays():
         cb_dist = np.linalg.norm(cb_xyz[gi] - cb_xyz[gj], axis=1)   # Cbeta-Cbeta -> dashing
         styles = np.where((cb_dist >= 3.5) & (cb_dist <= 6.5), 'solid', 'dashed')
 
-        fo.write(f'draw color {color}\n')
+        to_write.append(f'draw color {color}\n')
         for a, b, style in zip(gi, gj, styles):
             p1, p2 = ca_xyz[a], ca_xyz[b]
-            fo.write("draw line {%.3f %.3f %.3f} {%.3f %.3f %.3f} style %s width 2\n"
+            to_write.append("draw line {%.3f %.3f %.3f} {%.3f %.3f %.3f} style %s width 2\n"
                      % (p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], style))
 
-    fo.write('''mol delrep top 0
+        to_write.append('''mol delrep top 0
             mol color Beta
             mol representation NewCartoon 0.300000 10.000000 4.100000 0
             mol selection all
@@ -1035,7 +1036,7 @@ def process_dense_arrays():
             ''')
 
     if movie_name:
-        fo.write('''axes location Off
+            to_write.append('''axes location Off
             color Display Background white
             display resize 800 800
             display projection Orthographic
@@ -1090,11 +1091,18 @@ def process_dense_arrays():
             exit
         ''')
     elif still_image_name:
-        fo.write(f'set output "{still_image_name}"' + '''
+        to_write.append(f'set output "{still_image_name}"' + '''
             render snapshot $output
             exit
         ''')
-    fo.close()
+    # Write the tcl script once at the end
+    try:
+        with open(tcl_script, 'w+') as _fo:
+            _fo.write(''.join(to_write))
+    except Exception:
+        # Fallback: ensure we raise the original error if file write fails
+        with open(str(tcl_script), 'w+') as _fo:
+            _fo.write(''.join(to_write))
     return tcl_script
 
 

From 07b555d5ec93e2567b0d7f10566e734091d545ad Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Sun, 21 Jun 2026 17:10:26 -0500
Subject: [PATCH 09/13] update Frustratometer.vmd() to actually support
 np.array-type args for single and pair, as promised in the function type
 hints

---
 frustratometer/classes/Frustratometer.py | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/frustratometer/classes/Frustratometer.py b/frustratometer/classes/Frustratometer.py
index 32d66b4..35ac5c1 100644
--- a/frustratometer/classes/Frustratometer.py
+++ b/frustratometer/classes/Frustratometer.py
@@ -389,7 +389,7 @@ def auc(self):
         """
         return frustration.compute_auc(self.roc())
 
-    def vmd(self, sequence: str = None, single:Union[str,np.array] = 'singleresidue', pair:Union[str,np.array] = 'mutational',
+    def vmd(self, sequence: str = None, single:Union[str,np.array,None] = 'singleresidue', pair:Union[str,np.array] = 'mutational',
              aa_freq:np.array = None, correction:int = 0, max_connections:Union[int,None] = None, movie_name=None, still_image_name=None,
              call_vmd: bool = True):
         """
@@ -416,11 +416,17 @@ def vmd(self, sequence: str = None, single:Union[str,np.array] = 'singleresidue'
                     from the sequence that was passed to this vmd function. Proceeding further may not\n\
                     perform the computation that you intend to perform.")
         
+        if isinstance(single, str):
+            single = -self.frustration(kind=single, sequence=sequence, aa_freq=aa_freq)
 
-        tcl_script = frustration.write_tcl_script(self.pdb_file, self.chain, self.mask, self.distance_matrix, self.distance_cutoff,
-                                      -self.frustration(kind=single, sequence=sequence, aa_freq=aa_freq),
-                                      -self.frustration(kind=pair, sequence=sequence, aa_freq=aa_freq),
-                                      max_connections=max_connections, movie_name=movie_name, still_image_name=still_image_name)
+        if isinstance(pair, str):
+            pair = -self.frustration(kind=pair, sequence=sequence, aa_freq=aa_freq)
+
+        tcl_script = frustration.write_tcl_script(
+            self.pdb_file, self.chain, self.mask, self.distance_matrix, self.distance_cutoff,
+            single, pair,
+            max_connections=max_connections, movie_name=movie_name, still_image_name=still_image_name)
+        
         if call_vmd:
             frustration.call_vmd(self.pdb_file, tcl_script)
 

From 611fbb2ca64c0e37e2db18a23fea1f5e702f6b1b Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Sun, 21 Jun 2026 17:18:46 -0500
Subject: [PATCH 10/13] support None for single_frustration argument to
 write_tcl_script

---
 frustratometer/frustration/frustration.py | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/frustratometer/frustration/frustration.py b/frustratometer/frustration/frustration.py
index 044e081..3be967a 100644
--- a/frustratometer/frustration/frustration.py
+++ b/frustratometer/frustration/frustration.py
@@ -890,8 +890,13 @@ def plot_singleresidue_decoy_energy(decoy_energy, native_energy, method='cluster
     return g
 
 
-def write_tcl_script(pdb_file: Union[Path,str], chain: str, mask: np.array, distance_matrix: np.array, distance_cutoff: float, single_frustration: np.array,
-                    pair_frustration: np.array, tcl_script: Union[Path, str] ='frustration.tcl',max_connections: int =None, movie_name: Union[Path, str] =None, still_image_name: Union[Path, str] =None) -> Union[Path, str]:
+def write_tcl_script(
+    pdb_file: Union[Path,str], chain: str, mask: np.array, 
+    distance_matrix: np.array, distance_cutoff: float, 
+    single_frustration: Union[np.array,None], pair_frustration: np.array, 
+    tcl_script: Union[Path, str] ='frustration.tcl',max_connections: int =None, 
+    movie_name: Union[Path, str] =None, still_image_name: Union[Path, str] =None,
+    ) -> Union[Path, str]:
     """
     Writes a tcl script that can be run with VMD to superimpose the frustration patterns onto the corresponding PDB structure. 
     
@@ -952,13 +957,9 @@ def write_tcl_script(pdb_file: Union[Path,str], chain: str, mask: np.array, dist
 
     to_write.append(f'[atomselect top all] set beta 0\n')
     # Single residue frustration
-    for r, f in zip(residues, single_frustration):
-        to_write.append(f'[atomselect top "residue {int(r)}"] set beta {f}\n')
-
-    # if we're really worried about memory,
-    # we should loop over elements of mask
-    # and distance_matrix and check their indices
-    # and values, but distance_matrix shouldn't
+    if not (single_frustration is None):
+        for r, f in zip(residues, single_frustration):
+            to_write.append(f'[atomselect top "residue {int(r)}"] set beta {f}\n')
     # be very big (10^4 residues -> 10^8
     # elements in distance matrix -> 800 MB,
     # assuming 8-byte floats)

From 37e438054aaa539add9eb91aa544327f6bc271dc Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Sun, 21 Jun 2026 17:28:00 -0500
Subject: [PATCH 11/13] finish sparse matrix updates to
 frustration.write_tcl_script()

---
 frustratometer/frustration/frustration.py | 108 ++++++++++------------
 1 file changed, 48 insertions(+), 60 deletions(-)

diff --git a/frustratometer/frustration/frustration.py b/frustratometer/frustration/frustration.py
index 3be967a..fefdd47 100644
--- a/frustratometer/frustration/frustration.py
+++ b/frustratometer/frustration/frustration.py
@@ -960,6 +960,11 @@ def write_tcl_script(
     if not (single_frustration is None):
         for r, f in zip(residues, single_frustration):
             to_write.append(f'[atomselect top "residue {int(r)}"] set beta {f}\n')
+
+    # if we're really worried about memory,
+    # we should loop over elements of mask
+    # and distance_matrix and check their indices
+    # and values, but distance_matrix shouldn't
     # be very big (10^4 residues -> 10^8
     # elements in distance matrix -> 800 MB,
     # assuming 8-byte floats)
@@ -1037,73 +1042,56 @@ def process_dense_arrays():
             ''')
 
     if movie_name:
-            to_write.append('''axes location Off
-            color Display Background white
-            display resize 800 800
-            display projection Orthographic
-            display depthcue off
-            display resetview
-            display resize [expr [lindex [display get size] 0]/2*2] [expr [lindex [display get size] 1]/2*2] ;#Resize display to even height and width
-            display update ui
-
-            # Set up the movie directory and base file name
-            mkdir movie_tmp
-            set workdir "movie_tmp"
-            ''' + f'set basename "{movie_name}"' + '''
-            set numframes 360
-            set framerate 25
-
-            # Function to rotate the molecule and capture frames
-            proc captureFrames {} {
-                global workdir basename numframes
-                for {set i 0} {$i < $numframes} {incr i} {
-                    # Rotate the molecule around the Y-axis
-                    rotate y by 1
-                    
-                    # Capture the frame
-                    set output [format "%s/$basename.%05d.tga" $workdir $i]
-                    render snapshot $output
-                }
-            }
-
-            # Function to convert frames to MP4
-            proc convertToMP4 {} {
-                global workdir basename numframes framerate
-
-                set mybasefilename [format "%s/%s" $workdir $basename]
-                set outputFile [format "%s.mp4" $basename]
-                
-                # Construct and execute the ffmpeg command
-                
-                set command "ffmpeg -y -framerate $framerate -i $mybasefilename.%05d.tga -c:v libx264 -profile:v high -crf 20 -pix_fmt yuv420p $outputFile"
-                puts "Executing: $command"
-                exec ffmpeg -y -framerate $framerate -i $mybasefilename.%05d.tga -c:v libx264 -profile:v high -crf 20 -pix_fmt yuv420p $outputFile >&@ stdout
-            }
-
-            # Main script execution
-            captureFrames
-            convertToMP4
-
-            # Cleanup the TGA files if desired
+        to_write.append('''axes location Off
+        color Display Background white
+        display resize 800 800
+        display projection Orthographic
+        display depthcue off
+        display resetview
+        display resize [expr [lindex [display get size] 0]/2*2] [expr [lindex [display get size] 1]/2*2]
+        display update ui
+
+        mkdir movie_tmp
+        set workdir "movie_tmp"
+        ''' + f'set basename "{movie_name}"' + '''
+        set numframes 360
+        set framerate 25
+
+        proc captureFrames {} {
+            global workdir basename numframes
             for {set i 0} {$i < $numframes} {incr i} {
+                rotate y by 1
                 set output [format "%s/$basename.%05d.tga" $workdir $i]
-                exec rm $output
+                render snapshot $output
             }
-            exit
+        }
+
+        proc convertToMP4 {} {
+            global workdir basename numframes framerate
+            set mybasefilename [format "%s/%s" $workdir $basename]
+            set outputFile [format "%s.mp4" $basename]
+            set command "ffmpeg -y -framerate $framerate -i $mybasefilename.%05d.tga -c:v libx264 -profile:v high -crf 20 -pix_fmt yuv420p $outputFile"
+            puts "Executing: $command"
+            exec ffmpeg -y -framerate $framerate -i $mybasefilename.%05d.tga -c:v libx264 -profile:v high -crf 20 -pix_fmt yuv420p $outputFile >&@ stdout
+        }
+
+        captureFrames
+        convertToMP4
+
+        for {set i 0} {$i < $numframes} {incr i} {
+            set output [format "%s/$basename.%05d.tga" $workdir $i]
+            exec rm $output
+        }
+        exit
         ''')
     elif still_image_name:
         to_write.append(f'set output "{still_image_name}"' + '''
-            render snapshot $output
-            exit
+        render snapshot $output
+        exit
         ''')
-    # Write the tcl script once at the end
-    try:
-        with open(tcl_script, 'w+') as _fo:
-            _fo.write(''.join(to_write))
-    except Exception:
-        # Fallback: ensure we raise the original error if file write fails
-        with open(str(tcl_script), 'w+') as _fo:
-            _fo.write(''.join(to_write))
+
+    with open(tcl_script, 'w+') as _fo:
+        _fo.write(''.join(to_write))
     return tcl_script
 
 

From 249c994311d5a682c83245afd983aca5ca9e8982 Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Sun, 21 Jun 2026 17:29:43 -0500
Subject: [PATCH 12/13] clean up write_tcl_script()

---
 frustratometer/frustration/frustration.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/frustratometer/frustration/frustration.py b/frustratometer/frustration/frustration.py
index fefdd47..c1b1c66 100644
--- a/frustratometer/frustration/frustration.py
+++ b/frustratometer/frustration/frustration.py
@@ -936,11 +936,9 @@ def write_tcl_script(
         tcl script file
     """
     to_write = []
-    #fo = open(tcl_script, 'w+')
     single_frustration = np.nan_to_num(single_frustration,nan=0,posinf=0,neginf=0)
     pair_frustration = np.nan_to_num(pair_frustration,nan=0,posinf=0,neginf=0)
     
-    
     structure = prody.parsePDB(str(pdb_file))
     selection = structure.select('protein', chain=chain)
     residues = np.unique(selection.getResindices())

From 11bb77d52ae542151c3005c3333dabd76fd2baaf Mon Sep 17 00:00:00 2001
From: Finley Clark <60708336+ftclark3@users.noreply.github.com>
Date: Sun, 21 Jun 2026 17:31:01 -0500
Subject: [PATCH 13/13] update frustration.write_tcl_script() to draw lines
 very quickly and anchor them to CA positions rather than fixed coordinates;
 source the script in the tk console following alignment to move the lines
 along with the molecule

---
 frustratometer/frustration/frustration.py | 66 ++++++++++++++++-------
 1 file changed, 48 insertions(+), 18 deletions(-)

diff --git a/frustratometer/frustration/frustration.py b/frustratometer/frustration/frustration.py
index c1b1c66..92a8ffe 100644
--- a/frustratometer/frustration/frustration.py
+++ b/frustratometer/frustration/frustration.py
@@ -950,8 +950,10 @@ def write_tcl_script(
     cb_sel = selection.select('name CB or (resname GLY and name CA)')
     ca_xyz = np.full((len(residues), 3), np.nan)
     cb_xyz = np.full((len(residues), 3), np.nan)
+    ca_indices = np.full(len(residues), -1, dtype=int)
     ca_xyz[np.searchsorted(residues, ca_sel.getResindices())] = ca_sel.getCoords()
     cb_xyz[np.searchsorted(residues, cb_sel.getResindices())] = cb_sel.getCoords()
+    ca_indices[np.searchsorted(residues, ca_sel.getResindices())] = ca_sel.getIndices()
 
     to_write.append(f'[atomselect top all] set beta 0\n')
     # Single residue frustration
@@ -1001,8 +1003,24 @@ def process_dense_arrays():
         frust = pair_frustration[ii,jj] # pair_frustration is always dense
         dist = distance_matrix.lookup(ii,jj)
 
-    # Green = minimally frustrated (most negative first); red = frustrated (most positive
-    # first). For each: sort, cap at max_connections, then draw the in-range contacts.
+    # --- START OF HIGH-PERFORMANCE REDRAW Tcl PROCEDURE ---
+    to_write.append('\nproc redraw_bound_line {} {\n')
+    to_write.append('    set molid top\n')
+    to_write.append('    graphics $molid delete all\n\n')
+    
+    # 1. Grab all CA atoms AT ONCE
+    to_write.append('    set sel [atomselect $molid "name CA"]\n')
+    to_write.append('    set atom_data [$sel get {index x y z}]\n')
+    to_write.append('    $sel delete\n\n')
+
+    # 2. Build a Tcl array mapping index to {x y z} for O(1) lookups
+    to_write.append('    array unset coords\n')
+    to_write.append('    foreach record $atom_data {\n')
+    to_write.append('        lassign $record idx x y z\n')
+    to_write.append('        set coords($idx) [list $x $y $z]\n')
+    to_write.append('    }\n\n')
+
+    # Green = minimally frustrated; red = frustrated. Sort, cap, draw.
     for color, group, reverse in (
         ('green', frust < -0.78, False),
         ('red',   frust > 1,     True),
@@ -1016,28 +1034,40 @@ def process_dense_arrays():
             gi, gj, gd = gi[:max_connections], gj[:max_connections], gd[:max_connections]
 
         draw = (gd >= 3.5) & (gd <= 9.5)
-        draw &= np.abs(residues[gi] - residues[gj]) != 1 # Skip adjacent residues, which are often connected by a covalent bond and thus always close in space, but not necessarily frustrated.
+        draw &= np.abs(residues[gi] - residues[gj]) != 1 
         draw &= ~(np.isnan(ca_xyz[gi]).any(1) | np.isnan(ca_xyz[gj]).any(1) |
                   np.isnan(cb_xyz[gi]).any(1) | np.isnan(cb_xyz[gj]).any(1))
+        # Ensure mapping to valid ProDy/VMD indices
+        draw &= (ca_indices[gi] != -1) & (ca_indices[gj] != -1)
         gi, gj = gi[draw], gj[draw]
 
-        cb_dist = np.linalg.norm(cb_xyz[gi] - cb_xyz[gj], axis=1)   # Cbeta-Cbeta -> dashing
+        cb_dist = np.linalg.norm(cb_xyz[gi] - cb_xyz[gj], axis=1)   
         styles = np.where((cb_dist >= 3.5) & (cb_dist <= 6.5), 'solid', 'dashed')
 
-        to_write.append(f'draw color {color}\n')
-        for a, b, style in zip(gi, gj, styles):
-            p1, p2 = ca_xyz[a], ca_xyz[b]
-            to_write.append("draw line {%.3f %.3f %.3f} {%.3f %.3f %.3f} style %s width 2\n"
-                     % (p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], style))
-
-        to_write.append('''mol delrep top 0
-            mol color Beta
-            mol representation NewCartoon 0.300000 10.000000 4.100000 0
-            mol selection all
-            mol material Opaque
-            mol addrep top
-            color scale method GWR
-            ''')
+        # 3. Draw lines instantly using the array
+        if len(gi) > 0:
+            to_write.append(f'    graphics $molid color {color}\n')
+            for a, b, style in zip(gi, gj, styles):
+                id1, id2 = ca_indices[a], ca_indices[b]
+                to_write.append(f'    graphics $molid line $coords({id1}) $coords({id2}) style {style} width 2\n')
+            to_write.append('\n')
+
+    to_write.append('    puts "Lines redrawn successfully."\n')
+    to_write.append('}\n\n')
+    
+    # Execute the procedure once initially to render the graphic
+    to_write.append('redraw_bound_line\n\n')
+    # --- END OF Tcl PROCEDURE ---
+
+    # Visual representation setup
+    to_write.append('''mol delrep top 0
+        mol color Beta
+        mol representation NewCartoon 0.300000 10.000000 4.100000 0
+        mol selection all
+        mol material Opaque
+        mol addrep top
+        color scale method GWR
+        ''')
 
     if movie_name:
         to_write.append('''axes location Off