Repository: cctbx/cctbx_project (language: python)

Issue:
Matplotlib colormap deprecation
According to [matplotlib issue #20853](https://github.com/matplotlib/matplotlib/issues/20853) the interface for colormaps in matplotlib (`mpl`) changes. Method `mpl.cm.get_cmap('name')` is deprecated in version 3.7 and will be removed in version 3.9, probably in a year or so. Instead, `mpl.colormap.get_cmap('name')` or `mpl.colormap['name']` is encouraged. `plt.get_cmap()` is expected to remain available.

This is but a very minor issue for `cctbx.xfel` and other areas in active development. As they continue to use the most recent matplotlib version, the problem can be simply circumvented by refactoring all code to use the new syntax.

Nonetheless, this made me think about legacy software and codebase supported by cctbx. In their case, refactoring the code might break any functionality plotting with gradients, as the new syntax was introduced in matplotlib version 3.4.0, barely two years ago (Mar 26, 2021). If the current codebase should indeed support older applications such as the ones using python2, it needs to allow access to the current syntax existing back from matplotlib [version 2.2](https://matplotlib.org/2.2.4/api/cm_api.html?highlight=get_cmap#matplotlib.cm.get_cmap) from June 28, 2018.


Discussion / hints:
Yep -- I think we have a fair number of try-excepts built into cctbx for cases like what you describe, and that's what I'd recommend here as well. It's a little clunky but does seem to do the job.
I agree, although quite frankly I'd just prefer the Python2 support to be dropped. @irisdyoung, are you aware of any place in the code, for example in libtbx, where such aliases are stored centrally? Just thinking about scattering try-excepts in every place in cctbx where colormaps are used makes me feel physically ill.
I guess the closest to what you're hoping for would be how many files import from `six.moves` (as a result of the py2 -> py3 move and hypothetically maintaining backward compatibility). We don't have a centralized location of our own for this, to my knowledge. I think it's a good idea though, with the caveat that I wouldn't want to end up with a ton of unrelated code in one place. What about defining some plotting functionality in a central location, including colormaps? Some plotting utilities could be really handy anyway.
I feel like libtbx already acts as such centralised location; for example see the code for `libtbx.mpi4py` whose behavior changes with environment. I had an idea to define colormaps there as `from libtbx.matplotlib import colormaps` where colormap is conditionally `cm` or `colormaps`. This is still a rough patch, albeit maybe a little better one. If `libtbx` in not an appropriate location for this, then maybe we could have `libtbx.moves`? Just an idea.
Great, I endorse this. Will defer to someone closer to `libtbx` to formally approve adding this.
Hi all, cctbx.xfel depends on packages which have no python 2 support so I believe the argument about maintaining python2 is moot.  At this point we are mostly keeping python3 syntax out of cctbx_project/xfel for convenience in passing our code checkers, and indeed most of the exceptions in [py2_syntax_exceptions.txt](https://github.com/cctbx/cctbx_project/blob/master/.azure-pipelines/py2_syntax_exceptions.txt) are from cctbx.xfel.

Outside of cctbx_project/xfel, I see only `plt.get_cmap` is used.  Therefore, @Baharis rather than implementing a libtbx.matplotlib shim, I recommend updating any XFEL code to the new matplotlib API.  Thanks!
AFAIK the `plt.get_cmap` is a syntax discouraged by matplotlib compared to `mpl.colormap['name']`, but it is safe across all matplotlib versions. I'll see about updating all the code in `xfel` to either one of those supported versions, probably the first one (though I find the latter more readable).
There are only a few places where `cm.get_cmap` needs to be updated, see #897. However, I have found multiple places where colormaps are accessed directly, e.g. `matplotlib.cm.jet`. I tried inquiring about it but it remains unclear to me whether this syntax is encouraged and expected to stay. It might be better in a long run to refactor those as well.

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 3628f06ffd..231eb19ac2 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -1,3 +1,11 @@
+2023.6
+======
+
+* Basic space group check when reading mmCIF
+* Fix hybrid36 output in mmCIF
+* Updates to lists of modified amino and nucleic acids
+* DataManager updates for handling multiple uses (types) for generating fmodel
+
 2023.5
 ======
 
diff --git a/cctbx/maptbx/prepare_map_for_docking.py b/cctbx/maptbx/prepare_map_for_docking.py
index 29b6f78dac..543ecdab27 100644
--- a/cctbx/maptbx/prepare_map_for_docking.py
+++ b/cctbx/maptbx/prepare_map_for_docking.py
@@ -1428,7 +1428,6 @@ def local_mean_density(mm, radius):
 
 def assess_cryoem_errors(
     mmm, d_min,
-    map_1_id="map_manager_1", map_2_id="map_manager_2",
     determine_ordered_volume=True,
     ordered_mask_id='ordered_volume_mask', sphere_points=500,
     sphere_cent=None, radius=None,
@@ -1444,9 +1443,6 @@ def assess_cryoem_errors(
     target region
 
   Optional arguments:
-  map_1_id: identifier of first half-map, if different from default of
-    map_manager_1
-  map_2_id: same for second half-map
   determine_ordered_volume: flag for whether ordered volume should be assessed
   ordered_mask_id: identifier for mask defining ordered volume
   sphere_cent: center of sphere defining target region for analysis
@@ -1479,7 +1475,16 @@ def assess_cryoem_errors(
     raise Sorry("Provided maps must be full-size")
   spacings = get_grid_spacings(unit_cell,unit_cell_grid)
   ucpars = unit_cell.parameters()
-  # Keep track of shifted origin
+
+  # Keep track of shifted origin, only allowing origin_shift_grid_units method
+  # that is specified by NXSTART, NYSTART, NZSTART in map header.
+  # Files in .mrc format can have non-zero ORIGIN (Angstrom units), referred to
+  # here as external_origin, but this must be redefined in grid units outside
+  # this program, e.g. with method implemented in em_placement_script.py.
+  external_origin = mmm.map_manager().external_origin
+  if (external_origin is not None) and (flex.sum(flex.abs(flex.double(external_origin))) > 0):
+    raise Sorry("Non-zero origin must be specified by origin_shift_grid_units, not Angstrom units (.mrc ORIGIN)")
+
   origin_shift = mmm.map_manager().origin_shift_grid_units
   if flex.sum(flex.abs(flex.double(origin_shift))) > 0:
     shifted_origin = True
@@ -1590,8 +1595,8 @@ def assess_cryoem_errors(
   minimum_possible_d_min = 2.*max(spacings)
   d_min = max(d_min, minimum_possible_d_min)
   d_min_extended = 1./(1./d_min + r_star)
-  mc1 = working_mmm.map_as_fourier_coefficients(d_min=d_min_extended, d_max=d_max, map_id=map_1_id)
-  mc2 = working_mmm.map_as_fourier_coefficients(d_min=d_min_extended, d_max=d_max, map_id=map_2_id)
+  mc1 = working_mmm.map_manager_1().map_as_fourier_coefficients(d_min=d_min_extended, d_max=d_max)
+  mc2 = working_mmm.map_manager_2().map_as_fourier_coefficients(d_min=d_min_extended, d_max=d_max)
 
   success = False
   while guess_d_min and not success:
@@ -1607,8 +1612,8 @@ def assess_cryoem_errors(
       else:
         success = True # d_min is already highest possible value
     d_min_extended = 1./(1./d_min + r_star)
-    mc1 = working_mmm.map_as_fourier_coefficients(d_min=d_min_extended, d_max=d_max, map_id=map_1_id)
-    mc2 = working_mmm.map_as_fourier_coefficients(d_min=d_min_extended, d_max=d_max, map_id=map_2_id)
+    mc1 = working_mmm.map_manager_1().map_as_fourier_coefficients(d_min=d_min_extended, d_max=d_max)
+    mc2 = working_mmm.map_manager_2().map_as_fourier_coefficients(d_min=d_min_extended, d_max=d_max)
 
   f1 = flex.abs(mc1.data())
   f2 = flex.abs(mc2.data())
@@ -1864,8 +1869,8 @@ def assess_cryoem_errors(
   wEmean2 = lweight*expectE
   working_mmm.add_map_from_fourier_coefficients(map_coeffs=wEmean2, map_id='map_manager_lwtd')
   # working_mmm.write_map(map_id='map_manager_wtd',file_name='prepmap.map')
-  working_mmm.remove_map_manager_by_id(map_1_id)
-  working_mmm.remove_map_manager_by_id(map_2_id)
+  working_mmm.remove_map_manager_by_id('map_manager_1')
+  working_mmm.remove_map_manager_by_id('map_manager_2')
 
   shift_cart = working_mmm.shift_cart()
   if shift_map_origin:
diff --git a/cctbx/maptbx/segment_and_split_map.py b/cctbx/maptbx/segment_and_split_map.py
index 957c399e43..dc9c1e0302 100644
--- a/cctbx/maptbx/segment_and_split_map.py
+++ b/cctbx/maptbx/segment_and_split_map.py
@@ -2899,12 +2899,17 @@ def find_symmetry_center(map_data, crystal_symmetry = None, out = sys.stdout):
   return xyz_cart
 
 def get_center_of_map(map_data, crystal_symmetry,
+    round_down = True,
     place_on_grid_point = True):
   all = list(map_data.all())
   origin = list(map_data.origin())
   if place_on_grid_point:
-    sx, sy, sz = [int(all[0]/2)+origin[0], int(all[1]/2)+origin[1],
-       int(all[2]/2)+origin[2]]
+    if round_down:
+      sx, sy, sz = [int(all[0]/2)+origin[0], int(all[1]/2)+origin[1],
+         int(all[2]/2)+origin[2]]
+    else:
+      sx, sy, sz = [int(1+all[0]/2)+origin[0], int(1+all[1]/2)+origin[1],
+         int(1+all[2]/2)+origin[2]]
   else:
     sx, sy, sz = [all[0]/2+origin[0], all[1]/2+origin[1],
   all[2]/2+origin[2]]
@@ -2983,19 +2988,35 @@ def run_get_ncs_from_map(params = None,
       wrapping=False)
     mm.resolution_filter(d_min=params.crystal_info.resolution)
     map_data = mm.map_data()
+  sym_cen = params.reconstruction_symmetry.symmetry_center
 
   ncs_obj_to_check = None
+
   if params.reconstruction_symmetry.symmetry and (
      not ncs_obj or ncs_obj.max_operators()<2):
     if params.reconstruction_symmetry.optimize_center:
-      center_try_list = [True, False]
+                       # [use_center_of_map, round_down, symmetry_center]
+      center_try_list = [[True,True,sym_cen],
+                         [False,True,sym_cen],
+                         [True,True,None],
+                         [True,False,None],
+                        ]
     else:
-      center_try_list = [True]
+      center_try_list = [[True,True,sym_cen],
+                         [True,True,None],
+                         [True,False,None],
+                        ]
   elif ncs_obj:
-    center_try_list = [True]
+    center_try_list = [[True,True,sym_cen],
+                         [True,True,None],
+                         [True,False,None],
+                        ]
     ncs_obj_to_check = ncs_obj
   elif params.reconstruction_symmetry.optimize_center:
-    center_try_list = [None]
+    center_try_list = [[True,True,sym_cen],
+                         [True,True,None],
+                         [True,False,None],
+                        ]
   else:
     return None, None, None # did not even try
   # check separately for helical symmetry
@@ -3010,9 +3031,10 @@ def run_get_ncs_from_map(params = None,
     helical_list = [False]
 
   new_ncs_obj, ncs_cc, ncs_score = None, None, None
-  for use_center_of_map in center_try_list:
+  for [use_center_of_map, round_down, symmetry_center] in center_try_list:
    for include_helical in helical_list:
     local_params = deepcopy(params)
+    local_params.reconstruction_symmetry.symmetry_center=None
     local_params.reconstruction_symmetry.include_helical_symmetry = \
        include_helical
     new_ncs_obj, ncs_cc, ncs_score = get_ncs_from_map(params = local_params,
@@ -3021,6 +3043,8 @@ def run_get_ncs_from_map(params = None,
       use_center_of_map_as_center = use_center_of_map,
       crystal_symmetry = crystal_symmetry,
       ncs_obj_to_check = ncs_obj_to_check,
+      symmetry_center = symmetry_center,
+      round_down = round_down,
       out = out
       )
     if new_ncs_obj:
@@ -3046,6 +3070,7 @@ def get_ncs_from_map(params = None,
       identify_ncs_id = None,
       ncs_obj_to_check = None,
       ncs_in_cell_only = False,
+      round_down = True,
       out = sys.stdout):
 
 
@@ -3058,6 +3083,8 @@ def get_ncs_from_map(params = None,
   #  If center is not supplied and use_center_of_map_as_center, try that
   #  and return None if it fails to achieve a map cc of min_ncs_cc
 
+  # round_down defines where to guess center if gridding has even number of pts
+
   if ncs_in_cell_only is None:
     ncs_in_cell_only = (not params.crystal_info.use_sg_symmetry)
   if symmetry is None:
@@ -3088,10 +3115,9 @@ def get_ncs_from_map(params = None,
   if ncs_obj_to_check and ncs_obj_to_check.max_operators()>1:
     symmetry = "SUPPLIED NCS"
 
-
   if map_symmetry_center is None:
-    map_symmetry_center = get_center_of_map(map_data, crystal_symmetry)
-
+    map_symmetry_center = get_center_of_map(map_data, crystal_symmetry,
+      round_down = round_down)
   if optimize_center is None:
     if symmetry_center is None and (not use_center_of_map_as_center):
       optimize_center = True
@@ -3157,20 +3183,27 @@ def get_ncs_from_map(params = None,
   if n_rescore and not ncs_obj_to_check:
     print("Rescoring top %d results" %(min(n_rescore, len(results_list))), file = out)
     rescore_list = results_list[n_rescore:]
+    all_rescore = []
     new_sites_orth = get_points_in_map(
       map_data, n = 10*random_points, crystal_symmetry = crystal_symmetry)
     new_sites_orth.extend(sites_orth)
+    found_ok = False
     for orig_score, orig_cc_avg, ncs_obj, symmetry in results_list[:n_rescore]:
       score, cc_avg = score_ncs_in_map(map_data = map_data, ncs_object = ncs_obj,
         identify_ncs_id = identify_ncs_id,
         ncs_in_cell_only = ncs_in_cell_only,
         sites_orth = new_sites_orth, crystal_symmetry = crystal_symmetry, out = out)
+      all_rescore.append([score, cc_avg, ncs_obj, symmetry])
       if cc_avg is None or cc_avg < min_ncs_cc:
         score = 0. # Do not allow low CC values to be used
+      else:
+        found_ok = True
       if score is None:
         print("symmetry:", symmetry, " no score", ncs_obj.max_operators(), file = out)
       else:
         rescore_list.append([score, cc_avg, ncs_obj, symmetry])
+    if not found_ok:
+      rescore_list = all_rescore
     rescore_list.sort(key=itemgetter(0))
     rescore_list.reverse()
     results_list = rescore_list
@@ -3197,6 +3230,24 @@ def get_ncs_from_map(params = None,
        score, cc_avg, ncs_obj.max_operators(), symmetry.strip(), ), file = out)
 
   score, cc_avg, ncs_obj, ncs_info = results_list[0]
+  # check for offset by gridding
+  if  hasattr(params.reconstruction_symmetry,'check_grid_offset') and \
+     params.reconstruction_symmetry.check_grid_offset:
+    symmetry_center, cc_avg, score, ncs_obj = optimize_center_position(
+       map_data, sites_orth,
+       crystal_symmetry,
+       ncs_info, symmetry_center, ncs_obj, score, cc_avg,
+       params = params,
+       helical_rot_deg = helical_rot_deg,
+       two_fold_along_x = two_fold_along_x,
+       op_max = op_max,
+       min_ncs_cc = min_ncs_cc,
+       identify_ncs_id = identify_ncs_id,
+       ncs_obj_to_check = ncs_obj_to_check,
+       ncs_in_cell_only = ncs_in_cell_only,
+       helical_trans_z_angstrom = helical_trans_z_angstrom,
+       check_grid_offset = True,
+       out = out)
 
   # Optimize center if necessary
   if optimize_center:
@@ -3215,7 +3266,7 @@ def get_ncs_from_map(params = None,
        helical_trans_z_angstrom = helical_trans_z_angstrom, out = out)
     print("New center: (%7.3f, %7.3f, %7.3f)" %(tuple(symmetry_center)), file = out)
 
-  if cc_avg < min_ncs_cc:
+  if (not cc_avg) or (cc_avg < min_ncs_cc):
     print("No suitable symmetry found", file = out)
     return None, None, None
 
@@ -3237,12 +3288,20 @@ def optimize_center_position(map_data, sites_orth, crystal_symmetry,
      ncs_obj_to_check = None,
      ncs_in_cell_only = None,
      min_ncs_cc = None,
-     helical_trans_z_angstrom = None, out = sys.stdout):
+     helical_trans_z_angstrom = None,
+     check_grid_offset = False,
+     out = sys.stdout):
 
   if ncs_info is None:
     ncs_info = "None"
+
   symmetry = ncs_info.split()[0]
-  print("Optimizing center position...type is %s" %(ncs_info), file = out)
+  if check_grid_offset:
+    print(
+      "Checking for improvement with offset by +/-1 grid unit...type is %s" %(
+       ncs_info), file = out)
+  else:
+    print("Optimizing center position...type is %s" %(ncs_info), file = out)
 
   if len(ncs_info.split())>1 and ncs_info.split()[1] == '(a)':
     two_fold_along_x = True
@@ -3255,14 +3314,36 @@ def optimize_center_position(map_data, sites_orth, crystal_symmetry,
   best_ncs_obj = ncs_obj
   best_score = score
   best_cc_avg = cc_avg
+  if (ncs_info.find("Helical") > -1) and check_grid_offset:
+    return best_center, best_cc_avg, best_score, best_ncs_obj
   print("Starting center: (%7.3f, %7.3f, %7.3f)" %(tuple(best_center)), file = out)
+  if check_grid_offset:
+    n_range = 1
+    k_range = list(range(-1,2))
+    i_range = list(range(-1,2))
+    j_range = list(range(-1,2))
+    abc = crystal_symmetry.unit_cell().parameters()[:3]
+    N_ = map_data.all()
+
+    scale_x = abc[0]/N_[0]
+    scale_y = abc[1]/N_[1]
+    scale_z = abc[2]/N_[2]
+  else:
+    n_range = 6
+    k_range = range(0,1)
+    i_range = range(-4,5)
+    j_range = range(-4,5)
+    scale_x = 1
+    scale_y = 1
+    scale_z = 1
   from libtbx.utils import null_out
   scale = 5.
-  for itry in range(6):
-    scale = scale/5.
-    for i in range(-4, 5):
-     for j in range(-4, 5):
-      local_center = matrix.col(symmetry_center)+matrix.col((scale*i, scale*j, 0., ))
+  for itry in range(n_range):
+   scale = scale/5.
+   for k in k_range:
+    for i in i_range:
+     for j in j_range:
+      local_center = matrix.col(symmetry_center)+matrix.col((scale_x*i, scale_y*j, scale_z*k, ))
       ncs_list = get_ncs_list(params = params, symmetry = symmetry,
        symmetry_center = local_center,
        helical_rot_deg = helical_rot_deg,
@@ -3295,6 +3376,7 @@ def optimize_center_position(map_data, sites_orth, crystal_symmetry,
   cc_avg = best_cc_avg
   score = best_score
   ncs_obj = best_ncs_obj
+  print("Optimized center: (%7.3f, %7.3f, %7.3f)" %(tuple(best_center)), file = out)
   return best_center, best_cc_avg, best_score, best_ncs_obj
 
 
diff --git a/crys3d/hklviewer/hklview_frame.py b/crys3d/hklviewer/hklview_frame.py
index 67bcd8eb61..9cdbace2ba 100644
--- a/crys3d/hklviewer/hklview_frame.py
+++ b/crys3d/hklviewer/hklview_frame.py
@@ -168,7 +168,7 @@ def thread_process_arguments(self, **kwds):
         self.mprint("Processing PHIL file: %s" %fname)
         with open(fname, "r") as f:
           philstr = f.read()
-          self.update_from_philstr(philstr)
+          self.update_from_philstr(philstr, msgtype="preset_philstr", postrender=True)
           self.mprint("thread_process_arguments() waiting for run_external_sem.acquire", verbose="threadingmsg")
           if not self.run_external_sem.acquire(timeout=jsview_3d.lock_timeout):
             self.mprint("Timed out getting run_external_sem semaphore within %s seconds" %jsview_3d.lock_timeout, verbose=1)
@@ -177,7 +177,7 @@ def thread_process_arguments(self, **kwds):
         if 'image_file' in kwds: # save displayed reflections to an image file
           time.sleep(1)
           fname = kwds.get('image_file', "testimage.png" )
-          self.update_from_philstr('save_image_name = "%s"' %fname)
+          self.update_from_philstr('save_image_name = "%s"' %fname, msgtype="preset_philstr", postrender=True)
       # if we are invoked using Qtgui close us gracefully if requested
       if 'closing_time' in kwds:
         time.sleep(10)
@@ -462,11 +462,11 @@ def show_current_phil(self, useful_for_preset_button=True):
     return "\nCurrent non-default phil parameters:\n\n" + diffphil.as_str()
 
 
-  def update_from_philstr(self, philstr):
+  def update_from_philstr(self, philstr, msgtype="philstr", postrender=False):
     # Convenience function for scripting HKLviewer that mostly superseedes other functions for
     # scripting such as ExpandAnomalous(True), SetScene(0) etc.
     new_phil = libtbx.phil.parse(philstr)
-    self.guarded_process_PHIL_parameters(new_phil, msgtype="preset_philstr", postrender=True)
+    self.guarded_process_PHIL_parameters(new_phil, msgtype=msgtype, postrender=postrender)
 
 
   def guarded_process_PHIL_parameters(self, new_phil=None, msgtype="philstr",
diff --git a/iotbx/bioinformatics/structure.py b/iotbx/bioinformatics/structure.py
index e319fe5fae..27093e2e33 100644
--- a/iotbx/bioinformatics/structure.py
+++ b/iotbx/bioinformatics/structure.py
@@ -68,10 +68,8 @@ def summarize_blast_output(blast_out=None, blast_file=None,
   blast = next(parsed)
   if (len(blast.alignments) == 0):
     if stop_if_no_alignment:
-      raise StopIteration("No matching sequences!")
+      raise Sorry("No matching sequences!")
     else: return list()
-  #if (len(blast.alignments) == 0):
-  #  raise StopIteration("No matching sequences!")
   results = []
   for i_hit, hit in enumerate(blast.alignments):
     pdb_chain_id = str(hit.accession)
diff --git a/iotbx/map_manager.py b/iotbx/map_manager.py
index f5522a08e5..40e2efab23 100644
--- a/iotbx/map_manager.py
+++ b/iotbx/map_manager.py
@@ -48,6 +48,10 @@ class map_manager(map_reader, write_ccp4_map):
    shifted to place the origin at (0, 0, 0) and the original origin is
    recorded as self.origin_shift_grid_units.
 
+   NOTE: do not set self.origin_shift_grid_units directly; instead use
+   the method set_original_origin_and_gridding() or supply the value when
+   initializing the map_manager.
+
    You can also create a map_manager with a map_data object (3D flex.double()
    array) along with the meta-data below.
 
@@ -221,6 +225,15 @@ def __init__(self,
         scattering_table (electron n_gaussian wk1995 it1992 neutron)
         resolution (nominal resolution of map)
 
+     Note on external origin: If input map had ORIGIN specified,
+       so that the value of self.external_origin is not (0,0,0) and not None,
+       and apply_external_origin_if_present is set, and shift_origin() is
+       called, then:
+       determine if self.external_origin is on a grid point and if so, convert
+        and use negative of it as origin .
+       NOTE: apply_external_origin_if_present will be ignored if
+       origin_shift_grid_units is also set and is not (0,0,0).
+
       NOTE on "crystal_symmetry" objects
       There are two objects that are "crystal_symmetry" objects:
       A.  unit_cell_crystal_symmetry():  This is the symmetry of the
@@ -251,6 +264,8 @@ def __init__(self,
     if origin_shift_grid_units is not None:
       origin_shift_grid_units = tuple(origin_shift_grid_units)
       assert len(origin_shift_grid_units) ==3
+    else:
+      origin_shift_grid_units = (0, 0, 0)
 
     # Initialize log filestream
     self.set_log(log)
@@ -327,6 +342,7 @@ def __init__(self,
       self.unit_cell_grid = unit_cell_grid
       self._unit_cell_crystal_symmetry = unit_cell_crystal_symmetry
       self._wrapping = wrapping
+      self.external_origin = (0, 0, 0)
 
       # Calculate values for self._crystal_symmetry
       # Must always run this method after changing
@@ -334,8 +350,6 @@ def __init__(self,
       self.set_crystal_symmetry_of_partial_map()
 
       # Optional initialization information
-      if origin_shift_grid_units is None:
-        origin_shift_grid_units = (0, 0, 0)
       self.origin_shift_grid_units = origin_shift_grid_units
 
     # Initialization steps always done:
@@ -539,10 +553,13 @@ def origin_is_zero(self):
     else:
       return False
 
-  def shift_origin(self, desired_origin = (0, 0, 0)):
+  def shift_origin(self, desired_origin = (0, 0, 0),
+     apply_external_origin_if_present = True,):
+
     '''
     Shift the origin of the map to desired_origin
-        (normally desired_origin = (0, 0, 0) and update origin_shift_grid_units
+      (normally desired_origin is (0, 0, 0), so just update
+      origin_shift_grid_units
 
     Origin is the value of self.map_data().origin()
     origin_shift_grid_units is the shift to apply to self.map_data() to
@@ -563,17 +580,27 @@ def shift_origin(self, desired_origin = (0, 0, 0)):
 
     the new origin of map_data will be (d, e, f)
 
-    '''
+     Note on external origin: If input map had ORIGIN specified,
+       so that the value of self.external_origin is not (0,0,0) and not None,
+       and apply_external_origin_if_present is set, then:
+       determine if self.external_origin is on a grid point and if so, convert
+        and use negative of it as origin. Then self.external_origin to zero.
+       Does not apply if origin is already not (0,0,0).
 
+    '''
     if(self.map_data() is None): return
 
     # Figure out what the shifts are (in grid units)
-    shift_info = self._get_shift_info(desired_origin = desired_origin)
+    shift_info = self._get_shift_info(desired_origin = desired_origin,
+      apply_external_origin_if_present = apply_external_origin_if_present)
 
     # Update the value of origin_shift_grid_units
     #  This is position of the origin of the new map in the full unit cell grid
     self.origin_shift_grid_units = shift_info.new_origin_shift_grid_units
 
+    # Set external_origin to zero as it has now been used
+    self.external_origin = (0, 0, 0)
+
     # Shift map_data if necessary
     if shift_info.shift_to_apply !=  (0, 0, 0):
       # map will start at desired_origin and have current size:
@@ -597,7 +624,8 @@ def shift_origin(self, desired_origin = (0, 0, 0)):
       self._ncs_object=self._ncs_object.coordinate_offset(
         shift_info.shift_to_apply_cart)
 
-  def _get_shift_info(self, desired_origin = None):
+  def _get_shift_info(self, desired_origin = None,
+    apply_external_origin_if_present = None):
     '''
       Utility to calculate the shift necessary (grid units)
       map to place the origin of the current map
@@ -615,6 +643,34 @@ def _get_shift_info(self, desired_origin = None):
     # Current origin and shift to apply
     current_origin = self.map_data().origin()
 
+    self._warning_message = ""
+    if apply_external_origin_if_present and \
+         tuple(self.external_origin) != (0,0,0): # check for external origin
+      if self.external_origin_is_compatible_with_gridding():
+         external_origin_as_grid_units = self.external_origin_as_grid_units()
+      else:
+        external_origin_as_grid_units = (0,0,0)
+        self._warning_message="External origin is not on a grid point" +\
+         "...ignoring external origin" +\
+         "\n***Please contact the Phenix "+\
+          "developers if you need Phenix to use this external_origin***\n"
+    else:
+      external_origin_as_grid_units = (0,0,0)
+
+    if self.external_origin_as_grid_units and \
+        (external_origin_as_grid_units != (0,0,0)):
+      if current_origin and \
+          (current_origin != (0,0,0)):
+        self._warning_message="Map has external origin as well as existing " +\
+         "origin shift...ignoring external origin" +\
+         "\n***Please contact the Phenix "+\
+          "developers if you need Phenix to use this external_origin***\n"
+      else:  # take it
+        self._warning_message="Map has external origin " +\
+         "...using external origin as origin shift after "+\
+         "conversion to grid units"
+        current_origin = external_origin_as_grid_units
+
     # Original location of first element of map
     map_corner_original_location = add_tuples_int(current_origin,
          self.origin_shift_grid_units)
@@ -644,6 +700,38 @@ def _get_shift_info(self, desired_origin = None):
        )
     return shift_info
 
+  def external_origin_is_compatible_with_gridding(self):
+    value = self.external_origin_as_grid_units()
+    if value is not None:
+      return True
+    else:
+      return False
+
+  def external_origin_as_grid_units(self, as_inverse = False):
+    unit_cell = self.unit_cell_crystal_symmetry().unit_cell()
+    unit_cell_grid = self.unit_cell_grid
+    spacings = [(a/n) for a,n in zip(unit_cell.parameters()[:3],
+        unit_cell_grid)]
+    if self.external_origin:
+      external_origin = flex.double(self.external_origin)
+    else:
+      external_origin = flex.double((0.,0.,0.))
+    if flex.sum(flex.abs(external_origin)) > 0:
+      import math
+      origin_shift = tuple([round(o/s) for o,s in zip(external_origin, spacings)])
+      origin_check = [(g*s) for g,s in zip(origin_shift, spacings)]
+      origin_distance_to_grid = math.sqrt(flex.sum(
+            flex.pow2(flex.double(external_origin)-flex.double(origin_check))))
+      if origin_distance_to_grid > 0.001:
+        return None # not compatible with grid
+      else:
+        if as_inverse:
+          return tuple([-a for a in origin_shift])
+        else: # as-is:
+          return origin_shift
+    else: # no shift
+      return (0,0,0)
+
   def shift_origin_to_match_original(self):
     '''
      Shift origin by self.origin_shift_grid_units to place origin in its
diff --git a/iotbx/map_model_manager.py b/iotbx/map_model_manager.py
index af3386a3cb..f4ea3c0ab6 100644
--- a/iotbx/map_model_manager.py
+++ b/iotbx/map_model_manager.py
@@ -287,6 +287,7 @@ def __init__(self,
         absolute_angle_tolerance = absolute_angle_tolerance,
         absolute_length_tolerance = absolute_length_tolerance,
         ignore_symmetry_conflicts = ignore_symmetry_conflicts)
+    mmmn.set_log(self.log)
     mmmn.add_map_manager(any_map_manager)
     for m in extra_model_list: # Check all extra models
       mmmn.add_model(m.deep_copy(),
@@ -513,7 +514,7 @@ def set_log(self, log = sys.stdout):
     else:
       self.log = log
 
-  def _print(self, m):
+  def _print(self, m, force = False):
     '''
       Print to log if it is present
     '''
@@ -521,6 +522,8 @@ def _print(self, m):
     if (self.log is not None) and hasattr(self.log, 'closed') and (
         not self.log.closed):
       print(m, file = self.log)
+    elif force:
+      print(m)
 
   # Methods for obtaining models, map_managers, symmetry, ncs_objects
 
@@ -8670,10 +8673,12 @@ def set_log(self, log = sys.stdout):
     else:
       self.log = log
 
-  def _print(self, m):
+  def _print(self, m, force = False):
     if (self.log is not None) and hasattr(self.log, 'closed') and (
         not self.log.closed):
-      self._print(m, file = self.log)
+      print(m, file = self.log)
+    elif force:
+      print(m)
 
   def write_map(self, file_name = None):
     if not self._map_manager:
@@ -8873,8 +8878,11 @@ def shift_origin(self, desired_origin = (0, 0, 0)):
 
     # Apply shift to model, map and ncs object
 
-    # Shift origin of map_manager
+    # Shift origin of map_manager. Note if there was any problem
+    self._map_manager._warning_message = ''
     self._map_manager.shift_origin(desired_origin = desired_origin)
+    if self._map_manager.warning_message():
+      self._print("\n%s\n" %(self._map_manager.warning_message()), force=True)
 
     # Shift origin of model  Note this sets model shift_cart
     if self._model:
diff --git a/iotbx/merging_statistics.py b/iotbx/merging_statistics.py
index daaf1371fa..355e34e444 100644
--- a/iotbx/merging_statistics.py
+++ b/iotbx/merging_statistics.py
@@ -823,7 +823,7 @@ def as_cif_block(self, cif_block=None):
 
     for k, v in six.iteritems(mmcif_to_name_reflns):
       value = self.overall.__getattribute__(v)
-      if 'percent' in v:
+      if 'percent' in k:
         value *= 100
       cif_block['_reflns.' + k] = value
 
diff --git a/iotbx/mrcfile/__init__.py b/iotbx/mrcfile/__init__.py
index f881ac24ff..6104f600d4 100644
--- a/iotbx/mrcfile/__init__.py
+++ b/iotbx/mrcfile/__init__.py
@@ -51,8 +51,9 @@ class map_reader:
      NOTE phenix calls "origin" the position of the lower left corner
        of the map.
      mrcfile calls "origin" the value of the field "origin" which is 3
-       real numbers indicating the placement of the grid point (0,0,0) relative
-       to an external reference frame (typically that of a model)
+       real numbers (A units, not grid units) indicating the placement
+       of the grid point (0,0,0) relative to an external reference frame
+       (typically that of a big map that this map is cut out from).
 
      Here we will use "external_origin" to refer to the mrc origin and
      "origin" for nxstart_nystart_nzstart.
diff --git a/iotbx/pdb/common_residue_names.h b/iotbx/pdb/common_residue_names.h
index 3a8197f328..9d20da50df 100644
--- a/iotbx/pdb/common_residue_names.h
+++ b/iotbx/pdb/common_residue_names.h
@@ -356,139 +356,117 @@ namespace iotbx { namespace pdb { namespace common_residue_names {
   inline
   void
   initialize_set(
-    std::set<str3>& name_set,
+    std::set<std::string>& name_set,
     const char** names)
   {
     if (name_set.size() != 0) return;
-    for(const char** n=names; *n; n++) name_set.insert(str3(*n));
+    for(const char** n=names; *n; n++) name_set.insert(std::string(*n));
   }
 
   inline
-  str3
-  std_string_as_str3(std::string const& name)
-  {
-    if (name.size() > 3) {
-      throw std::runtime_error(
-        "residue name with more than 3 characters: \""+name+"\"");
-    }
-    return str3(name.data(), name.size(), /*i_begin*/ 0, /*pad_with*/ ' ');
-  }
-
-  inline
-  void
-  initialize_set(
-    std::set<str3>& name_set,
-    scitbx::af::const_ref<std::string> const& names)
-  {
-    for(std::size_t i=0;i<names.size();i++) {
-      name_set.insert(std_string_as_str3(names[i]));
-    }
-  }
-
-  inline
-  const std::set<str3>&
+  const std::set<std::string>&
   amino_acid_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, amino_acid);
     return result;
   }
 
   inline
-  const std::set<str3>&
+  const std::set<std::string>&
   d_amino_acid_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, d_amino_acid);
     return result;
   }
 
   inline
-  const std::set<str3>&
+  const std::set<std::string>&
   modified_amino_acid_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, modified_amino_acid);
     return result;
   }
 
   inline
-  const std::set<str3>&
+  const std::set<std::string>&
   rna_dna_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, rna_dna);
     return result;
   }
 
   inline
-  const std::set<str3>&
+  const std::set<std::string>&
   modified_rna_dna_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, modified_rna_dna);
     return result;
   }
 
   inline
-  const std::set<str3>&
+  const std::set<std::string>&
   ccp4_mon_lib_rna_dna_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, ccp4_mon_lib_rna_dna);
     return result;
   }
 
   inline
-  const std::set<str3>&
+  const std::set<std::string>&
   water_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, water);
     return result;
   }
 
   inline
-  const std::set<str3>&
+  const std::set<std::string>&
   small_molecule_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, small_molecule);
     return result;
   }
 
   inline
-  const std::set<str3>&
+  const std::set<std::string>&
   saccharide_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, common_saccharide);
     return result;
   }
 
   inline
-  const std::set<str3>&
+  const std::set<std::string>&
   element_set()
   {
-    static std::set<str3> result;
+    static std::set<std::string> result;
     initialize_set(result, element);
     return result;
   }
 
   inline
   std::string const&
-  get_class(str3 const& name, bool consider_ccp4_mon_lib_rna_dna=false)
+  get_class(std::string const& name, bool consider_ccp4_mon_lib_rna_dna=false)
   {
-    static const std::set<str3>& aa_set = amino_acid_set();
-    static const std::set<str3>& d_aa_set = d_amino_acid_set();
-    static const std::set<str3>& modified_aa_set = modified_amino_acid_set();
-    static const std::set<str3>& na_set = rna_dna_set();
-    static const std::set<str3>& modified_na_set = modified_rna_dna_set();
-    static const std::set<str3>& ml_na_set = ccp4_mon_lib_rna_dna_set();
-    static const std::set<str3>& w_set = water_set();
-    static const std::set<str3>& sm_set = small_molecule_set();
-    static const std::set<str3>& s_set = saccharide_set();
-    static const std::set<str3>& e_set = element_set();
+    static const std::set<std::string>& aa_set = amino_acid_set();
+    static const std::set<std::string>& d_aa_set = d_amino_acid_set();
+    static const std::set<std::string>& modified_aa_set = modified_amino_acid_set();
+    static const std::set<std::string>& na_set = rna_dna_set();
+    static const std::set<std::string>& modified_na_set = modified_rna_dna_set();
+    static const std::set<std::string>& ml_na_set = ccp4_mon_lib_rna_dna_set();
+    static const std::set<std::string>& w_set = water_set();
+    static const std::set<std::string>& sm_set = small_molecule_set();
+    static const std::set<std::string>& s_set = saccharide_set();
+    static const std::set<std::string>& e_set = element_set();
     static const std::string common_amino_acid("common_amino_acid");
     static const std::string d_amino_acid("d_amino_acid");
     static const std::string modified_amino_acid("modified_amino_acid");
@@ -500,47 +478,44 @@ namespace iotbx { namespace pdb { namespace common_residue_names {
     static const std::string common_saccharide("common_saccharide");
     static const std::string common_element("common_element");
     static const std::string other("other");
-    if (aa_set.find(name) != aa_set.end()) {
+    std::string padded(name);
+    if (padded.size() < 3) {
+      padded.insert(padded.begin(), 3-padded.size(), ' ');
+    }
+    if (aa_set.find(padded) != aa_set.end()) {
       return common_amino_acid;
     }
-    if (d_aa_set.find(name) != d_aa_set.end()) {
+    if (d_aa_set.find(padded) != d_aa_set.end()) {
       return d_amino_acid;
     }
-    if (modified_aa_set.find(name) != modified_aa_set.end()) {
+    if (modified_aa_set.find(padded) != modified_aa_set.end()) {
       return modified_amino_acid;
     }
-    if (na_set.find(name) != na_set.end()) {
+    if (na_set.find(padded) != na_set.end()) {
       return common_rna_dna;
     }
-    if (modified_na_set.find(name) != modified_na_set.end()) {
+    if (modified_na_set.find(padded) != modified_na_set.end()) {
       return modified_rna_dna;
     }
     if (   consider_ccp4_mon_lib_rna_dna
-        && ml_na_set.find(name) != ml_na_set.end()) {
+        && ml_na_set.find(padded) != ml_na_set.end()) {
       return ccp4_mon_lib_rna_dna;
     }
-    if (w_set.find(name) != w_set.end()) {
+    if (w_set.find(padded) != w_set.end()) {
       return common_water;
     }
-    if (sm_set.find(name) != sm_set.end()) {
+    if (sm_set.find(padded) != sm_set.end()) {
       return common_small_molecule;
     }
-    if (s_set.find(name) != s_set.end()) {
+    if (s_set.find(padded) != s_set.end()) {
       return common_saccharide;
     }
-    if (e_set.find(name) != e_set.end()) {
+    if (e_set.find(padded) != e_set.end()) {
       return common_element;
     }
     return other;
   }
 
-  inline
-  std::string const&
-  get_class(std::string const& name, bool consider_ccp4_mon_lib_rna_dna=false)
-  {
-    return get_class(std_string_as_str3(name), consider_ccp4_mon_lib_rna_dna);
-  }
-
 }}} // namespace iotbx::pdb::common_residue_names
 
 #endif // IOTBX_PDB_COMMON_RESIDUE_NAMES_H
diff --git a/iotbx/pdb/hierarchy.cpp b/iotbx/pdb/hierarchy.cpp
index eb6181a328..02aa4cd386 100644
--- a/iotbx/pdb/hierarchy.cpp
+++ b/iotbx/pdb/hierarchy.cpp
@@ -628,15 +628,16 @@ namespace {
     return std::string(result, 15U);
   }
 
-  small_str<19>
+  std::string
   atom::pdb_label_columns_segid_small_str() const
+  // Only used in overall_counts in find_duplicate_atom_labels()
   {
     char blank = ' ';
     small_str<19> result(small_str_no_init);
     atom_label_columns_formatter().format(result.elems, *this);
     data->segid.copy_left_justified(result.elems+15, 4U, blank);
     result.elems[19] = '\0';
-    return result;
+    return std::string(result.elems);
   }
 
   std::string
@@ -1343,7 +1344,7 @@ namespace {
     std::string ag_class;
     // Using correct overloaded function so residue name would be padded.
     // Crucial for mmCIF.
-    ag_class = common_residue_names::get_class(std::string(data->resname.c_str()));
+    ag_class = common_residue_names::get_class(data->resname);
     if (ag_class == "common_rna_dna" || ag_class == "modified_rna_dna") {
       if (!get_atom("N9")) {
         std::sort(
diff --git a/iotbx/pdb/hierarchy.h b/iotbx/pdb/hierarchy.h
index 446b6d2412..634d3c0b66 100644
--- a/iotbx/pdb/hierarchy.h
+++ b/iotbx/pdb/hierarchy.h
@@ -709,7 +709,7 @@ namespace hierarchy {
       pdb_label_columns() const;
 
       //! Not available in Python.
-      small_str<19>
+      std::string
       pdb_label_columns_segid_small_str() const;
 
       std::string
diff --git a/iotbx/pdb/overall_counts.cpp b/iotbx/pdb/overall_counts.cpp
index e241bff45d..1790802280 100644
--- a/iotbx/pdb/overall_counts.cpp
+++ b/iotbx/pdb/overall_counts.cpp
@@ -9,10 +9,10 @@ namespace detail {
 
   struct cmp_atom_labels_functor
   {
-    const small_str<19>* labels;
+    std::vector<std::string> const& labels;
 
     cmp_atom_labels_functor(
-      const small_str<19>* labels_)
+      std::vector<std::string> const& labels_)
     :
       labels(labels_)
     {}
@@ -22,7 +22,7 @@ namespace detail {
       unsigned i,
       unsigned j) const
     {
-      return (std::memcmp(labels[i].elems, labels[j].elems, 19U) < 0);
+      return (labels[i] < labels[j]);
     }
   };
 
@@ -39,7 +39,7 @@ namespace detail {
     af::shared<af::shared<atom> >& duplicate_atom_labels,
     hierarchy::model const& model,
     unsigned model_atoms_size,
-    const small_str<19>* model_atom_labels)
+    std::vector<std::string> const& model_atom_labels)
   {
     if (model_atoms_size == 0) return 0;
     boost::scoped_array<unsigned> indices(new unsigned[model_atoms_size]);
@@ -55,9 +55,7 @@ namespace detail {
     unsigned j_start = 0;
     for(unsigned j=1;j<model_atoms_size+1U;j++) {
       if (j != model_atoms_size) {
-        if (std::memcmp(
-              model_atom_labels[indices[j_start]].elems,
-              model_atom_labels[indices[j      ]].elems, 19U) == 0) continue;
+        if (model_atom_labels[indices[j_start]] == model_atom_labels[indices[j]]) continue;
       }
       if (j_start+1U == j) {
         j_start++;
@@ -125,8 +123,7 @@ namespace detail {
       model_ids[model.data->id]++;
       std::map<std::string, unsigned> model_chain_ids;
       unsigned model_atoms_size = model.atoms_size();
-      boost::scoped_array<small_str<19> > model_atom_labels(
-        new small_str<19>[model_atoms_size]);
+      std::vector<std::string > model_atom_labels;
       unsigned i_model_atom = 0;
       unsigned n_ch = model.chains_size();
       n_chains += n_ch;
@@ -173,8 +170,9 @@ namespace detail {
             for(unsigned i_at=0;i_at<n_ats;i_at++) {
               hierarchy::atom const& atom = ag.atoms()[i_at];
               if (atom.uij_is_defined()) n_anisou++;
-              model_atom_labels[i_model_atom++] =
-                atom.pdb_label_columns_segid_small_str();
+              model_atom_labels.push_back(atom.pdb_label_columns_segid_small_str());
+              i_model_atom++;
+              // std::cout << "'" << atom.pdb_label_columns_segid_small_str() << "'\n";
               element_charge_types[atom.pdb_element_charge_columns()]++;
             }
           }
@@ -264,7 +262,7 @@ namespace detail {
         duplicate_atom_labels,
         model,
         model_atoms_size,
-        model_atom_labels.get());
+        model_atom_labels);
     }
     {
       typedef std::map<std::string, unsigned>::const_iterator it;
diff --git a/iotbx/pdb/tst_hierarchy.py b/iotbx/pdb/tst_hierarchy.py
index fcf7f83c4e..1356476772 100644
--- a/iotbx/pdb/tst_hierarchy.py
+++ b/iotbx/pdb/tst_hierarchy.py
@@ -479,6 +479,10 @@ def exercise_atom_group():
   assert ag.altloc == ""
   assert ag.resname == ""
   assert ag.confid() == "  ", "'%s'" % ag.confid()
+  ag = pdb.hierarchy.atom_group(altloc="", resname="5char")
+  assert ag.altloc == ""
+  assert ag.resname == "5char"
+  assert ag.confid() == " 5char", "'%s'" % ag.confid()
   #
   ag.altloc = "l"
   ag.resname = "res"
@@ -773,6 +777,8 @@ def exercise_chain():
   assert c.id == ""
   c = pdb.hierarchy.chain(id="a")
   assert c.id == "a"
+  c = pdb.hierarchy.chain(id="long_chain_id")
+  assert c.id == "long_chain_id"
   c.id = "x"
   assert c.id == "x"
   c.id = None
@@ -4502,6 +4508,65 @@ def exercise_is_similar_hierarchy():
       .is_similar_hierarchy(
         other=h2.models()[0].only_chain().residue_groups()[0]) == (an == "N ")
 
+def exercise_is_similar_hierarchy_long():
+  s0 = """\
+loop_
+_atom_site.group_PDB
+_atom_site.id
+_atom_site.type_symbol
+_atom_site.label_atom_id
+_atom_site.label_alt_id
+_atom_site.label_comp_id
+_atom_site.label_asym_id
+_atom_site.label_entity_id
+_atom_site.label_seq_id
+_atom_site.pdbx_PDB_ins_code
+_atom_site.Cartn_x
+_atom_site.Cartn_y
+_atom_site.Cartn_z
+_atom_site.occupancy
+_atom_site.B_iso_or_equiv
+_atom_site.Cartn_x_esd
+_atom_site.Cartn_y_esd
+_atom_site.Cartn_z_esd
+_atom_site.occupancy_esd
+_atom_site.B_iso_or_equiv_esd
+_atom_site.pdbx_formal_charge
+_atom_site.auth_seq_id
+_atom_site.auth_comp_id
+_atom_site.auth_asym_id
+_atom_site.auth_atom_id
+_atom_site.pdbx_PDB_model_num
+ATOM   9414   C CA   . %s A-2  1 1   ? 71.805  71.447  63.447  1.00 78.79  ? ? ? ? ? ? 2   SER A-2  CA   1
+ATOM   9414   C CA   . %s A-2  2 2   ? 71.805  71.447  63.447  1.00 78.79  ? ? ? ? ? ? 2   SER A-2  CA   1
+ATOM   9414   C CA   . %s A-2  3 3   ? 71.805  71.447  63.447  1.00 78.79  ? ? ? ? ? ? 2   SER A-2  CA   1
+"""
+
+  i1 = pdb.input(source_info=None, lines=flex.split_lines(
+    s0 % ("SERine", "SERine", "SERine")))
+  h1 = i1.construct_hierarchy()
+  assert h1.is_similar_hierarchy(other=h1)
+  assert h1.is_similar_hierarchy(other=h1.deep_copy())
+  assert h1.models()[0].is_similar_hierarchy(
+    other=h1.models()[1])
+  assert h1.models()[0].only_chain().is_similar_hierarchy(
+    other=h1.models()[1].only_chain())
+  assert h1.models()[0].only_chain().residue_groups()[0].is_similar_hierarchy(
+    other=h1.models()[1].only_chain().residue_groups()[0])
+  for rn in ["SER","Alanine"]:
+    i2 = pdb.input(source_info=None, lines=flex.split_lines(
+      s0 % (rn, rn, rn)))
+    h2 = i2.construct_hierarchy()
+    assert not h1.is_similar_hierarchy(other=h2)
+    assert not h2.is_similar_hierarchy(other=h1)
+    assert not h1.models()[0].is_similar_hierarchy(
+      other=h2.models()[0]) == (rn == "GLY")
+    assert h1.models()[0].only_chain().is_similar_hierarchy(
+      other=h2.models()[0].only_chain()) == (an == "N ")
+    assert h1.models()[0].only_chain().residue_groups()[0] \
+      .is_similar_hierarchy(
+        other=h2.models()[0].only_chain().residue_groups()[0]) == (an == "N ")
+
 def exercise_atoms():
   pdb_inp = pdb.input(source_info=None, lines=flex.split_lines("""\
 ATOM      1  N   GLN A   3      35.299  11.075  19.070  1.00 36.89           N
@@ -7028,6 +7093,8 @@ def exercise_rename_chain_id():
   ph_in = pi.construct_hierarchy()
   ph_in.rename_chain_id(old_id="C", new_id="A")
   assert [c.id.strip() for c in ph_in.chains()] == ["B","A"]
+  ph_in.rename_chain_id(old_id="A", new_id="long_id")
+  assert [c.id.strip() for c in ph_in.chains()] == ["B","long_id"]
 
 def exercise_shift_to_origin():
   pdb_inp = pdb.input(source_info=None, lines="""\
@@ -7166,6 +7233,7 @@ def exercise(args):
     exercise_residue()
     exercise_is_identical_hierarchy()
     exercise_is_similar_hierarchy()
+    # exercise_is_similar_hierarchy_long()
     exercise_atoms()
     exercise_atoms_interleaved_conf()
     exercise_as_pdb_string(
diff --git a/iotbx/regression/tst_hierarchy_copy_select.py b/iotbx/regression/tst_hierarchy_copy_select.py
new file mode 100644
index 0000000000..9efee524be
--- /dev/null
+++ b/iotbx/regression/tst_hierarchy_copy_select.py
@@ -0,0 +1,52 @@
+from __future__ import absolute_import, division, print_function
+import time
+import iotbx.pdb
+
+pdb_str = """
+ATOM      1  N   GLY A   1      -9.009   4.612   6.102  1.00 16.77           N
+ATOM      2  CA  GLY A   1      -9.052   4.207   4.651  1.00 16.57           C
+ATOM      3  C   GLY A   1      -8.015   3.140   4.419  1.00 16.16           C
+ATOM      4  O   GLY A   1      -7.523   2.521   5.381  1.00 16.78           O
+ATOM      5  N   ASN A   2      -7.656   2.923   3.155  1.00 15.02           N
+ATOM      6  CA  ASN A   2      -6.522   2.038   2.831  1.00 14.10           C
+ATOM      7  C   ASN A   2      -5.241   2.537   3.427  1.00 13.13           C
+ATOM      8  O   ASN A   2      -4.978   3.742   3.426  1.00 11.91           O
+ATOM      9  CB  ASN A   2      -6.346   1.881   1.341  1.00 15.38           C
+ATOM     10  CG  ASN A   2      -7.584   1.342   0.692  1.00 14.08           C
+ATOM     11  OD1 ASN A   2      -8.025   0.227   1.016  1.00 17.46           O
+ATOM     12  ND2 ASN A   2      -8.204   2.155  -0.169  1.00 11.72           N
+ATOM     13  N   ASN A   3      -4.438   1.590   3.905  1.00 12.26           N
+ATOM     14  CA  ASN A   3      -3.193   1.904   4.589  1.00 11.74           C
+ATOM     15  C   ASN A   3      -1.955   1.332   3.895  1.00 11.10           C
+ATOM     16  O   ASN A   3      -1.872   0.119   3.648  1.00 10.42           O
+ATOM     17  CB  ASN A   3      -3.259   1.378   6.042  1.00 12.15           C
+ATOM     18  CG  ASN A   3      -2.006   1.739   6.861  1.00 12.82           C
+ATOM     19  OD1 ASN A   3      -1.702   2.925   7.072  1.00 15.05           O
+"""
+
+def test1():
+  """mimic from_ca
+  """
+  inp = iotbx.pdb.input(lines=pdb_str.split("\n"), source_info=None)
+  h = inp.construct_hierarchy()
+  assert h.atoms()[0].id_str() == 'pdb=" N   GLY A   1 "', "'%s'" % h.atoms()[0].id_str()
+  asc=h.atom_selection_cache()
+  sel = asc.selection(string = "protein")
+  h1 = h.deep_copy().select(sel) # deep copy required for failure
+  oc_h = h.overall_counts()
+  oc_h1 = h1.overall_counts()
+  # for a in h1.atoms():
+  #   print (a.id_str())
+  # resname, chain id and resnum are gone!
+  # 'pdb=" N             "'
+  #
+  assert h1.atoms()[0].id_str() == 'pdb=" N   GLY A   1 "', "'%s'" % h1.atoms()[0].id_str()
+  assert oc_h.n_duplicate_atom_labels == 0
+  # Overall counts are also broken since atom labels are mostly the same
+  assert oc_h1.n_duplicate_atom_labels == 0, oc_h1.n_duplicate_atom_labels # This is the problem
+  # exactly labels are not in h1 anymore
+
+if (__name__ == "__main__"):
+  t0 = time.time()
+  test1()
+  print("OK. Time: %8.3f"%(time.time()-t0))
diff --git a/iotbx/regression/tst_hierarchy_long_chain_ids_1.py b/iotbx/regression/tst_hierarchy_long_chain_ids_1.py
index 98ef5a19e7..e3c45a7033 100644
--- a/iotbx/regression/tst_hierarchy_long_chain_ids_1.py
+++ b/iotbx/regression/tst_hierarchy_long_chain_ids_1.py
@@ -248,12 +248,22 @@ def test1():
     """)
 
 def test2():
+  """ Overall counts duplicate atom labels due to incorrect long chain processing
+  """
   inp = iotbx.pdb.input(lines=mmcif_str.split("\n"), source_info=None)
   h = inp.construct_hierarchy()
-  for c in h.chains():
-    for rg in c.residue_groups():
-      for ag in rg.atom_groups():
-        print(c.id, rg.resseq, ag.resname)
+  print("here")
+  oc = h.overall_counts()
+  print("here1")
+  oc.show()
+  print("here2")
+  print(oc.errors())
+  print("here3")
+  # h = inp.construct_hierarchy()
+  # for c in h.chains():
+  #   for rg in c.residue_groups():
+  #     for ag in rg.atom_groups():
+  #       print(c.id, rg.resseq, ag.resname)
 
 if (__name__ == "__main__"):
   t0 = time.time()
diff --git a/iotbx/regression/tst_hierarchy_long_resname_1.py b/iotbx/regression/tst_hierarchy_long_resname_1.py
new file mode 100644
index 0000000000..5b94dfc648
--- /dev/null
+++ b/iotbx/regression/tst_hierarchy_long_resname_1.py
@@ -0,0 +1,121 @@
+from __future__ import absolute_import, division, print_function
+import time
+import iotbx.pdb
+from libtbx.test_utils import assert_lines_in_text
+
+# ------------------------------------------------------------------------------
+
+# from https://github.com/wwPDB/extended-wwPDB-identifier-examples
+# https://github.com/wwPDB/extended-wwPDB-identifier-examples/blob/main/Models/7fgz-extended_CCD_code-model.cif
+mmcif_str = '''
+data_XXXX
+#
+loop_
+_atom_site.group_PDB
+_atom_site.id
+_atom_site.type_symbol
+_atom_site.label_atom_id
+_atom_site.label_alt_id
+_atom_site.label_comp_id
+_atom_site.label_asym_id
+_atom_site.label_entity_id
+_atom_site.label_seq_id
+_atom_site.pdbx_PDB_ins_code
+_atom_site.Cartn_x
+_atom_site.Cartn_y
+_atom_site.Cartn_z
+_atom_site.occupancy
+_atom_site.B_iso_or_equiv
+_atom_site.pdbx_formal_charge
+_atom_site.auth_seq_id
+_atom_site.auth_comp_id
+_atom_site.auth_asym_id
+_atom_site.auth_atom_id
+_atom_site.pdbx_PDB_model_num
+ATOM   2140 N  N   . LYS A 1 261 ? 0.399   -10.171 39.802 1.00 40.11 ? 279 LYS A N   1
+ATOM   2141 C  CA  . LYS A 1 261 ? -0.169  -9.988  41.173 1.00 43.86 ? 279 LYS A CA  1
+ATOM   2142 C  C   . LYS A 1 261 ? 0.687   -9.011  41.991 1.00 41.94 ? 279 LYS A C   1
+ATOM   2143 O  O   . LYS A 1 261 ? 1.044   -7.920  41.556 1.00 39.32 ? 279 LYS A O   1
+ATOM   2144 C  CB  . LYS A 1 261 ? -0.260  -11.336 41.902 1.00 46.47 ? 279 LYS A CB  1
+ATOM   2145 C  CG  . LYS A 1 261 ? -1.583  -12.074 41.713 1.00 49.13 ? 279 LYS A CG  1
+ATOM   2146 C  CD  . LYS A 1 261 ? -1.611  -13.468 42.315 1.00 51.03 ? 279 LYS A CD  1
+ATOM   2147 C  CE  . LYS A 1 261 ? -2.923  -13.799 42.993 1.00 52.86 ? 279 LYS A CE  1
+ATOM   2148 N  NZ  . LYS A 1 261 ? -3.209  -12.856 44.100 1.00 54.19 ? 279 LYS A NZ  1
+HETATM 2149 CA CA  . CA  B 2 .   ? -17.362 -22.385 28.047 1.00 15.20 ? 301 CA  A CA  1
+HETATM 2150 N  N1  . 7ZTVU C 3 .   ? -7.743  -6.355  8.243  1.00 18.72 ? 302 7ZTVU A N1  1
+HETATM 2151 C  C2  . 7ZTVU C 3 .   ? -8.462  -5.534  9.265  1.00 16.68 ? 302 7ZTVU A C2  1
+HETATM 2152 C  C3  . 7ZTVU C 3 .   ? -8.092  -5.865  10.711 1.00 17.35 ? 302 7ZTVU A C3  1
+HETATM 2153 N  N4  . 7ZTVU C 3 .   ? -7.975  -7.334  10.767 1.00 17.04 ? 302 7ZTVU A N4  1
+HETATM 2154 C  C5  . 7ZTVU C 3 .   ? -6.781  -7.689  10.027 1.00 17.11 ? 302 7ZTVU A C5  1
+HETATM 2155 C  C6  . 7ZTVU C 3 .   ? -7.363  -7.730  8.633  1.00 16.97 ? 302 7ZTVU A C6  1
+HETATM 2156 C  C7  . 7ZTVU C 3 .   ? -8.071  -7.951  12.064 1.00 16.77 ? 302 7ZTVU A C7  1
+HETATM 2157 C  C8  . 7ZTVU C 3 .   ? -8.490  -9.440  12.025 1.00 16.60 ? 302 7ZTVU A C8  1
+HETATM 2158 O  O8  . 7ZTVU C 3 .   ? -8.391  -10.140 10.781 1.00 14.14 ? 302 7ZTVU A O8  1
+HETATM 2159 C  C9  . 7ZTVU C 3 .   ? -8.438  -6.311  6.935  1.00 20.13 ? 302 7ZTVU A C9  1
+HETATM 2160 C  C10 . 7ZTVU C 3 .   ? -7.646  -6.965  5.796  1.00 22.62 ? 302 7ZTVU A C10 1
+HETATM 2161 S  S   . 7ZTVU C 3 .   ? -8.372  -6.708  4.279  1.00 24.90 ? 302 7ZTVU A S   1
+HETATM 2162 O  O1S . 7ZTVU C 3 .   ? -7.579  -7.371  3.224  1.00 25.47 ? 302 7ZTVU A O1S 1
+HETATM 2163 O  O2S . 7ZTVU C 3 .   ? -8.533  -5.211  4.078  1.00 25.01 ? 302 7ZTVU A O2S 1
+HETATM 2164 O  O3S . 7ZTVU C 3 .   ? -9.691  -7.326  4.282  1.00 26.42 ? 302 7ZTVU A O3S 1
+HETATM 2165 N  N1  . 7ZTVU D 3 .   ? -11.233 1.732   10.446 1.00 74.61 ? 303 7ZTVU A N1  1
+HETATM 2166 C  C2  . 7ZTVU D 3 .   ? -11.682 2.495   9.264  1.00 77.03 ? 303 7ZTVU A C2  1
+HETATM 2167 C  C3  . 7ZTVU D 3 .   ? -11.907 3.979   9.617  1.00 76.82 ? 303 7ZTVU A C3  1
+HETATM 2168 N  N4  . 7ZTVU D 3 .   ? -12.335 4.282   11.022 1.00 73.77 ? 303 7ZTVU A N4  1
+HETATM 2169 C  C5  . 7ZTVU D 3 .   ? -12.654 3.123   11.908 1.00 69.90 ? 303 7ZTVU A C5  1
+HETATM 2170 C  C6  . 7ZTVU D 3 .   ? -12.333 1.719   11.409 1.00 69.07 ? 303 7ZTVU A C6  1
+HETATM 2171 C  C7  . 7ZTVU D 3 .   ? -11.415 5.232   11.725 1.00 71.46 ? 303 7ZTVU A C7  1
+HETATM 2172 C  C8  . 7ZTVU D 3 .   ? -12.126 5.976   12.871 1.00 70.74 ? 303 7ZTVU A C8  1
+HETATM 2173 O  O8  . 7ZTVU D 3 .   ? -11.360 5.921   14.096 1.00 61.89 ? 303 7ZTVU A O8  1
+HETATM 2174 C  C9  . 7ZTVU D 3 .   ? -10.756 0.357   10.142 1.00 74.40 ? 303 7ZTVU A C9  1
+HETATM 2175 C  C10 . 7ZTVU D 3 .   ? -10.121 -0.289  11.396 1.00 75.92 ? 303 7ZTVU A C10 1
+HETATM 2176 S  S   . 7ZTVU D 3 .   ? -8.757  -1.214  11.091 1.00 77.53 ? 303 7ZTVU A S   1
+HETATM 2177 O  O1S . 7ZTVU D 3 .   ? -7.735  -0.305  10.525 1.00 74.01 ? 303 7ZTVU A O1S 1
+HETATM 2178 O  O2S . 7ZTVU D 3 .   ? -9.042  -2.313  10.141 1.00 83.78 ? 303 7ZTVU A O2S 1
+HETATM 2179 O  O3S . 7ZTVU D 3 .   ? -8.264  -1.864  12.321 1.00 70.44 ? 303 7ZTVU A O3S 1
+HETATM 2180 O  O   . HOH E 4 .   ? 10.567  -6.539  28.064 1.00 30.11 ? 401 HOH A O   1
+HETATM 2181 O  O   . HOH E 4 .   ? 2.856   -8.848  40.692 1.00 32.14 ? 402 HOH A O   1
+HETATM 2182 O  O   . HOH E 4 .   ? -10.577 -29.762 22.833 1.00 22.07 ? 403 HOH A O   1
+HETATM 2183 O  O   . HOH E 4 .   ? 14.705  -5.437  13.988 1.00 28.06 ? 404 HOH A O   1
+'''
+
+def test1():
+  """
+  Test correct reading of long residue name mmCIF file into hierarchy.
+  Testing output to pdb and cif.
+  """
+  inp = iotbx.pdb.input(lines=mmcif_str.split("\n"), source_info=None)
+  h = inp.construct_hierarchy()
+  chains = list(h.overall_counts().chain_ids.keys())
+  chains.sort()
+  answer = ['A']
+  assert (chains == answer), '%s %s' % (chains, answer)
+  # print(h.overall_counts().resnames)
+  resnames = sorted(h.overall_counts().resnames.keys())
+  assert resnames == [' CA', '7ZTVU', 'HOH', 'LYS'], resnames
+
+  o_pdb_str = h.as_pdb_string()
+  # print(o_pdb_str)
+  # Note here incorrect/trimmed residue names
+  # There's no way to correctly output resnames longer than 3 char in PDB format
+  assert_lines_in_text(o_pdb_str, """\
+ATOM      8  CE  LYS A 279      -2.923 -13.799  42.993  1.00 52.86           C
+ATOM      9  NZ  LYS A 279      -3.209 -12.856  44.100  1.00 54.19           N
+HETATM   10 CA    CA A 301     -17.362 -22.385  28.047  1.00 15.20          CA
+HETATM   11  C10 7ZT A 302      -7.646  -6.965   5.796  1.00 22.62           C
+HETATM   12  C2  7ZT A 302      -8.462  -5.534   9.265  1.00 16.68           C
+    """)
+
+  o_cif_str = "%s" % h.as_cif_block()
+  # print(o_cif_str)
+  assert_lines_in_text(o_cif_str, """\
+ATOM     8  CE   .  LYS    A  279  ?   -2.92300  -13.79900  42.99300  1.000  52.86000  C   ?  A  ?  1  1
+ATOM     9  NZ   .  LYS    A  279  ?   -3.20900  -12.85600  44.10000  1.000  54.19000  N   ?  A  ?  1  1
+HETATM  10  CA   .  CA     A  301  ?  -17.36200  -22.38500  28.04700  1.000  15.20000  CA  ?  B  ?  .  1
+HETATM  11  C10  .  7ZTVU  A  302  ?   -7.64600   -6.96500   5.79600  1.000  22.62000  C   ?  C  ?  .  1
+HETATM  12  C2   .  7ZTVU  A  302  ?   -8.46200   -5.53400   9.26500  1.000  16.68000  C   ?  C  ?  .  1
+    """)
+
+if (__name__ == "__main__"):
+  t0 = time.time()
+  test1()
+  print("OK. Time: %8.3f"%(time.time()-t0))
diff --git a/iotbx/regression/tst_merging_statistics.py b/iotbx/regression/tst_merging_statistics.py
index 236243ff71..c4d064a083 100644
--- a/iotbx/regression/tst_merging_statistics.py
+++ b/iotbx/regression/tst_merging_statistics.py
@@ -38,6 +38,7 @@ def exercise(debug=False):
   assert "_reflns_shell" in cif_block
   assert approx_equal(float(cif_block["_reflns.pdbx_Rpim_I_all"]), result.overall.r_pim)
   assert approx_equal(float(cif_block["_reflns.pdbx_CC_half"]), result.overall.cc_one_half)
+  assert approx_equal(float(cif_block["_reflns.percent_possible_obs"]), result.overall.completeness * 100.0)
   assert approx_equal(
     flex.int(cif_block["_reflns_shell.number_measured_obs"]),
     [15737, 15728, 15668, 15371, 14996, 14771, 13899, 13549, 13206, 12528])
diff --git a/iotbx/regression/tst_reflection_file_reader.py b/iotbx/regression/tst_reflection_file_reader.py
index b589e5e9fd..3167240e7a 100644
--- a/iotbx/regression/tst_reflection_file_reader.py
+++ b/iotbx/regression/tst_reflection_file_reader.py
@@ -65,7 +65,7 @@ def exercise_sigma_filtering():
   hkl_in = reflection_file_reader.any_reflection_file(
     "tst_iotbx_hkl_reader.hkl")
   ma = hkl_in.as_miller_arrays(merge_equivalents=False)
-  assert (ma[0].size() == 6)
+  assert (ma[0].size() == 3)
   ma = hkl_in.as_miller_arrays(merge_equivalents=False,
     enforce_positive_sigmas=True)
   assert (ma[0].size() == 3)
diff --git a/iotbx/xds/read_ascii.py b/iotbx/xds/read_ascii.py
index d700c98362..7d3830ccbf 100644
--- a/iotbx/xds/read_ascii.py
+++ b/iotbx/xds/read_ascii.py
@@ -77,10 +77,13 @@ def __init__(self, file_handle, header_only=False, allow_unmerged=True):
           break
         data = line.split()
         assert len(data) == self.number_of_items_in_each_data_record
+        sigma_value = float(data[self.sigma_iobs_column])
+        if sigma_value < 0:
+          continue # negative sigma values should be ignored (outliers)
         h = [int(data[i]) for i in self.miller_index_columns]
         self.miller_indices.append(h)
         self.iobs.append(float(data[self.iobs_column]))
-        self.sigma_iobs.append(float(data[self.sigma_iobs_column]))
+        self.sigma_iobs.append(sigma_value)
         if (self.zd_column is not None):
           self.zd.append(float(data[self.zd_column]))
 
diff --git a/libtbx/env_config.py b/libtbx/env_config.py
index 408407ae99..4b2693cac1 100644
--- a/libtbx/env_config.py
+++ b/libtbx/env_config.py
@@ -582,11 +582,14 @@ def under_root(self, path, return_relocatable_path=False):
     return abs(self.build_path / '..' / path)
 
   def under_base(self, path, return_relocatable_path=False):
-    if self.build_options.use_conda:
-      return abs(self.build_path / '..' / 'conda_base' / path)
-    else:
-      return abs(self.build_path / '..' / 'base' / path)
-    # return abs(os.path.join(abs(self.build_path), '..', 'base', path))
+    result = self.build_path / '..' / 'conda_base' / path
+    if self.installed:
+      result = self.bin_path.anchor / path
+    elif not self.build_options.use_conda:
+      result = self.build_path / '..' / 'base' / path
+    if not return_relocatable_path:
+      result = abs(result)
+    return result
 
   def under_build(self, path, return_relocatable_path=False):
     result = self.as_relocatable_path(path)
diff --git a/mmtbx/model/model.py b/mmtbx/model/model.py
index eb168207e8..59268799d2 100644
--- a/mmtbx/model/model.py
+++ b/mmtbx/model/model.py
@@ -3341,7 +3341,7 @@ def restraints_manager_energies_sites(self,
     )
     return result
 
-  def solvent_selection(self, include_ions=False):
+  def solvent_selection(self, include_ions=False, offset=None):
     result = flex.bool()
     get_class = iotbx.pdb.common_residue_names_get_class
     for a in self.get_hierarchy().atoms():
@@ -3351,6 +3351,9 @@ def solvent_selection(self, include_ions=False):
       elif (a.segid.strip() == "ION") and (include_ions):
         result.append(True)
       else: result.append(False)
+    if(offset is not None):
+      for i in offset:
+        result[i] = False
     return result
 
   def xray_structure_macromolecule(self):
diff --git a/mmtbx/probe/AtomTypes.py b/mmtbx/probe/AtomTypes.py
index c761939032..c5f9b8caa9 100644
--- a/mmtbx/probe/AtomTypes.py
+++ b/mmtbx/probe/AtomTypes.py
@@ -66,7 +66,7 @@ def IsSpecialAminoAcidCarbonyl(resName, atomName):
 # a list of residue names with trailing spaces trimmed.  The second is a list of atoms that qualify
 # for all of the entries in the residue names.  In both cases, the strings are stripped of all
 # spaces to the left and right.
-_AromaticTable = [
+_AromaticAcceptorTable = [
   # Note: Some atoms from these residues are listed in other sections.  The combination of
   # reside and atom name is not duplicated, but there are multiple entries for some residues --
   # this is not a set.
@@ -122,14 +122,15 @@ def IsSpecialAminoAcidCarbonyl(resName, atomName):
 ]
 
 # Is a carbon or nitrogen or hydrogen atom part of an aromatic ring?
-def IsAromatic(resName, atomName):
-  """Given a residue and atom name, determine whether that atom is part of an aromatic ring.
+def IsAromaticAcceptor(resName, atomName):
+  """Given a residue and atom name, determine whether that atom is and acceptor because it is
+  part of an aromatic ring.
   :param resName: String containing the 1-3-character residue name in all caps, including leading space.
   :param atomName: String containing the 1-4-character atom name in all caps, including leading space.
   :returns True if the atom is aromatic in a standard residue, False if not.  Does not handle HET atoms.
   """
 
-  for e in _AromaticTable:
+  for e in _AromaticAcceptorTable:
     if resName.strip() in e[0] and atomName.strip() in e[1]:
       return True
   return False
@@ -138,14 +139,14 @@ def IsAromatic(resName, atomName):
 def Test(inFileName = None):
 
   #========================================================================
-  # Check IsAromatic() to ensure it gives results when expected and not when not.
+  # Check IsAromaticAcceptor() to ensure it gives results when expected and not when not.
   aromaticChecks = [
       ['PHE', 'CE2', True],
       ['  U', 'HN3', True],
       ['ASN',   'O', False]
     ]
   for a in aromaticChecks:
-    assert IsAromatic(a[0],a[1]) == a[2], "AtomTypes.Test(): {} {} not marked as aromatic {}".format(a[0],a[1],a[2])
+    assert IsAromaticAcceptor(a[0],a[1]) == a[2], "AtomTypes.Test(): {} {} not marked as aromatic {}".format(a[0],a[1],a[2])
 
 if __name__ == '__main__':
 
diff --git a/mmtbx/probe/Helpers.py b/mmtbx/probe/Helpers.py
index 0860fe6b0f..7049cd0f67 100644
--- a/mmtbx/probe/Helpers.py
+++ b/mmtbx/probe/Helpers.py
@@ -95,7 +95,7 @@
   gap_weight = 0.25
     .type = float
     .short_caption = Gap weight
-    .help = Weight applied to gap score (-gapweight in probe)
+    .help = Weight applied to gap score (-gapweight in probe). Should be no smaller than probe_radius.
 
   allow_weak_hydrogen_bonds = False
     .type = bool
@@ -146,15 +146,17 @@ def createDotSphereCache(probePhil):
 def createDotScorer(extraAtomInfo, probePhil):
   """
     Helper function to create a DotScorer object, passing it all required Phil parameters
-    in addition to the other constructor parameters.
+    in addition to the other constructor parameters. Also enforces the constraint that the
+    gap_weight is at least as large as the probe_radius.
     :param extraAtomInfo: Can be obtained from getExtraAtomInfo().  Holds extra atom information
     needed for scoring atoms.
-    :param probePhi: Subobject of PHIL parameters for probe.  Can be obtained using
+    :param probePhil: Subobject of PHIL parameters for probe.  Can be obtained using
     self.params.probe from a Program Template program that includes the probe_phil_parameters
     from above in its master PHIL parameters string.
     :returns A mmtbx_probe_ext.DotScorer object.
   """
-  return probeExt.DotScorer(extraAtomInfo, probePhil.gap_weight,
+  gap_weight = max(probePhil.gap_weight, probePhil.probe_radius)
+  return probeExt.DotScorer(extraAtomInfo, gap_weight,
         probePhil.bump_weight, probePhil.hydrogen_bond_weight,
         probePhil.uncharged_hydrogen_cutoff, probePhil.charged_hydrogen_cutoff,
         probePhil.clash_cutoff, probePhil.worse_clash_cutoff,
@@ -408,7 +410,7 @@ def getExtraAtomInfo(model, bondedNeighborLists, useNeutronDistances = False, pr
                 # as an acceptor, perhaps making it a cylinder or a sphere in the center
                 # of the ring.
                 if a.element in ['C','N']:
-                  if AtomTypes.IsAromatic(ag.resname, a.name):
+                  if AtomTypes.IsAromaticAcceptor(ag.resname, a.name):
                     extra.isAcceptor = True
                     warnings += "Marking "+a.name.strip()+" as an aromatic-ring acceptor\n"
 
@@ -680,11 +682,11 @@ def __init__(self, atom, overlap):
         # and if we find one that is on the same aromatic ring then we either ignore this new
         # atom (if it is further) or replace the existing one (if it is closer).
         skip = False
-        if AtomTypes.IsAromatic(a.parent().resname.strip().upper(), a.name.strip().upper()):
+        if AtomTypes.IsAromaticAcceptor(a.parent().resname.strip().upper(), a.name.strip().upper()):
           for c in candidates:
             # See if we belong to the same atom group and are both ring acceptors.  If so, we need to replace
             # or else squash this atom.
-            if (AtomTypes.IsAromatic(c._atom.parent().resname.strip().upper(), c._atom.name.strip().upper()) and
+            if (AtomTypes.IsAromaticAcceptor(c._atom.parent().resname.strip().upper(), c._atom.name.strip().upper()) and
                 a.parent() == c._atom.parent()):
               if overlap < c._overlap:
                 # Replace the further atom with this atom.
diff --git a/mmtbx/probe/README.md b/mmtbx/probe/README.md
index f3051bea85..3001183e18 100644
--- a/mmtbx/probe/README.md
+++ b/mmtbx/probe/README.md
@@ -46,8 +46,8 @@ As mentioned above, the C++ classes can be imported using the **mmtbx_probe_ext*
 There are also a set of Python classes that provide support for applications that want to use the C++ classes.
 
 * **AtomTypes.py:** This contains classes and functions that help determine the characteristics of atoms, which is needed
-to know how to handle them during Probe2 calculations.  The **IsAromatic()** function tells whether
-a specified atom from a specified residue is part of an aromatic ring in a standard residue.
+to know how to handle them during Probe2 calculations.  The **IsAromaticAcceptor()** function tells whether
+a specified atom from a specified residue is an acceptor because it is part of an aromatic ring in a standard residue.
 
 * **Helpers.py.** This contains helper functions needed by both Probe2 and Reduce2.  See the file itself for a complete
 list of functions and parameters.  Notable ones include:
diff --git a/mmtbx/probe/Scoring.cpp b/mmtbx/probe/Scoring.cpp
index 4f002c12b8..9ba504d509 100644
--- a/mmtbx/probe/Scoring.cpp
+++ b/mmtbx/probe/Scoring.cpp
@@ -105,6 +105,19 @@ DotScorer::CheckDotResult DotScorer::check_dot(
   Point absoluteDotLocation = sourceAtom.data->xyz + dotOffset;
   Point probLoc = sourceAtom.data->xyz + dotOffset.normalize() * (dotOffset.length() + probeRadius);
 
+  // Check to see if the dot should be removed from consideration because it is inside an excluded atom.
+  // Doing this test ahead of the neighbor-interaction test makes things faster for the long-running
+  // 4fen test case.
+  for (scitbx::af::shared<iotbx::pdb::hierarchy::atom>::const_iterator e = exclude.begin();
+    e != exclude.end(); e++) {
+    // The original Probe code does internal checks for Phantom Hydrogens, but we handle that
+    // in the calling routine by properly adjusting the list of atoms to be excluded.
+    double vdwe = m_extraInfoMap.getMappingFor(*e).getVdwRadius();
+    if ((absoluteDotLocation - e->data->xyz).length_sq() < vdwe * vdwe) {
+      return ret;
+    }
+  }
+
   bool isHydrogenBond = false;          ///< Are we looking at a hydrogen bond to our neighbor?
   bool tooCloseHydrogenBond = false;    ///< Are we too close to be a hydrogen bond?
   double hydrogenBondMinDist = 0;       ///< Hydrogen bond minimum distance based on the atom types (will be set below).
@@ -175,20 +188,6 @@ DotScorer::CheckDotResult DotScorer::check_dot(
     }
   }
 
-  // Check to see if the dot should be removed from consideration because it is also inside an excluded atom.
-  if (keepDot) {
-    for (scitbx::af::shared<iotbx::pdb::hierarchy::atom>::const_iterator e = exclude.begin();
-         e != exclude.end(); e++) {
-      // The original Probe code does internal checks for Phantom Hydrogens, but we handle that
-      // in the calling routine by properly adjusting the list of atoms to be excluded.
-      double vdwe = m_extraInfoMap.getMappingFor(*e).getVdwRadius();
-      if ((absoluteDotLocation - e->data->xyz).length_sq() < vdwe * vdwe) {
-        keepDot = false;
-        break;
-      }
-    }
-  }
-
   // If this dot is a keeper, fill in non-default return values.
   if (keepDot) {
     // Determine the overlap type and amount of overlap.
diff --git a/mmtbx/programs/reduce2.py b/mmtbx/programs/reduce2.py
index 2d4932cfeb..c2b93c1d24 100644
--- a/mmtbx/programs/reduce2.py
+++ b/mmtbx/programs/reduce2.py
@@ -1,5 +1,5 @@
 ##################################################################################
-# Copyright(c) 2021, Richardson Lab at Duke
+# Copyright(c) 2021-2023, Richardson Lab at Duke
 # Licensed under the Apache 2 license
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
@@ -34,14 +34,20 @@
 from mmtbx.programs import probe2
 import copy
 import tempfile
+from iotbx.data_manager import DataManager
+import csv
 
-version = "0.6.0"
+version = "1.0.0"
 
 master_phil_str = '''
 approach = *add remove
   .type = choice
   .short_caption = Add or remove Hydrogens
   .help = Determines whether Reduce will add (and optimize) or remove Hydrogens from the model
+keep_existing_H = False
+  .type = bool
+  .short_caption = Do not remove Hydrogens in the original model
+  .help = Keep existing H atoms in the model
 n_terminal_charge = *residue_one first_in_chain no_charge
   .type = choice(multi=False)
   .short_caption = N terminal charge approach
@@ -82,6 +88,22 @@
   .type = bool
   .short_caption = Profile the entire run
   .help = Profile the performance of the entire run
+comparison_file = None
+  .type = str
+  .short_caption = Compare the Mover scores from this run with those in comparison_file
+  .help = Points to a comparison_file that is the result of running Hydrogenate or Reduce or Reduce2 or some other hydrogen-placement program. The Probe2 scores for the Movers found in the current run are compared against the scores for comparison_file and stored in a file with the same name as output.file_name with _comparison.csv appended. If None, no comparison is done.
+verbosity = 3
+  .type = int
+  .short_caption = Level of detail in description file
+  .help = Level of detail in description file.
+bonded_neighbor_depth = 4
+  .type = int
+  .short_caption = How many neighbors to consider bonded (>3 accepted only if hydrogen)
+  .help = When looking for interactions between atoms, this specifies how many hops we should take when looking for atoms that are excluded because we share a common chain of bonds. Lengths more than 3 only happen when one end is a hydrogen. This value should not be changed from the default except for regression tests against earlier versions.
+stop_on_any_missing_hydrogen = False
+  .type = bool
+  .short_caption = Emit a Sorry and stop when any hydrogen in the model has insufficient restraints.
+  .help = Emit a Sorry and stop when any hydrogen in the model has insufficient restraints. It will always emit when there are one or more residues without restraints.
 
 output
   .style = menu_item auto_align
@@ -525,11 +547,11 @@ def _AddFlipkinBase(states, views, fileName, fileBaseName, model, alts, bondedNe
 
     # Find out the type of the residue, used to determine the type of flip.
     type = '?'
-    if s.resName == 'ASN':
+    if s.resName[-3:] == 'ASN':
       type = 'N'
-    elif s.resName == 'GLN':
+    elif s.resName[-3:] == 'GLN':
       type = 'Q'
-    elif s.resName == 'HIS':
+    elif s.resName[-3:] == 'HIS':
       type = 'H'
 
     if i > 0:
@@ -865,7 +887,7 @@ def _MakeProbePhil(self, movers_to_check, temp_output_file_name):
     newParams.__inject__('source_selection', source_selection)      # Not default
     newParams.__inject__('target_selection', target_selection)      # Not default
     newParams.approach = 'both'                                     # Not default
-    newParams.__inject__('excluded_bond_chain_length', 4)
+    newParams.__inject__('excluded_bond_chain_length', self._bondedNeighborDepth) # Not default
     newParams.__inject__('minimum_water_hydrogen_occupancy', 0.66)  # Not default
     newParams.__inject__('maximum_water_hydrogen_b', 40.0)          # Not default
     newParams.__inject__('include_mainchain_mainchain', True)
@@ -913,8 +935,8 @@ def _AddHydrogens(self):
       use_neutron_distances=self.params.use_neutron_distances,
       n_terminal_charge=self.params.n_terminal_charge,
       exclude_water=True,
-      stop_for_unknowns=True,
-      keep_existing_H=False
+      stop_for_unknowns=False,
+      keep_existing_H=self.params.keep_existing_H
     )
     reduce_add_h_obj.run()
     reduce_add_h_obj.show(None)
@@ -924,6 +946,12 @@ def _AddHydrogens(self):
       for res in missed_residues:
         bad += " " + res
       raise Sorry("Restraints were not found for the following residues:"+bad)
+    insufficient_restraints = list(reduce_add_h_obj.site_labels_no_para)
+    if self.params.stop_on_any_missing_hydrogen and len(insufficient_restraints) > 0:
+      bad = insufficient_restraints[0]
+      for res in insufficient_restraints[1:]:
+        bad += "," + res
+      raise Sorry("Insufficient restraints were found for the following atoms:"+bad)
     self.model = reduce_add_h_obj.get_model()
 
 # ------------------------------------------------------------------------------
@@ -993,16 +1021,31 @@ def _DescribeLockdown(self, flipMover, invertFlip, fixedUp):
       flipStateString = 'Flipped'
     else:
       flipStateString = 'Unflipped'
+    # Look up the residue that is described and find Nitrogens in it whose names are
+    # longer than one character. These will be part of the flipping part of the residue.
+    # See if any of these atoms residue groups have the same altid as the flipMover; if
+    # so, use that one.  If not, use ''.  If we use '', replace it with '.' so that we
+    # can properly parse the string.
+    # @todo Atom selection here is will not necessarily work for future flip Movers.
     # @todo How to handle cases where the different alternates resulted in different placement
     # for the flip Movers?  Presumably, we want to do them in backwards order and replace all of
     # the alternates with '.' so that they always get handled and the last one is the one that
     # is set.
-    '''
-    altId = flipMover.altId
+    altId = ''
+    flipAlt = flipMover.altId.strip()
+    resAlt = ''
+    for chain in self.model.chains():
+      if flipMover.chain.strip() == chain.id.strip():
+        for rg in chain.residue_groups():
+          if int(flipMover.resId) == rg.resseq_as_int():
+            for a in rg.atoms():
+              if (a.element == 'N') and (len(a.name.strip()) > 1):
+                if a.parent().altloc.lower() == flipAlt.lower():
+                  resAlt = flipAlt
+    if flipAlt == resAlt:
+      altId = flipAlt
     if altId.strip() == '':
       altId = '.'
-    '''
-    altId = '.'
     return '{} {} {} {} {}{} {} {}'.format(flipMover.modelId+1, altId.lower(), flipMover.chain,
       flipMover.resName, flipMover.resId, flipMover.iCode, flipStateString, adjustedString)
 
@@ -1019,7 +1062,7 @@ def validate(self):
         pad = 'H'
       base = os.path.splitext(os.path.basename(inName))[0] + pad
       self.params.output.filename = base + suffix
-      print('Writing model output to', self.params.output.filename)
+      print('Writing model output to', self.params.output.filename, file=self.logger)
 
     self.data_manager.has_models(raise_sorry=True)
     if self.params.output.description_file_name is None:
@@ -1039,6 +1082,9 @@ def validate(self):
 
   def run(self):
 
+    # Set our bonded-neighbor depth
+    self._bondedNeighborDepth = self.params.bonded_neighbor_depth
+
     # String describing the run that will be output to the specified file.
     outString = 'reduce2 v.{}, run {}\n'.format(version, datetime.now().strftime("%Y-%m-%d %H:%M:%S"))
     for a in sys.argv:
@@ -1079,18 +1125,19 @@ def run(self):
       opt = Optimizers.FastOptimizer(self.params.probe, self.params.add_flip_movers,
         self.model, probeRadius=0.25, altID=self.params.alt_id, modelIndex=self.params.model_id,
         preferenceMagnitude=self.params.preference_magnitude,
+        bondedNeighborDepth = self._bondedNeighborDepth,
         nonFlipPreference=self.params.non_flip_preference,
         skipBondFixup=self.params.skip_bond_fix_up,
         flipStates = self.params.set_flip_states,
-        verbosity=3)
+        verbosity=self.params.verbosity)
       doneOpt = work_clock()
       outString += opt.getInfo()
       outString += 'Time to Add Hydrogen = '+str(doneAdd-startAdd)+'\n'
       outString += 'Time to Interpret = '+str(doneInt-startInt)+'\n'
       outString += 'Time to Optimize = '+str(doneOpt-startOpt)+'\n'
       if self.params.output.print_atom_info:
-        print('Atom information used during calculations:')
-        print(opt.getAtomDump())
+        print('Atom information used during calculations:', file=self.logger)
+        print(opt.getAtomDump(), file=self.logger)
 
     else: # Removing Hydrogens from the model rather than adding them.
       make_sub_header('Removing Hydrogens', out=self.logger)
@@ -1122,8 +1169,117 @@ def run(self):
     print('Wrote', self.params.output.filename,'and',
       self.params.output.description_file_name, file = self.logger)
 
+    # If we've been asked to do a comparison with another program's output, do it.
+    if self.params.comparison_file is not None:
+      make_sub_header('Comparing with other model', out=self.logger)
+
+      # Construct the file names we'll be using.
+      compareOutName = self.params.output.filename + "_comparison.csv"
+
+      # Find the list of all Movers in the model, which will be used to select
+      # each in turn for evaluation.
+      moverLocations = _FindMoversInOutputString(outString)
+
+      # Make the first line in our table be the header line
+      table = [ ["Mover",
+                 "Score from Reduce2",
+                 "Score from {}".format(self.params.comparison_file),
+                 "Difference",
+                 "Star if other is higher"] ]
+
+      # Run Probe2 on each of the Movers for our current model and for
+      # comparison_file and determine the summary score of the Mover against
+      # the rest of the model in each case. Accumulate these into a table.
+      for i, m in enumerate(moverLocations):
+        print ('Computing scores for Mover',i,'of',len(moverLocations))
+
+        #============================================================
+        # Find the score for our output.
+
+        # Make the Probe2 Phil parameters, then overwrite the ones that were
+        # filled in with values that we want for our summaries.
+        source = [ m ]
+        tempName = tempfile.mktemp()
+        newParams = self._MakeProbePhil(source, tempName)
+        newParams.approach = 'once'       # Reduce2/Reduce only check Mover atom to other
+        newParams.output.format = 'raw'
+        newParams.output.contact_summary = True
+        newParams.output.condensed = True
+        newParams.output.count_dots = True
+
+        # Run Probe2
+        p2 = probe2.Program(self.data_manager, newParams, master_phil=probe2.master_phil_str,
+          logger=self.logger)
+        p2.overrideModel(self.model)
+        dots, output = p2.run()
+
+        # Parse the output to find the scores for each of the two directions
+        # and sum them up. There is one per line and the number ends with #.
+        lines = output.strip().split('\n')
+        values = []
+        for line in lines:
+          stripped_line = line.strip()
+          if stripped_line:
+            number = float(stripped_line.split('#')[0].strip())
+            values.append(number)
+        myScore = sum(values)
+
+        print('My values for Mover', str(m), 'are', values, 'sum is', myScore, file=self.logger)
+        os.unlink(tempName)
+
+        #============================================================
+        # Find the score for the comparison file.
+
+        # Read the file using a new datamanager
+        dm = DataManager()
+        dm.process_model_file(self.params.comparison_file)
+        #otherModel = dm.get_model(self.params.comparison_file)
+
+        # Make the Probe2 Phil parameters, then overwrite the ones that were
+        # filled in with values that we want for our summaries.
+        tempName = tempfile.mktemp()
+        newParams = self._MakeProbePhil(source, tempName)
+        newParams.approach = 'once'       # Reduce2/Reduce only check Mover atom to other
+        newParams.output.format = 'raw'
+        newParams.output.contact_summary = True
+        newParams.output.condensed = True
+        newParams.output.count_dots = True
+
+        # Run Probe2
+        p2 = probe2.Program(dm, newParams, master_phil=probe2.master_phil_str,
+          logger=self.logger)
+        p2.overrideModel(dm.get_model(self.params.comparison_file))
+        dots, output = p2.run()
+
+        # Parse the output to find the scores for each of the two directions
+        # and sum them up. There is one per line and the number ends with #.
+        lines = output.strip().split('\n')
+        values = []
+        for line in lines:
+          stripped_line = line.strip()
+          if stripped_line:
+            number = float(stripped_line.split('#')[0].strip())
+            values.append(number)
+        otherScore = sum(values)
+
+        print('Other values for Mover', str(m), 'are', values, 'sum is', otherScore, file=self.logger)
+        os.unlink(tempName)
+
+        #============================================================
+        # Add the line to the table, indicating if the other is better.
+        mark = ''
+        if otherScore > myScore:
+          mark = '*'
+        table.append( [str(m), myScore, otherScore, myScore - otherScore, mark] )
+
+      # Write the table to our output CSV file
+      with open(compareOutName, 'w', newline='') as csvfile:
+          writer = csv.writer(csvfile)
+          writer.writerows(table)
+
     # If we've been asked to make Flipkins, then make each of them.
     if self.params.add_flip_movers and self.params.output.flipkin_directory is not None:
+      make_sub_header('Producing flipkins', out=self.logger)
 
       # Find the base name of the two output files we will produce.
       inName = self.data_manager.get_default_model_name()
@@ -1226,6 +1382,8 @@ def run(self):
             skipBondFixup=self.params.skip_bond_fix_up,
             flipStates = flipStates,
             verbosity=3)
+          print('Results of optimization:', file=self.logger)
+          print(opt.getInfo(), file=self.logger)
           self._ReinterpretModel()
 
           # Get the other characteristics we need to know about each atom to do our work.
@@ -1293,7 +1451,7 @@ def run(self):
 
         # Write the base information in the Flipkin, not including the moving atoms in
         # the Movers that will be placed, or atoms bonded to the moving atoms.
-        flipkinText = _AddFlipkinBase(amides, views, self.params.output.filename, base, self.model,
+        flipkinText = _AddFlipkinBase(hists, views, self.params.output.filename, base, self.model,
           alts, bondedNeighborLists, moverLocations, inSideChain, inWater, inHet)
 
         # Make two configurations, the one that Reduce picked and the one
@@ -1315,7 +1473,7 @@ def run(self):
           # @todo Remove this reinterpretation once place_hydrogens() does all the interpretation we need.
           self._ReinterpretModel()
 
-          # Run optimization, locking the specified Amides into each configuration.
+          # Run optimization, locking the specified Histidines into each configuration.
           # Don't do fixup on the ones that are locked down.  Make sure that we can
           # avoid adding a comma on the first flip state that is added.
           flipStates = self.params.set_flip_states
@@ -1347,6 +1505,8 @@ def run(self):
             skipBondFixup=self.params.skip_bond_fix_up,
             flipStates = flipStates,
             verbosity=3)
+          print('Results of optimization:', file=self.logger)
+          print(opt.getInfo(), file=self.logger)
           self._ReinterpretModel()
 
           # Get the other characteristics we need to know about each atom to do our work.
@@ -1392,7 +1552,7 @@ def run(self):
 
     # Report profiling info if we've been asked to in the Phil parameters
     if self.params.profile:
-      print('Profile results:')
+      print('Profile results:', file=self.logger)
       import pstats
       profile_params = {'sort_by': 'time', 'num_entries': 20}
       self._pr.disable()
diff --git a/mmtbx/reduce/InteractionGraph.py b/mmtbx/reduce/InteractionGraph.py
index 550f8f5834..19393be7a8 100644
--- a/mmtbx/reduce/InteractionGraph.py
+++ b/mmtbx/reduce/InteractionGraph.py
@@ -198,7 +198,7 @@ def _PairsOverlap(mover1, atoms1, positions1,
   :param mover1: The first Mover
   :param atoms1: Atom list for the first Mover
   :param positions1: probe.PositionReturn.positions holding possible positions for each.
-  :param mover2: The first Mover
+  :param mover2: The second Mover
   :param atoms2: Atom list for the second Mover
   :param positions2: probe.PositionReturn.positions holding possible positions for each.
   :param extraAtomInfoMap: probe.ExtraAtomInfoMap that can be used to look
@@ -211,24 +211,41 @@ def _PairsOverlap(mover1, atoms1, positions1,
   :returns True if a pair of atoms with one from each overlap, False if not.
   """
 
-  for i1, p1 in enumerate(positions1):
+  # Construct look-up tables for the radii of each atom to pull these calculations
+  # outside of the loops.
+  radii1 = []
+  for ai1 in range(len(atoms1)):
+    radii1.append(extraAtomInfoMap.getMappingFor(atoms1[ai1]).vdwRadius)
+  radii2 = []
+  for ai2 in range(len(atoms2)):
+    radii2.append(extraAtomInfoMap.getMappingFor(atoms2[ai2]).vdwRadius)
+
+  ret = False
+  for p1 in positions1:
     for ai1 in range(len(p1)):
-      r1 = extraAtomInfoMap.getMappingFor(atoms1[ai1]).vdwRadius
-      for i2, p2 in enumerate(positions2):
+      r1 = radii1[ai1]
+      x1 = p1[ai1][0]
+      y1 = p1[ai1][1]
+      z1 = p1[ai1][2]
+      limit1 = 2*probeRad + r1
+      for p2 in positions2:
         for ai2 in range(len(p2)):
-          r2 = extraAtomInfoMap.getMappingFor(atoms2[ai2]).vdwRadius
-          dx = p1[ai1][0] - p2[ai2][0]
-          dy = p1[ai1][1] - p2[ai2][1]
-          dz = p1[ai1][2] - p2[ai2][2]
+          r2 = radii2[ai2]
+          dx = x1 - p2[ai2][0]
+          dy = y1 - p2[ai2][1]
+          dz = z1 - p2[ai2][2]
           dSquared = dx*dx + dy*dy + dz*dz
-          limit = r1 + r2 + 2*probeRad
+          limit = limit1 + r2
           limitSquared = limit*limit
           if dSquared <= limitSquared:
             # Add the opposite Mover to each atom; they interact
             atomMoverSets[atoms1[ai1]].add(mover2)
             atomMoverSets[atoms2[ai2]].add(mover1)
-            return True
-  return False
+            # Set the return value to True but do not quit looking
+            # because we need to catalog all of the atomMoverSets
+            # interactions, not just the first one found.
+            ret = True
+  return ret
 
 #######################################################################################################
 # Test code and objects below here
@@ -341,7 +358,7 @@ def Test():
     if maxLen != e[2]:
       return "Expected max sized component of "+str(e[2])+", found "+str(maxLen)+" for case "+str(i)
 
-    # Check atom/Mover overlaps by finding the set of lengths that are present accross all atoms.
+    # Check atom/Mover overlaps by finding the set of lengths that are present across all atoms.
     lengths = set()
     for a in atoms:
       lengths.add(len(am[a]))
diff --git a/mmtbx/reduce/Movers.py b/mmtbx/reduce/Movers.py
index f22f48a58f..843a91ec8e 100644
--- a/mmtbx/reduce/Movers.py
+++ b/mmtbx/reduce/Movers.py
@@ -172,7 +172,7 @@ def FixUp(self, coarseIndex):
     # No fixups for any coarse index.
     return FixUpReturn([], [], [], [])
   def PoseDescription(self, coarseIndex, fineIndex, fixedUp):
-    if coarseIndex >= len(self.CoarsePositions().positions) or (
+    if coarseIndex >= len(self.CoarsePositions().positions) or fineIndex is not None and (
         fineIndex > 0 and fineIndex >= len(self.FinePositions(0).positions)):
       return "Unrecognized state . ."
     else:
@@ -180,7 +180,7 @@ def PoseDescription(self, coarseIndex, fineIndex, fixedUp):
 
 ##################################################################################
 class _MoverRotator(object):
-  def __init__(self, atoms, axis, dihedral, offset, coarseRange, coarseStepDegrees = 30.0,
+  def __init__(self, atoms, axis, dihedral, offset, coarseRange, coarseStepDegrees = 15.0,
                 doFineRotations = True, fineStepDegrees = 1.0,
                 preferenceFunction = None, preferredOrientationScale = 1.0):
     """Constructs a Rotator, to be called by a derived class or test program but
@@ -364,11 +364,14 @@ def FixUp(self, coarseIndex):
     return FixUpReturn([], [], [], [])
 
   def PoseDescription(self, coarseIndex, fineIndex, fixedUp):
-    if coarseIndex >= len(self.CoarsePositions().positions) or (
+    if coarseIndex >= len(self.CoarsePositions().positions) or fineIndex is not None and (
         fineIndex > 0 and fineIndex >= len(self.FinePositions(0).positions)):
       return "Unrecognized state . ."
     else:
-      angle = self._offset + self._coarseAngles[coarseIndex] + self._fineAngles[fineIndex]
+      fineOffset = 0
+      if fineIndex is not None:
+        fineOffset = self._fineAngles[fineIndex]
+      angle = self._offset + self._coarseAngles[coarseIndex] + fineOffset
       while angle > 180: angle -= 360
       while angle < -180: angle += 360
       return "Angle {:.1f} deg .".format(angle)
@@ -426,6 +429,9 @@ def __init__(self, atom, bondedNeighborLists, hParameters, potentialAcceptors =
         friends.append(b)
 
     # Determine the preference function (180 or 120) based on friends bonding structure
+    # NOTE: We replace the preferenceFunction with None below to match original reduce
+    # behavior. We leave the math in here so that if we later decide to have a preference
+    # it will be correctly determined based on the number of neighbors.
     # @todo Consider parameterizing the magic constant of 0.1 for the preference magnitude
     if len(friends) == 2:
       # There are two neighbors, so our function repeats every 180 degrees
@@ -436,6 +442,8 @@ def preferenceFunction(degrees): return 0.1 + 0.1 * math.cos(degrees * (math.pi/
     else:
       raise ValueError("MoverSingleHydrogenRotator(): Atom's bonded neighbor's neighbor does not have 2-3 other bonds "+
       "it has "+str(len(friends)))
+    # Original Reduce did not have a preference, so we remove that here.
+    preferenceFunction = None
 
     # Determine the axis to rotate around, which starts at the partner atom and points at the neighbor.
     normal = (rvec3(neighbor.xyz) - rvec3(partner.xyz)).normalize()
@@ -484,7 +492,7 @@ def preferenceFunction(degrees): return 0.1 + 0.1 * math.cos(degrees * (math.pi/
 
 ##################################################################################
 class MoverNH3Rotator(_MoverRotator):
-  def __init__(self, atom, bondedNeighborLists, hParameters, coarseStepDegrees = 30.0,
+  def __init__(self, atom, bondedNeighborLists, hParameters, coarseStepDegrees = 15.0,
                 fineStepDegrees = 1.0, preferredOrientationScale = 1.0):
     """ A Mover that rotates three Hydrogens around an axis from their bonded Nitrogen neighbor
        to the single bonded partner of its partner.  This is designed for use with NH3+,
@@ -534,13 +542,19 @@ def __init__(self, atom, bondedNeighborLists, hParameters, coarseStepDegrees = 3
     for b in bonded:
       if b.i_seq != neighbor.i_seq:
         friends.append(b)
-    if len(friends) != 3:
-      raise ValueError("MoverNH3Rotator(): Partner does not have three bonded friends")
+    if len(friends) < 3:
+      raise ValueError("MoverNH3Rotator(): Partner does not have at least three bonded friends")
 
-    # Set the preference function to like 120-degree rotations away from the starting location.
-    # @todo Consider parameterizing the magic constant of 0.1 for the preference magnitude
-    # (for example, we might use preferredOrientationScale for this and not set a separate one?)
-    def preferenceFunction(degrees): return 0.1 + 0.1 * math.cos(degrees * (math.pi/180) * (360/120))
+    preferenceFunction = None
+
+    # If we have exactly three neighbors, then we set a preference function that prefers
+    # the staggered orientation. If there are more, we just leave things alone.
+    if len(friends) == 3:
+
+      # Set the preference function to prefer 120-degree rotations away from the starting location.
+      # @todo Consider parameterizing the magic constant of 0.1 for the preference magnitude
+      # (for example, we might use preferredOrientationScale for this and not set a separate one?)
+      def preferenceFunction(degrees): return 0.1 + 0.1 * math.cos(degrees * (math.pi/180) * (360/120))
 
     # Determine the axis to rotate around, which starts at the partner and points at the neighbor.
     normal = (rvec3(neighbor.xyz) - rvec3(partner.xyz)).normalize()
@@ -562,7 +576,7 @@ def preferenceFunction(degrees): return 0.1 + 0.1 * math.cos(degrees * (math.pi/
     atoms.extend(hydrogens)
 
     # Construct our parent class, which will do all of the actual work based on our inputs.
-    _MoverRotator.__init__(self, atoms, axis, dihedral, offset, 180, coarseStepDegrees,
+    _MoverRotator.__init__(self, atoms, axis, dihedral, offset, 60, coarseStepDegrees,
       fineStepDegrees = fineStepDegrees, preferenceFunction = preferenceFunction,
       preferredOrientationScale = preferredOrientationScale)
 
@@ -643,7 +657,7 @@ def __init__(self, atom, bondedNeighborLists, hParameters):
 
 ##################################################################################
 class MoverTetrahedralMethylRotator(_MoverRotator):
-  def __init__(self, atom, bondedNeighborLists, hParameters, coarseStepDegrees = 30.0,
+  def __init__(self, atom, bondedNeighborLists, hParameters, coarseStepDegrees = 15.0,
                   fineStepDegrees = 1.0, preferredOrientationScale = 1.0):
     """ A Mover that rotates three Hydrogens around an axis from their bonded Carbon neighbor
        to the single bonded partner of its partner.  This is designed for use with tetrahedral
@@ -912,7 +926,7 @@ def PoseDescription(self, coarseIndex, fineIndex, fixedUp):
         fString = 'AnglesNotAdjusted'
     else:
       fString = '.'
-    if coarseIndex >= len(self.CoarsePositions().positions) or (
+    if coarseIndex >= len(self.CoarsePositions().positions) or fineIndex is not None and (
         fineIndex > 0 and fineIndex >= len(self.FinePositions(0).positions)):
       return "Unrecognized state ."
     elif coarseIndex == 0:
@@ -1208,7 +1222,7 @@ def FixUp(self, coarseIndex):
 
   def PoseDescription(self, coarseIndex, fineIndex, fixedUp):
     if self._enabledFlipStates == 3:
-      if coarseIndex >= len(self.CoarsePositions().positions) or (
+      if coarseIndex >= len(self.CoarsePositions().positions) or fineIndex is not None and (
           fineIndex > 0 and fineIndex >= len(self.FinePositions(0).positions)):
         return "Unrecognized state . ."
       elif coarseIndex < 4:
@@ -1283,9 +1297,9 @@ def _rotateAroundAxis(atom, axis, degrees):
      :param atom: iotbx.pdb.hierarchy.atom or scitbx::vec3<double> or
      scitbx.matrix.rec(xyz, (3,1)) to be moved.
      :param axis: flex array with two scitbx::vec3<double> points, the first
-     of which is the origin and the second is a vector pointing in the direction
-     of the axis of rotation.  Positive rotations will be right-handed rotations
-     around this axis.
+     of which is the origin in space around which to rotate and the second is
+     a vector pointing from the origin in the direction of the axis of rotation.
+     Positive rotations will be right-handed rotations around this axis.
      :param degrees: How much to rotate the atom around the axis.
      Positive rotation is right-handed around the axis.
      :returns the new location for the atom.
@@ -1297,27 +1311,12 @@ def _rotateAroundAxis(atom, axis, degrees):
   except Exception:
     pos = rvec3(atom)
 
-  # Project the atom position onto the axis, finding its closest point on the axis.
-  # The position lies on a plane whose normal points along the axis vector.  The
-  # point we seek is the intersection of the axis with this plane.
-  # The plane equation will be the normalized axis direction vector and the offset
-  # from the origin such that N.point + d = 0 defines the plane.
-  # Solve for the time at which the line's ray crosses the plane and then solve for
-  # the location along the line at that time.  t = - (d + (lineOrigin * planeNormal)) /
-  # (lineDirection * planeNormal).  Because the line direction and normal are the
-  # same, the divisor is 1.
-  normal = lvec3(axis[1]).normalize()
-  d = -normal*pos
-  t = - (d + (normal * rvec3(axis[0])))
-  nearPoint = lvec3(axis[0]) + t * normal
-
-  # Find the vector from the closest point towards the atom, which is its offset
-  offset = lvec3(pos) - nearPoint
-
-  # Rotate the offset vector around the axis by the specified angle.  Add the new
-  # offset to the closest point. Store this as the new location for this atom and angle.
-  newOffset = offset.rotate_around_origin(lvec3(axis[1]), degrees*math.pi/180)
-  return nearPoint + newOffset
+  # The axis of rotation for this function is specified as the two ends of the axis.
+  # The axis passed in has the point around which to rotate and the direction vector
+  # from the origin, so we need to add those.
+  return lvec3(scitbx.matrix.rotate_point_around_axis(
+      axis_point_1 = axis[0], axis_point_2 = rvec3(axis[0]) + rvec3(axis[1]),
+      point = pos, angle = degrees, deg = True))
 
 def _rotateHingeDock(movableAtoms, hingeIndex, firstDockIndex, secondDockIndex, alphaCarbon):
   '''Perform the three-step rotate-hinge-dock calculation described in
@@ -1567,8 +1566,8 @@ def prefFunc2(x):
     # Test coarseStepDegrees default behavior.
     rot = _MoverRotator(atoms,axis, 0, 0, 180)
     coarse = rot.CoarsePositions()
-    if len(coarse.positions) != 12:
-      return "Movers.Test() _MoverRotator Default coarse step: Expected 12, got "+str(len(coarse.positions))
+    if len(coarse.positions) != 24:
+      return "Movers.Test() _MoverRotator Default coarse step: Expected 24, got "+str(len(coarse.positions))
 
     # Test doFineRotations = False and 180 degree coarseStepDegrees.
     rot = _MoverRotator(atoms,axis, 0, 0, 180, 180, False)
@@ -1582,18 +1581,18 @@ def prefFunc2(x):
     # Test fineStepDegrees setting.
     rot = _MoverRotator(atoms,axis, 0, 0, 180, fineStepDegrees = 2)
     fine = rot.FinePositions(0)
-    # +/- 15 degrees in 1-degree steps, but we don't do the +15 because it will be handled by the next
-    # rotation up.
-    if len(fine.positions) != 14:
-      return "Movers.Test() _MoverRotator setting fine step: Expected 14, got "+str(len(fine.positions))
+    # +/- 7.5 degrees in 2-degree steps, but we wouldn't do the +7.5 because it will be handled by the next
+    # rotation up. So we get +/- 2, 4, and 6
+    if len(fine.positions) != 6:
+      return "Movers.Test() _MoverRotator setting fine step: Expected 6, got "+str(len(fine.positions))
 
     # Test the PoseDescription
-    if rot.PoseDescription(1,1, False) != "Angle -28.0 deg .":
+    if rot.PoseDescription(1,1, False) != "Angle -13.0 deg .":
       return "Movers.Test() _MoverRotator: Unexpected results for PoseDescription, got "+rot.PoseDescription(1,1, False)
 
     # Test setting an offset and counterbalancing dihedral.
     rot = _MoverRotator(atoms,axis, -10, 10, 180, fineStepDegrees = 2)
-    if rot.PoseDescription(1,1, False) != "Angle -18.0 deg .":
+    if rot.PoseDescription(1,1, False) != "Angle -3.0 deg .":
       return "Movers.Test() _MoverRotator: Unexpected results for offset, got "+rot.PoseDescription(1,1, False)
 
     # Verify that the coarse and fine results don't change when the atom position is moved after
@@ -1677,6 +1676,7 @@ def __init__(self, paramN, paramA2):
     if h.xyz[2] != 1 or abs(-h.xyz[0]-math.sqrt(2)) > 1e-5:
       return "Movers.Test() MoverSingleHydrogenRotator pair: bad H placement"
 
+    ''' Preference function has been removed from this class
     # Check fitness function preferring 0 and 180 rotations
     zero = mover._preferenceFunction(0)
     ninety = mover._preferenceFunction(90)
@@ -1685,6 +1685,7 @@ def __init__(self, paramN, paramA2):
       return "Movers.Test() MoverSingleHydrogenRotator pair: bad preference function"
     if zero - ninety < 1e-5:
       return "Movers.Test() MoverSingleHydrogenRotator pair: bad preference function"
+    '''
 
     # Check that one of the orientations has a dihedral angle of 0 compared to the original
     # 13-degree-off potential acceptor atom.
@@ -1779,6 +1780,7 @@ def __init__(self, paramN, paramA2):
     if not (h.xyz[2] == 1 and h.xyz[0]+math.sqrt(2) < 1e-5):
       return "Movers.Test() MoverSingleHydrogenRotator triple: bad H placement"
 
+    ''' Preference function has been removed from this class
     # Check fitness function preferring 180 and +/- 120 from there rotations away from
     # the angle away from 180 degrees.
     zero = mover._preferenceFunction(angle+0)
@@ -1791,6 +1793,7 @@ def __init__(self, paramN, paramA2):
       return "Movers.Test() MoverSingleHydrogenRotator triple: bad preference function"
     if zero - oneEighty < 1e-5:
       return "Movers.Test() MoverSingleHydrogenRotator triple: bad preference function"
+    '''
 
   except Exception as e:
     return "Movers.Test() MoverSingleHydrogenRotator triple: Exception during test: "+str(e)+"\n"+traceback.format_exc()
@@ -2760,13 +2763,6 @@ def __init__(self, paramN, paramA2):
     if 'AnglesNotAdjusted' != mover.PoseDescription(0, 0, False).split()[-1]:
       return "Movers.Test() MoverHisFlip: Unexpected angle description Flipped, got "+mover.PoseDescription(0, 0, False)
 
-    # Test whether Hydrogen placement optimization is happening when we do only flipped or only
-    # not flipped.
-    # @todo
-
-    # Test whether we can control fixup behavior when locking down in flipped state.
-    # @todo
-
   except Exception as e:
     return "Movers.Test() MoverHisFlip: Exception during test: "+str(e)+"\n"+traceback.format_exc()
 
diff --git a/mmtbx/reduce/Optimizers.py b/mmtbx/reduce/Optimizers.py
index a1b2408857..e0ce7c0a91 100644
--- a/mmtbx/reduce/Optimizers.py
+++ b/mmtbx/reduce/Optimizers.py
@@ -16,7 +16,7 @@
 from __future__ import print_function, nested_scopes, generators, division
 from __future__ import absolute_import
 
-import argparse, time, itertools
+import argparse, itertools, re
 
 from mmtbx.reduce import Movers
 from mmtbx.reduce import InteractionGraph
@@ -25,7 +25,6 @@
 
 from iotbx.map_model_manager import map_model_manager
 from iotbx.data_manager import DataManager
-from iotbx import pdb
 from iotbx.pdb import common_residue_names_get_class
 import mmtbx
 from scitbx.array_family import flex
@@ -108,9 +107,10 @@ def _ResNameAndID(a):
   chainID = a.parent().parent().parent().id
   resName = a.parent().resname.strip().upper()
   resID = str(a.parent().parent().resseq_as_int())
+  altLoc = a.parent().altloc
   # Don't print the code if it is a space (blank).
   insertionCode = a.parent().parent().icode.strip()
-  return "chain "+str(chainID)+" "+resName+" "+resID+insertionCode
+  return "chain "+str(chainID)+" "+altLoc+resName+" "+resID+insertionCode
 
 ##################################################################################
 # Helper classes
@@ -430,8 +430,8 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex, altID,
         # Compute the interaction graph, of which each connected component is a Clique.
         # Get a list of singleton Cliques and a list of other Cliques.  Keep separate lists
         # of the singletons and the groups.
-        self._interactionGraph, self._atomMoverSets = InteractionGraph.InteractionGraphAllPairs(self._movers, self._extraAtomInfo,
-          probeRadius=probeRadius)
+        self._interactionGraph, self._atomMoverSets = InteractionGraph.InteractionGraphAllPairs(self._movers,
+          self._extraAtomInfo, probeRadius=probeRadius)
         components = cca.connected_components( graph = self._interactionGraph )
         maxLen = 0
         singletonCliques = []   # Each entry is a list of integer indices into models with one entry
@@ -460,11 +460,12 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex, altID,
           for a in m.FixUp(0).atoms:
             moverAtoms.add(a)
 
-        # Get the excluded list for each atom in the set, making a dictionary
+        # Get the excluded list for each atom in the set, making a dictionary.
+        # We go at most 3 hops unless one end of the chain has a hydrogen.
         self._excludeDict = {}
         for a in moverAtoms:
           self._excludeDict[a] = mmtbx.probe.Helpers.getAtomsWithinNBonds(a,
-            bondedNeighborLists, self._extraAtomInfo, probeRadius, self._bondedNeighborDepth)
+            bondedNeighborLists, self._extraAtomInfo, probeRadius, self._bondedNeighborDepth, 3)
         self._infoString += _ReportTiming(self._verbosity, "determine excluded atoms")
 
         ################################################################################
@@ -568,7 +569,7 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex, altID,
         # Compute and store the initial score for each Mover in its info
         for m in self._movers:
           coarse = m.CoarsePositions()
-          score = coarse.preferenceEnergies[0] * self._preferenceMagnitude
+          score = self._preferenceMagnitude * coarse.preferenceEnergies[0]
           self._setMoverState(coarse, 0)
           for a in coarse.atoms:
             score += self._scoreAtom(a)
@@ -576,7 +577,7 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex, altID,
 
         ################################################################################
         # Call internal methods to optimize the single-element Cliques and then to optimize
-        # the multi-element Cliques and then to do indepenedent fine adjustment of all
+        # the multi-element Cliques and then to do independent fine adjustment of all
         # Cliques.  Subclasses should overload the called routines, but the global approach
         # taken here will be the same for all of them.  If we want to change the recipe
         # so that we can do global fine optimization, we'll do that here rather than in the
@@ -605,7 +606,7 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex, altID,
         # each of them, whether they are part of a multi-Mover Clique or not.
         self._fineLocations = {}
         for m in self._movers:
-          self._fineLocations[m] = 0
+          self._fineLocations[m] = None
         self._infoString += _VerboseCheck(self._verbosity, 1,"Fine optimization on all Movers\n")
         for m in self._movers:
           self._optimizeSingleMoverFine(m)
@@ -615,7 +616,7 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex, altID,
         # Print the final state and score for all Movers
         def _scoreMoverReportClash(self, m, index):
           coarse = m.CoarsePositions()
-          score = coarse.preferenceEnergies[index] * self._preferenceMagnitude
+          score = self._preferenceMagnitude * coarse.preferenceEnergies[index]
           clash = False
           for atom in coarse.atoms:
             maxRadiusWithoutProbe = self._extraAtomInfo.getMappingFor(atom).vdwRadius + self._maximumVDWRadius
@@ -803,7 +804,7 @@ def _optimizeSingleMoverCoarse(self, mover):
     # :side_effect: self._setMoverState() is called to put the Mover's atoms into its best state.
     # :side_effect: self._coarseLocations is set to the Mover's best state.
     # :side_effect: Changes the value of self._highScores[mover] to the score at the coarse position
-    # selected
+    # selected.
     coarse = mover.CoarsePositions()
     scores = coarse.preferenceEnergies[:]
     for i in range(len(scores)):
@@ -887,8 +888,8 @@ def _optimizeSingleMoverFine(self, mover):
 
       # Put the Mover into its final position (which may be back to its initial position)
       # and update the high score.
-      self._fineLocations[mover] = maxIndex
       if maxScore > self._highScores[mover]:
+        self._fineLocations[mover] = maxIndex
         self._infoString += _VerboseCheck(self._verbosity, 3,"Setting single Mover to fine orientation {}".format(maxIndex)+
           ", max score = {:.2f} (coarse score {:.2f})\n".format(maxScore,self._highScores[mover]))
         self._setMoverState(fine, maxIndex)
@@ -897,6 +898,7 @@ def _optimizeSingleMoverFine(self, mover):
         self._highScores[mover] = maxScore
       else:
         # Put us back to the initial coarse location and don't change the high score.
+        self._fineLocations[mover] = None
         self._infoString += _VerboseCheck(self._verbosity, 3,"Leaving single Mover at coarse orientation\n")
         self._setMoverState(coarse, self._coarseLocations[mover])
     return maxScore
@@ -1035,7 +1037,7 @@ def _PlaceMovers(self, atoms, rotatableHydrogenIDs, bondedNeighborLists, hParame
         if numH == 3:
           # See if the Carbon's other neighbor is attached to two other atoms (3 total).  If so,
           # then insert a MoverAromaticMethylRotator and if not, generate a MoverTetrahedralMethylRotator
-          # so that they Hydrogens will be staggered but do not add it to those to be optimized.
+          # so that the Hydrogens will be staggered but do not add it to those to be optimized.
           if len(bondedNeighborLists[neighbor]) == 3:
             try:
               self._movers.append(Movers.MoverAromaticMethylRotator(a, bondedNeighborLists, hParameters))
@@ -1087,7 +1089,7 @@ def _PlaceMovers(self, atoms, rotatableHydrogenIDs, bondedNeighborLists, hParame
           try:
             # Check to see if the state of this Mover has been specified. If so, place it in
             # the requested state and don't insert the Mover. If not, then insert the Mover.
-            s = _FindFlipState(a.parent().parent(), fs)
+            s = _FindFlipState(a, fs)
             if s is not None:
               self._infoString += _VerboseCheck(self._verbosity, 1,"Setting MoverAmideFlip for "+resNameAndID+": flipped = "+
                 str(s.flipped)+", angles adjusted = "+str(s.fixedUp)+"\n")
@@ -1204,7 +1206,7 @@ def _modifyIfNeeded(nitro, coarseNitroPos, hydro):
             # Check to see if the state of this Mover has been specified. If so, place it in
             # the requested state and insert the Mover as a non-flipping Histidine.
             # If not, then insert the Histidine flip Mover.
-            s = _FindFlipState(a.parent().parent(), fs)
+            s = _FindFlipState(a, fs)
             if s is not None:
               if s.flipped:
                 enabledFlips = 2
@@ -1332,16 +1334,17 @@ def _optimizeCliqueCoarse(self, clique):
     for curStateValues in _generateAllStates(numStates):
 
       # Set all movers to match the state list.
-      # @todo Optimize this so that it only changes states that differed from last time.
-      for i in range(len(movers)):
-        self._setMoverState(states[i], curStateValues[i])
-        self._coarseLocations[movers[i]] = curStateValues[i]
+      for i,m in enumerate(movers):
+        # Only change this Mover if it is different from the last time.
+        if self._coarseLocations[m] != curStateValues[i]:
+          self._setMoverState(states[i], curStateValues[i])
+          self._coarseLocations[m] = curStateValues[i]
 
       # Compute the score over all atoms in all Movers and see if it is the best.  If so,
       # update the best.
       score = 0
       for i in range(len(movers)):
-        score += states[i].preferenceEnergies[curStateValues[i]]
+        score += self._preferenceMagnitude * states[i].preferenceEnergies[curStateValues[i]]
         for a in states[i].atoms:
           score += self._scoreAtom(a)
       self._infoString += _VerboseCheck(self._verbosity, 5,"Score is {:.2f} at {}\n".format(score, curStateValues))
@@ -1350,14 +1353,14 @@ def _optimizeCliqueCoarse(self, clique):
         bestScore = score
         bestState = curStateValues[:]
 
-    # Put each Mover into its best state and compute its high-score value.
-    # Compute the best individual scores for these Movers for use in later fine-motion
+    # Put each Mover into its state in the best configuration and compute its high-score value.
+    # Store the individual scores for these Movers in the best config for use in later fine-motion
     # processing.  Return the total score.
     ret = 0.0
     for i,m in enumerate(movers):
       self._setMoverState(states[i], bestState[i])
       self._coarseLocations[m] = bestState[i]
-      self._highScores[m] = states[i].preferenceEnergies[curStateValues[i]]
+      self._highScores[m] = self._preferenceMagnitude * states[i].preferenceEnergies[curStateValues[i]]
       for a in states[i].atoms:
         self._highScores[m] += self._scoreAtom(a)
       self._infoString += _VerboseCheck(self._verbosity, 3,"Setting Mover in clique to coarse orientation {}".format(bestState[i])+
@@ -1380,7 +1383,7 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex, altID,
               ):
     """Constructor for _CliqueOptimizer.  This uses a recursive algorithm to break down the total
     clique into sets of smaller cliques.  It looks for a vertex cut in the Clique it is called with
-    that will separate the remaining vertices into two more more connected subcomponents.  It then tests
+    that will separate the remaining vertices into two or more more connected subcomponents.  It then tests
     each combined state of the set of Movers in the vertex cut to find the one with the best overall
     maximum score.  For each state, it first recursively optimizes all of the connected subcomponents
     and then (with each of the subcomponents in its optimal state) computes the score for the Movers in
@@ -1446,7 +1449,7 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex, altID,
 
   def _optimizeCliqueCoarse(self, clique):
     """
-    Looks for a vertex cut in the Clique that will separate the remaining vertices into two more
+    Looks for a vertex cut in the Clique that will separate the remaining vertices into two or more
     more connected subcomponents.  Test each combined state of the set of Movers in the vertex cut
     to find the one with the best overall maximum score.
       For each state, recursively optimize all of the connected subcomponents and then (with each
@@ -1491,10 +1494,11 @@ def _optimizeCliqueCoarse(self, clique):
         numStates.append(len(states[m].positions))
       for curStateValues in _generateAllStates(numStates):
         # Set all cutMovers to match the state list.
-        # @todo Optimize this so that it only changes states that differed from last time.
         for i,m in enumerate(cutMovers):
-          self._setMoverState(states[m], curStateValues[i])
-          self._coarseLocations[m] = curStateValues[i]
+          # Only change this Mover if it is different from the last time.
+          if self._coarseLocations[m] != curStateValues[i]:
+            self._setMoverState(states[m], curStateValues[i])
+            self._coarseLocations[m] = curStateValues[i]
 
         # Recursively compute the best score across all connected components in the cutGraph.
         # This will leave each subgraph in its best state for this set of cutMovers states.
@@ -1503,12 +1507,13 @@ def _optimizeCliqueCoarse(self, clique):
         for g in components:
           subMovers = [cutGraph.vertex_label(i) for i in g]
           subset = _subsetGraph(cutGraph, subMovers)
-          score += self._optimizeCliqueCoarse(subset)
+          # Call this exact same method, not a derived class method
+          score += _CliqueOptimizer._optimizeCliqueCoarse(self, subset)
 
         # Add the score over all atoms in the vertex-cut Movers and see if it is the best.  If so,
         # update the best.
         for i, m in enumerate(cutMovers):
-          score += states[m].preferenceEnergies[curStateValues[i]]
+          score += self._preferenceMagnitude * states[m].preferenceEnergies[curStateValues[i]]
           for a in states[m].atoms:
             score += self._scoreAtom(a)
         self._infoString += _VerboseCheck(self._verbosity, 5,"Score is {:.2f} at {}\n".format(score, curStateValues))
@@ -1525,7 +1530,7 @@ def _optimizeCliqueCoarse(self, clique):
       for i,m in enumerate(movers):
         self._setMoverState(states[m], bestState[i])
         self._coarseLocations[m] = bestState[i]
-        self._highScores[m] = states[m].preferenceEnergies[bestState[i]]
+        self._highScores[m] = self._preferenceMagnitude * states[m].preferenceEnergies[bestState[i]]
         for a in states[m].atoms:
           self._highScores[m] += self._scoreAtom(a)
         self._infoString += _VerboseCheck(self._verbosity, 3,"Setting Mover in clique to coarse orientation {}".format(bestState[i])+
@@ -1552,7 +1557,7 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex = 0, altID = None
                 verbosity = 1
               ):
     """Constructor for FastOptimizer.  This uses the same algorithm as the
-    parent-class but first constructs a cache for every atom in every Mover
+    parent class but first constructs a cache for every atom in every Mover
     of all the Movers whose positions can affect its answer.  The _scoreAtom() method
     is overridden to use this cached value in clique optimization (but not in singleton
     or fine optimization) when it has already been computed for a given configuration
@@ -1615,6 +1620,10 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex = 0, altID = None
     # it calls.
     self._doScoreCaching = False
 
+    # Keep track of the ratio of calculated vs. cached results.
+    self._numCalculated = 0
+    self._numCached = 0
+
     super(FastOptimizer, self).__init__(probePhil, addFlipMovers, model, modelIndex = modelIndex, altID = altID,
                 bondedNeighborDepth = bondedNeighborDepth,
                 probeRadius = probeRadius, useNeutronDistances = useNeutronDistances, probeDensity = probeDensity,
@@ -1622,6 +1631,13 @@ def __init__(self, probePhil, addFlipMovers, model, modelIndex = 0, altID = None
                 nonFlipPreference = nonFlipPreference, skipBondFixup = skipBondFixup,
                 flipStates = flipStates, verbosity = verbosity)
 
+    if self._numCalculated > 0:
+      self._infoString += _VerboseCheck(self._verbosity, 1,
+          'Calculated : cached atom scores: {} : {}; fraction calculated {:.2f}\n'.format(
+          self._numCalculated, self._numCached,
+          self._numCalculated/(self._numCalculated + self._numCached)))
+
+
   def _initializeAlternate(self):
     # Ensure that we have a per-atom _scoreCache dictionary that will store already-computed
     # results for a given atom based on the configurations of the Movers that can affect its
@@ -1644,11 +1660,15 @@ def _scoreAtom(self, atom):
       # compute and store it and then return that value.
       state = tuple([self._coarseLocations[m] for m in self._atomMoverSets[atom]])
       try:
+        self._numCached += 1
         return self._scoreCache[atom][state]
       except Exception:
+        self._numCached -= 1      # Undo the increment above
+        self._numCalculated += 1
         self._scoreCache[atom][state] = super(FastOptimizer, self)._scoreAtom(atom)
         return self._scoreCache[atom][state]
     else:
+      self._numCalculated += 1
       return super(FastOptimizer, self)._scoreAtom(atom)
 
   def _optimizeCliqueCoarse(self, clique):
@@ -1704,19 +1724,22 @@ def _ParseFlipStates(fs):
   return ret
 
 
-def _FindFlipState(rg, flipStates):
-  '''Is a residue group in the list of residues that have flip states? If so, return the associated state
-  :param rg: residue group
+def _FindFlipState(a, flipStates):
+  '''Is an atom in the list of residues that have flip states? If so, return the associated state
+  :param a: atom
   :param flipStates: List of FlipMoverState objects
   :return: FlipMoverState if the residue is in the list, None if not
   '''
+  ag = a.parent()
+  rg = ag.parent()
+  chain = rg.parent()
   for fs in flipStates:
-    chain = rg.parent()
     modelId = chain.parent().id
     # We must offset the index of the model by 1 to get to the 1-based model ID
     if ( (modelId == fs.modelId + 1 or modelId == '') and
           (chain.id == fs.chain) and
           (rg.resseq_as_int() == fs.resId) and
+          (ag.altloc.strip() == '' or fs.altId == '.' or ag.altloc.strip().lower() == fs.altId.lower()) and
           (rg.icode.strip() == fs.iCode.strip())
         ):
       return fs
@@ -1837,6 +1860,75 @@ def _generateAllStates(numStates):
 ##################################################################################
 # Test function and associated data and helpers to verify that all functions behave properly.
 
+def _optimizeFragment(pdb_raw):
+  """Returns an optimizer constructed based on the raw PDB data snippet passed in.
+  :param pdb_raw: A string that includes a snippet of a PDB file. Should have no alternates.
+  :return: FastOptimizer initialized based on the snippet.
+  """
+  dm = DataManager(['model'])
+  dm.process_model_str("my_data.pdb", pdb_raw)
+  model = dm.get_model()
+
+  # Make sure we have a valid unit cell.  Do this before we add hydrogens to the model
+  # to make sure we have a valid unit cell.
+  model = shift_and_box_model(model = model)
+
+  # Add Hydrogens to the model
+  reduce_add_h_obj = reduce_hydrogen.place_hydrogens(model = model)
+  reduce_add_h_obj.run()
+  model = reduce_add_h_obj.get_model()
+
+  # Interpret the model after adding Hydrogens to it so that
+  # all of the needed fields are filled in when we use them below.
+  # @todo Remove this once place_hydrogens() does all the interpretation we need.
+  p = mmtbx.model.manager.get_default_pdb_interpretation_params()
+  p.pdb_interpretation.allow_polymer_cross_special_position=True
+  p.pdb_interpretation.clash_guard.nonbonded_distance_threshold=None
+  p.pdb_interpretation.proceed_with_excessive_length_bonds=True
+  model.process(make_restraints=True,pdb_interpretation_params=p) # make restraints
+
+  # Get the first model in the hierarchy.
+  firstModel = model.get_hierarchy().models()[0]
+
+  # Get the list of alternate conformation names present in all chains for this model.
+  alts = AlternatesInModel(firstModel)
+
+  # Get the atoms from the first conformer in the first model (the empty string is the name
+  # of the first conformation in the model; if there is no empty conformation, then it will
+  # pick the first available conformation for each atom group.
+  atoms = GetAtomsForConformer(firstModel, "")
+
+  # Get the Cartesian positions of all of the atoms we're considering for this alternate
+  # conformation.
+  carts = flex.vec3_double()
+  for a in atoms:
+    carts.append(a.xyz)
+
+  # Get the bond proxies for the atoms in the model and conformation we're using and
+  # use them to determine the bonded neighbor lists.
+  bondProxies = model.get_restraints_manager().geometry.get_all_bond_proxies(sites_cart = carts)[0]
+  bondedNeighborLists = Helpers.getBondedNeighborLists(atoms, bondProxies)
+
+  # Get the probeExt.ExtraAtomInfo needed to determine which atoms are potential acceptors.
+  class philLike:
+    def __init__(self, useImplicitHydrogenDistances = False):
+      self.implicit_hydrogens = useImplicitHydrogenDistances
+      self.set_polar_hydrogen_radius = True
+  probePhil = philLike(False)
+  ret = Helpers.getExtraAtomInfo(model = model, bondedNeighborLists = bondedNeighborLists,
+      useNeutronDistances=False,probePhil=probePhil)
+  extra = ret.extraAtomInfo
+
+  # Also compute the maximum VDW radius among all atoms.
+  maxVDWRad = 1
+  for a in atoms:
+    maxVDWRad = max(maxVDWRad, extra.getMappingFor(a).vdwRadius)
+
+  # Optimization will place the movers, which should be none because the Histidine flip
+  # will be constrained by the ionic bonds.
+  probePhil = philLike(False)
+  return FastOptimizer(probePhil, True, model)
+
 def Test(inFileName = None, dumpAtoms = False):
   """Test function for all functions provided above.
   :param inFileName: Name of a PDB or CIF file to load (default makes a small molecule)
@@ -1876,6 +1968,124 @@ def Test(inFileName = None, dumpAtoms = False):
   #========================================================================
   # Unit tests for each type of Optimizer.
 
+  ################################################################################
+  # Test using a snippet from 7c31 to make sure the single-hydrogen rotator sets
+  # the orientation of the Hydrogen to make good contact with a nearby Oxygen. The B
+  # alternate has been removed and the A moved to the main alternate. This is not
+  # close enough for a good Hydrogen bond, so it does not point straight at the Oxygen.
+
+  pdb_7c31_two_residues = """\
+CRYST1   27.854   27.854   99.605  90.00  90.00  90.00 P 43          8
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.035901  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.035901  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.010040        0.00000
+ATOM     68  N   SER A   5     -31.155  49.742   0.887  1.00 10.02           N
+ATOM     69  CA  SER A   5     -32.274  48.937   0.401  1.00  9.76           C
+ATOM     70  C   SER A   5     -33.140  49.851  -0.454  1.00  9.47           C
+ATOM     71  O   SER A   5     -33.502  50.939  -0.012  1.00 10.76           O
+ATOM     72  CB  SER A   5     -33.086  48.441   1.599  1.00 12.34           C
+ATOM     73  OG  SER A   5     -34.118  47.569   1.179  1.00 19.50           O
+ATOM    758  N   VAL B   2     -34.430  42.959   3.043  1.00 19.95           N
+ATOM    759  CA  VAL B   2     -32.994  42.754   2.904  0.50 18.09           C
+ATOM    760  C   VAL B   2     -32.311  44.083   2.629  1.00 17.65           C
+ATOM    761  O   VAL B   2     -32.869  44.976   1.975  1.00 18.97           O
+ATOM    762  CB  VAL B   2     -32.703  41.738   1.778  0.50 19.70           C
+ATOM    763  CG1 VAL B   2     -33.098  42.307   0.419  0.50 19.61           C
+ATOM    764  CG2 VAL B   2     -31.224  41.334   1.755  0.50 16.49           C
+TER    1447      CYS B  47
+END
+"""
+  opt = _optimizeFragment(pdb_7c31_two_residues)
+  movers = opt._movers
+  if len(movers) != 1:
+    return "Optimizers.Test(): Incorrect number of Movers for single-hydrogen rotator test: " + str(len(movers))
+
+  # See what the pose angle is on the Mover. It should be 111 degrees, and is reported
+  # after 'pose Angle '.
+  angle = int(re.search('(?<=pose Angle )\d+', opt.getInfo()).group(0))
+  if angle != 62:
+    return "Optimizers.Test(): Unexpected angle ("+str(angle)+") for single-hydrogen rotator, expected 62"
+
+  ################################################################################
+  # Test using a modified snippet from 7c31 to make sure the single-hydrogen rotator sets
+  # the orientation of the Hydrogen to make a good hydrogen bond with a nearby Oxygen.
+  # The above example has been stripped down and has the Oxygen moved closer
+
+  pdb_7c31_two_residues_close_O = """\
+CRYST1   27.854   27.854   99.605  90.00  90.00  90.00 P 43          8
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.035901  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.035901  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.010040        0.00000
+ATOM     68  N   SER A   5     -31.155  49.742   0.887  1.00 10.02           N
+ATOM     69  CA  SER A   5     -32.274  48.937   0.401  1.00  9.76           C
+ATOM     70  C   SER A   5     -33.140  49.851  -0.454  1.00  9.47           C
+ATOM     71  O   SER A   5     -33.502  50.939  -0.012  1.00 10.76           O
+ATOM     72  CB  SER A   5     -33.086  48.441   1.599  1.00 12.34           C
+ATOM     73  OG  SER A   5     -34.118  47.569   1.179  1.00 19.50           O
+ATOM    761  O   VAL B   2     -33.862  46.385   1.577  1.00 18.97           O
+TER    1447      CYS B  47
+END
+"""
+  opt = _optimizeFragment(pdb_7c31_two_residues_close_O)
+  movers = opt._movers
+  if len(movers) != 1:
+    return "Optimizers.Test(): Incorrect number of Movers for single-hydrogen rotator H-bond test: " + str(len(movers))
+
+  # See what the pose angle is on the Mover. It should be 111 degrees, and is reported
+  # after 'pose Angle '.
+  angle = int(re.search('(?<=pose Angle )\d+', opt.getInfo()).group(0))
+  if angle != 109:
+    return "Optimizers.Test(): Unexpected angle ("+str(angle)+") for single-hydrogen rotator H-bond, expected 109"
+
+  ################################################################################
+  # Test using a modified snippet from 1ubq to make sure the NH3 rotator sets
+  # the orientation as expected to align with a nearby Oxygen and waters.
+
+  pdb_1ubq_two_residues_close_O = """\
+CRYST1   50.840   42.770   28.950  90.00  90.00  90.00 P 21 21 21    4
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.019670  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.023381  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.034542        0.00000
+ATOM      1  N   MET A   1      27.340  24.430   2.614  1.00  9.67           N
+ATOM      2  CA  MET A   1      26.266  25.413   2.842  1.00 10.38           C
+ATOM      3  C   MET A   1      26.913  26.639   3.531  1.00  9.62           C
+ATOM      4  O   MET A   1      27.886  26.463   4.263  1.00  9.62           O
+ATOM      5  CB  MET A   1      25.112  24.880   3.649  1.00 13.77           C
+ATOM      6  CG  MET A   1      25.353  24.860   5.134  1.00 16.29           C
+ATOM      7  SD  MET A   1      23.930  23.959   5.904  1.00 17.17           S
+ATOM      8  CE  MET A   1      24.447  23.984   7.620  1.00 16.11           C
+ATOM    127  N   VAL A  17      30.310  25.458   5.384  1.00  8.99           N
+ATOM    128  CA  VAL A  17      30.288  24.245   6.193  1.00  8.85           C
+ATOM    129  C   VAL A  17      29.279  23.227   5.641  1.00  8.04           C
+ATOM    130  O   VAL A  17      28.478  23.522   4.725  1.00  8.99           O
+ATOM    131  CB  VAL A  17      29.903  24.590   7.665  1.00  9.78           C
+ATOM    132  CG1 VAL A  17      30.862  25.496   8.389  1.00 12.05           C
+ATOM    133  CG2 VAL A  17      28.476  25.135   7.705  1.00 10.54           C
+TER     603      GLY A  76
+HETATM  625  O   HOH A  98      25.928  21.774   2.325  1.00 13.70           O
+HETATM  637  O   HOH A 110      28.824  25.094   0.886  0.77 36.99           O
+END
+"""
+  opt = _optimizeFragment(pdb_1ubq_two_residues_close_O)
+  movers = opt._movers
+  if len(movers) != 1:
+    return "Optimizers.Test(): Incorrect number of Movers for NH3 rotator test: " + str(len(movers))
+
+  # See what the pose angle is on the Mover. It should be 163 degrees, and is reported
+  # after 'pose Angle '.
+  angle = int(re.search('(?<=pose Angle )\d+', opt.getInfo()).group(0))
+  if angle != 163:
+    return "Optimizers.Test(): Unexpected angle ("+str(angle)+") for NH3 rotator, expected 163"
+
   ################################################################################
   # Test using snippet from 1xso to ensure that the Histidine placement code will lock down the
   # Histidine, set its Nitrogen states, and mark its Hydrogens for deletion.  To
@@ -1980,9 +2190,6 @@ def __init__(self, useImplicitHydrogenDistances = False):
     if a.element.upper() == "ZN":
       a.xyz = origPositionZN
 
-  # Get the spatial-query information needed to quickly determine which atoms are nearby
-  sq = probeExt.SpatialQuery(atoms)
-
   # Optimization will place the movers, which should be none because the Histidine flip
   # will be constrained by the ionic bonds.
   probePhil = philLike(False)
@@ -2003,8 +2210,165 @@ def __init__(self, useImplicitHydrogenDistances = False):
         return 'Optimizers.Test(): '+name+' in 1xso Histidine test was not set for deletion'
 
   #========================================================================
-  # Check a case where a HisFlip would be flipped and locked down and have its Hydrogen removed.
-  # @todo
+  # Check a clique with multiple elements to be sure that it was properly
+  # globally optimized. This is a set of ACT residues that are offset
+  # such that they want to line up the same way.  A carbon is placed to
+  # force the orientaion away from the initial solution.
+  pdb_multi_act = (
+"""
+CRYST1   93.586  127.886  251.681  90.00  90.00  90.00 I 2 2 2
+SCALE1      0.010685  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.007819  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.003973        0.00000
+HETATM    1  C   ACT A   1       6.494 -47.273 -37.006  1.00 16.65           C
+HETATM    2  O   ACT A   1       5.654 -47.981 -37.645  1.00 15.33           O
+HETATM    3  CH3 ACT A   1       6.790 -47.410 -35.548  1.00 17.13           C
+HETATM    4  OXT ACT A   1       7.110 -46.425 -37.699  1.00 17.05           O
+HETATM  101  C   ACT A 101       8.994 -49.773 -37.606  1.00 16.65           C
+HETATM  102  O   ACT A 101       8.154 -50.481 -38.245  1.00 15.33           O
+HETATM  103  CH3 ACT A 101       9.290 -49.910 -36.148  1.00 17.13           C
+HETATM  104  OXT ACT A 101       9.610 -48.925 -38.299  1.00 17.05           O
+HETATM  111  C   ACT A 111      11.494 -52.273 -38.206  1.00 16.65           C
+HETATM  112  O   ACT A 111      10.654 -52.981 -38.845  1.00 15.33           O
+HETATM  113  CH3 ACT A 111      11.790 -52.410 -36.748  1.00 17.13           C
+HETATM  114  OXT ACT A 111      12.110 -51.425 -38.899  1.00 17.05           O
+HETATM  121  C   ACT A 121      13.994 -54.773 -38.806  1.00 16.65           C
+HETATM  122  O   ACT A 121      13.154 -55.481 -39.445  1.00 15.33           O
+HETATM  123  CH3 ACT A 121      14.290 -54.910 -37.348  1.00 17.13           C
+HETATM  124  OXT ACT A 121      14.610 -53.925 -39.499  1.00 17.05           O
+HETATM  131  C   ACT A 131      16.494 -57.273 -39.406  1.00 16.65           C
+HETATM  132  O   ACT A 131      15.654 -57.981 -40.045  1.00 15.33           O
+HETATM  133  CH3 ACT A 131      16.790 -57.410 -37.948  1.00 17.13           C
+HETATM  134  OXT ACT A 131      17.110 -56.425 -40.099  1.00 17.05           O
+HETATM  141  C   ACT A 141      18.994 -59.773 -40.006  1.00 16.65           C
+HETATM  142  O   ACT A 141      18.154 -60.481 -40.645  1.00 15.33           O
+HETATM  143  CH3 ACT A 141      19.290 -59.910 -38.548  1.00 17.13           C
+HETATM  144  OXT ACT A 141      19.610 -58.925 -40.699  1.00 17.05           O
+HETATM  200  C   ACT A 200       4.500 -45.500 -35.000  1.00 16.65           C
+END
+"""
+    )
+  dm = DataManager(['model'])
+  dm.process_model_str("pdb_multi_act.pdb",pdb_multi_act)
+  model = dm.get_model()
+
+  # Make sure we have a valid unit cell.  Do this before we add hydrogens to the model
+  # to make sure we have a valid unit cell.
+  model = shift_and_box_model(model = model)
+
+  # Add Hydrogens to the model
+  reduce_add_h_obj = reduce_hydrogen.place_hydrogens(model = model)
+  reduce_add_h_obj.run()
+  model = reduce_add_h_obj.get_model()
+
+  # Interpret the model after adding Hydrogens to it so that
+  # all of the needed fields are filled in when we use them below.
+  # @todo Remove this once place_hydrogens() does all the interpretation we need.
+  p = mmtbx.model.manager.get_default_pdb_interpretation_params()
+  p.pdb_interpretation.allow_polymer_cross_special_position=True
+  p.pdb_interpretation.clash_guard.nonbonded_distance_threshold=None
+  p.pdb_interpretation.proceed_with_excessive_length_bonds=True
+  model.process(make_restraints=True,pdb_interpretation_params=p) # make restraints
+
+  # Get the first model in the hierarchy.
+  firstModel = model.get_hierarchy().models()[0]
+
+  # Get the list of alternate conformation names present in all chains for this model.
+  alts = AlternatesInModel(firstModel)
+
+  # Get the atoms from the first conformer in the first model (the empty string is the name
+  # of the first conformation in the model; if there is no empty conformation, then it will
+  # pick the first available conformation for each atom group.
+  atoms = GetAtomsForConformer(firstModel, "")
+
+  # Get the Cartesian positions of all of the atoms we're considering for this alternate
+  # conformation.
+  carts = flex.vec3_double()
+  for a in atoms:
+    carts.append(a.xyz)
+
+  # Get the bond proxies for the atoms in the model and conformation we're using and
+  # use them to determine the bonded neighbor lists.
+  bondProxies = model.get_restraints_manager().geometry.get_all_bond_proxies(sites_cart = carts)[0]
+  bondedNeighborLists = Helpers.getBondedNeighborLists(atoms, bondProxies)
+
+  # Get the probeExt.ExtraAtomInfo needed to determine which atoms are potential acceptors.
+  class philLike:
+    def __init__(self, useImplicitHydrogenDistances = False):
+      self.implicit_hydrogens = useImplicitHydrogenDistances
+      self.set_polar_hydrogen_radius = True
+  probePhil = philLike(False)
+  ret = Helpers.getExtraAtomInfo(model = model, bondedNeighborLists = bondedNeighborLists,
+      useNeutronDistances=False,probePhil=probePhil)
+  extra = ret.extraAtomInfo
+
+  # Also compute the maximum VDW radius among all atoms.
+  maxVDWRad = 1
+  for a in atoms:
+    maxVDWRad = max(maxVDWRad, extra.getMappingFor(a).vdwRadius)
+
+  # Optimization will place the movers. Make sure we got as many as we expected.
+  # Make sure that the orientation for all of the movers is correct.
+  # Test with each type of optimizer, from the base to the more derived, so
+  # that we find out about failures on the base classes first.
+  probePhil = philLike(False)
+
+  print('Testing _BruteForceOptimizer')
+  opt = _BruteForceOptimizer(probePhil, True, model, probeRadius=0.25, modelIndex = 0, altID = None,
+                bondedNeighborDepth = 4,
+                useNeutronDistances = False,
+                probeDensity = 16.0,
+                minOccupancy = 0.02,
+                preferenceMagnitude = 1.0,
+                nonFlipPreference = 0.5,
+                skipBondFixup = False,
+                flipStates = '',
+                verbosity = 1)
+  movers = opt._movers
+  if len(movers) != 6:
+    return "Optimizers.Test(): Incorrect number of Movers for _BruteForceOptimizer multi-ACT test: " + str(len(movers))
+  res = re.findall('pose Angle [-+]?\d+', opt.getInfo())
+  for r in res:
+    if not 'pose Angle 90' in r:
+      return "Optimizers.Test(): Unexpected angle ("+str(r)+") for _BruteForceOptimizer multi-ACT test"
+
+  print('Testing _CliqueOptimizer')
+  opt = _CliqueOptimizer(probePhil, True, model, probeRadius=0.25, modelIndex = 0, altID = None,
+                bondedNeighborDepth = 4,
+                useNeutronDistances = False,
+                probeDensity = 16.0,
+                minOccupancy = 0.02,
+                preferenceMagnitude = 1.0,
+                nonFlipPreference = 0.5,
+                skipBondFixup = False,
+                flipStates = '',
+                verbosity = 1)
+  movers = opt._movers
+  if len(movers) != 6:
+    return "Optimizers.Test(): Incorrect number of Movers for _CliqueOptimizer multi-ACT test: " + str(len(movers))
+  res = re.findall('pose Angle [-+]?\d+', opt.getInfo())
+  for r in res:
+    if not 'pose Angle 90' in r:
+      return "Optimizers.Test(): Unexpected angle ("+str(r)+") for _CliqueOptimizer multi-ACT test"
+
+  print('Testing FastOptimizer')
+  opt = FastOptimizer(probePhil, True, model, probeRadius=0.25, modelIndex = 0, altID = None,
+                bondedNeighborDepth = 4,
+                useNeutronDistances = False,
+                probeDensity = 16.0,
+                minOccupancy = 0.02,
+                preferenceMagnitude = 1.0,
+                nonFlipPreference = 0.5,
+                skipBondFixup = False,
+                flipStates = '',
+                verbosity = 1)
+  movers = opt._movers
+  if len(movers) != 6:
+    return "Optimizers.Test(): Incorrect number of Movers for FastOptimizer multi-ACT test: " + str(len(movers))
+  res = re.findall('pose Angle [-+]?\d+', opt.getInfo())
+  for r in res:
+    if not 'pose Angle 90' in r:
+      return "Optimizers.Test(): Unexpected angle ("+str(r)+") for FastOptimizer multi-ACT test"
 
   #========================================================================
   # Check a case where an AmideFlip would be locked down and have its Hydrogen removed.
@@ -2014,7 +2378,100 @@ def __init__(self, useImplicitHydrogenDistances = False):
   # Check that the occupancy and B-factor cut-offs for water Oxygens are causing them
   # to be ignored in the calculations by putting an atom in the way and making sure it
   # is ignored.
-  # @todo
+  pdb_b_factor = (
+"""
+CRYST1   93.586  127.886  251.681  90.00  90.00  90.00 I 2 2 2
+SCALE1      0.010685  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.007819  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.003973        0.00000
+HETATM    1  C   ACT A   1       6.494 -47.273 -37.006  1.00 16.65           C
+HETATM    2  O   ACT A   1       5.654 -47.981 -37.645  1.00 15.33           O
+HETATM    3  CH3 ACT A   1       6.790 -47.410 -35.548  1.00 17.13           C
+HETATM    4  OXT ACT A   1       7.110 -46.425 -37.699  1.00 17.05           O
+HETATM  200  O   HOH A 200       8.200 -48.610 -35.148  1.00 50.00           O
+HETATM  201  O   HOH A 201       5.800 -46.810 -34.548  1.00 20.00           O
+END
+"""
+    )
+  dm = DataManager(['model'])
+  dm.process_model_str("pdb_b_factor.pdb",pdb_b_factor)
+  model = dm.get_model()
+
+  # Make sure we have a valid unit cell.  Do this before we add hydrogens to the model
+  # to make sure we have a valid unit cell.
+  model = shift_and_box_model(model = model)
+
+  # Add Hydrogens to the model
+  reduce_add_h_obj = reduce_hydrogen.place_hydrogens(model = model)
+  reduce_add_h_obj.run()
+  model = reduce_add_h_obj.get_model()
+
+  # Interpret the model after adding Hydrogens to it so that
+  # all of the needed fields are filled in when we use them below.
+  # @todo Remove this once place_hydrogens() does all the interpretation we need.
+  p = mmtbx.model.manager.get_default_pdb_interpretation_params()
+  p.pdb_interpretation.allow_polymer_cross_special_position=True
+  p.pdb_interpretation.clash_guard.nonbonded_distance_threshold=None
+  p.pdb_interpretation.proceed_with_excessive_length_bonds=True
+  model.process(make_restraints=True,pdb_interpretation_params=p) # make restraints
+
+  # Get the first model in the hierarchy.
+  firstModel = model.get_hierarchy().models()[0]
+
+  # Get the list of alternate conformation names present in all chains for this model.
+  alts = AlternatesInModel(firstModel)
+
+  # Get the atoms from the first conformer in the first model (the empty string is the name
+  # of the first conformation in the model; if there is no empty conformation, then it will
+  # pick the first available conformation for each atom group.
+  atoms = GetAtomsForConformer(firstModel, "")
+
+  # Get the Cartesian positions of all of the atoms we're considering for this alternate
+  # conformation.
+  carts = flex.vec3_double()
+  for a in atoms:
+    carts.append(a.xyz)
+
+  # Get the bond proxies for the atoms in the model and conformation we're using and
+  # use them to determine the bonded neighbor lists.
+  bondProxies = model.get_restraints_manager().geometry.get_all_bond_proxies(sites_cart = carts)[0]
+  bondedNeighborLists = Helpers.getBondedNeighborLists(atoms, bondProxies)
+
+  # Get the probeExt.ExtraAtomInfo needed to determine which atoms are potential acceptors.
+  class philLike:
+    def __init__(self, useImplicitHydrogenDistances = False):
+      self.implicit_hydrogens = useImplicitHydrogenDistances
+      self.set_polar_hydrogen_radius = True
+  probePhil = philLike(False)
+  ret = Helpers.getExtraAtomInfo(model = model, bondedNeighborLists = bondedNeighborLists,
+      useNeutronDistances=False,probePhil=probePhil)
+  extra = ret.extraAtomInfo
+
+  # Also compute the maximum VDW radius among all atoms.
+  maxVDWRad = 1
+  for a in atoms:
+    maxVDWRad = max(maxVDWRad, extra.getMappingFor(a).vdwRadius)
+
+  # Optimization will place the movers. Make sure we got as many as we expected.
+  # Make sure that the orientation for all of the movers is correct.
+  probePhil = philLike(False)
+  opt = FastOptimizer(probePhil, True, model, probeRadius=0.25, modelIndex = 0, altID = None,
+                bondedNeighborDepth = 4,
+                useNeutronDistances = False,
+                probeDensity = 16.0,
+                minOccupancy = 0.02,
+                preferenceMagnitude = 1.0,
+                nonFlipPreference = 0.5,
+                skipBondFixup = False,
+                flipStates = '',
+                verbosity = 1)
+  movers = opt._movers
+  if len(movers) != 1:
+    return "Optimizers.Test(): Incorrect number of Movers for B-factor test: " + str(len(movers))
+  res = re.findall('pose Angle [-+]?\d+', opt.getInfo())
+  for r in res:
+    if not 'pose Angle 90' in r:
+      return "Optimizers.Test(): Unexpected angle ("+str(r)+") for B-factor test"
 
   #========================================================================
   # Generate an example data model with a small molecule in it or else read
@@ -2047,11 +2504,9 @@ def __init__(self, useImplicitHydrogenDistances = False):
   p.pdb_interpretation.proceed_with_excessive_length_bonds=True
   model.process(make_restraints=True, pdb_interpretation_params=p) # make restraints
 
-  # Run each type of optimizer on the model.
-  # @todo Verify that they all get equivalent results
-  print('Testing _BruteForceOptimizer')
-  opt = _BruteForceOptimizer(probePhil, True, model,probeRadius=0.25, modelIndex = 0, altID = None,
-                bondedNeighborDepth = 3,
+  # Run each type of optimizer on the model to make sure they don't crash.
+  opt = _BruteForceOptimizer(probePhil, True, model, probeRadius=0.25, modelIndex = 0, altID = None,
+                bondedNeighborDepth = 4,
                 useNeutronDistances = False,
                 probeDensity = 16.0,
                 minOccupancy = 0.02,
@@ -2060,9 +2515,8 @@ def __init__(self, useImplicitHydrogenDistances = False):
                 skipBondFixup = False,
                 flipStates = '',
                 verbosity = 1)
-  print('Testing _CliqueOptimizer')
-  opt = _CliqueOptimizer(probePhil, True, model,probeRadius=0.25, modelIndex = 0, altID = None,
-                bondedNeighborDepth = 3,
+  opt = _CliqueOptimizer(probePhil, True, model, probeRadius=0.25, modelIndex = 0, altID = None,
+                bondedNeighborDepth = 4,
                 useNeutronDistances = False,
                 probeDensity = 16.0,
                 minOccupancy = 0.02,
@@ -2071,11 +2525,7 @@ def __init__(self, useImplicitHydrogenDistances = False):
                 skipBondFixup = False,
                 flipStates = '',
                 verbosity = 1)
-  print('Testing FastOptimizer')
-  opt = FastOptimizer(probePhil, True, model,probeRadius=0.25)
-
-  # Test specifying the flip states.
-  # @todo
+  opt = FastOptimizer(probePhil, True, model, probeRadius=0.25)
 
   # Write debugging output if we've been asked to
   if dumpAtoms:
diff --git a/mmtbx/regression/tst_reduce.py b/mmtbx/regression/tst_reduce.py
index e7538b9cc4..8fd57f3d0f 100644
--- a/mmtbx/regression/tst_reduce.py
+++ b/mmtbx/regression/tst_reduce.py
@@ -2,7 +2,7 @@
 # This is a test program to validate that the Python wrapping of Reduce worked.
 #
 
-#                Copyright 2021  Richardson Lab at Duke University
+#                Copyright 2021-2023  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.
@@ -21,6 +21,12 @@
 from mmtbx.reduce import Movers
 from mmtbx.reduce import InteractionGraph
 from mmtbx.reduce import Optimizers
+from mmtbx.programs import reduce2
+from iotbx.cli_parser import run_program
+from six.moves import cStringIO as StringIO
+import os.path as op
+import os
+import math
 
 def RunReduceTests():
 
@@ -47,6 +53,454 @@ def RunReduceTests():
 
   return ret
 
+# Each test case has a name, a raw PDB file, a chain ID, a residue ID, an atom name,
+# a list of possible locations for the atom (for sets of 3 hydrogens, any 60-degree
+# rotation is equivalent in score), a maximum distance threshold, and a list of extra
+# arguments to make to Reduce2 for this case.
+testCases = [
+
+  ["7c31_single_hydrogen_rotator",
+   """\
+CRYST1   27.854   27.854   99.605  90.00  90.00  90.00 P 43          8
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.035901  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.035901  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.010040        0.00000
+ATOM     68  N   SER A   5     -31.155  49.742   0.887  1.00 10.02           N
+ATOM     69  CA  SER A   5     -32.274  48.937   0.401  1.00  9.76           C
+ATOM     70  C   SER A   5     -33.140  49.851  -0.454  1.00  9.47           C
+ATOM     71  O   SER A   5     -33.502  50.939  -0.012  1.00 10.76           O
+ATOM     72  CB  SER A   5     -33.086  48.441   1.599  1.00 12.34           C
+ATOM     73  OG  SER A   5     -34.118  47.569   1.179  1.00 19.50           O
+ATOM    758  N   VAL B   2     -34.430  42.959   3.043  1.00 19.95           N
+ATOM    759  CA  VAL B   2     -32.994  42.754   2.904  0.50 18.09           C
+ATOM    760  C   VAL B   2     -32.311  44.083   2.629  1.00 17.65           C
+ATOM    761  O   VAL B   2     -32.869  44.976   1.975  1.00 18.97           O
+ATOM    762  CB  VAL B   2     -32.703  41.738   1.778  0.50 19.70           C
+ATOM    763  CG1 VAL B   2     -33.098  42.307   0.419  0.50 19.61           C
+ATOM    764  CG2 VAL B   2     -31.224  41.334   1.755  0.50 16.49           C
+TER    1447      CYS B  47
+END
+""",
+   "A",
+   5,
+   "HG",
+   [ (-33.788, 46.891, 0.809)
+   ],
+   0.1,
+   []
+  ],
+
+  ["1xso_amide_unflipped",
+   """\
+CRYST1   73.450   68.940   58.760  90.00  90.00  90.00 P 21 21 21    8
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.013615  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.014505  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.017018        0.00000
+MTRIX1   1  0.987830 -0.010290 -0.155220        1.43777    1
+MTRIX2   1 -0.008360 -0.999880  0.013070       33.15571    1
+MTRIX3   1 -0.155330 -0.011610 -0.987790       17.07902    1
+ATOM    147  N   GLN A  22      32.268  17.697  -8.014  1.00 15.32           N
+ATOM    148  CA  GLN A  22      33.617  18.261  -8.189  1.00 14.46           C
+ATOM    149  C   GLN A  22      34.546  17.130  -8.669  1.00 18.34           C
+ATOM    150  O   GLN A  22      34.525  16.056  -8.055  1.00 20.98           O
+ATOM    151  CB  GLN A  22      34.120  18.672  -6.771  1.00 15.81           C
+ATOM    152  CG  GLN A  22      35.452  19.443  -6.847  1.00 17.90           C
+ATOM    153  CD  GLN A  22      35.878  19.981  -5.501  1.00 14.30           C
+ATOM    154  OE1 GLN A  22      35.529  19.359  -4.472  1.00 18.17           O
+ATOM    155  NE2 GLN A  22      36.611  21.049  -5.457  1.00 20.19           N
+ATOM    173  N   GLU A  25      40.443  19.021  -7.568  1.00 47.81           N
+ATOM    174  CA  GLU A  25      40.160  20.153  -6.631  1.00 42.14           C
+ATOM    175  C   GLU A  25      39.675  21.334  -7.469  1.00 36.95           C
+ATOM    176  O   GLU A  25      39.955  22.496  -7.223  1.00 45.86           O
+ATOM    177  CB  GLU A  25      41.460  20.522  -5.912  0.00 41.53           C
+ATOM    178  CG  GLU A  25      42.518  19.432  -5.947  0.00 41.54           C
+ATOM    179  CD  GLU A  25      43.919  19.983  -6.125  0.00 41.41           C
+ATOM    180  OE1 GLU A  25      44.311  20.257  -7.279  0.00 41.37           O
+ATOM    181  OE2 GLU A  25      44.629  20.145  -5.110  0.00 41.34           O
+ATOM    182  N   GLY A  25A     38.940  21.006  -8.547  1.00 35.16           N
+ATOM    183  CA  GLY A  25A     38.466  22.062  -9.438  1.00 30.84           C
+ATOM    184  C   GLY A  25A     37.135  22.649  -9.008  1.00 26.52           C
+ATOM    185  O   GLY A  25A     36.678  22.555  -7.874  1.00 24.92           O
+ATOM    191  N   VAL A  27      33.241  23.494  -8.492  1.00 13.48           N
+ATOM    192  CA  VAL A  27      32.081  22.686  -8.152  1.00 14.56           C
+ATOM    193  C   VAL A  27      30.869  23.108  -9.008  1.00 14.24           C
+ATOM    194  O   VAL A  27      30.552  24.277  -9.114  1.00 15.24           O
+ATOM    195  CB  VAL A  27      31.716  22.923  -6.662  1.00 13.92           C
+ATOM    196  CG1 VAL A  27      30.493  22.048  -6.301  1.00 19.33           C
+ATOM    197  CG2 VAL A  27      32.934  22.549  -5.797  1.00 15.63           C
+ATOM    749  N   SER A 103      34.060  24.070  -2.055  1.00 11.00           N
+ATOM    750  CA  SER A 103      34.714  22.821  -2.367  1.00 11.78           C
+ATOM    751  C   SER A 103      33.932  21.688  -1.675  1.00 11.18           C
+ATOM    752  O   SER A 103      33.182  21.910  -0.755  1.00 10.91           O
+ATOM    753  CB  SER A 103      36.150  22.787  -1.786  1.00 14.12           C
+ATOM    754  OG  SER A 103      36.858  21.681  -2.392  1.00 15.49           O
+ATOM    755  N   LEU A 104      34.225  20.450  -2.085  1.00 12.60           N
+ATOM    756  CA  LEU A 104      33.727  19.278  -1.440  1.00 13.49           C
+ATOM    757  C   LEU A 104      34.844  18.562  -0.681  1.00 14.51           C
+ATOM    758  O   LEU A 104      34.606  17.475  -0.159  1.00 16.82           O
+ATOM    759  CB  LEU A 104      33.019  18.367  -2.402  1.00 13.27           C
+ATOM    760  CG  LEU A 104      31.817  18.935  -3.135  1.00 13.24           C
+ATOM    761  CD1 LEU A 104      31.200  17.904  -4.040  1.00 16.51           C
+ATOM    762  CD2 LEU A 104      30.821  19.574  -2.182  1.00 22.63           C
+ATOM    763  N   LYS A 105      36.033  19.165  -0.594  1.00 14.92           N
+ATOM    764  CA  LYS A 105      37.139  18.617   0.145  1.00 16.30           C
+ATOM    765  C   LYS A 105      38.020  19.765   0.678  1.00 17.78           C
+ATOM    766  O   LYS A 105      37.916  20.873   0.200  1.00 20.07           O
+ATOM    767  CB  LYS A 105      37.943  17.627  -0.649  1.00 26.20           C
+ATOM    768  CG  LYS A 105      38.587  18.132  -1.914  1.00 38.16           C
+ATOM    769  CD  LYS A 105      38.295  17.264  -3.120  0.00 45.15           C
+ATOM    770  CE  LYS A 105      38.612  15.811  -2.887  1.00 51.58           C
+ATOM    771  NZ  LYS A 105      39.344  15.183  -4.045  1.00 74.91           N
+TER    1093      PRO A 151
+END
+""",
+   "A",
+   22,
+   "OE1",
+   [ (35.529, 19.359, -4.472)
+   ],
+   0.1,
+   []
+  ],
+
+  ["1xso_amide_flipped",
+   """\
+CRYST1   73.450   68.940   58.760  90.00  90.00  90.00 P 21 21 21    8
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.013615  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.014505  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.017018        0.00000
+MTRIX1   1  0.987830 -0.010290 -0.155220        1.43777    1
+MTRIX2   1 -0.008360 -0.999880  0.013070       33.15571    1
+MTRIX3   1 -0.155330 -0.011610 -0.987790       17.07902    1
+ATOM    924  N   GLY A 127      23.876  47.136  -0.439  1.00 18.38           N
+ATOM    925  CA  GLY A 127      22.692  47.748   0.105  1.00 17.97           C
+ATOM    926  C   GLY A 127      22.873  49.006   0.862  1.00 18.02           C
+ATOM    927  O   GLY A 127      21.908  49.714   1.225  1.00 24.20           O
+ATOM    932  N   ASN A 129      24.271  49.439   3.964  1.00 14.98           N
+ATOM    933  CA  ASN A 129      24.135  49.362   5.428  1.00 15.72           C
+ATOM    934  C   ASN A 129      24.971  48.198   5.914  1.00 13.04           C
+ATOM    935  O   ASN A 129      25.459  47.398   5.129  1.00 14.62           O
+ATOM    936  CB  ASN A 129      22.691  49.260   5.878  1.00 15.14           C
+ATOM    937  CG  ASN A 129      21.911  48.176   5.160  1.00 19.11           C
+ATOM    938  OD1 ASN A 129      22.371  47.040   5.114  1.00 15.87           O
+ATOM    939  ND2 ASN A 129      20.811  48.560   4.517  1.00 30.20           N
+ATOM    948  N   GLU A 131      24.275  45.568   7.628  1.00 14.32           N
+ATOM    949  CA  GLU A 131      23.665  44.242   7.404  1.00 15.35           C
+ATOM    950  C   GLU A 131      23.963  43.701   6.022  1.00 13.39           C
+ATOM    951  O   GLU A 131      24.150  42.499   5.830  1.00 13.14           O
+ATOM    952  CB  GLU A 131      22.201  44.278   7.753  1.00 18.42           C
+ATOM    953  CG  GLU A 131      21.474  42.964   7.698  1.00 23.30           C
+ATOM    954  CD  GLU A 131      21.905  42.014   8.811  1.00 27.07           C
+ATOM    955  OE1 GLU A 131      22.005  42.457   9.969  1.00 37.65           O
+ATOM    956  OE2 GLU A 131      22.106  40.834   8.513  1.00 32.51           O
+ATOM    957  N   SER A 132      24.016  44.584   5.024  1.00 13.86           N
+ATOM    958  CA  SER A 132      24.353  44.193   3.664  1.00 12.64           C
+ATOM    959  C   SER A 132      25.696  43.468   3.641  1.00 12.67           C
+ATOM    960  O   SER A 132      25.800  42.483   2.891  1.00 13.76           O
+ATOM    961  CB  SER A 132      24.320  45.391   2.737  1.00 14.43           C
+ATOM    962  OG  SER A 132      24.685  45.044   1.405  1.00 13.58           O
+ATOM    991  N   ASN A 137      20.989  41.281   3.329  1.00 12.56           N
+ATOM    992  CA  ASN A 137      19.607  41.651   3.597  1.00 13.15           C
+ATOM    993  C   ASN A 137      18.639  40.676   2.991  1.00 10.75           C
+ATOM    994  O   ASN A 137      17.511  41.040   2.611  1.00 15.89           O
+ATOM    995  CB  ASN A 137      19.280  43.097   3.159  1.00 13.87           C
+ATOM    996  CG  ASN A 137      19.949  44.125   4.040  1.00 16.63           C
+ATOM    997  OD1 ASN A 137      19.570  44.245   5.254  1.00 19.77           O
+ATOM    998  ND2 ASN A 137      20.916  44.887   3.583  1.00 16.13           N
+TER    1093      PRO A 151
+END
+""",
+   "A",
+   129,
+   "OD1",
+   [ (22.371, 47.040, 5.114)
+   ],
+   0.1,
+   []
+  ],
+
+  ["1xso_histidine_unflipped",
+   """\
+CRYST1   73.450   68.940   58.760  90.00  90.00  90.00 P 21 21 21    8
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.013615  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.014505  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.017018        0.00000
+MTRIX1   1  0.987830 -0.010290 -0.155220        1.43777    1
+MTRIX2   1 -0.008360 -0.999880  0.013070       33.15571    1
+MTRIX3   1 -0.155330 -0.011610 -0.987790       17.07902    1
+ATOM    103  N   VAL A  17      15.777  20.067  -3.917  1.00 11.00           N
+ATOM    104  CA  VAL A  17      16.936  19.211  -3.933  1.00 12.41           C
+ATOM    105  C   VAL A  17      18.055  19.934  -4.702  1.00  9.37           C
+ATOM    106  O   VAL A  17      17.836  20.513  -5.742  1.00 13.60           O
+ATOM    107  CB  VAL A  17      16.635  17.858  -4.617  1.00 16.16           C
+ATOM    108  CG1 VAL A  17      17.933  16.983  -4.557  1.00 21.33           C
+ATOM    109  CG2 VAL A  17      15.541  17.118  -3.839  1.00 19.65           C
+ATOM    117  N   HIS A  19      22.094  19.284  -6.044  1.00 11.05           N
+ATOM    118  CA  HIS A  19      23.225  18.356  -6.269  1.00 10.79           C
+ATOM    119  C   HIS A  19      24.527  19.153  -6.473  1.00 11.77           C
+ATOM    120  O   HIS A  19      24.547  20.205  -7.093  1.00 13.93           O
+ATOM    121  CB  HIS A  19      22.906  17.530  -7.527  1.00 18.16           C
+ATOM    122  CG  HIS A  19      21.858  16.494  -7.322  1.00 26.17           C
+ATOM    123  ND1 HIS A  19      22.149  15.200  -6.888  1.00 28.19           N
+ATOM    124  CD2 HIS A  19      20.525  16.560  -7.452  1.00 32.51           C
+ATOM    125  CE1 HIS A  19      21.056  14.541  -6.758  1.00 32.39           C
+ATOM    126  NE2 HIS A  19      20.028  15.310  -7.102  1.00 34.60           N
+HETATM 2433  O   HOH B 215      21.490  11.727  -5.002  1.00 34.66           O
+HETATM 2436  O   HOH B 218      16.653  14.508   9.188  1.00 35.04           O
+HETATM 2498  O   HOH B 280       8.420  -1.071  13.530  1.00 50.43           O
+HETATM 2525  O   HOH B 307      37.167  -8.451   4.794  1.00 59.87           O
+HETATM 2538  O   HOH B 320      14.429  18.062  16.293  1.00 67.01           O
+TER    1093      PRO A 151
+END
+""",
+   "A",
+   19,
+   "ND1",
+   [ (22.149, 15.200, -6.888)
+   ],
+   0.1,
+   []
+  ],
+
+  ["1xso_histidine_flipped",
+   """\
+CRYST1   73.450   68.940   58.760  90.00  90.00  90.00 P 21 21 21    8
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.013615  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.014505  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.017018        0.00000
+MTRIX1   1  0.987830 -0.010290 -0.155220        1.43777    1
+MTRIX2   1 -0.008360 -0.999880  0.013070       33.15571    1
+MTRIX3   1 -0.155330 -0.011610 -0.987790       17.07902    1
+ATOM   1115  N   VAL B   5      26.079  16.885  14.997  1.00  8.83           N
+ATOM   1116  CA  VAL B   5      24.686  16.588  15.374  1.00  8.78           C
+ATOM   1117  C   VAL B   5      23.969  16.005  14.172  1.00  9.25           C
+ATOM   1118  O   VAL B   5      24.335  16.257  13.044  1.00 10.86           O
+ATOM   1119  CB  VAL B   5      23.982  17.879  15.867  1.00 10.17           C
+ATOM   1120  CG1 VAL B   5      23.763  18.899  14.740  1.00 14.21           C
+ATOM   1121  CG2 VAL B   5      22.685  17.620  16.633  1.00 14.81           C
+ATOM   1196  N   VAL B  17      17.298  13.044  18.182  1.00 10.96           N
+ATOM   1197  CA  VAL B  17      18.457  13.879  18.007  1.00  9.04           C
+ATOM   1198  C   VAL B  17      19.658  13.118  18.600  1.00 11.44           C
+ATOM   1199  O   VAL B  17      19.557  12.560  19.692  1.00 14.98           O
+ATOM   1200  CB  VAL B  17      18.288  15.251  18.678  1.00 13.50           C
+ATOM   1201  CG1 VAL B  17      19.502  16.123  18.558  1.00 17.18           C
+ATOM   1202  CG2 VAL B  17      17.090  16.002  18.108  1.00 20.25           C
+ATOM   1210  N   HIS B  19      23.916  13.702  19.282  1.00  9.93           N
+ATOM   1211  CA  HIS B  19      25.062  14.561  19.425  1.00  9.29           C
+ATOM   1212  C   HIS B  19      26.366  13.749  19.441  1.00  8.50           C
+ATOM   1213  O   HIS B  19      26.409  12.668  20.031  1.00 12.76           O
+ATOM   1214  CB  HIS B  19      24.996  15.427  20.688  1.00 12.81           C
+ATOM   1215  CG  HIS B  19      23.810  16.312  20.701  1.00 19.55           C
+ATOM   1216  ND1 HIS B  19      23.897  17.622  20.251  1.00 24.11           N
+ATOM   1217  CD2 HIS B  19      22.556  16.137  21.208  1.00 25.44           C
+ATOM   1218  CE1 HIS B  19      22.723  18.181  20.377  1.00 27.18           C
+ATOM   1219  NE2 HIS B  19      21.897  17.308  20.959  1.00 29.89           N
+HETATM 2413  O   HOH B 195      10.371  -1.926  15.288  1.00 30.69           O
+HETATM 2486  O   HOH B 268      10.982   5.736  29.100  1.00 46.97           O
+HETATM 2528  O   HOH B 310      19.641  17.849  22.382  1.00 61.47           O
+TER    1093      PRO A 151
+END
+""",
+   "B",
+   19,
+   "ND1",
+   [ (22.666, 16.182, 21.033)
+   ],
+   0.1,
+   []
+  ],
+
+  ["1xso_histidine_large_unflip_preference",
+   """\
+CRYST1   73.450   68.940   58.760  90.00  90.00  90.00 P 21 21 21    8
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.013615  0.000000  0.000000        0.00000
+SCALE2      0.000000  0.014505  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.017018        0.00000
+MTRIX1   1  0.987830 -0.010290 -0.155220        1.43777    1
+MTRIX2   1 -0.008360 -0.999880  0.013070       33.15571    1
+MTRIX3   1 -0.155330 -0.011610 -0.987790       17.07902    1
+ATOM   1115  N   VAL B   5      26.079  16.885  14.997  1.00  8.83           N
+ATOM   1116  CA  VAL B   5      24.686  16.588  15.374  1.00  8.78           C
+ATOM   1117  C   VAL B   5      23.969  16.005  14.172  1.00  9.25           C
+ATOM   1118  O   VAL B   5      24.335  16.257  13.044  1.00 10.86           O
+ATOM   1119  CB  VAL B   5      23.982  17.879  15.867  1.00 10.17           C
+ATOM   1120  CG1 VAL B   5      23.763  18.899  14.740  1.00 14.21           C
+ATOM   1121  CG2 VAL B   5      22.685  17.620  16.633  1.00 14.81           C
+ATOM   1196  N   VAL B  17      17.298  13.044  18.182  1.00 10.96           N
+ATOM   1197  CA  VAL B  17      18.457  13.879  18.007  1.00  9.04           C
+ATOM   1198  C   VAL B  17      19.658  13.118  18.600  1.00 11.44           C
+ATOM   1199  O   VAL B  17      19.557  12.560  19.692  1.00 14.98           O
+ATOM   1200  CB  VAL B  17      18.288  15.251  18.678  1.00 13.50           C
+ATOM   1201  CG1 VAL B  17      19.502  16.123  18.558  1.00 17.18           C
+ATOM   1202  CG2 VAL B  17      17.090  16.002  18.108  1.00 20.25           C
+ATOM   1210  N   HIS B  19      23.916  13.702  19.282  1.00  9.93           N
+ATOM   1211  CA  HIS B  19      25.062  14.561  19.425  1.00  9.29           C
+ATOM   1212  C   HIS B  19      26.366  13.749  19.441  1.00  8.50           C
+ATOM   1213  O   HIS B  19      26.409  12.668  20.031  1.00 12.76           O
+ATOM   1214  CB  HIS B  19      24.996  15.427  20.688  1.00 12.81           C
+ATOM   1215  CG  HIS B  19      23.810  16.312  20.701  1.00 19.55           C
+ATOM   1216  ND1 HIS B  19      23.897  17.622  20.251  1.00 24.11           N
+ATOM   1217  CD2 HIS B  19      22.556  16.137  21.208  1.00 25.44           C
+ATOM   1218  CE1 HIS B  19      22.723  18.181  20.377  1.00 27.18           C
+ATOM   1219  NE2 HIS B  19      21.897  17.308  20.959  1.00 29.89           N
+HETATM 2413  O   HOH B 195      10.371  -1.926  15.288  1.00 30.69           O
+HETATM 2486  O   HOH B 268      10.982   5.736  29.100  1.00 46.97           O
+HETATM 2528  O   HOH B 310      19.641  17.849  22.382  1.00 61.47           O
+TER    1093      PRO A 151
+END
+""",
+   "B",
+   19,
+   "ND1",
+   [ (23.897, 17.622, 20.251)
+   ],
+   0.1,
+   ["non_flip_preference=1000"]
+  ],
+
+["1ehz_aromatic_methyl_rotator",
+   """\
+CRYST1   54.981   33.389   61.921  90.00  90.20  90.00 P 1 21 1      2
+ORIGX1      1.000000  0.000000  0.000000        0.00000
+ORIGX2      0.000000  1.000000  0.000000        0.00000
+ORIGX3      0.000000  0.000000  1.000000        0.00000
+SCALE1      0.018188  0.000000  0.000063        0.00000
+SCALE2      0.000000  0.029950  0.000000        0.00000
+SCALE3      0.000000  0.000000  0.016150        0.00000
+ATOM    133  P     U A   7      67.463  47.074  35.969  1.00 44.37           P
+ATOM    134  OP1   U A   7      68.318  46.756  34.822  1.00 48.09           O
+ATOM    135  OP2   U A   7      67.945  47.948  37.077  1.00 45.68           O
+ATOM    136  O5'   U A   7      66.104  47.724  35.455  1.00 40.88           O
+ATOM    137  C5'   U A   7      65.285  47.024  34.459  1.00 37.89           C
+ATOM    138  C4'   U A   7      64.055  47.852  34.101  1.00 35.74           C
+ATOM    139  O4'   U A   7      63.297  48.107  35.326  1.00 38.13           O
+ATOM    140  C3'   U A   7      64.317  49.197  33.459  1.00 36.87           C
+ATOM    141  O3'   U A   7      63.402  49.394  32.378  1.00 37.45           O
+ATOM    142  C2'   U A   7      64.097  50.171  34.624  1.00 36.55           C
+ATOM    143  O2'   U A   7      63.595  51.417  34.246  1.00 35.54           O
+ATOM    144  C1'   U A   7      63.015  49.475  35.442  1.00 37.23           C
+ATOM    145  N1    U A   7      63.056  49.858  36.864  1.00 36.91           N
+ATOM    146  C2    U A   7      62.011  50.628  37.343  1.00 34.52           C
+ATOM    147  O2    U A   7      61.087  50.966  36.653  1.00 34.65           O
+ATOM    148  N3    U A   7      62.112  50.993  38.659  1.00 37.03           N
+ATOM    149  C4    U A   7      63.131  50.684  39.541  1.00 40.15           C
+ATOM    150  O4    U A   7      63.105  51.143  40.699  1.00 36.62           O
+ATOM    151  C5    U A   7      64.179  49.865  38.971  1.00 36.52           C
+ATOM    152  C6    U A   7      64.106  49.490  37.691  1.00 36.25           C
+HETATM 1061  P   5MC A  49      65.638  54.464  33.461  1.00 31.61           P
+HETATM 1062  OP1 5MC A  49      66.252  55.251  34.624  1.00 26.54           O
+HETATM 1063  OP2 5MC A  49      65.643  52.989  33.540  1.00 33.43           O
+HETATM 1064  O5' 5MC A  49      64.126  54.972  33.441  1.00 31.88           O
+HETATM 1065  C5' 5MC A  49      63.204  54.605  32.392  1.00 32.68           C
+HETATM 1066  C4' 5MC A  49      61.796  55.006  32.810  1.00 30.80           C
+HETATM 1067  O4' 5MC A  49      61.292  54.110  33.848  1.00 32.40           O
+HETATM 1068  C3' 5MC A  49      61.745  56.381  33.468  1.00 31.44           C
+HETATM 1069  O3' 5MC A  49      61.666  57.398  32.483  1.00 31.25           O
+HETATM 1070  C2' 5MC A  49      60.485  56.287  34.312  1.00 34.62           C
+HETATM 1071  O2' 5MC A  49      59.414  56.469  33.417  1.00 35.41           O
+HETATM 1072  C1' 5MC A  49      60.568  54.849  34.825  1.00 33.10           C
+HETATM 1073  N1  5MC A  49      61.331  54.787  36.102  1.00 32.34           N
+HETATM 1074  C2  5MC A  49      60.737  55.340  37.271  1.00 32.42           C
+HETATM 1075  O2  5MC A  49      59.568  55.832  37.197  1.00 28.24           O
+HETATM 1076  N3  5MC A  49      61.420  55.348  38.428  1.00 29.50           N
+HETATM 1077  C4  5MC A  49      62.626  54.785  38.503  1.00 32.97           C
+HETATM 1078  N4  5MC A  49      63.188  54.744  39.702  1.00 28.75           N
+HETATM 1079  C5  5MC A  49      63.270  54.215  37.358  1.00 34.81           C
+HETATM 1080  C6  5MC A  49      62.583  54.243  36.161  1.00 33.69           C
+HETATM 1081  CM5 5MC A  49      64.655  53.581  37.465  1.00 34.52           C
+HETATM 1692  O   HOH A 131      65.908  53.779  40.692  1.00 47.68           O
+HETATM 1783  O   HOH A 675      68.361  55.414  38.613  1.00 58.70           O
+TER    1653        A A  76
+END
+""",
+   "A",
+   49,
+   "HM51",
+   [ (64.602, 52.635, 37.674), (65.188, 53.997, 38.160), (65.153, 53.663, 36.637)
+   ],
+   0.1,
+   []
+  ]
+
+]
+
+def RunRegressionTests():
+  for tc in testCases:
+    name = tc[0]
+    pdb_raw = tc[1]
+    chain = tc[2]
+    resID = tc[3]
+    atomName = tc[4]
+    positions = tc[5]
+    maxDist = tc[6]
+    extraArgs = tc[7]
+
+    print('Testing regression on', name)
+
+    pdb_file = "./deleteme.pdb"
+    with open(pdb_file, "w") as f:
+      f.write(pdb_raw)
+    if (op.exists("./deletemeFH.pdb")):
+      os.remove("./deletemeFH.pdb")
+    if (op.exists("./deleteme_description.txt")):
+      os.remove("./deleteme_description.txt")
+    out = StringIO()
+    try:
+      # Run the program
+      args = [pdb_file, "add_flip_movers=True", "output.description_file=./deleteme_description.txt"]
+      args.extend(extraArgs)
+      results = run_program(program_class=reduce2.Program, logger=out, args=args)
+      # Check the position of the atom to see if it is near enough to one of the expected locations.
+      found = False
+      for c in results.model.chains():
+        if c.id == chain:
+          for rg in c.residue_groups():
+            if rg.resseq_as_int() == resID:
+              for atom in rg.atoms():
+                if atom.name.strip().upper() == atomName:
+                  found = True
+                  loc = atom.xyz
+                  closeEnough = False
+                  for pos in positions:
+                    dist = math.sqrt( (loc[0]-pos[0])*(loc[0]-pos[0]) +
+                      (loc[1]-pos[1])*(loc[1]-pos[1]) +
+                      (loc[2]-pos[2])*(loc[2]-pos[2])
+                    )
+                    if dist <= maxDist:
+                      closeEnough = True
+                  if not closeEnough:
+                    return "Atom "+chain+" "+str(resID)+" "+atomName+" in "+name+" too far from expected locations"
+      if not found:
+        return "Did not find atom "+atomName+" in chain "+chain+" residue "+str(resID)+" of "+name
+    except Exception as e:
+      return "Exception when running reduce2: "+str(e)
+    if (op.exists(pdb_file)):
+      os.remove(pdb_file)
+    if (op.exists("./deletemeFH.pdb")):
+      os.remove("./deletemeFH.pdb")
+    if (op.exists("./deleteme_description.txt")):
+      os.remove("./deleteme_description.txt")
+
+  return ""
+
 if __name__ == '__main__':
 
   ret = RunReduceTests()
@@ -56,3 +510,11 @@ def RunReduceTests():
     print(ret)
 
   assert (len(ret) == 0)
+
+  ret = RunRegressionTests()
+  if len(ret) == 0:
+    print('Success!')
+  else:
+    print(ret)
+
+  assert (len(ret) == 0)
diff --git a/mmtbx/solvent/ordered_solvent.py b/mmtbx/solvent/ordered_solvent.py
index 751469952a..257d723f72 100644
--- a/mmtbx/solvent/ordered_solvent.py
+++ b/mmtbx/solvent/ordered_solvent.py
@@ -101,6 +101,9 @@
             after ferst few macro-cycles, filter_only - remove water only, \
             every_macro_cycle - do water update every macro-cycle
     .short_caption = Mode
+  keep_existing = False
+    .type = bool
+    .help = Keep existing in the model water
   n_cycles = 3
     .type = int
     .short_caption = Number of cycles
@@ -235,6 +238,9 @@ def __init__(self, fmodel,
     self._peaks     = None
     self.n_water    = None
     self.model_size_init = self.model.size()
+    self.existing_solvent = None
+    if(self.params.keep_existing):
+      self.existing_solvent = self.model.solvent_selection().iselection()
     #
     self._call(msg="Start")
     self._call(msg="Compute maps",     func=self._get_maps)
@@ -277,7 +283,7 @@ def _get_maps(self):
       map_2_type = self.params.secondary_map_and_map_cc_filter.cc_map_2_type)
 
   def _filter_solvent(self):
-    sol_sel   = self.model.solvent_selection()
+    sol_sel   = self.model.solvent_selection(offset = self.existing_solvent)
     hd_sel    = self.model.get_hd_selection()
     n_sol_start = self.n_water
     mfp = self.params.secondary_map_and_map_cc_filter
@@ -310,6 +316,8 @@ def _filter_solvent(self):
             for atom in r.atoms():
               if(atom.element_is_hydrogen()): continue
               i_seq = atom.i_seq
+              if(self.params.keep_existing and i_seq in self.existing_solvent):
+                continue
               if(not distance_selection[i_seq]): keep = False
               if(atom.occ > self.params.occupancy_max or
                  atom.occ < self.params.occupancy_min): keep = False
@@ -321,7 +329,8 @@ def _filter_solvent(self):
                 atom      = atom,
                 min_cc    = mfp.poor_cc_threshold,
                 min_value = mfp.poor_map_value_threshold)
-              if(not good_map): keep = False
+              if(not good_map):
+                keep = False
               #
               # XXX It does not work for 1f8t. WHY? XXX
               #
@@ -412,7 +421,7 @@ def refine_adp(self):
        ):
       hd_sel     = self.model.get_hd_selection()
       not_hd_sel = ~hd_sel
-      sol_sel    = self.model.solvent_selection()
+      sol_sel    = self.model.solvent_selection(offset=self.existing_solvent)
       not_sol_sel= ~sol_sel
       selection_aniso = self.model.get_xray_structure().use_u_aniso().deep_copy()
       if(self.params.new_solvent == "anisotropic"):
@@ -447,7 +456,8 @@ def refine_occupancies(self):
        ):
       self.fmodels.fmodel_xray().xray_structure.scatterers().flags_set_grads(
         state = False)
-      i_selection = self.model.solvent_selection().iselection()
+      i_selection = self.model.solvent_selection(
+        offset=self.existing_solvent).iselection()
       self.fmodels.fmodel_xray().xray_structure.scatterers(
         ).flags_set_grad_occupancy(iselection = i_selection)
       lbfgs_termination_params = scitbx.lbfgs.termination_parameters(
@@ -467,7 +477,7 @@ def refine_occupancies(self):
   def _correct_drifted_waters(self):
     if(self.params.mode != "filter_only"): return
     if(not self.params.correct_drifted_waters): return
-    sol_sel = self.model.solvent_selection()
+    sol_sel = self.model.solvent_selection(offset=self.existing_solvent)
     hd_sel  = self.model.get_hd_selection()
     hd_sol = sol_sel & hd_sel
     if(hd_sol.count(True)>0): return
diff --git a/scitbx/matrix/__init__.py b/scitbx/matrix/__init__.py
index 6c6a8ccd76..12d4fb4998 100644
--- a/scitbx/matrix/__init__.py
+++ b/scitbx/matrix/__init__.py
@@ -39,7 +39,7 @@ class rec(object):
 
   Examples
   --------
-  >>> from scitbx.matix import rec
+  >>> from scitbx.matrix import rec
   >>> m = rec((1, 2, 3, 4), (2, 2))
   >>> print m.trace()
   5
diff --git a/simtbx/command_line/hopper.py b/simtbx/command_line/hopper.py
index 89f396454d..cfd0e8ac75 100644
--- a/simtbx/command_line/hopper.py
+++ b/simtbx/command_line/hopper.py
@@ -126,6 +126,7 @@ def run(self):
             refl_names_already = COMM.bcast(refl_names_already)
 
         exp_gatheredRef_spec = []  # optional list of expt, refls, spectra
+        trefs = []
         for i_exp, line in enumerate(input_lines):
             if i_exp == self.params.max_process:
                 break
@@ -248,10 +249,29 @@ def run(self):
             if SIM.refining_Fhkl:
                 nparam += SIM.Num_ASU*SIM.num_Fhkl_channels
             x0 = [1] * nparam
+            tref = time.time()
+            MAIN_LOGGER.info("Beginning refinement of shot %d / %d" % (i_exp+1, len(input_lines)))
             try:
                 x = Modeler.Minimize(x0, SIM, i_exp=i_exp)
             except StopIteration:
                 x = Modeler.target.x0
+            tref = time.time()-tref
+            sigz = niter = None
+            try:
+                niter = len(Modeler.target.all_hop_id)
+                sigz = np.mean(Modeler.target.all_sigZ)
+            except Exception:
+                pass
+
+            trefs.append(tref)
+            print_s = "Finished refinement of shot %d / %d in %.4f sec. (rank mean t/im=%.4f sec.)" \
+                        % (i_exp+1, len(input_lines), tref, np.mean(trefs))
+            if sigz is not None and niter is not None:
+                print_s += " Ran %d iterations. Final sigmaZ = %.1f," % (niter, sigz)
+            if COMM.rank==0:
+                MAIN_LOGGER.info(print_s)
+            else:
+                MAIN_LOGGER.debug(print_s)
             if self.params.profile:
                 SIM.D.show_timings(COMM.rank)
 
diff --git a/simtbx/command_line/integrate.py b/simtbx/command_line/integrate.py
index f910c5dbed..c23e6011e9 100644
--- a/simtbx/command_line/integrate.py
+++ b/simtbx/command_line/integrate.py
@@ -437,24 +437,25 @@ def df_iter():
         utils.refls_to_hkl(Rindexed, data_expt.detector, data_expt.beam, data_expt.crystal, update_table=True)
         try:
             int_expt, int_refl = integrate(args.procPhil, data_exptList, Rindexed, pred)
-        except RuntimeError:
-            print("Integration failed for %s" % pred_file)
-            continue
-        int_expt_name = "%s/%s_%d_integrated.expt" % ( args.outdir,tag, i_f)
-        int_expt.as_file(int_expt_name)
-        int_refl['bbox'] = int_refl['shoebox'].bounding_boxes()
-        int_refl_name = int_expt_name.replace(".expt", ".refl")
-        int_refl.as_file(int_refl_name)
+            int_expt_name = "%s/%s_%d_integrated.expt" % ( args.outdir,tag, i_f)
+            int_expt.as_file(int_expt_name)
+            int_refl['bbox'] = int_refl['shoebox'].bounding_boxes()
+            int_refl_name = int_expt_name.replace(".expt", ".refl")
+            int_refl.as_file(int_refl_name)
+        except RuntimeError as err:
+            print("Integration failed for %s because: '%s'" % (pred_file, str(err)))
+
         all_dfs.append(df)
         all_pred_names.append(pred_file)
         spec_name = df.spectrum_filename.values[0]
         if spec_name is None:
             spec_name = ""
-        exp_ref_spec_lines.append("%s %s %s\n" % (new_expt_name, int_refl_name, spec_name))
+        exp_ref_spec_lines.append("%s %s %s\n" % (new_expt_name, pred_file, spec_name))
 
     if all_dfs:
         all_dfs = pandas.concat(all_dfs)
         all_dfs["predicted_refls"] = all_pred_names
+        all_dfs["predictions"] = all_pred_names
     else:
         all_dfs = None
 
diff --git a/simtbx/diffBragg/hopper_io.py b/simtbx/diffBragg/hopper_io.py
index 5f833760f1..da9316e37e 100644
--- a/simtbx/diffBragg/hopper_io.py
+++ b/simtbx/diffBragg/hopper_io.py
@@ -83,6 +83,10 @@ def save_to_pandas(x, Mod, SIM, orig_exp_name, params, expt, rank_exp_idx, stg1_
         diff_gam_b = diff_gam_c = diff_gam_a
     if params.isotropic.diffuse_sigma:
         diff_sig_b = diff_sig_c = diff_sig_a
+
+    if params.simulator.crystal.has_isotropic_ncells:
+        Nb = Nc = Na
+
     eta_a, eta_b, eta_c = hopper_utils.get_mosaicity_from_x(x, Mod, SIM)
     a_init, b_init, c_init, al_init, be_init, ga_init = SIM.crystal.dxtbx_crystal.get_unit_cell().parameters()
 
diff --git a/simtbx/diffBragg/hopper_utils.py b/simtbx/diffBragg/hopper_utils.py
index 8ac7709cab..8d44789a8b 100644
--- a/simtbx/diffBragg/hopper_utils.py
+++ b/simtbx/diffBragg/hopper_utils.py
@@ -150,7 +150,6 @@ def __init__(self, params):
         self.pids=None  # panel id (per spot)
         self.tilt_abc=None  # background plane constants (per spot), a,b are fast,slow scan components, c is offset
         self.selection_flags=None  # whether the spot was selected for refinement (sometimes poorly conditioned spots are rejected)
-        self.background=None  # background for entire image (same shape as the detector)
         self.tilt_cov = None  # covariance estimates from background fitting (not used)
         self.simple_weights = None  # not used
         self.refls_idx = None  # position of modeled spot in original refl array
@@ -432,7 +431,7 @@ def GatherFromExperiment(self, exp, ref, remove_duplicate_hkl=True, sg_symbol=No
         if roi_packet is None:
             return False
 
-        self.rois, self.pids, self.tilt_abc, self.selection_flags, self.background, self.tilt_cov = roi_packet
+        self.rois, self.pids, self.tilt_abc, self.selection_flags, background, self.tilt_cov = roi_packet
 
         if remove_duplicate_hkl and not self.no_rlp_info:
             is_not_a_duplicate = ~self.is_duplicate_hkl(refls)
@@ -475,7 +474,7 @@ def GatherFromExperiment(self, exp, ref, remove_duplicate_hkl=True, sg_symbol=No
         if not self.no_rlp_info:
             self.Q = [np.linalg.norm(refls[i_roi]["rlp"]) for i_roi in range(len(refls)) if self.selection_flags[i_roi]]
 
-        self.data_to_one_dim(img_data, is_trusted, self.background)
+        self.data_to_one_dim(img_data, is_trusted, background)
         return True
 
     def data_to_one_dim(self, img_data, is_trusted, background):
@@ -1404,6 +1403,9 @@ def model(x, Mod, SIM,  compute_grad=True, dont_rescale_gradient=False, update_s
 
     Ndef_params = [Mod.P["Ndef%d" % (i_n,)] for i_n in range(3)]
     Nd, Ne, Nf = [n_param.get_val(x[n_param.xpos]) for n_param in Ndef_params]
+    if SIM.D.isotropic_ncells:
+        Ne = Nd
+        Nf = Nd
     SIM.D.Ncells_def = Nd, Ne, Nf
 
     # diffuse signals
diff --git a/simtbx/diffBragg/phil.py b/simtbx/diffBragg/phil.py
index f1a0f837b7..806780f47c 100644
--- a/simtbx/diffBragg/phil.py
+++ b/simtbx/diffBragg/phil.py
@@ -1117,7 +1117,7 @@
     .type = str
     .help = panel list for refinement as a string, e.g. "0-8,10,32-40" . The ranges are inclusive,
     .help = e.g. 0-8 means panels 0,1,2,3,4,5,6,7,8
-  fit_tilt_using_weights = True
+  fit_tilt_using_weights = False
     .type = bool
     .help = if not using robust estimation for background, and instead using tilt plane fit,
     .help = then this parameter will toggle the use of weights. Weights are the estimated
diff --git a/simtbx/diffBragg/refiners/stage_two_refiner.py b/simtbx/diffBragg/refiners/stage_two_refiner.py
index 974daf81a2..b877c40f70 100644
--- a/simtbx/diffBragg/refiners/stage_two_refiner.py
+++ b/simtbx/diffBragg/refiners/stage_two_refiner.py
@@ -797,7 +797,7 @@ def _save_Zscore_data(self):
         if not os.path.exists(outdir):
             os.makedirs(outdir)
         fname = os.path.join(outdir, "sigZ_iter%d_rank%d" % (self.iterations, self.rank))
-        np.save(fname, self._shot_Zscores)
+        np.save(fname, np.array(self._shot_Zscores, object))
 
     def _sanity_check_grad(self):
         pass
diff --git a/simtbx/diffBragg/utils.py b/simtbx/diffBragg/utils.py
index 283fdb5447..3692633de0 100644
--- a/simtbx/diffBragg/utils.py
+++ b/simtbx/diffBragg/utils.py
@@ -641,7 +641,7 @@ def fit_plane_equation_to_background_pixels(shoebox_img, fit_sel, sigma_rdout=3,
     try:
         AWA_inv = np.linalg.inv(AWA)
     except np.linalg.LinAlgError:
-        MAIN_LOGGER.warning("WARNING: Fit did not work.. investigate reflection, number of pixels used in fit=%d" % len(fast))
+        MAIN_LOGGER.debug("WARNING: Fit did not work.. investigate reflection, number of pixels used in fit=%d" % len(fast))
         return None
     AtW = np.dot(A.T, W)
     t1, t2, t3 = np.dot(np.dot(AWA_inv, AtW), rho_bg)
@@ -650,7 +650,8 @@ def fit_plane_equation_to_background_pixels(shoebox_img, fit_sel, sigma_rdout=3,
     # vector of residuals
     r = rho_bg - np.dot(A, (t1, t2, t3))
     Nbg = len(rho_bg)
-    r_fact = np.dot(r.T, np.dot(W, r)) / (Nbg - 3)  # 3 parameters fit
+    with np.errstate(invalid='ignore'):
+        r_fact = np.dot(r.T, np.dot(W, r)) / (Nbg - 3)  # 3 parameters fit
     var_covar = AWA_inv * r_fact  # TODO: check for correlations in the off diagonal elems
 
     return (t1, t2, t3), var_covar
diff --git a/simtbx/nanoBragg/utils.py b/simtbx/nanoBragg/utils.py
index e160bde818..c901eea5d1 100644
--- a/simtbx/nanoBragg/utils.py
+++ b/simtbx/nanoBragg/utils.py
@@ -160,7 +160,7 @@ def flexBeam_sim_colors(CRYSTAL, DETECTOR, BEAM, Famp, energies, fluxes,
 
 
 def sim_background(DETECTOR, BEAM, wavelengths, wavelength_weights, total_flux, pidx=0, beam_size_mm=0.001,
-                   Fbg_vs_stol=None, sample_thick_mm=100, density_gcm3=1, molecular_weight=18):
+                   Fbg_vs_stol=None, sample_thick_mm=100, density_gcm3=1, molecular_weight=18, roi=None):
   """
   :param DETECTOR:
   :param BEAM: see sim_spots
@@ -173,6 +173,9 @@ def sim_background(DETECTOR, BEAM, wavelengths, wavelength_weights, total_flux,
   :param sample_thick_mm: path length of background that is exposed by the beam
   :param density_gcm3: density of background  (defaults to water)
   :param molecular_weight: molecular weight of background (defaults to water)
+  :param roi: region of interst on detector ((fast_min, fast_max), (slow_min, slow_max)),
+              e.g., ((10,20), (30,35)) would simulate in the 10x5 pixel region of interest
+              starting at pixel 10,30 (fast, slow). See also region_of_interest in nanoBragg_ext.cpp
   :return: raw_pixels as flex array, these can be passed to sim_spots function below
   """
   wavelength_weights = np.array(wavelength_weights)
@@ -192,6 +195,8 @@ def sim_background(DETECTOR, BEAM, wavelengths, wavelength_weights, total_flux,
   SIM.amorphous_density_gcm3 = density_gcm3
   SIM.amorphous_molecular_weight_Da = molecular_weight
   SIM.progress_meter = False
+  if roi is not None:
+    SIM.region_of_interest=roi
   SIM.add_background()
   background_raw_pixels = SIM.raw_pixels.deep_copy()
   SIM.free_all()
diff --git a/xfel/command_line/cspad_detector_congruence.py b/xfel/command_line/cspad_detector_congruence.py
index 9de82127e8..5196db4354 100644
--- a/xfel/command_line/cspad_detector_congruence.py
+++ b/xfel/command_line/cspad_detector_congruence.py
@@ -170,10 +170,7 @@ def detector_plot_dict(params, detector, data, title, units_str, show=True, reve
   # initialize the color map
   values = flex.double(list(data.values()))
   norm = Normalize(vmin=flex.min(values), vmax=flex.max(values))
-  if reverse_colormap:
-    cmap = plt.cm.get_cmap(params.colormap + "_r")
-  else:
-    cmap = plt.cm.get_cmap(params.colormap)
+  cmap = plt.get_cmap(params.colormap + ("_r" if reverse_colormap else ''))
   sm = cm.ScalarMappable(norm=norm, cmap=cmap)
   if len(values) == 0:
     print("no values")
diff --git a/xfel/command_line/cspad_detector_congruence2.py b/xfel/command_line/cspad_detector_congruence2.py
index cdde743624..5782eb4148 100644
--- a/xfel/command_line/cspad_detector_congruence2.py
+++ b/xfel/command_line/cspad_detector_congruence2.py
@@ -836,10 +836,7 @@ def detector_plot_dict(self, detector, data, title, units_str, show=True, revers
     # initialize the color map
     values = flex.double(data.values())
     norm = Normalize(vmin=flex.min(values), vmax=flex.max(values))
-    if reverse_colormap:
-      cmap = plt.cm.get_cmap(self.params.colormap + "_r")
-    else:
-      cmap = plt.cm.get_cmap(self.params.colormap)
+    cmap = plt.get_cmap(self.params.colormap + ("_r" if reverse_colormap else ''))
     sm = cm.ScalarMappable(norm=norm, cmap=cmap)
     if len(values) == 0:
       print("no values")
diff --git a/xfel/command_line/detector_residuals.py b/xfel/command_line/detector_residuals.py
index d9072390a1..8ca4533c46 100644
--- a/xfel/command_line/detector_residuals.py
+++ b/xfel/command_line/detector_residuals.py
@@ -473,7 +473,7 @@ def get_normalized_colors(self, data, vmin=None, vmax=None):
 
     # initialize the color map
     norm = Normalize(vmin=vmin, vmax=vmax)
-    cmap = plt.cm.get_cmap(self.params.colormap)
+    cmap = plt.get_cmap(self.params.colormap)
     sm = cm.ScalarMappable(norm=norm, cmap=cmap)
     color_vals = np.linspace(vmin, vmax, 11)
     sm.set_array(color_vals) # needed for colorbar
diff --git a/xfel/metrology/panel_fitting.py b/xfel/metrology/panel_fitting.py
index 1ba594595e..835ee34383 100644
--- a/xfel/metrology/panel_fitting.py
+++ b/xfel/metrology/panel_fitting.py
@@ -90,8 +90,8 @@ def plot_contours(self, ax, show=False):
       mahal_emp_cov = mahal_emp_cov.reshape(xx.shape)
       emp_cov_contour = ax.contour(xx, yy, np.sqrt(mahal_emp_cov),
                                   levels=[1.,2.,3.,4.,5.],
-                                  #cmap=plt.cm.PuBu_r,
-                                  cmap=plt.cm.cool_r,
+                                  #cmap=plt.get_cmap('PuBu_r'),
+                                  cmap=plt.get_cmap('cool_r'),
                                   linestyles='dashed')
 
     COV = self.rob_cov
@@ -111,8 +111,8 @@ def plot_contours(self, ax, show=False):
       mahal_robust_cov = mahal_robust_cov.reshape(xx.shape)
       robust_contour = ax.contour(xx, yy, np.sqrt(mahal_robust_cov),
                                  levels=[1.,2.,3.,4.,5.],
-                                 #cmap=plt.cm.YlOrBr_r,
-                                 cmap=plt.cm.spring_r,
+                                 #cmap=plt.get_cmap('YlOrBr_r'),
+                                 cmap=plt.get_cmap('spring_r'),
                                  linestyles='dotted')
 
 class three_feature_fit(Panel_MCD_Filter):


Test patch:
diff --git a/iotbx/run_tests.py b/iotbx/run_tests.py
index b45a0d359a..6f57a069a0 100644
--- a/iotbx/run_tests.py
+++ b/iotbx/run_tests.py
@@ -145,6 +145,7 @@
   "$D/regression/tst_group_rounding.py",
   "$D/regression/tst_hierarchy_occupancies_rounding.py",
   "$D/regression/tst_hierarchy_long_chain_ids_1.py",
+  "$D/regression/tst_hierarchy_long_resname_1.py",
   ]
 
 # failing tests on Windows, Python 2.7
@@ -152,7 +153,9 @@
   "$D/detectors/tst_debug_write.py",
 ]
 
-tst_list_fail = ["$D/regression/ncs/tst_ncs_reordered_chains.py",
+tst_list_fail = [
+  "$D/regression/ncs/tst_ncs_reordered_chains.py",
+  "$D/regression/tst_hierarchy_copy_select.py",
   ]
 if sys.platform == 'win32':
   tst_list_fail += tst_list_windows_fail