Repository: cctbx/cctbx_project (language: python)

Issue:
Explicit template instantiation of smtbx' geometrical_hydrogen_sites<N>
Anyone see any problems/have any objection to this? It's fixing a problem that isn't currently exhibited by the *standard* build, but I believe is technically is more standards-correct.

Effectively, it fixes a problem compiling `libsmtbx_refinement_constraints` as a shared library
on macOS (which is currently explicitly disabled for default builds).

The Problem
-------------
  This template class is defined privately in `libsmtbx_refinement_constraints`. On linux systems the implicit instantiation via usage in the rest of the library results in a weak-symbol export that can be linked to (by `smtbx_refinement_constraints_ext`), but it looks like the macOS linker is stricter about symbol visibility and exporting - the symbol ends up in a private (non-external) ELF section.

The Cause of Problem:
-------------
Relying on linking to visible implicit declarations seems to be, at best, implementation-specific behaviour. See standards:

>   14 Templates
>   8 ... A ... template must be defined in every translation unit in which it is implicitly instantiated (14.7.1), unless the corresponding specialization is explicitly instantiated (14.7.2) in some translation unit; no diagnostic is required.

The Solution
-------------
  Explicitly declare the template instantiations
Registry: pass along objects from one understand method to the next
Currently if one understand method from a parent Format class wishes to pass information to children Format classes, there is no way to do this.  For example, if FormatCBF's understand method reads the header and returns True, then each of the children classes (mini, full, etc) re-opens the file and re-reads the header to see if it can understand it.  Instead, if the understand method returned a tuple:

True, anyobject

or

False, None

Then the registry could pass on anyobject to the children.  Perhaps anyobject is a header dictionary from FormatSMV, or a detector address interpreted by FormatXTC.  It's implementation specific.

This would reduce file reads and allow communication of derived information to children classes.

It could be done with a minimum amount of intrusion by updating registry.py to look at what the understand method returned.  If it's a bool, just call the children's understand methods.  If it's a tuple, pass along the second object too.

Thoughts?
cluster.unit_cell stdevs bug
It looks like the unit cell standard deviations are calculated using the input unit cell parameters, rather than after being mapped to the niggli cell, as demonstrated by the following example.
 Note the different ordering of the unit cell parameters on the last line - all unit cells appear to be mapped to the correct niggli cell.

```
$ cat symmetry.txt
37.96099853515625 79.125 79.125 90.0 90.0 90.0 P1
37.96200180053711 79.1259994506836 79.1259994506836 90.0 90.0 90.0 P1
37.9630012512207 79.12000274658203 79.12000274658203 90.0 90.0 90.0 P1
37.95899963378906 79.12348359159503 79.12348359159503 90.01280173257825 90.0 90.0 P1
37.941001892089844 79.10900115966797 79.10900115966797 90.0 90.0 90.0 P1
37.957000732421875 79.12699890136719 79.12699890136719 90.0 90.0 90.0 P1
37.96500015258789 79.14099884033203 79.14099884033203 90.0 90.0 90.0 P1
37.93000030517578 79.13600158691406 79.13600158691406 90.0 90.0 90.0 P1
37.95500183105469 79.12899780273438 79.12899780273438 90.0 90.0 90.0 P1
79.02999878  79.02700043 37.91400146  90.          90.01300049  90. P1

$ cluster.unit_cell symmetry.txt --text --noplot
There are 10 lines in the input file
Unit cell: (37.961, 79.125, 79.125, 90, 90, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.961, 79.125, 79.125, 90, 90, 90)
   niggli-->Space group: P 1 (No. 1)
Unit cell: (37.962, 79.126, 79.126, 90, 90, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.962, 79.126, 79.126, 90, 90, 90)
   niggli-->Space group: P 1 (No. 1)
Unit cell: (37.963, 79.12, 79.12, 90, 90, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.963, 79.12, 79.12, 90, 90, 90)
   niggli-->Space group: P 1 (No. 1)
Unit cell: (37.959, 79.1235, 79.1235, 90.0128, 90, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.959, 79.1235, 79.1235, 90.0128, 90, 90)
   niggli-->Space group: P 1 (No. 1)
Unit cell: (37.941, 79.109, 79.109, 90, 90, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.941, 79.109, 79.109, 90, 90, 90)
   niggli-->Space group: P 1 (No. 1)
Unit cell: (37.957, 79.127, 79.127, 90, 90, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.957, 79.127, 79.127, 90, 90, 90)
   niggli-->Space group: P 1 (No. 1)
Unit cell: (37.965, 79.141, 79.141, 90, 90, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.965, 79.141, 79.141, 90, 90, 90)
   niggli-->Space group: P 1 (No. 1)
Unit cell: (37.93, 79.136, 79.136, 90, 90, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.93, 79.136, 79.136, 90, 90, 90)
   niggli-->Space group: P 1 (No. 1)
Unit cell: (37.955, 79.129, 79.129, 90, 90, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.955, 79.129, 79.129, 90, 90, 90)
   niggli-->Space group: P 1 (No. 1)
Unit cell: (79.03, 79.027, 37.914, 90, 90.013, 90)
Space group: P 1 (No. 1)
   niggli-->Unit cell: (37.914, 79.027, 79.03, 90, 90.013, 90)
   niggli-->Space group: P 1 (No. 1)
10 lattices will be analyzed
INFO cluster.ab_cluster: Hierarchical clustering of unit cells
INFO cluster.ab_cluster: Using Andrews-Bernstein distance from Andrews & Bernstein J Appl Cryst 47:346 (2014)
INFO cluster.ab_cluster: Distances have been calculated

0 singletons:

Point group    a           b           c          alpha        beta         gamma       

1 clusters:

Cluster_id       N_xtals  Med_a         Med_b         Med_c         Med_alpha    Med_beta     Med_gamma   
cluster_1        10       37.96 (12.32) 79.13 (0.03 ) 79.13 (12.36) 90.00 (0.00) 90.00 (0.00) 90.00 (0.00)
10 in P1.
Standard deviations are in brackets.
```


Discussion / hints:
Since nobody seems to have any objections...

Got it.  Yet another bug to work out.  Will try to do something sensible
this weekend.

Nick

Nicholas K. Sauter, Ph. D.
Senior Scientist, Molecular Biophysics & Integrated Bioimaging Division
Lawrence Berkeley National Laboratory
1 Cyclotron Rd., Bldg. 33R0345
Berkeley, CA 94720
(510) 486-5713

On Fri, Nov 17, 2017 at 2:21 PM, Richard Gildea <notifications@github.com>
wrote:

> It looks like the unit cell standard deviations are calculated using the
> input unit cell parameters, rather than after being mapped to unit cell
> parameters, as demonstrated by the following example.
> Note the different ordering of the unit cell parameters on the last line -
> all unit cells appear to be mapped to the correct niggli cell.
>
> $ cat symmetry.txt
> 37.96099853515625 79.125 79.125 90.0 90.0 90.0 P1
> 37.96200180053711 79.1259994506836 79.1259994506836 90.0 90.0 90.0 P1
> 37.9630012512207 79.12000274658203 79.12000274658203 90.0 90.0 90.0 P1
> 37.95899963378906 79.12348359159503 79.12348359159503 90.01280173257825 90.0 90.0 P1
> 37.941001892089844 79.10900115966797 79.10900115966797 90.0 90.0 90.0 P1
> 37.957000732421875 79.12699890136719 79.12699890136719 90.0 90.0 90.0 P1
> 37.96500015258789 79.14099884033203 79.14099884033203 90.0 90.0 90.0 P1
> 37.93000030517578 79.13600158691406 79.13600158691406 90.0 90.0 90.0 P1
> 37.95500183105469 79.12899780273438 79.12899780273438 90.0 90.0 90.0 P1
> 79.02999878  79.02700043 37.91400146  90.          90.01300049  90. P1
>
> $ cluster.unit_cell symmetry.txt --text --noplot
> There are 10 lines in the input file
> Unit cell: (37.961, 79.125, 79.125, 90, 90, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.961, 79.125, 79.125, 90, 90, 90)
>    niggli-->Space group: P 1 (No. 1)
> Unit cell: (37.962, 79.126, 79.126, 90, 90, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.962, 79.126, 79.126, 90, 90, 90)
>    niggli-->Space group: P 1 (No. 1)
> Unit cell: (37.963, 79.12, 79.12, 90, 90, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.963, 79.12, 79.12, 90, 90, 90)
>    niggli-->Space group: P 1 (No. 1)
> Unit cell: (37.959, 79.1235, 79.1235, 90.0128, 90, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.959, 79.1235, 79.1235, 90.0128, 90, 90)
>    niggli-->Space group: P 1 (No. 1)
> Unit cell: (37.941, 79.109, 79.109, 90, 90, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.941, 79.109, 79.109, 90, 90, 90)
>    niggli-->Space group: P 1 (No. 1)
> Unit cell: (37.957, 79.127, 79.127, 90, 90, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.957, 79.127, 79.127, 90, 90, 90)
>    niggli-->Space group: P 1 (No. 1)
> Unit cell: (37.965, 79.141, 79.141, 90, 90, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.965, 79.141, 79.141, 90, 90, 90)
>    niggli-->Space group: P 1 (No. 1)
> Unit cell: (37.93, 79.136, 79.136, 90, 90, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.93, 79.136, 79.136, 90, 90, 90)
>    niggli-->Space group: P 1 (No. 1)
> Unit cell: (37.955, 79.129, 79.129, 90, 90, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.955, 79.129, 79.129, 90, 90, 90)
>    niggli-->Space group: P 1 (No. 1)
> Unit cell: (79.03, 79.027, 37.914, 90, 90.013, 90)
> Space group: P 1 (No. 1)
>    niggli-->Unit cell: (37.914, 79.027, 79.03, 90, 90.013, 90)
>    niggli-->Space group: P 1 (No. 1)
> 10 lattices will be analyzed
> INFO cluster.ab_cluster: Hierarchical clustering of unit cells
> INFO cluster.ab_cluster: Using Andrews-Bernstein distance from Andrews & Bernstein J Appl Cryst 47:346 (2014)
> INFO cluster.ab_cluster: Distances have been calculated
>
> 0 singletons:
>
> Point group    a           b           c          alpha        beta         gamma
>
> 1 clusters:
>
> Cluster_id       N_xtals  Med_a         Med_b         Med_c         Med_alpha    Med_beta     Med_gamma
> cluster_1        10       37.96 (12.32) 79.13 (0.03 ) 79.13 (12.36) 90.00 (0.00) 90.00 (0.00) 90.00 (0.00)
> 10 in P1.
> Standard deviations are in brackets.
>
> —
> You are receiving this because you are subscribed to this thread.
> Reply to this email directly, view it on GitHub
> <https://github.com/cctbx/cctbx_project/issues/97>, or mute the thread
> <https://github.com/notifications/unsubscribe-auth/AOGuVsScLa4QjNc_2D4xNZXHsxeNUbHFks5s3gbvgaJpZM4Qisv1>
> .
>

Richard,
I've committed a reasonably sensible fix.  It won't work for corner cases
where neighboring unit cells cross polytope boundaries, but for the simple
tetragonal input (lysozyme) there are two lines of output:
1) mean and std dev cells, calculated in Niggli setting.
2) Median cell back-converted to the input setting, in this case tetragonal.

1 clusters:


Cluster_id       N_xtals  Med_a         Med_b         Med_c
Med_alpha    Med_beta     Med_gamma   Delta(deg)

10 in P1.

cluster_1        10       37.96 (0.02 ) 79.13 (0.03 ) 79.13 (0.03 ) 90.00
(0.00) 90.00 (0.00) 90.00 (0.00)

     P 4/m m m (No. 123)  79.13         79.13         37.96         90.00
      90.00        90.00         0.0



Standard deviations are in brackets.

Each cluster:

Input lattice count, with integration Bravais setting space group.

Cluster median with Niggli cell parameters (std dev in brackets).

Highest possible metric symmetry and unit cell using LePage (J Appl Cryst
1982, 15:255) method, maximum delta 3deg.

Nick

Nicholas K. Sauter, Ph. D.
Senior Scientist, Molecular Biophysics & Integrated Bioimaging Division
Lawrence Berkeley National Laboratory
1 Cyclotron Rd., Bldg. 33R0345
Berkeley, CA 94720
(510) 486-5713

On Fri, Nov 17, 2017 at 3:21 PM, Nicholas Sauter <nksauter@lbl.gov> wrote:

