text: ! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  ! H0 X
  ! H0 X   libAtoms+QUIP: atomistic simulation library
  ! H0 X
  ! H0 X   Portions of this code were written by
  ! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
  ! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
  ! H0 X
  ! H0 X   Copyright 2006-2010.
  ! H0 X
  ! H0 X   These portions of the source code are released under the GNU General
  ! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
  ! H0 X
  ! H0 X   If you would like to license the source code under different terms,
  ! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
  ! H0 X
  ! H0 X   Portions of this code were written by Noam Bernstein as part of
  ! H0 X   his employment for the U.S. Government, and are not subject
  ! H0 X   to copyright in the USA.
  ! H0 X
  ! H0 X
  ! H0 X   When using this software, please cite the following reference:
  ! H0 X
  ! H0 X   http://www.libatoms.org
  ! H0 X
  ! H0 X  Additional contributions by
  ! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
  ! H0 X
  ! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  
  #include "error.inc"
  module clusters_module
    use error_module
    use system_module
    use linearalgebra_module
    use Table_module
    use periodictable_module
    use Atoms_types_module
    use Atoms_module
    use dictionary_module
    use ParamReader_module
    use group_module
    use constraints_module
    use dynamicalsystem_module
    use cinoutput_module
  
  #ifdef _MPI
  #ifndef _OLDMPI
    use mpi
  #endif
  #endif
  implicit none
  #ifdef _MPI
  #ifdef _OLDMPI
  include 'mpif.h'
  #endif
  #endif
  private
  
  integer, parameter :: &
     HYBRID_ACTIVE_MARK = 1, &
     HYBRID_BUFFER_MARK = 2, &
     HYBRID_TRANS_MARK = 3, &
     HYBRID_TERM_MARK = 4, &
     HYBRID_BUFFER_OUTER_LAYER_MARK = 5, &
     HYBRID_ELECTROSTATIC_MARK = 6, &
     HYBRID_NO_MARK = 0
  !   HYBRID_FIT_MARK = 6, &
  
  public :: HYBRID_ACTIVE_MARK, HYBRID_BUFFER_MARK, HYBRID_TRANS_MARK, HYBRID_TERM_MARK, &
     HYBRID_BUFFER_OUTER_LAYER_MARK, HYBRID_ELECTROSTATIC_MARK, HYBRID_NO_MARK
  ! HYBRID_FIT_MARK
  
  integer, parameter :: MAX_CUT_BONDS = 6
  public :: MAX_CUT_BONDS
  
  character(len=TABLE_STRING_LENGTH), parameter :: hybrid_mark_name(0:6) = &
    (/ "h_none    ", &
       "h_active  ", &
       "h_buffer  ", &
       "h_trans   ", &
       "h_term    ", &
       "h_outer_l ", &
       "h_fit     " /)
  
  public :: create_cluster_info_from_mark, carve_cluster, create_cluster_simple, create_hybrid_weights, &
      bfs_grow, bfs_step, multiple_images, discard_non_min_images, make_convex, create_embed_and_fit_lists, &
      create_embed_and_fit_lists_from_cluster_mark, &
      add_cut_hydrogens, construct_hysteretic_region, &
      create_pos_or_list_centred_hybrid_region, get_hybrid_list, &
      estimate_origin_extent
      !, construct_region, select_hysteretic_quantum_region
  
  !% Grow a selection list by bond hopping.
  interface bfs_grow
     module procedure bfs_grow_single
     module procedure bfs_grow_list
  end interface
  
  contains
  
    !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    !X
    !X   Cluster carving routines
    !X
    !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  
    !% On exit, 'list' will contain 'atom' (with shift '000') 
    !% plus the atoms within 'n' bonds hops of it.
    subroutine bfs_grow_single(this, list, atom, n, nneighb_only, min_images_only, alt_connect)
      type(Atoms), intent(in)  :: this
      type(Table), intent(out) :: list
      integer, intent(in) :: atom, n
      logical, optional, intent(in)::nneighb_only, min_images_only
      type(Connection), intent(in), optional :: alt_connect
  
      call append(list, (/atom, 0,0,0/)) ! Add atom with shift 000
      call bfs_grow_list(this, list, n, nneighb_only, min_images_only, alt_connect)
  
    end subroutine bfs_grow_single
  
  
    !% On exit, 'list' will have been grown by 'n' bond hops.
    subroutine bfs_grow_list(this, list, n, nneighb_only, min_images_only, alt_connect)
      type(Atoms), intent(in)   ::this
      type(Table), intent(inout)::list
      integer, intent(in) :: n
      logical, optional, intent(in)::nneighb_only, min_images_only
      type(Connection), intent(in), optional :: alt_connect
  
      type(Table)::tmplist
      integer::i
  
      do i=1,n
         call bfs_step(this, list, tmplist, nneighb_only, min_images_only, alt_connect=alt_connect)
         call append(list, tmplist)
      end do
  
      if (n >= 1) call finalise(tmplist)
    end subroutine bfs_grow_list
  
  
    !% Execute one Breadth-First-Search move on the atomic connectivity graph.
    subroutine bfs_step(this,input,output,nneighb_only, min_images_only, max_r, alt_connect, property, debugfile, error)
      type(Atoms),        intent(in), target      :: this  !% The atoms structure to perform the step on.
      type(Table),        intent(in)      :: input !% Table with intsize 4. First integer column is indices of atoms
                                                   !% already in the region, next 3 are shifts.
      type(Table),        intent(out)     :: output !% Table with intsize 4, containing the new atomic 
                                                    !% indices and shifts.
  
  
      logical, optional,  intent(in)      :: nneighb_only
      !% If present and true, sets whether only neighbours
      !% within the sum of the two respective covalent radii (multiplied by the atom's nneightol) are included,
      !% irrespective of the cutoff in the atoms structure
      !% (default is true).
  
      logical, optional, intent(in)       :: min_images_only 
      !% If true, there will be no repeated atomic indices in final list - only the
      !% minimum shift image of those found will be included. Default is false.
  
      real(dp), optional, intent(in)      :: max_r
      !% if present, only neighbors within this range will be included
  
      type(Connection), intent(in), optional, target :: alt_connect
      integer, intent(in), optional:: property(:)
  
      type(inoutput), optional :: debugfile
  
      integer, optional, intent(out) :: error
  
      !local
      logical                             :: do_nneighb_only, do_min_images_only
      integer                             :: i, j, n, m, jshift(3), ishift(3)
  
      integer :: n_i, keep_row(4), in_i, min_image
      integer, allocatable, dimension(:) :: repeats
      real(dp), allocatable, dimension(:) :: normsqshift
      type(Connection), pointer :: use_connect
  
      INIT_ERROR(error)
  
      if (present(debugfile)) call print("bfs_step", file=debugfile)
      if (present(alt_connect)) then
         use_connect => alt_connect
      else
        use_connect => this%connect
      endif
  
      if (present(debugfile)) call print("bfs_step cutoff " // this%cutoff, file=debugfile)
  
      if (.not.use_connect%initialised) then
           RAISE_ERROR('BFS_Step: Atomic structure has no connectivity data', error)
      endif
  
      do_nneighb_only = optional_default(.true., nneighb_only)
  
      do_min_images_only = optional_default(.false., min_images_only)
  
      if (present(debugfile)) call print('bfs_step: do_nneighb_only = ' // do_nneighb_only // ' do_min_images_only = '//do_min_images_only, file=debugfile)
      call print('bfs_step: do_nneighb_only = ' // do_nneighb_only // ' do_min_images_only = '//do_min_images_only, PRINT_NERD)
  
      if(input%intsize /= 4 .or. input%realsize /= 0) then
           RAISE_ERROR("bfs_step: input table must have intsize=4.", error)
      endif
  
      call allocate(output, 4, 0, 0, 0)
  
      ! Now go though the atomic indices
      do m = 1, input%N
         i = input%int(1,m)
         ishift = input%int(2:4,m)
  
         if (present(debugfile)) call print("bfs_step check atom " // m // " " // i // " with n_neigh " // n_neighbours(this, i, alt_connect=use_connect), file=debugfile)
         ! Loop over i's neighbours
         do n = 1, n_neighbours(this,i, alt_connect=use_connect)
  	  j = neighbour(this,i,n,shift=jshift,max_dist=max_r, alt_connect=use_connect)
  	  if (present(debugfile)) call print("bfs_step   check neighbour " // n // " " // j, file=debugfile)
  	  if (j == 0) cycle
  
            ! Look at next neighbour if j with correct shift is already in the cluster
            ! Must check input AND output tables
            if (find(input,(/j,ishift+jshift/)) > 0) cycle
  	  if (present(debugfile)) call print("bfs_step   not in input", file=debugfile)
            if (find(output,(/j,ishift+jshift/)) > 0) cycle
  	  if (present(debugfile)) call print("bfs_step   not in output", file=debugfile)
  
            if (do_nneighb_only .and. .not. is_nearest_neighbour(this, i, n, alt_connect=use_connect)) cycle
  	  if (present(debugfile)) call print("bfs_step   acceptably near neighbor", file=debugfile)
  
            if (present(property)) then
               if (property(j) == 0) cycle
            endif
  	  if (present(debugfile)) call print("bfs_step   property matches OK", file=debugfile)
  
            ! Everything checks out ok, so add j to the output table
            ! with correct shift
  	  if (present(debugfile)) call print("bfs_step   appending", file=debugfile)
            call append(output,(/j,ishift+jshift/)) 
         end do
      end do
  
      if (present(debugfile)) call print("bfs_step   raw output", file=debugfile)
      if (present(debugfile)) call print(output, file=debugfile)
  
      if (do_min_images_only) then
  
         ! If there are repeats of any atomic index, 
         ! we want to keep the one with smallest normsq(shift)
  
         ! must check input and output
  
         n = 1
         do while (n <= output%N)
  
            ! How many occurances in the output list?
            n_i = count(int_part(output,1) == output%int(1,n))
  
            ! Is there one in the input list as well?
            in_i = find_in_array(int_part(input,1),output%int(1,n))
  
            ! If there's only one, and we've not got one already then move on
            if (n_i == 1 .and. in_i == 0) then 
               n = n + 1
               cycle
            end if
  
            ! otherwise, things are more complicated...
            ! we want to keep the one with the smallest shift, bearing
            ! in mind that it could well be the one in the input list
  
            allocate(repeats(n_i), normsqshift(n_i))
  
            ! Get indices of repeats of this atomic index
            repeats = pack((/ (j, j=1,output%N) /), &
                 int_part(output,1) == output%int(1,n))
  
            if (in_i /= 0) then
               ! atom is in input list, remove all new occurances
               call delete_multiple(output, repeats)
            else
               ! Find row with minimum normsq(shift)
               normsqshift = normsq(real(output%int(2:4,repeats),dp),1)
               min_image = repeats(minloc(normsqshift,dim=1))
               keep_row = output%int(:,min_image)
  
               ! keep the minimum image
               call delete_multiple(output, &
                    pack(repeats, repeats /= min_image))
            end if
  
            ! don't increment n, since delete_multiple copies items from
            ! end of list over deleted items, so we need to retest
  
            deallocate(repeats, normsqshift)
  
         end do
  
      end if
  
      if (present(debugfile)) call print("bfs_step   final output", file=debugfile)
      if (present(debugfile)) call print(output, file=debugfile)
  
    end subroutine bfs_step
  
  
    !% Check if the list of indices and shifts 'list' contains
    !% any repeated atomic indices.
    function multiple_images(list)
      type(Table), intent(in) :: list
      logical :: multiple_images
  
  
      integer :: i
      multiple_images = .false.
  
      if (list%N == 0) return
  
      do i = 1, list%N
         if (count(int_part(list,1) == list%int(1,i)) /= 1) then
            multiple_images = .true.
            return
         end if
      end do
  
    end function multiple_images
  
    !% Given an input list with 4 integer columns for atomic indices
    !% and shifts, keep only the minimum images of each atom
    subroutine discard_non_min_images(list)
      type(Table), intent(inout) :: list
  
      integer :: n, n_i, j, min_image, keep_row(4)
      integer, allocatable, dimension(:) :: repeats
      real(dp), allocatable, dimension(:) :: normsqshift
  
      n = 1
      do while (n <= list%N)
  
         ! How many occurances in the list list?
         n_i = count(int_part(list,1) == list%int(1,n))
  
         ! If there's only one, and we've not got one already then move on
         if (n_i == 1) then 
            n = n + 1
            cycle
         end if
  
         ! otherwise, things are more complicated...
         ! we want to keep the one with the smallest shift
  
         allocate(repeats(n_i), normsqshift(n_i))
  
         ! Get indices of repeats of this atomic index
         repeats = pack((/ (j, j=1,list%N) /), &
              int_part(list,1) == list%int(1,n))
  
         ! Find row with minimum normsq(shift)
         normsqshift = normsq(real(list%int(2:4,repeats),dp),1)
         min_image = repeats(minloc(normsqshift,dim=1))
         keep_row = list%int(:,min_image)
  
         ! keep the minimum image
         call delete_multiple(list, &
              pack(repeats, repeats /= min_image))
  
         ! don't increment n, since delete_multiple copies items from
         ! end of list over deleted items, so we need to retest
  
         deallocate(repeats, normsqshift)
  
      end do
  
    end subroutine discard_non_min_images
  
  
    !% Add atoms to 'list' to make the selection region convex, i.e. if $i$ and
    !% $j$ are nearest neighbours, with $i$ in the list and not $j$ then $j$ will be added
    !% if more than half its nearest neighbours are in the list.
    subroutine make_convex(this, list)
      type(Atoms), intent(in) :: this
      type(Table), intent(inout) :: list
      type(Table)::tmplist
      do while(make_convex_step(this, list, tmplist) /= 0)
         call append(list, tmplist)
      end do
      call finalise(tmplist)
    end subroutine make_convex
  
    !% OMIT
    ! Given an input list, what needs to be added to make the selection region convex?
    function make_convex_step(this, input, output, error) result(newatoms)
      type(Atoms), intent(in) :: this
      type(Table), intent(in) :: input
      type(Table), intent(out) :: output
      integer, optional, intent(out) :: error
      integer::newatoms
  
      integer :: n, i, j, k, p, m, n_in, nn, ishift(3), jshift(3), kshift(3)
      real(dp) :: r_ij, r_kj
  
      INIT_ERROR(error)
  
      ! Check table size
      if(input%intsize /= 4 .or. input%realsize /= 0) then
           RAISE_ERROR("make_convex_step: input table must have intsize=4", error)
      endif
  
      if(input%intsize /= 4 .or. input%realsize /= 0) then
           RAISE_ERROR("bfs_step: input table must have intsize=4.", error)
      endif
  
      call allocate(output, 4, 0, 0, 0)
  
      do n=1,input%N
         i = input%int(1,n)
         ishift = input%int(2:4,n)
  
         !Loop over neighbours
         do m = 1, n_neighbours(this,i)
            j = neighbour(this,i,m, r_ij,shift=jshift)
  
            ! Look for nearest neighbours of i not in input list
            if (find(input,(/j,ishift+jshift/)) == 0 .and. is_nearest_neighbour(this, i, m)) then
  
               n_in = 0
               nn = 0
               ! Count number of nearest neighbours of j, and how many
               ! of them are in input list
               do p = 1, n_neighbours(this,j)
                  k = neighbour(this,j,p, r_kj, shift=kshift)
                  if (is_nearest_neighbour(this, j, p)) then
                     nn = nn + 1
                     if (find(input,(/k,ishift+jshift+kshift/)) /= 0) n_in = n_in + 1
                  end if
               end do
  
               !If more than half of  j's nearest neighbours are in then add it to output
               if (find(output, (/j,ishift+jshift/)) == 0 .and. real(n_in,dp)/real(nn,dp) > 0.5_dp) &
                    call append(output, (/j,ishift+jshift/))
  
            end if
  
         end do
      end do
  
      newatoms = output%N
  
    end function  make_convex_step
  
  
  ! Gotcha 1: Hollow sections 
  ! NB equivalent to reduce_n_cut_bonds when new number of bonds is 0 
  ! JRK: reduce_n_cut_bonds() is not equivalent to cluster_in_out_in()
  ! for silica, where adding IN-OUT-IN atoms doesn't change number of
  ! cut bonds OUT and IN refers to the list in cluster_info Look at the
  ! OUT nearest neighbours of IN atoms. If mode='all', then if ALL the
  ! nearest neighbours of the OUT atom are IN, then make the OUT atom
  ! IN. If mode='any', then add the OUT atom if ANY of it's nearest
  ! neighours (apart from the original IN atom)!  are IN. This stronger
  ! condition is required for cp2k.
  
  function cluster_in_out_in(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, mode) result(cluster_changed)
    type(Atoms), intent(in) :: this
    type(Table), intent(inout) :: cluster_info
    logical, intent(in) :: connectivity_just_from_connect
    type(Connection), intent(in) :: use_connect
    logical, intent(in) :: atom_mask(6)
    character(len=*), intent(in) :: mode
    logical :: cluster_changed
  
    integer :: n, i, ishift(3), m, j, jshift(3), p, k, kshift(3), n_nearest, n_in
    logical :: add_atom
  
    cluster_changed = .false.
    if (trim(mode) /= 'all' .and. trim(mode) /= 'any') &
         call system_abort('cluster_in_out_in: bad mode '//mode//' - should be "all" or "any"')
  
    n = 1
    ! Loop over cluster atoms (including ones that may get added in this loop)
    call print('cluster_in_out_in: Checking for hollow sections', PRINT_NERD)
    do while (n <= cluster_info%N)
      i = cluster_info%int(1,n)
      ishift = cluster_info%int(2:4,n)
      call print('cluster_in_out_in: i = '//i//' ['//ishift//'] Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
  
      !Loop over neighbours
      do m = 1, n_neighbours(this,i,alt_connect=use_connect)
      j = neighbour(this,i,m, shift=jshift,alt_connect=use_connect)
  
      if (find(cluster_info,(/j,ishift+jshift,this%Z(j),0/), atom_mask) == 0 .and. &
          (connectivity_just_from_connect .or. is_nearest_neighbour(this, i, m, alt_connect=use_connect)) ) then
        ! j is out and is nearest neighbour
  
        call print('cluster_in_out_in:   checking j = '//j//" ["//jshift//"]",PRINT_ANALYSIS)
  
        ! We have an OUT nearest neighbour, loop over its nearest neighbours to see if they are IN
  
        n_nearest = 0
        n_in = 0
        do p = 1, n_neighbours(this,j,alt_connect=use_connect)
           k = neighbour(this,j,p, shift=kshift,alt_connect=use_connect)
           if (k == i) cycle
           if (connectivity_just_from_connect .or. is_nearest_neighbour(this, j, p,alt_connect=use_connect)) then
              n_nearest = n_nearest + 1
              if (find(cluster_info,(/k,ishift+jshift+kshift,this%Z(k),0/), atom_mask) /= 0) then
                 n_in = n_in + 1
              end if
           end if
        end do
  
        if (trim(mode) == "all") then
           !If all of j's nearest neighbours are IN then add it
           add_atom = n_nearest == n_in
        else
           !If any of j's nearest neighbours (apart from i) are in then add it
           add_atom = n_in /= 0
        end if
  
        if (add_atom) then
          call append(cluster_info, (/j,ishift+jshift,this%Z(j),0/), (/this%pos(:,j), 1.0_dp/), (/ "inoutin   "/) )
          cluster_changed = .true.
          call print('cluster_in_out_in:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
        end if
  
      end if
  
      end do ! m
      n = n + 1
    end do ! while (n <= cluster_info%N)
  
    call print('cluster_in_out_in: Finished checking',PRINT_NERD)
    call print("cluster_in_out_in: cluster list:", PRINT_NERD)
    call print(cluster_info, PRINT_NERD)
  end function cluster_in_out_in
  
    !% Find cases where two IN atoms have a common
    !% OUT nearest neighbour, and see if termination would cause the hydrogen
    !% atoms to be too close together. If so, include the OUT nearest neighbour
    !% in the cluster
    !% returns true if cluster was changed
    function cluster_fix_termination_clash(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, termination_clash_factor) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      real(dp), intent(in) :: termination_clash_factor
      logical :: cluster_changed
  
      integer :: n, i, ishift(3), m, j, jshift(3), p, k, kshift(3)
      real(dp) :: dhat_ij(3), dhat_jk(3), r_ij, r_jk, H1(3), H2(3), diff_ik(3)
      ! real(dp) :: t_norm
  
      cluster_changed = .false.
  
      call print('doing cluster_fix_termination_clash', PRINT_NERD)
  
      !Loop over atoms in the cluster
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)     ! index of atom in the cluster
        ishift = cluster_info%int(2:4,n)
        call print('cluster_fix_termination_clash: i = '//i//'. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
        !Loop over atom i's neighbours
        do m = 1, n_neighbours(this,i,alt_connect=use_connect)
  	j = neighbour(this,i,m, shift=jshift, diff=dhat_ij, distance=r_ij,alt_connect=use_connect)
  	dhat_ij = dhat_ij/r_ij
  
  	!If j is IN the cluster, or not a nearest neighbour then try the next neighbour
  	if(find(cluster_info,(/j,ishift+jshift,this%Z(j),0/), atom_mask) /= 0) then
  	  call print('cluster_fix_termination_clash:   j = '//j//" ["//jshift//"] is in cluster",PRINT_ANALYSIS)
  	  cycle
  	end if
  	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	  call print('cluster_fix_termination_clash:   j = '//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	  cycle
  	end if
  
  	! So j is an OUT nearest neighbour of i.
  	call print('cluster_fix_termination_clash:   checking j = '//j//" ["//jshift//"]",PRINT_ANALYSIS)
  
  	!Determine the position of the would-be hydrogen along i--j
  	call print('cluster_fix_termination_clash:  Finding i--j hydrogen position',PRINT_ANALYSIS)
  
  	H1 = this%pos(:,i) + (this%lattice .mult. (ishift)) + &
  	     termination_bond_rescale(this%Z(i), this%Z(j)) * r_ij * dhat_ij
  
  	!Do a loop over j's nearest neighbours
  	call print('cluster_fix_termination_clash:  Looping over '//n_neighbours(this,j, alt_connect=use_connect)//' neighbours of j', PRINT_ANALYSIS)
  
  	do p = 1, n_neighbours(this,j, alt_connect=use_connect)
  	  k = neighbour(this,j,p, shift=kshift, diff=dhat_jk, distance=r_jk, alt_connect=use_connect)
  	  dhat_jk = dhat_jk/r_jk
  
  	  !If k is OUT of the cluster or k == i or it is not a nearest neighbour of j
  	  !then try the next neighbour
  
  	  if(find(cluster_info,(/k,ishift+jshift+kshift,this%Z(k),0/), atom_mask) == 0) then
  	    call print('cluster_fix_termination_clash:   k = '//k//" ["//kshift//"] not in cluster",PRINT_ANALYSIS)
  	    cycle
  	  end if
  	  if(k == i .and. all( jshift+kshift == 0 )) cycle
  	  if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,j, p, alt_connect=use_connect))) then
  	    call print('cluster_fix_termination_clash:   k = '//k//" ["//kshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	    cycle
  	  end if
  
  	  call print('cluster_fix_termination_clash: testing k = '//k//" ["//kshift//"]", PRINT_ANALYSIS)
  	  !Determine the position of the would-be hydrogen along k--j
  	  call print('cluster_fix_termination_clash:   Finding k--j hydrogen position',PRINT_ANALYSIS)
  
  	  diff_ik = r_ij * dhat_ij + r_jk * dhat_jk
  
  	  H2 = this%pos(:,i) + (this%lattice .mult. (ishift)) + diff_ik - &
  	    termination_bond_rescale(this%Z(k), this%Z(j)) * r_jk * dhat_jk
  	  call print('cluster_fix_termination_clash:   Checking i--k distance and hydrogen distance',PRINT_ANALYSIS)
  
  	  call print("cluster_fix_termination_clash: proposed hydrogen positions:", PRINT_ANALYSIS)
  	  call print(H1, PRINT_ANALYSIS)
  	  call print(H2, PRINT_ANALYSIS)
            call print('norm(H1-H2) '//norm(H1-H2)//' tol'//(bond_length(1,1)*this%nneightol*termination_clash_factor), PRINT_ANALYSIS)
  	  !NB workaround for pgf90 bug (as of 9.0-1)
  	  ! t_norm = norm(H1-H2); call print("cluster_fix_termination_clash: hydrogen distance would be "//t_norm, PRINT_ANALYSIS)
  	  !NB end of workaround for pgf90 bug (as of 9.0-1)
  	  ! If i and k are nearest neighbours, or the terminating hydrogens would be very close, then
  	  ! include j in the cluster. The H--H checking is conservative, hence the extra factor of 1.2
  	  if ((norm(diff_ik) < bond_length(this%Z(i),this%Z(k))*this%nneightol) .or. &
  	      (norm(H1-H2) < bond_length(1,1)*this%nneightol*termination_clash_factor) ) then
  
  	    call append(cluster_info,(/j,ishift+jshift,this%Z(j),0/),(/this%pos(:,j),1.0_dp/), (/ "clash     "/) )
  	    cluster_changed = .true.
  	    call print('cluster_fix_termination_clash:  Atom '//j//'with shift '//(ishift+jshift)//' added to cluster. Atoms = '//cluster_info%N, PRINT_NERD)
  	    ! j is now included in the cluster, so we can exit this do loop (over p)
  	    exit
  	  end if 
  
  	end do ! p
  
        end do ! m
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_fix_termination_clash: Finished checking',PRINT_NERD)
      call print("cluster_fix_termination_clash: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
    end function cluster_fix_termination_clash
  
    !% keep each whole residue (using atom_res_number(:) property) if any bit of it is already included
    function cluster_keep_whole_residues(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, present_keep_whole_residues) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical, intent(in) :: present_keep_whole_residues !% if true, keep_whole_residues was specified by user, not just a default value
      logical :: cluster_changed
      
      integer :: n, i, ishift(3), m, j, jshift(3)
      integer, pointer :: atom_res_number(:)
  
      call print('doing cluster_keep_whole_residues', PRINT_NERD)
  
      cluster_changed = .false.
      if (.not. assign_pointer(this, 'atom_res_number', atom_res_number)) then
        if (present_keep_whole_residues) then
  	call print("WARNING: cluster_keep_whole_residues got keep_whole_residues requested explicitly, but no proper atom_res_number property available", PRINT_ALWAYS)
        endif
        return
      endif
  
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
        if (atom_res_number(i) < 0) cycle
        call print('cluster_keep_whole_residues: i = '//i//' residue # = ' // atom_res_number(i) //'. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
        do m=1, n_neighbours(this, i, alt_connect=use_connect)
  	j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  
  	if (atom_res_number(i) /= atom_res_number(j)) then
  	  call print("cluster_keep_whole_residues:   j = "//j//" ["//jshift//"] has different res number " // atom_res_number(j), PRINT_ANALYSIS)
  	  cycle
  	endif
  	if(find(cluster_info,(/j,ishift+jshift,this%Z(j),0/), atom_mask) /= 0) then
  	  call print("cluster_keep_whole_residues:   j = "//j//" ["//jshift//"] is in cluster",PRINT_ANALYSIS)
  	  cycle
  	end if
  	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	  call print("cluster_keep_whole_residues:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	  cycle
  	end if
  
  	call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"res_num   "/) )
  	cluster_changed = .true.
  	call print('cluster_keep_whole_residues:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  
        end do ! m
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_keep_whole_residues: Finished checking',PRINT_NERD)
      call print("cluster_keep_whole_residues: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
  
    end function cluster_keep_whole_residues
  
    !% keep each whole proline (using atom_subgroup_number(:) and atom_res_number(:) property) if any bit of it is already included
    function cluster_keep_whole_prolines(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, present_keep_whole_prolines) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical, intent(in) :: present_keep_whole_prolines !% if true, keep_whole_prolines was specified by user, not just a default value
      logical :: cluster_changed
      
      integer :: n, i, ishift(3), m, j, jshift(3)
      integer, pointer :: atom_res_number(:), atom_subgroup_number(:)
      logical :: is_proline
  
      call print('doing cluster_keep_whole_prolines', PRINT_NERD)
  
      cluster_changed = .false.
      if (.not. assign_pointer(this, 'atom_res_number', atom_res_number) .or. .not. assign_pointer(this, 'atom_subgroup_number', atom_subgroup_number)) then
        if (present_keep_whole_prolines) then
  	call print("WARNING: cluster_keep_whole_prolines got keep_whole_prolines requested explicitly, but no proper atom_res_number or atom_subgroup_number properties available", PRINT_ALWAYS)
        endif
        return
      endif
  
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
        if (atom_res_number(i) < 0) cycle
        call print('cluster_keep_whole_prolines: i = '//i//' residue # = ' // atom_res_number(i) //' subgroup # = '// atom_subgroup_number(i) // &
                   '. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
        do m=1, n_neighbours(this, i, alt_connect=use_connect)
  	j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  	is_proline = .false.
  
  	if (atom_res_number(i) /= atom_res_number(j)) then
  	  call print("cluster_keep_whole_prolines:   j = "//j//" ["//jshift//"] has different res number " // atom_res_number(j), PRINT_ANALYSIS)
  	  cycle
  	endif
  	if (atom_subgroup_number(i) <= -10 .or. atom_subgroup_number(j) <= -10) then
  	  call print("cluster_keep_whole_prolines:   j = "//j//" ["//jshift//"] has proline subgroup number " // atom_subgroup_number(j), PRINT_ANALYSIS)
  	  is_proline = .true.
  	endif
  	if(find(cluster_info,(/j,ishift+jshift,this%Z(j),0/), atom_mask) /= 0) then
  	  call print("cluster_keep_whole_prolines:   j = "//j//" ["//jshift//"] is in cluster",PRINT_ANALYSIS)
  	  cycle
  	end if
  	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	  call print("cluster_keep_whole_prolines:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	  cycle
  	end if
  
  	if (is_proline) then
  	  call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"subgp_num "/) )
  	  cluster_changed = .true.
  	  call print('cluster_keep_whole_prolines:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  	endif
  
        end do ! m
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_keep_whole_prolines: Finished checking',PRINT_NERD)
      call print("cluster_keep_whole_prolines: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
  
    end function cluster_keep_whole_prolines
  
    !% keep each proline sidechain and directly connected bit of backbone (using atom_subgroup_number and atom_res_number(:) property) if any bit of it is already included
    function cluster_keep_whole_proline_sidechains(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, present_keep_whole_proline_sidechains) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical, intent(in) :: present_keep_whole_proline_sidechains !% if true, keep_whole_proline_sidechains was specified by user, not just a default value
      logical :: cluster_changed
      
      integer :: n, i, ishift(3), m, j, jshift(3)
      integer, pointer :: atom_res_number(:), atom_subgroup_number(:)
      logical :: is_proline
  
      call print('doing cluster_keep_whole_proline_sidechains', PRINT_NERD)
  
      cluster_changed = .false.
      if (.not. assign_pointer(this, 'atom_res_number', atom_res_number) .or. .not. assign_pointer(this, 'atom_subgroup_number', atom_subgroup_number)) then
        if (present_keep_whole_proline_sidechains) then
  	call print("WARNING: cluster_keep_whole_proline_sidechains got keep_whole_proline_sidechains requested explicitly, but no proper atom_res_number or atom_subgroup_number properties available", PRINT_ALWAYS)
        endif
        return
      endif
  
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
        if (atom_res_number(i) < 0) cycle
        call print('cluster_keep_whole_proline_sidechains: i = '//i//' residue # = ' // atom_res_number(i) //' subgroup # = '// atom_subgroup_number(i) // &
                   '. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
        do m=1, n_neighbours(this, i, alt_connect=use_connect)
  	j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  	is_proline = .false.
  
  	if (atom_res_number(i) /= atom_res_number(j)) then
  	  call print("cluster_keep_whole_proline_sidechains:   j = "//j//" ["//jshift//"] has different res number " // atom_res_number(j), PRINT_ANALYSIS)
  	  cycle
  	endif
  	if (atom_subgroup_number(i) == -1000 .or. atom_subgroup_number(j) == -1000) then
  	  call print("cluster_keep_proline_sidechains:   i or j = "//j//" ["//jshift//"] has proline sidechain subgroup number " // atom_subgroup_number(j), PRINT_ANALYSIS)
  	  is_proline = .true.
  	endif
  	if(find(cluster_info,(/j,ishift+jshift,this%Z(j),0/), atom_mask) /= 0) then
  	  call print("cluster_keep_proline_sidechains:   j = "//j//" ["//jshift//"] is in cluster",PRINT_ANALYSIS)
  	  cycle
  	end if
  	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	  call print("cluster_keep_whole_proline_sidechains:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	  cycle
  	end if
  
  	if (is_proline) then
  	  call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"subgp_num "/) )
  	  cluster_changed = .true.
  	  call print('cluster_keep_whole_proline_sidechains:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  	endif
  
        end do ! m
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_keep_whole_proline_sidechains: Finished checking',PRINT_NERD)
      call print("cluster_keep_whole_proline_sidechains: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
  
    end function cluster_keep_whole_proline_sidechains
  
    !% keep each whole subgroup (using atom_subgroup_number and atom_res_number(:) property) if any bit of it is already included
    function cluster_keep_whole_subgroups(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, present_keep_whole_subgroups) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical, intent(in) :: present_keep_whole_subgroups !% if true, keep_whole_subgroups was specified by user, not just a default value
      logical :: cluster_changed
      
      integer :: n, i, ishift(3), m, j, jshift(3)
      integer, pointer :: atom_res_number(:), atom_subgroup_number(:)
  
      call print('doing cluster_keep_whole_subgroups', PRINT_NERD)
  
      cluster_changed = .false.
      if (.not. assign_pointer(this, 'atom_res_number', atom_res_number) .or. .not. assign_pointer(this, 'atom_subgroup_number', atom_subgroup_number)) then
        if (present_keep_whole_subgroups) then
  	call print("WARNING: cluster_keep_whole_subgroups got keep_whole_subgroups requested explicitly, but no proper atom_res_number or atom_subgroup_number properties available", PRINT_ALWAYS)
        endif
        return
      endif
  
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
        if (atom_res_number(i) < 0) cycle
        if (atom_subgroup_number(i) == 0) cycle
        call print('cluster_keep_whole_subgroups: i = '//i//' residue # = ' // atom_res_number(i) //' subgroup # = '// atom_subgroup_number(i) // '. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
        do m=1, n_neighbours(this, i, alt_connect=use_connect)
  	j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  	call print("cluster_keep_whole_subgroups:    m = " // m //" j = " // j, PRINT_ANALYSIS)
  	if (atom_res_number(j) < 0) cycle
  	if (atom_subgroup_number(j) == 0) cycle
  
  	if (atom_res_number(i) /= atom_res_number(j)) then
  	  call print("cluster_keep_whole_subgroups:   j = "//j//" ["//jshift//"] has different res number " // atom_res_number(j), PRINT_ANALYSIS)
  	  cycle
  	endif
  	if (atom_subgroup_number(i) /= atom_subgroup_number(j)) then
  	  call print("cluster_keep_whole_subgroups:   j = "//j//" ["//jshift//"] has different subgroup number " // atom_subgroup_number(j), PRINT_ANALYSIS)
  	  cycle
  	endif
  	if(find(cluster_info,(/j,ishift+jshift,this%Z(j),0/), atom_mask) /= 0) then
  	  call print("cluster_keep_whole_subgroups:   j = "//j//" ["//jshift//"] is in cluster",PRINT_ANALYSIS)
  	  cycle
  	end if
  	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	  call print("cluster_keep_whole_subgroups:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	  cycle
  	end if
  
  	call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"subgp_num "/) )
  	cluster_changed = .true.
  	call print('cluster_keep_whole_subgroups:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  
        end do ! m
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_keep_whole_subgroups: Finished checking',PRINT_NERD)
      call print("cluster_keep_whole_subgroups: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
  
    end function cluster_keep_whole_subgroups
  
    !% keep whole silica tetrahedra -- that is, for each silicon atom, keep all it's oxygen nearest neighbours
    function cluster_keep_whole_silica_tetrahedra(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical :: cluster_changed
      
      integer :: n, i, ishift(3), m, j, jshift(3)
  
      call print('doing cluster_keep_whole_silica_tetrahedra', PRINT_NERD)
  
      cluster_changed = .false.
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
        if (this%z(i) /= 14) then
           ! Consider only silicon atoms, which form centres of tetrahedra
           cycle  
        end if
  
        call print('cluster_keep_whole_silica_tetrahedra: i = '//i//'. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
        do m=1, n_neighbours(this, i, alt_connect=use_connect)
  	j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  
  	if (this%z(j) /= 8) then
  	  call print("cluster_keep_whole_silica_tetrahedra:   j = "//j//" ["//jshift//"] is not oxygen ", PRINT_ANALYSIS)
  	  cycle
  	endif
  	if(find(cluster_info,(/j,ishift+jshift,this%Z(j),0/), atom_mask) /= 0) then
  	  call print("cluster_keep_whole_silica_tetrahedra:   j = "//j//" ["//jshift//"] is in cluster",PRINT_ANALYSIS)
  	  cycle
  	end if
  	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	  call print("cluster_keep_whole_silica_tetrahedra:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	  cycle
  	end if
  
  	call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"tetra     "/) )
  	cluster_changed = .true.
  	call print('cluster_keep_whole_silica_tetrahedra:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  
        end do ! m
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_keep_whole_silica_tetrahedra: Finished checking',PRINT_NERD)
      call print("cluster_keep_whole_silica_tetrahedra: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
  
    end function cluster_keep_whole_silica_tetrahedra
  
    !% keep whole octahedra for titania that is, for each Ti atom, keep all it's oxygen nearest neighbours but the O atoms with only 1 neighbour.
    function cluster_keep_whole_titania_octahedra(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical :: cluster_changed
      
      integer :: n, i, ishift(3), m, j, jshift(3)
  
      call print('doing cluster_keep_whole_titania_octahedra', PRINT_NERD)
  
      cluster_changed = .false.
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
        if (this%z(i) /= 22) then
           ! Consider only Ti atoms, which form centres of tetrahedra
           cycle  
        end if
  
        call print('cluster_keep_whole_titania_octahedra: i = '//i//'. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
        do m=1, n_neighbours(this, i, alt_connect=use_connect)
          j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  
          if (this%z(j) /= 8) then
            call print("cluster_keep_whole_titania_octahedra:   j = "//j//" ["//jshift//"] is not oxygen ", PRINT_ANALYSIS)
            cycle
          endif
          if(find(cluster_info,(/j,ishift+jshift,this%Z(j),0/), atom_mask) /= 0) then
            call print("cluster_keep_whole_titania_octahedra:   j = "//j//" ["//jshift//"] is in cluster",PRINT_ANALYSIS)
            cycle
          end if
          if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
            call print("cluster_keep_whole_titania_octahedra:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
            cycle
          end if
  
          !check if the oxygen has more than 1 neighbour in the cluster.
  !       neigh_O = 0
  !       do l=1,n_neighbours(this, j, alt_connect=use_connect)
  !           k = neighbour(this, j, l, distance=rdist, shift=kshift, alt_connect=use_connect)
  !           if(find(cluster_info,(/k,ishift+jshift+kshift,this%Z(k),0/), atom_mask) == 0) then
  !             call print(" j : "// j // " Z(j) :" // this%Z(j) // ' shift ' // jshift)
  !             call print(" k : "// k // " Z(k) :" // this%Z(k) // ' shift ' // kshift)
  !             call print(" i : "// i // " Z(i) :" // this%Z(i) // ' shift ' // ishift)
  !             call print('distance between k and j '// rdist)
  !             call print("cluster_keep_whole_titania_octahedra:   k = "//k//" ["//kshift//"] is not in cluster",PRINT_ANALYSIS)
  !             cycle
  !           end if
  !           if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,j, l, alt_connect=use_connect))) then
  !             call print(" j : "// j // " Z(j) :" // this%Z(j) // ' shift ' // jshift)
  !             call print(" k : "// k // " Z(k) :" // this%Z(k) // ' shift ' // kshift)
  !             call print(" i : "// i // " Z(i) :" // this%Z(i) // ' shift ' // ishift)
  !             call print('distance between k and j '// rdist)
  !             call print("cluster_keep_whole_titania_octahedra:   k = "//k//" ["//kshift//"] not nearest neighbour",PRINT_ANALYSIS)
  !             cycle
  !           end if
  !           neigh_O = neigh_O + 1 
  !       end do ! l
  
          call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"octa      "/) )
          cluster_changed = .true.
          call print('cluster_keep_whole_titania_octahedra  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
          
        end do ! m
  
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_keep_whole_titania_octahedra: Finished checking',PRINT_NERD)
      call print("cluster_keep_whole_titania_octahedra: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
  
    end function cluster_keep_whole_titania_octahedra
  
    !% adding each neighbour outside atom if doing so immediately reduces the number of cut bonds
    function cluster_reduce_n_cut_bonds(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical :: cluster_changed
  
      integer :: n, i, ishift(3), m, j, jshift(3), p, k, kshift(3)
      integer :: n_bonds_in, n_bonds_out
  
      call print('doing cluster_reduce_n_cut_bonds', PRINT_NERD)
  
      cluster_changed = .false.
  
      ! loop over each atom in the cluster already
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
        call print('cluster_reduce_n_cut_bonds: i = '//i//'. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
  
        ! loop over every neighbour, looking for outside neighbours
        do m=1, n_neighbours(this, i, alt_connect=use_connect)
  	j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  	if (find(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), atom_mask ) /= 0) cycle
  	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	  call print("cluster_reduce_n_cut_bonds:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	  cycle
  	end if
  	! if we're here, j must be an outside neighbour of i
  	call print('cluster_reduce_n_cut_bonds: j = '//j//'. Looping over '//n_neighbours(this,j,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
  
  	n_bonds_in = 0
  	n_bonds_out = 0
  	do p=1, n_neighbours(this, j, alt_connect=use_connect)
  	  k = neighbour(this, j, p, shift=kshift, alt_connect=use_connect)
  	  if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,j, p, alt_connect=use_connect))) then
  	    call print("cluster_reduce_n_cut_bonds:   k = "//k//" ["//kshift//"] not nearest neighbour of j = " // j,PRINT_ANALYSIS)
  	    cycle
  	  end if
  	  ! count how many bonds point in vs. out
  	  if (find(cluster_info, (/ k, ishift+jshift+kshift, this%Z(k), 0 /), atom_mask ) /= 0) then
  	    n_bonds_in = n_bonds_in + 1
  	  else
  	    n_bonds_out = n_bonds_out + 1
  	  endif
  	end do ! p
  
  	if (n_bonds_out < n_bonds_in) then ! adding this one would reduce number of cut bonds
  	  call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"n_cut_bond" /) )
  	  cluster_changed = .true.
  	  call print('cluster_reduce_n_cut_bonds:  Added atom ' //j//' ['//(ishift+jshift)//'] n_bonds_in ' // n_bonds_in // &
  		     ' out ' // n_bonds_out // '  to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  	endif
        end do ! m
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_reduce_n_cut_bonds: Finished checking',PRINT_NERD)
      call print("cluster_reduce_n_cut_bonds: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
    end function cluster_reduce_n_cut_bonds
  
    !% Go through the cluster atoms and find cut bonds. If the bond is to
    !% hydrogen, then include the hydrogen in the cluster list (since AMBER
    !% doesn't cap a bond to hydrogen with a link atom)
    function cluster_protect_X_H_bonds(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical :: cluster_changed
  
      integer :: n, i, ishift(3), m, j, jshift(3)
  
      call print('doing cluster_protect_X_H_bonds', PRINT_NERD)
      cluster_changed = .false.
  
      ! loop over each atom in the cluster already
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
        call print('cluster_protect_X_H_bonds: i = '//i//'. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
  
        ! loop over every neighbour, looking for outside neighbours
        do m=1, n_neighbours(this, i, alt_connect=use_connect)
  	j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  
  	! if j is in, or not a nearest neighbour, go on to next
  	if (find(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), atom_mask ) /= 0) cycle
  	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	  call print("cluster_protect_X_H_bonds:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	  cycle
  	end if
  	! if we're here, j must be an outside neighbour of i
  
  	if (this%Z(i) == 1 .or. this%Z(j) == 1) then
  	  call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"X_H_bond  " /) )
  	  cluster_changed = .true.
  	  call print('cluster_protect_X_H_bonds:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  	endif
        end do ! m
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_protect_X_H_bonds: Finished checking',PRINT_NERD)
      call print("cluster_protect_X_H_bonds: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
    end function cluster_protect_X_H_bonds
  
    !% Go through the cluster atoms and find cut bonds. 
    !% Include the bonded-to atom in the cluster list.
    !% cheap alternative to keep_residues_whole
    function cluster_keep_whole_molecules(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical :: cluster_changed
  
      integer :: n, i, ishift(3), m, j, jshift(3)
  
      call print('doing cluster_keep_whole_molecules', PRINT_NERD)
      cluster_changed = .false.
  
      ! loop over each atom in the cluster already
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
        call print('cluster_keep_whole_molecules: i = '//i//'. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
  
        ! loop over every neighbour, looking for outside neighbours
        do m=1, n_neighbours(this, i, alt_connect=use_connect)
  	j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  
  	! if j is in, or not a nearest neighbour, go on to next
  	if (find(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), atom_mask ) /= 0) cycle
  	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	  call print("cluster_keep_whole_molecules:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	  cycle
  	end if
  	! if we're here, j must be an outside neighbour of i
  
  	call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"molec     " /) )
  	cluster_changed = .true.
  	call print('cluster_keep_whole_molecules:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
        end do ! m
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_keep_whole_molecules: Finished checking',PRINT_NERD)
      call print("cluster_keep_whole_molecules: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
    end function cluster_keep_whole_molecules
  
    !if by accident a single atom is included, that is part of a larger entity (not ion)
    !then delete it from the list as it will cause error in SCF (e.g. C=O oxygen)
    !call allocate(to_delete,4,0,0,0,1)
    !do n = 1, cluster%N
    !        i = cluster%int(1,n)
    !        is_alone = .true.
    !        do j = 1, n_neighbours(at,i)
    !                k = neighbour(at,i,j)
    !                if(find(cluster_info,(/k,shift/)) /= 0) then
    !                   if(is_nearest_neighbour(at,i,j)) GOTO 120
    !                else
    !                   if(is_nearest_neighbour(at,i,j)) is_alone = .false.
    !                end if
    !        end do
    !        if(.not. is_alone) call append(to_delete,(/i,shift/))
    !   120  cycle
    !end do
  
    !do i = 1, to_delete%N
    ! call delete(cluster_info,to_delete%int(:,i))
    !end do
  
  !OUTDATED   function cluster_biochem_in_out_in(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask) result(cluster_changed)
  !OUTDATED     type(Atoms), intent(in) :: this
  !OUTDATED     type(Table), intent(inout) :: cluster_info
  !OUTDATED     logical, intent(in) :: connectivity_just_from_connect
  !OUTDATED     type(Connection), intent(in) :: use_connect
  !OUTDATED     logical, intent(in) :: atom_mask(6)
  !OUTDATED     logical :: cluster_changed
  !OUTDATED 
  !OUTDATED     integer :: n, i, ishift(3), m, j, jshift(3), p, k, kshift(3)
  !OUTDATED 
  !OUTDATED     call print('doing cluster_biochem_in_out_in', PRINT_NERD)
  !OUTDATED     cluster_changed = .false.
  !OUTDATED 
  !OUTDATED     ! loop over each atom in the cluster already
  !OUTDATED     n = 1
  !OUTDATED     do while (n <= cluster_info%N)
  !OUTDATED       i = cluster_info%int(1,n)
  !OUTDATED       ishift = cluster_info%int(2:4,n)
  !OUTDATED       call print('cluster_reduce_n_cut_bonds: i = '//i//'. Looping over '//n_neighbours(this,i,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
  !OUTDATED 
  !OUTDATED       ! loop over every neighbour, looking for outside neighbours
  !OUTDATED       do m=1, n_neighbours(this, i, alt_connect=use_connect)
  !OUTDATED 	j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  !OUTDATED 
  !OUTDATED 	! if j is in, or not a neareste neighbour, go on to next
  !OUTDATED 	if (find(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), atom_mask ) /= 0) cycle
  !OUTDATED 	if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  !OUTDATED 	  call print("cluster_reduce_n_cut_bonds:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  !OUTDATED 	  cycle
  !OUTDATED 	end if
  !OUTDATED 	! if we're here, j must be an outside neighbour of i
  !OUTDATED 	call print('cluster_reduce_n_cut_bonds: j = '//j//'. Looping over '//n_neighbours(this,j,alt_connect=use_connect)//' neighbours...',PRINT_ANALYSIS)
  !OUTDATED 
  !OUTDATED 	do p=1, n_neighbours(this, j, alt_connect=use_connect)
  !OUTDATED 	  k = neighbour(this, j, p, shift=kshift, alt_connect=use_connect)
  !OUTDATED 	  if (k == i .and. all(jshift + kshift == 0)) cycle
  !OUTDATED 	  if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,j, p, alt_connect=use_connect))) then
  !OUTDATED 	    call print("cluster_reduce_n_cut_bonds:   k = "//k//" ["//kshift//"] not nearest neighbour of j = " // j,PRINT_ANALYSIS)
  !OUTDATED 	    cycle
  !OUTDATED 	  end if
  !OUTDATED 	  ! if k is in, then j has 2 in neighbours, and so we add it
  !OUTDATED 	  if (find(cluster_info, (/ k, ishift+jshift+kshift, this%Z(k), 0 /), atom_mask ) /= 0) then
  !OUTDATED 	    call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"bio_IOI   " /) )
  !OUTDATED 	    cluster_changed = .true.
  !OUTDATED 	    call print('cluster_biochem_in_out_in:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  !OUTDATED 	    exit
  !OUTDATED 	  end if 
  !OUTDATED 	end do ! p
  !OUTDATED       end do ! m
  !OUTDATED 
  !OUTDATED       n = n + 1
  !OUTDATED     end do ! while (n <= cluster_info%N)
  !OUTDATED 
  !OUTDATED     call print('cluster_biochem_in_out_in: Finished checking',PRINT_NERD)
  !OUTDATED     call print("cluster_biochem_in_out_in: cluster list:", PRINT_NERD)
  !OUTDATED     call print(cluster_info, PRINT_NERD)
  !OUTDATED   end function cluster_biochem_in_out_in
  
    !% if an in non-H atom is in a residue and not at right coordination, add all of its neighbours
    function cluster_protect_double_bonds(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, present_protect_double_bonds) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical, intent(in) :: present_protect_double_bonds !% if true, protect_double_bonds was specified by user, not just a default value
      logical :: cluster_changed
  
      integer :: n, i, ishift(3), m, j, jshift(3)
      integer ::  n_nearest_neighbours
      integer, pointer :: atom_res_number(:)
  
      call print('doing cluster_protect_double_bonds', PRINT_NERD)
      cluster_changed = .false.
  
      if (.not. assign_pointer(this, 'atom_res_number', atom_res_number)) then
        if (present_protect_double_bonds) then
  	call print("WARNING: cluster_protect_double_bonds got protect_double_bonds requested explicitly, but no proper atom_res_number property available", PRINT_ALWAYS)
        endif
        return
      endif
  
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
  
        if (this%Z(i) /= 0 .and. atom_res_number(i) >= 0) then
  
  	! count nearest neighbours of i
  	n_nearest_neighbours = 0
  	do m=1, n_neighbours(this, i, alt_connect=use_connect)
  	  j = neighbour(this, i, m, alt_connect=use_connect)
  
  	  ! if j is not nearest neighbour, go on to next one
  	  if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	    call print("cluster_protect_double_bonds:   j = "//j//" not nearest neighbour",PRINT_ANALYSIS)
  	    cycle
  	  end if
  	  ! if we're here, j must be an outside neighbour of i
  	  n_nearest_neighbours = n_nearest_neighbours + 1
  	end do ! m
  
  	if (ElementValence(this%Z(i)) /= -1 .and. (ElementValence(this%Z(i)) /= n_nearest_neighbours)) then ! i has a valence, and it doesn't match nn#
  	  do m=1, n_neighbours(this, i)
  	    j = neighbour(this, i, m, shift=jshift, alt_connect=use_connect)
  
  	    ! if j is in, or isn't nearest neighbour, go on to next
  	    if (find(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), atom_mask ) /= 0) cycle
  	    if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, m, alt_connect=use_connect))) then
  	      call print("cluster_protect_double_bonds:   j = "//j//" ["//jshift//"] not nearest neighbour",PRINT_ANALYSIS)
  	      cycle
  	    end if
  
  	    ! if we're here, j must be an outside neighbour of i
  	    call append(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), (/ this%pos(:,j), 1.0_dp /), (/"dbl_bond  " /) )
  	    cluster_changed = .true.
  	    call print('cluster_protect_double_bonds:  Added atom ' //j//' ['//(ishift+jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  	  end do ! m
  	endif ! valence defined
        end if !  atom_res_number(i) >= 0
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_protect_double_bonds: Finished checking',PRINT_NERD)
      call print("cluster_protect_double_bonds: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
    end function cluster_protect_double_bonds
  
    !% protect bond between R-C ... -N-R
    !%                        |      |
    !%                        O      R (usually H or C (proline))
    function cluster_protect_peptide_bonds(this, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, present_protect_peptide_bonds) result(cluster_changed)
      type(Atoms), intent(in) :: this !% atoms structure 
      type(Table), intent(inout) :: cluster_info !% table of cluster info, modified if necessary on output
      logical, intent(in) :: connectivity_just_from_connect !% if true, we're doing hysterestic connect and should rely on the connection object completely
      type(Connection), intent(in) :: use_connect !% connection object to use for connectivity info
      logical, intent(in) :: atom_mask(6) !% which fields in int part of table to compare when checking for identical atoms
      logical, intent(in) :: present_protect_peptide_bonds !% if true, protect_peptide_bonds was specified by user, not just a default value
      logical :: cluster_changed
  
      integer :: n, i, ishift(3), ji, j, ki, k, jshift(3)
      integer, pointer :: atom_res_number(:)
      integer :: found_other_j, found_other_jshift(3)
      integer :: this_neighbour_Z, other_neighbour_Z, other_Z
      logical :: found_this_neighbour, found_other_neighbour
  
      call print('doing cluster_protect_peptide_bonds', PRINT_NERD)
      cluster_changed = .false.
  
      if (.not. assign_pointer(this, 'atom_res_number', atom_res_number)) then
        if (present_protect_peptide_bonds) then
  	call print("WARNING: cluster_protect_peptide_bonds got protect_peptide_bonds requested explicitly, but no proper atom_res_number property available", PRINT_ALWAYS)
        endif
        return
      endif
  
      n = 0
      do while (n < cluster_info%N)
        n = n + 1
        i = cluster_info%int(1,n)
        ishift = cluster_info%int(2:4,n)
  
        if (this%Z(i) /= 0 .and. atom_res_number(i) >= 0) then
  
  	call print("cluster_protect_peptide_bonds: i = "//i//" Z(i) " // this%Z(i), PRINT_ANALYSIS)
  
  	! peptide bond is C-N
  	!                 |
  	!                 O
  	if (this%Z(i) == 6) then
  	  other_Z = 7
  	  this_neighbour_Z = 8
  	  other_neighbour_Z = 0
  	else if (this%Z(i) == 7) then
  	  other_Z = 6
  	  this_neighbour_Z = 0
  	  other_neighbour_Z = 8
  	else
  	  call print("cluster_protect_peptide_bonds:   i = "//i//" not C or N",PRINT_ANALYSIS)
  	  cycle
  	endif
  
  	found_other_j = 0
  	found_this_neighbour = (this_neighbour_Z == 0)
  	found_other_neighbour = (other_neighbour_Z == 0)
  
  	call print("cluster_protect_peptide_bonds: other_Z " // other_Z // " this_neighbour_Z " // this_neighbour_Z // " other_neighbour_z " // other_neighbour_Z, PRINT_ANALYSIS)
  	call print("cluster_protect_peptide_bonds: found_this_neighbour " // found_this_neighbour // " found_other_neighbour " // found_other_neighbour, PRINT_ANALYSIS)
  
  	! must have 3 neighbours
  	if (n_neighbours(this, i, alt_connect=use_connect) /= 3) then
  	  call print("cluster_protect_peptide_bonds:   i = "//i//" doesn't have 3 neighbours",PRINT_ANALYSIS)
  	  cycle
  	endif
  	do ji=1, n_neighbours(this, i, alt_connect=use_connect) ! look for correct neighbour
  	  j=neighbour(this, i, ji, shift=jshift, alt_connect=use_connect)
  	  call print("cluster_protect_peptide_bonds:   j = " //j//" Z(j) " // this%Z(j), PRINT_ANALYSIS)
  	  ! only nearest neighbours count, usually
  	  if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,i, ji, alt_connect=use_connect))) then
  	    call print("cluster_protect_peptide_bonds:   neighbour j = "//j//" not nearest neighbour",PRINT_ANALYSIS)
  	    cycle
  	  end if
  	  if (this%Z(j) == this_neighbour_Z .and. n_neighbours(this, j) == 1) then ! neighbour has right Z and coordination
  	    found_this_neighbour = .true.
  	    call print("cluster_protect_peptide_bonds:   neighbour j = "//j//" is neighbor (O) we're looking for",PRINT_ANALYSIS)
  	  endif
  	  if (this%Z(j) == other_Z) then ! found other
  	    ! if j is already in, this isn't a peptide bond in need of protection from being cut
  	    if (find(cluster_info, (/ j, ishift+jshift, this%Z(j), 0 /), atom_mask) /= 0) then
  	      call print("cluster_protect_peptide_bonds:   other_Z neighbour j = "//j//" already in",PRINT_ANALYSIS)
  	      cycle
  	    endif
  	    ! save info on other
  	    call print("cluster_protect_peptide_bonds: found_other_j " // j // " Z(found_other_j) " // this%Z(j), PRINT_ANALYSIS)
  	    found_other_j = j
  	    found_other_jshift = jshift
  
  	    ! if (.not. found_this_neighbour .or. found_other_j == 0) then ! some part of motif missing
  	      ! call print("cluster_protect_peptide_bonds:   failed to find this_neighbour_Z or other_j",PRINT_ANALYSIS)
  	      ! cycle
  	    ! endif
  
  	    ! other must have 3 neighbours
  	    if (n_neighbours(this, found_other_j, alt_connect=use_connect) /= 3) then
  	      call print("cluster_protect_peptide_bonds:   found_other_j = "//found_other_j//" doesn't have 3 neighbours", PRINT_ANALYSIS)
  	      cycle
  	    endif
  	    ! perhaps we should require that other be in a different residue, but not doing that now
  	    ! if (atoms_res_number(i) == atom_res_number(found_other_j)) then
  	      ! call print("cluster_protect_peptide_bonds:   i and found_other_j are in the same residue", PRINT_ANALYSIS)
  	      ! cycle
  	    ! endif
  
  	    ! check neighbours of other
  	    do ki=1, n_neighbours(this, found_other_j, alt_connect=use_connect)
  	      k=neighbour(this, found_other_j, ki, alt_connect=use_connect)
  	      ! only nearest neighbours count
  	      if(.not. (connectivity_just_from_connect .or. is_nearest_neighbour(this,found_other_j, ki, alt_connect=use_connect))) then
  		call print("cluster_protect_peptide_bonds:   found_other_j's neighbour k = "//k//" not nearest neighbour",PRINT_ANALYSIS)
  		cycle
  	      end if
  
  	      if (this%Z(k) == other_neighbour_Z .and. n_neighbours(this, k) == 1) then ! neighbour has right Z and coordination
  		found_other_neighbour = .true.
  		call print("cluster_protect_peptide_bonds:   neighbour k = "//k//" is neighbor (O) we're looking for",PRINT_ANALYSIS)
  	      endif
  	    end do ! ki
  
  	    ! some part of motif is missing
  	    if (.not. found_other_neighbour) then
  	      call print("cluster_protect_peptide_bonds:   couldn't find other_neighbour",PRINT_ANALYSIS)
  	      cycle
  	    endif
  
  	    call append(cluster_info, (/ found_other_j, ishift+found_other_jshift, this%Z(found_other_j), 0 /), &
  			(/ this%pos(:,found_other_j), 1.0_dp /), (/"pep_bond  " /) )
  	    cluster_changed = .true.
  	    call print('cluster_protect_peptide_bonds:  Added atom ' //found_other_j//' ['//(ishift+found_other_jshift)//'] to cluster. Atoms = ' // cluster_info%N, PRINT_NERD)
  	  endif ! Z(j) == other_Z
  	end do ! ji
  
        end if ! Z(i) /= 0 .and. atom_res_number(i) /= 0
  
      end do ! while (n <= cluster_info%N)
  
      call print('cluster_protect_peptide_bonds: Finished checking',PRINT_NERD)
      call print("cluster_protect_peptide_bonds: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
    end function cluster_protect_peptide_bonds
  
    subroutine create_cluster_simple(at, args_str, cluster, mark_name, error)
       type(Atoms), intent(inout) :: at
       character(len=*), intent(in) :: args_str
       type(Atoms), intent(out) :: cluster
       character(len=*), intent(in), optional :: mark_name
       integer, intent(out), optional :: ERROR
  
       type(Table) :: cluster_info
  
       INIT_ERROR(error)
  
       call calc_connect(at)
       cluster_info = create_cluster_info_from_mark(at, args_str, mark_name=mark_name, error=error)
       PASS_ERROR(error)
       call carve_cluster(at, args_str, cluster_info, cluster=cluster, mark_name=mark_name, error=error)
       PASS_ERROR(error)
       call finalise(cluster_info)
       
    end subroutine create_cluster_simple
  
    !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    !
    !% Create cluster from atoms and cluster information table (which
    !% which should be generated by a call to :func:`create_cluster_info_from_hybrid_mark`.
    !%
    !% The output cluster contains all properties of the initial atoms object, and
    !% some additional columns, which are:
    !%
    !%  ``index``
    !%      index of the cluster atoms into the initial atoms object.
    !%
    !%  ``termindex`` 
    !%     nonzero for termination atoms, and is an index into
    !%     the cluster atoms specifiying which atom is being terminated,
    !%     it is used in collecting the forces.
    !%
    !%  ``rescale`` 
    !%      a real number which for nontermination atoms is 1.0,
    !%      for termination atoms records the scaling applied to
    !%      termination bond lengths
    !%  
    !%  ``shift``
    !%      the shift of each atom
    !
    !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    subroutine carve_cluster(at, args_str, cluster_info, cluster, mark_name, error)
      type(Atoms), intent(in), target :: at
      character(len=*), intent(in) :: args_str
      type(Table), intent(in) :: cluster_info
      type(Atoms), intent(out) :: cluster
      character(len=*), optional, intent(in) :: mark_name
      integer, optional, intent(out) :: error
  
      type(Dictionary) :: params
      logical :: do_rescale_r
      real(dp) :: r_scale, cluster_vacuum, cluster_fix_lattice(3), cluster_round_lattice
      logical :: terminate, randomise_buffer, print_clusters, do_calc_connect, do_same_lattice, do_fix_lattice
      logical :: hysteretic_connect
      integer :: i, j, k, m
      real(dp) :: maxlen(3), sep(3), lat_maxlen(3), lat_sep(3)
      logical :: do_periodic(3)
      integer, pointer :: hybrid_mark(:), cluster_index(:), cluster_hybrid_mark(:), cluster_shift(:,:)
      type(CInoutput)                    :: clusterfile
      character(len=255)                :: clusterfilename
      character(len=255)                :: my_mark_name
      character(STRING_LENGTH)          :: run_suffix
      type(Table) :: outer_layer
      logical :: in_outer_layer, do_map_into_cell
  
  #ifdef _MPI
      integer::mpi_size, mpi_rank, PRINT_ALWAYS
      integer :: mpi_force_size
      real(dp), allocatable, dimension(:)  :: mpi_force
  
      call get_mpi_size_rank(MPI_COMM_WORLD, mpi_size, mpi_rank)
  #endif /* _MPI */
  
      INIT_ERROR(error)
  
      call print('carve_cluster got args_str "'//trim(args_str)//'"', PRINT_VERBOSE)
  
      call initialise(params)
      call param_register(params, 'terminate', 'T', terminate,&
        help_string="do hydrogen termination")
      call param_register(params, 'do_rescale_r', 'F', do_rescale_r,&
        help_string="do rescale cluster lattice positions and lattice (for matching lattice constants of different models")
      call param_register(params, 'r_scale', '1.0', r_scale,&
        help_string="cluster position/lattice rescaling factor")
      call param_register(params, 'randomise_buffer', 'T', randomise_buffer,&
        help_string="randomly perturb positions of buffer atoms")
      call param_register(params, 'print_clusters', 'F', print_clusters,&
        help_string="Print out clusters into clusters.xyz file (clusters.MPI_RANK.xyz if in MPI)")
      call param_register(params, 'cluster_calc_connect', 'F', do_calc_connect,&
        help_string="run calc_connect before doing stuff (passivation, etc)")
      call param_register(params, 'cluster_same_lattice', 'F', do_same_lattice,&
        help_string="Don't modify cluster cell vectors from incoming cell vectors")
      call param_register(params, 'cluster_fix_lattice', '0. 0. 0.', cluster_fix_lattice, &
        has_value_target=do_fix_lattice, &
        help_string="fix the cluster lattice to be cubic with given 3 lengths. useful with DFT codes where one doesn't pay for vacuum.")
      call param_register(params, 'cluster_round_lattice', '3.0', cluster_round_lattice, &
        help_string="distance to which to round cluster lattice (default is 3.0 A)")
      call param_register(params, 'cluster_periodic_x', 'F', do_periodic(1),&
        help_string="do not add vaccum along x, so as to make cluster's cell vectors periodic")
      call param_register(params, 'cluster_periodic_y', 'F', do_periodic(2),&
        help_string="do not add vaccum along x, so as to make cluster's cell vectors periodic")
      call param_register(params, 'cluster_periodic_z', 'F', do_periodic(3),&
        help_string="do not add vaccum along x, so as to make cluster's cell vectors periodic")
      call param_register(params, 'cluster_vacuum', '10.0', cluster_vacuum,&
        help_string="Amount of vacuum to add around cluster (minimum - lattice vectors are rounded up to integers)")
      call param_register(params, 'hysteretic_connect', 'F', hysteretic_connect,&
        help_string="Use hysteretic connection object for identifying bonds")
      call param_register(params, 'map_into_cell', 'T', do_map_into_cell,&
        help_string="Call map_into_cell (lattice coordinates=[-0.5,0.5)) on cluster")
      call param_register(params, 'run_suffix', '', run_suffix, &
           help_string="string to append to names of extra properties added to cluster. NOTE: does not override mark_name argument")
      if (.not. param_read_line(params, args_str, ignore_unknown=.true.,task='carve_cluster arg_str') ) then
        RAISE_ERROR("carve_cluster failed to parse args_str='"//trim(args_str)//"'", error)
      endif
      call finalise(params)
  
      ! first pick up an atoms structure with the right number of atoms and copy the properties
      !Now turn the cluster_temp table into an atoms structure
      call print('carve_cluster: Copying atomic data to output object',PRINT_NERD)
      call print('List of atoms in cluster:', PRINT_NERD)
      call print(int_part(cluster_info,1), PRINT_NERD)
  
      call select(cluster, at, list=int_part(cluster_info,1), error=error)
      PASS_ERROR(error)
      ! then reset the positions species and Z (latter two needed because termination atoms have Z=1)
      ! unfold the positions to real positions using the stored shifts, at is neede because
      ! next we might make the unit cell much smaller
      do i=1,cluster_info%N
         cluster%pos(:,i) = cluster_info%real(1:3, i)+(at%lattice .mult. cluster_info%int(2:4, i))
         cluster%Z(i) = cluster_info%int(5,i)
         cluster%species(:,i) = ' '
         cluster%species(1:3,i) = s2a(ElementName(cluster_info%int(5,i)))
      end do
      ! add properties to cluster
      call add_property(cluster, 'index'//trim(run_suffix), int_part(cluster_info,1))
      call add_property(cluster, 'shift'//trim(run_suffix), 0, n_cols=3, ptr2=cluster_shift)
      cluster_shift(:,1:cluster%N) = cluster_info%int(2:4,1:cluster_info%N)
      call add_property(cluster, 'termindex'//trim(run_suffix), int_part(cluster_info,6))
      call add_property(cluster, 'rescale'//trim(run_suffix), real_part(cluster_info,4))
      call add_property(cluster, 'cluster_ident'//trim(run_suffix), cluster_info%str(1,1:cluster_info%N))
  
      ! Find smallest bounding box for cluster
      ! Find boxes aligned with xyz (maxlen) and with a1 a2 a3 (lat_maxlen)
      maxlen = 0.0_dp
      lat_maxlen = 0.0_dp
      do i=1,cluster%N
         do j=1,cluster%N
  	  sep = cluster%pos(:,i)-cluster%pos(:,j)
  	  lat_sep = cluster%g .mult. sep
  	  do k=1,3
  	     if (abs(sep(k)) > maxlen(k)) maxlen(k) = abs(sep(k))
  	     if (abs(lat_sep(k)) > lat_maxlen(k)) lat_maxlen(k) = abs(lat_sep(k))
  	  end do
         end do
      end do
  
      ! renormalize lat_maxlen to real dist units
      do k=1,3
        lat_maxlen(k) = lat_maxlen(k) * norm(cluster%lattice(:,k))
      end do
  
      ! Round up maxlen to be divisible by cluster_round_lattice (default 3 A), so that it does not fluctuate too much
      forall (k=1:3) maxlen(k) = cluster_round_lattice*ceiling(maxlen(k)/cluster_round_lattice)
      forall (k=1:3) lat_maxlen(k) = cluster_round_lattice*ceiling(lat_maxlen(k)/cluster_round_lattice)
  
      ! vacuum pad cluster (if didn't set do_same_lattice or cluster_fix_lattice)
      ! if not periodic at all, just do vacuum padding
      ! if cluster_fix_lattice was given, just set the diagonal lattice elements to those values
      ! if periodic along some dir, keep supercell vector directions, and set
      !    extent in each direction to lesser of cluster extent + vacuum or original extent
      if (do_same_lattice) then
        cluster%lattice = at%lattice
      else if (do_fix_lattice) then
        cluster%lattice = 0.0_dp
        do k=1,3
           cluster%lattice(k,k) = cluster_fix_lattice(k)
        end do
      else
        if (any(do_periodic)) then
  	do k=1,3
  	  if (do_periodic(k)) then
  	    if (lat_maxlen(k)+cluster_vacuum >= norm(at%lattice(:,k))) then
  	      cluster%lattice(:,k) = at%lattice(:,k)
  	    else
  	      cluster%lattice(:,k) = (lat_maxlen(k)+cluster_vacuum)*at%lattice(:,k)/norm(at%lattice(:,k))
  	    endif
  	  else
  	    cluster%lattice(:,k) = (lat_maxlen(k)+cluster_vacuum)*at%lattice(:,k)/norm(at%lattice(:,k))
  	  endif
  	end do
        else
  	cluster%lattice = 0.0_dp
  	do k=1,3
  	  cluster%lattice(k,k) = maxlen(k) + cluster_vacuum
  	end do
        endif
      endif
  
      call set_lattice(cluster, cluster%lattice, scale_positions=.false.)
  
      ! Remap positions so any image atoms end up inside the cell
      if (do_map_into_cell) then
         call map_into_cell(cluster)
      end if
  
      call print ('carve_cluster: carved cluster with '//cluster%N//' atoms', PRINT_VERBOSE)
  
      write (line, '(a,f10.3,f10.3,f10.3)') &
           'carve_cluster: Cluster dimensions are ', cluster%lattice(1,1), &
           cluster%lattice(2,2), cluster%lattice(3,3)
      call print(line, PRINT_VERBOSE)
  
      my_mark_name = optional_default('hybrid_mark', mark_name)
      ! reassign pointers
      if (.not. assign_pointer(at, trim(my_mark_name), hybrid_mark)) then
           RAISE_ERROR('cannot reassign "'//trim(my_mark_name)//'" property', error)
      endif
  
      ! rescale cluster positions and lattice 
      if (do_rescale_r) then
         call print('carve_cluster: rescaling cluster positions (and lattice if not fixed) by factor '//r_scale, PRINT_VERBOSE)
         cluster%pos = r_scale * cluster%pos
         if (.not. do_fix_lattice) then
            call set_lattice(cluster, r_scale * cluster%lattice, scale_positions=.false.)
         end if
      end if
  
      if (randomise_buffer .and. .not. any(hybrid_mark == HYBRID_BUFFER_OUTER_LAYER_MARK)) &
           do_calc_connect = .true.
  
      if (do_calc_connect) then
         call print('carve_cluster: doing calc_connect', PRINT_VERBOSE)
         ! Does QM force model need connectivity information?
         call set_cutoff(cluster, at%cutoff)
         call calc_connect(cluster)
      end if
  
      if (randomise_buffer) then
         ! If outer buffer layer is not marked do so now. This will be
         ! the case if we're using hysteretic_buffer selection option.
         ! In this case we consider any atoms connected to terminating
         ! hydrogens to be in the outer layer - obviously this breaks down 
         ! if we're not terminating so we abort if that's the case.
         if (.not. assign_pointer(cluster, 'index', cluster_index)) then
              RAISE_ERROR('carve_cluster: cluster is missing index property', error)
         endif
  
         if (.not. any(hybrid_mark == HYBRID_BUFFER_OUTER_LAYER_MARK)) then
            if (.not. terminate) then
  	    RAISE_ERROR('cannot determine which buffer atoms to randomise if terminate=F and hysteretic_buffer=T', error)
  	  endif
  
            if (.not. assign_pointer(cluster, trim(my_mark_name), cluster_hybrid_mark)) then
                 RAISE_ERROR(trim(my_mark_name)//' property not found in cluster', error)
  	  endif
  
            do i=1,cluster%N
               if (hybrid_mark(cluster_info%int(1,i)) /= HYBRID_BUFFER_MARK) cycle
               in_outer_layer = .false.
               do m = 1, cluster_info%N
                  if (cluster_info%int(6,m).eq. i) then
                     in_outer_layer = .true.
                     exit
                  end if
               enddo
               if (in_outer_layer) then
                  hybrid_mark(cluster_info%int(1,i)) = HYBRID_BUFFER_OUTER_LAYER_MARK
                  cluster_hybrid_mark(i) = HYBRID_BUFFER_OUTER_LAYER_MARK
               end if
            end do
         end if
  
        if (any(hybrid_mark == HYBRID_BUFFER_OUTER_LAYER_MARK)) then
  	 ! Slightly randomise the positions of the outermost layer
  	 ! of the buffer region in order to avoid systematic errors
  	 ! in QM forces.
  
  	 call initialise(outer_layer, 1,0,0,0)
  	 call append(outer_layer, find(hybrid_mark == HYBRID_BUFFER_OUTER_LAYER_MARK))
  
  	 do i=1,outer_layer%N
  	    ! Shift atom by randomly distributed vector of magnitude 0.05 A
  	    j = find_in_array(cluster_index, outer_layer%int(1,i))
  	    cluster%pos(:,j) = cluster%pos(:,j) + 0.05_dp*random_unit_vector()
  	 end do
  
  	 call finalise(outer_layer)
        end if
  
      end if
  
      if (value(mainlog%verbosity_stack) >= PRINT_VERBOSE .or. print_clusters) then
  #ifdef _MPI
         write (clusterfilename, '(a,i3.3,a)') 'clusters.',mpi_rank,'.xyz'
  #else
         clusterfilename = 'clusters.xyz'
  #endif /* _MPI */
         call initialise(clusterfile, clusterfilename, append=.true., action=OUTPUT)
         call write(clusterfile, cluster)
         call finalise(clusterfile)
      end if
  
   end subroutine carve_cluster
  
    !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    !
    !% Create a cluster using the mark_name (optional arg, default 'hybrid_mark') 
    !% property and options in 'args_str'.
    !% All atoms that are marked with anything other than 'HYBRID_NO_MARK' will
    !% be included in the cluster; this includes active, transition and buffer
    !% atoms.   Atoms added (to satisfy heuristics) will be marked in new property
    !% modified_//trim(mark_name) as 'HYBRID_BUFFER_MARK'.
    !% Returns a Table object (cluster_info) which contains info on atoms whose
    !% indices are given in atomlist, possibly with some extras for consistency,
    !% and optionally terminated with Hydrogens, that can be used by carve_cluster().
    !
    !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  
    function create_cluster_info_from_mark(at, args_str, cut_bonds, mark_name, error) result(cluster_info)
      type(Atoms), intent(inout), target :: at
      character(len=*), intent(in) :: args_str
      type(Table), optional, intent(out)   :: cut_bonds !% Return a list of the bonds cut when making
      !% the cluster.  See create_cluster() documentation.
      character(len=*), intent(in), optional :: mark_name
      integer, optional, intent(out) :: error
      type(Table) :: cluster_info
  
      type(Dictionary) :: params
      logical :: terminate, periodic_x, periodic_y, periodic_z, &
           even_electrons, do_periodic(3), cluster_hopping_nneighb_only, cluster_heuristics_nneighb_only, &
           cluster_allow_modification, hysteretic_connect, same_lattice, &
           fix_termination_clash, keep_whole_residues, keep_whole_subgroups, keep_whole_prolines, keep_whole_proline_sidechains, &
  	 keep_whole_silica_tetrahedra, keep_whole_titania_octahedra, terminate_octahedra,reduce_n_cut_bonds, in_out_in, &
           protect_X_H_bonds, protect_double_bonds, protect_peptide_bonds, keep_whole_molecules, has_termination_rescale, &
  	 combined_protein_heuristics, force_no_fix_termination_clash
      character(STRING_LENGTH) :: in_out_in_mode
      logical :: keep_whole_residues_has_value, keep_whole_subgroups_has_value, keep_whole_prolines_has_value, keep_whole_proline_sidechains_has_value, &
                 protect_double_bonds_has_value, protect_peptide_bonds_has_value, keep_whole_molecules_has_value
      real(dp) :: r, r_min, centre(3), termination_rescale, termination_clash_factor
      type(Table) :: cluster_list, currentlist, nextlist, activelist, bufferlist
      integer :: i, j, jj, first_active, old_n, n_cluster
      integer, pointer :: hybrid_mark(:), modified_hybrid_mark(:)
      integer :: prev_cluster_info_n
      integer, allocatable, dimension(:) :: uniqed, tmp_index
  
      type(Table)                              :: n_term, in_n_term, sorted_n_term
      integer                                  :: m, n, p, p_in, in_n
  
      real(dp)                                 :: H1(3)
      real(dp)                                 :: dhat_ij(3), d(3)
      real(dp)                                 :: r_ij, rescale
      real(dp)                                 :: orig_length
      real(dp), dimension(3)                   :: bond1, bond2, new_bond 
      integer                                  :: ishift(3), jshift(3), oldN, most_hydrogens
      logical                                  :: atom_mask(6)
      integer, allocatable, dimension(:) :: idx
  
      type(Connection), pointer :: use_connect
      logical :: connectivity_just_from_connect
      logical :: cluster_changed, cluster_hopping
      character(len=STRING_LENGTH) :: my_mark_name
  
      INIT_ERROR(error)
  
      my_mark_name = optional_default('hybrid_mark', mark_name)
  
      call print('create_cluster_info_from_mark got args_str "'//trim(args_str)//'"', PRINT_VERBOSE)
  
      call initialise(params)
      call param_register(params, 'terminate', 'T', terminate, &
        help_string="Do hydrogen termination for dangling bonds")
      call param_register(params, 'cluster_periodic_x', 'F', periodic_x,&
        help_string="do not add vaccum along x, so as to make cluster's cell vectors periodic")
      call param_register(params, 'cluster_periodic_y', 'F', periodic_y,&
        help_string="do not add vaccum along y, so as to make cluster's cell vectors periodic")
      call param_register(params, 'cluster_periodic_z', 'F', periodic_z,&
        help_string="do not add vaccum along z, so as to make cluster's cell vectors periodic")
      call param_register(params, 'even_electrons', 'F', even_electrons,&
        help_string="Remove passivation hydrogen to keep total number of electrons even")
      call param_register(params, 'cluster_hopping_nneighb_only', 'T', cluster_hopping_nneighb_only,&
        help_string="Apply is_nearest_neighbor test to bonds when doing bond hops to find cluster.  If false, uses connectivity object as is")
      call param_register(params, 'cluster_heuristics_nneighb_only', 'T', cluster_heuristics_nneighb_only,&
        help_string="Apply is_nearest_neighbor test to bonds when deciding what's bonded for heuristics.  If false, uses connectivity object as is")
      call param_register(params, 'cluster_allow_modification', 'T', cluster_allow_modification,&
        help_string="Allow cluster to be modified using various heuristics")
      call param_register(params, 'hysteretic_connect', 'F', hysteretic_connect,&
        help_string="Use hysteretic connect object to decide what's bonded.  Usually goes with cluster_*_nneighb_only=F.")
      call param_register(params, 'cluster_same_lattice', 'F', same_lattice,&
        help_string="Don't modify cluster cell vectors from incoming cell vectors")
      call param_register(params, 'fix_termination_clash','T', fix_termination_clash,&
        help_string="Apply termination clash (terimnation H too close together) heureistic")
      call param_register(params, 'force_no_fix_termination_clash','F', force_no_fix_termination_clash,&
        help_string="Overrides values of (terminate .or. termination_clash) to explicitly disable termination clash heuristic")
      call param_register(params, 'combined_protein_heuristics','F', combined_protein_heuristics, &
        help_string="Apply heuristics for proteins: overrides keep_whole_molecules=F, keep_whole_residues=F, keep_whole_subgroups=T, keep_whole_prolines=T, keep_whole_proline_sidechains=F, protect_double_bonds=F, protect_peptide_bonds=T")
      call param_register(params, 'keep_whole_residues','T', keep_whole_residues, has_value_target=keep_whole_residues_has_value,&
        help_string="Apply heuristic keeping identified residues whole")
      call param_register(params, 'keep_whole_prolines','F', keep_whole_prolines, has_value_target=keep_whole_prolines_has_value,&
        help_string="Apply heuristic keeping identified prolines whole")
      call param_register(params, 'keep_whole_proline_sidechains','F', keep_whole_proline_sidechains, has_value_target=keep_whole_proline_sidechains_has_value,&
        help_string="Apply heuristic keeping identified prolines sidechain+directly connected backbone part whole")
      call param_register(params, 'keep_whole_subgroups','F', keep_whole_subgroups, has_value_target=keep_whole_subgroups_has_value,&
        help_string="Apply heuristic keeping identified subgroups within a residue whole")
      call param_register(params, 'keep_whole_silica_tetrahedra','F', keep_whole_silica_tetrahedra,&
        help_string="Apply heureistic keeping SiO2 tetrahedra whole")
      call param_register(params, 'keep_whole_titania_octahedra','F', keep_whole_titania_octahedra,&
        help_string="Apply heureistic keeping TiO2 tetragonal structure whole")
      call param_register(params, 'terminate_octahedra','F', terminate_octahedra,&
        help_string="Apply heureistic terminating TiO2 octahedra")
      call param_register(params, 'reduce_n_cut_bonds','T', reduce_n_cut_bonds,&
        help_string="Apply heuristic that adds atoms if doing so reduces number of broken bonds")
      call param_register(params, 'in_out_in', 'F', in_out_in, &
        help_string="Apply heuristic than adds atoms to avoid IN-OUT-IN configurations. Usually equivalent to reduce_n_cut_bonds.")
      call param_register(params, 'in_out_in_mode', 'all', in_out_in_mode, &
        help_string="Mode to use for in_out_in heuristic - can be either 'all' or 'any'. Default is 'all'.")
      call param_register(params, 'protect_X_H_bonds','T', protect_X_H_bonds,&
        help_string="Apply heuristic protecting X-H bonds - no point H passivating bonds with an H")
      call param_register(params, 'protect_double_bonds','T', protect_double_bonds, has_value_target=protect_double_bonds_has_value,&
        help_string="Apply heuristic protecting double bonds from being cut (based one some heuristic idea of what's a double bond")
      call param_register(params, 'protect_peptide_bonds','F', protect_peptide_bonds, has_value_target=protect_peptide_bonds_has_value,&
        help_string="Apply heuristic protecting peptide bonds from being cut (RO-C...N-RR)")
      call param_register(params, 'keep_whole_molecules','F', keep_whole_molecules, has_value_target=keep_whole_molecules_has_value, &
        help_string="Apply heuristic keeping identified molecules (i.e. contiguous bonded groups of atoms) whole")
      call param_register(params, 'termination_rescale', '0.0', termination_rescale, has_value_target=has_termination_rescale,&
        help_string="rescale factor for X-H passivation bonds")
      call param_register(params, 'termination_clash_factor', '1.2', termination_clash_factor, &
        help_string="Factor to multiply H-H bond-length by when deciding if two Hs are too close. Default 1.2")
      call param_register(params, 'cluster_hopping', 'T', cluster_hopping,&
        help_string="Identify cluster region by hopping from a seed atom.  Needed so that shifts will be correct.")
  
      if (.not. param_read_line(params, args_str, ignore_unknown=.true.,task='create_cluster_info_from_mark args_str') ) then
           RAISE_ERROR("create_cluster_info_from_mark failed to parse args_str='"//trim(args_str)//"'", error)
      endif
      call finalise(params)
  
      if (combined_protein_heuristics) then
        if (keep_whole_molecules_has_value) call print("WARNING: create_cluster_info_from_mark got combined_protein_heuristics, overriding keep_whole_molecules")
        if (keep_whole_residues_has_value) call print("WARNING: create_cluster_info_from_mark got combined_protein_heuristics, overriding keep_whole_residues")
        if (keep_whole_subgroups_has_value) call print("WARNING: create_cluster_info_from_mark got combined_protein_heuristics, overriding keep_whole_subgroups")
        if (keep_whole_prolines_has_value) call print("WARNING: create_cluster_info_from_mark got combined_protein_heuristics, overriding keep_whole_prolines")
        if (keep_whole_proline_sidechains_has_value) call print("WARNING: create_cluster_info_from_mark got combined_protein_heuristics, overriding keep_whole_proline_sidechains")
        if (protect_double_bonds_has_value) call print("WARNING: create_cluster_info_from_mark got combined_protein_heuristics, overriding protect_double_bonds_has_value")
        if (protect_peptide_bonds_has_value) call print("WARNING: create_cluster_info_from_mark got combined_protein_heuristics, overriding protect_peptide_bonds_has_value")
        keep_whole_molecules = .false.
        keep_whole_molecules_has_value = .true.
        keep_whole_residues = .false.
        keep_whole_residues_has_value = .true.
        keep_whole_subgroups = .true.
        keep_whole_subgroups_has_value = .true.
        keep_whole_prolines = .true.
        keep_whole_prolines_has_value = .true.
        keep_whole_proline_sidechains = .true.
        keep_whole_proline_sidechains_has_value = .true.
        protect_double_bonds = .false.
        protect_peptide_bonds = .true.
        protect_peptide_bonds_has_value = .true.
      endif
  
      do_periodic = (/periodic_x,periodic_y,periodic_z/)
  
      if (.not. has_property(at, trim(my_mark_name))) then
         RAISE_ERROR('create_cluster_info_from_mark: atoms structure has no "'//trim(my_mark_name)//'" property', error)
      endif
  
      call assign_property_pointer(at, trim(my_mark_name), hybrid_mark, error=error)
      PASS_ERROR_WITH_INFO('create_cluster_info_from_mark passed atoms structure with no hybrid_mark property', error)
  
      ! Calculate centre of cluster
      call allocate(cluster_list, 1,0,0,0)
      call append(cluster_list, find(hybrid_mark /= HYBRID_NO_MARK))
      ! unreliable when cluster extends across periodic boundaries
      ! centre = 0.0_dp
      ! do i=1,cluster_list%N
      ! centre = centre + at%pos(:,cluster_list%int(1,i)) + (at%lattice .mult. cluster_list%int(2:4,i))
      ! end do
      ! centre = centre / cluster_list%N
      centre = at%pos(:,cluster_list%int(1,1))
  !!$    call print('centre = '//centre)
  
      ! Find ACTIVE atom closest to centre of cluster, using min image convention  
      r_min = huge(1.0_dp)
      do i=1,cluster_list%N
         if (hybrid_mark(cluster_list%int(1,i)) /= HYBRID_ACTIVE_MARK) cycle
         r = distance_min_image(at, centre, cluster_list%int(1,i))
  !!$       call print(cluster_list%int(1,i)//' '//r)
         if (r < r_min) then
            first_active = cluster_list%int(1,i)
            r_min = r
         end if
      end do
  
      if (cluster_hopping) then
         call system_timer('cluster_hopping')
  
         n_cluster = cluster_list%N
         call wipe(cluster_list)
         call allocate(cluster_list, 4,0,0,0)
  
         ! Add first marked atom to cluster_list. shifts will be relative to this atom
         call print('Growing cluster starting from atom '//first_active//', n_cluster='//n_cluster, PRINT_VERBOSE)
         call append(cluster_list, (/first_active,0,0,0/))
         call append(currentlist, cluster_list)
  
         ! Add other active atoms using bond hopping from the central cluster atom
         ! to find the other cluster atoms and hence to determine the correct 
         ! periodic shifts. 
         !
         ! This will fail if marked atoms do not form a single connected cluster
         old_n = cluster_list%N
         do 
            call wipe(nextlist)
            if (hysteretic_connect) then
               call BFS_step(at, currentlist, nextlist, nneighb_only=cluster_hopping_nneighb_only, &
                    min_images_only = any(do_periodic) .or. same_lattice , alt_connect=at%hysteretic_connect)
            else
               call BFS_step(at, currentlist, nextlist, nneighb_only=cluster_hopping_nneighb_only, &
                    min_images_only = any(do_periodic) .or. same_lattice)
            endif
            do j=1,nextlist%N
               jj = nextlist%int(1,j)
               !          shift = nextlist%int(2:4,j)
               if (hybrid_mark(jj) /= HYBRID_NO_MARK .and. find(cluster_list, nextlist%int(:,j)) == 0) &
                    call append(cluster_list, nextlist%int(:,j))
            end do
            call append(currentlist, nextlist)
  
            ! check exit condition
            allocate(tmp_index(cluster_list%N))
            tmp_index = int_part(cluster_list,1)
            call sort_array(tmp_index)
            call uniq(tmp_index, uniqed)
            call print('cluster hopping: got '//cluster_list%N//' atoms, of which '//size(uniqed)//' are unique.', PRINT_VERBOSE)
            if (size(uniqed) == n_cluster) exit !got them all
            deallocate(uniqed, tmp_index)
  
            ! check that cluster is still growing
            if (cluster_list%N == old_n) then
               call write(at, 'create_cluster_abort.xyz')
               call print(cluster_list)
               RAISE_ERROR('create_cluster_info_from_mark: cluster stopped growing before all marked atoms found - check for split QM region', error)
            end if
            old_n = cluster_list%N
         end do
         deallocate(tmp_index, uniqed)
         call finalise(nextlist)
         call finalise(currentlist)
         call system_timer('cluster_hopping')
      else
         call system_timer('cluster_diff_min_image')
         call allocate(cluster_list, 4,0,0,0)
  
         call append(cluster_list, (/first_active,0,0,0/))
         do i=1, at%n
            if (hybrid_mark(i) == HYBRID_NO_MARK .or. i == first_active) cycle
            ! shifts relative to first_active
            d = distance_min_image(at, first_active, i, shift=ishift)
            call append(cluster_list, (/i,ishift/))
         end do
         call system_timer('cluster_diff_min_image')
      end if
      ! partition cluster_list so that active atoms come first
      call wipe(activelist)
      call wipe(bufferlist)
      do i=1,cluster_list%N
         if (hybrid_mark(cluster_list%int(1,i)) == HYBRID_ACTIVE_MARK) then
            call append(activelist, cluster_list%int(:,i))
         else
            call append(bufferlist, cluster_list%int(:,i))
         end if
      end do
  
      call wipe(cluster_list)
      call append(cluster_list, activelist)
      call append(cluster_list, bufferlist)
      call finalise(activelist)
      call finalise(bufferlist)
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
      ! 
      ! Validate arguments
      !
  
      ! check for consistency in optional arguments
  
      if (.not. (count(do_periodic) == 0 .or. count(do_periodic) == 1 .or. count(do_periodic) == 3)) then
           RAISE_ERROR('count(periodic) must be zero, one or three.', error)
      endif
  
      if (same_lattice) do_periodic = .true.
  
      if (any(do_periodic) .and. multiple_images(cluster_list)) then
           RAISE_ERROR("create_cluster: can't make a periodic cluster since cluster_list contains repeats", error)
      endif
  
      ! check for empty list
  
      if(cluster_list%N == 0) then
         call print('create_cluster: empty cluster_list', PRINT_NORMAL)
         return
      end if
  
      call print('create_cluster: Entering create_cluster', PRINT_NERD)
  
      if(.not.(cluster_list%intsize == 1 .or. cluster_list%intsize == 4) .or. cluster_list%realsize /= 0) then
           RAISE_ERROR("create_cluster: cluster_list table must have intsize=1 or 4 and realsize=0.", error)
      endif
  
  
      ! Cluster_info is extensible storage for the cluster
      ! It stores atomic indices and shifts (4 ints)
      ! atomic number (1 int)
      ! termination index (1 int): for termination atoms, which atom is being terminated?
      ! and atomic positions (3 reals)
      ! It's length will be at least cluster_list%N
  
      call print('create_cluster: Creating temporary cluster table', PRINT_NERD)
      call allocate(cluster_info,6,4,1,0,cluster_list%N)
  
      ! First, put all the marked atoms into cluster_info, storing their positions and shifts
      call print('create_cluster: Adding specified atoms to the cluster', PRINT_NERD)
      do i = 1, cluster_list%N
         if(cluster_list%intsize == 4) then
            ! we have shifts
            ishift = cluster_list%int(2:4,i)
         else
            ! no incoming shifts
            ishift = (/0,0,0/)
         end if
         call append(cluster_info, (/cluster_list%int(1,i),ishift,at%Z(cluster_list%int(1,i)),0/),&
  	    (/at%pos(:,cluster_list%int(1,i)),1.0_dp/), (/ hybrid_mark_name(hybrid_mark(cluster_list%int(1,i))) /) )
      end do
  
      call print("create_cluster: cluster list:", PRINT_NERD)
      call print(cluster_info, PRINT_NERD)
  
      call system_timer('cluster_consistency')
      ! Next, check for various gotchas
  
      ! at mask is used to match with atoms already in the cluster_info table
      ! if we are periodic in a direction, we don't care about the shifts in that direction when matching
      atom_mask = (/.true.,.not.do_periodic, .true., .true./)
  
      connectivity_just_from_connect = .not. cluster_heuristics_nneighb_only
      if (hysteretic_connect) then
         use_connect => at%hysteretic_connect
         ! will also pass true for connectivity_just_from_connect to cluster_...() routines, so that
         ! is_nearest_neighbour() won't be used
         connectivity_just_from_connect = .true.
      else
         use_connect => at%connect
      endif
  
      if (cluster_allow_modification) then
         call add_property(at, 'modified_'//trim(my_mark_name), 0, ptr=modified_hybrid_mark)
         cluster_changed = .true.
         modified_hybrid_mark = hybrid_mark
         prev_cluster_info_n = cluster_info%N
  
         do while (cluster_changed) 
            cluster_changed = .false.
            call print("fixing up cluster according to heuristics keep_whole_residues " // keep_whole_residues // &
                 ' keep_whole_prolines ' // keep_whole_prolines // &
                 ' keep_whole_proline_sidechains ' // keep_whole_proline_sidechains // &
                 ' keep_whole_subgroups ' // keep_whole_subgroups // &
                 ' keep_whole_silica_tetrahedra ' // keep_whole_silica_tetrahedra // &
                 ' keep_whole_titania_octahedra ' // keep_whole_titania_octahedra // &
                 ' terminate_octahedra ' // terminate_octahedra // &
                 ' reduce_n_cut_bonds ' // reduce_n_cut_bonds // &
                 ' in_out_in ' // in_out_in // &
                 ' in_out_in_mode ' // trim(in_out_in_mode) // &
                 ' protect_X_H_bonds ' // protect_X_H_bonds // &
                 ' protect_double_bonds ' // protect_double_bonds // &
                 ' protect_peptide_bonds ' // protect_peptide_bonds // &
                 ' terminate .or. fix_termination_clash ' // (terminate .or. fix_termination_clash) // &
                 ' force_no_fix_termination_clash '// force_no_fix_termination_clash, verbosity=PRINT_NERD)
            if (keep_whole_proline_sidechains) then
               cluster_changed = cluster_changed .or. cluster_keep_whole_proline_sidechains(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, keep_whole_proline_sidechains_has_value)
            endif
            if (keep_whole_prolines) then
               cluster_changed = cluster_changed .or. cluster_keep_whole_prolines(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, keep_whole_prolines_has_value)
            endif
            if (keep_whole_subgroups) then
               cluster_changed = cluster_changed .or. cluster_keep_whole_subgroups(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, keep_whole_subgroups_has_value)
            endif
            if (keep_whole_residues) then
               cluster_changed = cluster_changed .or. cluster_keep_whole_residues(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, keep_whole_residues_has_value)
            endif
            if (keep_whole_silica_tetrahedra) then
               cluster_changed = cluster_changed .or. cluster_keep_whole_silica_tetrahedra(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask)
            end if
            if (reduce_n_cut_bonds) then
               cluster_changed = cluster_changed .or. cluster_reduce_n_cut_bonds(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask)
            endif
            if (keep_whole_titania_octahedra) then
               cluster_changed = cluster_changed .or. cluster_keep_whole_titania_octahedra(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask)
            end if
            if (in_out_in) then
               cluster_changed = cluster_changed .or. cluster_in_out_in(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, in_out_in_mode)
            endif
            if (protect_X_H_bonds) then
               cluster_changed = cluster_changed .or. cluster_protect_X_H_bonds(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask)
            endif
            if (protect_double_bonds) then
               cluster_changed = cluster_changed .or. cluster_protect_double_bonds(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, protect_double_bonds_has_value)
            endif
            if (protect_peptide_bonds) then
               cluster_changed = cluster_changed .or. cluster_protect_peptide_bonds(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, protect_peptide_bonds_has_value)
            endif
            if ((terminate .or. fix_termination_clash) .and. .not. force_no_fix_termination_clash) then
               cluster_changed = cluster_changed .or. cluster_fix_termination_clash(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask, termination_clash_factor)
            endif
            if (keep_whole_molecules) then
               cluster_changed = cluster_changed .or. cluster_keep_whole_molecules(at, cluster_info, connectivity_just_from_connect, use_connect, atom_mask)
            endif
         end do ! while cluster_changed
  
         if (prev_cluster_info_n < cluster_info%N) modified_hybrid_mark(cluster_info%int(1,prev_cluster_info_n+1:cluster_info%N)) = HYBRID_BUFFER_MARK
  
         call print('create_cluster: Finished fixing cluster for various heuristic pathologies',PRINT_NERD)
         call print("create_cluster: cluster list:", PRINT_NERD)
         call print(cluster_info, PRINT_NERD)
      end if ! allow_cluster_mod
  
      call system_timer('cluster_consistency')
  
      call system_timer('cluster_termination')
  
      !So now cluster_info contains all the atoms that are going to be in the cluster.
      !If terminate is set, we need to add terminating hydrogens along nearest neighbour bonds
      if (terminate.or.terminate_octahedra) then
         call print('create_cluster: Terminating cluster with hydrogens',PRINT_NERD)
  
         call allocate(n_term, 5, 0, 0, 0)
         call allocate(in_n_term, 5, 0, 0, 0)
         oldN = cluster_info%N
  
         !Loop over atoms in the cluster
         do n = 1, oldN
  
            i = cluster_info%int(1,n)
            ishift = cluster_info%int(2:4,n)
            if(terminate_octahedra .and. at%z(i).ne.8 ) cycle 
            !Loop over atom i's neighbours
            do m = 1, n_neighbours(at,i, alt_connect=use_connect)
  
               j = neighbour(at,i,m, r_ij, diff=dhat_ij, shift=jshift, alt_connect=use_connect)
               dhat_ij = dhat_ij / r_ij
  
               if (find(cluster_info,(/j,ishift+jshift,at%Z(j),0/), atom_mask) == 0 .and. &
                    (hysteretic_connect .or. is_nearest_neighbour(at, i, m, alt_connect=use_connect))) then
  
                  ! If j is an OUT atom, and it is close enough, put a terminating hydrogen
                  ! at the scaled distance between i and j
  
                  if (.not. has_termination_rescale) then
                     rescale = termination_bond_rescale(at%Z(i), at%Z(j))
                  else
                     rescale = termination_rescale
                  end if
                  H1 = at%pos(:,i) + rescale * r_ij * dhat_ij
  
                  ! Label term atom with indices into original atoms structure.
                  ! j is atom it's generated from and n is index into cluster table of atom it's attached to
                  call append(cluster_info,(/j,ishift,1,n/),(/H1, rescale/), (/ "term      " /)) 
                   
                  ! Keep track of how many termination atoms each cluster atom has
                  p = find_in_array(int_part(n_term,(/1,2,3,4/)),(/n,ishift/))
                  if (p == 0) then
                     call append(n_term, (/n,ishift,1/))
                  else
                     n_term%int(5,p) = n_term%int(5,p) + 1
                  end if
  
                  ! optionally keep a record of the bonds that we have cut
                  if (present(cut_bonds)) &
                       call append(cut_bonds, (/i,j,ishift,jshift/))
  
                  if(current_verbosity() .ge. PRINT_NERD) then
                     write(line,'(a,i0,a,i0,a)')'create_cluster: Replacing bond ',i,'--',j,' with hydrogen'
                     call print(line, PRINT_NERD)
                  end if
               end if
  
               ! Keep track of how many other neighbours has
               if (terminate_octahedra.and.find(cluster_info,(/j,ishift+jshift,at%Z(j),0/), atom_mask) /= 0 .and. &
                    (hysteretic_connect .or. is_nearest_neighbour(at, i, m, alt_connect=use_connect))) then
                  p_in = find_in_array(int_part(in_n_term,(/1,2,3,4/)),(/n,ishift/))
                  if (p_in == 0) then
                     call append(in_n_term, (/n,ishift,1/))
                  else
                     in_n_term%int(5,p_in) = in_n_term%int(5,p_in) + 1
                  end if
               endif
            end do
  
         end do
  
         !If terminate octahedra is set we have to adjust how many H atoms are added to the oxygens.
         !Since with Ti, the oxygens are threefold coordinated, we have to remove some of the added O
         !to restore the correct twofold coordination.
         if(terminate_octahedra) then
  
            j = cluster_info%N 
            do i=n_term%n,1,-1
               n = n_term%int(1,i) !Reference atom
               in_n = find_in_array(int_part(in_n_term,(/1/)),(/n/))
               ishift = n_term%int(2:4,i)
  
               if (n_term%int(5,i) .eq. 2) then
                 if (all(cluster_info%int(2:6,j) == (/ishift,1,n/))) then
                   bond1 = cluster_info%real(1:3,j)   - cluster_info%real(1:3,n)
                   bond2 = cluster_info%real(1:3,j-1) - cluster_info%real(1:3,n)
                   orig_length = sqrt(sum(bond1 * bond1))
                   
                   new_bond = bond1 + bond2 
                   new_bond = new_bond / sqrt(sum(new_bond*new_bond))
                   new_bond = new_bond * orig_length 
                   call delete(cluster_info, j) !Delete second H atom
                   call print('create_cluster: removed one of atom '//cluster_info%int(1,n)//" "//maxval(int_part(n_term,5))// &
                       ' terminating hydrogens ', PRINT_VERBOSE)
                   cluster_info%real(1:3,j-1) = cluster_info%real(1:3,n) + new_bond 
                 endif
                 j = j - 2
               elseif (n_term%int(5,i) .eq. 1 .and. in_n_term%int(5,in_n) .le. 1) then  !Do not remove the H atom
                 j = j - 1
               elseif (n_term%int(5,i) .eq. 1 .and. in_n_term%int(5,in_n) .gt. 1) then  !Remove the H atom
                 if (all(cluster_info%int(2:6,j) == (/ishift,1,n/))) then
                   call delete(cluster_info, j)
                   call print('create_cluster: removed one of atom '//cluster_info%int(1,n)//" "//maxval(int_part(n_term,5))// &
                       ' terminating hydrogens ', PRINT_VERBOSE)
                 endif
                 j = j - 1
               endif
            end do
  
         endif
  
         ! Do we need to remove a hydrogen atom to ensure equal n_up and n_down electrons?
         if (even_electrons .and. mod(sum(int_part(cluster_info,5)),2) == 1) then
  
            ! Find first atom with a maximal number of terminating hydrogens
  
            do i=1,n_term%n
               call append(sorted_n_term, (/n_term%int(5,i), n_term%int(1,i)/))
            end do
            allocate(idx(n_term%n))
            call sort(sorted_n_term, idx)
  
            n = sorted_n_term%int(2, sorted_n_term%n)
            most_hydrogens = idx(sorted_n_term%n)
            ishift = n_term%int(2:4,most_hydrogens)
  
            ! Loop over termination atoms
            do j=oldN,cluster_info%N
               ! Remove first H atom attached to atom i
               if (all(cluster_info%int(2:6,j) == (/ishift,1,n/))) then
                  call delete(cluster_info, j)
                  call print('create_cluster: removed one of atom '//cluster_info%int(1,n)//" "//maxval(int_part(n_term,5))// &
                       ' terminating hydrogens to zero total spin', PRINT_VERBOSE)
                  exit 
               end if
            end do
  
            deallocate(idx)
         end if
  
         call finalise(n_term)
  
         call print('create_cluster: Finished terminating cluster',PRINT_NERD)
      end if
      call system_timer('cluster_termination')
  
      call print ('Exiting create_cluster_info', PRINT_NERD)
  
      call finalise(cluster_list)
  
    end function create_cluster_info_from_mark
  
  
    !% Given an atoms structure with a 'hybrid_mark' property, set to 
    !% 'HYBRID_ACTIVE_MARK' for some atoms, this routine adds transition 
    !% and buffer regions, and creates a 'weight_region1' property, 
    !% whose values are between 0 and 1. 
    !%
    !% Atoms marked with 'HYBRID_ACTIVE_MARK' in 'hybrid_mark' get weight
    !% 1.  Neighbour hopping is done up to transition_hops times, over
    !% which the weight linearly decreases to zero with hop count if
    !% 'weight_interpolation=hop_ramp'.  If
    !% 'weight_interpolation=distance_ramp', weight is 1 up to
    !% 'distance_ramp_inner_radius', and goes linearly to 0 by
    !% 'distance_ramp_outer_radius', where distance is calculated from
    !% 'distance_ramp_center(1:3)'.  'distance_ramp' makes the most sense if
    !% the 'HYBRID_ACTIVE_MARK' is set on a sphere with radius
    !% distance_ramp_inner_radius from distance_ramp_center, no
    !% hysteresis, and with 'transition_hops' < 0 so hops continue until
    !% no more atoms are added.  Transition atoms are marked with
    !% 'HYBRID_TRANS_MARK'. If 'hysteretic_buffer=F' (default), buffer
    !% atoms are selected by hopping 'buffer_hops' times, otherwise
    !% hysteretically with radii 'hysteretic_buffer_inner_radius' and
    !% 'hysteretic_buffer_outer_radius'.  Buffer atoms are marked with
    !% 'HYBRID_BUFFER_MARK' and given weight 0.  If 'hysteretic_connect' is
    !% set, all hops are done on bonds in the hysteretic_connect
    !% structure.
  
    subroutine create_hybrid_weights(at, args_str, error)
      type(Atoms), intent(inout) :: at
      character(len=*), intent(in) :: args_str
      integer, optional, intent(out) :: error
  
      type(Dictionary) :: params
      logical :: has_distance_ramp_inner_radius, has_distance_ramp_outer_radius, has_distance_ramp_center
      real(dp) :: distance_ramp_inner_radius, distance_ramp_outer_radius, distance_ramp_center(3)
      logical :: min_images_only, mark_buffer_outer_layer, cluster_hopping_nneighb_only, hysteretic_buffer, hysteretic_connect, hysteretic_buffer_nneighb_only
      real(dp) :: hysteretic_buffer_inner_radius, hysteretic_buffer_outer_radius
      real(dp) :: hysteretic_connect_cluster_radius, hysteretic_connect_inner_factor, hysteretic_connect_outer_factor
      real(dp) :: hysteretic_connect_inner_cutoff, hysteretic_connect_outer_cutoff
      integer :: buffer_hops, transition_hops
      character(STRING_LENGTH) :: weight_interpolation, run_suffix
      logical :: construct_buffer_use_only_heavy_atoms, cluster_hopping_skip_unreachable
  
      integer, pointer :: hybrid_mark(:)
      real(dp), pointer :: weight_region1(:)
      integer :: n_region1, n_trans, n_region2 !, n_term
      integer :: i, j, jj, first_active, shift(3)
      logical :: dummy
      type(Table) :: activelist, currentlist, nextlist, distances, oldbuffer, bufferlist
      real(dp) :: core_CoM(3), core_mass, mass
      integer :: list_1, hybrid_number 
      type(Table) :: total_embedlist 
  
      real(dp) :: origin(3), extent(3,3)
      integer :: old_n
      integer, allocatable :: uniq_Z(:)
  
      logical :: distance_ramp, hop_ramp
      logical :: add_this_atom, more_hops
      integer :: cur_trans_hop
      real(dp) :: bond_len, d
      logical :: ignore_silica_residue
      integer :: res_num_silica
  
      INIT_ERROR(error)
  
  ! only one set of defaults now, not one in args_str and one in arg list 
      call initialise(params)
      call param_register(params, 'run_suffix', '', run_suffix, help_string="string to append to hyrbid_mark for proper mark property name")
      call param_register(params, 'transition_hops', '0', transition_hops, help_string="Number of hops to do transition region over")
      call param_register(params, 'buffer_hops', '3', buffer_hops, help_string="Number of hops for buffer region")
      call param_register(params, 'weight_interpolation', 'hop_ramp', weight_interpolation, help_string="How to ramp down weights for transition region: hop_ramp or distance_ramp")
      call param_register(params, 'distance_ramp_inner_radius', '0', distance_ramp_inner_radius, has_value_target=has_distance_ramp_inner_radius, help_string="If distance_ramp, inner radius to start reducing weight from 1")
      call param_register(params, 'distance_ramp_outer_radius', '0', distance_ramp_outer_radius, has_value_target=has_distance_ramp_outer_radius, help_string="If distance_ramp, outer radius by which to reach weight of 0")
      call param_register(params, 'distance_ramp_center', '0 0 0', distance_ramp_center, has_value_target=has_distance_ramp_center, help_string="If present, origin for distances of distance ramp (otherwise cluster center of mass)")
      call param_register(params, 'cluster_hopping_nneighb_only', 'T', cluster_hopping_nneighb_only, help_string="If true, only hop to atom pairs that are is_nearest_neighbor().")
      call param_register(params, 'cluster_hopping_skip_unreachable', 'F', cluster_hopping_skip_unreachable, help_string='If true, remove buffer atoms no longer reachable by bond hopping')
      call param_register(params, 'min_images_only', 'F', min_images_only, help_string="If true, consider only minimum images, not multiple periodic images")
      call param_register(params, 'mark_buffer_outer_layer', 'T', mark_buffer_outer_layer, help_string="If true, mark outermost buffer layer")
      call param_register(params, 'hysteretic_buffer', 'F', hysteretic_buffer, help_string="If true, do hysteretic buffer")
      call param_register(params, 'hysteretic_buffer_inner_radius', '5.0', hysteretic_buffer_inner_radius, help_string="Inner radius for hysteretic buffer")
      call param_register(params, 'hysteretic_buffer_outer_radius', '7.0', hysteretic_buffer_outer_radius, help_string="Outer radius for hysteretic buffer")
      call param_register(params, 'hysteretic_buffer_nneighb_only', 'F', hysteretic_buffer_nneighb_only, help_string="If true, only hop to nearest neighbour pairs when constructing hysteretic buffer")
      call param_register(params, 'hysteretic_connect', 'F', hysteretic_connect, help_string="If true, use hysteretic connect object")
      call param_register(params, 'hysteretic_connect_cluster_radius', '1.2', hysteretic_connect_cluster_radius, help_string="If using hysteretic connect, use as cluster radius for origin and extent estimates")
      call param_register(params, 'hysteretic_connect_inner_factor', '1.2', hysteretic_connect_inner_factor, help_string="Inner factor (multiplier for covalent bond distance) for hysteretic connect")
      call param_register(params, 'hysteretic_connect_outer_factor', '1.5', hysteretic_connect_outer_factor, help_string="Outer factor (multiplier for covalent bond distance) for hysteretic connect")
      call param_register(params, 'construct_buffer_use_only_heavy_atoms', 'F', construct_buffer_use_only_heavy_atoms, help_string="If true, use only non-H atoms for constructing buffer")
      call param_register(params, 'ignore_silica_residue', 'F', ignore_silica_residue, help_string="If true, do special things for silica") !lam81
      call param_register(params, 'res_num_silica', '1', res_num_silica, help_string="Residue number for silica") !lam81
      if (.not. param_read_line(params, args_str, ignore_unknown=.true.,task='create_hybrid_weights args_str') ) then
        RAISE_ERROR("create_hybrid_weights failed to parse args_str='"//trim(args_str)//"'", error)
      endif
      call finalise(params)
  
      hop_ramp = .false.
      distance_ramp = .false.
      if (trim(weight_interpolation) == 'hop_ramp') then
         hop_ramp = .true.
      else if (trim(weight_interpolation) == 'distance_ramp') then
         distance_ramp = .true.
      else
         RAISE_ERROR('create_hybrid_weights: unknown weight_interpolation value: '//trim(weight_interpolation), error)
      end if
  
      call print('create_hybrid_weights: transition_hops='//transition_hops//' buffer_hops='//buffer_hops//' weight_interpolation='//weight_interpolation, PRINT_VERBOSE)
      call print('  cluster_hopping_nneighb_only='//cluster_hopping_nneighb_only//'  min_images_only='//min_images_only//' mark_buffer_outer_layer='//mark_buffer_outer_layer, PRINT_VERBOSE)
      call print('  hysteretic_buffer='//hysteretic_buffer//' hysteretic_buffer_inner_radius='//hysteretic_buffer_inner_radius, PRINT_VERBOSE)
      call print('  hysteretic_buffer_outer_radius='//hysteretic_buffer_outer_radius//' hysteretic_buffer_nneighb_only='//hysteretic_buffer_nneighb_only, PRINT_VERBOSE)
      call print('  hysteretic_connect='//hysteretic_connect//' hysteretic_connect_cluster_radius='//hysteretic_connect_cluster_radius, PRINT_VERBOSE)
      call print('  hysteretic_connect_inner_factor='//hysteretic_connect_inner_factor //' hysteretic_connect_outer_factor='//hysteretic_connect_outer_factor, PRINT_VERBOSE)
  
  
      ! check to see if atoms has a 'weight_region1' property already, if so, check that it is compatible, if not present, add it
      if(assign_pointer(at, 'weight_region1'//trim(run_suffix), weight_region1)) then
         weight_region1 = 0.0_dp
      else
         call add_property(at, 'weight_region1'//trim(run_suffix), 0.0_dp)
         dummy = assign_pointer(at, 'weight_region1'//trim(run_suffix), weight_region1)
      end if
  
      ! check for a compatible hybrid_mark property. it must be present
      if(.not.assign_pointer(at, 'hybrid_mark'//trim(run_suffix), hybrid_mark)) then
         RAISE_ERROR('create_hybrid_weights: atoms structure has no "hybrid_mark'//trim(run_suffix)//'" property', error)
      endif
  
      ! Fast implementation of trivial case where buffer_hops=0 and transition_hops=0
      if (buffer_hops == 0 .and. transition_hops == 0) then
         call print('WARNING: buffer_hops=0 and transition_hops=0, aborting create_hybrid_weights() early')
  
         where (hybrid_mark /= HYBRID_ACTIVE_MARK)
            hybrid_mark = HYBRID_NO_MARK
         end where
  
         weight_region1 = 0.0_dp
         where (hybrid_mark == HYBRID_ACTIVE_MARK)
            weight_region1 = 1.0_dp
         end where
  
         return
      end if
  
      ! Add first marked atom to activelist. shifts will be relative to this atom
      first_active = find_in_array(hybrid_mark, HYBRID_ACTIVE_MARK)
  
      n_region1 = count(hybrid_mark == HYBRID_ACTIVE_MARK)
      list_1 = 0  
      call allocate(activelist, 4,0,0,0)
      call allocate(total_embedlist, 4,0,0,0) 
  
      do while (total_embedlist%N < n_region1)
  
         call wipe(currentlist)
         call wipe(activelist)
         call wipe(nextlist)    
  
         call append(activelist, (/first_active,0,0,0/))
         call append(currentlist, activelist)
         weight_region1(first_active) = 1.0_dp
  
         if (distance_ramp) then ! we might need ACTIVE_MARK center of mass
            if (has_property(at, 'mass')) then
               core_mass = at%mass(first_active)
            else
               core_mass = ElementMass(at%Z(first_active))
            end if
            core_CoM = core_mass*at%pos(:,first_active) ! we're taking this as reference atom so no shift needed here
         end if
  
         if (hysteretic_connect) then
           call system_timer('hysteretic_connect')
           call estimate_origin_extent(at, hybrid_mark == HYBRID_ACTIVE_MARK, hysteretic_connect_cluster_radius, origin, extent)
           call calc_connect_hysteretic(at, hysteretic_connect_inner_cutoff, hysteretic_connect_outer_cutoff, at%hysteretic_connect, origin, extent)
           call system_timer('hysteretic_connect')
         endif
  
         ! Add other active atoms using bond hopping from the first atom
         ! in the cluster to find the other active atoms and hence to determine the correct 
         ! periodic shifts
         
         hybrid_number = 1
         do while (hybrid_number .ne. 0)
            if (hysteretic_connect) then
              call BFS_step(at, currentlist, nextlist, nneighb_only=cluster_hopping_nneighb_only, min_images_only = min_images_only, alt_connect=at%hysteretic_connect)
            else
              call BFS_step(at, currentlist, nextlist, nneighb_only=cluster_hopping_nneighb_only, min_images_only = min_images_only, property=hybrid_mark)
            endif
            hybrid_number = 0 
            do j=1,nextlist%N
               jj = nextlist%int(1,j) 
               shift = nextlist%int(2:4,j)
               if (hybrid_mark(jj) == HYBRID_ACTIVE_MARK .and. find_in_array(activelist%int(1,1:activelist%N), jj) == 0) then
                  hybrid_number = hybrid_number+1 
                  call append(activelist, nextlist%int(:,j))
                  weight_region1(jj) = 1.0_dp
  
                  if (distance_ramp) then ! we might need ACTIVE_MARK center of mass
                     if (has_property(at, 'mass')) then
                        mass = at%mass(jj)
                     else
                        mass = ElementMass(at%Z(jj))
                     end if
                     core_CoM = core_CoM + mass*(at%pos(:,jj) + (at%lattice .mult. (shift)))
                     core_mass = core_mass + mass
                  end if
  
               end if
            end do
            call append(currentlist, nextlist)
         enddo
         list_1 = list_1 + activelist%N
         call append(total_embedlist, activelist)    
  
         if (distance_ramp) then ! calculate actual distance ramp parameters 
  	 ! distance_ramp center is as specified, otherwise ACTIVE_MARK center of mass
  	 core_CoM = core_CoM/core_mass
  	 if (.not. has_distance_ramp_center) distance_ramp_center = core_CoM
  
  	 ! distance_ramp_inner_radius is as specified, otherwise distance of atom furthest from distance_ramp_center
  	 if (.not. has_distance_ramp_inner_radius) then
  	   call initialise(distances, 1, 1, 0, 0)
  	   do i=1, activelist%N
  	     jj = activelist%int(1,i)
  	     call append(distances, jj, distance_min_image(at, distance_ramp_center, jj))
  	   end do
  	   distance_ramp_inner_radius = maxval(distances%real(1,1:distances%N))
  	   call finalise(distances)
  	 endif
  
  	 if (.not. has_distance_ramp_outer_radius) then
  	   distance_ramp_outer_radius = 0.0_dp
  	   call uniq(at%Z, uniq_Z)
  	   do i=1, size(uniq_Z)
  	   do j=1, size(uniq_Z)
  	     bond_len = bond_length(uniq_Z(i), uniq_Z(j))
  	     if (bond_len > distance_ramp_outer_radius) distance_ramp_outer_radius = bond_len 
  	   end do
  	   end do
  	   if (transition_hops <= 0) &
  	     call print("WARNING: using transition_hops for distance_ramp outer radius, but transition_hops="//transition_hops,PRINT_ALWAYS)
  	   distance_ramp_outer_radius = distance_ramp_inner_radius + distance_ramp_outer_radius*transition_hops
  	 endif
         endif ! distance_ramp
  
         call wipe(currentlist)
         call append(currentlist, activelist)
  
         ! create transition region
         call print('create_hybrid_mark: creating transition region',PRINT_VERBOSE)
         n_trans = 0
  
         cur_trans_hop = 1
         if (distance_ramp) transition_hops = -1 ! distance ramp always does as many hops as needed to get every atom within outer radius
         more_hops = (transition_hops < 0 .or. transition_hops >= 2)
         do while (more_hops)
  	 more_hops = .false.
           if (hysteretic_connect) then
             call BFS_step(at, currentlist, nextlist, nneighb_only = cluster_hopping_nneighb_only .and. (transition_hops > 0), min_images_only = min_images_only, alt_connect=at%hysteretic_connect)
           else
             call BFS_step(at, currentlist, nextlist, nneighb_only = cluster_hopping_nneighb_only .and. (transition_hops > 0), min_images_only = min_images_only)
           endif
  
           call wipe(currentlist)
           do j = 1,nextlist%N
              jj = nextlist%int(1,j)
              if(hybrid_mark(jj) == HYBRID_NO_MARK) then
  	       add_this_atom = .true.
  	       if (distance_ramp) then
                    d = distance_min_image(at, distance_ramp_center, jj)
  		  if (d >= distance_ramp_outer_radius) add_this_atom = .false.
  	       endif
  
  	       if (add_this_atom) then
  		 if (hop_ramp) weight_region1(jj) = 1.0_dp - real(cur_trans_hop,dp)/real(transition_hops,dp) ! linear transition
  		 if (distance_ramp) then        ! Save distance, weight will be calculated later
  		   if (d <= distance_ramp_inner_radius) then
  		     weight_region1(jj) = 1.0_dp
  		   else
  		     weight_region1(jj) = 1.0_dp - (d-distance_ramp_inner_radius)/(distance_ramp_outer_radius-distance_ramp_inner_radius)
  		   endif
  		 endif
  		 call append(currentlist, nextlist%int(:,j))
  		 hybrid_mark(jj) = HYBRID_TRANS_MARK
  		 n_trans = n_trans+1
  		 more_hops = .true.
  	       end if
              end if
           end do ! do j=1,nextlist%N
  
  	 cur_trans_hop = cur_trans_hop + 1
  	 if (transition_hops >= 0 .and. cur_trans_hop >= transition_hops) more_hops = .false.
         end do ! more_hops
  
         if (list_1 < n_region1) then
            call print('searching for a new quantum zone as found '//list_1//' atoms, need to get to '//n_region1, PRINT_VERBOSE)
            do i =1, at%N
               if (hybrid_mark(i) == HYBRID_ACTIVE_MARK .and. .not. Is_in_Array(total_embedlist%int(1,1:total_embedlist%N), i)) then
                  first_active = i
                  exit
               endif
            enddo
         endif
  
      enddo ! while (total_embedlist%N < n_region1)
  
      ! create region2 (buffer region) 
      if (.not. hysteretic_buffer) then
         ! since no hysteresis, safe to reset non ACTIVE/TRANS atoms to HYBRID_NO_MARK
         where (hybrid_mark /= HYBRID_ACTIVE_MARK .and. hybrid_mark /= HYBRID_TRANS_MARK) hybrid_mark = HYBRID_NO_MARK
  
         n_region2 = 0
         do i = 0,buffer_hops-1
  	  if (hysteretic_connect) then
  	    call BFS_step(at, currentlist, nextlist, nneighb_only = cluster_hopping_nneighb_only, min_images_only = min_images_only, alt_connect=at%hysteretic_connect)
  	  else
  	    call BFS_step(at, currentlist, nextlist, nneighb_only = cluster_hopping_nneighb_only, min_images_only = min_images_only)
  	  endif
            call wipe(currentlist)
            do j = 1,nextlist%N
               jj = nextlist%int(1,j)
               if(hybrid_mark(jj) == HYBRID_NO_MARK) then
                  call append(currentlist, nextlist%int(:,j))
                  weight_region1(jj) = 0.0_dp
                  if (i==buffer_hops-1 .and. mark_buffer_outer_layer) then
                     hybrid_mark(jj) = HYBRID_BUFFER_OUTER_LAYER_MARK
                  else
                     hybrid_mark(jj) = HYBRID_BUFFER_MARK
                  end if
                  n_region2 = n_region2+1
               end if
            end do
         end do
      else 
         ! hysteretic buffer here
  
         call initialise(oldbuffer, 1,0,0,0)
         !call wipe(oldbuffer)
         call append(oldbuffer, find(hybrid_mark /= HYBRID_NO_MARK))
  
         call wipe(currentlist)
  
         ! Add first marked atom to embedlist. shifts will be relative to this atom
         first_active = find_in_array(hybrid_mark, HYBRID_ACTIVE_MARK)
         call append(bufferlist, (/first_active,0,0,0/))
         call append(currentlist, bufferlist)
  
         n_region2 = count(hybrid_mark /= HYBRID_NO_MARK)
  
         ! Find old embed + buffer atoms using bond hopping from the first atom
         ! in the cluster to find the other buffer atoms and hence to determine the correct 
         ! periodic shifts
         !
         ! This will fail if marked atoms do not form connected clusters around the active atoms
         old_n = bufferlist%N
         do while (bufferlist%N < n_region2)
            if (hysteretic_connect) then
               call BFS_step(at, currentlist, nextlist, nneighb_only = cluster_hopping_nneighb_only, min_images_only = min_images_only, alt_connect=at%hysteretic_connect)
            else
               call BFS_step(at, currentlist, nextlist, nneighb_only = cluster_hopping_nneighb_only, min_images_only = min_images_only, property =hybrid_mark)
            endif
            do j=1,nextlist%N
               jj = nextlist%int(1,j)
               shift = nextlist%int(2:4,j)
               if (hybrid_mark(jj) /= HYBRID_NO_MARK) call append(bufferlist, nextlist%int(:,j))
            end do
            call append(currentlist, nextlist)
  
            ! check that cluster is still growing
            if (bufferlist%N == old_n) then
               exit
               if (cluster_hopping_skip_unreachable) then
                  call print('WARNING: skipping buffer atoms no longer connected')
                  exit
               else
                  RAISE_ERROR('create_hybrid_weights: buffer cluster stopped growing before all marked atoms found - check for split QM or buffer region', error)
               end if
            end if
            old_n = bufferlist%N
         end do
         n_region2 = bufferlist%N ! possibly excluding some no-longer-reachable atoms
  
         ! Remove marks on all buffer atoms
         do i=1,oldbuffer%N
            if (hybrid_mark(oldbuffer%int(1,i)) == HYBRID_BUFFER_MARK .or. &
                hybrid_mark(oldbuffer%int(1,i)) == HYBRID_BUFFER_OUTER_LAYER_MARK) &
                hybrid_mark(oldbuffer%int(1,i)) = HYBRID_NO_MARK
         end do
  
         !construct the hysteretic buffer region:
         if (hysteretic_connect) then
  	 call print("create_hybrid_weights calling construct_hysteretic_region", verbosity=PRINT_NERD)
  	 call construct_hysteretic_region(region=bufferlist,at=at,core=total_embedlist,loop_atoms_no_connectivity=.false., &
  	   inner_radius=hysteretic_buffer_inner_radius,outer_radius=hysteretic_buffer_outer_radius,use_avgpos=.false., &
  	   add_only_heavy_atoms=construct_buffer_use_only_heavy_atoms, cluster_hopping_nneighb_only=cluster_hopping_nneighb_only, &
             force_hopping_nneighb_only=hysteretic_buffer_nneighb_only, min_images_only=min_images_only, &
  	   alt_connect=at%hysteretic_connect, &
             ignore_silica_residue=ignore_silica_residue, res_num_silica=res_num_silica, error=error) !NB, debugfile=mainlog) lam81
         else
  	 call print("create_hybrid_weights calling construct_hysteretic_region", verbosity=PRINT_NERD)
  	 call construct_hysteretic_region(region=bufferlist,at=at,core=total_embedlist,loop_atoms_no_connectivity=.false., &
  	   inner_radius=hysteretic_buffer_inner_radius,outer_radius=hysteretic_buffer_outer_radius,use_avgpos=.false., &
  	   add_only_heavy_atoms=construct_buffer_use_only_heavy_atoms, cluster_hopping_nneighb_only=cluster_hopping_nneighb_only, &
             force_hopping_nneighb_only=hysteretic_buffer_nneighb_only, min_images_only=min_images_only, &
             ignore_silica_residue=ignore_silica_residue, res_num_silica=res_num_silica, error=error) !NB, debugfile=mainlog) lam81
         endif
  
         call print('bufferlist=',PRINT_VERBOSE)
         call print(bufferlist,PRINT_VERBOSE)
  
         ! Mark new buffer region, leaving core QM region alone
         ! at the moment only ACTIVE and NO marks are present, because hybrid_mark was set to =hybrid
         do i=1,bufferlist%N
            if (hybrid_mark(bufferlist%int(1,i)) == HYBRID_NO_MARK) & 
                 hybrid_mark(bufferlist%int(1,i)) = HYBRID_BUFFER_MARK
  
            ! Marking the outer layer with  HYBRID_BUFFER_OUTER_LAYER_MARK is
            ! dealt with in create_cluster_from_mark.
  
         end do
  
         call finalise(bufferlist)
         call finalise(oldbuffer)
  !       call finalise(embedlist)
  
      end if
  
      ! this is commented out for now, terminations are controlled in create_cluster only
      ! create terminations
      !n_term = 0
      !if(terminate) then
      !   call BFS_step(at, currentlist, nextlist, nneighb_only = .true., min_images_only = .true.)
      !   do j = 1,nextlist%N
      !      jj = nextlist%int(1,j)
      !      if(hybrid_mark(jj) == HYBRID_NO_MARK) then
      !         weight_region1(jj) = 0.0_dp
      !         hybrid_mark(jj) = HYBRID_TERM_MARK
      !         n_term = n_term+1
      !      end if
      !   end do
      !end do
  
      call print('create_hybrid_weights: '//list_1//' region 1, '//n_trans//' transition, '//n_region2//&
           ' region 2, '//count(hybrid_mark /= HYBRID_NO_MARK)//' in total', PRINT_VERBOSE)
      !call print('create_hybrid_weights: '//n_region1//' region 1, '//n_trans//' transition, '//n_region2//&
      !     ' region 2, '//count(hybrid_mark /= HYBRID_NO_MARK)//' in total', PRINT_VERBOSE)
      !    call sort(total_embedlist)
  
      call finalise(activelist)
      call finalise(currentlist)
      call finalise(nextlist)
      call finalise(distances)
      call finalise(total_embedlist)
  
    end subroutine create_hybrid_weights
  
    !% Return estimated `(origin, extent)` for hysteretic connectivity 
    !% calculator to include all atoms within a distance `cluster_radius`
    !% of an atom in the `active` array.
    subroutine estimate_origin_extent(at, active, cluster_radius, origin, extent)
      type(Atoms), intent(in) :: at
      logical, intent(in) :: active(:)
      real(dp), intent(in) :: cluster_radius
      real(dp), intent(out) :: origin(3), extent(3,3)
  
      real(dp) :: center(3), low_corner(3), high_corner(3), dr(3)
      integer :: i, n_active, first_active
      logical :: found_first_active
  
      found_first_active = .false.
      n_active = 0
      do i=1, at%N
        if (.not. active(i)) cycle
        n_active = n_active + 1
        if (found_first_active) then
  	center = center + diff_min_image(at, first_active, i)
        else
  	center = 0.0_dp
  	first_active = i
  	found_first_active = .true.
        endif
      end do
      center = center / real(n_active, dp)
      center = center + at%pos(:,first_active)
  
      call print("estimate_origin_extent: got center" // center, verbosity=PRINT_VERBOSE)
  
      low_corner = 1.0e38_dp
      high_corner = -1.0e38_dp
      do i=1, at%N
        if (.not. active(i)) cycle
        dr = diff_min_image(at, at%pos(:,i), center)
        low_corner = min(low_corner, dr)
        high_corner = max(high_corner, dr)
      end do
      call print("estimate_origin_extent: got relative low_corner" // low_corner, verbosity=PRINT_NERD)
      call print("estimate_origin_extent: got relative high_corner" // high_corner, verbosity=PRINT_NERD)
      low_corner = low_corner + center
      high_corner = high_corner + center
  
      call print("estimate_origin_extent: got low_corner" // low_corner, verbosity=PRINT_NERD)
      call print("estimate_origin_extent: got high_corner" // high_corner, verbosity=PRINT_NERD)
  
      origin = low_corner - cluster_radius
      extent = 0.0_dp
      extent(1,1) = (high_corner(1)-low_corner(1))+2.0_dp*cluster_radius
      extent(2,2) = (high_corner(2)-low_corner(2))+2.0_dp*cluster_radius
      extent(3,3) = (high_corner(3)-low_corner(3))+2.0_dp*cluster_radius
  
      call print("estimate_origin_extent: got origin" // origin, verbosity=PRINT_VERBOSE)
      call print("estimate_origin_extent: got extent(1,:) " // extent(1,:), verbosity=PRINT_VERBOSE)
      call print("estimate_origin_extent: got extent(2,:) " // extent(2,:), verbosity=PRINT_VERBOSE)
      call print("estimate_origin_extent: got extent(3,:) " // extent(3,:), verbosity=PRINT_VERBOSE)
  
    end subroutine estimate_origin_extent
  
    !% Given an Atoms structure with an active region marked in the mark_name (default 'hybrid_mark')
    !% property using 'HYBRID_ACTIVE_MARK', grow the embed region by 'fit_hops'
    !% bond hops to form a fit region. Returns the embedlist and fitlist with correct
    !% periodic shifts.
    subroutine create_embed_and_fit_lists(at, fit_hops, embedlist, fitlist, cluster_hopping_nneighb_only, min_images_only, mark_name, error)
  
      type(Atoms), intent(inout) :: at
      integer :: fit_hops
      type(Table), intent(out) :: embedlist, fitlist
      logical, intent(in), optional :: cluster_hopping_nneighb_only, min_images_only
      character(len=*), intent(in), optional :: mark_name
      integer, optional, intent(out) :: error
  
      character(len=255)                :: my_mark_name
      integer, pointer :: hybrid_mark(:)
      integer :: n_region1, n_region2 !, n_term
      integer :: i, j, jj, first_active, shift(3)
      logical :: do_hopping_nneighb_only, do_min_images_only
      type(Table) :: currentlist, nextlist, tmpfitlist
      integer :: n, hybrid_number
      type(Table) :: totallist
      integer :: list_1
  
      integer :: old_n
  
      INIT_ERROR(error)
  
      my_mark_name = optional_default('hybrid_mark', mark_name)
  
      call print('Entered create_embed_and_fit_lists.',PRINT_VERBOSE)
      do_hopping_nneighb_only = optional_default(.false., cluster_hopping_nneighb_only)
      do_min_images_only = optional_default(.true., min_images_only)
  
      ! check for a compatible hybrid_mark property. it must be present
      if(.not.assign_pointer(at, trim(my_mark_name), hybrid_mark)) then
         RAISE_ERROR('create_fit_region: atoms structure has no "'//trim(my_mark_name)//'"', error)
      endif
  
      call wipe(embedlist)
      call wipe(fitlist)
  
      ! Add first marked atom to embedlist. shifts will be relative to this atom
      first_active = find_in_array(hybrid_mark, HYBRID_ACTIVE_MARK)
  
      n_region1 = count(hybrid_mark == HYBRID_ACTIVE_MARK)
      list_1 = 0
      call allocate(totallist,4,0,0,0)
      call allocate(currentlist,4,0,0,0)
  
      ! Add other active atoms using bond hopping from the first atom
      ! in the cluster to find the other active atoms and hence to determine the correct 
      ! periodic shifts
      !
      old_n = embedlist%N
      do while (embedlist%N < n_region1)
  
        call wipe(currentlist)
        call wipe(nextlist)
        call wipe(totallist)
  
        call append(totallist, (/first_active,0,0,0/))
        call print('create_embed_and_fit_lists: expanding quantum region starting from '//first_active//' atom', PRINT_VERBOSE)
        call append(currentlist, totallist)
  
        hybrid_number = 1 
        do while (hybrid_number .ne. 0) 
  !!!CHECK THIS
         call BFS_step(at, currentlist, nextlist, nneighb_only = .false., min_images_only = do_min_images_only, property =hybrid_mark)
  
         hybrid_number = 0 
         do j=1,nextlist%N
            jj = nextlist%int(1,j)
            shift = nextlist%int(2:4,j)
            if (hybrid_mark(jj) == HYBRID_ACTIVE_MARK) then
                  hybrid_number = hybrid_number + 1 
                  call append(totallist, nextlist%int(:,j))
            endif
         end do
         call append(currentlist, nextlist)
        enddo
  
        list_1 = list_1 + totallist%N
        call append(embedlist, totallist)
        call print('create_embed_and_fit_lists: number of atoms in the embedlist is '//embedlist%N, PRINT_VERBOSE)
  
        if (list_1 .lt. n_region1) then
           n = 0
           do i =1, at%N
              if (hybrid_mark(i) == HYBRID_ACTIVE_MARK) then
                 n = n+1
                 if (.not. Is_in_Array(int_part(embedlist,1), i)) then
                    first_active = i
                    exit
                 endif
              endif
           enddo
        endif
  
         ! check that cluster is still growing
         !if (embedlist%N == old_n) then
         !     RAISE_ERROR('create_embed_and_fit_lists: (embedlist) cluster stopped growing before all marked atoms found - check for split QM region', error)
         !endif
         !old_n = embedlist%N
      end do
  
  
      call wipe(currentlist)
      call append(currentlist, embedlist)
  
  
      ! create region2 (fit region)
      call initialise(tmpfitlist,4,0,0,0,0)
      n_region2 = 0
      do i = 0,fit_hops-1
         call BFS_step(at, currentlist, nextlist, nneighb_only = do_hopping_nneighb_only, min_images_only = do_min_images_only)
         do j = 1,nextlist%N
            jj = nextlist%int(1,j)
            if(hybrid_mark(jj) /= HYBRID_ACTIVE_MARK) then
              call append(tmpfitlist, nextlist%int(:,j))
              n_region2 = n_region2+1
            end if
         end do
         call append(currentlist, nextlist)
      end do
  
      call print('create_embed_and_fit_lists: '//list_1//' embed, '//n_region2//' fit', PRINT_VERBOSE)
  
      ! Sort order to we are stable to changes in neighbour ordering introduced
      ! by calc_connect. 
      call sort(embedlist)
      call sort(tmpfitlist)
  
      ! fitlist consists of sorted embedlist followed by sorted list of remainder of fit atoms
      call append(fitlist, embedlist)
      call append(fitlist, tmpfitlist)
  
      if (do_min_images_only) then
         if (multiple_images(embedlist)) then
            call discard_non_min_images(embedlist)
            call print('create_embed_and_fits_lists: multiple images discarded from embedlist', PRINT_VERBOSE)
         endif
  
         if (multiple_images(fitlist)) then
            call discard_non_min_images(fitlist)
            call print('create_embed_and_fits_lists: multiple images discarded from fitlist', PRINT_VERBOSE)
         endif
      endif
  
      call finalise(currentlist)
      call finalise(nextlist)
      call finalise(tmpfitlist)
  
    end subroutine create_embed_and_fit_lists
  
    !% Given an Atoms structure with an active region marked in the mark_name (default 'cluster_mark')
    !% property using 'HYBRID_ACTIVE_MARK', and buffer region marked using 'HYBRID_BUFFER_MARK',
    !% 'HYBRID_TRANS_MARK' or 'HYBRID_BUFFER_OUTER_LAYER_MARK', simply returns the embedlist and fitlist
    !% according to 'cluster_mark'. It does not take into account periodic shifts.
    subroutine create_embed_and_fit_lists_from_cluster_mark(at,embedlist,fitlist, mark_name, error)
  
      type(Atoms), intent(in)  :: at
      type(Table), intent(out) :: embedlist, fitlist
      character(len=*), intent(in), optional :: mark_name
      integer, optional, intent(out) :: error
  
      character(len=255) :: my_mark_name
      type(Table)              :: tmpfitlist
  
      INIT_ERROR(error)
  
      my_mark_name = optional_default('cluster_mark', mark_name)
  
      call print('Entered create_embed_and_fit_lists_from_cluster_mark.',PRINT_VERBOSE)
      call wipe(embedlist)
      call wipe(fitlist)
      call wipe(tmpfitlist)
  
      !build embed list from ACTIVE atoms
      call list_matching_prop(at,embedlist,trim(my_mark_name),HYBRID_ACTIVE_MARK)
  
      !build fitlist from BUFFER and TRANS atoms
      call list_matching_prop(at,tmpfitlist,trim(my_mark_name),HYBRID_BUFFER_MARK)
      call append(fitlist,tmpfitlist)
      call list_matching_prop(at,tmpfitlist,trim(my_mark_name),HYBRID_TRANS_MARK)
      call append(fitlist,tmpfitlist)
      call list_matching_prop(at,tmpfitlist,trim(my_mark_name),HYBRID_BUFFER_OUTER_LAYER_MARK)
      call append(fitlist,tmpfitlist)
  
      call wipe(tmpfitlist)
      call append(tmpfitlist,fitlist)
  
      ! Sort in order to we are stable to changes in neighbour ordering introduced
      ! by calc_connect. 
      call sort(embedlist)
      call sort(tmpfitlist)
  
      ! fitlist consists of sorted embedlist followed by sorted list of remainder of fit atoms
      call wipe(fitlist)
      call append(fitlist, embedlist)
      call append(fitlist, tmpfitlist)
  
      call print('Embedlist:',PRINT_ANALYSIS)
      call print(int_part(embedlist,1),PRINT_ANALYSIS)
      call print('Fitlist:',PRINT_ANALYSIS)
      call print(int_part(fitlist,1),PRINT_ANALYSIS)
  
      call finalise(tmpfitlist)
      call print('Leaving create_embed_and_fit_lists_from_cluster_mark.',PRINT_VERBOSE)
  
    end subroutine create_embed_and_fit_lists_from_cluster_mark
  
    !% Return the atoms in a hysteretic region:
    !% To become part of the 'list' region, atoms must drift within the
    !% 'inner' list. To leave the 'list' region, an atom
    !% must drift further than 'outer'.
    !% Optionally use time averaged positions
    !
    subroutine update_hysteretic_region(at,inner,outer,list,min_images_only,error)
  
      type(Atoms),       intent(in)    :: at
      type(Table),       intent(inout) :: inner, outer
      type(Table),       intent(inout) :: list
      logical, optional, intent(in)    :: min_images_only
      integer, optional, intent(out)   :: error
  
      integer n, my_verbosity
      logical do_min_images_only
      integer, allocatable :: cols(:)
  
      INIT_ERROR(error)
  
      my_verbosity = current_verbosity()
      do_min_images_only = optional_default(.false., min_images_only)
      call print('update_hysteretic_region: do_min_images_only='//do_min_images_only, PRINT_VERBOSE)
  
      if (do_min_images_only) then
         ! only compare atom indices - first column in table
         allocate(cols(1))
         cols(1) = 1
      else
         ! compare atom indices and shifts - all four columns in table
         allocate(cols(4))
         cols(:) = (/1, 2, 3, 4/)
      end if
  
      if (my_verbosity == PRINT_ANALYSIS) then
         call print('In Select_Hysteretic_Quantum_Region:')
         call print('List currently contains '//list%N//' atoms')
      end if
  
      if ((list%intsize /= inner%intsize) .or. (list%intsize /= outer%intsize)) then
         RAISE_ERROR('update_hysteretic_region: inner, outer and list must have the same intsize', error)
      endif
  
      ! Speed up searching
      call sort(outer)
      call sort(list)
  
      !Check for atoms in 'list' and not in 'outer'
      do n = list%N, 1, -1
         if (search(outer,list%int(cols,n))==0) then
            call delete(list,n)
            if (my_verbosity > PRINT_NORMAL) call print('Removed atom ('//list%int(cols,n)//') from quantum list')
         end if
      end do
  
      call sort(list)
  
      !Check for new atoms in 'inner' cluster and add them to list
      do n = 1, inner%N
         if (search(list,inner%int(cols,n))==0) then
            call append(list,inner%int(:,n)) ! copy shifts columns as well, even if we're not comparing them
            call sort(list)
            if (my_verbosity > PRINT_NORMAL) call print('Added atom ('//inner%int(cols,n)//') to quantum list')
         end if
      end do
  
      if (my_verbosity >= PRINT_NORMAL) call print(list%N//' atoms selected for quantum treatment')
      deallocate(cols)
  
    end subroutine update_hysteretic_region
  
  
    !
    !% Given an atoms object, and a point 'centre' or a list of atoms 'core',
    !% fill the 'region' table with all atoms hysteretically within 'inner_radius -- outer_radius'
    !% of the 'centre' or any 'core' atom, which can be reached by connectivity hopping.
    !% In the case of building the 'region' around 'centre', simply_loop_over_atoms loops instead of
    !% bond hopping.
    !% Optionally use the time averaged positions.
    !% Optionally use only heavy atom selection.
    !
    subroutine construct_hysteretic_region(region,at,core,centre,loop_atoms_no_connectivity,inner_radius,outer_radius,use_avgpos,add_only_heavy_atoms,cluster_hopping_nneighb_only, &
                                           force_hopping_nneighb_only, min_images_only,alt_connect,debugfile, ignore_silica_residue, res_num_silica, error) !lam81
  
      type(Table),           intent(inout) :: region
      type(Atoms),           intent(in)    :: at
      type(Table), optional, intent(in)    :: core
      real(dp),    optional, intent(in)    :: centre(3)
      logical,     optional, intent(in)    :: loop_atoms_no_connectivity
      real(dp),              intent(in)    :: inner_radius
      real(dp),              intent(in)    :: outer_radius
      logical,     optional, intent(in)    :: use_avgpos
      logical,     optional, intent(in)    :: add_only_heavy_atoms
      logical,     optional, intent(in)  :: cluster_hopping_nneighb_only
      logical,     optional, intent(in)  :: force_hopping_nneighb_only
      logical,     optional, intent(in)  :: min_images_only
      type(Connection), optional,intent(in) :: alt_connect
      type(inoutput), optional :: debugfile
      integer, optional, intent(out) :: error
  
      type(Table)                          :: inner_region
      type(Table)                          :: outer_region
      logical                              :: no_hysteresis
  
      logical,     optional, intent(in)    :: ignore_silica_residue ! lam81
      integer,     optional, intent(in)    :: res_num_silica        ! lam81
  
      INIT_ERROR(error)
  
      if (present(debugfile)) call print("construct_hysteretic_region radii " // inner_radius // " " // outer_radius, file=debugfile)
      !check the input arguments only in construct_region.
      if (inner_radius.lt.0._dp) then
        RAISE_ERROR('inner_radius must be > 0 and it is '//inner_radius, error)
      endif
      if (outer_radius.lt.inner_radius) then
        RAISE_ERROR('outer radius ('//outer_radius//') must not be smaller than inner radius ('//inner_radius//').', error)
      endif
  
      no_hysteresis = .false.
      if ((outer_radius-inner_radius).lt.epsilon(0._dp)) no_hysteresis = .true.
      if (inner_radius .feq. 0.0_dp) call print('WARNING: hysteretic region with inner_radius .feq. 0.0', PRINT_ALWAYS)
      if (no_hysteresis) call print('WARNING! construct_hysteretic_region: inner_buffer=outer_buffer. no hysteresis applied: outer_region = inner_region used.',PRINT_ALWAYS)
      if (present(debugfile)) call print("   no_hysteresis " // no_hysteresis, file=debugfile)
  
      if (present(debugfile)) call print("   constructing inner region", file=debugfile)
      call construct_region(region=inner_region,at=at,core=core,centre=centre,loop_atoms_no_connectivity=loop_atoms_no_connectivity, &
        radius=inner_radius,use_avgpos=use_avgpos,add_only_heavy_atoms=add_only_heavy_atoms,cluster_hopping_nneighb_only=cluster_hopping_nneighb_only, &
        force_hopping_nneighb_only=force_hopping_nneighb_only, &
        min_images_only=min_images_only,alt_connect=alt_connect,debugfile=debugfile, &
        ignore_silica_residue=ignore_silica_residue, res_num_silica=res_num_silica) !lam81
      if (no_hysteresis) then
         call initialise(outer_region,4,0,0,0,0)
         call append(outer_region,inner_region)
      else
      if (present(debugfile)) call print("   constructing outer region", file=debugfile)
         call construct_region(region=outer_region,at=at,core=core,centre=centre,loop_atoms_no_connectivity=loop_atoms_no_connectivity, &
  	radius=outer_radius, use_avgpos=use_avgpos,add_only_heavy_atoms=add_only_heavy_atoms,cluster_hopping_nneighb_only=cluster_hopping_nneighb_only, &
          force_hopping_nneighb_only=force_hopping_nneighb_only, &
  	min_images_only=min_images_only,alt_connect=alt_connect,debugfile=debugfile, &
          ignore_silica_residue=ignore_silica_residue, res_num_silica=res_num_silica) ! lam81
      endif
  
      if (present(debugfile)) call print("   orig inner_region list", file=debugfile)
      if (present(debugfile)) call print(inner_region, file=debugfile)
      if (present(debugfile)) call print("   orig outer_region list", file=debugfile)
      if (present(debugfile)) call print(outer_region, file=debugfile)
      if (present(debugfile)) call print("   old region list", file=debugfile)
      if (present(debugfile)) call print(region, file=debugfile)
  
       call print('construct_hysteretic_region: old region list:',PRINT_VERBOSE)
       call print(region,PRINT_VERBOSE)
  
      call update_hysteretic_region(at,inner_region,outer_region,region, min_images_only=min_images_only)
  
      call print('construct_hysteretic_region: inner_region list:',PRINT_VERBOSE)
      call print(inner_region,PRINT_VERBOSE)
      call print('construct_hysteretic_region: outer_region list:',PRINT_VERBOSE)
      call print(outer_region,PRINT_VERBOSE)
  
      if (present(debugfile)) call print("   new inner region list", file=debugfile)
      if (present(debugfile)) call print(inner_region, file=debugfile)
      if (present(debugfile)) call print("   new outer region list", file=debugfile)
      if (present(debugfile)) call print(outer_region, file=debugfile)
  
      call print('construct_hysteretic_region: new region list:',PRINT_VERBOSE)
      call print(region,PRINT_VERBOSE)
  
      if (present(debugfile)) call print("   new region list", file=debugfile)
      if (present(debugfile)) call print(region, file=debugfile)
  
      call finalise(inner_region)
      call finalise(outer_region)
  
    end subroutine construct_hysteretic_region
  
    !
    !% Given an atoms object, and a point or a list of atoms in the first integer of
    !% the 'core' table, fill the 'buffer' table with all atoms within 'radius'
    !% of any core atom (which can be reached by connectivity hopping) or the point
    !%(with bond hopping or simply looping once over the atoms).
    !% Optionally use the time averaged positions (only for the radius).
    !% Optionally use a heavy atom based selection (applying to both the core and the region atoms).
    !% Alternatively use the hysteretic connection, only nearest neighbours and/or min_images (only for the n_connectivity_hops).
    !
    subroutine construct_region(region,at,core,centre,loop_atoms_no_connectivity,radius,n_connectivity_hops,use_avgpos,add_only_heavy_atoms,cluster_hopping_nneighb_only,&
                                force_hopping_nneighb_only, min_images_only,alt_connect,debugfile, ignore_silica_residue, res_num_silica, error) !lam81
  
      type(Table),           intent(out) :: region
      type(Atoms),           intent(in)  :: at
      type(Table), optional, intent(in)  :: core
      real(dp),    optional, intent(in)  :: centre(3)
      logical,     optional, intent(in)  :: loop_atoms_no_connectivity
      real(dp),    optional, intent(in)  :: radius
      integer,     optional, intent(in)  :: n_connectivity_hops
      logical,     optional, intent(in)  :: use_avgpos
      logical,     optional, intent(in)  :: add_only_heavy_atoms
      logical,     optional, intent(in)  :: cluster_hopping_nneighb_only
      logical,     optional, intent(in)  :: force_hopping_nneighb_only
      logical,     optional, intent(in)  :: min_images_only
      type(Connection), optional,intent(in) :: alt_connect
  type(inoutput), optional :: debugfile
      integer, optional, intent(out) :: error
  
      logical                             :: do_use_avgpos
      logical                             :: do_loop_atoms_no_connectivity
      logical                             :: do_add_only_heavy_atoms
      logical                             :: do_hopping_nneighb_only, do_force_hopping_nneighb_only
      logical                             :: do_min_images_only
      integer                             :: i, j, ii, ij
      type(Table)                         :: nextlist
      logical                             :: more_hops, add_i
      integer                             :: cur_hop
      real(dp), pointer                   :: use_pos(:,:)
      integer                             ::  shift_i(3)
      logical :: do_ignore_silica_residue, got_atom_res_number
  
      integer, pointer :: atom_res_number(:) !lam81
      logical,     optional, intent(in)  :: ignore_silica_residue !lam81
      integer,     optional, intent(in)  :: res_num_silica
  
      INIT_ERROR(error)
  
      do_loop_atoms_no_connectivity = optional_default(.false.,loop_atoms_no_connectivity)
      do_ignore_silica_residue = optional_default(.false., ignore_silica_residue)
  
      do_add_only_heavy_atoms = optional_default(.false.,add_only_heavy_atoms)
      if (do_add_only_heavy_atoms .and. .not. has_property(at,'Z')) then
        RAISE_ERROR("construct_region: atoms has no Z property", error)
      endif
  
      got_atom_res_number = assign_pointer(at, 'atom_res_number', atom_res_number)
  
      do_use_avgpos = optional_default(.false.,  use_avgpos)
  
      if (count((/ present(centre), present(core) /)) /= 1) then
        RAISE_ERROR("Need either centre or core, but not both present(centre) " // present(centre) // " present(core) "// present(core), error)
      endif
  
      ! if we have radius, we'll have to decide what positions to use
      if (present(radius)) then
        if (do_use_avgpos) then
          if (.not. assign_pointer(at, "avgpos", use_pos)) then
  	  RAISE_ERROR("do_use_avgpos is true, but no avgpos property", error)
  	endif
        else
          if (.not. assign_pointer(at, "pos", use_pos)) then
  	  RAISE_ERROR("do_use_avgpos is false, but no pos property", error)
  	endif
        endif
      endif
  
      call initialise(region,4,0,0,0,0)
  
      if (do_loop_atoms_no_connectivity) then
        if (present(debugfile)) call print("   loop over atoms", file=debugfile)
        if (.not. present(radius)) then
  	RAISE_ERROR("do_loop_atoms_no_connectivity=T requires radius", error)
        endif
        if (present(n_connectivity_hops) .or. present(min_images_only) .or. present(cluster_hopping_nneighb_only)) &
  	call print("WARNING: do_loop_atoms_no_connectivity, but specified unused arg n_connectivity_hops " // present(n_connectivity_hops) // &
  	  " min_images_only " // present(min_images_only) // " cluster_hopping_nneighb_only " // present(cluster_hopping_nneighb_only), PRINT_ALWAYS)
         ! jrk33: disabled this warning - should be fine with small cells since distance_min_image is used.
         ! call print('WARNING: check if your cell is greater than the radius, looping only works in that case.',PRINT_ALWAYS)
         !if (any((/at%lattice(1,1),at%lattice(2,2),at%lattice(3,3)/) < radius)) then
         !RAISE_ERROR('too small cell', error)
         !endifx
         do i = 1,at%N
  	 if (present(debugfile)) call print("check atom i " // i // " Z " // at%Z(i) // " pos " // at%pos(:,i), file=debugfile)
  	 if (do_add_only_heavy_atoms .and. at%Z(i) == 1) cycle
  	 if (present(centre)) then ! distance from centre is all that matters
  	   if (present(debugfile)) call print(" distance_min_image use_pos i " // use_pos(:,i) // " centre " // centre // " dist " // distance_min_image(at, use_pos(:,i), centre(1:3), shift_i(1:3)), file=debugfile)
  	   if (distance_min_image(at,use_pos(:,i),centre(1:3),shift_i(1:3)) < radius) then
  	     if (present(debugfile)) call print("   adding i " // i  // " at " // use_pos(:,i) // " center " // centre // " dist " // distance_min_image(at, use_pos(:,i), centre(1:3), shift_i(1:3)), file=debugfile)
  	     call append(region,(/i,shift_i(1:3)/))
  	    endif
  	  else ! no centre, check distance from each core list atom
  	    do ij=1, core%N
  	      j = core%int(1,ij)
  	      if (present(debugfile)) call print(" core atom ij " // ij // " j " // j // " distance_min_image use_pos i " // use_pos(:,i) // " use_pos j " // use_pos(:,j) // " dist " // distance_min_image(at, use_pos(:,i), use_pos(:,j), shift_i(1:3)), file=debugfile)
  	      if (distance_min_image(at,use_pos(:,i),use_pos(:,j),shift_i(1:3)) < radius) then
  		if (present(debugfile)) call print("   adding i " // i // " at " // use_pos(:,i) // " near j " // j // " at " // use_pos(:,j) // " dist " // distance_min_image(at, use_pos(:,i), use_pos(:,j), shift_i(1:3)), file=debugfile)
  	        call append(region,(/i,shift_i(1:3)/))
  		exit
  	      endif
  	    end do ! j
  	  endif ! no centre, use core
         enddo ! i
  
      else ! do_loop_atoms_no_connectivity
  
        if (.not. present(core)) then
  	RAISE_ERROR("do_loop_atoms_no_connectivity is false, trying connectivity hops, but no core list is specified", error)
        endif
  
        if (present(debugfile)) call print("   connectivity hopping", file=debugfile)
        if (present(debugfile) .and. present(radius)) call print("    have radius " // radius, file=debugfile)
        if (present(debugfile)) call print("   present cluster_hopping_nneighb_only " // present(cluster_hopping_nneighb_only), file=debugfile)
        if (present(debugfile) .and. present(cluster_hopping_nneighb_only)) call print("   cluster_hopping_nneighb_only " // cluster_hopping_nneighb_only, file=debugfile)
        do_hopping_nneighb_only = optional_default(.true., cluster_hopping_nneighb_only)
        do_force_hopping_nneighb_only = optional_default(.false., force_hopping_nneighb_only)
        do_min_images_only = optional_default(.true., min_images_only)
        if (do_use_avgpos) then
  	RAISE_ERROR("can't use avgpos with connectivity hops - make sure your connectivity is based on pos instead", error)
        endif
  
        if (do_ignore_silica_residue) then
           call print_warning('Overriding min_images_only since ignore_silica_residue=T')
           do_min_images_only = .true. !lam81 
        end if
        ! start with core
        call append(region,core)
  
        more_hops = .true.
        cur_hop = 1
        do while (more_hops)
  	if (present(debugfile)) call print('   construct_region do_hopping_nneighb_only = ' // do_hopping_nneighb_only // ' do_min_images_only = '//do_min_images_only, file=debugfile)
  	if (present(debugfile)) call print('   doing hop ' // cur_hop, file=debugfile)
  	if (present(debugfile)) call print('   cutoffs ' // at%cutoff, file=debugfile)
  	more_hops = .false.
  	call bfs_step(at, region, nextlist, nneighb_only=(do_hopping_nneighb_only .and. .not. present(radius)) .or. do_force_hopping_nneighb_only, &
                        min_images_only=do_min_images_only, max_r=radius, alt_connect=alt_connect, debugfile=debugfile)
  	if (present(debugfile)) call print("   bfs_step returned nextlist%N " // nextlist%N, file=debugfile)
  	if (present(debugfile)) call print(nextlist, file=debugfile)
  	if (nextlist%N /= 0) then ! go over things in next hop
  	  do ii=1, nextlist%N
  	    i = nextlist%int(1,ii)
  	    add_i = .true.
  	    if (do_add_only_heavy_atoms) then
  	      if (at%Z(i) == 1) cycle
  	    endif
              if (do_ignore_silica_residue .and. got_atom_res_number) then     ! lam81
                 if (atom_res_number(i) == res_num_silica) cycle ! lam81
              endif
  	    if (present(debugfile)) call print("  i " // i // " is heavy or all atoms requested", file=debugfile)
  	    if (present(radius)) then ! check to make sure we're close enough to core list
  	      if (present(debugfile)) call print("  radius is present, check distance from core list", file=debugfile)
  	      add_i = .false.
  	      do ij=1, core%N
  		j = core%int(1,ij)
  		if (present(debugfile)) call print("  distance of ii " // ii // " i " // i // " at  " // use_pos(:,i) // " from ij " // ij // " j " // j // " at " // use_pos(:,j) // " is " // distance_min_image(at,use_pos(:,i),use_pos(:,j),shift_i(1:3)), file=debugfile)
  		if (distance_min_image(at,use_pos(:,i),use_pos(:,j),shift_i(1:3)) <= radius) then
  		  if (present(debugfile)) call print("  decide to add i", file=debugfile)
  		  add_i = .true.
  		  exit ! from ij loop
  		endif
  	      end do ! ij
  	    endif ! present(radius)
  	    if (add_i) then
  	      if (present(debugfile)) call print("  actually adding i", file=debugfile)
  	      more_hops = .true.
  	      call append(region, nextlist%int(1:4,ii))
  	    endif
  	  end do
  	else ! nextlist%N == 0
  	  if (present(n_connectivity_hops)) then
  	    if (n_connectivity_hops > 0) &
  	      call print("WARNING: n_connectivity_hops = " // n_connectivity_hops // " cur_hop " // cur_hop // " but no atoms added")
  	  endif
  	endif ! nextlist%N
  
  	cur_hop = cur_hop + 1
  	if (present(n_connectivity_hops)) then
  	  if (cur_hop > n_connectivity_hops) more_hops = .false.
  	endif
        end do ! while more_hops
  
      endif ! do_loop_atoms_no_connectivity
      if (present(debugfile)) call print("leaving construct_region()", file=debugfile)
  
    end subroutine construct_region
  
  
    subroutine update_active(this, nneightol, avgpos, reset, error)
      type(Atoms) :: this
      real(dp), optional :: nneightol
      logical, optional :: avgpos, reset
      integer, optional, intent(out) :: error
  
      type(Atoms), save :: nn_atoms
      integer, pointer, dimension(:) :: nn, old_nn, active
      integer  :: i
      real(dp) :: use_nn_tol
      logical  :: use_avgpos, do_reset
  
      INIT_ERROR(error)
  
      use_nn_tol = optional_default(this%nneightol, nneightol)
      use_avgpos = optional_default(.true., avgpos)
      do_reset   = optional_default(.false., reset)
  
      ! First time copy entire atoms structure into nn_atoms and
      ! add properties for nn, old_nn and active
      if (reset .or. .not. is_initialised(nn_atoms)) then
         nn_atoms = this
         call add_property(this, 'nn', 0)
         call add_property(this, 'old_nn', 0)
         call add_property(this, 'active', 0)
  
  
      end if
  
      if (this%N /= nn_atoms%N) then
           RAISE_ERROR('update_actives: Number mismatch between this%N ('//this%N//') and nn_atoms%N ('//nn_atoms%N//')', error)
      endif
  
      if (.not. assign_pointer(this, 'nn', nn)) then
           RAISE_ERROR('update_actives: Atoms is missing "nn" property', error)
      endif
  
      if (.not. assign_pointer(this, 'old_nn', old_nn)) then
           RAISE_ERROR('update_actives: Atoms is missing "old_nn" property', error)
      endif
  
      if (.not. assign_pointer(this, 'active', active)) then
           RAISE_ERROR('update_actives: Atoms is missing "active" property', error)
      endif
  
      call print('update_actives: recalculating nearest neighbour table', PRINT_VERBOSE)
      if (use_avgpos .and. associated(this%avgpos)) then
         call print('update_actives: using time averaged atomic positions', PRINT_VERBOSE)
         nn_atoms%pos = this%avgpos
      else
         call print('update_actives: using instantaneous atomic positions', PRINT_VERBOSE)
         nn_atoms%pos = this%pos
      end if
      !% Use hysteretic connect so we can work with relative cutoffs
      call calc_connect_hysteretic(nn_atoms, cutoff_factor=use_nn_tol, &
                                   cutoff_break_factor=use_nn_tol)
  
      nn = 0
      do i = 1,nn_atoms%N
         nn(i) = n_neighbours(nn_atoms, i)
      end do
  
      if (all(old_nn == 0)) old_nn = nn ! Special case for first time
  
      ! Decrement the active counts
      where (active /= 0) active = active -1
  
      ! Find newly active atoms
      where (nn /= old_nn) active = 1
  
      if (count(active == 1) /= 0) then
         call print('update_actives: '//count(active == 1)//' atoms have become active.')
      end if
  
      old_nn = nn
  
    end subroutine update_active
  
  
    !
    !% Given an atoms structure and a list of quantum atoms, find X-H
    !% bonds which have been cut and include the other atom of 
    !% the pair in the quantum list.
    !
    subroutine add_cut_hydrogens(this,qmlist,heuristics_nneighb_only,verbosity,alt_connect)
  
      type(Atoms),       intent(in),          target :: this
      type(Table),       intent(inout)               :: qmlist
      integer, optional, intent(in)                  :: verbosity
      logical, optional, intent(in)                  :: heuristics_nneighb_only
      type(Connection), intent(in), optional, target :: alt_connect
  
      type(Table)                :: neighbours, bonds, centre
      logical                    :: more_atoms
      integer                    :: i, j, n, nn, added
      type(Connection), pointer :: use_connect
      logical :: my_heuristics_nneighb_only
  
      ! Check for atomic connectivity
      if (present(alt_connect)) then
        use_connect => alt_connect
      else
        use_connect => this%connect
      endif
  
      my_heuristics_nneighb_only = optional_default(.true., heuristics_nneighb_only)
  
      more_atoms = .true.
      added = 0
      call allocate(centre,4,0,0,0,1)
  
      !Repeat while the search trigger is true
      do while(more_atoms)
  
         more_atoms = .false.
  
         !Find nearest neighbours of the cluster
         call bfs_step(this,qmlist,neighbours,nneighb_only=my_heuristics_nneighb_only,min_images_only=.true.,alt_connect=use_connect)
  
         !Loop over neighbours
         do n = 1, neighbours%N
  
            i = neighbours%int(1,n)
  
            call wipe(centre)
            call append(centre,(/i,0,0,0/))
  
            ! Find atoms bonded to this neighbour
            call bfs_step(this,centre,bonds,nneighb_only=my_heuristics_nneighb_only,min_images_only=.true.,alt_connect=use_connect)
  
            !Loop over these bonds
            do nn = 1, bonds%N
  
               j = bonds%int(1,nn)
  
               !Try to find j in the qmlist
               if (find_in_array(int_part(qmlist,1),j)/=0) then
                  !We have a cut bond. If i or j are hydrogen then add i to the 
                  !quantum list and trigger another search after this one
                  if (this%Z(i)==1 .or. this%Z(j)==1) then
                     call append(qmlist,(/i,0,0,0/))
                     if (present(verbosity)) then
                        if (verbosity >= PRINT_NORMAL) call print('Add_Cut_Hydrogens: Added atom '//i//', neighbour of atom '//j)
                     end if
                     more_atoms = .true.
                     added = added + 1
                     exit !Don't add the same atom more than once
                  end if
               end if
  
            end do
  
         end do
  
      end do
  
      !Report findings
      if (present(verbosity)) then
         if (verbosity >= PRINT_NORMAL) then
            write(line,'(a,i0,a)')'Add_Cut_Hydrogens: Added ',added,' atoms to quantum region'
            call print(line)
         end if
      end if
  
      call finalise(centre)
      call finalise(bonds)
      call finalise(neighbours)
  
    end subroutine add_cut_hydrogens
  
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  !X
  !X GROUP CONVERSION
  !X
  !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
  
    !
    !% Convert a constrained atom's group to a normal group prior to quantum treatment
    !% The constraints are not deleted, just ignored. The number of degrees of freedom
    !% of the system is also updated.
    !
    subroutine constrained_to_quantum(this,i)
  
      type(DynamicalSystem), intent(inout) :: this
      integer,               intent(in)    :: i
      integer                              :: g!, n, j, a1, a2
  
      g = this%group_lookup(i)
  
      if (this%group(g)%type == TYPE_CONSTRAINED) then
         write(line,'(2(a,i0),a)')' INFO: Converting group ',g,' from Constrained to Quantum'
         call print(line)
  
         !Change the group type
         call set_type(this%group(g),TYPE_ATOM)
  
         !Update the number of degrees of freedom (add 1 per removed constraint)
         this%Ndof = this%Ndof + Group_N_Objects(this%group(g))
  
         !Thermalize the released constraints
  !       do n = 1, group_n_objects(this%group(g))
  !
  !          j = group_nth_object(this%group(g),n)
  !
  !          !If not a two body constraint then go to the next one.
  !          if (this%constraint(j)%N /= 2) cycle
  !
  !          !Get the indices of the atoms in this constraint
  !          a1 = this%constraint(j)%atom(1)
  !          a2 = this%constraint(j)%atom(2)
  !
  !          call thermalize_bond(this,a1,a2,this%sim_temp)
  !
  !       end do
  
  !       call thermalize_group(this,g,this%sim_temp)
  
      end if
  
    end subroutine constrained_to_quantum
  
    !
    !% When a constraint is released as a molecule drifts into the quantum region the
    !% new, unconstrained degree of freedom has zero temperature. This routine adds
    !% velocities in the directions between atoms which make up two body constraints
    !% at a specified temperature
    !
    subroutine thermalize_bond(this,i,j,T)
  
      type(DynamicalSystem), intent(inout) :: this         !% The dynamical system
      integer,               intent(in)    :: i,j          !% The two atoms
      real(dp),              intent(in)    :: T            !% Temperature
      real(dp)                             :: meff         !  effective mass
      real(dp), dimension(3)               :: u_hat        !  unit vector between the atoms
      real(dp)                             :: w, wi, wj, & ! magnitude of velocity component in bond direction
                                              mi, mj       ! masses of the atoms
      real(dp)                             :: E0,E1        ! kinetic energy in bond before and after
      real(dp), save                       :: Etot = 0.0_dp ! cumulative energy added to the system
      real(dp), save                       :: Eres = 0.0_dp ! cumulative energy already in bonds (QtoC error)
      real(dp)                             :: r, v         ! relative distance and velocity
  
      call print('Thermalize_Bond: Thermalizing bond '//i//'--'//j)
  
      ! Find unit vector along bond and current relative velocity
      call distance_relative_velocity(this%atoms,i,j,r,v)
  
      ! Find unit vector between atoms
      u_hat = diff_min_image(this%atoms,i,j)
      u_hat = u_hat / norm(u_hat)
  
      ! Calculate effective mass
      mi = ElementMass(this%atoms%Z(i))
      mj = ElementMass(this%atoms%Z(j))
      meff = mi * mj / (mi + mj)
  
      ! Calculate energy currently stored in this bond (should be zero)
      E0 = 0.5_dp * meff * v * v
      Eres = Eres + E0
  
      call print('Thermalize_Bond: Bond currently contains '//round(E0,8)//' eV')
      call print('Thermalize_Bond: Residual energy encountered so far: '//round(Eres,8)//' eV')
  
      ! Draw a velocity for this bond
      w = gaussian_velocity_component(meff,T)
      w = w - v ! Adjust for any velocity which the bond already had
  
      ! Distribute this between the two atoms, conserving momentum
      wi = (meff / mi) * w
      wj = -(meff / mj) * w
  
      this%atoms%velo(:,i) = this%atoms%velo(:,i) + wi * u_hat
      this%atoms%velo(:,j) = this%atoms%velo(:,j) + wj * u_hat
  
      call distance_relative_velocity(this%atoms,i,j,r,v)
      E1 = 0.5_dp * meff * v * v
  
      Etot = Etot + E1 - E0
  
      call print('Thermalize_Bond: Added '//round(E1-E0,8)//' eV to bond')
      call print('Thermalize_Bond: Added '//round(Etot,8)//' eV to the system so far')
  
    end subroutine thermalize_bond
  
    !
    !% Convert a quantum atom's group to a constrained group. The group is searched for
    !% two body constraints, which are assumed to be bond length or time dependent bond length
    !% constraints. The required final length is extracted, and a new time dependent bond length
    !% constraint is imposed, taking the current length and relative velocity along the bond
    !% to the required values, in the time given by smoothing_time.
    !%
    !% Any other constraints are left untouched.
    !
    subroutine quantum_to_constrained(this,i,smoothing_time)
  
      type(DynamicalSystem), intent(inout) :: this
      integer,               intent(in)    :: i
      real(dp),              intent(in)    :: smoothing_time
      integer, save                        :: CUBIC_BONDLENGTH_SQ_FUNC
      logical, save                        :: first_call = .true.
      integer                              :: g, j, n, a1, a2, datalength
      real(dp)                             :: x0,y0,y0p,x1,y1,y1p, coeffs(4), t, dE, m1, m2, meff
      real(dp), save                       :: Etot = 0.0_dp
  
      !Register the constraint function if this is the first call
      if (first_call) then
         CUBIC_BONDLENGTH_SQ_FUNC = Register_Constraint(CUBIC_BONDLENGTH_SQ)
         first_call = .false.
      end if
  
      g = this%group_lookup(i)
  
      if (this%group(g)%type == TYPE_ATOM) then
         write(line,'(2(a,i0),a)')' INFO: Converting group ',g,' from Quantum to Constrained'
         call print(line)
  
         !Change the group type
         call set_type(this%group(g),TYPE_CONSTRAINED)
  
         !Update the number of degrees of freedom (subtract 1 per added constraint)
         this%Ndof = this%Ndof - Group_N_Objects(this%group(g))
  
         !Loop over constraints
         do n = 1, Group_N_Objects(this%group(g))
  
            j = Group_Nth_Object(this%group(g),n)
  
            !If not a two body constraint then go to the next one.
            if (this%constraint(j)%N /= 2) cycle
  
            !Get the indices of the atoms in this constraint
            a1 = this%constraint(j)%atom(1)
            a2 = this%constraint(j)%atom(2)
  
            !Detect the type of constraint (simple bondlength or time dependent bond length) 
            !from the size of the data array (1 or 6)
            datalength = size(this%constraint(j)%data)
            select case(datalength)
            case(1) 
               !Simple bond length: the required final length is the only entry
               y1 = this%constraint(j)%data(1)
            case(6)
               !Time dependent bond length: the final length is the polynomial evaluated at the end point
               coeffs = this%constraint(j)%data(1:4)
               t = this%constraint(j)%data(6)
               y1 = ((coeffs(1)*t + coeffs(2))*t + coeffs(3))*t + coeffs(4)
            case default
               !This isn't a known bond length constraint. Set y1 to a negative value.
               y1 = -1.0_dp
            end select
  
            !If an existing bond length constraint has been found then set up a time dependent constraint
            if (y1 > 0.0_dp) then
  
               call print('Quantum_to_Constrained: Reconstraining bond '//a1//'--'//a2)
  
               x0 = this%t
               call distance_relative_velocity(this%atoms,a1,a2,y0,y0p)
               x1 = this%t + smoothing_time
               y1p = 0.0_dp
  
               !Find the coefficients of the cubic which fits these boundary conditions
               call fit_cubic(x0,y0,y0p,x1,y1,y1p,coeffs)
  
               !Change the constraint
               call ds_amend_constraint(this,j,CUBIC_BONDLENGTH_SQ_FUNC,(/coeffs,x0,x1/))
  
               !Work out the amount of energy which will be removed
               m1 = ElementMass(this%atoms%Z(a1))
               m2 = ElementMass(this%atoms%Z(a2))
               meff = m1 * m2 / (m1 + m2)
  
               dE = 0.5_dp * meff * y0p * y0p
               Etot = Etot + dE
  
               call print('Quantum_to_Constrained: Gradually removing '//round(dE,8)//' eV from bond')
               call print('Quantum_to_Constrained: Will have removed approximately '//round(Etot,8)//' eV in total')
  
            end if
  
         end do
  
      end if
  
    end subroutine quantum_to_constrained
  
  
    subroutine thermalize_group(this,g,T)
  
      type(DynamicalSystem), intent(inout) :: this
      integer,               intent(in)    :: g
      real(dp),              intent(in)    :: T
  
      integer  :: i, j, k
      real(dp) :: mi, mj, mk, mij, mik, mjk, vij, vik, vjk, rij, rik, rjk, gij, gik, gjk
      real(dp) :: ri(3), rj(3), rk(3), Rcom(3)
      real(dp) :: A(9,9), b(9), A_inv(9,9), x(9)
  
      real(dp) :: dEij, dEik, dEjk
      real(dp), save :: Etot = 0.0_dp
  
      ! Get atoms
      i = this%group(g)%atom(1)
      j = this%group(g)%atom(2)
      k = this%group(g)%atom(3)
  
      ! Get shifted positions
      ri = 0.0_dp
      rj = diff_min_image(this%atoms,i,j)
      rk = diff_min_image(this%atoms,i,k)  
  
      ! Get separations and current bond velocities
      call distance_relative_velocity(this%atoms,i,j,rij,vij)
      call distance_relative_velocity(this%atoms,i,k,rik,vik)
      call distance_relative_velocity(this%atoms,j,k,rjk,vjk)
  
      ! Calculate masses
      mi = ElementMass(this%atoms%Z(i))
      mj = ElementMass(this%atoms%Z(j))
      mk = ElementMass(this%atoms%Z(k))   
      mij = mi * mj / (mi + mj)
      mik = mi * mk / (mi + mk)
      mjk = mj * mk / (mj + mk)
  
      ! Calculate centre of mass
      Rcom = (mi*ri + mj*rj + mk*rk) / (mi+mj+mk)
  
      ! Select a gaussian distributed velocity adjustment for each bond
      gij = gaussian_velocity_component(mij,T)
      gik = gaussian_velocity_component(mik,T)
      gjk = gaussian_velocity_component(mjk,T)
      vij = gij - vij
      vik = gik - vik
      vjk = gjk - vjk
  
      ! Build matrix...
      A = 0.0_dp
      b = 0.0_dp
  
      ! Momentum conservation
      A(1,1) = mi; A(2,2) = mi; A(3,3) = mi
      A(1,4) = mj; A(2,5) = mj; A(3,6) = mj
      A(1,7) = mk; A(2,8) = mk; A(3,9) = mk
  
      ! Angular momentum conservation
      A(4,2) = -mi * (ri(3) - Rcom(3)); A(4,3) = mi * (ri(2) - Rcom(2))
      A(4,5) = -mj * (rj(3) - Rcom(3)); A(4,6) = mj * (rj(2) - Rcom(2))
      A(4,8) = -mk * (rk(3) - Rcom(3)); A(4,9) = mk * (rk(2) - Rcom(2))
  
      A(5,1) = mi * (ri(3) - Rcom(3)); A(5,3) = -mi * (ri(1) - Rcom(1))
      A(5,4) = mj * (rj(3) - Rcom(3)); A(5,6) = -mj * (rj(1) - Rcom(1))
      A(5,7) = mk * (rk(3) - Rcom(3)); A(5,9) = -mk * (rk(1) - Rcom(1))
  
      A(6,1) = -mi * (ri(2) - Rcom(2)); A(6,2) = mi * (ri(1) - Rcom(1))
      A(6,4) = -mj * (rj(2) - Rcom(2)); A(6,5) = mj * (rj(1) - Rcom(1))
      A(6,7) = -mk * (rk(2) - Rcom(2)); A(6,8) = mk * (rk(1) - Rcom(1))
  
      ! Bond length velocities
      A(7,1) = ri(1) - rj(1); A(7,2) = ri(2) - rj(2); A(7,3) = ri(3) - rj(3)
      A(7,4) = rj(1) - ri(1); A(7,5) = rj(2) - ri(2); A(7,6) = rj(3) - ri(3)
      b(7) = rij * vij
  
      A(8,1) = ri(1) - rk(1); A(8,2) = ri(2) - rk(2); A(8,3) = ri(3) - rk(3)
      A(8,7) = rk(1) - ri(1); A(8,8) = rk(2) - ri(2); A(8,9) = rk(3) - ri(3)
      b(8) = rik * vik
  
      A(9,4) = rj(1) - rk(1); A(9,5) = rj(2) - rk(2); A(9,6) = rj(3) - rk(3)
      A(9,7) = rk(1) - rj(1); A(9,8) = rk(2) - rj(2); A(9,9) = rk(3) - rj(3)
      b(9) = rjk * vjk
  
      ! Invert
      call inverse(A,A_inv)
      x = A_inv .mult. b
  
      call print('Momentum before update = '//momentum(this,(/i,j,k/)))
  
      call print('Angular momentum before update = '//&
           (mi*((ri-Rcom) .cross. this%atoms%velo(:,i)) + &
            mj*((rj-Rcom) .cross. this%atoms%velo(:,j)) + &
            mk*((rk-Rcom) .cross. this%atoms%velo(:,k))))
  
      dEij = -bond_energy(this,i,j)
      dEik = -bond_energy(this,i,k)
      dEjk = -bond_energy(this,j,k)
  
      this%atoms%velo(:,i) = this%atoms%velo(:,i) + x(1:3)
      this%atoms%velo(:,j) = this%atoms%velo(:,j) + x(4:6)
      this%atoms%velo(:,k) = this%atoms%velo(:,k) + x(7:9)
  
      dEij = dEij + bond_energy(this,i,j)
      dEik = dEik + bond_energy(this,i,k)
      dEjk = dEjk + bond_energy(this,j,k)
  
      call print('Momentum after update = '//momentum(this,(/i,j,k/)))
  
      call print('Angular momentum after update = '//&
           (mi*((ri-Rcom) .cross. this%atoms%velo(:,i)) + &
            mj*((rj-Rcom) .cross. this%atoms%velo(:,j)) + &
            mk*((rk-Rcom) .cross. this%atoms%velo(:,k))))
  
      ! Now check the result
  
      call print('Required velocity for bond '//i//'--'//j//': '//gij)
      call print('Required velocity for bond '//i//'--'//k//': '//gik)
      call print('Required velocity for bond '//j//'--'//k//': '//gjk)
  
      call distance_relative_velocity(this%atoms,i,j,rij,vij)
      call distance_relative_velocity(this%atoms,i,k,rik,vik)
      call distance_relative_velocity(this%atoms,j,k,rjk,vjk)
  
      call print('')
  
      call print('Actual velocity for bond '//i//'--'//j//': '//vij)
      call print('Actual velocity for bond '//i//'--'//k//': '//vik)
      call print('Actual velocity for bond '//j//'--'//k//': '//vjk)
  
      call print('Energy added to bond '//i//'--'//j//': '//dEij)
      call print('Energy added to bond '//i//'--'//k//': '//dEik)
      call print('Energy added to bond '//j//'--'//k//': '//dEjk)
  
      Etot = Etot + dEij + dEik + dEjk
  
      call print('Total energy added so far = '//Etot)
  
    end subroutine thermalize_group
  
    function bond_energy(this,i,j)
  
      type(DynamicalSystem), intent(in) :: this
      integer,               intent(in) :: i, j
      real(dp)                          :: bond_energy
  
      real(dp)                          :: mi, mj, vij, rij
  
      mi = ElementMass(this%atoms%Z(i))
      mj = ElementMass(this%atoms%Z(j))
  
      call distance_relative_velocity(this%atoms,i,j,rij,vij)
  
      bond_energy = 0.5_dp * mi * mj * vij * vij / (mi + mj)
  
    end function bond_energy
  
    !% Updates the core QM flags saved in 'hybrid' and 'hybrid_mark'
    !% properties.  Do this hysteretically, from 'R_inner' to 'R_outer'
    !% around 'origin' or 'atomlist', that is the centre of the QM
    !% region (a position in space or a list of atoms).  Also saves the
    !% old 'hybrid_mark' property in 'old_hybrid_mark' property.
    !% optionally correct selected region with heuristics (as coded in
    !% :func:`create_cluster_info_from_mark`())
    !
    subroutine create_pos_or_list_centred_hybrid_region(my_atoms,R_inner,R_outer,origin, atomlist,use_avgpos,add_only_heavy_atoms, &
  	     cluster_hopping_nneighb_only,cluster_heuristics_nneighb_only,min_images_only,use_create_cluster_info, create_cluster_info_args, list_changed, mark_postfix, &
               ignore_silica_residue, res_num_silica, error) ! lam81
  
      type(Atoms),        intent(inout) :: my_atoms
      real(dp),           intent(in)    :: R_inner
      real(dp),           intent(in)    :: R_outer
      real(dp), optional, intent(in)    :: origin(3)
      type(Table), optional, intent(in)    :: atomlist !the seed of the QM region
      logical,  optional, intent(in)   :: use_avgpos, add_only_heavy_atoms, cluster_hopping_nneighb_only, cluster_heuristics_nneighb_only, min_images_only, use_create_cluster_info
      character(len=*), optional, intent(in) :: create_cluster_info_args
      logical,  optional, intent(out)   :: list_changed
      character(len=*),  optional, intent(in)   :: mark_postfix
      integer, optional, intent(out) :: error
  
      logical ::  do_hopping_nneighb_only, do_heuristics_nneighb_only
      logical :: my_use_create_cluster_info
      type(Atoms) :: atoms_for_add_cut_hydrogens
      type(Table) :: core, old_core, old_all_but_term
      integer, pointer :: hybrid_p(:), hybrid_mark_p(:), old_hybrid_mark_p(:)
      character(len=STRING_LENGTH) :: my_create_cluster_info_args, my_mark_postfix
      integer, pointer :: hybrid_region_core_tmp_p(:)
      type(Table) :: padded_cluster_info
  
      logical,  optional, intent(in)   :: ignore_silica_residue ! lam81
      integer,  optional, intent(in)   :: res_num_silica        ! lam81
  
      INIT_ERROR(error)
  
      my_use_create_cluster_info = optional_default(.false., use_create_cluster_info)
  
      do_hopping_nneighb_only = optional_default(.false., cluster_hopping_nneighb_only)
      do_heuristics_nneighb_only = optional_default(.true., cluster_heuristics_nneighb_only)
  
      if (count((/present(origin),present(atomlist)/))/=1) then
        RAISE_ERROR('create_pos_or_list_centred_hybrid_mark: Exactly 1 of origin and atomlist must be present.', error)
      endif
  
      my_mark_postfix = optional_default("", trim(mark_postfix))
  
      call map_into_cell(my_atoms)
  
      ! save various subsets of atoms based on old hybrid marks
      call allocate(core,4,0,0,0,0)
      call allocate(old_core,4,0,0,0,0)
      call allocate(old_all_but_term,4,0,0,0,0)
      call get_hybrid_list(my_atoms,old_core,active_trans_only=.true., int_property='hybrid_mark'//trim(my_mark_postfix), error=error)
      PASS_ERROR_WITH_INFO("create_pos_or_list_centred_hybrid_region getting old core", error)
      call get_hybrid_list(my_atoms,core,active_trans_only=.true., int_property='hybrid_mark'//trim(my_mark_postfix), error=error)
      PASS_ERROR_WITH_INFO("create_pos_or_list_centred_hybrid_region getting core", error)
      call get_hybrid_list(my_atoms,old_all_but_term, all_but_term=.true., int_property='hybrid_mark'//trim(my_mark_postfix), error=error)
      PASS_ERROR_WITH_INFO("create_pos_or_list_centred_hybrid_region getting all but term", error)
  
     ! Build the new hysteretic QM core based on distances from point/list
     if (present(atomlist)) then
       call print("create_pos_or_list_centred_hybrid_region calling construct_hysteretic_region", verbosity=PRINT_NERD)
       call construct_hysteretic_region(region=core,at=my_atoms,core=atomlist,loop_atoms_no_connectivity=.false., &
         inner_radius=R_inner,outer_radius=R_outer, use_avgpos=use_avgpos, add_only_heavy_atoms=add_only_heavy_atoms, &
         cluster_hopping_nneighb_only=cluster_hopping_nneighb_only, min_images_only=min_images_only, &
         ignore_silica_residue=ignore_silica_residue, res_num_silica=res_num_silica, error=error) !NB , debugfile=mainlog) lam81
     else !present origin
       call print("create_pos_or_list_centred_hybrid_region calling construct_hysteretic_region", verbosity=PRINT_NERD)
       call construct_hysteretic_region(region=core,at=my_atoms,centre=origin,loop_atoms_no_connectivity=.true., &
         inner_radius=R_inner,outer_radius=R_outer, use_avgpos=use_avgpos, add_only_heavy_atoms=add_only_heavy_atoms, &
         cluster_hopping_nneighb_only=cluster_hopping_nneighb_only, min_images_only=min_images_only, &
         ignore_silica_residue=ignore_silica_residue, res_num_silica=res_num_silica, error=error) !NB , debugfile=mainlog) lam81
     endif
     PASS_ERROR_WITH_INFO("create_pos_or_list_centred_hybrid_region constructing hysteretic region", error)
  
     ! call create_cluster_info_from_mark, if requested, for various chemical intuition fixes to core
     if (my_use_create_cluster_info) then
        call add_property(my_atoms, "hybrid_region_core_tmp", HYBRID_NO_MARK, ptr=hybrid_region_core_tmp_p)
        hybrid_region_core_tmp_p(int_part(core,1)) = HYBRID_ACTIVE_MARK
        if (present(create_cluster_info_args) .and. (present(cluster_hopping_nneighb_only) .or. present(cluster_heuristics_nneighb_only))) then
  	 RAISE_ERROR("Got both create_cluster_info_args and (do_hopping_nneighb_only or do_heuristics_nneighb_only), but these conflict", error)
        endif
        my_create_cluster_info_args = optional_default("terminate=F cluster_hopping_nneighb_only="//do_hopping_nneighb_only// &
  	 " cluster_heuristics_nneighb_only="//do_heuristics_nneighb_only//" cluster_allow_modification", create_cluster_info_args)
        padded_cluster_info = create_cluster_info_from_mark(my_atoms, my_create_cluster_info_args, mark_name="hybrid_region_core_tmp")
        call finalise(core)
        call initialise(core, Nint=4, Nreal=0, Nstr=0, Nlogical=0, max_length=padded_cluster_info%N)
        call append(core, blank_rows=padded_cluster_info%N)
        core%int(1:4,1:padded_cluster_info%N) = padded_cluster_info%int(1:4,1:padded_cluster_info%N)
        call remove_property(my_atoms,"hybrid_region_core_tmp")
     endif
  
  !TO BE OPTIMIZED : add avgpos to add_cut_hydrogen
     ! add cut hydrogens, according to avgpos
      atoms_for_add_cut_hydrogens = my_atoms
      atoms_for_add_cut_hydrogens%oldpos = my_atoms%avgpos
      atoms_for_add_cut_hydrogens%avgpos = my_atoms%avgpos
      atoms_for_add_cut_hydrogens%pos = my_atoms%avgpos
  
      !% Use hysteretic connectivity calculator so we can work with relative cutoffs
      call calc_connect_hysteretic(atoms_for_add_cut_hydrogens, cutoff_factor=DEFAULT_NNEIGHTOL, &
           cutoff_break_factor=DEFAULT_NNEIGHTOL)
     
  
      call add_cut_hydrogens(atoms_for_add_cut_hydrogens,core)
      call finalise(atoms_for_add_cut_hydrogens)
  
      ! check changes in QM list and set the new QM list
      if (present(list_changed)) then
         list_changed = check_list_change(old_list=old_core,new_list=core)
         if (list_changed)  call print('QM list around the origin  has changed')
      endif
  
      ! sadd old hybrid_mark property
      call add_property(my_atoms,'old_hybrid_mark'//trim(my_mark_postfix),HYBRID_NO_MARK)
      if (.not. assign_pointer(my_atoms,'old_hybrid_mark'//trim(my_mark_postfix),old_hybrid_mark_p)) then
        RAISE_ERROR("create_pos_or_list_centred_hybrid_region couldn't get old_hybrid_mark"//trim(my_mark_postfix)//" property", error)
      endif
      ! update QM_flag of my_atoms
      if (.not. assign_pointer(my_atoms,'hybrid_mark'//trim(my_mark_postfix),hybrid_mark_p)) then
        RAISE_ERROR("create_pos_or_list_centred_hybrid_region couldn't get hybrid_mark"//trim(my_mark_postfix)//" property", error)
      endif
      ! save old_hybrid_mark
      old_hybrid_mark_p(1:my_atoms%N) = hybrid_mark_p(1:my_atoms%N)
      ! default to NO_MARK
      hybrid_mark_p(1:my_atoms%N) = HYBRID_NO_MARK
      ! anything which was marked before (except termination) is now BUFFER (so hysteretic buffer, later, has correct previous values)
      hybrid_mark_p(int_part(old_all_but_term,1)) = HYBRID_BUFFER_MARK
      ! anything returned in the core list is now ACTIVE
      hybrid_mark_p(int_part(core,1)) = HYBRID_ACTIVE_MARK
  
     ! update hybrid property of my_atoms
      if (.not. assign_pointer(my_atoms,'hybrid'//trim(my_mark_postfix),hybrid_p)) then
        RAISE_ERROR("create_pos_or_list_centred_hybrid_region couldn't get hybrid"//trim(my_mark_postfix)//" property", error)
      endif
      hybrid_p(1:my_atoms%N) = 0
      hybrid_p(int_part(core,1)) = 1
  
      call finalise(core)
      call finalise(old_core)
      call finalise(old_all_but_term)
  
    end subroutine create_pos_or_list_centred_hybrid_region
  
  !  !% Returns a $hybridlist$ table with the atom indices whose $cluster_mark$
  !  !% (or optionally any $int_property$) property takes no greater than $hybridflag$ positive value.
  !  !
  !  subroutine get_hybrid_list(my_atoms,hybridflag,hybridlist,int_property,get_up_to_mark_value)
  !
  !    type(Atoms), intent(in)  :: my_atoms
  !    integer,     intent(in)  :: hybridflag
  !    type(Table), intent(out) :: hybridlist
  !    character(len=*), optional, intent(in) :: int_property
  !    logical, intent(in), optional :: get_up_to_mark_value
  !
  !    integer, pointer :: mark_p(:)
  !    integer              :: i
  !    logical :: do_get_up_to_mark_value
  !
  !    do_get_up_to_mark_value = optional_default(.false., get_up_to_mark_value)
  !
  !    if (present(int_property)) then
  !       if (.not. assign_pointer(my_atoms, trim(int_property), mark_p)) then
  !	 RAISE_ERROR("get_hybrid_list_int couldn't get int_property='"//trim(int_property)//"'", error)
  !       endif
  !    else
  !       if (.not. assign_pointer(my_atoms, 'cluster_mark', mark_p)) then
  !	 RAISE_ERROR("get_hybrid_list_int couldn't get default int_property='cluster_mark'", error)
  !       endif
  !    endif
  !
  !    call initialise(hybridlist,4,0,0,0,0)      !1 int, 0 reals, 0 str, 0 log, num_hybrid_atoms entries
  !    if (do_get_up_to_mark_value) then
  !       do i=1,my_atoms%N
  !          ! if (my_atoms%data%int(hybrid_flag_index,i).gt.0.and. &
  !             ! my_atoms%data%int(hybrid_flag_index,i).le.hybridflag) &
  !          if (mark_p(i) > 0 .and.  mark_p(i) <= hybridflag) &
  !               call append(hybridlist,(/i,0,0,0/))
  !       enddo
  !    else
  !       do i=1,my_atoms%N
  !          ! if (my_atoms%data%int(hybrid_flag_index,i).eq.hybridflag) &
  !          if (mark_p(i) == hybridflag) call append(hybridlist,(/i,0,0,0/))
  !       enddo
  !    endif
  !
  !    if (hybridlist%N.eq.0) call print('Empty QM list with cluster_mark '//hybridflag,verbosity=PRINT_SILENT)
  !
  !  end subroutine get_hybrid_list
  
    !% Checks and reports the changes between two tables, $old_qmlist$ and $new_qmlist$.
    !
    function check_list_change(old_list,new_list) result(list_changed)
  
      type(Table), intent(in) :: old_list, new_list
      integer :: i
      logical :: list_changed
  
      list_changed = .false.
      if (old_list%N.ne.new_list%N) then
         call print ('list has changed: new number of atoms is: '//new_list%N//', was: '//old_list%N)
         list_changed = .true.
      else
         if (any(old_list%int(1,1:old_list%N).ne.new_list%int(1,1:new_list%N))) then
             do i=1,old_list%N
                if (.not.find_in_array(int_part(old_list,1),(new_list%int(1,i))).gt.0) then
                   call print('list has changed: atom '//new_list%int(1,i)//' has entered the region')
                   list_changed = .true.
                endif
                if (.not.find_in_array(int_part(new_list,1),(old_list%int(1,i))).gt.0) then
                   call print('list has changed: atom '//old_list%int(1,i)//' has left the region')
                   list_changed = .true.
                endif
             enddo
         endif
      endif
  
    end function check_list_change
  
    !% return list of atoms that have various subsets of hybrid marks set
    subroutine get_hybrid_list(at,hybrid_list,all_but_term,active_trans_only,int_property, error)
      type(Atoms), intent(in)  :: at !% object to scan for marked atoms
      type(Table), intent(out) :: hybrid_list !% on return, list of marked atoms
      logical, intent(in), optional :: all_but_term !% if present and true, select all marked atoms that aren't TERM
      logical, intent(in), optional :: active_trans_only !% if present and true, select all atoms marked ACTIVE or TRANS
      !% exactly one of all_but_term and active_trans_only must be present and true
      character(len=*), optional, intent(in) :: int_property !% if present, property to check, default cluster_mark
      integer, optional, intent(out) :: error
  
      integer :: i
      integer, pointer :: hybrid_mark(:)
      logical              :: my_all_but_term, my_active_trans_only
      character(STRING_LENGTH) :: my_int_property
  
      INIT_ERROR(error)
  
      if (.not. present(all_but_term) .and. .not. present(active_trans_only)) then
        RAISE_ERROR("get_hybrid_list called with neither all_but_term nor active_trans_only present", error)
      endif
  
      my_all_but_term = optional_default(.false., all_but_term)
      my_active_trans_only = optional_default(.false., active_trans_only)
  
      if ((my_all_but_term .and. my_active_trans_only) .or. (.not. my_all_but_term .and. .not. my_active_trans_only)) then
        RAISE_ERROR("get_hybrid_list needs exactly one of all_but_term=" // all_but_term // " and active_trans_only="//my_active_trans_only, error)
      endif
  
      my_int_property = ''
      if (present(int_property)) then
         my_int_property = trim(int_property)
      else
         my_int_property = "cluster_mark"
      endif
      if (.not.(assign_pointer(at, trim(my_int_property), hybrid_mark))) then
        RAISE_ERROR("get_hybrid_list couldn't find "//trim(my_int_property)//" field", error)
      endif
  
      call initialise(hybrid_list,4,0,0,0,0)      !1 int, 0 reals, 0 str, 0 log, num_qm_atoms entries
      do i=1, at%N
        if (my_all_but_term) then
  	if (hybrid_mark(i) /= HYBRID_NO_MARK .and. hybrid_mark(i) /= HYBRID_TERM_MARK) call append(hybrid_list,(/i,0,0,0/))
        else if (my_active_trans_only) then
  	if (hybrid_mark(i) == HYBRID_ACTIVE_MARK .or. hybrid_mark(i) == HYBRID_TRANS_MARK) call append(hybrid_list,(/i,0,0,0/))
        else
  	RAISE_ERROR("impossible! get_hybrid_list has no selection mode set", error)
        endif
      end do
  
      if (hybrid_list%N.eq.0) call print('get_hybrid_list returns empty hybrid list with field '//trim(my_int_property)// &
                                     ' all_but_term ' // my_all_but_term // ' active_trans_only ' // my_active_trans_only ,PRINT_ALWAYS)
    end subroutine get_hybrid_list
  
  
  end module clusters_module
  