> Got it.  Yet another bug to work out.  Will try to do something sensible
> this weekend.
>
> Nick
>
> Nicholas K. Sauter, Ph. D.
> Senior Scientist, Molecular Biophysics & Integrated Bioimaging Division
> Lawrence Berkeley National Laboratory
> 1 Cyclotron Rd., Bldg. 33R0345
> Berkeley, CA 94720
> (510) 486-5713
>
> On Fri, Nov 17, 2017 at 2:21 PM, Richard Gildea <notifications@github.com>
> wrote:
>
>> It looks like the unit cell standard deviations are calculated using the
>> input unit cell parameters, rather than after being mapped to unit cell
>> parameters, as demonstrated by the following example.
>> Note the different ordering of the unit cell parameters on the last line
>> - all unit cells appear to be mapped to the correct niggli cell.
>>
>> $ cat symmetry.txt
>> 37.96099853515625 79.125 79.125 90.0 90.0 90.0 P1
>> 37.96200180053711 79.1259994506836 79.1259994506836 90.0 90.0 90.0 P1
>> 37.9630012512207 79.12000274658203 79.12000274658203 90.0 90.0 90.0 P1
>> 37.95899963378906 79.12348359159503 79.12348359159503 90.01280173257825 90.0 90.0 P1
>> 37.941001892089844 79.10900115966797 79.10900115966797 90.0 90.0 90.0 P1
>> 37.957000732421875 79.12699890136719 79.12699890136719 90.0 90.0 90.0 P1
>> 37.96500015258789 79.14099884033203 79.14099884033203 90.0 90.0 90.0 P1
>> 37.93000030517578 79.13600158691406 79.13600158691406 90.0 90.0 90.0 P1
>> 37.95500183105469 79.12899780273438 79.12899780273438 90.0 90.0 90.0 P1
>> 79.02999878  79.02700043 37.91400146  90.          90.01300049  90. P1
>>
>> $ cluster.unit_cell symmetry.txt --text --noplot
>> There are 10 lines in the input file
>> Unit cell: (37.961, 79.125, 79.125, 90, 90, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.961, 79.125, 79.125, 90, 90, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> Unit cell: (37.962, 79.126, 79.126, 90, 90, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.962, 79.126, 79.126, 90, 90, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> Unit cell: (37.963, 79.12, 79.12, 90, 90, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.963, 79.12, 79.12, 90, 90, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> Unit cell: (37.959, 79.1235, 79.1235, 90.0128, 90, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.959, 79.1235, 79.1235, 90.0128, 90, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> Unit cell: (37.941, 79.109, 79.109, 90, 90, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.941, 79.109, 79.109, 90, 90, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> Unit cell: (37.957, 79.127, 79.127, 90, 90, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.957, 79.127, 79.127, 90, 90, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> Unit cell: (37.965, 79.141, 79.141, 90, 90, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.965, 79.141, 79.141, 90, 90, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> Unit cell: (37.93, 79.136, 79.136, 90, 90, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.93, 79.136, 79.136, 90, 90, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> Unit cell: (37.955, 79.129, 79.129, 90, 90, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.955, 79.129, 79.129, 90, 90, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> Unit cell: (79.03, 79.027, 37.914, 90, 90.013, 90)
>> Space group: P 1 (No. 1)
>>    niggli-->Unit cell: (37.914, 79.027, 79.03, 90, 90.013, 90)
>>    niggli-->Space group: P 1 (No. 1)
>> 10 lattices will be analyzed
>> INFO cluster.ab_cluster: Hierarchical clustering of unit cells
>> INFO cluster.ab_cluster: Using Andrews-Bernstein distance from Andrews & Bernstein J Appl Cryst 47:346 (2014)
>> INFO cluster.ab_cluster: Distances have been calculated
>>
>> 0 singletons:
>>
>> Point group    a           b           c          alpha        beta         gamma
>>
>> 1 clusters:
>>
>> Cluster_id       N_xtals  Med_a         Med_b         Med_c         Med_alpha    Med_beta     Med_gamma
>> cluster_1        10       37.96 (12.32) 79.13 (0.03 ) 79.13 (12.36) 90.00 (0.00) 90.00 (0.00) 90.00 (0.00)
>> 10 in P1.
>> Standard deviations are in brackets.
>>
>> —
>> You are receiving this because you are subscribed to this thread.
>> Reply to this email directly, view it on GitHub
>> <https://github.com/cctbx/cctbx_project/issues/97>, or mute the thread
>> <https://github.com/notifications/unsubscribe-auth/AOGuVsScLa4QjNc_2D4xNZXHsxeNUbHFks5s3gbvgaJpZM4Qisv1>
>> .
>>
>
>

Can this be closed yet? Thanks.

Fix patch:
diff --git a/.gitignore b/.gitignore
index 29990d223a..014d44be79 100644
--- a/.gitignore
+++ b/.gitignore
@@ -65,3 +65,8 @@ scitbx/array_family/detail/operator_functors.h
 scitbx/array_family/detail/std_imports.h
 sphinx/build
 xfel/clustering/test_cluster.py
+mmtbx/suitename/test/
+mmtbx/suitename/suitename.stderr.log
+mmtbx/suitename/assist/
+mmtbx/suitename/*.show.txt
+*.txt
diff --git a/libtbx/auto_build/bootstrap.py b/libtbx/auto_build/bootstrap.py
index d277c39fb6..13e5fb5fae 100644
--- a/libtbx/auto_build/bootstrap.py
+++ b/libtbx/auto_build/bootstrap.py
@@ -1001,10 +1001,6 @@ class probe_module(SourceModule):
   module = 'probe'
   anonymous = ['svn', 'https://github.com/rlabduke/probe.git/trunk']
 
-class suitename_module(SourceModule):
-  module = 'suitename'
-  anonymous = ['svn', 'https://github.com/rlabduke/suitename.git/trunk']
-
 class reduce_module(SourceModule):
   module = 'reduce'
   anonymous = ['git', '-b v3.7.2',
@@ -1483,7 +1479,7 @@ def run(self):
 
   def _add_svn(self, module, url):
     update_list = ['update']
-    if module in ["reduce", "probe", "king", "suitename"]:
+    if module in ["reduce", "probe", "king"]:
       pass
     elif self.revert:
       update_list = ['update', '-r', self.revert]
@@ -1969,8 +1965,7 @@ class MOLPROBITYBuilder(Builder):
     'molprobity',
     #'chem_data', #chem_data removed from molprobity builder until accessible outside cci, -CJW
     'reduce',
-    'probe',
-    'suitename'
+    'probe'
   ]
   # Copy these sources from cci.lbl.gov
   HOT = [
@@ -2314,7 +2309,6 @@ class PhenixBuilder(CCIBuilder):
     'reduce',
     'probe',
     'king',
-    'suitename',
     'dials',
     'xia2',
     'phaser',
diff --git a/mmtbx/command_line/suitename.py b/mmtbx/command_line/suitename.py
new file mode 100644
index 0000000000..00d5095195
--- /dev/null
+++ b/mmtbx/command_line/suitename.py
@@ -0,0 +1,11 @@
+# LIBTBX_SET_DISPATCHER_NAME phenix.suitename
+# LIBTBX_SET_DISPATCHER_NAME molprobity.suitename
+# LIBTBX_SET_DISPATCHER_NAME cctbx.suitename
+
+from __future__ import absolute_import, division, print_function
+import sys
+from  mmtbx.programs import suitename
+
+if __name__ == '__main__':
+  suitename.run(args=sys.argv[1:])
+
diff --git a/mmtbx/programs/suitename.py b/mmtbx/programs/suitename.py
new file mode 100644
index 0000000000..be397385f4
--- /dev/null
+++ b/mmtbx/programs/suitename.py
@@ -0,0 +1,236 @@
+#        Copyright 2021  Richardson Lab at Duke University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import division
+from __future__ import nested_scopes, generators, absolute_import
+from __future__ import with_statement, print_function
+import sys, os
+
+from mmtbx.suitename import dualparse, suites
+from mmtbx.suitename.suitename import main, version
+
+from iotbx.cli_parser import CCTBXParser
+from libtbx.program_template import ProgramTemplate
+from libtbx.utils import multi_out, show_total_time
+
+#import libtbx.load_env
+#from libtbx.utils import Usage
+
+import os, sys
+
+def run(args):
+  "The main program, if run from CCTBX / PHENIX."
+  logger = multi_out()
+  logger.register('stderr', sys.stderr)
+  logger2 = multi_out()
+  logger2.register('stdout', sys.stdout)
+  logger3 = multi_out()
+  # TO DIAGNOSE OPTIONS TROUBLES:
+  # logger3.register('verbiage', open("suitename.stderr.log", "w"))
+
+  parser = dualparse.parseArgs(Program, logger3)
+  working_phil = parser.working_phil
+  options = working_phil.extract().suitename
+
+  # now we call into the core of suitename itself
+  if options.version:
+      print(version, file=logger2)
+      return
+  if options.infile == "" or options.residuein or options.suitein:
+      # let the core figure out the input
+      main(optionsIn=options, outFile=logger2, errorFile=logger)
+  else:
+      type, ext = analyzeFileType(options.infile)
+      if type=="":
+        logger.write("File extension "+str(ext)+" not recognized\n")
+        return
+      if type == "pdb":
+          suites.main(options=options, outFile=logger2, errorFile=logger)
+      else:
+        # help the core figure out the input file type
+        if type == "kinemage":
+          options.suitein = True
+        elif type == "dangle":
+          options.residuein = True
+        main(optionsIn=options, outFile=logger2, errorFile=logger)
+
+
+extensionList={
+    "pdb": "pdb",
+    "cif": "pdb",
+    "kin": "kinemage",
+    "dangle": "dangle",
+    "suitegeom": "dangle",
+}
+
+def analyzeFileType(filename):
+    base, extension = os.path.splitext(filename)
+    extension = extension[1:]
+    type = extensionList.get(extension, "")
+    return type, extension
+
+
+class Program(ProgramTemplate):
+  prog = os.getenv('LIBTBX_DISPATCHER_NAME')
+  program_name = "suitename"
+  description="""
+   < insert help text here>
+""" % locals()
+
+  # The substructure below has been commented out because it would make it more
+  # difficult to integrate with suitename's ability to run independently
+  # of CCTBX
+  master_phil_str = """
+    suitename {
+      # input
+        infile=""
+          .type=str
+          .help="the file to process"
+        anglefields = 9
+          .type=int
+          .help="number of angle fields provided, for textual input only"
+        pointidfields = 7
+          .type=int
+          .help="number of point id fields before the angle fields"
+        ptid=0
+          .type=int
+          .help="number of point id fields before the angle fields"
+        residuein=false
+          .type=bool
+          .help="expect dangle format giving residues"
+        suitein=false
+          .type=bool
+          .help="expect kinemage format giving suites directly"
+      # output
+        string=False
+          .type=bool
+          .help="output in string format, 3 characters per suite"
+        kinemage=False
+          .type=bool
+          .help="output in kinemage format, useful for visualization"
+        report=true
+          .type=bool
+          .help="output as a report, giving statistical details"
+        chart=False
+          .type=bool
+          .help="modifier to standard report, output without statistical summary"
+        nosequence = False
+          .type=bool
+          .help="modifier to string format, do not include base letters"
+        causes=False
+          .type=bool
+          .help="output extra details concerning the causes of each assignment made"
+        test=False
+          .type=bool
+          .help="display a lat of additional information about program internals"
+      # compute
+        satellites=False
+          .type=bool
+          .help="use the special satelliteWidths values for satellites"
+        nowannabe=False
+          .type=bool
+          .help="do not consider 'wannabe' clusters"
+        noinc=False
+          .type=bool
+          .help="do not display incomplete suites"
+        etatheta=False
+          .type=bool
+        altid="A"
+          .type=str
+          .help="which alternate conformer to use (A, B, etc)"
+        altidfield = 6
+          .type=int
+          .help="which field (1-based) gives the alternate conformer code"
+        version=False
+          .type=bool
+          .help="give the version number of suite name"
+      # deprecated
+        oneline=False
+          .type=bool
+      }
+"""
+
+# might add:
+        # altidval="A"
+        #   .type=str
+        #   .help="which alternate conformer to use (A, B, etc)"
+#
+  datatypes = ['model', 'phil']  # also
+  data_manager_options = ['model_skip_expand_with_mtrix']
+  known_article_ids = ['molprobity']
+
+  def validate(self):
+    pass
+
+  def run(self, args):
+    pass
+
+
+  # end of class Program
+
+
+# def rest_of_old_run(args):  -- an intermediate state, outdated technique
+#   namespace, others = parser.parse_args(sys.argv[1:])
+#   # whatever the old fashioned parser won't use becomes part of <others>
+#   # and will be given to the phil parser
+#
+#   parser = CCTBXParser(
+#     program_class=Program,
+#     logger=logger)
+#   namespace = parser.parse_args(others)
+#
+#   # start program
+#   print('Starting job', file=logger)
+#   print('='*79, file=logger)
+#   phil1 = parser.working_phil.extract()
+#   args2 = parseCommandLine()
+#   task = Program(
+#     parser.data_manager, phil1, logger=logger2)
+#   main()
+
+#=============================================================================
+# def run(args):
+#
+#   # create parser
+#   logger = multi_out()
+#   logger.register('stderr', sys.stderr)
+#   logger2 = multi_out()
+#   logger2.register('stdout', sys.stdout)
+#
+#   parser = CCTBXParser(
+#     program_class=cablam.Program,
+#     logger=logger)
+#   namespace = parser.parse_args(sys.argv[1:])
+#
+#   # start program
+#   print('Starting job', file=logger)
+#   print('='*79, file=logger)
+#   task = cablam.Program(
+#     parser.data_manager, parser.working_phil.extract(), logger=logger2)
+#
+#   # validate inputs
+#   task.validate()
+#
+#   # run program
+#   task.run()
+#
+#   # stop timer
+#   print('', file=logger)
+#   print('='*79, file=logger)
+#   print('Job complete', file=logger)
+#   show_total_time(out=logger)
+
+# =============================================================================
+if __name__ == '__main__':
+  run(sys.argv[1:])
diff --git a/mmtbx/suitename/LICENSE.txt b/mmtbx/suitename/LICENSE.txt
new file mode 100644
index 0000000000..0074ce01f1
--- /dev/null
+++ b/mmtbx/suitename/LICENSE.txt
@@ -0,0 +1,54 @@
+Apache License
+
+Version 2.0, January 2004
+
+http://www.apache.org/licenses/
+
+TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION
+
+1. Definitions.
+
+"License" shall mean the terms and conditions for use, reproduction, and distribution as defined by Sections 1 through 9 of this document.
+
+"Licensor" shall mean the copyright owner or entity authorized by the copyright owner that is granting the License.
+
+"Legal Entity" shall mean the union of the acting entity and all other entities that control, are controlled by, or are under common control with that entity. For the purposes of this definition, "control" means (i) the power, direct or indirect, to cause the direction or management of such entity, whether by contract or otherwise, or (ii) ownership of fifty percent (50%) or more of the outstanding shares, or (iii) beneficial ownership of such entity.
+
+"You" (or "Your") shall mean an individual or Legal Entity exercising permissions granted by this License.
+
+"Source" form shall mean the preferred form for making modifications, including but not limited to software source code, documentation source, and configuration files.
+
+"Object" form shall mean any form resulting from mechanical transformation or translation of a Source form, including but not limited to compiled object code, generated documentation, and conversions to other media types.
+
+"Work" shall mean the work of authorship, whether in Source or Object form, made available under the License, as indicated by a copyright notice that is included in or attached to the work (an example is provided in the Appendix below).
+
+"Derivative Works" shall mean any work, whether in Source or Object form, that is based on (or derived from) the Work and for which the editorial revisions, annotations, elaborations, or other modifications represent, as a whole, an original work of authorship. For the purposes of this License, Derivative Works shall not include works that remain separable from, or merely link (or bind by name) to the interfaces of, the Work and Derivative Works thereof.
+
+"Contribution" shall mean any work of authorship, including the original version of the Work and any modifications or additions to that Work or Derivative Works thereof, that is intentionally submitted to Licensor for inclusion in the Work by the copyright owner or by an individual or Legal Entity authorized to submit on behalf of the copyright owner. For the purposes of this definition, "submitted" means any form of electronic, verbal, or written communication sent to the Licensor or its representatives, including but not limited to communication on electronic mailing lists, source code control systems, and issue tracking systems that are managed by, or on behalf of, the Licensor for the purpose of discussing and improving the Work, but excluding communication that is conspicuously marked or otherwise designated in writing by the copyright owner as "Not a Contribution."
+
+"Contributor" shall mean Licensor and any individual or Legal Entity on behalf of whom a Contribution has been received by Licensor and subsequently incorporated within the Work.
+
+2. Grant of Copyright License. Subject to the terms and conditions of this License, each Contributor hereby grants to You a perpetual, worldwide, non-exclusive, no-charge, royalty-free, irrevocable copyright license to reproduce, prepare Derivative Works of, publicly display, publicly perform, sublicense, and distribute the Work and such Derivative Works in Source or Object form.
+
+3. Grant of Patent License. Subject to the terms and conditions of this License, each Contributor hereby grants to You a perpetual, worldwide, non-exclusive, no-charge, royalty-free, irrevocable (except as stated in this section) patent license to make, have made, use, offer to sell, sell, import, and otherwise transfer the Work, where such license applies only to those patent claims licensable by such Contributor that are necessarily infringed by their Contribution(s) alone or by combination of their Contribution(s) with the Work to which such Contribution(s) was submitted. If You institute patent litigation against any entity (including a cross-claim or counterclaim in a lawsuit) alleging that the Work or a Contribution incorporated within the Work constitutes direct or contributory patent infringement, then any patent licenses granted to You under this License for that Work shall terminate as of the date such litigation is filed.
+
+4. Redistribution. You may reproduce and distribute copies of the Work or Derivative Works thereof in any medium, with or without modifications, and in Source or Object form, provided that You meet the following conditions:
+
+You must give any other recipients of the Work or Derivative Works a copy of this License; and
+You must cause any modified files to carry prominent notices stating that You changed the files; and
+You must retain, in the Source form of any Derivative Works that You distribute, all copyright, patent, trademark, and attribution notices from the Source form of the Work, excluding those notices that do not pertain to any part of the Derivative Works; and
+If the Work includes a "NOTICE" text file as part of its distribution, then any Derivative Works that You distribute must include a readable copy of the attribution notices contained within such NOTICE file, excluding those notices that do not pertain to any part of the Derivative Works, in at least one of the following places: within a NOTICE text file distributed as part of the Derivative Works; within the Source form or documentation, if provided along with the Derivative Works; or, within a display generated by the Derivative Works, if and wherever such third-party notices normally appear. The contents of the NOTICE file are for informational purposes only and do not modify the License. You may add Your own attribution notices within Derivative Works that You distribute, alongside or as an addendum to the NOTICE text from the Work, provided that such additional attribution notices cannot be construed as modifying the License.
+
+You may add Your own copyright statement to Your modifications and may provide additional or different license terms and conditions for use, reproduction, or distribution of Your modifications, or for any such Derivative Works as a whole, provided Your use, reproduction, and distribution of the Work otherwise complies with the conditions stated in this License.
+5. Submission of Contributions. Unless You explicitly state otherwise, any Contribution intentionally submitted for inclusion in the Work by You to the Licensor shall be under the terms and conditions of this License, without any additional terms or conditions. Notwithstanding the above, nothing herein shall supersede or modify the terms of any separate license agreement you may have executed with Licensor regarding such Contributions.
+
+6. Trademarks. This License does not grant permission to use the trade names, trademarks, service marks, or product names of the Licensor, except as required for reasonable and customary use in describing the origin of the Work and reproducing the content of the NOTICE file.
+
+7. Disclaimer of Warranty. Unless required by applicable law or agreed to in writing, Licensor provides the Work (and each Contributor provides its Contributions) on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied, including, without limitation, any warranties or conditions of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A PARTICULAR PURPOSE. You are solely responsible for determining the appropriateness of using or redistributing the Work and assume any risks associated with Your exercise of permissions under this License.
+
+8. Limitation of Liability. In no event and under no legal theory, whether in tort (including negligence), contract, or otherwise, unless required by applicable law (such as deliberate and grossly negligent acts) or agreed to in writing, shall any Contributor be liable to You for damages, including any direct, indirect, special, incidental, or consequential damages of any character arising as a result of this License or out of the use or inability to use the Work (including but not limited to damages for loss of goodwill, work stoppage, computer failure or malfunction, or any and all other commercial damages or losses), even if such Contributor has been advised of the possibility of such damages.
+
+9. Accepting Warranty or Additional Liability. While redistributing the Work or Derivative Works thereof, You may choose to offer, and charge a fee for, acceptance of support, warranty, indemnity, or other liability obligations and/or rights consistent with this License. However, in accepting such obligations, You may act only on Your own behalf and on Your sole responsibility, not on behalf of any other Contributor, and only if You agree to indemnify, defend, and hold each Contributor harmless for any liability incurred by, or claims asserted against, such Contributor by reason of your accepting any such warranty or additional liability.
+
+END OF TERMS AND CONDITIONS
+
diff --git a/mmtbx/suitename/README.md b/mmtbx/suitename/README.md
new file mode 100644
index 0000000000..56f9c9f3a8
--- /dev/null
+++ b/mmtbx/suitename/README.md
@@ -0,0 +1,56 @@
+# suitename
+Suitename - a tool for classifying backbone "suite" conformations (linkages between ribose sugars) in an RNA molecule.
+
+## Requirements
+Suitename requires Python 2.7 or later and it requires numpy to have been installed.
+
+## Running
+
+Suitename is integrated into the CCTBX library and can be run from within it as `phenix.suitename`, `molprobity.suitename` or `cctbx.suitename`.
+
+Suitename can also be run from the command line using a Python 3 interpreter: `python3 suitename.py [arguments]`
+
+Suitename takes an input file as its first argument, or if none provided, reads from standard input. It writes its results to standard output. It classifies the conformation of each suite into one of several dozen predefined clusters that have been determined by years of study. 
+
+Two forms of input are supported on standard input:
+  1. A list of the 6 dihedral angles $/alpha$, $/beta$, $/gamma$, $/delta$, $/epsilon$, $/zeta$) for each residue of the RNA molecule. Suitename will re-parse the dihedral angles in the residues to obtain the 7 dihedral angles in each suite ($/delta/-1, $/epsilon/-1, $/zeta/-1, $/alpha$, $/beta$, $/gamma$, $/delta$), and then operate on those. This is the default input.
+
+    Each line of this format describes one residue, with fields separated by colons. The first several (default 6) are ID information, the remainder are angles. Sample:
+1: A:   7: : :  U:-75.533:-154.742:48.162:80.895:-148.423:-159.688
+    If the number of fields before the first angle is not 6, use
+    -- pointIDfields <n> 
+    to tell Suitename the true number.
+
+  2. A kinemage file providing a list of 7 or 9 dihedral angles in each suite. Mark this by using the --suitein command line flag. 9 angles include $/chi$-1 and $/chi$.
+
+A third form of input is supported when running phenix.suitename or molprobity.suitename and passing the file name on the command line: a molecule file in CIF or PDB format.  When run on this type of file, Suitename will compute the dihedral angles for all RNA bases and compute the resulting scores.
+
+Three forms of output are supported:
+  1. A text report, showing the classification of each suite into a cluster or outlier, and the neatness of fit ("suiteness") of the suite into that cluster. Suiteness is the cosine of the normalized distance of a suite datapoint from the power 3 hyperellipsoid boundary of the cluster in toward its center. A statistical summary is included at the end. This is the default output format.
+  2. A kinemage file, which will display a 7D data point for each suite in the data. Points are color-coded according to the clusters to which they have been assigned. Each cluster is displayed as a colored ring surrounding its defined center. Specify this by using the --kinemage command line flag.
+  3. A brief string showing only the cluster assignments. It consists of three characters per suite - base identity (uc) and 2-character number-letter name of the suite cluster (e.g., C1aG1gU1aA1aA1cG). Specify this by using the --string command line flag.
+
+Many other command line options are available; type `python3 suitename.py --help` to display them.
+
+## Use as a library
+
+See the suites.rtf document for more information about using suites as a library
+rather than a program.
+
+## Testing
+
+In addition to the tests found in the unit-test direcgtory, this Python code was run
+across all files in the PDB that have RNA and its output was
+compared to the original C code output for -report, -kinemage, and -string outputs.
+It was tested with both Python 2.7 and Python 3.6.
+It was tested on the dangle-format output of mmtbx.mp_geo and found to produce the same
+output modulo the following slight differences:
+- Python 3 uses Banker's rounding, so the low-order digit may be off by one.
+- Python performs floating-point calculations at double resolution, so some of
+the calculations are slightly different for numbers right near the cutoffs, sometimes
+making corner cases fall into a neighboring bin.
+
+## Directories
+
+**unit-test** Contains unit-test and regression test code for SuiteName.
+
diff --git a/mmtbx/suitename/__init__.py b/mmtbx/suitename/__init__.py
new file mode 100644
index 0000000000..e32cc32f4d
--- /dev/null
+++ b/mmtbx/suitename/__init__.py
@@ -0,0 +1,9 @@
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+import sys, os
+
+currentdir = os.path.dirname(os.path.abspath(__file__))
+sys.path.append(currentdir)
+# print("init:", sys.path)
+
+# Don't know why I should need this... but it helps
diff --git a/mmtbx/suitename/diangle.py b/mmtbx/suitename/diangle.py
new file mode 100644
index 0000000000..a8fde58fcf
--- /dev/null
+++ b/mmtbx/suitename/diangle.py
@@ -0,0 +1,361 @@
+#        Copyright 2021  Richardson Lab at Duke University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+import os, sys
+
+from iotbx.data_manager import DataManager    #   Load in the DataManager
+from scitbx.matrix import dihedral_angle
+from suitenamedefs import Residue, findBase
+
+
+n = "alpha beta gamma delta epsilon zeta"
+a = "O3' P O5' C5' C4' C3'"
+bone = (" O3'", " P  ", " O5'", " C5'", " C4'", " C3'")
+
+names = n.split(" ")
+bones = 2 * bone
+
+dihedrals = (
+  ("alpha", (" O3'", " P  ", " O5'", " C5'", )),
+  ("beta", (" P  ", " O5'", " C5'", " C4'", )),
+  ("gamma", (" O5'", " C5'", " C4'", " C3'", )),
+  ("delta", (" C5'", " C4'", " C3'", " O3'", )),
+  ("epsilon", (" C4'", " C3'", " O3'", " P  ", )),
+  ("zeta", (" C3'", " O3'", " P  ", " O5'", ))
+)
+
+
+class Failure(Exception):
+  message = ""
+
+  def __init__(self, msg):
+    self.message=msg
+
+
+def main(inFile):
+  # for testing only
+  alt = 'A'
+  mgr = loadModel(inFile)
+
+  residues = getResidueDihedrals(mgr, alt,
+                                 name=os.path.splitext(inFile)[0])
+  for r in residues:
+    #print(residueString(r))
+    pass
+
+
+def loadModel(filename):
+  dm = DataManager()             #   Initialize the DataManager and call it dm
+  dm.set_overwrite(True)         #   tell the DataManager to overwrite files with the same name
+  #print("Reading file")
+  manager = dm.get_model(filename)
+  return manager
+
+
+def getResidueDihedrals(mgr, altcode='', name='', errorFile = sys.stderr):
+  # use altcode='A', 'B', etc. for real alt work
+  # if name is provided produce a debug file <name>.suites
+  global manager, hierarchy, backbone_hierarchy, chain
+
+  if altcode=='':  altcode = 'A'
+  manager = mgr
+  hierarchy = manager.get_hierarchy()
+
+  # dboutput = open("hierarchy.show.txt", "w")
+  selector = "name P or name O5' or name C5' or name C4' or name C3' or name O3'"
+  selection = manager.selection(selector)
+  #print (len(selection), "atoms")
+  backbone_hierarchy = hierarchy.select(selection)
+  # backbone_hierarchy.show(dboutput)
+  chains = backbone_hierarchy.chains()
+  all_residues = []
+  cchains = list(chains)
+  nchains = len(cchains)
+  # print(cchains)
+  if nchains == 0:
+    print("Model contains no RNA", file=errorFile)
+    return []
+
+  for chain in cchains:
+    # print("chain ", chain.id)
+    conf = get_matching_conformer(chain, altcode)
+
+    names = conf.get_residue_names_padded(pad=False)
+    if all((n == "HOH" for n in names)):
+      # print("ignoring water chain")
+      continue
+    # print(list(names))
+    ids = conf.get_residue_ids(pad=False)
+    # #print(list(ids))
+    backbone_atoms = conf.atoms()
+
+    try:
+      residues = build_dihedrals(backbone_atoms, chain.id, conf.altloc)
+      all_residues.extend(residues)
+    except Failure as e:
+      print("chain ", chain.id, e.message, file=errorFile)
+      continue
+
+  # if name != "":  # useful for seeing what suites were generated
+  #   writeSuitesFile(residues, name)
+  if len(all_residues) == 0:
+    errorFile.write("read no residues: perhaps wrong alternate code\n")
+  return all_residues
+
+
+def get_matching_conformer(chain, altcode):
+  confo = [c for c in chain.conformers() if c.altloc in (altcode, '')]
+  if len(confo)>0:
+    conf = confo[0]
+  else:
+    conf = chain.only_conformer()
+  return conf
+
+
+def build_dihedrals(mainchain, chain_id, alt):
+  residues = []
+  # print (list(mainchain.extract_name()))  # gives a list of atom names
+
+  backbone = [atom for atom in mainchain if atom.name in bone]
+  residue = None
+
+  k = 0
+  i = find_start(backbone)
+  while k < len(backbone) - 3:
+    pivot = backbone[k + 1]
+    labels = pivot.fetch_labels()
+    res_number = int(labels.resseq)
+    if residue is None or res_number != residue.sequence:
+      # we are seeing our first angle of a new residue
+      residue = new_residue(residues, pivot)
+    name, angle = easy_make_dihedral(backbone, i, k)
+    # ----------------- on failure: new strategy ----------------------
+    if angle == 9999:
+      # this residue is strangely ordered due to partial alt conformers;
+      # build dihedrals the hard way hereafter
+      k = harder_build_chain(backbone, i, k, residue, residues,
+          chain_id, alt)
+      return residues
+
+    else:
+      residue.angle[i] = angle
+      k = k + 1
+      i = (i + 1) % 6
+  return residues
+
+
+def find_start(backbone):
+  for i in range(len(bones)):
+    if backbone[0].name == bones[i]:
+      return i
+  assert False, "backbone atom not recognized"
+
+
+def easy_make_dihedral(backbone, i, k):
+  dh = dihedrals[i % 6]
+  if backbone[k].name == dh[1][0]:
+    group = backbone[k:k+4]
+    group_names = tuple((a.name for a in group))
+    if group_names == dh[1]:
+      clump = backbone[k:k + 4]
+      points = [atom.xyz for atom in clump]
+      name = dh[0]
+      angle = dihedral_angle(sites=points, deg=True)
+      if angle < 0:  # suitename standard: all angles are positive
+        angle += 360.0
+      return name, angle
+    return "failure: " + str(group_names), 9999
+  return "failure: " + backbone[k].name + "!= " + dh[1][0], 9999
+
+
+def harder_build_chain(backbone, i, k, residue, residues, chain_id, alt):
+  if k - i - 2 < 0:
+    syn_chain, i = jump_start()
+    k = 0
+  else:
+    syn_chain = backbone[k:k + 3]
+  ksyn = k  # point on backbone where syn_chain shadowing begins
+  # print("syn_chain")
+  # for atom in syn_chain:
+  #   print(atom.name, atom.fetch_labels().resseq, atom.serial) #just for debugging
+
+  ii = i  # indicates which dihedral angle we are working on
+  while k < len(backbone) - 3:
+    k2 = harder_build_residue(syn_chain, ii, k-ksyn, residue)
+    # k2 is an index into syn_chain
+    if k2 > k-ksyn:
+      k = k2+ksyn
+      ii = 0
+      residue = new_residue(residues, syn_chain[k2+1])
+    else:
+      # things are horribly broken, just keep on looking for sanity
+      ok = resync(residues, syn_chain, k2)
+      if not ok:  # we can't fix it, terminate here
+        break
+      k += 1
+  # delete a possible final dead one
+  if residues[-1].is_dead() or \
+      (len(residues) > 1 and residues[-1].sequence == residues[-2].sequence):
+    del residues[-1]
+  return k
+
+
+def harder_build_residue(syn_chain, i, kk, residue):
+  # kk is an index to syn_chain
+  assert len(syn_chain) >= kk + 3, "We have run out of syn_chain"
+
+  res = syn_chain[-1].fetch_labels().resseq
+  res_atoms, res_names = get_one_residue(res)
+  res2_atoms, res2_names = get_one_residue(str(int(res)+1))
+
+  for j in range(i, 6):
+    atoms = res_atoms if j < 4 else res2_atoms
+    names = res_names if j < 4 else res2_names
+    if bones[j+3] in names:  # add atom 3 ahead of dihedral start
+      # so that we have 4 atoms to work with
+      x = names.index(bones[j+3])
+      atom = atoms[x]
+      syn_chain.append(atom)
+      # print("appending", atom.name, atom.fetch_labels().resseq, atom.serial)
+      name, angle = easy_make_dihedral(syn_chain, j, kk)
+      # print(" ", syn_chain[kk+1].fetch_labels().resseq, name, angle, kk, j)
+      residue.angle[j] = angle  # if second failure, we get 9999
+      kk = kk + 1
+    # syn_name = [a.name for a in syn_chain]
+  assert len(syn_chain) >= kk + 3, "We have run out of syn_chain"
+  return kk
+
+
+def resync(residues, syn_chain, k2):
+  "previous residue is not working for us, find next useful one"
+  oldres = syn_chain[-1].fetch_labels().resseq
+  for r in range(1, 21):
+    res = str(int(oldres) + r)
+    selection=manager.selection("resseq " + res + " and (name P or name O5' or name C5' or name C4' or name C3' or name O3')")
+    res_hierarchy = hierarchy.select(selection)
+    atoms = res_hierarchy.atoms()
+    if len(atoms) > 0:  # yes! we have found a useful residue
+      syn_chain.append(atoms[0])
+      return True
+  return False
+
+
+def jump_start():
+  """Used when a hiccup occurs at the beginning of the chain.
+  Get three items on the start of syn_chain so that harder_build_residue
+  has an adequate starting point to work with.
+  """
+  syn_chain = []
+  res = list(hierarchy.residue_groups())[0].resseq
+  res_atoms, res_names = get_one_residue(res)
+  res2_atoms, res2_names = get_one_residue(str(int(res)+1))
+
+  i=1
+  while i<7:
+    if bones[i] in res_names:
+      # we have found our starting point
+      atom = res_atoms[res_names.index(bones[i])]
+      syn_chain.append(atom)
+      break
+    i += 1
+  j = i+1; n = 1
+  jmax = 10  # if we haven't got them by now, give up
+  while n<3 and j<jmax:
+    atoms = res_atoms if j < 6 else res2_atoms
+    names = res_names if j < 6 else res2_names
+    if bones[j] in names:  # add atom 3 ahead of residue start
+      # so that we have 4 atoms to work with
+      atom = atoms[names.index(bones[j])]
+      syn_chain.append(atom)
+      n += 1
+    j += 1
+  if j>=jmax:
+    raise Failure(" unable to find a backbone")
+  return syn_chain, i
+
+
+def new_residue(residues, pivot):
+  """Start a new Residue object, whose angles will be filled in later as we
+     continue down the chain."""
+  residue = Residue("", "", [9999, 9999, 9999, 9999, 9999, 9999])
+  labels = pivot.fetch_labels()
+  res_number = int(labels.resseq)
+  alt = labels.altloc
+  chainID = labels.chain_id
+  residue.sequence = res_number
+  grandpa = pivot.parent().parent()
+  residue_name = grandpa.unique_resnames()[0].rjust(3)
+  # CIF files sometimes give shorter names, but PDB files always give 3 chars
+  # rjust bridges the gap
+  id = (" ", "1", "{:>2}".format(chainID), labels.resseq, " ", "{:1}".format(alt), residue_name)
+  residue.pointIDs = id
+  base = findBase(residue_name)
+  residue.base = base
+  if residue.base is not None:
+    residues.append(residue)
+  return residue
+
+
+def get_one_residue(res):
+  selection=manager.selection("resseq " + res)
+  res_hierarchy = hierarchy.select(selection)
+  atoms = res_hierarchy.atoms()
+  names = [a.name for a in atoms]
+  # print("residue ", res)
+  # for a in atoms:
+  #   print("  ", a.name, a.parent().parent().resseq)
+  return atoms, names
+
+
+# ----------------------------  tool room  -----------------------------
+
+# useful for seeing what suites were generated
+def writeSuitesFile(r, name):
+  file = open(name + ".suites", "w")
+  for r in all_residues:
+    id = ":".join(r.pointIDs)
+    angles = "".join([":{:.3f}".format(a) for a in r.angle])
+    # print(id + angles, file=file)
+  file.close()
+
+
+# for debug printouts
+def nameList(atoms):
+  list = [a.name for a in atoms]
+  return list
+
+
+def residueString(r):
+  sizes=[1, 1, 1, 3, 1, 1, 3]
+  id = "".join(["{num:{width}}:".format(num=x,width=y) for x, y in zip(r.pointIDs, sizes)])
+  angles = "".join(["{:8.3f}:".format(a) for a in r.angle])
+  return id + angles[:-1]
+
+
+# The following code was used to build the dihedrals list above
+def make_groups(output):
+  i = 0
+  for name in names:
+    output.write('("{}" ('.format(name))
+    for j in range(4):
+      output.write('"{}".format(a[i+j]), ')
+    output.write("))\n")
+    i = i+1
+
+
+# The following is for debugging and not for production use:
+if __name__ == '__main__':
+  main(sys.argv[1])
diff --git a/mmtbx/suitename/dualparse.py b/mmtbx/suitename/dualparse.py
new file mode 100644
index 0000000000..7a4be4a0cc
--- /dev/null
+++ b/mmtbx/suitename/dualparse.py
@@ -0,0 +1,120 @@
+#        Copyright 2021  Richardson Lab at Duke University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+import argparse
+from suiteninit import buildParser
+
+from iotbx.cli_parser import CCTBXParser
+import sys
+
+
+# suitename first
+def parseArgs(programClass, logger):
+  for i, arg in enumerate(sys.argv):  # make case insensitive
+    if arg[0]=="-":
+      sys.argv[i] = arg.lower()
+  parser = argparse.ArgumentParser()
+  buildParser(parser)  # set up the arguments for suite name
+  args, others = parser.parse_known_args()
+  args2 = argsToPhilForm(args, others)
+  args2.extend(others)
+
+  # create the cctbx aspect of the parser
+  parser2 = CCTBXParser(
+    program_class = programClass,
+    # custom_process_arguments = custom_args_proc,
+    logger=logger
+    )
+  # do it
+  namespace = parser2.parse_args(args2)
+  extracted = parser2.working_phil.extract()
+  # for name, value in extracted.suitename.__dict__.items():
+    # print(name, ":", value)
+  return parser2
+
+
+class TableItem(object):
+    main_name = ""
+    alias_name = ""
+    main_default = "false"
+    alias_default = "false"
+
+    def __init__(self, name1, name2, default1="", default2=""):
+        self.main_name = name1
+        self.alias_name = name2
+        self.main_default = default1
+        self.alias_default = default2
+
+
+aliasList=(
+('residuein', 'residuesin', False, False) ,
+('suitein', 'suitesin', False, False) ,
+('altid', 'altidval', 'A', 'A') ,
+('pointidfields', 'ptid', 7, 0) ,
+('etatheta', 'thetaeta', False, False),
+(None, 'angles', None, 9),  # deprecated items
+(None, 'resAngles', None, 6)
+)
+
+def argsToPhilForm(namespace, others):
+  """
+  Converting the parseArgs data to phil input format:
+  The raw namespace as an entry for every argument that could possibly appear,
+  including some that are aliases for others, some that are deprecated that
+  we will not support in this context. We determine that an alias has been used
+  if its value is different from its default. If so, the alias overrides the
+  main item. In any case, only one of a main/alias pair will be passed on.
+
+  The output of this function is a list of command line arguments suitable
+  for phil.
+  """
+  aliasTable = [TableItem(*entry) for entry in aliasList]
+  aliases = [item.alias_name for item in aliasTable]
+  aliasDict = {aliases[i]: aliasTable[i] for i in range(len(aliasTable))}
+  ###### res = {test_keys[i]: test_values[i] for i in range(len(test_keys))}
+  dictionary = vars(namespace)
+  for  key,value in dictionary.items(): # vars turns namespace to dict
+    if key in aliasDict:  # if key is an alias or deprecated
+      item = aliasDict[key]
+      if item.main_name is None:
+        pass
+      if value != item.alias_default:
+        # If the caller has actually specified a value for the alias,
+        # lift the alias's value to the main entry
+        dictionary[item.main_name] = value
+      dictionary[key] = None
+      # -- all aliases and deprecated items are suppressed from the dictionary
+
+  argsOut = []
+  for  key,value in dictionary.items(): # vars turns namespace to dict
+    if value is not None and value != "":
+      argsOut.append(str(key)+"="+str(value))
+  return argsOut
+
+# CCTBX first
+def parseArgs1(Program, logger):
+  # create the cc tbx aspect of the parser
+  parser = CCTBXParser(
+    program_class = Program,
+    # custom_process_arguments = custom_args_proc,
+    logger=logger)
+
+  # add the old suitename aspect of the parser
+  buildParser(parser)
+  # do it
+  namespace = parser.parse_args(sys.argv[1:])
+  return parser.working_phil
+
diff --git a/mmtbx/suitename/mainchain.py b/mmtbx/suitename/mainchain.py
new file mode 100644
index 0000000000..f2163d58fb
--- /dev/null
+++ b/mmtbx/suitename/mainchain.py
@@ -0,0 +1,24 @@
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+
+# This file has been merged into myangle.py
+
+def build_dihedrals1(mainchain):
+    print (mainchain.extract_name())
+    backbone = []
+    for atom in mainchain:
+        if atom.name in bones:
+            backbone.append(atom)
+    #help(mainchain)
+    backbone2 = [atom for atom in mainchain if atom.name in bones]
+    sanity = (backbone == backbone2)
+    for atom in backbone2:
+      print(atom.i_seq, atom.name)
+    i = find_start(backbone)
+    for k in range(len(backbone) - 3):
+        name, angle = easy_make_dihedral(backbone, i, k)
+        labels = backbone[k + 1].fetch_labels()
+        residue = int(labels.resseq)
+        print(residue, name, angle,)
+        i = (i + 1) % 6
+
diff --git a/mmtbx/suitename/suitename.bat b/mmtbx/suitename/suitename.bat
new file mode 100644
index 0000000000..60418dd29a
--- /dev/null
+++ b/mmtbx/suitename/suitename.bat
@@ -0,0 +1,1 @@
+python suitename.py %*
diff --git a/mmtbx/suitename/suitename.py b/mmtbx/suitename/suitename.py
new file mode 100644
index 0000000000..6d11680b02
--- /dev/null
+++ b/mmtbx/suitename/suitename.py
@@ -0,0 +1,583 @@
+#        Copyright 2021  Richardson Lab at Duke University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import division
+from __future__ import nested_scopes, generators, absolute_import
+from __future__ import with_statement, print_function
+
+"""
+Suitename is a program to aid model-building and perform validation of RNA
+backbone conformation. Rather than focusing on the 6 dihedral angles of each
+phosphate-to-phosphate nucleotide, it parses the RNA backbone into the more
+conformationally diagnostic unit of 7 dihedrals from delta of one nucleotide to
+delta of the next. This unit is known as a "suite", because it runs between the
+ribose sugars (and also between adjacent bases).
+
+These seven dihedral angles define a 7-dimensional space, analogous to the 2-D
+Ramachandran phi, psi space but vastly larger. Suitename's analysis encapsulates
+the results of research that used the well-ordered parts of a high-accuracy
+reference structure dataset to classify feasible conformations of these angles
+into several dozen well-populated clusters of observed cases
+(Richardson 2008, RNA 14: 465). Each "conformer" cluster has been given a
+2-character number-letter name, such as 1a for A-form conformation. The outer
+limit of each cluster is a "super-ellipse" shape (Gridgeman 1970, Math Gaz 54:31)
+which fits this data better than either a radially symmetric or a box shape.
+Its extent in each direction uses the overall standard deviation of all clusters
+in that dihedral angle. Some clusters are not cleanly separated, but are
+"satellites" of a "dominant" larger cluster such as 1a. Those are distinguished
+by a plane closer to the satellite center, at a distance split proportional to
+their relative populations.
+
+Delta-1 and delta are cleanly bimodal by C2'-endo vs C3'-endo ribose pucker and
+gamma is cleanly trimodal, so each input suite is initially placed into one of 12
+delta-delta-gamma "bins" or else is triaged as an outlier if outside all bins or
+outside the allowable ranges in any dihedral. Then, most of the detailed work of
+cluster definition is done in the 4-dimensional space of epsilon-1, zeta-1, alpha,
+and beta, with a final refinement of the fit calculation in all 7 parameters. Each
+suite within an input RNA model file is thus assigned the 2-character name of the
+cluster it falls into and a 0-1.0 "suiteness" value for how centrally it fits
+within that cluster. Suites that lie outside any cluster are marked as outliers,
+tagged with '!!' (bang-bang) as their 2-character code.
+
+
+Suitename can take input in several formats:
+1. A PDB or CIF format file, from which it will extract the dihedral angles.
+2. A "dangle" format file of nucleotides with the 6 dihedral angles already
+   calculated, most often provided by the mp_geo utility.
+3. A kinemage format file of suites with their 7 dihedral angles already calculated.
+The latter two formats may be provided as a file, or as standard input.
+
+Suitename can output in several formats:
+1. A report showing the classification and suiteness of each suite, and
+   optionally a statistical summary.
+2. A compact string showing only the 2-character conformer names in lower case,
+   alternating with the single upper-case letter for the base type.
+3. A kinemage file showing all the suites of that model as points in
+   7-dimensional space, colored according to their classification.
+Each suite is tagged with the ID information of the SECOND of the two
+nucleotides that make it up, because that one contains 4 of the 7 dihedrals.
+The output is typically written to standard output.
+
+The various operations of Suitename are also provided as a library, suites.py,
+for use by other CCTBX programs.
+
+"""
+from suitenamedefs import Cluster, Issue, failMessages
+from suitenamedefs import Holder, globals
+# suitenamedefs must come first!
+import suiteninit
+from suiteninit import bins, MAX_CLUSTERS, normalWidths, satelliteWidths
+from suiteninput import readResidues, readKinemageFile, buildSuites
+
+import sys
+import numpy as np
+from math import cos, pi
+
+
+version = "suitename.1.1.081021"
+dbCounter = 0
+dbTarget = -99  # triggers extra output on this suite
+
+# A collection of variables used for output
+outNote = Holder()
+outNote.version = version
+outNote.comment = ""
+outNote.wannabes = 0
+outNote.outliers = 0
+
+
+# ***main()******************************************************************
+"The main program, if run without using CCTBX or ProgramTemplate (old style)."""
+def main(inStream=None, outFile=None, errorFile=None, optionsIn=None):
+  # inStream is used in internal testing
+  global options
+  if optionsIn is None:
+    # suitename is being run in standalone mode
+    options = suiteninit.parseCommandLine()
+  else:
+    options = optionsIn
+  globals.options = options
+  # makes them available to other parts of this program
+
+  global dbCounter  # for debugging KPB 210222
+  if options.version:
+    print(version)
+    return
+
+  # 1. read the input
+  if inStream:
+    inFile = inStream
+  elif options.infile != "":
+    inFile = open(options.infile)
+  else:
+    inFile = sys.stdin
+
+  if not outFile:  outFile = sys.stdout
+  if not errorFile: errorFile = sys.stderr
+
+  suites = read(inFile, errorFile)
+  suites = compute(suites)
+  finalStats()
+  write(outFile, suites)
+
+
+def loadOptions(optionsIn):
+  # A service routine, just to be called by suites.py
+  global options
+  options = optionsIn
+
+
+def read(inFile, errorFile=sys.stderr):
+  "Read input as residues or as suites, resulting in suites."
+  if options.suitein:
+    suites = readKinemageFile(inFile)
+    if len(suites) == 0:
+      errorFile.write("read no suites: perhaps wrong type of kinemage file\n")
+  else:
+    residues = readResidues(inFile)
+    if len(residues) == 0:
+      errorFile.write("read no residues: perhaps wrong alternate code\n")
+      suites = []
+    else:
+      suites = buildSuites(residues)
+      suites = suites[:-1]
+  return suites
+
+
+def compute(suites):
+  """ Process the suites, resulting in various attributes added to each
+     suite object, and statistics accrued in cluster objects."""
+  global dbCounter
+  for s in suites:
+    if not s.validate():
+      if options.test:
+        sys.stderr.write("! failed validation: {}\n".format(s.pointID))
+      annotate(s, bins[13], bins[13].cluster[0], 0, 0, " tangled ",
+            "", "", "", "")
+      continue  # makes sense but does not ?? match C version
+
+    # At this point we have a complete suite
+    bin, issue, text, pointMaster = evaluateSuite(s)
+    pointColor = "white"
+    if bin is None:
+      s.cluster = bins[0].cluster[0]
+      bins[0].cluster[0].count += 1
+      annotate(s, bins[0], bins[0].cluster[0], 0, 0, text, issue,
+          "", pointMaster, pointColor)
+    else:
+      memberPack = membership(bin, s)
+      (cluster, distance, suiteness, situation, comment,
+          pointMaster, pointColor) = memberPack
+      annotate(s, bin, cluster, distance, suiteness, situation, issue,
+          comment, pointMaster, pointColor)
+    dbCounter += 1
+  return suites
+
+
+def annotate(suite, bin, cluster, distance, suiteness, situation,
+                issue, comment, pointMaster, pointColor):
+  """ Add attributes to suite objects, thus recording the results of
+  Suitename computations."""
+  suite.bin = bin
+  suite.cluster = cluster
+  suite.distance = distance
+  suite.suiteness = suiteness
+  suite.situation = situation
+  suite.issue = issue
+  suite.comment = comment
+  suite.pointMaster = pointMaster
+  suite.pointColor = pointColor
+
+
+def write(outFile, suites):
+  "Write the output in the requested format."
+  from suitenout import output
+  # late import is required to ensure that globals.options is in place
+  # for the suitenout module
+  output(outFile, suites, outNote)
+
+
+# *** evaluateSuite and its tools ***************************************
+
+# Boundaries of various angle ranges
+epsilonmin = 155
+epsilonmax = 310  # 070130
+delta3min = 60
+delta3max = 105  #  changed by S.J. on 06/06/2011
+delta2min = 125
+delta2max = 165
+gammapmin = 20
+gammapmax = 95  # max 070326
+gammatmin = 140
+gammatmax = 215  # max 070326
+gammammin = 260
+gammammax = 335  # max 070326
+alphamin = 25
+alphamax = 335
+betamin = 50
+betamax = 290
+zetamin = 25
+zetamax = 335
+
+# triage table for yes-no angles:
+# each of these filters, applied to a suite, will provide a
+# true or false answer as to whether this angle is in a reasonable range.
+# pointMaster is for grouping points in kinemage display
+
+# data per line: angle index, min, max, code, text, pointMaster
+triageFilters = {
+  "epsilon": (2, epsilonmin, epsilonmax, Issue.EPSILON_M, "e out", "E"),
+  "alpha":   (4, alphamin, alphamax, Issue.ALPHA, " a out", "T"),
+  "beta":    (5, betamin, betamax, Issue.BETA, " b out", "T"),
+  "zeta":    (3, zetamin, zetamax, Issue.ZETA_M, " z out", "T"),
+}
+
+def triage(selector, suite):
+  # if angle lies outside the acceptable range, triage immediately
+  filter = triageFilters[selector]
+  index, min, max, failCode, failText, pointMaster = filter
+  if suite.angle[index] < min or suite.angle[index] > max:
+    return False, failCode, failText, pointMaster
+  else:
+    return True, None, None, ""
+
+
+# The more complex angles are handled by a "sieve".
+# A sieve will determine whether an angle is within one of several ranges
+# and provide an appropriate code indicating the range.
+# This is handled by the sift() function.
+sieveDelta = (
+  (delta3min, delta3max, 3),
+  (delta2min, delta2max, 2),
+)
+
+sieveGamma = (
+  (gammatmin, gammatmax, "t"),
+  (gammapmin, gammapmax, "p"),
+  (gammammin, gammammax, "m"),
+)
+
+def sift(sieve, angle, failCode):
+  for filter in sieve:
+    min, max, code = filter
+    if min <= angle <= max:
+      return code, "", ""
+  failMessage = failMessages[failCode]
+  return None, failCode, failMessage
+
+
+def evaluateSuite(suite):
+  '''Determine whether suite falls into one of 12 predefined bins,
+   and if so, which.'''
+  global bins
+
+  # The order of triage operations, though it may seem arbitrary,
+  # was carefully chosen by the scientists.
+  ok, failCode, notes, pointMaster = triage("epsilon", suite)
+  if not ok:
+    return None, failCode, notes, pointMaster
+
+  # Angles with several meaningful ranges:
+  # for each angle, find out which range it lies in, or none
+  # this becomes a selector to help choose a bin
+  puckerdm, failCode, notes = sift(sieveDelta, suite.deltaMinus, Issue.DELTA_M)
+  if not puckerdm:
+    return None, failCode, notes, "D"
+
+  puckerd, failCode, notes = sift(sieveDelta, suite.delta, Issue.DELTA)
+  if not puckerd:
+    return None, failCode, notes, "D"
+
+  gammaname, failCode, notes = sift(sieveGamma, suite.gamma, Issue.GAMMA)
+  if not gammaname:
+    return None, failCode, notes, "T"
+
+  ok, failCode, notes, pointMaster = triage("alpha", suite)
+  if not ok:
+    return None, failCode, notes, pointMaster
+
+  ok, failCode, notes, pointMaster = triage("beta", suite)
+  if not ok:
+    return None, failCode, notes, pointMaster
+
+  ok, failCode, notes, pointMaster = triage("zeta", suite)
+  if not ok:
+    return None, failCode, notes, pointMaster
+
+  # We have passed the test: now use this information to select a bin
+  bin = bins[(puckerdm, puckerd, gammaname)]
+  # Bins is an associated dictionary indexed by the triplet of three angle classifiers
+  # Each unique triplet of classifiers selects one unique bin, for a total of 12 bins.
+  return bin, None, None, ""
+
+
+# ***membership()***************************************************************
+
+def membership(bin, suite):
+  """
+  Having selected a bin, we look for the best cluster match within that bin
+
+  Cluster membership:
+  Three of the seven dimensions (delta-1, gamma, and delta) have already been
+  used to select a bin. These angles are fairly cut and dried, in-or-out type
+  decisions. The other four are used in this function to select the closest
+  cluster in four-dimensional space. Much research has gone into determining
+  meaningful boundaries for clusters in this four dimensional space. After the
+  nearest cluster is selected, we go back and use seven-dimensional space to
+  refine the distance measurement; this may change whether our suite is in or
+  out of the cluster, but it will not change the selection of cluster.
+
+  New research would be required to determine appropriate cluster boundaries in
+  seven dimensional space.
+  """
+  matches = np.full(MAX_CLUSTERS, 999.9)
+  matchCount = 0
+  comment = ""
+  pointMaster = ""
+  pointColor = "white"
+  lDominant = False
+
+  if bin.dominant > 0:  # this bin has a dominant cluster, note it
+    dominantJ = bin.dominant
+    domCluster = bin.cluster[bin.dominant]
+
+  # find the closest cluster
+  # search every cluster in the bin except cluster 0, which is for outliers
+  closestD = 999
+  closestCluster = bin.cluster[0]  # default, representing an outlier
+  for j, c in enumerate(bin.cluster[1:], 1):
+    if c.status == "wannabe" and options.nowannabe:
+      continue
+    distance = hyperEllipsoidDistance(
+        suite.angle, bin.cluster[j].angle, 4, normalWidths
+    )
+    if distance < closestD:
+      closestD = distance
+      closestJ = j
+      closestCluster = c
+    matches[j] = distance
+    if distance < 1:  # suite could be a member of this cluster
+      matchCount += 1
+      if c.dominance == "dom":
+        lDominant = True
+        # there is a close dominant cluster to consider
+
+  if matchCount == 1:
+    theCluster = closestCluster
+    situation = "1-only-one"
+
+  elif matchCount > 1 and not lDominant:
+    # dominant cluster is not a possible cluster
+    # just output than minimum distance match
+    theCluster = closestCluster
+    situation = "{}-None-dom".format(matchCount)
+
+  elif matchCount > 1:  # and lDominant
+    # find the closest cluster that is not the dominant cluster
+    closestNonD = 999
+    for j, c in enumerate(bin.cluster[1:], 1):
+      if c.status == "wannabe" and options.nowannabe:
+        continue
+      if matches[j] < closestNonD and c.dominance != "dom":
+        closestNonD = matches[j]
+        closestJ = j
+        theCluster = c
+
+    if theCluster.dominance == "sat":
+      # We need to distinguish carefully whether our suite
+      # is in the dominant or satellite cluster
+      theCluster, closestJ, situation = domSatDistinction(
+          suite, domCluster, theCluster, matches, matchCount
+      )
+    else:
+      if matches[dominantJ] < matches[closestJ]:
+        closestJ = dominantJ
+        theCluster = domCluster
+      situation = "{}-not-sat".format(matchCount)
+  else:
+    # no match, it's an outlier
+    closestJ = 0
+    theCluster = closestCluster
+    if closestCluster.name != "!!":
+      situation = "outlier distance {:.3}".format(closestD)
+    else:
+      situation = "vacant bin"
+    outNote.outliers += 1
+    pointMaster = "O"
+    pointColor = "white"
+
+  # final computation of suiteness
+  # this time we use all 7 dimensions
+  if dbCounter >= dbTarget and dbCounter <= dbTarget + 1:  # KPB debug tool
+    print(suite.pointID)
+    print(suite.angle)
+    print(theCluster.name)
+
+  distance = hyperEllipsoidDistance(suite.angle, theCluster.angle, 7, normalWidths)
+  # this calculation can assign or deassign a cluster
+  if distance <= 1:
+    suiteness = (cos(pi * distance) + 1) / 2
+    if suiteness < 0.01:
+      suiteness = 0.01
+  else:
+    if closestJ != 0:
+      # 7D distance forces this suite to be an outlier
+      # so we deassign it here
+      closestJ = 0
+      comment = "7D dist {}".format(theCluster.name)
+    theCluster = bin.cluster[0]  # outlier
+    suiteness = 0
+
+  theCluster.count += 1
+  suite.cluster = theCluster
+  if theCluster.status == "wannabe" and not options.nowannabe:
+    outNote.wannabes = 1  # once set, stays set
+  pointColor = 0  # will be handled later!!
+  if options.test:
+    print(" [suite: %s %s 4Ddist== %f, 7Ddist== %f, suiteness==%f] \n" % \
+        (theCluster.name, suite.pointID[:11], closestD, distance, suiteness))
+  return theCluster, distance, suiteness, situation, comment, pointMaster, pointColor
+
+
+def domSatDistinction(suite, domCluster, satCluster, matches, matchCount):
+  """
+  The special case where the two best matches are the dominant cluster in the
+  bin, and its satellite. Since dominant clusters are typically much larger
+  than their satellites, we don't just accept whatever cluster is closest.
+  First, we do a vector test to see if our suite actually lies in between the
+  two. If so, we use specially chosen "width" parameters in each dimension to
+  do a comparison. In a few cases, even more specially chosen width parameters
+  (satelliteInfo) apply to this exact dominant/satellite pair. Thus we test
+  our suite's position against two hyperellipsoids of unequal size and shape,
+  to see which is the better fit.
+  """
+  closestCluster = satCluster
+  closestJ = satCluster.ordinal
+  dominantJ = domCluster.ordinal
+
+  # use vector properties of numpy.array to determine difference vectors
+  domToPoint = domCluster.angle - suite.angle
+  satToPoint = satCluster.angle - suite.angle
+  domToSat = domCluster.angle - satCluster.angle
+  satToDom = -domToSat
+
+  dps = narrowDotProduct(domToPoint, domToSat, 4)
+  spd = narrowDotProduct(satToPoint, satToDom, 4)
+
+  if dps > 0 and spd > 0:
+    # the trickiest case: point is between dom and sat
+    domWidths = normalWidths.copy()
+    if options.satellites:
+      satWidths = satelliteWidths.copy()
+    else:
+      satWidths = normalWidths.copy()
+    if satCluster.satelliteInfo is not None:
+      modifyWidths(domWidths, satWidths, satCluster.satelliteInfo)
+    disttodom = hyperEllipsoidDistance(suite.angle, domCluster.angle, 4, domWidths)
+    disttosat = hyperEllipsoidDistance(suite.angle, satCluster.angle, 4, satWidths)
+    if disttodom < disttosat:
+      closestJ = dominantJ
+      closestCluster = domCluster
+    situation = "{}-BETWEEN-dom-sat({:7.3}|{:7.3})".format(matchCount,disttodom,disttosat)
+    # else the satellite cluster remains the chosen cluster
+
+  else:
+    # the point is not in between
+    # just assign by closest standard distance evaluation
+    if matches[dominantJ] < matches[closestJ]:
+      closestJ = dominantJ
+      closestCluster = domCluster
+    if dps <= 0:
+      situation = "{}-OUTSIDE-dom".format(matchCount)
+    else:
+      situation = "{}-OUTSIDE-sat".format(matchCount)
+
+  return closestCluster, closestJ, situation
+
+
+# *** Gathering some statistics *********************************************
+
+def finalStats():
+  "Compute summary statistics for each bin"
+  for bin in bins.values():
+    for c in bin.cluster:
+      bin.count += c.count
+
+
+def clearStats():
+  "Remove all accrued statistics, appropriate for a fresh new run"
+  from suitenout import clearStatistics
+  clearStatistics()
+
+
+# *** The fancy math ********************************************************
+
+# This variable was experimental but we have settled on 3:
+power = 3
+
+def hyperEllipsoidDistance(suiteAngles, clusterAngles, nAngles, widthArray):
+  global dbCounter
+  if nAngles == 4:
+    workRange = range(2, 6)
+  else:
+    workRange = range(1, 8)
+
+  summation = 0
+  for k in workRange:
+    delta = abs(suiteAngles[k] - clusterAngles[k])
+    delta = delta / widthArray[k]
+    delToPower = pow(delta, power)
+    summation = summation + delToPower
+    # if dbCounter >= dbTarget and nAngles > 4:  # KPB debug 120221
+    #     sys.stderr.write("db=%3d, k=%d, del=%8.4f, delpower=%10.6f, dpower=%10.6f\n" %
+    #                 (dbCounter, k, delta, delToPower, summation) )
+  result = pow(summation, 1 / power)
+  # if dbCounter == dbTarget and nAngles > 4:
+  #     sys.stderr.write("final = %7.3f\n" % result)
+  return result
+
+
+def narrowDotProduct(a, b, nAngles):
+  """The narrow dot product involves only a subset of the dimensions,
+  either 4 or 7. In practice, only 4 is in use."""
+  if nAngles == 4:
+    return np.dot(a[2:6], b[2:6])
+  else:
+    return np.dot(a[1:8], b[1:8])
+
+
+def modifyWidths(dom, sat, satInfo):
+  for m in range(9):
+    if satInfo.satelliteWidths[m] > 0:
+      sat[m] = satInfo.satelliteWidths[m]
+    if satInfo.dominantWidths[m] > 0:
+      dom[m] = satInfo.dominantWidths[m]
+
+
+if (__name__ == "__main__"):
+  main()
+
+
+# CHANGE LOG (yymmdd datestamp format):
+# 0.2.070524 preserve chi-1 and chi, so could preserve eta, theta
+# 0.3.070525 general read dangle record for, e.g.,  eta, theta
+# 0.3.070628 triage reports zeta-1, epsilon-1, delta-1,... Ltriage codes
+# 0.3.070803 notes: rearranged suitenhead.h/janesviews ...
+#                  put something in to say what veiws mean
+# 0.3.070919 3g wannabe (tRNA TpseudoUC loop)
+# 0.3.110606 range of delta updated by S.J.
+# 01/07/2014 S.J. updated so that it can take input with alternate conformations,
+#            and by default will calculate the suite for altA
+# 09/18/2014 S.J. updated so that suitename will ignore DNA residues
+# 03/01/2021 Ken Brooks converted to Python, major reorganization
diff --git a/mmtbx/suitename/suitename.sh b/mmtbx/suitename/suitename.sh
new file mode 100644
index 0000000000..de09f24f61
--- /dev/null
+++ b/mmtbx/suitename/suitename.sh
@@ -0,0 +1,3 @@
+#!/bin/sh
+
+python suitename $*
diff --git a/mmtbx/suitename/suitenamedefs.py b/mmtbx/suitename/suitenamedefs.py
new file mode 100644
index 0000000000..147237fa79
--- /dev/null
+++ b/mmtbx/suitename/suitenamedefs.py
@@ -0,0 +1,378 @@
+"Class definitions used throughout suitename"
+
+#        Copyright 2021  Richardson Lab at Duke University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+import numpy as np
+from numpy import array
+from enum import Enum
+
+
+class Holder(object):
+  "these objects exist only to have case specific attributes added to them"
+  pass
+
+globals = Holder()
+# globals.options will be created in the main program
+
+
+# reasons why a suite may fail to be classified:
+Issue = Enum('Issue', 'DELTA_M EPSILON_M ZETA_M ALPHA BETA GAMMA DELTA')
+reasons = {
+  Issue.DELTA_M:    "delta-1",
+  Issue.EPSILON_M:  "epsilon-1",
+  Issue.ZETA_M:     "zeta-1",
+  Issue.ALPHA:      "alpha",
+  Issue.BETA:       "beta",
+  Issue.GAMMA:      "gamma",
+  Issue.DELTA:      "delta"
+}
+
+failMessages = {
+  Issue.DELTA_M:    "bad deltam",
+  Issue.GAMMA:      "g out",
+  Issue.DELTA:      "bad delta"
+}
+
+
+class Bin(object):
+  """ primary (coarse grained) classification of suites """
+  # permanence properties
+  name = ""
+  ordinal = 0
+  cluster = ()
+  # a tuple of cluster objects
+  dominant = -1
+  # statistics gathered during the run
+  count = 0
+  active = False
+
+  def __init__(self, ordinal, name, clusters=()):
+    self.ordinal = ordinal
+    self.name = name
+    self.cluster = clusters
+    self.dominant = -1
+    self.active = False
+
+    for i, c in enumerate(clusters):
+      if c.dominance == "dom":
+        self.dominant = i
+        break
+
+
+class Cluster(object):
+  """secondary (fine grained) classification of suites"""
+  # intrinsic data:
+  ordinal = 0       # its place in the bin
+  name = ""         # the name of the cluster
+  status = ""       # certain, wannabe, triaged, outlier, nothing, incomplete
+  clusterColor = "" # kinemage color names
+  dominance = ""    # dom, sat, ord, out, tri, inc
+  satelliteInfo = None    # present only if this cluster is a satellite
+  angle = ()        # tuple of 9 angles:
+  # chiMinus, deltaMinus, epsilon, zeta, alpha, beta, gamma, delta, chi
+
+  # gathered statistics:
+  count = 0  # number of data points found in this cluster
+  suitenessSum = 0
+  suitenessCounts = None
+
+  def __init__(self, ordinal, name, status, color, dominance, angles):
+    self.ordinal = ordinal
+    self.name = name
+    self.LOK = (name != "!!")
+    self.status = status
+    self.clusterColor = color
+    self.dominance = dominance
+    self.angle = array(angles)
+    self.satelliteInfo = None
+    # sometimes modified in suiteninit.buildBin
+    self.suitenessCounts = np.zeros(12)
+    self.suitenessSum = 0
+
+
+class SatelliteInfo(object):
+  # numbers used when suite is between satellite and dominant centers
+  name = ""
+  satelliteWidths = ()   # vector of 9 angles
+  dominantWidths = ()    # vector of 9 angles
+
+  def __init__(self, name, satelliteWidths, dominantWidths):
+    self.name = name
+    self.satelliteWidths = satelliteWidths
+    self.dominantWidths = dominantWidths
+
+
+class Residue(object):
+  '''
+  # A residue as normally read in, consisting of its six dihedral angles
+  Used only briefly as input.
+  '''
+  sequence = -1  # sequence number in PDB file
+  pointIDs = []
+  base = " "    # A, C, G, U, ...
+  angle = np.empty(0)  # will have 6 or 7 elements:
+  # alpha, beta, gamma, delta, epsilon, zeta [, chi]
+
+  def __init__(self, ID, base, angles):
+    self.pointIDs = ID
+    self.base = base
+    self.angle = angles
+
+  def is_dead(self):
+    dead = [a == 9999 for a in self.angle]
+    return all(dead)
+
+
+  # def __str__(self):
+  #     string = "%s:%s:%s:%4s:%s:%s:%s:%s:%s:%s:%s:%s:%s" \
+  #        % (" ",
+  #           "1",
+  #           chainID,
+  #           resnum,
+  #           i_code,
+  #           altloc,
+  #           resname,
+  #           alpha,
+  #           beta,
+  #           gamma,
+  #           delta,
+  #           epsilon,
+  #           zeta)
+  #     return string
+
+
+  # nicknames: for ease of reading the code, each angle is given
+  # a meaningful alias. Here they are:
+  # 0   alpha
+  # 1   beta
+  # 2   gamma
+  # 3   delta
+  # 4   epsilon
+  # 5   zeta
+  # 7   chi
+  @property
+  def alpha(self):
+    return self.angle[0]
+
+  @alpha.setter
+  def alpha(self, value):
+    self.angle[0] = value
+
+  @property
+  def beta(self):
+    return self.angle[1]
+
+  @beta.setter
+  def beta(self, value):
+    self.angle[1] = value
+
+  @property
+  def gamma(self):
+    return self.angle[2]
+
+  @gamma.setter
+  def gamma(self, value):
+    self.angle[2] = value
+
+  @property
+  def delta(self):
+    return self.angle[3]
+
+  @delta.setter
+  def delta(self, value):
+    self.angle[3] = value
+
+  @property
+  def epsilon(self):
+    return self.angle[4]
+
+  @epsilon.setter
+  def epsilon(self, value):
+    self.angle[4] = value
+
+  @property
+  def zeta(self):
+    return self.angle[5]
+
+  @zeta.setter
+  def zeta(self, value):
+    self.angle[5] = value
+
+  @property
+  def chi(self):
+    if len(self.angle) > 6:
+      return self.angle[6]
+
+  @chi.setter
+  def chi(self, value):
+    if len(self.angle) > 6:
+      self.angle[6] = value
+
+
+class Suite(object):
+  '''
+  The set of angles forming the linkage BETWEEN residues.
+  This is the core data structure used in most operations of the program.
+  '''
+  pointID = ()
+  base = " "    # A, C, G, U, ...
+  angle = np.empty(0) # will become an np.array of 9 angles:
+  # chiMinus, deltaMinus, epsilon, zeta, alpha, beta, gamma, delta, chi
+
+  # fields computed during analysis:
+  valid = False  # False means an incomplete, malformed suite
+  cluster = None  # The cluster to which it is assigned
+  suiteness = 0.0
+  distance = 0.0
+  situation = ""  # by what logical path this cluster was assigned
+  pointMaster = ""
+  pointColor = ""
+
+  def __init__(self, ID, base, angles=None):
+    self.pointID = ID
+    self.base = base
+    if angles is None:
+      self.angle = np.full(9, 0.0)
+    else:
+      self.angle = angles
+    self.suiteness = 0.0
+    self.distance = 0.0
+    self.notes = ""
+    self.dbflag = False
+
+  def validate(self):
+    # make sure that angles deltaMinus through delta are reasonable
+    self.valid = True
+    for i in range(1, 8):
+      if self.angle[i] < 0 or self.angle[i] > 360:
+        self.valid = False
+    return self.valid
+
+  # nicknames: for ease of reading the code, each angle is given
+  # a meaningful alias. Here they are:
+  # 0   chiMinus
+  # 1   deltaMinus
+  # 2   epsilon
+  # 3   zeta
+  # 4   alpha
+  # 5   beta
+  # 6   gamma
+  # 7   delta
+  # 8   chi
+
+  @property
+  def chiMinus(self):
+    return self.angle[0]
+
+  @chiMinus.setter
+  def chiMinus(self, value):
+    self.angle[0] = value
+
+  @property
+  def deltaMinus(self):
+    return self.angle[1]
+
+  @deltaMinus.setter
+  def deltaMinus(self, value):
+    self.angle[1] = value
+
+  @property
+  def epsilon(self):
+    return self.angle[2]
+
+  @epsilon.setter
+  def epsilon(self, value):
+    self.angle[2] = value
+
+  @property
+  def zeta(self):
+    return self.angle[3]
+
+  @zeta.setter
+  def zeta(self, value):
+    self.angle[3] = value
+
+  @property
+  def alpha(self):
+    return self.angle[4]
+
+  @alpha.setter
+  def alpha(self, value):
+    self.angle[4] = value
+
+  @property
+  def beta(self):
+    return self.angle[5]
+
+  @beta.setter
+  def beta(self, value):
+    self.angle[5] = value
+
+  @property
+  def gamma(self):
+    return self.angle[6]
+
+  @gamma.setter
+  def gamma(self, value):
+    self.angle[6] = value
+
+  @property
+  def delta(self):
+    return self.angle[7]
+
+  @delta.setter
+  def delta(self, value):
+    self.angle[7] = value
+
+  @property
+  def chi(self):
+    return self.angle[8]
+
+  @chi.setter
+  def chi(self, value):
+    self.angle[8] = value
+
+
+# The great variety of codes that may represent each base in the input file
+NAListA = ":ADE:  A:A  : Ar:ATP:ADP:AMP:T6A:1MA:RIA:  I:I  :"
+NAListG = ":GUA:  G:G  : Gr:GTP:GDP:GMP:GSP:1MG:2MG:M2G:OMG: 7MG:"
+NAListC = ":CYT:  C:C  : Cr:CTP:CDP:CMP:5MC:OMC:"
+NAListU = ":URA:URI:  U: Ur:U  :UTP:UDP:UMP:5MU:H2U:PSU:4SU:"
+NAListY = ": YG:YG :  Y:Y  :"
+#NAListT = ":THY:  T:T  : Tr:TTP:TDP:TMP:"
+IgnoreDNAList = ": DA: DG: DC: DT:THY:  T:T  : Tr:TTP:TDP:TMP:"
+
+# out of the noise, determine the base
+def findBase(baseCode):
+  if IgnoreDNAList.find(baseCode) >= 0:
+    return None  # we ignore DNA residues
+  elif len(baseCode) != 3:
+    return "Z"
+
+  if NAListA.find(baseCode) >= 0:
+    base = "A"
+  elif NAListG.find(baseCode) >= 0:
+    base = "G"
+  elif NAListC.find(baseCode) >= 0:
+    base = "C"
+  elif NAListU.find(baseCode) >= 0:
+    base = "U"
+  elif NAListY.find(baseCode) >= 0:
+    base = "Y"
+  else:
+    base = "?"
+  return base
diff --git a/mmtbx/suitename/suiteninit.py b/mmtbx/suitename/suiteninit.py
new file mode 100644
index 0000000000..fb59bdfeb5
--- /dev/null
+++ b/mmtbx/suitename/suiteninit.py
@@ -0,0 +1,426 @@
+"""
+This is a self initializing module that embodies the data around which
+this program is built. It exports its primary data structure:
+  bins    the bin and cluster definitions
+"""
+
+#        Copyright 2021  Richardson Lab at Duke University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+import suitenamedefs
+from suitenamedefs import Bin, Cluster, SatelliteInfo
+
+from numpy import array
+import argparse, sys
+
+MAX_CLUSTERS = 16  # practical, observed limit of clusters in a bin
+
+
+# ***parseCommandLine()*******************************************************
+def parseCommandLine():
+    "Parse a command line, returning a parseargs.Namespace"
+    for i, arg in enumerate(sys.argv):
+        sys.argv[i] = arg.lower()
+
+    parser = argparse.ArgumentParser()
+    buildParser(parser)
+
+    # now actually parse them
+    args = parser.parse_args()
+    if args.ptid:
+        args.pointidfields = args.ptid
+    pass
+    return args
+
+
+def buildParser(parser):
+    "Set up the expected arguments, ready for parsing"
+
+    # the input file may be given as an argument or as a redirect
+    parser.add_argument("infile", nargs="?", default="")
+
+    # input styles
+    parser.add_argument("--residuein", "-residuein", action="store_true")
+    parser.add_argument("--residuesin", "-residuesin", action="store_true")
+    parser.add_argument("--suitein", "-suitein", action="store_true")
+    parser.add_argument("--suitesin", "-suitesin", action="store_true")
+
+    # output styles (default is --report)
+    outputStyle = parser.add_mutually_exclusive_group()
+    outputStyle.add_argument("--report", "-report", action="store_true")
+    outputStyle.add_argument("--string", "-string", action="store_true")
+    outputStyle.add_argument("--kinemage", "-kinemage", action="store_true")
+    # output modifiers
+    parser.add_argument("--chart", "-chart", action="store_true")
+      # a modifier to --report, suppress as the statistical summary
+    parser.add_argument("--causes", "-causes", action="store_true")
+      # a modifier to --report, reveals algorithm details
+    parser.add_argument("--nosequence", "-nosequence", action="store_true")
+      # a modifier to --string, places ':' instead of residue code
+
+    # additional options
+    parser.add_argument("--satellites", "-satellites", action="store_true")
+    parser.add_argument("--nowannabe", "-nowannabe", action="store_true")
+    parser.add_argument("--noinc", "-noinc", action="store_true")
+    parser.add_argument("--thetaeta", "-thetaeta", action="store_true")
+    parser.add_argument("--etatheta", "-etatheta", action="store_true")
+    parser.add_argument("--test", "-test", action="store_true")
+    parser.add_argument("--version", "-version", action="store_true")
+
+    # numerical options
+    parser.add_argument("--anglefields", "-anglefields", type=int, default=9)
+    parser.add_argument("--pointidfields", "-pointidfields", type=int, default=7)
+    parser.add_argument("--ptid", "-ptid", type=int, default=0)
+    parser.add_argument("--altid", "-altid", type=str, default="A")
+    parser.add_argument("--altidval", "-altidval", type=str, default="A")
+    parser.add_argument("--altidfield", "-altidfield", type=int, default=6)
+
+    # the following are deprecated:
+    parser.add_argument("--angles", type=int, default=9)
+    parser.add_argument("--resAngles", type=int, default=6)
+    parser.add_argument("--oneline", "-oneline", action="store_true")
+    #   "--help" is automatically available, it summarizes this list.
+    return parser
+
+
+# *** codes to match various residues ****************************
+
+idFields = 5
+match_list = (
+    (":ADE:  A:A  : Ar:ATP:ADP:AMP:T6A:1MA:RIA:  I:I  :", "A"),
+    (":GUA:  G:G  : Gr:GTP:GDP:GMP:GSP:1MG:2MG:M2G:OMG: YG: 7MG:YG :", "G"),
+    (":CYT:  C:C  : Cr:CTP:CDP:CMP:5MC:OMC:", "C"),
+    (":URA:URI:  U: Ur:U  :UTP:UDP:UMP:5MU:H2U:PSU:4SU:", "U"),
+    (":THY:  T:T  : Tr:TTP:TDP:TMP:", "T"),
+)
+
+
+# *** Special data relating to satellite clusters ****************
+
+# This function operates on the satelliteData list below
+# it creates an associated dictionary based on the name
+def buildSatelliteTable():
+    global satelliteTable
+    satelliteTable = {}
+    for item in satelliteData:
+        name = item[0]
+        satWidths = item[1]
+        domWidths = item[2]
+        satelliteTable[name] = SatelliteInfo(name, satWidths, domWidths)
+
+
+# The satellite data:
+# The widths below are used to determine multidimensional hyperellipsoidal distances.
+# A distance <= 1  is considered "in" the cluster
+
+# There are three tiers of widths here:
+# 1. The normal widths, deltamw etc.
+# 2. The general satellite widths, epsilonsatw etc.
+# 3. The special satellite widths in the table farther below, satelliteData
+# The normalWidths are used for a typical cluster.
+# The satelliteWidths are an exception, used for satellite clusters.
+# The satelliteData are exceptions to the satelliteWidths, for certain specific
+# satellite clusters.
+
+# SITUATION as of 210213:
+#   The satelliteWidths are used ONLY if the --satellites arg is used
+
+clusterhalfwidthsversion = "070328"
+deltamw = 28
+epsilonw = 60
+epsilonsatw = 50  # satw 070328
+zetaw = 55
+zetasatw = 50  # satw 070328
+alphaw = 50
+alphasatw = 45  # satw 070328
+betaw = 70
+betasatw = 60  # satw 070328
+gammaw = 35
+deltaw = 28
+
+
+# width arrays set the widths of clusters in the various dimensions
+# the zeroes on either end may someday be replaced with widths for the
+# chi angles
+normalWidths = array((0, deltamw, epsilonw, zetaw, alphaw, betaw, gammaw, deltaw, 0))
+satelliteWidths = array(
+    (0, deltamw, epsilonsatw, zetasatw, alphasatw, betasatw, gammaw, deltaw, 0)
+)
+
+satelliteData = (
+    #  sat         9 angles sat widths                    9 angles dom width
+    ("1m", (0, 0, 0, 0, 0, 32, 0, 0, 0), (0, 0, 0, 0, 0, 64, 0, 0, 0)),
+    ("1L", (0, 0, 18, 0, 0, 18, 0, 0, 0), (0, 0, 70, 0, 0, 70, 0, 0, 0)),
+    ("&a", (0, 0, 20, 20, 0, 0, 0, 0, 0), (0, 0, 60, 60, 0, 0, 0, 0, 0)),
+    ("1f", (0, 0, 0, 0, 0, 47, 0, 0, 0), (0, 0, 0, 0, 0, 65, 0, 0, 0)),
+    ("1[", (0, 0, 0, 0, 0, 34, 0, 0, 0), (0, 0, 0, 0, 0, 56, 0, 0, 0)),
+    ("4a", (0, 0, 40, 40, 0, 0, 0, 0, 0), (0, 0, 50, 50, 0, 0, 0, 0, 0)),
+    ("#a", (0, 0, 26, 26, 0, 0, 0, 0, 0), (0, 0, 36, 36, 0, 0, 0, 0, 0)),
+    ("0i", (0, 0, 0, 0, 0, 60, 0, 0, 0), (0, 0, 0, 0, 0, 60, 0, 0, 0)),
+    ("6j", (0, 0, 0, 0, 0, 60, 0, 0, 0), (0, 0, 0, 0, 0, 60, 0, 0, 0)),
+)
+# note on the satelliteData:
+# The satellite cluster is stated. The dominant cluster can be found by looking
+# in the bin containing the satellite cluster, it will be the cluster with a
+# dominance of "dom"
+
+
+def getSatelliteInfo(name):
+  if name in satelliteTable:
+    return satelliteTable[name]
+  else:
+    return None
+
+
+# *** Cluster data  *******************************************************
+
+# The cluster data: centers of each cluster in 7 dimensions
+#   (number, name, status, color, dominance ... the 7 angles)
+bin0data = (0, "trig",
+    ( 0 , "!!", "triaged", "white      ", "tri",
+        (0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0)),
+)
+
+bin1data = (1, "33 p",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0)),
+    ( 1 , "1a", "certain", "yellowtint ", "dom",
+        (180.0,  81.495,  212.25,  288.831,  294.967,  173.99,  53.55,  81.035,  180.0)),
+    ( 2 , "1m", "certain", "blue       ", "sat",
+        (180.0,  83.513,  218.12,  291.593,  292.247,  222.3,  58.067,  86.093,  180.0)),
+    ( 3 , "1L", "certain", "green      ", "sat",
+        (180.0,  85.664,  245.014,  268.257,  303.879,  138.164,  61.95,  79.457,  180.0)),
+    ( 4 , "&a", "certain", "cyan       ", "sat",
+        (180.0,  82.112,  190.682,  264.945,  295.967,  181.839,  51.455,  81.512,  180.0)),
+    ( 5 , "7a", "certain", "pink       ", "ord",
+        (180.0,  83.414,  217.4,  222.006,  302.856,  160.719,  49.097,  82.444,  180.0)),
+    ( 6 , "3a", "certain", "magenta    ", "ord",
+        (180.0,  85.072,  216.324,  173.276,  289.32,  164.132,  45.876,  84.956,  180.0)),
+    ( 7 , "9a", "certain", "hotpink    ", "ord",
+        (180.0,  83.179,  210.347,  121.474,  288.568,  157.268,  49.347,  81.047,  180.0)),
+    ( 8 , "1g", "certain", "sea        ", "ord",
+        (180.0,  80.888,  218.636,  290.735,  167.447,  159.565,  51.326,  85.213,  180.0)),
+    ( 9 , "7d", "certain", "purple     ", "ord",
+        (180.0,  83.856,  238.75,  256.875,  69.562,  170.2,  52.8,  85.287,  180.0)),
+    ( 10 , "3d", "certain", "peach      ", "ord",
+        (180.0,  85.295,  244.085,  203.815,  65.88,  181.13,  54.68,  86.035,  180.0)),
+    ( 11 , "5d", "certain", "yellow     ", "ord",
+        (180.0,  79.671,  202.471,  63.064,  68.164,  143.45,  49.664,  82.757,  180.0)),
+    ( 12 , "3g", "wannabe", "gray       ", "ord",
+        (180.0,  84.0,  195.0,  146.0,  170.0,  170.0,  52.0,  84.0,  180.0)),
+)
+
+bin2data = (2, "33 t",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0)),
+    ( 1 , "1e", "certain", "red        ", "ord",
+        (180.0,  80.514,  200.545,  280.51,  249.314,  82.662,  167.89,  85.507,  180.0)),
+    ( 2 , "1c", "certain", "gold       ", "dom",
+        (180.0,  80.223,  196.591,  291.299,  153.06,  194.379,  179.061,  83.648,  180.0)),
+    ( 3 , "1f", "certain", "lime       ", "sat",
+        (180.0,  81.395,  203.03,  294.445,  172.195,  138.54,  175.565,  84.47,  180.0)),
+    ( 4 , "5j", "certain", "sky        ", "ord",
+        (180.0,  87.417,  223.558,  80.175,  66.667,  109.15,  176.475,  83.833,  180.0)),
+    ( 5 , "5n", "wannabe", "gray       ", "ord",
+        (180.0,  86.055,  246.502,  100.392,  73.595,  213.752,  183.395,  85.483,  180.0)),
+)
+
+bin3data = (3, "33 m",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0)),
+    # ( 1 , "!!", "nothing", "white      ", "out",
+    #     (0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0)),
+    # viewed as highly dubious, KPB 210308
+)
+
+bin4data = (4, "32 p",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+    ( 1 , "1b", "certain", "cyan       ", "dom",
+        (180.000, 084.215, 215.014, 288.672, 300.420, 177.476, 058.307, 144.841, 180.000)),
+    ( 2 , "1[", "certain", "pink       ", "sat",
+        (180.000, 082.731, 220.463, 288.665, 296.983, 221.654, 054.213, 143.771, 180.000)),
+    ( 3 , "3b", "certain", "lilac      ", "ord",
+        (180.000, 084.700, 226.400, 168.336, 292.771, 177.629, 048.629, 147.950, 180.000)),
+    ( 4 , "1z", "certain", "peach      ", "ord",
+        (180.000, 083.358, 206.042, 277.567, 195.700, 161.600, 050.750, 145.258, 180.000)),
+    ( 5 , "5z", "certain", "purple     ", "ord",
+        (180.000, 082.614, 206.440, 052.524, 163.669, 148.421, 050.176, 147.590, 180.000)),
+    ( 6 , "7p", "certain", "sea        ", "ord",
+        (180.000, 084.285, 236.600, 220.400, 068.300, 200.122, 053.693, 145.730, 180.000)),
+    ( 7 , "5p", "wannabe", "gray       ", "ord",
+        (180.000, 084.457, 213.286, 069.086, 075.500, 156.671, 057.486, 147.686, 180.000)),
+)
+
+bin5data = (5, "32 t",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+    ( 1 , "1t", "certain", "red        ", "ord",
+        (180.000, 081.200, 199.243, 288.986, 180.286, 194.743, 178.200, 147.386, 180.000)),
+    ( 2 , "5q", "certain", "yellow     ", "ord",
+        (180.000, 082.133, 204.933, 069.483, 063.417, 115.233, 176.283, 145.733, 180.000)),
+)
+
+bin6data = (6, "32 m",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+    ( 1 , "1o", "certain", "sky        ", "ord",
+        (180.000, 083.977, 216.508, 287.192, 297.254, 225.154, 293.738, 150.677, 180.000)),
+    ( 2 , "7r", "certain", "lilactint  ", "ord",
+        (180.000, 084.606, 232.856, 248.125, 063.269, 181.975, 295.744, 149.744, 180.000)),
+    ( 3 , "5r", "wannabe", "gray       ", "ord",
+        (180.000, 083.000, 196.900, 065.350, 060.150, 138.425, 292.550, 154.275, 180.000)),
+)
+
+bin7data = (7, "23 p",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+    ( 1 , "2a", "certain", "cyan       ", "ord",
+        (180.000, 145.399, 260.339, 288.756, 288.444, 192.733, 053.097, 084.067, 180.000)),
+    ( 2 , "4a", "certain", "yellow     ", "sat",
+        (180.000, 146.275, 259.783, 169.958, 298.450, 169.583, 050.908, 083.967, 180.000)),
+    ( 3 , "0a", "certain", "green      ", "dom",
+        (180.000, 149.286, 223.159, 139.421, 284.559, 158.107, 047.900, 084.424, 180.000)),
+    ( 4 , "#a", "certain", "hotpink    ", "sat",
+        (180.000, 148.006, 191.944, 146.231, 289.288, 150.781, 042.419, 084.956, 180.000)),
+    ( 5 , "4g", "certain", "greentint  ", "ord",
+        (180.000, 148.028, 256.922, 165.194, 204.961, 165.194, 049.383, 082.983, 180.000)),
+    ( 6 , "6g", "certain", "gold       ", "ord",
+        (180.000, 145.337, 262.869, 079.588, 203.863, 189.688, 058.000, 084.900, 180.000)),
+    ( 7 , "8d", "certain", "red        ", "ord",
+        (180.000, 148.992, 270.596, 240.892, 062.225, 176.271, 053.600, 087.262, 180.000)),
+    ( 8 , "4d", "certain", "sky        ", "ord",
+        (180.000, 149.822, 249.956, 187.678, 080.433, 198.133, 061.000, 089.378, 180.000)),
+    ( 9 , "6d", "certain", "orange     ", "ord",
+        (180.000, 146.922, 241.222, 088.894, 059.344, 160.683, 052.333, 083.417, 180.000)),
+    ( 10 , "2g", "wannabe", "gray       ", "ord",
+        (180.000, 141.900, 258.383, 286.517, 178.267, 165.217, 048.350, 084.783, 180.000)),
+)
+
+bin8data = (8, "23 t",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+    ( 1 , "2h", "certain", "sea        ", "ord",
+        (180.000, 147.782, 260.712, 290.424, 296.200, 177.282, 175.594, 086.565, 180.000)),
+    ( 2 , "4n", "certain", "peach      ", "ord",
+        (180.000, 143.722, 227.256, 203.789, 073.856, 216.733, 194.444, 080.911, 180.000)),
+    ( 3 , "0i", "certain", "lilactint  ", "sat",
+        (180.000, 148.717, 274.683, 100.283, 080.600, 248.133, 181.817, 082.600, 180.000)),
+    ( 4 , "6n", "certain", "lilac      ", "dom",
+        (180.000, 150.311, 268.383, 084.972, 063.811, 191.483, 176.644, 085.600, 180.000)),
+    ( 5 , "6j", "certain", "purple     ", "sat",
+        (180.000, 141.633, 244.100, 066.056, 071.667, 122.167, 182.200, 083.622, 180.000)),
+)
+
+bin9data = (9, "23 m",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+    ( 1 , "0k", "wannabe", "gray       ", "ord",
+        (180.000, 149.070, 249.780, 111.520, 278.370, 207.780, 287.820, 086.650, 180.000)),
+)
+
+bin10data = (10, "22 p",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+    ( 1 , "2[", "certain", "sea        ", "ord",
+        (180.000, 146.383, 259.402, 291.275, 291.982, 210.048, 054.412, 147.760, 180.000)),
+    ( 2 , "4b", "certain", "gold       ", "ord",
+        (180.000, 145.256, 244.622, 162.822, 294.159, 171.630, 045.900, 145.804, 180.000)),
+    ( 3 , "0b", "certain", "red        ", "ord",
+        (180.000, 147.593, 248.421, 112.086, 274.943, 164.764, 056.843, 146.264, 180.000)),
+    ( 4 , "4p", "certain", "purple     ", "ord",
+        (180.000, 150.077, 260.246, 213.785, 071.900, 207.638, 056.715, 148.131, 180.000)),
+    ( 5 , "6p", "certain", "sky        ", "ord",
+        (180.000, 146.415, 257.831, 089.597, 067.923, 173.051, 055.513, 147.623, 180.000)),
+    ( 6 , "2z", "wannabe", "gray       ", "ord",
+        (180.000, 142.900, 236.550, 268.800, 180.783, 185.133, 054.467, 143.350, 180.000)),
+)
+
+bin11data = (11, "22 t",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+    ( 1 , "4s", "certain", "lime       ", "ord",
+        (180.000, 149.863, 247.562, 170.488, 277.938, 084.425, 176.413, 148.087, 180.000)),
+    ( 2 , "2u", "wannabe", "gray       ", "ord",
+        (180.000, 143.940, 258.200, 298.240, 279.640, 183.680, 183.080, 145.120, 180.000)),
+)
+
+bin12data = (12, "22 m",
+    ( 0 , "!!", "outlier", "white      ", "out",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+    ( 1 , "2o", "certain", "hotpink    ", "ord",
+        (180.000, 147.342, 256.475, 295.508, 287.408, 194.525, 293.725, 150.458, 180.000)),
+)
+
+bin13data = (13, "inc ",
+    ( 0 , "__", "incompl", "white      ", "inc",
+        (000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000, 000.000)),
+)
+
+
+def buildBin(data):
+    ordinal = data[0]
+    name = data[1]
+    clusters = []
+    for item in data[2:]:
+        c = suitenamedefs.Cluster(*item)
+        if c.dominance == "sat":
+            c.satelliteInfo = getSatelliteInfo(c.name)
+        clusters.append(c)
+    bin = Bin(ordinal, name, clusters)
+    return bin
+
+
+# The bins become an associative table that maps the bin selector information
+# directly to the bin
+# selector = (puckerdm, puckerd, gammaname)
+# bins 0 and 13 are catchbasins for outliers, they are not indexed by selectors
+def buildTheBins():
+    bins = {}
+    bins[0] = buildBin(bin0data)
+    bins[(3, 3, "p")] = buildBin(bin1data)
+    bins[(3, 3, "t")] = buildBin(bin2data)
+    bins[(3, 3, "m")] = buildBin(bin3data)
+    bins[(3, 2, "p")] = buildBin(bin4data)
+    bins[(3, 2, "t")] = buildBin(bin5data)
+    bins[(3, 2, "m")] = buildBin(bin6data)
+    bins[(2, 3, "p")] = buildBin(bin7data)
+    bins[(2, 3, "t")] = buildBin(bin8data)
+    bins[(2, 3, "m")] = buildBin(bin9data)
+    bins[(2, 2, "p")] = buildBin(bin10data)
+    bins[(2, 2, "t")] = buildBin(bin11data)
+    bins[(2, 2, "m")] = buildBin(bin12data)
+    bins[13] = buildBin(bin13data)
+
+    # build aliases so that bins can be indexed by number during output
+    bins[1] = bins[(3, 3, "p")]
+    bins[2] = bins[(3, 3, "t")]
+    bins[3] = bins[(3, 3, "m")]
+    bins[4] = bins[(3, 2, "p")]
+    bins[5] = bins[(3, 2, "t")]
+    bins[6] = bins[(3, 2, "m")]
+    bins[7] = bins[(2, 3, "p")]
+    bins[8] = bins[(2, 3, "t")]
+    bins[9] = bins[(2, 3, "m")]
+    bins[10] = bins[(2, 2, "p")]
+    bins[11] = bins[(2, 2, "t")]
+    bins[12] = bins[(2, 2, "m")]
+    return bins
+
+
+# parseCommandLine()
+buildSatelliteTable()
+bins = buildTheBins()
diff --git a/mmtbx/suitename/suiteninput.py b/mmtbx/suitename/suiteninput.py
new file mode 100644
index 0000000000..bac62f1691
--- /dev/null
+++ b/mmtbx/suitename/suiteninput.py
@@ -0,0 +1,219 @@
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+from suitenamedefs import Suite, Residue, findBase, globals
+
+"""
+This module handles reading suites from "dangle" format files
+reading residues from kinemage format files and regrouping them into suites.
+Extraction of suites from loaded cctbx models is handled elsewhere.
+"""
+
+#        Copyright 2021  Richardson Lab at Duke University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import numpy as np
+import math, sys
+
+
+def stringToFloat(string):
+  try:
+    n = float(string)
+  except ValueError:
+    n = 9999.0  # or maybe math.nan?
+  return n
+
+
+def readResidues(inFile):
+  options = globals.options
+
+  lines = inFile.readlines()
+  if not lines:
+    return []
+  # catch a specific error notation from mp_geo:
+  while lines[0].strip().startswith("Atom pair"):
+    del lines[0]
+  residues = []
+  try:
+    i = 0
+    line = ""
+    for line in lines:
+      i += 1
+      if len(line.strip()) == 0 or line[0] == "#":  # blank or comment line
+        continue
+      fields = line.split(":")
+      ids = fields[: options.pointidfields]
+      baseCode = fields[options.pointidfields - 1]
+      angleStrings = fields[options.pointidfields :]
+      if (
+          ids[options.altidfield-1].strip() != ""
+          and ids[options.altidfield-1] != options.altid
+      ):    # -1 converts 1-based to 0-based counting
+        continue  # lines for the wrong alternative conformation are ignored
+
+      base = findBase(baseCode)
+      if not base:  # ignore DNA bases
+        continue
+      angles = np.array([stringToFloat(s) for s in angleStrings])
+      for i in range(len(angles)):
+        if angles[i] < 0:
+          angles[i] += 360.0
+      residue = Residue(ids, base, angles)
+      residues.append(residue)
+  except IndexError:
+    print("Suitename found malformed input on line {}, reading no further:".format(i),
+        file=sys.stderr)
+    print("    ", line, file=sys.stderr)
+  return residues
+
+
+def readKinemageFile(inFile):
+  """
+  We glean the following information from a kinemage file:
+  The @dimension command gives us the number of dimensions in the data
+  Anything between a @balllist command and a subsequent @ command
+  is a data line.
+  """
+  options = globals.options
+
+  lines = inFile.readlines()
+  goodLines = []
+  place, line = findPrefixInList(lines, "@dimension")
+  if place > 0:
+    items = line.split()
+    dimension = len(items) - 1
+  else:
+    dimension = options.anglefields
+    place = 0
+  while place >= 0:
+    begin, line = findPrefixesInList(lines, "@balllist", "@dotlist", place)
+    if begin > 0:
+      end, line = findPrefixInList(lines, "@", begin + 1)
+      place = end
+      if end < 0:
+        end = len(lines)
+      goodLines += lines[begin + 1 : end]
+    else:
+      break
+  if len(goodLines) == 0:
+    goodLines = lines  # assume a pure data file
+  return readKinemageSuites(goodLines, dimension)
+
+
+def readKinemageSuites(lines, dimension):
+  """Read a list of kinemage data lines to yield a suite."""
+  suites = []
+  for line in lines:
+    if len(line.strip()) == 0 or line[0] == "#":  # blank or comment line
+      continue
+    # A meaningful line begins with an id string enclosed in braces
+    if line[0] == "{":
+      mark = line.find("}")
+      if mark > 0:
+        idString = line[1:mark]
+        ids = idString.split(":")
+
+        # there may be some miscellaneous markers after the id string
+        k = mark + 1
+        while k < len(line) and not line[k].isdigit():
+          k = k + 1
+        mark2 = k
+
+        # once we see a number, everything else is angles
+        angleText = line[mark2:]
+        angleStrings = angleText.split(" ")
+        angleStrings2 = angleText.split(",")
+        if len(angleStrings2) > len(angleStrings):
+          angleStrings = angleStrings2
+        angleList = [stringToFloat(s) for s in angleStrings]
+        if len(angleList) != dimension:
+          continue  # wrong number of dimensions means probably not a data point
+        if dimension == 9:
+          angles = np.array(angleList)
+        else:  # given only 7 angles,skipping the chi angles on the ends
+          angles = np.array([180.0] + angleList + [180.0])
+        for i in range(len(angles)):
+          if angles[i] < 0:
+            angles[i] += 360.0
+
+        suite = Suite(ids, ids[9][2], angles)
+        suites.append(suite)
+  return suites
+
+
+def findPrefixInList(list, prefix, start=0):
+  for i, s in enumerate(list[start:]):
+      if s.startswith(prefix):
+        return i + start, s
+  return -1, None
+
+
+def findPrefixesInList(list, prefix1, prefix2, start=0):
+  for i, s in enumerate(list[start:]):
+      if s.startswith(prefix1) or s.startswith(prefix2):
+        return i + start, s
+  return -1, None
+
+
+def buildSuiteBetweenResidues(r1, r2):
+  suite = Suite(r2.pointIDs, r2.base)
+  if len(r1.angle) > 6:
+    suite.chiMinus = r1.chi
+  suite.deltaMinus = r1.delta
+  suite.epsilon = r1.epsilon
+  suite.zeta = r1.zeta
+  suite.alpha = r2.alpha
+  suite.beta = r2.beta
+  suite.gamma = r2.gamma
+  suite.delta = r2.delta
+  if len(r2.angle) > 6:
+    suite.chi = r2.chi
+  return suite
+
+
+def buildSuiteFirst(r2):
+  suite = Suite(r2.pointIDs, r2.base)
+  suite.alpha = r2.alpha
+  suite.beta = r2.beta
+  suite.gamma = r2.gamma
+  suite.delta = r2.delta
+  if len(r2.angle) > 6:
+    suite.chi = r2.chi
+  suite.epsilon = 999
+  suite.zeta = 999
+  suite.chiMinus = 999
+  suite.deltaMinus = 999
+  return suite
+
+
+def buildSuiteLast(r1):
+  suite = Suite((), "")
+  if len(r1.angle) > 6:
+    suite.chiMinus = r1.chi
+  suite.deltaMinus = r1.delta
+  suite.epsilon = r1.epsilon
+  suite.zeta = r1.zeta
+  suite.alpha = 999
+  suite.beta = 999
+  suite.gamma = 999
+  suite.delta = 999
+  suite.chi = 999
+  return suite
+
+
+def buildSuites(residues):
+  suites = [buildSuiteFirst(residues[0])]
+  for i in range(len(residues) - 1):
+    suites.append(buildSuiteBetweenResidues(residues[i], residues[i + 1]))
+  suites.append(buildSuiteLast(residues[-1]))
+  return suites
diff --git a/mmtbx/suitename/suitenout.py b/mmtbx/suitename/suitenout.py
new file mode 100644
index 0000000000..b801eb6df1
--- /dev/null
+++ b/mmtbx/suitename/suitenout.py
@@ -0,0 +1,422 @@
+#        Copyright 2021  Richardson Lab at Duke University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function, unicode_literals
+from suitenamedefs import globals, reasons
+from suiteninit import bins
+
+import numpy as np
+
+options = None  # global,  will be copied from suitename.py global
+
+reportCountAll = 0
+trigCountAll = 0
+suitenessSumAll = 0
+binSuiteCountAll = 0
+
+
+def output(outFile, suites, outNote):
+  global options
+  options = globals.options
+  for s in suites:
+    outSuite(outFile, s)
+  writeFinalOutput(outFile, suites, outNote)
+
+
+def outSuite(outFile, s):
+  if options.string:
+    string1Suite(outFile, s)
+  elif options.kinemage:
+    pass
+  else:
+    reportSuite(outFile, s)
+
+
+def writeFinalOutput(outFile, suites, outNote):
+  if options.satellites:
+    outNote.comment = " special general case satellite widths, power = 3.00"
+  else:
+    outNote.comment = " all general case widths, power = 3.00"
+  if options.string:
+    outFile.write("\n")
+  elif options.kinemage:
+    kinemageFinal(outFile, suites, outNote)
+  else:
+    reportFinal(outFile, outNote)
+
+
+def string1Suite(outFile, suite):
+  if options.nosequence:
+    basestring = ":"
+  else:
+    basestring = suite.base
+  if suite.cluster:
+    name = suite.cluster.name
+  else:
+    name = "!!"
+  outFile.write("{}{}".format(name,basestring))
+
+
+def reportSuite(outFile, suite):
+  global reportCountAll, trigCountAll, suitenessSumAll, binSuiteCountAll
+
+  if options.noinc and not suite.valid:
+    return
+
+  # 1. write one line of output for this suite
+  reason = ""; note=""
+  outIDs = ":".join(suite.pointID)
+  if suite.issue:
+    reason = " " + reasons[suite.issue]
+  elif suite.comment:
+    reason = " " + suite.comment
+  elif suite.situation and options.causes:
+    reason += " " + suite.situation
+  if suite.cluster.status == "wannabe":
+    note = " wannabe"
+  output = (
+    "{} {} {} {:5.3f}{}{}\n".format(outIDs, suite.bin.name, suite.cluster.name, float(suite.suiteness), reason, note)
+  )
+  outFile.write(output)
+
+  # 2. gather statistics
+  reportCountAll += 1
+  bin = suite.bin
+  cluster = suite.cluster
+  suiteness = suite.suiteness
+  if bin.ordinal == 0:
+    trigCountAll += 1
+  elif bin.ordinal < 13:
+    suitenessSumAll += suiteness
+    binSuiteCountAll += 1
+
+  if cluster.ordinal == 0:
+    cluster.suitenessCounts[11] += 1
+  else:
+    cluster.suitenessSum += suiteness
+    # report in statistical baskets at intervals of 0.1:
+    # everything from 0 to 0.1 goes in bucket 1
+    # ... everything from 0.9 to 1.o goes into bucket 10
+    if suiteness == 0:
+      bucket = 0
+    else:
+      bucket = 1 + int(suiteness * 10)
+    cluster.suitenessCounts[bucket] += 1
+
+
+def reportFinal(outFile, outNote):
+  if not options.chart:
+    outFile.write(outNote.comment + "\n")
+    suitenessAverage(outFile, 0)
+    if bins[1].cluster[1].count > 0:  # Aform 1a    070325
+      suitenessAverage(outFile, 1)
+      suitenessAverage(outFile, 2)
+
+
+def suitenessAverage(outFile,  mode):
+  """
+  Gather statistics on suiteness.
+  12 buckets:
+    one for suiteness=0
+    one each for divisions by tenths from 0 to 1
+    one for outliers
+  """
+  bucket = np.zeros(12, dtype=int)
+  sum = 0
+  average = 0
+  allCount = 0
+  excludedCluster = None
+
+  if mode == 1:
+    # cluster 1a all by itself
+    comment = " A form (1a)"
+    cluster = bins[1].cluster[1]
+    sum = cluster.suitenessSum
+    for k in range(12):
+      bucket[k] = cluster.suitenessCounts[k]
+  else:
+    if mode == 0:
+      # all complete suites
+      startWith = 0
+      comment = "For all"
+      outFile.write(
+          "Found {} complete suites derived from {} entries\n".format(
+          binSuiteCountAll + trigCountAll, reportCountAll
+        )
+      )
+      outFile.write(
+          "{} suites were triaged, leaving {} assigned to bins\n".format(
+          trigCountAll, binSuiteCountAll
+        )
+      )
+    elif mode == 2:
+      # all complete suites except cluster 1a
+      startWith = 1
+      comment = " non-1a  has"
+      excludedCluster = bins[1].cluster[1]
+
+    for i in range(1, 13):
+      bin = bins[i]
+      # the final bin 13 is for the pseudo-suites with incomplete
+      # angles and is ignored
+      for cluster in bin.cluster[startWith:]:
+        # cluster 0 in every bin is for outliers and is ignored
+        if cluster is excludedCluster:
+          continue  # mode 2: ignore cluster 1a
+        sum += cluster.suitenessSum
+        for k in range(12):
+          bucket[k] += cluster.suitenessCounts[k]
+
+  allCount = np.sum(bucket)
+  if allCount > 1:
+    average = sum / allCount
+  else:
+    average = 0
+  outFile.write(
+    "{} {} suites: average suiteness== {:5.3f} (power==3.00)\n".format(
+        comment, allCount, average
+    )
+  )
+  if mode == 0:
+    outFile.write("{:6d} suites are  outliers\n".format(bucket[11]))
+  outFile.write("{:6d} suites have suiteness == 0    \n".format(bucket[0]))
+  outFile.write("{:6d} suites have suiteness >  0 <.1\n".format(bucket[1]))
+  for k in range(2, 10):
+    outFile.write(
+      "{:6d} suites have suiteness >=.{} <.{}\n".format(bucket[k], k - 1, k)
+    )
+  outFile.write("{:6d} suites have suiteness >=.9    \n".format(bucket[10]))
+
+
+def clearStatistics():
+  global reportCountAll, trigCountAll, suitenessSumAll, binSuiteCountAll
+  reportCountAll = 0
+  trigCountAll = 0
+  suitenessSumAll = 0
+  binSuiteCountAll = 0
+
+  for bin in bins.values():
+    for cluster in bin.cluster:
+      cluster.count = 0
+      cluster.suitenessSum = 0
+      for k in range(12):
+        cluster.suitenessCounts[k] = 0
+
+
+# ***** kinemage output format *****************************************************
+
+def kinemage1Suite(suite, bin, cluster, notes, distance, suiteness, issue, comment,
+                    pointMaster, pointColor):
+    suite.pointMaster = pointMaster
+    suite.pointColor = pointColor
+    suite.notes = notes
+
+
+# static text: viewing parameters
+janesviews = """
+@viewid {d e z}
+@zoom 1.00
+@zslab 200
+@ztran 0
+@center 197.500 172.300 178.300
+@axischoice 2 3 4
+@matrix
+0.07196 0.11701 -0.99052 -0.00336 0.99312 0.11707 0.99740 -0.00509 0.07186
+@2viewid {zag front}
+@2zoom 1.00
+@2zslab 200
+@2ztran 0
+@2center 174.091 194.887 207.768
+@2axischoice 4 5 7
+@2matrix
+0.99508 -0.00018 -0.09905 -0.00135 -0.99993 -0.01172 -0.09904 0.0118 -0.99501
+@3viewid {a b g}
+@3zoom 1.00
+@3zslab 200
+@3ztran 0
+@3center 175.700 189.600 64.100
+@3axischoice 5 6 7
+@3matrix
+0.99955 0.000101 0.030002 0.0002 0.99995 -0.010012 -0.030001 0.010013 0.9995
+
+"""
+
+# static text: static items in the display
+kinemageFrame = ("\n"+
+  "@group {frame} dominant \n"
+  "@vectorlist {frame} color= white \n"
+  "P   0.000   0.000   0.000   5.000   0.000   0.000 \n"+
+  "P  35.000   0.000   0.000  40.000   0.000   0.000 \n"+
+  "P  80.000   0.000   0.000 160.000   0.000   0.000 \n"+
+  "P   0.000   0.000   0.000   0.000   5.000   0.000 \n"+
+  "P   0.000  35.000   0.000   0.000  40.000   0.000 \n"+
+  "P   0.000  80.000   0.000   0.000 160.000   0.000 \n"+
+  "P   0.000   0.000   0.000   0.000   0.000   5.000 \n"+
+  "P   0.000   0.000  35.000   0.000   0.000  40.000 \n"+
+  "P   0.000   0.000  80.000   0.000   0.000 160.000 \n"+
+  "P 200.000   0.000   0.000 280.000   0.000   0.000 \n"+
+  "P 320.000   0.000   0.000 360.000   0.000   0.000 \n"+
+  "P   0.000 200.000   0.000   0.000 280.000   0.000 \n"+
+  "P   0.000 320.000   0.000   0.000 360.000   0.000 \n"+
+  "P   0.000   0.000 200.000   0.000   0.000 280.000 \n"+
+  "P   0.000   0.000 320.000   0.000   0.000 360.000 \n"+
+  "@labellist {XYZ} color= white \n"+
+  "{X}  20.000  -5.000  -5.000 \n"+
+  "{X} 380.000  -5.000  -5.000 \n"+
+  "{Y}  -5.000  20.000  -5.000 \n"+
+  "{Y}  -5.000 380.000  -5.000 \n"+
+  "{Z}  -5.000  -5.000  20.000 \n"+
+  "{Z}  -5.000  -5.000 380.000 \n"+
+  "@labellist {mtp} color= green \n"+
+  "{p}  60.000   0.000   0.000 \n"+
+  "{t} 180.000   0.000   0.000 \n"+
+  "{m} 300.000   0.000   0.000 \n"+
+  "{p}   0.000  60.000   0.000 \n"+
+  "{t}   0.000 180.000   0.000 \n"+
+  "{m}   0.000 300.000   0.000 \n"+
+  "{p}   0.000   0.000  60.000 \n"+
+  "{t}   0.000   0.000 180.000 \n"+
+  "{m}   0.000   0.000 300.000 \n"+
+  "\n"
+)
+
+def kinemageFinal(outFile, suites, outNote):
+  """
+  Output the content of a kinemage file
+  The 3, 2 order may seem odd, but it is a standard
+  """
+  if not suites or all([not s.valid for s in suites]):
+    # empty input gives empty output
+    return
+  kinemageHeader(outFile, outNote)
+  for deltaMinus in (3, 2):
+    for delta in (3, 2):
+      binGroupOut(outFile, deltaMinus, delta, suites)
+  triaged = bins[0]
+  if triaged.count > 0:
+    outFile.write("@group {triaged} dominant dimension=9 wrap=360 select off\n")
+    binOut(outFile, triaged, suites)
+
+
+def kinemageHeader(outFile, outNote):
+  """ The invariant portion of a kinemage file """
+  outFile.write("@text\n {}\n {}\n".format(outNote.version, outNote.comment))
+  outFile.write("@kinemage 1\n")
+  outFile.write("@onewidth\n")
+  if options.etatheta: # 070524
+    outFile.write(
+        "@dimension {theta} {delta-1} {epsilon-1} {zeta-1} {alpha} "
+        "{beta} {gamma} {delta} {eta}\n"
+    )
+  else:
+    outFile.write(
+        "@dimension {chi-1} {delta-1} {epsilon-1} {zeta-1} {alpha} "
+        "{beta} {gamma} {delta} {chi}\n"
+    )
+  outFile.write(
+      "@dimminmax 0.000 360.000 0.000 360.000 0.000 360.000 0.000 "
+      "360.000 0.000 360.000 0.000 360.000 0.000 360.000 0.000 360.000 "
+      "0.000 360.000\n"
+    )
+  if outNote.outliers:
+    outFile.write("@pointmaster 'O' {outliers}\n")
+  if outNote.wannabes:
+    outFile.write("@master {wannabees}\n")
+  outFile.write(janesviews)
+  outFile.write(kinemageFrame)
+
+
+def binGroupOut(outFile, deltaMinus, delta, suites):
+  # If any bin in the group has data, generate a group header
+  groupCount = 0
+  for gamma in ("p", "t", "m"):
+    bin = bins[(deltaMinus, delta, gamma)]
+    groupCount += bin.count
+  if groupCount > 0:
+    outFile.write(
+      "@group {{{}{}}} recessiveon dimension=9".format(deltaMinus,delta)+
+      " wrap=360 select animate off\n")
+
+  # generate the data
+  for gamma in ("p", "t", "m"):
+    bin = bins[(deltaMinus, delta, gamma)]
+    if bin.count > 0:
+      binOut(outFile, bin, suites)
+
+
+def binOut(outFile, bin, suites):
+  if any([c.count > 0 for c in bin.cluster]):
+    outFile.write("@subgroup {{{}}} recessiveon \n".format(bin.name))
+  for cluster in bin.cluster[1:]:
+      # the first cluster, for outliers, will be handled later
+      if cluster.count > 0:  # empty clusters vanish
+        extras = ""
+        if cluster.status == "wannabe":
+          extras = " master= {wannabees}"
+        # display a ball for each for each suite in this cluster
+        ballList = (
+            "@balllist {{{} {}}} color= {} radius= 1 "
+            "nohilite master= {{data}}{}\n"
+        ).format(bin.name, cluster.name, cluster.clusterColor, extras)
+        # display a ring surrounding the center of the cluster
+        ringList = (
+            "@ringlist {{{} {}}} color= {} radius= 10 width= 1 "
+            "nobutton master= {{avsigma}}{}\n"
+        ).format(bin.name, cluster.name, cluster.clusterColor, extras)
+        angleList = formatAngles(cluster.angle[1:-1], ' ')
+        ringList2 = "{{{} {}}} 180 {} 180\n".format(
+            bin.name, cluster.name, angleList
+        )
+        labelList = (
+            "@labellist {{{} {}}} color= {} nobutton "
+            "master= {{labels}}{}\n"
+        ).format(bin.name, cluster.name, cluster.clusterColor, extras)
+        outFile.write(ballList)
+        outPoints(outFile, bin, cluster, suites, "")
+        outFile.write(ringList)
+        outFile.write(ringList2)
+        outFile.write(labelList)
+        outFile.write(ringList2)
+
+  # handle outliers if there are any:
+  if bin.cluster[0].count > 0:
+    cluster = bin.cluster[0]
+    ballList = (
+      "@balllist {{{} {}}} color= {} radius= 1 "
+      "nohilite master= {{data}}\n"
+    ).format(bin.name, cluster.name, cluster.clusterColor)
+    outFile.write(ballList)
+    outPoints(outFile, bin, cluster, suites, "'O' white")
+
+
+def outPoints(outFile, bin, cluster, suites, extra1):
+  for s in suites:
+    if s.cluster is cluster:
+      if bin.name == "trig":  # the triage bin is specially handled
+        extra = "'{}'".format(s.pointMaster)
+      else:
+        extra = extra1
+      ids = ':'.join(s.pointID)
+      line = \
+        ("{{{} {} :D=={:5.3f}".format(bin.name, cluster.name, s.distance) \
+       + ":S=={:5.3f}: {}}} {} ,".format(s.suiteness, ids, extra)) \
+       + formatAngles(s.angle, ",") + "\n"
+      outFile.write(line)
+
+
+def formatAngles(angles, separator):
+    strings = ["{:7.2f}".format(a) for a in angles]
+    out = separator.join(strings)
+    return out
diff --git a/mmtbx/suitename/suites.py b/mmtbx/suitename/suites.py
new file mode 100644
index 0000000000..673c1874fc
--- /dev/null
+++ b/mmtbx/suitename/suites.py
@@ -0,0 +1,186 @@
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+import sys, os
+
+#        Copyright 2021  Richardson Lab at Duke University
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from suitenamedefs import globals
+from iotbx.data_manager import DataManager    #   Load in the DataManager
+from libtbx import phil
+from libtbx.utils import Sorry
+# from mmtbx.validation import utils
+# from cctbx import geometry_restraints
+# from collections import defaultdict
+from diangle import getResidueDihedrals
+
+
+# IMPORT TO EXPORT:
+from mmtbx.suitename.suitename import compute, write, \
+    finalStats, clearStats
+
+
+# The following are the options available, in Phil format,
+# for human and computer comprehension.
+philOptions = """
+  suitename {
+    # input
+      infile=""
+        .type=str
+        .help="the file to process"
+      anglefields = 9
+        .type=int
+        .help="number of angle fields provided, for textual input only"
+      pointidfields = 7
+        .type=int
+        .help="number of point id fields before the angle fields"
+      ptid=0
+        .type=int
+        .help="number of point id fields before the angle fields"
+      residuein=false
+        .type=bool
+        .help="expect dangle format giving residues"
+      suitein=false
+        .type=bool
+        .help="expect kinemage format giving suites directly"
+    # output
+      string=False
+        .type=bool
+        .help="output in string format, 3 characters per suite"
+      kinemage=False
+        .type=bool
+        .help="output in kinemage format, useful for visualization"
+      report=true
+        .type=bool
+        .help="output as a report, giving statistical details"
+      chart=False
+        .type=bool
+        .help="modifier to standard report, output without statistical summary"
+      nosequence = False
+        .type=bool
+        .help="modifier to string format, do not include base letters"
+      causes=False
+        .type=bool
+        .help="output extra details concerning the causes of each assignment made"
+      test=False
+        .type=bool
+        .help="display a lat of additional information about program internals"
+    # compute
+      satellites=False
+        .type=bool
+        .help="use the special satelliteWidths values for satellites"
+      nowannabe=False
+        .type=bool
+        .help="do not consider 'wannabe' clusters"
+      noinc=False
+        .type=bool
+        .help="do not display incomplete suites"
+      etatheta=False
+        .type=bool
+      altid="A"
+        .type=str
+        .help="which alternate conformer to use (A, B, etc)"
+      altidfield = 6
+        .type=int
+        .help="which field (1-based) gives the alternate conformer code"
+      version=false
+        .type=bool
+        .help="give the version number of suite name"
+    # deprecated and automatically true:
+      oneline=false
+        .type=bool
+    }
+"""
+
+def main(options, outFile=None, errorFile=None):
+  """The main track for handling PDB and CIF input formats, which will involve
+     parsing the model hierarchy to get the dihedral angles for ourselves"""
+  setOptions(options)
+  import suiteninput  # must be AFTER setOptions
+
+  if not outFile: outFile = sys.stdout
+  if not errorFile: errorFile = sys.stderr
+  inFile = options.infile
+  model = loadModel(inFile)
+
+  residues = getResidueDihedrals(model, options.altid,
+                                 name=os.path.splitext(inFile)[0],
+                                 errorFile=errorFile)
+  ### to print mp_geo-like output:
+  # for r in residues:
+  #   print(residueString(r))
+  # useful for seeing what suites were generated
+
+  if residues is not None and len(residues) > 0:
+    suiteList = suiteninput.buildSuites(residues)
+    suiteList = suiteList[:-1]
+    suiteList = compute(suiteList)
+    finalStats()
+    write(outFile, suiteList)
+    clearStats()
+
+
+def parseOptions(optionString, errorFile=None):
+  """ Use optionString to modify the defaults given in philOptions above.
+  Returns a Python object that has an attribute for every option listed
+  in philOptions.  Example: "chart=true noinc=true causes=true"
+  The values in optionString are case insensitive.
+  """
+  opt2 = """  # use this for more complex option types e.g. multiples
+  suitename {
+    report=true
+    chart=true
+}  """
+  # user_phil = phil.parse(opt2)
+
+  master_phil = phil.parse(philOptions)
+  interp = master_phil.command_line_argument_interpreter()
+  optionList = optionString.split()
+  try:
+    user_phil = interp.process(args=optionList)
+  except Sorry as e:
+    if errorFile is None:  errorFile = sys.stderr
+    print(e, file=errorFile)
+
+  working_phil = master_phil.fetch(sources=user_phil)
+  full_options = working_phil.extract()
+  return full_options.suitename
+
+
+def setOptions(optionsIn):
+  """optionsIn may be the result of parseOptions above
+     or the result of an argparse parse_args operation"""
+  from mmtbx.suitename.suitename import loadOptions
+  globals.options = optionsIn
+  loadOptions(optionsIn)
+
+
+def loadModel(filename):
+  dm = DataManager()             #   Initialize the DataManager and call it dm
+  dm.set_overwrite(True)         #   tell the DataManager to overwrite files with the same name
+  #print("Reading file")
+  model = dm.get_model(filename)
+  #print("Processing model")
+  #model.process_input_model(make_restraints=True)
+  # removed because Restraints Manager will not operate
+  # on unfamiliar residues  KPB 6/10/2021
+  return model
+
+
+def testResidues(model):
+  #print("computing dihedrals")
+  residues = getResidueDihedrals(model)
+  for r in residues:
+    print(r.pointIDs, " : ", r.angle)
+
diff --git a/mmtbx/suitename/suites.rtf b/mmtbx/suitename/suites.rtf
new file mode 100644
index 0000000000..4bdf4162fe
Binary files /dev/null and b/mmtbx/suitename/suites.rtf differ
diff --git a/mmtbx/suitename/transdent.py b/mmtbx/suitename/transdent.py
new file mode 100644
index 0000000000..812c7d69a7
--- /dev/null
+++ b/mmtbx/suitename/transdent.py
@@ -0,0 +1,103 @@
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+import re, os, sys
+
+N = 2   # indentation standard
+white = "([ \t]*)"
+simple = re.compile(".*(#|:)")
+stringy = re.compile("[^\"#:]*(\".*\"|'.*'|#|:[ \t]*(#.*)?$)")
+
+class Line(object):
+    content = ""
+    level = 0
+    actual = 0  # the actual whitespace at the front
+
+
+def analyzeWhitespace(white):
+    count = 0
+    for char in white:
+        if char == "\t":
+            count += N
+        else:
+            count += 1
+    return count
+
+
+def analyzeLine(raw):
+    # discover indentation of raw line and whether it ends with a colon
+    content = raw.strip()
+    match1 = re.match(white, raw)
+    mark = match1.end()
+    spaces = analyzeWhitespace(match1[1])
+    if content == "":  # blank line
+        spaces = 0
+    while True:
+        # this loop is exited only by returning
+        match2 = stringy.match(raw, mark)
+        if match2 is None:
+            return content, spaces, False
+
+        what = match2.group(1)
+        char = what[0]
+        if char == "#":
+            return content, spaces, False
+        elif char == ":":
+            return content, spaces, True
+        elif char == '"' or char == "'":
+            mark = match2.end()
+
+
+def main():
+    file = sys.argv[1] + ".py"
+    #root, ext = os.path.splitext(file)
+    outFile = sys.argv[1] + ".ind.py"
+    # file = "target.py"
+    stream = open(file)
+    lines = stream.readlines()
+    oStream = open(outFile, "w")
+    #oStream = sys.stderr
+
+    indents = []
+    prevSpaces = 0
+    level = 0
+    newLevel = False
+    freeze = False
+    count=-1
+    for rawLine in lines:
+        count += 1
+        content, spaces, indenting = analyzeLine(rawLine)
+        if spaces == 0 and rawLine.startswith("def "):
+          # sanity reset: no matter what, restart at left margin
+          # when defining a top level function
+          level = 0; newLevel = False
+        elif newLevel:
+          if spaces > prevSpaces:
+            # A true syntactic indent (else a stray colon, ignore)
+            level += 1
+            indents.append(prevSpaces)
+          newLevel = False
+        elif spaces > prevSpaces:
+          # indent without : is a line-continuation or data-formatting situation
+          freeze = True
+          indents.append(prevSpaces)
+        if content == "":  # blank lines influence nothing
+          print("", file=oStream)
+          continue
+#        if spaces < prevSpaces and level > 0:
+        if spaces < prevSpaces and spaces <= indents[-1]:
+            oldSpaces = indents.pop()
+            while spaces < oldSpaces:
+                oldSpaces = indents.pop()
+            level = len(indents)
+            freeze = False
+        prevSpaces = spaces
+        if freeze:
+          print((level * N + (spaces-indents[-1])) * " " + content, file=oStream)
+        else:
+          print(level * N * " " + content, file=oStream)
+        if indenting:
+            newLevel = True
+    oStream.close()
+
+main()
+


Test patch:
diff --git a/mmtbx/run_tests.py b/mmtbx/run_tests.py
index 5badc565df..59e664aeae 100644
--- a/mmtbx/run_tests.py
+++ b/mmtbx/run_tests.py
@@ -5,6 +5,8 @@
 import libtbx.load_env
 
 general_tests = [
+  "$D/suitename/unit-test/FullTest.py",
+  #
   "$D/regression/tst_sampled_model_density.py",
   # pair interaction
   "$D/pair_interaction/tst_00.py",
diff --git a/mmtbx/suitename/unit-test/FullTest.py b/mmtbx/suitename/unit-test/FullTest.py
new file mode 100644
index 0000000000..efe89d8ce7
--- /dev/null
+++ b/mmtbx/suitename/unit-test/FullTest.py
@@ -0,0 +1,16 @@
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function, unicode_literals
+import UnitTest, regression
+import os, sys
+
+workDir = os.path.dirname(sys.argv[0]).strip()
+if workDir:
+  os.chdir(workDir)
+  # else we must be there already
+
+def fullTest():
+  UnitTest.testAll()
+  regression.test()
+
+if __name__ == '__main__':
+  fullTest()
diff --git a/mmtbx/suitename/unit-test/UnitTest.py b/mmtbx/suitename/unit-test/UnitTest.py
new file mode 100644
index 0000000000..05dbe4fce3
--- /dev/null
+++ b/mmtbx/suitename/unit-test/UnitTest.py
@@ -0,0 +1,332 @@
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function
+import os, sys, inspect
+from io import StringIO
+currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
+parentdir = os.path.dirname(currentdir)
+sys.path.insert(0, parentdir)
+
+# here, before importing suitename, is an opportunity to set command line
+# options and to redirect output.
+# sys.argv.extend(["--noinc", "--chart", "--causes"])
+# sys.stdout = open("output.txt", "w")
+import suitename, suites, regression
+
+
+# The test data:
+#
+# The first portion is to test membership in every cluster.
+# The second portion is to test triage by out of bounds angle
+# for each dihedral angle. The unnecessary part of each residue is damaged so
+# the extra suite generated will be incomplete and will vanish.
+input1=u''' :1a: : : : :  Z:  9999.000: 9999.000: 9999.000:   81.495:  212.250:  288.831:  180.000
+ :1a: : : : :  Z:   294.967:  173.990:   53.550:   81.035: 9999.000: 9999.000:  180.000
+ :1m: : : : :  Z:  9999.000: 9999.000: 9999.000:   83.513:  218.120:  291.593:  180.000
+ :1m: : : : :  Z:   292.247:  222.300:   58.067:   86.093: 9999.000: 9999.000:  180.000
+ :1L: : : : :  Z:  9999.000: 9999.000: 9999.000:   85.664:  245.014:  268.257:  180.000
+ :1L: : : : :  Z:   303.879:  138.164:   61.950:   79.457: 9999.000: 9999.000:  180.000
+ :&a: : : : :  Z:  9999.000: 9999.000: 9999.000:   82.112:  190.682:  264.945:  180.000
+ :&a: : : : :  Z:   295.967:  181.839:   51.455:   81.512: 9999.000: 9999.000:  180.000
+ :7a: : : : :  Z:  9999.000: 9999.000: 9999.000:   83.414:  217.400:  222.006:  180.000
+ :7a: : : : :  Z:   302.856:  160.719:   49.097:   82.444: 9999.000: 9999.000:  180.000
+ :3a: : : : :  Z:  9999.000: 9999.000: 9999.000:   85.072:  216.324:  173.276:  180.000
+ :3a: : : : :  Z:   289.320:  164.132:   45.876:   84.956: 9999.000: 9999.000:  180.000
+ :9a: : : : :  Z:  9999.000: 9999.000: 9999.000:   83.179:  210.347:  121.474:  180.000
+ :9a: : : : :  Z:   288.568:  157.268:   49.347:   81.047: 9999.000: 9999.000:  180.000
+ :1g: : : : :  Z:  9999.000: 9999.000: 9999.000:   80.888:  218.636:  290.735:  180.000
+ :1g: : : : :  Z:   167.447:  159.565:   51.326:   85.213: 9999.000: 9999.000:  180.000
+ :7d: : : : :  Z:  9999.000: 9999.000: 9999.000:   83.856:  238.750:  256.875:  180.000
+ :7d: : : : :  Z:    69.562:  170.200:   52.800:   85.287: 9999.000: 9999.000:  180.000
+ :3d: : : : :  Z:  9999.000: 9999.000: 9999.000:   85.295:  244.085:  203.815:  180.000
+ :3d: : : : :  Z:    65.880:  181.130:   54.680:   86.035: 9999.000: 9999.000:  180.000
+ :5d: : : : :  Z:  9999.000: 9999.000: 9999.000:   79.671:  202.471:   63.064:  180.000
+ :5d: : : : :  Z:    68.164:  143.450:   49.664:   82.757: 9999.000: 9999.000:  180.000
+ :3g: : : : :  Z:  9999.000: 9999.000: 9999.000:   84.000:  195.000:  146.000:  180.000
+ :3g: : : : :  Z:   170.000:  170.000:   52.000:   84.000: 9999.000: 9999.000:  180.000
+ :1e: : : : :  Z:  9999.000: 9999.000: 9999.000:   80.514:  200.545:  280.510:  180.000
+ :1e: : : : :  Z:   249.314:   82.662:  167.890:   85.507: 9999.000: 9999.000:  180.000
+ :1c: : : : :  Z:  9999.000: 9999.000: 9999.000:   80.223:  196.591:  291.299:  180.000
+ :1c: : : : :  Z:   153.060:  194.379:  179.061:   83.648: 9999.000: 9999.000:  180.000
+ :1f: : : : :  Z:  9999.000: 9999.000: 9999.000:   81.395:  203.030:  294.445:  180.000
+ :1f: : : : :  Z:   172.195:  138.540:  175.565:   84.470: 9999.000: 9999.000:  180.000
+ :5j: : : : :  Z:  9999.000: 9999.000: 9999.000:   87.417:  223.558:   80.175:  180.000
+ :5j: : : : :  Z:    66.667:  109.150:  176.475:   83.833: 9999.000: 9999.000:  180.000
+ :5n: : : : :  Z:  9999.000: 9999.000: 9999.000:   86.055:  246.502:  100.392:  180.000
+ :5n: : : : :  Z:    73.595:  213.752:  183.395:   85.483: 9999.000: 9999.000:  180.000
+ :!!: : : : :  Z:  9999.000: 9999.000: 9999.000:    0.000:    0.000:    0.000:    0.000
+ :!!: : : : :  Z:     0.000:    0.000:    0.000:    0.000: 9999.000: 9999.000:    0.000
+ :1b: : : : :  Z:  9999.000: 9999.000: 9999.000:   84.215:  215.014:  288.672:  180.000
+ :1b: : : : :  Z:   300.420:  177.476:   58.307:  144.841: 9999.000: 9999.000:  180.000
+ :1[: : : : :  Z:  9999.000: 9999.000: 9999.000:   82.731:  220.463:  288.665:  180.000
+ :1[: : : : :  Z:   296.983:  221.654:   54.213:  143.771: 9999.000: 9999.000:  180.000
+ :3b: : : : :  Z:  9999.000: 9999.000: 9999.000:   84.700:  226.400:  168.336:  180.000
+ :3b: : : : :  Z:   292.771:  177.629:   48.629:  147.950: 9999.000: 9999.000:  180.000
+ :1z: : : : :  Z:  9999.000: 9999.000: 9999.000:   83.358:  206.042:  277.567:  180.000
+ :1z: : : : :  Z:   195.700:  161.600:   50.750:  145.258: 9999.000: 9999.000:  180.000
+ :5z: : : : :  Z:  9999.000: 9999.000: 9999.000:   82.614:  206.440:   52.524:  180.000
+ :5z: : : : :  Z:   163.669:  148.421:   50.176:  147.590: 9999.000: 9999.000:  180.000
+ :7p: : : : :  Z:  9999.000: 9999.000: 9999.000:   84.285:  236.600:  220.400:  180.000
+ :7p: : : : :  Z:    68.300:  200.122:   53.693:  145.730: 9999.000: 9999.000:  180.000
+ :5p: : : : :  Z:  9999.000: 9999.000: 9999.000:   84.457:  213.286:   69.086:  180.000
+ :5p: : : : :  Z:    75.500:  156.671:   57.486:  147.686: 9999.000: 9999.000:  180.000
+ :1t: : : : :  Z:  9999.000: 9999.000: 9999.000:   81.200:  199.243:  288.986:  180.000
+ :1t: : : : :  Z:   180.286:  194.743:  178.200:  147.386: 9999.000: 9999.000:  180.000
+ :5q: : : : :  Z:  9999.000: 9999.000: 9999.000:   82.133:  204.933:   69.483:  180.000
+ :5q: : : : :  Z:    63.417:  115.233:  176.283:  145.733: 9999.000: 9999.000:  180.000
+ :1o: : : : :  Z:  9999.000: 9999.000: 9999.000:   83.977:  216.508:  287.192:  180.000
+ :1o: : : : :  Z:   297.254:  225.154:  293.738:  150.677: 9999.000: 9999.000:  180.000
+ :7r: : : : :  Z:  9999.000: 9999.000: 9999.000:   84.606:  232.856:  248.125:  180.000
+ :7r: : : : :  Z:    63.269:  181.975:  295.744:  149.744: 9999.000: 9999.000:  180.000
+ :5r: : : : :  Z:  9999.000: 9999.000: 9999.000:   83.000:  196.900:   65.350:  180.000
+ :5r: : : : :  Z:    60.150:  138.425:  292.550:  154.275: 9999.000: 9999.000:  180.000
+ :2a: : : : :  Z:  9999.000: 9999.000: 9999.000:  145.399:  260.339:  288.756:  180.000
+ :2a: : : : :  Z:   288.444:  192.733:   53.097:   84.067: 9999.000: 9999.000:  180.000
+ :4a: : : : :  Z:  9999.000: 9999.000: 9999.000:  146.275:  259.783:  169.958:  180.000
+ :4a: : : : :  Z:   298.450:  169.583:   50.908:   83.967: 9999.000: 9999.000:  180.000
+ :0a: : : : :  Z:  9999.000: 9999.000: 9999.000:  149.286:  223.159:  139.421:  180.000
+ :0a: : : : :  Z:   284.559:  158.107:   47.900:   84.424: 9999.000: 9999.000:  180.000
+ :#a: : : : :  Z:  9999.000: 9999.000: 9999.000:  148.006:  191.944:  146.231:  180.000
+ :#a: : : : :  Z:   289.288:  150.781:   42.419:   84.956: 9999.000: 9999.000:  180.000
+ :4g: : : : :  Z:  9999.000: 9999.000: 9999.000:  148.028:  256.922:  165.194:  180.000
+ :4g: : : : :  Z:   204.961:  165.194:   49.383:   82.983: 9999.000: 9999.000:  180.000
+ :6g: : : : :  Z:  9999.000: 9999.000: 9999.000:  145.337:  262.869:   79.588:  180.000
+ :6g: : : : :  Z:   203.863:  189.688:   58.000:   84.900: 9999.000: 9999.000:  180.000
+ :8d: : : : :  Z:  9999.000: 9999.000: 9999.000:  148.992:  270.596:  240.892:  180.000
+ :8d: : : : :  Z:    62.225:  176.271:   53.600:   87.262: 9999.000: 9999.000:  180.000
+ :4d: : : : :  Z:  9999.000: 9999.000: 9999.000:  149.822:  249.956:  187.678:  180.000
+ :4d: : : : :  Z:    80.433:  198.133:   61.000:   89.378: 9999.000: 9999.000:  180.000
+ :6d: : : : :  Z:  9999.000: 9999.000: 9999.000:  146.922:  241.222:   88.894:  180.000
+ :6d: : : : :  Z:    59.344:  160.683:   52.333:   83.417: 9999.000: 9999.000:  180.000
+ :2g: : : : :  Z:  9999.000: 9999.000: 9999.000:  141.900:  258.383:  286.517:  180.000
+ :2g: : : : :  Z:   178.267:  165.217:   48.350:   84.783: 9999.000: 9999.000:  180.000
+ :2h: : : : :  Z:  9999.000: 9999.000: 9999.000:  147.782:  260.712:  290.424:  180.000
+ :2h: : : : :  Z:   296.200:  177.282:  175.594:   86.565: 9999.000: 9999.000:  180.000
+ :4n: : : : :  Z:  9999.000: 9999.000: 9999.000:  143.722:  227.256:  203.789:  180.000
+ :4n: : : : :  Z:    73.856:  216.733:  194.444:   80.911: 9999.000: 9999.000:  180.000
+ :0i: : : : :  Z:  9999.000: 9999.000: 9999.000:  148.717:  274.683:  100.283:  180.000
+ :0i: : : : :  Z:    80.600:  248.133:  181.817:   82.600: 9999.000: 9999.000:  180.000
+ :6n: : : : :  Z:  9999.000: 9999.000: 9999.000:  150.311:  268.383:   84.972:  180.000
+ :6n: : : : :  Z:    63.811:  191.483:  176.644:   85.600: 9999.000: 9999.000:  180.000
+ :6j: : : : :  Z:  9999.000: 9999.000: 9999.000:  141.633:  244.100:   66.056:  180.000
+ :6j: : : : :  Z:    71.667:  122.167:  182.200:   83.622: 9999.000: 9999.000:  180.000
+ :0k: : : : :  Z:  9999.000: 9999.000: 9999.000:  149.070:  249.780:  111.520:  180.000
+ :0k: : : : :  Z:   278.370:  207.780:  287.820:   86.650: 9999.000: 9999.000:  180.000
+ :2[: : : : :  Z:  9999.000: 9999.000: 9999.000:  146.383:  259.402:  291.275:  180.000
+ :2[: : : : :  Z:   291.982:  210.048:   54.412:  147.760: 9999.000: 9999.000:  180.000
+ :4b: : : : :  Z:  9999.000: 9999.000: 9999.000:  145.256:  244.622:  162.822:  180.000
+ :4b: : : : :  Z:   294.159:  171.630:   45.900:  145.804: 9999.000: 9999.000:  180.000
+ :0b: : : : :  Z:  9999.000: 9999.000: 9999.000:  147.593:  248.421:  112.086:  180.000
+ :0b: : : : :  Z:   274.943:  164.764:   56.843:  146.264: 9999.000: 9999.000:  180.000
+ :4p: : : : :  Z:  9999.000: 9999.000: 9999.000:  150.077:  260.246:  213.785:  180.000
+ :4p: : : : :  Z:    71.900:  207.638:   56.715:  148.131: 9999.000: 9999.000:  180.000
+ :6p: : : : :  Z:  9999.000: 9999.000: 9999.000:  146.415:  257.831:   89.597:  180.000
+ :6p: : : : :  Z:    67.923:  173.051:   55.513:  147.623: 9999.000: 9999.000:  180.000
+ :2z: : : : :  Z:  9999.000: 9999.000: 9999.000:  142.900:  236.550:  268.800:  180.000
+ :2z: : : : :  Z:   180.783:  185.133:   54.467:  143.350: 9999.000: 9999.000:  180.000
+ :4s: : : : :  Z:  9999.000: 9999.000: 9999.000:  149.863:  247.562:  170.488:  180.000
+ :4s: : : : :  Z:   277.938:   84.425:  176.413:  148.087: 9999.000: 9999.000:  180.000
+ :2u: : : : :  Z:  9999.000: 9999.000: 9999.000:  143.940:  258.200:  298.240:  180.000
+ :2u: : : : :  Z:   279.640:  183.680:  183.080:  145.120: 9999.000: 9999.000:  180.000
+ :2o: : : : :  Z:  9999.000: 9999.000: 9999.000:  147.342:  256.475:  295.508:  180.000
+ :2o: : : : :  Z:   287.408:  194.525:  293.725:  150.458: 9999.000: 9999.000:  180.000
+ :epsilon: : : : :  Z:   294.967:  173.990:   53.550:   81.495:  154.000:  288.831:  180.000
+ :epsilon: : : : :  Z:   294.967:  173.990:   53.550:   81.495:  212.250:  288.831:  180.000
+ :alpha: : : : :  Z:   294.967:  173.990:   53.550:   81.495:  212.250:  288.831:  180.000
+ :alpha: : : : :  Z:    24.000:  173.990:   53.550:   81.495:  212.250:  288.831:  180.000
+ :beta: : : : :  Z:   294.967:  173.990:   53.550:   81.495:  212.250:  288.831:  180.000
+ :beta: : : : :  Z:   294.967:   49.000:   53.550:   81.495:  212.250:  288.831:  180.000
+ :zeta: : : : :  Z:   294.967:  173.990:   53.550:   81.495:  212.250:   24.000:  180.000
+ :zeta: : : : :  Z:   294.967:  173.990:   53.550:   81.495:  212.250:  288.831:  180.000
+ :delta-1: : : : :  Z:   294.967:  173.990:   53.550:   59.000:  212.250:  288.831:  180.000
+ :delta-1: : : : :  Z:   294.967:  173.990:   53.550:   81.495:  212.250:  288.831:  180.000
+ :gamma: : : : :  Z:   294.967:  173.990:   53.550:   81.495:  212.250:  288.831:  180.000
+ :gamma: : : : :  Z:   294.967:  173.990:  139.000:   81.495:  212.250:  288.831:  180.000
+ :delta: : : : :  Z:   294.967:  173.990:   53.550:   81.495:  212.250:  288.831:  180.000
+ :delta: : : : :  Z:   294.967:  173.990:   53.550:   59.000:  212.250:  288.831:  180.000
+'''
+
+# A manually selected group of test cases from real files
+# designed to test each code path through membership(). The first residue from
+# each pair has intentionally been damaged so that it will not produce separate
+# output from a report. We use the causes option to illustrate what code path
+# is being used.
+input2 = u'''2xLk:1: C:  11: : :  G:__?__:__?__:__?__:81.132:-127.583:-70.677
+2xLk:1: C:  12: : :  U:169.008:153.891:51.391:80.277:-135.347:-70.614
+3cgp:1: B:  19: : :  U:__?__:__?__:__?__:82.839:-147.528:-179.087
+3cgp:1: B:  20: : :  A:139.983:-154.445:63.134:88.055:-145.599:70.874
+4pco:1: B:   3: : :  U:__?__:__?__:__?__:77.659:-165.227:-68.525
+4pco:1: B:   4: : :  G:151.914:-179.903:176.058:83.039:-148.171:-66.728
+5b2q:1: B:  62: : :  G:__?__:__?__:__?__:83.537:-131.816:-116.417
+5b2q:1: B:  63: : :  U:-69.320:-146.615:47.107:147.038:-148.815:45.665
+6gc5:1: F:   2: : :  U:__?__:__?__:__?__:144.610:-116.227:152.694
+6gc5:1: F:   3: : :  U:-66.167:162.580:41.697:145.644:-122.673:127.881
+3bns:1: A:  21: : :  C:__?__:__?__:__?__:76.224:-166.174:-73.594
+3bns:1: A:  22: : :  G:150.784:-158.788:175.706:87.605:-146.172:-63.516
+3gm7:1: H:   5: : :  U:__?__:__?__:__?__:68.910:-153.989:-56.381
+3gm7:1: H:   6: : :  G:-105.747:164.057:92.120:74.597:-150.523:-79.724
+6qit:1: A:   2: : :  C:__?__:__?__:__?__:82.169:-138.695:-63.417
+6qit:1: A:   3: : :  A:-71.504:-131.618:54.061:144.409:-95.827:-140.754
+3rer:1: K:   7: : :  U:__?__:__?__:__?__:87.510:-99.276:-118.108
+3rer:1: K:   8: : :  A:-66.924:-158.118:48.287:81.250:__?__:__?__
+3diL:1: A:  59: : :  C:__?__:__?__:__?__:80.668:-145.667:-36.026
+3diL:1: A:  60: : :  G:-143.441:115.188:149.951:86.379:-141.567:-69.901
+5ho4:1: B:   3: : :  G:__?__:__?__:__?__:160.213:-123.685:-174.677
+5ho4:1: B:   4: : :  G:-107.676:163.883:39.081:85.911:-157.392:-71.638
+4mcf:1: E:   4: : :  U:__?__:__?__:__?__:78.239:-156.881:-70.399
+4mcf:1: E:   5: : :  G:-91.794:163.594:87.552:70.675:-141.886:-72.556
+3pdr:1: A:  59: : :  C:__?__:__?__:__?__:80.441:-149.674:-76.690
+3pdr:1: A:  60: : :  A:-62.415:171.383:47.537:79.461:-145.680:-71.359
+3gm7:1: G:   1: : :  C:__?__:__?__:__?__:84.065:-128.784:-61.905
+3gm7:1: G:   2: : :  U:-76.914:-166.398:55.279:74.218:-157.766:-64.720
+6h0r:1: B:  15: : :  U:__?__:__?__:__?__:83.971:-122.349:-103.636
+6h0r:1: B:  16: : :  U:-30.804:145.657:33.314:81.109:-141.719:-75.527
+2zko:1: C:  13: : :  G:__?__:__?__:__?__:76.629:-150.027:-67.298
+2zko:1: C:  14: : :  C:-70.016:164.567:71.735:76.499:-160.106:-73.474
+3pdr:1: X: 138: : :  U:__?__:__?__:__?__:77.324:-177.192:-105.412
+3pdr:1: X: 139: : :  A:-46.950:179.570:49.599:71.442:-143.233:-61.461
+4jah:1: B:  10: : :  U:__?__:__?__:__?__:85.890:-164.804:-95.055
+4jah:1: B:  11: : :  G:-64.134:178.767:49.773:77.067:-152.496:-70.128
+3diL:1: A:  13: : :  C:__?__:__?__:__?__:135.303:-125.074:-69.725
+3diL:1: A:  14: : :  G:75.452:147.741:32.719:83.048:-146.012:-75.223
+3pdr:1: X: 132: : :  U:__?__:__?__:__?__:77.469:-157.795:-115.458
+3pdr:1: X: 133: : :  U:47.309:136.943:-25.259:83.460:-150.210:-61.763
+'''
+
+output1='''
+ :1a: : : : :  Z 33 p 1a 1.000
+ :1m: : : : :  Z 33 p 1m 1.000
+ :1L: : : : :  Z 33 p 1L 1.000
+ :&a: : : : :  Z 33 p &a 1.000
+ :7a: : : : :  Z 33 p 7a 1.000
+ :3a: : : : :  Z 33 p 3a 1.000
+ :9a: : : : :  Z 33 p 9a 1.000
+ :1g: : : : :  Z 33 p 1g 1.000
+ :7d: : : : :  Z 33 p 7d 1.000
+ :3d: : : : :  Z 33 p 3d 1.000
+ :5d: : : : :  Z 33 p 5d 1.000
+ :3g: : : : :  Z 33 p 3g 1.000 wannabe
+ :1e: : : : :  Z 33 t 1e 1.000
+ :1c: : : : :  Z 33 t 1c 1.000
+ :1f: : : : :  Z 33 t 1f 1.000
+ :5j: : : : :  Z 33 t 5j 1.000
+ :5n: : : : :  Z 33 t 5n 1.000 wannabe
+ :!!: : : : :  Z trig !! 0.000 epsilon-1
+ :1b: : : : :  Z 32 p 1b 1.000
+ :1[: : : : :  Z 32 p 1[ 1.000
+ :3b: : : : :  Z 32 p 3b 1.000
+ :1z: : : : :  Z 32 p 1z 1.000
+ :5z: : : : :  Z 32 p 5z 1.000
+ :7p: : : : :  Z 32 p 7p 1.000
+ :5p: : : : :  Z 32 p 5p 1.000 wannabe
+ :1t: : : : :  Z 32 t 1t 1.000
+ :5q: : : : :  Z 32 t 5q 1.000
+ :1o: : : : :  Z 32 m 1o 1.000
+ :7r: : : : :  Z 32 m 7r 1.000
+ :5r: : : : :  Z 32 m 5r 1.000 wannabe
+ :2a: : : : :  Z 23 p 2a 1.000
+ :4a: : : : :  Z 23 p 4a 1.000
+ :0a: : : : :  Z 23 p 0a 1.000
+ :#a: : : : :  Z 23 p #a 1.000
+ :4g: : : : :  Z 23 p 4g 1.000
+ :6g: : : : :  Z 23 p 6g 1.000
+ :8d: : : : :  Z 23 p 8d 1.000
+ :4d: : : : :  Z 23 p 4d 1.000
+ :6d: : : : :  Z 23 p 6d 1.000
+ :2g: : : : :  Z 23 p 2g 1.000 wannabe
+ :2h: : : : :  Z 23 t 2h 1.000
+ :4n: : : : :  Z 23 t 4n 1.000
+ :0i: : : : :  Z 23 t 0i 1.000
+ :6n: : : : :  Z 23 t 6n 1.000
+ :6j: : : : :  Z 23 t 6j 1.000
+ :0k: : : : :  Z 23 m 0k 1.000 wannabe
+ :2[: : : : :  Z 22 p 2[ 1.000
+ :4b: : : : :  Z 22 p 4b 1.000
+ :0b: : : : :  Z 22 p 0b 1.000
+ :4p: : : : :  Z 22 p 4p 1.000
+ :6p: : : : :  Z 22 p 6p 1.000
+ :2z: : : : :  Z 22 p 2z 1.000 wannabe
+ :4s: : : : :  Z 22 t 4s 1.000
+ :2u: : : : :  Z 22 t 2u 1.000 wannabe
+ :2o: : : : :  Z 22 m 2o 1.000
+ :epsilon: : : : :  Z trig !! 0.000 epsilon-1
+ :alpha: : : : :  Z 33 p 1a 0.999
+ :alpha: : : : :  Z trig !! 0.000 alpha
+ :beta: : : : :  Z 33 p 1a 0.999
+ :beta: : : : :  Z trig !! 0.000 beta
+ :zeta: : : : :  Z 33 p 1a 0.999
+ :zeta: : : : :  Z trig !! 0.000 zeta-1
+ :delta-1: : : : :  Z trig !! 0.000 delta
+ :delta-1: : : : :  Z trig !! 0.000 delta-1
+ :gamma: : : : :  Z 33 p 1a 0.999
+ :gamma: : : : :  Z trig !! 0.000 gamma
+ :delta: : : : :  Z 33 p 1a 0.999
+ :delta: : : : :  Z trig !! 0.000 delta
+'''
+
+output2 = '''2xLk:1: C:  12: : :  U 33 p 1g 0.839 1-only-one
+3cgp:1: B:  20: : :  A 33 p 3g 0.040 1-only-one wannabe
+4pco:1: B:   4: : :  G 33 t 1c 0.890 2-BETWEEN-dom-sat(   0.22|  0.913)
+5b2q:1: B:  63: : :  U 32 p 1[ 0.072 2-BETWEEN-dom-sat(  0.941|  0.829)
+6gc5:1: F:   3: : :  U 22 p 4b 0.889 2-None-dom
+3bns:1: A:  22: : :  G 33 t 1c 0.901 2-OUTSIDE-dom
+3gm7:1: H:   6: : :  G 33 p !! 0.000 7D dist 1a
+6qit:1: A:   3: : :  A 32 p 1[ 0.899 2-OUTSIDE-sat
+3rer:1: K:   8: : :  A 33 p 7a 0.047 2-None-dom
+3diL:1: A:  60: : :  G 33 t !! 0.000 7D dist 1e
+5ho4:1: B:   4: : :  G 23 p !! 0.000 7D dist 4a
+4mcf:1: E:   5: : :  G 33 p !! 0.000 7D dist 1a
+3pdr:1: A:  60: : :  A 33 p 1a 0.916 4-BETWEEN-dom-sat(    0.1|    1.2)
+3gm7:1: G:   2: : :  U 33 p 1a 0.589 4-BETWEEN-dom-sat(  0.428|  0.904)
+6h0r:1: B:  16: : :  U 33 p 1L 0.033 4-BETWEEN-dom-sat(  0.862|  0.655)
+2zko:1: C:  14: : :  C 33 p 1a 0.444 4-OUTSIDE-dom
+3pdr:1: X: 139: : :  A 33 p &a 0.555 4-OUTSIDE-sat
+4jah:1: B:  11: : :  G 33 p &a 0.912 5-BETWEEN-dom-sat(  0.442|  0.226)
+3diL:1: A:  14: : :  G 23 p !! 0.000 outlier distance 1.01
+3pdr:1: X: 133: : :  U 33 m !! 0.000 vacant bin
+'''
+
+def test(input, canonicalOutput, options, identity):
+  opt = suites.parseOptions(options)
+  stream = StringIO(input)
+  outFile=StringIO()
+  suitename.clearStats()
+  suitename.main(stream, outFile=outFile, optionsIn=opt)
+  suitename.clearStats()
+
+  output = outFile.getvalue()
+  assert output.strip() == canonicalOutput.strip(), identity
+
+
+def testAll():
+  test(input1, output1, "chart=true noinc=true", "cluster and triage test")
+  test(input2, output2, "chart=true noinc=true causes=true",
+      "code paths test")
+  test(regression.in_1ehz, regression.out_1ehz, "", "1ehz regression test")
+
+
+# Not normally used, but useful for diagnosing failures
+def testVerbose(input, canonicalOutput, options, identity):
+  opt = suites.parseOptions(options)
+  stream = StringIO(input)
+  outFile=StringIO("")
+  suitename.clearStats()
+  suitename.main(stream, outFile=outFile, optionsIn=opt)
+
+  output = outFile.getvalue()
+  if output.strip() == canonicalOutput.strip():
+   result = True
+   sys.stderr.write("Success\n")
+  else:
+    result = False
+    sys.stderr.write("Failed\n")
+    sys.stderr.write("========================================\n")
+    sys.stderr.write(canonicalOutput.strip())
+    sys.stderr.write("\n\n=========================================\n")
+    sys.stderr.write(output.strip())
+
+    out2 = open("UnitTest-output.txt", "w")
+    out2.write(output.strip())
+    out2.close()
+    return result
+
+testAll()
diff --git a/mmtbx/suitename/unit-test/gentest.py b/mmtbx/suitename/unit-test/gentest.py
new file mode 100644
index 0000000000..4172456589
--- /dev/null
+++ b/mmtbx/suitename/unit-test/gentest.py
@@ -0,0 +1,104 @@
+# This is an infrequently run file used to generate
+# unit test data from the data inherent in suitename.py
+# Its output was used to create the canonicalOutput string in UnitTest.py
+
+from __future__ import division
+import os, sys, inspect, copy
+
+currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))
+parentdir = os.path.dirname(currentdir)
+sys.path.insert(0, parentdir)
+
+from suitename import triageFilters, sieveDelta, sieveGamma
+from suiteninit import bins
+from suitenamedefs import Residue
+from suiteninput import readResidues
+
+import numpy as np
+from io import StringIO
+
+
+def run():
+    global pair
+    "generate a residue set that will fail each triage in turn"
+    s1a = " :1a: : : : Z:   294.967:  173.990:   53.550:   81.495:  212.250:  288.831:  180.000"
+    standard = readResidues(StringIO(s1a))
+    standard.append(copy.deepcopy(standard[0]))
+    # standard is now two copies of the "1a" cluster center
+    pair = 2 * standard
+    for key, filter in triageFilters.items():
+        testFromFilter(key, filter)
+    testFromSieve("delta-1", 1, sieveDelta)
+    testFromSieve("gamma", 6, sieveGamma)
+    testFromSieve("delta", 7, sieveDelta)
+
+
+def testFromFilter(id, filter):
+    index = filter[0]
+    min = filter[1]
+    testCore(id, index, min)
+
+
+def testFromSieve(id, index, sieve):
+    min = sieve[0][0]
+    testCore(id, index, min)
+
+
+def testCore(id, index, min):
+    basic = copy.deepcopy(pair)
+    if index >= 4:
+        m = index - 4
+        n = 1
+    else:
+        m = index + 2
+        n = 0
+    basic[n].angle[m] = min - 1.0
+    r1 = Residue(" :{}: : : : Z:".format(id), "Z", basic[0].angle)
+    r2 = Residue(" :{}: : : : Z:",format(id), "Z", basic[1].angle)
+    printResidue(r1)
+    printResidue(r2)
+
+
+def run1():
+    "Generate a residue set from the cluster definitions"
+    for i in range(1, 13):
+        b = bins[i]
+        for c in b.cluster[1:]:
+            residuesFromCluster(c.angle, c.name)
+
+
+def residuesFromCluster(angles, id):
+    angles1 = np.array(3 * [9999] + list(angles[1:4]) + list(angles[0:1]))
+    angles2 = np.array(list(angles[4:8]) + 2 * [9999] + list(angles[8:9]))
+    r1 = Residue(" :{}: : : : Z:".format(id), "Z", angles1)
+    r2 = Residue(" :{}: : : : Z:".format(id), "Z", angles2)
+    printResidue(r1)
+    printResidue(r2)
+
+
+def printResidue(r):
+    angleText = ":".join(["{:9.3f}".format(x) for x in r.angle])
+    print(r.pointIDs, angleText)
+
+
+temp = """chiMinus
+    # 1   deltaMinus
+    # 2   epsilon
+    # 3   zeta
+    # 4   alpha
+    # 5   beta
+    # 6   gamma
+    # 7   delta
+    # 8   chi    """
+
+
+def rearrange():
+    words = []
+    for line in temp.splitlines():
+        word = line.split()[-1]
+        words.append(word)
+    text = ", ".join(words)
+    return text
+
+
+run()
diff --git a/mmtbx/suitename/unit-test/regression.py b/mmtbx/suitename/unit-test/regression.py
new file mode 100644
index 0000000000..c10dd1a0cb
--- /dev/null
+++ b/mmtbx/suitename/unit-test/regression.py
@@ -0,0 +1,218 @@
+from __future__ import nested_scopes, generators, division, absolute_import
+from __future__ import  with_statement, print_function, unicode_literals
+from io import StringIO
+
+from mmtbx.suitename.suitename import main
+from mmtbx.suitename import suites
+
+in_1ehz = \
+""" :1: A:   1: : :  G:__?__:-128.053:67.793:82.863:-155.568:-68.623
+ :1: A:   2: : :  C:-67.449:-178.393:53.832:83.378:-145.149:-76.788
+ :1: A:   3: : :  G:-74.460:169.695:59.477:80.733:-148.283:-80.025
+ :1: A:   4: : :  G:-64.437:162.194:60.724:82.164:-157.372:-68.659
+ :1: A:   5: : :  A:-74.745:-176.472:53.356:84.944:-137.465:-81.718
+ :1: A:   6: : :  U:-48.819:157.567:55.308:81.262:-150.997:-77.008
+ :1: A:   7: : :  U:-59.498:-178.720:62.506:137.257:-105.934:-51.960
+ :1: A:   8: : :  U:-83.821:-145.649:55.418:78.644:-142.800:-118.577
+ :1: A:   9: : :  A:-69.711:-141.745:52.257:147.772:-106.213:-77.257
+ :1: A:  10: : :2MG:177.814:147.203:__?__:__?__:-126.196:-88.738
+ :1: A:  11: : :  C:-56.137:167.906:48.160:87.215:-150.469:-69.893
+ :1: A:  12: : :  U:-67.796:172.869:51.821:80.665:-158.532:-65.168
+ :1: A:  13: : :  C:166.574:-169.891:178.570:82.512:-153.131:-97.390
+ :1: A:  14: : :  A:83.353:-158.279:-114.591:92.038:-125.519:-57.251
+ :1: A:  15: : :  G:-55.064:162.460:51.892:79.826:-136.303:-143.947
+ :1: A:  16: : :H2U:-6.131:91.184:__?__:__?__:-61.796:-131.179
+ :1: A:  17: : :H2U:27.788:107.741:__?__:__?__:177.959:76.191
+ :1: A:  18: : :  G:45.383:-159.437:59.039:150.600:-95.226:-179.094
+ :1: A:  19: : :  G:-71.420:-178.863:53.824:153.803:-91.620:-83.735
+ :1: A:  20: : :  G:-81.332:-150.736:47.756:89.936:-122.273:-54.104
+ :1: A:  21: : :  A:-75.551:148.565:-176.602:78.228:-168.936:-75.552
+ :1: A:  22: : :  G:158.774:153.494:179.253:81.999:-144.965:-80.445
+ :1: A:  23: : :  A:-53.308:174.758:52.482:82.300:-155.347:-66.440
+ :1: A:  24: : :  G:-68.767:178.230:46.822:83.604:-144.254:-72.825
+ :1: A:  25: : :  C:-65.087:168.872:53.883:83.255:-145.086:-68.361
+ :1: A:  26: : :M2G:-53.808:170.843:__?__:__?__:-136.277:-76.875
+ :1: A:  27: : :  C:-52.975:166.853:43.623:83.427:-148.537:-73.433
+ :1: A:  28: : :  C:-72.439:178.326:49.309:80.115:-152.098:-66.986
+ :1: A:  29: : :  A:-66.553:173.956:55.560:81.361:-155.521:-78.260
+ :1: A:  30: : :  G:-53.966:165.929:56.930:83.563:-144.716:-62.268
+ :1: A:  31: : :  A:-69.947:177.767:52.310:83.748:-136.979:-75.532
+ :1: A:  32: : :OMC:-52.715:__?__:__?__:__?__:-145.941:-71.172
+ :1: A:  33: : :  U:-67.688:-177.012:47.001:82.123:-148.003:-53.686
+ :1: A:  34: : :OMG:171.113:148.145:__?__:__?__:-132.519:-71.830
+ :1: A:  35: : :  A:-47.663:163.713:40.211:80.882:-143.671:-59.500
+ :1: A:  36: : :  A:-52.433:165.690:51.257:72.220:-160.374:-85.207
+ :1: A:  37: : :YYG:-57.536:162.974:__?__:__?__:-148.148:-67.046
+ :1: A:  38: : :  A:-61.808:-179.981:46.930:82.458:-136.769:-76.362
+ :1: A:  39: : :PSU:-47.734:160.379:53.321:79.261:-140.089:-68.553
+ :1: A:  40: : :5MC:-67.366:172.001:__?__:__?__:-154.223:-74.895
+ :1: A:  41: : :  U:-68.226:-179.367:52.444:78.871:-137.302:-84.708
+ :1: A:  42: : :  G:-47.869:158.673:55.597:79.758:-160.327:-70.305
+ :1: A:  43: : :  G:-67.015:-178.303:55.596:81.592:-154.873:-76.392
+ :1: A:  44: : :  A:-59.692:162.125:60.042:85.307:-142.834:-57.171
+ :1: A:  45: : :  G:-71.930:-176.932:50.985:87.553:-135.100:-78.652
+ :1: A:  46: : :7MG:-56.771:-146.515:__?__:__?__:-102.653:-137.859
+ :1: A:  47: : :  U:62.398:-164.007:44.365:146.142:-93.662:-78.042
+ :1: A:  48: : :  C:-73.483:-174.257:161.496:145.555:-143.540:75.586
+ :1: A:  49: : :5MC:50.694:168.495:__?__:__?__:-145.028:-82.064
+ :1: A:  50: : :  U:-51.651:177.223:42.076:80.392:-150.600:-67.789
+ :1: A:  51: : :  G:-63.906:176.814:52.813:79.429:-150.411:-71.324
+ :1: A:  52: : :  U:-64.688:173.639:48.543:80.289:-156.462:-69.422
+ :1: A:  53: : :  G:-56.867:171.464:56.154:83.914:-159.389:-64.873
+ :1: A:  54: : :5MU:-79.665:__?__:__?__:__?__:-128.577:-70.746
+ :1: A:  55: : :PSU:-49.745:168.796:44.052:76.601:-140.782:-69.884
+ :1: A:  56: : :  C:166.416:171.805:53.305:83.405:-132.708:-70.648
+ :1: A:  57: : :  G:-65.658:167.070:57.463:81.657:-145.232:-67.602
+ :1: A:  58: : :1MA:-60.812:-146.102:__?__:__?__:-78.282:-169.333
+ :1: A:  59: : :  U:72.618:-158.840:63.678:84.575:-148.820:-53.724
+ :1: A:  60: : :  C:-72.162:179.454:66.018:148.304:-97.069:-66.356
+ :1: A:  61: : :  C:-84.288:179.827:38.237:82.964:-152.342:-74.495
+ :1: A:  62: : :  A:-60.070:179.638:46.853:80.531:-145.618:-74.140
+ :1: A:  63: : :  C:-62.017:167.262:50.948:80.715:-152.333:-70.724
+ :1: A:  64: : :  A:-66.852:179.977:44.052:75.809:-147.481:-76.520
+ :1: A:  65: : :  G:-43.958:164.226:49.876:79.755:-152.032:-73.335
+ :1: A:  66: : :  A:-57.891:178.545:52.046:81.724:-150.972:-73.549
+ :1: A:  67: : :  A:-61.953:164.096:54.240:83.200:-152.186:-78.339
+ :1: A:  68: : :  U:-59.846:175.336:47.283:82.185:-152.942:-65.390
+ :1: A:  69: : :  U:-63.754:168.082:55.094:79.119:-155.366:-85.643
+ :1: A:  70: : :  C:-61.653:164.642:53.122:78.973:-158.517:-64.505
+ :1: A:  71: : :  G:-78.355:173.623:60.313:80.276:-149.552:-68.409
+ :1: A:  72: : :  C:-73.183:176.236:62.079:82.989:-152.335:-67.878
+ :1: A:  73: : :  A:-63.297:177.662:50.401:81.559:-148.155:-66.150
+ :1: A:  74: : :  C:-66.885:-174.928:50.690:85.875:-145.004:-58.835
+ :1: A:  75: : :  C:-52.334:175.666:42.302:85.557:-131.936:163.904
+ :1: A:  76: : :  A:-70.955:130.178:164.639:160.907:__?__:__?__
+"""
+
+out_1ehz = (\
+  " :1: A:   1: : :  G inc  __ 0.000\n"+
+  " :1: A:   2: : :  C 33 p 1a 0.935\n"+
+  " :1: A:   3: : :  G 33 p 1a 0.868\n"+
+  " :1: A:   4: : :  G 33 p 1a 0.842\n"+
+  " :1: A:   5: : :  A 33 p 1a 0.847\n"+
+  " :1: A:   6: : :  U 33 p 1a 0.664\n"+
+  " :1: A:   7: : :  U 32 p 1b 0.803\n"+
+  " :1: A:   8: : :  U 23 p 2a 0.509\n"+
+  " :1: A:   9: : :  A 32 p 1[ 0.046\n"+
+  " :1: A:  10: : :2MG inc  __ 0.000\n"+
+  " :1: A:  11: : :  C inc  __ 0.000\n"+
+  " :1: A:  12: : :  U 33 p 1a 0.898\n"+
+  " :1: A:  13: : :  C 33 t 1c 0.824\n"+
+  " :1: A:  14: : :  A trig !! 0.000 gamma\n"+
+  " :1: A:  15: : :  G 33 p 1a 0.484\n"+
+  " :1: A:  16: : :H2U inc  __ 0.000\n"+
+  " :1: A:  17: : :H2U inc  __ 0.000\n"+
+  " :1: A:  18: : :  G inc  __ 0.000\n"+
+  " :1: A:  19: : :  G 22 p 4b 0.512\n"+
+  " :1: A:  20: : :  G 23 p 2a 0.623\n"+
+  " :1: A:  21: : :  A 33 t !! 0.000\n"+
+  " :1: A:  22: : :  G 33 t 1f 0.714\n"+
+  " :1: A:  23: : :  A 33 p 1a 0.840\n"+
+  " :1: A:  24: : :  G 33 p 1a 0.881\n"+
+  " :1: A:  25: : :  C 33 p 1a 0.967\n"+
+  " :1: A:  26: : :M2G inc  __ 0.000\n"+
+  " :1: A:  27: : :  C inc  __ 0.000\n"+
+  " :1: A:  28: : :  C 33 p 1a 0.923\n"+
+  " :1: A:  29: : :  A 33 p 1a 0.973\n"+
+  " :1: A:  30: : :  G 33 p 1a 0.838\n"+
+  " :1: A:  31: : :  A 33 p 1a 0.914\n"+
+  " :1: A:  32: : :OMC inc  __ 0.000\n"+
+  " :1: A:  33: : :  U inc  __ 0.000\n"+
+  " :1: A:  34: : :OMG inc  __ 0.000\n"+
+  " :1: A:  35: : :  A inc  __ 0.000\n"+
+  " :1: A:  36: : :  A 33 p 1a 0.670\n"+
+  " :1: A:  37: : :YYG inc  __ 0.000\n"+
+  " :1: A:  38: : :  A inc  __ 0.000\n"+
+  " :1: A:  39: : :PSU 33 p 1a 0.680\n"+
+  " :1: A:  40: : :5MC inc  __ 0.000\n"+
+  " :1: A:  41: : :  U inc  __ 0.000\n"+
+  " :1: A:  42: : :  G 33 p 1a 0.630\n"+
+  " :1: A:  43: : :  G 33 p 1a 0.882\n"+
+  " :1: A:  44: : :  A 33 p 1a 0.837\n"+
+  " :1: A:  45: : :  G 33 p 1a 0.749\n"+
+  " :1: A:  46: : :7MG inc  __ 0.000\n"+
+  " :1: A:  47: : :  U inc  __ 0.000\n"+
+  " :1: A:  48: : :  C 22 t 2u 0.283 wannabe\n"+
+  " :1: A:  49: : :5MC inc  __ 0.000\n"+
+  " :1: A:  50: : :  U inc  __ 0.000\n"+
+  " :1: A:  51: : :  G 33 p 1a 0.981\n"+
+  " :1: A:  52: : :  U 33 p 1a 0.945\n"+
+  " :1: A:  53: : :  G 33 p 1a 0.896\n"+
+  " :1: A:  54: : :5MU inc  __ 0.000\n"+
+  " :1: A:  55: : :PSU inc  __ 0.000\n"+
+  " :1: A:  56: : :  C 33 p 1g 0.894\n"+
+  " :1: A:  57: : :  G 33 p 1a 0.837\n"+
+  " :1: A:  58: : :1MA inc  __ 0.000\n"+
+  " :1: A:  59: : :  U inc  __ 0.000\n"+
+  " :1: A:  60: : :  C 32 p 1b 0.662\n"+
+  " :1: A:  61: : :  C 23 p 2a 0.553\n"+
+  " :1: A:  62: : :  A 33 p 1a 0.895\n"+
+  " :1: A:  63: : :  C 33 p 1a 0.964\n"+
+  " :1: A:  64: : :  A 33 p 1a 0.791\n"+
+  " :1: A:  65: : :  G 33 p 1a 0.586\n"+
+  " :1: A:  66: : :  A 33 p 1a 0.940\n"+
+  " :1: A:  67: : :  A 33 p 1a 0.941\n"+
+  " :1: A:  68: : :  U 33 p 1a 0.891\n"+
+  " :1: A:  69: : :  U 33 p 1a 0.951\n"+
+  " :1: A:  70: : :  C 33 p 1a 0.809\n"+
+  " :1: A:  71: : :  G 33 p 1a 0.761\n"+
+  " :1: A:  72: : :  C 33 p 1a 0.832\n"+
+  " :1: A:  73: : :  A 33 p 1a 0.965\n"+
+  " :1: A:  74: : :  C 33 p 1a 0.886\n"+
+  " :1: A:  75: : :  C 33 p 1a 0.639\n"+
+  " :1: A:  76: : :  A 32 t !! 0.000\n"+
+  " all general case widths, power = 3.00\n"+
+  "Found 52 complete suites derived from 76 entries\n"+
+  "1 suites were triaged, leaving 51 assigned to bins\n"+
+  "For all 51 suites: average suiteness== 0.750 (power==3.00)\n"+
+  "     2 suites are  outliers\n"+
+  "     0 suites have suiteness == 0    \n"+
+  "     1 suites have suiteness >  0 <.1\n"+
+  "     0 suites have suiteness >=.1 <.2\n"+
+  "     1 suites have suiteness >=.2 <.3\n"+
+  "     0 suites have suiteness >=.3 <.4\n"+
+  "     1 suites have suiteness >=.4 <.5\n"+
+  "     4 suites have suiteness >=.5 <.6\n"+
+  "     7 suites have suiteness >=.6 <.7\n"+
+  "     4 suites have suiteness >=.7 <.8\n"+
+  "    19 suites have suiteness >=.8 <.9\n"+
+  "    12 suites have suiteness >=.9    \n"+
+  " A form (1a) 38 suites: average suiteness== 0.838 (power==3.00)\n"+
+  "     0 suites have suiteness == 0    \n"+
+  "     0 suites have suiteness >  0 <.1\n"+
+  "     0 suites have suiteness >=.1 <.2\n"+
+  "     0 suites have suiteness >=.2 <.3\n"+
+  "     0 suites have suiteness >=.3 <.4\n"+
+  "     1 suites have suiteness >=.4 <.5\n"+
+  "     1 suites have suiteness >=.5 <.6\n"+
+  "     5 suites have suiteness >=.6 <.7\n"+
+  "     3 suites have suiteness >=.7 <.8\n"+
+  "    16 suites have suiteness >=.8 <.9\n"+
+  "    12 suites have suiteness >=.9    \n"+
+  " non-1a  has 11 suites: average suiteness== 0.584 (power==3.00)\n"+
+  "     0 suites have suiteness == 0    \n"+
+  "     1 suites have suiteness >  0 <.1\n"+
+  "     0 suites have suiteness >=.1 <.2\n"+
+  "     1 suites have suiteness >=.2 <.3\n"+
+  "     0 suites have suiteness >=.3 <.4\n"+
+  "     0 suites have suiteness >=.4 <.5\n"+
+  "     3 suites have suiteness >=.5 <.6\n"+
+  "     2 suites have suiteness >=.6 <.7\n"+
+  "     1 suites have suiteness >=.7 <.8\n"+
+  "     3 suites have suiteness >=.8 <.9\n"+
+  "     0 suites have suiteness >=.9    \n"
+)
+
+def regression(name, input, canonical_output):
+  inp = StringIO(input)
+  out = StringIO()
+  options2 = suites.parseOptions("report=True")
+
+  main(inStream=inp, outFile=out, optionsIn=options2)
+  output = out.getvalue()
+  assert output == canonical_output, name + " regression test failed"
+
+
+def test():
+  regression("1ehz", in_1ehz, out_1ehz)
+