From a48ed7c1567089adaceb8c85ea3a9fb2e4a5177f Mon Sep 17 00:00:00 2001
From: Axel Arnold <arnolda@icp.uni-stuttgart.de>
Date: Fri, 5 Nov 2004 23:16:45 +0000
Subject: [PATCH] Changes for the new exclusions code.

---
 RELEASE_NOTES      |  12 +-
 communication.c    | 211 ++++++-------
 communication.h    |  17 +-
 energy.h           |   2 +-
 forces.h           |   2 +-
 global.c           |  64 ++--
 global.h           |  64 ++--
 initialize.c       |  13 +-
 interaction_data.c |  84 -----
 interaction_data.h |  13 -
 particle_data.c    | 770 +++++++++++++++++++++++++++------------------
 particle_data.h    |  61 ++--
 pressure.c         |   2 +-
 pressure.h         |   2 +-
 rattle.c           |  64 ++--
 rattle.h           |   3 +
 verlet.c           |   4 +-
 17 files changed, 736 insertions(+), 652 deletions(-)

diff --git a/RELEASE_NOTES b/RELEASE_NOTES
index fe67141f5f0..dce7fc03364 100644
--- a/RELEASE_NOTES
+++ b/RELEASE_NOTES
@@ -8,7 +8,7 @@ What's new, what's been changed, what's been fixed
 Current working revision.
 
 
-- v1.8.1d -
+- v1.8.1e -
 -----------
 
 (e) November 5th, 2004.
@@ -18,6 +18,9 @@ Current working revision.
 (a) October 18th, 2004.
 
 New Features:
+(e) Completely rewritten exclusions code. Ready for testing.
+    Allows for general exclusions and is stored particle bound, so it
+    should be faster.
 (e) The Morse potential for nonbonded interactions has been added.
 (d) Added the famous soft-sphere potential.
 (c) Now the P3M limit for the mesh size is configurable and defaults to 128.
@@ -25,6 +28,9 @@ New Features:
 (a) Rattle is in, but probably not yet working.
 
 Changes:
+(e) Flags for the compilers changed. gcc is more aggressive now, icc 8.1 should
+    not issue the silly warnings anymore.
+(e) Arijit's exclusion code has more or less been removed.
 (c) Switched from POLYGAMMA_EPS to the commonly used ROUND_ERROR_PREC.
 (b) By request of too many users, grid.h now explicitely includes limits.h (which
     shouldn't be necessary on modern machines, but then you never know).
@@ -32,6 +38,10 @@ Changes:
     allocation problems (a 512 grid consumes roughly 3 GB memory!).
 
 Bugfixes:
+(e) fixed a mem leak in partCfg, bonds (updatePartCfg(WITH_BONDS)) where never freed.
+(e) fixed the particle print out, not handling particles by reference.
+(e) fixed try_delete_bond, deleting bonds probably did not work for a long time.
+(e) fixed COM with masses, was off by three
 (e) Added #ifdef statements for TABULATED and BUCKINGHAM potentials in the interprintall()   
 (b) Fixed the warning with the recent icc (Blade center).
 
diff --git a/communication.c b/communication.c
index 1b818471485..bd0f6e34e2c 100644
--- a/communication.c
+++ b/communication.c
@@ -65,9 +65,10 @@ typedef void (SlaveCallback)(int node, int param);
 /** Action number for \ref mpi_who_has. */
 #define REQ_WHO_HAS   2
 /** Action number for \ref mpi_bcast_event. */
-#define REQ_EVENT    3
+#define REQ_EVENT     3
 /** Action number for \ref mpi_place_particle. */
 #define REQ_PLACE     4
+
 /** Action number for \ref mpi_send_v. */
 #define REQ_SET_V     5
 /** Action number for \ref mpi_send_f. */
@@ -78,6 +79,7 @@ typedef void (SlaveCallback)(int node, int param);
 #define REQ_SET_TYPE  8
 /** Action number for \ref mpi_send_bond. */
 #define REQ_SET_BOND  9
+
 /** Action number for \ref mpi_recv_part. */
 #define REQ_GET_PART  10
 /** Action number for \ref mpi_integrate. */
@@ -88,6 +90,7 @@ typedef void (SlaveCallback)(int node, int param);
 #define REQ_BCAST_IA_SIZE  13
 /** Action number for \ref mpi_gather_stats. */
 #define REQ_GATHER    14
+
 /** Action number for \ref mpi_set_time_step. */
 #define REQ_SET_TIME_STEP  15
 /** Action number for \ref mpi_get_particles. */
@@ -98,6 +101,7 @@ typedef void (SlaveCallback)(int node, int param);
 #define REQ_SET_EXT 18
 /** Action number for \ref mpi_place_particle. */
 #define REQ_PLACE_NEW   19
+
 /** Action number for \ref mpi_remove_particle */
 #define REQ_REM_PART   20
 /** Action number for \ref mpi_bcast_constraint */
@@ -108,6 +112,7 @@ typedef void (SlaveCallback)(int node, int param);
 #define REQ_RANDOM_STAT 23
 /** Action number for \ref mpi_lj_cap_forces. */
 #define REQ_BCAST_LFC 24
+
 /** Action number for \ref mpi_tab_cap_forces. */
 #define REQ_BCAST_TFC 25
 /** Action number for \ref mpi_random_seed */
@@ -118,16 +123,18 @@ typedef void (SlaveCallback)(int node, int param);
 #define REQ_GET_CONSFOR 28
 /** Action number for \ref mpi_rescale_particles */
 #define REQ_RESCALE_PART 29
+
 /** Action number for \ref mpi_bcast_cell_structure */
-#define REQ_BCAST_CS 30
+#define REQ_BCAST_CS  30
 /** Action number for \ref mpi_send_quat. */
-#define REQ_SET_QUAT 31
+#define REQ_SET_QUAT  31
 /** Action number for \ref mpi_send_omega. */
 #define REQ_SET_OMEGA 32
 /** Action number for \ref mpi_send_torque. */
 #define REQ_SET_TORQUE 33
 /** Action number for \ref mpi_send_mol_id. */
 #define REQ_SET_MOLID 34
+
 /** Action number for \ref mpi_bcast_nptiso_geom */
 #define REQ_BCAST_NPTISO_GEOM 35
 /** Action number for \ref mpi_update_mol_ids  */
@@ -135,15 +142,16 @@ typedef void (SlaveCallback)(int node, int param);
 /** Action number for \ref mpi_sync_topo_part_info  */
 #define REQ_SYNC_TOPO 37
 /** Action number for \ref mpi_send_mass. */
-#define REQ_SET_MASS     38
-/** Action number for \ref mpi_bcast_bond_partners.*/
-#define REQ_BCAST_BOND_PARTNERS  39
+#define REQ_SET_MASS  38
 /** Action number for \ref mpi_buck_cap_forces. */
-#define REQ_BCAST_BFC 40
-/**Action number for \ref mpi_morse_cap_forces. */
-#define REQ_BCAST_MFC 41
+#define REQ_BCAST_BFC 39
+
 /** Action number for \ref mpi_gather_runtime_errors  */
-#define REQ_GET_ERRS 42
+#define REQ_GET_ERRS  40
+/** Action number for \ref mpi_send_exclusion. */
+#define REQ_SET_EXCL  41
+/**Action number for \ref mpi_morse_cap_forces. */
+#define REQ_BCAST_MFC 42
 /** Total number of action numbers. */
 #define REQ_MAXIMUM   43
 
@@ -179,7 +187,6 @@ void mpi_bcast_constraint_slave(int node, int parm);
 void mpi_random_seed_slave(int node, int parm);
 void mpi_random_stat_slave(int node, int parm);
 void mpi_lj_cap_forces_slave(int node, int parm);
-void mpi_morse_cap_forces_slave(int node, int parm);
 void mpi_tab_cap_forces_slave(int node, int parm);
 void mpi_bit_random_seed_slave(int node, int parm);
 void mpi_bit_random_stat_slave(int node, int parm);
@@ -194,10 +201,10 @@ void mpi_bcast_nptiso_geom_slave(int node, int parm);
 void mpi_update_mol_ids_slave(int node, int parm);
 void mpi_sync_topo_part_info_slave(int node, int parm);
 void mpi_send_mass_slave(int node, int parm);
-void mpi_bcast_bond_partners_slave(int node, int parm);
 void mpi_buck_cap_forces_slave(int node, int parm);
 void mpi_gather_runtime_errors_slave(int node, int parm);
-
+void mpi_send_exclusion_slave(int node, int parm);
+void mpi_morse_cap_forces_slave(int node, int parm);
 /*@}*/
 
 /** A list of which function has to be called for
@@ -242,10 +249,10 @@ SlaveCallback *callbacks[] = {
   mpi_update_mol_ids_slave,         /* 36: REQ_UPDATE_MOL_IDS */
   mpi_sync_topo_part_info_slave,    /* 37: REQ_SYNC_TOPO */
   mpi_send_mass_slave,              /* 38: REQ_SET_MASS */
-  mpi_bcast_bond_partners_slave,    /* 39: REQ_BCAST_BOND_PARTNERS*/
-  mpi_buck_cap_forces_slave,        /* 40: REQ_BCAST_BFC */
-  mpi_morse_cap_forces_slave,       /* 41: REQ_BCAST_MFC */
-  mpi_gather_runtime_errors_slave,  /* 42: REQ_GET_ERRS */
+  mpi_buck_cap_forces_slave,        /* 39: REQ_BCAST_LFC */
+  mpi_gather_runtime_errors_slave,  /* 40: REQ_GET_ERRS */
+  mpi_send_exclusion_slave,         /* 41: REQ_SET_EXCL */
+  mpi_morse_cap_forces_slave,       /* 42: REQ_BCAST_MFC */
 };
 
 /** Names to be printed when communication debugging is on. */
@@ -255,44 +262,52 @@ char *names[] = {
   "WHO_HAS"   ,     /*  2 */
   "EVENT"     ,     /*  3 */
   "SET_POS"   ,     /*  4 */
+
   "SET_V"     ,     /*  5 */
   "SET_F"     ,     /*  6 */
   "SET_Q"     ,     /*  7 */
   "SET_TYPE"  ,     /*  8 */
   "SET_BOND"  ,     /*  9 */
+
   "GET_PART"  ,     /* 10 */
   "INTEGRATE" ,     /* 11 */
   "BCAST_IA"  ,     /* 12 */
   "BCAST_IAS" ,     /* 13 */
   "GATHER"    ,     /* 14 */
+
   "TIME_STEP" ,     /* 15 */
   "GET_PARTS" ,     /* 16 */
   "BCAST_CIA" ,     /* 17 */
   "SEND_EXT"  ,     /* 18 */
   "PLACE_NEW" ,     /* 19 */
+
   "REM_PART"  ,     /* 20 */
   "BCAST_CON" ,     /* 21 */
   "RAND_SEED" ,     /* 22 */
   "RAND_STAT" ,     /* 23 */
   "BCAST_LFC" ,     /* 24 */
+
   "BCAST_TFC" ,     /* 25 */
   "BIT_RAND_SEED",  /* 26 */
   "BIT_RAND_STAT",  /* 27 */
   "GET_CONSTR",     /* 28 */
   "RESCALE_PART",   /* 29 */
+
   "BCAST_CS",       /* 30 */
   "SET_QUAT"  ,     /* 31 */
   "SET_OMEGA",      /* 32 */
   "SET_TORQUE",     /* 33 */
   "SET_MOLID",      /* 34 */
+
   "BCAST_NPT_GEOM", /* 35 */
   "UPDATE_MOL_IDS", /* 36 */
   "SYNC_TOPO",      /* 37 */
-  "SET_MASS"        /* 38 */
-  "BCAST_BOND_PARTNERS" /* 39*/
-  "BCAST_BFC" ,     /* 40 */
-  "BCAST_MFC" ,     /* 41 */
-  "GET_ERRS",       /* 42 */
+  "SET_MASS",       /* 38 */
+  "BCAST_BFC" ,     /* 39 */
+
+  "GET_ERRS",       /* 40 */
+  "SET_EXCL",       /* 41 */
+  "BCAST_MFC" ,     /* 42 */
 };
 
 /** the requests are compiled here. So after a crash you get the last issued request */
@@ -306,19 +321,7 @@ static int request[3];
 void mpi_core(MPI_Comm *comm, int *errcode,...) {
   fprintf(stderr, "%d: Aborting due to MPI error %d, forcing core dump\n", this_node, *errcode);
   fflush(stderr);
-  /* no longer necessary with core.* files.
-     Actually, it is nice to have all nodes...
-     #ifdef MPI_ERR_LOCALDEAD
-     #ifdef MPI_ERR_REMOTEDEAD
-     if (*errcode != MPI_ERR_LOCALDEAD &&
-     *errcode != MPI_ERR_REMOTEDEAD)
-     core();
-     #else
-  */
   core();
-  /* #endif
-     #endif
-  */
 }
 #endif
 
@@ -887,39 +890,48 @@ void mpi_send_bond_slave(int pnode, int part)
 void mpi_recv_part(int pnode, int part, Particle *pdata)
 {
   IntList *bl = &(pdata->bl);
+#ifdef EXCLUSIONS
+  IntList *el = &(pdata->el);
+#endif
   MPI_Status status;
 
-  if (pnode == this_node) {
-    Particle *p = local_particles[part];
-
-    memcpy(pdata, p, sizeof(Particle));
-
-    bl->max = bl->n = p->bl.n;
-
-    if (bl->n > 0) {
-      alloc_intlist(bl, bl->n);
-      memcpy(bl->e, p->bl.e, sizeof(int)*bl->n);
-    }
-    else
-      bl->e = NULL;
-  }
+  /* fetch fixed data */
+  if (pnode == this_node)
+    memcpy(pdata, local_particles[part], sizeof(Particle));
   else {
     mpi_issue(REQ_GET_PART, pnode, part);
-
-    /* add new properties here */
-
-    pdata->p.identity = part;
     MPI_Recv(pdata, sizeof(Particle), MPI_BYTE, pnode,
 	     REQ_GET_PART, MPI_COMM_WORLD, &status);
-    bl->max = bl->n;
-    if (bl->n > 0) {
-      alloc_intlist(bl, bl->n);
+  }
+
+  /* copy dynamic data */
+  /* bonds */
+  bl->max = bl->n;
+  if (bl->n > 0) {
+    alloc_intlist(bl, bl->n);
+    if (pnode == this_node)
+      memcpy(bl->e, local_particles[part]->bl.e, sizeof(int)*bl->n);
+    else
       MPI_Recv(bl->e, bl->n, MPI_INT, pnode,
-	       REQ_GET_PART, MPI_COMM_WORLD, &status);
-    }
+               REQ_GET_PART, MPI_COMM_WORLD, &status);
+  }
+  else
+    bl->e = NULL;
+
+#ifdef EXCLUSIONS
+  /* exclusions */
+  el->max = el->n;
+  if (el->n > 0) {
+    alloc_intlist(el, el->n);
+    if (pnode == this_node)
+      memcpy(el->e, local_particles[part]->el.e, sizeof(int)*el->n);
     else
-      pdata->bl.e = NULL;
+      MPI_Recv(el->e, el->n, MPI_INT, pnode,
+               REQ_GET_PART, MPI_COMM_WORLD, &status);
   }
+  else
+    el->e = NULL;
+#endif
 }
 
 void mpi_recv_part_slave(int pnode, int part)
@@ -935,6 +947,11 @@ void mpi_recv_part_slave(int pnode, int part)
   if (p->bl.n > 0)
     MPI_Send(p->bl.e, p->bl.n, MPI_INT, 0, REQ_GET_PART,
 	     MPI_COMM_WORLD);
+#ifdef EXCLUSIONS
+  if (p->el.n > 0)
+    MPI_Send(p->el.e, p->el.n, MPI_INT, 0, REQ_GET_PART,
+	     MPI_COMM_WORLD);
+#endif
 }
 
 /****************** REQ_REM_PART ************/
@@ -1719,63 +1736,6 @@ void mpi_bcast_nptiso_geom_slave(int node,int parm)
 
 }
 
-/*************** REQ_BCAST_BOND_PARTNERS *****************/
-
-void mpi_bcast_bond_partners(int width, int parm)
-{
-  mpi_issue(REQ_BCAST_BOND_PARTNERS, width, parm);
-  mpi_bcast_bond_partners_slave(width, parm);
-
-}
-
-void mpi_bcast_bond_partners_slave(int width, int parm)
-{
-  int dim = n_total_particles*width;
-  //printf("%d: width = %d parm = %d dimensions of partBondPartners %d\n",this_node, width, parm, dim);
-  if(parm==REALLOC_BOND_PARTNERS)
-  {
-     //free(partBondPartners);
-     //int *partBondPartners;
-     partBondPartners = (int *) realloc(partBondPartners, dim*sizeof(int));
-     //if(this_node==0)
-     //   partBondPartners = (int *) realloc(partBondPartners, n_total_particles*width*sizeof(int));
-     //MPI_Bcast(&partBondPartners, n_total_particles*width, MPI_BYTE, 0, MPI_COMM_WORLD);
-     MPI_Barrier(MPI_COMM_WORLD);
-     //printf("%d: after MPI_Bcast\n",this_node);
-  }
-  else if(parm==CURR_BOND_PARTNERS)
-  {
-    //printf("%d: inside CURR_BOND_PARTNERS, dimensions: %d\n",this_node, dim);
-    //partBondPartners = (int *) realloc(partBondPartners, n_total_particles*0*sizeof(int));
-    //int i;
-    //for(i=0;i<dim;i++)
-    //  MPI_Bcast(&partBondPartners[i], 1, MPI_INT, 0, MPI_COMM_WORLD);
-    //MPI_Bcast((int *)partBondPartners, dim*sizeof(int), MPI_INT, 0, MPI_COMM_WORLD);
-    MPI_Bcast((int *)partBondPartners, dim, MPI_INT, 0, MPI_COMM_WORLD);
-    MPI_Barrier(MPI_COMM_WORLD);
-    //printf("%d: got partBondPartners\n",this_node);
-    //printf("%d: partBondPartners[0]=%d\n",this_node,*(partBondPartners+0));
-    /*
-    int j;
-    for(i=0;i<n_total_particles;i=i+1)
-    {
-    printf("%d:\t",i);
-    for (j=0;j<width;j++)
-     printf("%d\t",partBondPartners[i*width+j]);
-    printf("\n");
-    }
-    */
-  }
-  else
-  {
-    fprintf(stderr,"ERROR: Wrong Parameter passed to 'mpi_bcast_bond_partners'\n");
-    errexit();
-  }
-
-  //printf("Comm Node: %d\n",this_node);
-
-}
-
 /***************REQ_UPDATE_MOL_IDS *********************/
 
 void mpi_update_mol_ids()
@@ -1913,6 +1873,27 @@ void mpi_gather_runtime_errors_slave(int node, int parm)
   }
 }
 
+/********************* REQ_SET_EXCL ********/
+void mpi_send_exclusion(int part1, int part2, int delete)
+{
+#ifdef EXCLUSIONS
+  mpi_issue(REQ_SET_EXCL, part1, part2);
+
+  MPI_Bcast(&delete, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  local_change_exclusion(part1, part2, delete);
+  on_particle_change();
+#endif
+}
+
+void mpi_send_exclusion_slave(int part1, int part2)
+{
+#ifdef EXCLUSIONS
+  int delete;
+  MPI_Bcast(&delete, 1, MPI_INT, 0, MPI_COMM_WORLD);  
+  local_change_exclusion(part1, part2, delete);
+  on_particle_change();
+#endif
+}
 
 /*********************** MAIN LOOP for slaves ****************/
 
diff --git a/communication.h b/communication.h
index c93f95886f3..82b2fe3e77f 100644
--- a/communication.h
+++ b/communication.h
@@ -194,6 +194,14 @@ void mpi_send_mol_id(int node, int part, int mid);
 */
 int mpi_send_bond(int pnode, int part, int *bond, int delete);
 
+/** Issue REQ_SET_EXCLUSION: send exclusions.
+    Also calls \ref on_particle_change.
+    \param part     identity of first particle of the exclusion.
+    \param part2    identity of secnd particle of the exclusion.
+    \param delete   if true, do not add the exclusion, rather delete it if found
+*/
+void mpi_send_exclusion(int part, int part2, int delete);
+
 /** Issue REQ_REM_PART: remove a particle.
     Also calls \ref on_particle_change.
     \param id   the particle to remove.
@@ -201,7 +209,9 @@ int mpi_send_bond(int pnode, int part, int *bond, int delete);
 */
 void mpi_remove_particle(int node, int id);
 
-/** Issue REQ_GET_PART: recv particle data.
+/** Issue REQ_GET_PART: recv particle data. The data has to be freed later
+    using \ref free_particle, otherwise the dynamically allocated parts, bonds
+    and exclusions are left over.
     \param part the particle.
     \param node the node it is attached to.
     \param part_data where to store the received data.
@@ -258,6 +268,7 @@ void mpi_gather_stats(int job, void *result);
     pointed to by il. The particles bonding information references this array,
     which is the only data you have to free later (besides the result array
     you allocated). YOU MUST NOT CALL \ref free_particle on any of these particles!
+    
   \param result where to store the gathered particles
   \param il if non-NULL, the integerlist where to store the bonding info
 */
@@ -323,10 +334,6 @@ void mpi_update_mol_ids(void);
 /** Issue REQ_SYNC_TOPO: Update the molecules ids to that they correspond to the topology */
 int mpi_sync_topo_part_info(void);
 
-/** Issue REQ_BCAST_BOND_PARTNERS: Broadcast the sorted array containing the mapping of
-    particle and its bond partners */
-void mpi_bcast_bond_partners(int width, int parm);
-
 /** Issue REQ_GET_ERRS: gather all error messages from all nodes and set the interpreter result
     to these error messages. This called only on the master node.
     The errors are append to the result, if ret_state == TCL_ERROR, otherwise the result is overwritten.
diff --git a/energy.h b/energy.h
index aab7073cc01..82181313224 100644
--- a/energy.h
+++ b/energy.h
@@ -137,7 +137,7 @@ MDINLINE void add_bonded_energy(Particle *p1)
   int i, type_num, type, n_partners, bond_broken;
   double ret, dx[3] = {0, 0, 0};
 
-  i=0;
+  i = 0;
   while(i<p1->bl.n) {
     type_num = p1->bl.e[i++];
     iaparams = &bonded_ia_params[type_num];
diff --git a/forces.h b/forces.h
index 5c029d0c2bf..78c1ad4286c 100644
--- a/forces.h
+++ b/forces.h
@@ -214,7 +214,7 @@ MDINLINE void add_bonded_force(Particle *p1)
   Bonded_ia_parameters *iaparams;
   int i, j, type_num, type, n_partners, bond_broken;
 
-  i=0;
+  i = 0;
   while(i<p1->bl.n) {
     type_num = p1->bl.e[i++];
     iaparams = &bonded_ia_params[type_num];
diff --git a/global.c b/global.c
index 989fb6aeb77..52c4d3cdc9a 100644
--- a/global.c
+++ b/global.c
@@ -31,6 +31,7 @@
 #include "domain_decomposition.h"
 #include "layered.h"
 #include "pressure.h"
+#include "rattle.h"
 
 /**********************************************
  * description of variables
@@ -51,38 +52,37 @@ const Datafield fields[] = {
   {&dpd_gamma,       TYPE_DOUBLE, 1, "dpd_gamma",     ro_callback,    5 },         /* 3  from thermostat.c */
   {&dpd_r_cut,       TYPE_DOUBLE, 1, "dpd_r_cut",     ro_callback,    5 },         /* 4  from thermostat.c */
   {&langevin_gamma,  TYPE_DOUBLE, 1, "gamma",         langevin_gamma_callback, 1 },/* 5  from thermostat.c */
-  {&ia_excl,            TYPE_INT, 1, "ia_excl",       ro_callback,    5 },         /* 6  from interaction_data.c */
-  {&integ_switch,       TYPE_INT, 1, "integ_switch",  ro_callback,    1 },         /* 7  from integrate.c */
-  {local_box_l,      TYPE_DOUBLE, 3, "local_box_l",   ro_callback,    2 },         /* 8  from global.c */
-  {&max_cut,         TYPE_DOUBLE, 1, "max_cut",       ro_callback,    5 },         /* 9  from interaction_data.c */
-  {&max_num_cells,      TYPE_INT, 1, "max_num_cells", max_num_cells_callback, 5 }, /* 10 from cells.c */
-  {&max_seen_particle,  TYPE_INT, 1, "max_part",      ro_callback,    5 },         /* 11 from particle_data.c */
-  {&max_range,       TYPE_DOUBLE, 1, "max_range",     ro_callback,    5 },         /* 12 from integrate.c */
-  {&max_skin,        TYPE_DOUBLE, 1, "max_skin",      ro_callback,    5 },         /* 13 from integrate.c */
-  {&min_num_cells,      TYPE_INT, 1, "min_num_cells", min_num_cells_callback, 5 }, /* 14  from cells.c */
-  {&n_layers,           TYPE_INT, 1, "n_layers",      ro_callback,    3 },         /* 15 from layered.c */
-  {&n_nodes,            TYPE_INT, 1, "n_nodes",       ro_callback,    3 },         /* 16 from communication.c */
-  {&n_total_particles,  TYPE_INT, 1, "n_part",        ro_callback,    6 },         /* 17 from particle.c */
-  {&n_particle_types,   TYPE_INT, 1, "n_part_types",  ro_callback,    8 },         /* 18 from interaction_data.c */
-  {&n_rigidbonds,       TYPE_INT, 1, "n_rigidbonds",  ro_callback,    5 },         /* 19 from interaction_data.c */
-  {node_grid,           TYPE_INT, 3, "node_grid",     node_grid_callback, 2 },     /* 20 from grid.c */
-  {&nptiso_gamma0,   TYPE_DOUBLE, 1, "nptiso_gamma0", ro_callback,    13 },        /* 21 from thermostat.c */
-  {&nptiso_gammav,   TYPE_DOUBLE, 1, "nptiso_gammav", ro_callback,    13 },        /* 22 from thermostat.c */
-  {&nptiso.p_ext,    TYPE_DOUBLE, 1, "npt_p_ext",     ro_callback,     7 },        /* 23 from pressure.c */
-  {&nptiso.p_inst,   TYPE_DOUBLE, 1, "npt_p_inst",    ro_callback,    10 },        /* 24 from pressure.c */
-  {&nptiso.p_inst_av,TYPE_DOUBLE, 1, "npt_p_inst_av", ro_callback,    10 },        /* 25 from pressure.c */
-  {&nptiso.p_diff,   TYPE_DOUBLE, 1, "npt_p_diff",    p_diff_callback, 7 },        /* 26 from pressure.c */
-  {&nptiso.piston,   TYPE_DOUBLE, 1, "npt_piston",    piston_callback, 6 },        /* 27 from pressure.c */
-  {&periodic,          TYPE_BOOL, 3, "periodicity",   per_callback,    1 },        /* 28 from grid.c */
-  {&skin,            TYPE_DOUBLE, 1, "skin",          skin_callback,   2 },        /* 29 from integrate.c */
-  {&temperature,     TYPE_DOUBLE, 1, "temperature",   temp_callback,   2 },        /* 30 from thermostat.c */
-  {&thermo_switch,      TYPE_INT, 1, "thermo_switch", ro_callback,     2 },        /* 31 from thermostat.c */
-  {&sim_time,        TYPE_DOUBLE, 1, "time",          time_callback,   4 },        /* 32 from integrate.c */
-  {&time_step,       TYPE_DOUBLE, 1, "time_step",     time_step_callback, 5 },     /* 33 from integrate.c */
-  {&timing_samples,     TYPE_INT, 1, "timings",       timings_callback, 4 },       /* 34 from tuning.c */
-  {&transfer_rate,      TYPE_INT, 1, "transfer_rate", ro_callback,     2 },        /* 35 from imd.c */
-  {&rebuild_verletlist,TYPE_BOOL, 1, "verlet_flag",   ro_callback,     8 },        /* 36 from verlet.c */
-  {&verlet_reuse,    TYPE_DOUBLE, 1, "verlet_reuse",  ro_callback,     8 },        /* 37 from integrate.c */
+  {&integ_switch,       TYPE_INT, 1, "integ_switch",  ro_callback,    1 },         /* 6  from integrate.c */
+  {local_box_l,      TYPE_DOUBLE, 3, "local_box_l",   ro_callback,    2 },         /* 7  from global.c */
+  {&max_cut,         TYPE_DOUBLE, 1, "max_cut",       ro_callback,    5 },         /* 8  from interaction_data.c */
+  {&max_num_cells,      TYPE_INT, 1, "max_num_cells", max_num_cells_callback, 5 }, /* 9 from cells.c */
+  {&max_seen_particle,  TYPE_INT, 1, "max_part",      ro_callback,    5 },         /* 10 from particle_data.c */
+  {&max_range,       TYPE_DOUBLE, 1, "max_range",     ro_callback,    5 },         /* 11 from integrate.c */
+  {&max_skin,        TYPE_DOUBLE, 1, "max_skin",      ro_callback,    5 },         /* 12 from integrate.c */
+  {&min_num_cells,      TYPE_INT, 1, "min_num_cells", min_num_cells_callback, 5 }, /* 13  from cells.c */
+  {&n_layers,           TYPE_INT, 1, "n_layers",      ro_callback,    3 },         /* 14 from layered.c */
+  {&n_nodes,            TYPE_INT, 1, "n_nodes",       ro_callback,    3 },         /* 15 from communication.c */
+  {&n_total_particles,  TYPE_INT, 1, "n_part",        ro_callback,    6 },         /* 16 from particle.c */
+  {&n_particle_types,   TYPE_INT, 1, "n_part_types",  ro_callback,    8 },         /* 17 from interaction_data.c */
+  {&n_rigidbonds,       TYPE_INT, 1, "n_rigidbonds",  ro_callback,    5 },         /* 18 from rattle.c */
+  {node_grid,           TYPE_INT, 3, "node_grid",     node_grid_callback, 2 },     /* 19 from grid.c */
+  {&nptiso_gamma0,   TYPE_DOUBLE, 1, "nptiso_gamma0", ro_callback,    13 },        /* 20 from thermostat.c */
+  {&nptiso_gammav,   TYPE_DOUBLE, 1, "nptiso_gammav", ro_callback,    13 },        /* 21 from thermostat.c */
+  {&nptiso.p_ext,    TYPE_DOUBLE, 1, "npt_p_ext",     ro_callback,     7 },        /* 22 from pressure.c */
+  {&nptiso.p_inst,   TYPE_DOUBLE, 1, "npt_p_inst",    ro_callback,    10 },        /* 23 from pressure.c */
+  {&nptiso.p_inst_av,TYPE_DOUBLE, 1, "npt_p_inst_av", ro_callback,    10 },        /* 24 from pressure.c */
+  {&nptiso.p_diff,   TYPE_DOUBLE, 1, "npt_p_diff",    p_diff_callback, 7 },        /* 25 from pressure.c */
+  {&nptiso.piston,   TYPE_DOUBLE, 1, "npt_piston",    piston_callback, 6 },        /* 26 from pressure.c */
+  {&periodic,          TYPE_BOOL, 3, "periodicity",   per_callback,    1 },        /* 27 from grid.c */
+  {&skin,            TYPE_DOUBLE, 1, "skin",          skin_callback,   2 },        /* 28 from integrate.c */
+  {&temperature,     TYPE_DOUBLE, 1, "temperature",   temp_callback,   2 },        /* 29 from thermostat.c */
+  {&thermo_switch,      TYPE_INT, 1, "thermo_switch", ro_callback,     2 },        /* 30 from thermostat.c */
+  {&sim_time,        TYPE_DOUBLE, 1, "time",          time_callback,   4 },        /* 31 from integrate.c */
+  {&time_step,       TYPE_DOUBLE, 1, "time_step",     time_step_callback, 5 },     /* 32 from integrate.c */
+  {&timing_samples,     TYPE_INT, 1, "timings",       timings_callback, 4 },       /* 33 from tuning.c */
+  {&transfer_rate,      TYPE_INT, 1, "transfer_rate", ro_callback,     2 },        /* 34 from imd.c */
+  {&rebuild_verletlist,TYPE_BOOL, 1, "verlet_flag",   ro_callback,     8 },        /* 35 from verlet.c */
+  {&verlet_reuse,    TYPE_DOUBLE, 1, "verlet_reuse",  ro_callback,     8 },        /* 36 from integrate.c */
   { NULL, 0, 0, NULL, NULL, 0 }
 };
 
diff --git a/global.h b/global.h
index 478f42c91a4..cef3544ae65 100644
--- a/global.h
+++ b/global.h
@@ -89,70 +89,68 @@ extern const Datafield fields[];
 #define FIELD_DPD_RCUT            4
 /** index of \ref langevin_gamma in  \ref #fields */
 #define FIELD_LANGEVIN_GAMMA      5
-/** index of \ref ia_excl in \ref #fields */
-#define FIELD_IA_EXCL             6
 /** index of \ref integ_switch in \ref #fields */
-#define FIELD_INTEG_SWITCH        7
+#define FIELD_INTEG_SWITCH        6
 /** index of \ref local_box_l in \ref #fields */
-#define FIELD_LBOXL               8
+#define FIELD_LBOXL               7
 /** index of \ref max_cut in \ref #fields */
-#define FIELD_MCUT                9
+#define FIELD_MCUT                8
 /** index of \ref max_num_cells  in \ref #fields */
-#define FIELD_MAXNUMCELLS         10
+#define FIELD_MAXNUMCELLS         9
 /** index of \ref max_seen_particle in \ref #fields */
-#define FIELD_MAXPART             11
+#define FIELD_MAXPART             10
 /** index of \ref max_range in \ref #fields */
-#define FIELD_MAXRANGE            12
+#define FIELD_MAXRANGE            11
 /** index of \ref max_skin in  \ref #fields */
-#define FIELD_MAXSKIN             13
+#define FIELD_MAXSKIN             12
 /** index of \ref min_num_cells  in \ref #fields */
-#define FIELD_MINNUMCELLS         14
+#define FIELD_MINNUMCELLS         13
 /** index of \ref n_layers in  \ref #fields */
-#define FIELD_NLAYERS             15
+#define FIELD_NLAYERS             14
 /** index of \ref n_nodes in \ref #fields */
-#define FIELD_NNODES              16
+#define FIELD_NNODES              15
 /** index of \ref n_total_particles in  \ref #fields */
-#define FIELD_NPART               17
+#define FIELD_NPART               16
 /** index of \ref n_particle_types in \ref #fields */
-#define FIELD_NPARTTYPE           18
+#define FIELD_NPARTTYPE           17
 /** index of \ref n_rigidbonds in \ref #fields */
-#define FIELD_RIGIDBONDS          19
+#define FIELD_RIGIDBONDS          18
 /** index of \ref node_grid in \ref #fields */
-#define FIELD_NODEGRID            20
+#define FIELD_NODEGRID            19
 /** index of \ref nptiso_gamma0 in \ref #fields */
-#define FIELD_NPTISO_G0           21
+#define FIELD_NPTISO_G0           20
 /** index of \ref nptiso_gammav in \ref #fields */
-#define FIELD_NPTISO_GV           22
+#define FIELD_NPTISO_GV           21
 /** index of \ref nptiso_struct::p_ext in \ref #fields */
-#define FIELD_NPTISO_PEXT         23      
+#define FIELD_NPTISO_PEXT         22      
 /** index of \ref nptiso_struct::p_inst in \ref #fields */
-#define FIELD_NPTISO_PINST        24
+#define FIELD_NPTISO_PINST        23
 /** index of \ref nptiso_struct::p_inst_av in \ref #fields */
-#define FIELD_NPTISO_PINSTAV      25
+#define FIELD_NPTISO_PINSTAV      24
 /** index of \ref nptiso_struct::p_diff in \ref #fields */
-#define FIELD_NPTISO_PDIFF        26
+#define FIELD_NPTISO_PDIFF        25
 /** index of \ref nptiso_struct::piston in \ref #fields */
-#define FIELD_NPTISO_PISTON       27
+#define FIELD_NPTISO_PISTON       26
 /** index of \ref #periodic in \ref #fields */
-#define FIELD_PERIODIC            28
+#define FIELD_PERIODIC            27
 /** index of \ref #skin in \ref #fields */
-#define FIELD_SKIN                29
+#define FIELD_SKIN                28
 /** index of \ref #temperature in \ref #fields */
-#define FIELD_TEMPERATURE         30
+#define FIELD_TEMPERATURE         29
 /** index of \ref thermo_switch in \ref #fields */
-#define FIELD_THERMO_SWITCH       31
+#define FIELD_THERMO_SWITCH       30
 /** index of \ref sim_time in  \ref #fields */
-#define FIELD_SIMTIME             32
+#define FIELD_SIMTIME             31
 /** index of \ref time_step in \ref #fields */
-#define FIELD_TIMESTEP            33
+#define FIELD_TIMESTEP            32
 /** index of \ref timing_samples in  \ref #fields */
-#define FIELD_TIMINGSAMP          34
+#define FIELD_TIMINGSAMP          33
 /** index of \ref transfer_rate  in \ref #fields */
-#define FIELD_TRANSFERRATE        35
+#define FIELD_TRANSFERRATE        34
 /** index of \ref rebuild_verletlist in \ref #fields */
-#define FIELD_VERLETFLAG          36
+#define FIELD_VERLETFLAG          35
 /** index of \ref verlet_reuse in  \ref #fields */
-#define FIELD_VERLETREUSE         37
+#define FIELD_VERLETREUSE         36
 
 /*@}*/
 
diff --git a/initialize.c b/initialize.c
index 4ad9e661501..961d95988fb 100644
--- a/initialize.c
+++ b/initialize.c
@@ -48,6 +48,7 @@
 #include "nemd.h"
 #include "domain_decomposition.h"
 #include "errorhandling.h"
+#include "rattle.h"
 
 /** whether before integration the thermostat has to be reinitialized */
 static int reinit_thermo = 1;
@@ -144,7 +145,7 @@ void on_integration_start()
 
   /* Update particle and observable information for routines in statistics.c */
   invalidate_obs();
-  free(partCfg); partCfg=NULL;
+  freePartCfg();
 
   /* Prepare particle structure: Communication step: number of ghosts and ghost information */
   if(resort_particles) {
@@ -163,7 +164,7 @@ void on_particle_change()
   invalidate_obs();
 
   /* the particle information is no longer valid */
-  free(partCfg); partCfg=NULL;
+  freePartCfg();
 }
 
 void on_coulomb_change()
@@ -351,8 +352,10 @@ void on_ghost_flags_change()
   if (thermo_switch & THERMO_DPD)
     /* DPD needs also ghost velocities */
     ghosts_have_v = 1;
-  else if (n_rigidbonds != 0)
+#ifdef CONSTRAINT_BOND
+  else if (n_rigidbonds)
     ghosts_have_v = 1;
+#endif
 #ifdef ELECTROSTATICS
   /* Maggs electrostatics needs ghost velocities too */
   else if(coulomb.method == COULOMB_MAGGS)
@@ -413,10 +416,6 @@ static void init_tcl(Tcl_Interp *interp)
   Tcl_CreateCommand(interp, "nemd", (Tcl_CmdProc *)nemd, 0, NULL);
   /* in thermostat.c */
   Tcl_CreateCommand(interp, "thermostat", (Tcl_CmdProc *)thermostat, 0, NULL);
-#ifdef EXCLUSIONS
-  /* in interaction_data.c */
-  Tcl_CreateCommand(interp, "exclusion", (Tcl_CmdProc *)exclusion, 0, NULL);
-#endif
 
   /* evaluate the Tcl initialization script */
   scriptdir = getenv("ESPRESSO_SCRIPTS");
diff --git a/interaction_data.c b/interaction_data.c
index 75f8c28fab6..744dbd9c342 100644
--- a/interaction_data.c
+++ b/interaction_data.c
@@ -42,8 +42,6 @@
  *****************************************/
 int n_particle_types = 0;
 int n_interaction_types = 0;
-int n_rigidbonds = 0;
-int ia_excl = -1;
 IA_parameters *ia_params = NULL;
 
 #ifdef ELECTROSTATICS
@@ -1387,85 +1385,3 @@ int inter(ClientData _data, Tcl_Interp *interp,
   /* check for background errors which have not been handled so far */
   return mpi_gather_runtime_errors(interp, err_code);
 }
-
-#ifdef EXCLUSIONS
-int exclusion_parse(Tcl_Interp * interp, int argc, char ** argv)
-{
-  int i;
-  if(n_total_particles==0)
-  {
-    Tcl_AppendResult(interp, "ERROR!! Non bonded Interaction Exclusions could not be set as there are no particles in the system.\n", (char *) NULL);
-    Tcl_AppendResult(interp, "USAGE: inter exclusion {set|delete}", (char *) NULL);
-    return TCL_ERROR;
-  }
-  if(argc==0)
-  {
-    Tcl_AppendResult(interp, "ERROR!! Non bonded Interaction Exclusions could not be set.", (char *) NULL);
-    Tcl_AppendResult(interp, "USAGE: inter exclusion {set|delete}", (char *) NULL);
-    return TCL_ERROR;
-  }
-  else
-  {
-    for (i=0;i<argc;i++)
-    {
-      if      (ARG0_IS_S("set"))
-      {
-         /**Excludes non bonded interactions 1-2 and 1-3 only*/
-         /**parameter 4 in the function assumes that each particle can have a maximum of 4 distinguishable bond partners if only 1-2 and 1-3 interaction exclusions are considered.*/
-	 ia_excl = create_bond_list(4);
-	 mpi_bcast_parameter(FIELD_IA_EXCL);
-	 mpi_bcast_bond_partners(ia_excl, CURR_BOND_PARTNERS);
-	 i++;argc--;argv++;break;
-      }
-      else if (ARG0_IS_S("delete"))
-      {
-	 ia_excl = -1;
-	 mpi_bcast_parameter(FIELD_IA_EXCL);
-	 mpi_bcast_bond_partners(0, REALLOC_BOND_PARTNERS);
-	 i++;argc--;argv++;break;
-      }
-      else
-      {
-  	 Tcl_AppendResult(interp, "ERROR: Unknown 'interaction exclusion' argument \"", argv[i], "\"", (char *) NULL);
-	 return TCL_ERROR;
-      }
-    }
-
-  }
-
-  return TCL_OK;
-}
-
-int exclusion_print(Tcl_Interp *interp)
-{
-  int start=1;
-  if(ia_excl != -1) {
-      if (start) {
-        Tcl_AppendResult(interp, "{", (char *)NULL);
-        start = 0;
-      }
-      else
-        Tcl_AppendResult(interp, " {", (char *)NULL);
-      Tcl_AppendResult(interp, "exclusion set");
-      Tcl_AppendResult(interp, "}", (char *)NULL);
-  }
-  return (TCL_OK);
-}
-#endif
-
-int exclusion(ClientData _data, Tcl_Interp *interp, int argc, char **argv)
-{
-#ifdef EXCLUSIONS
-  if(n_total_particles==0)
-  {
-    Tcl_AppendResult(interp, "ERROR!! Non bonded Interaction Exclusions could be set only after all particles are set.\n", (char *) NULL);
-    Tcl_AppendResult(interp, "USAGE: exclude_nonbonded_ia {set|delete}", (char *) NULL);
-    return TCL_ERROR;
-  }
-  if(argc==1) return exclusion_print(interp);
-  else        return exclusion_parse(interp, argc-1, argv+1);
-#else
-  Tcl_AppendResult(interp, "EXCLUSIONS not compiled in", (char *) NULL);
-  return TCL_ERROR;
-#endif
-}
diff --git a/interaction_data.h b/interaction_data.h
index 605c8ff8eb8..6ab927e0dbb 100644
--- a/interaction_data.h
+++ b/interaction_data.h
@@ -59,7 +59,6 @@
 
 /*@}*/
 
-
 /** \name Type codes for the type of Coulomb interaction
     Enumeration of implemented methods for the electrostatic
     interaction.
@@ -108,9 +107,6 @@
 
 /*@}*/
 
-#define REALLOC_BOND_PARTNERS 0
-#define CURR_BOND_PARTNERS    1
-
 /* Data Types */
 /************************************************************/
 
@@ -426,9 +422,6 @@ extern int n_particle_types;
 /* Number of nonbonded (short range) interactions. Not used so far.*/
 extern int n_interaction_types;
 
-/**Contains the number of rigid bonds set so far*/
-extern int n_rigidbonds;
-
 /** Structure containing the coulomb parameters. */
 extern Coulomb_parameters coulomb;
 
@@ -504,12 +497,6 @@ int inter(ClientData data, Tcl_Interp *interp,
 int constraint(ClientData _data, Tcl_Interp *interp,
 	       int argc, char **argv);
 
-/** Implementation of the Tcl function exclusion. This function
-    allows to set and delete exclusions of non bonded interactions between bonded particles
-    (currently supported exclusions are for only 1-2 and 1-3 interactions).
- */
-int exclusion(ClientData _data, Tcl_Interp *interp, int argc, char **argv);
-
 /** Callback for setmd niatypes. */
 int niatypes_callback(Tcl_Interp *interp, void *data);
 
diff --git a/particle_data.c b/particle_data.c
index 79cdc243162..141ffd64a35 100644
--- a/particle_data.c
+++ b/particle_data.c
@@ -55,7 +55,29 @@ int max_local_particles = 0;
 Particle **local_particles = NULL;
 Particle *partCfg = NULL;
 int partCfgSorted = 0;
-int *partBondPartners = NULL;
+
+/** bondlist for partCfg, if bonds are needed */
+IntList partCfg_bl = { NULL, 0, 0 };
+
+/************************************************
+ * local functions
+ ************************************************/
+
+/** Remove bond from particle if possible */
+int try_delete_bond(Particle *part, int *bond);
+
+/** Remove exclusion from particle if possible */
+void try_delete_exclusion(Particle *part, int part2);
+
+/** Insert an exclusion if not already set */
+void try_add_exclusion(Particle *part, int part2);
+
+/** Automatically add the next \<distance\> neighbors in each molecule to the exclusion list.
+    This uses the bond topology obtained directly from the particles, since only this contains
+    the full topology, in contrast to \ref topology. To easily setup the bonds, all data
+    should be on a single node, therefore the \ref partCfg array is used. With large amounts
+    of particles, you should avoid this function and setup exclusions manually. */
+void auto_exclusion(int distance);
 
 /************************************************
  * particle initialization functions
@@ -130,10 +152,16 @@ void init_particle(Particle *part)
 #endif
 
   init_intlist(&(part->bl));
+#ifdef EXCLUSIONS
+  init_intlist(&(part->el));
+#endif
 }
 
 void free_particle(Particle *part) {
   realloc_intlist(&(part->bl), 0);
+#ifdef EXCLUSIONS
+  realloc_intlist(&(part->el), 0);
+#endif
 }
 
 
@@ -141,22 +169,19 @@ void free_particle(Particle *part) {
  * organizational functions
  ************************************************/
 
-void updatePartCfg(int bonds_flag )
+void updatePartCfg(int bonds_flag)
 {
   int j;
-  IntList bl;
 
   if(partCfg)
     return;
 
-
   partCfg = malloc(n_total_particles*sizeof(Particle));
-  if ( bonds_flag != WITH_BONDS ){
+  if (bonds_flag != WITH_BONDS)
     mpi_get_particles(partCfg, NULL);
-  }
-  else {
-    mpi_get_particles(partCfg,&bl);
-  }
+  else
+    mpi_get_particles(partCfg,&partCfg_bl);
+ 
   for(j=0; j<n_total_particles; j++)
     unfold_position(partCfg[j].r.p,partCfg[j].l.i);
 
@@ -188,6 +213,13 @@ int sortPartCfg()
   return 1;
 }
 
+void freePartCfg()
+{
+  free(partCfg);
+  partCfg = NULL;
+  realloc_intlist(&partCfg_bl, 0);
+}
+
 /** resize \ref local_particles.
     \param part the highest existing particle
 */
@@ -269,37 +301,6 @@ void update_local_particles(ParticleList *pl)
     local_particles[p[i].p.identity] = &p[i];
 }
 
-int try_delete_bond(Particle *part, int *bond)
-{
-  IntList *bl = &part->bl;
-  int i, j, type, partners;
-  if (!bond) {
-    realloc_intlist(bl, bl->n = 0);
-    return TCL_OK;
-  }
-
-  for (i = 0; i < bl->n;) {
-    type = bond[i];
-    partners = bonded_ia_params[type].num;
-    if (type != bond[0])
-      i += 1 + partners;
-    else {
-      for(j = 1; j <= partners; j++)
-	if (bond[j] != bl->e[i + j])
-	  break;
-      if (j > partners) {
-	bl->n -= 1 + partners;
-	memcpy(bl->e + i, bl->e + i + 1 + partners,
-	       bl->n - i);
-	realloc_intlist(bl, bl->n);
-	return TCL_OK;
-      }
-      i += 1 + partners;
-    }
-  }
-  return TCL_ERROR;
-}
-
 Particle *got_particle(ParticleList *l, int id)
 {
   int i;
@@ -386,72 +387,67 @@ Particle *move_indexed_particle(ParticleList *dl, ParticleList *sl, int i)
 }
 
 #ifdef ROTATION
-void part_print_omega(Particle part, char *buffer, Tcl_Interp *interp)
+void part_print_omega(Particle *part, char *buffer, Tcl_Interp *interp)
 {
-  Tcl_PrintDouble(interp, part.m.omega[0], buffer);
+  Tcl_PrintDouble(interp, part->m.omega[0], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.m.omega[1], buffer);
+  Tcl_PrintDouble(interp, part->m.omega[1], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.m.omega[2], buffer);
+  Tcl_PrintDouble(interp, part->m.omega[2], buffer);
   Tcl_AppendResult(interp, buffer, (char *)NULL);
-  return;
 }
 
-void part_print_torque(Particle part, char *buffer, Tcl_Interp *interp)
+void part_print_torque(Particle *part, char *buffer, Tcl_Interp *interp)
 {
-  Tcl_PrintDouble(interp, part.f.torque[0], buffer);
+  Tcl_PrintDouble(interp, part->f.torque[0], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.f.torque[1], buffer);
+  Tcl_PrintDouble(interp, part->f.torque[1], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.f.torque[2], buffer);
+  Tcl_PrintDouble(interp, part->f.torque[2], buffer);
   Tcl_AppendResult(interp, buffer, (char *)NULL);
-  return;
 }
 
-void part_print_quat(Particle part, char *buffer, Tcl_Interp *interp)
+void part_print_quat(Particle *part, char *buffer, Tcl_Interp *interp)
 {
-  Tcl_PrintDouble(interp, part.r.quat[0], buffer);
+  Tcl_PrintDouble(interp, part->r.quat[0], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.r.quat[1], buffer);
+  Tcl_PrintDouble(interp, part->r.quat[1], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.r.quat[2], buffer);
+  Tcl_PrintDouble(interp, part->r.quat[2], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.r.quat[3], buffer);
+  Tcl_PrintDouble(interp, part->r.quat[3], buffer);
   Tcl_AppendResult(interp, buffer, (char *)NULL);
-  return;
 }
 #endif
 
-void part_print_v(Particle part, char *buffer, Tcl_Interp *interp)
+void part_print_v(Particle *part, char *buffer, Tcl_Interp *interp)
 {
   /* unscale velocities ! */
-  Tcl_PrintDouble(interp, part.m.v[0]/time_step, buffer);
+  Tcl_PrintDouble(interp, part->m.v[0]/time_step, buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.m.v[1]/time_step, buffer);
+  Tcl_PrintDouble(interp, part->m.v[1]/time_step, buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.m.v[2]/time_step, buffer);
+  Tcl_PrintDouble(interp, part->m.v[2]/time_step, buffer);
   Tcl_AppendResult(interp, buffer, (char *)NULL);
-  return;
 }
 
-void part_print_f(Particle part, char *buffer, Tcl_Interp *interp)
+void part_print_f(Particle *part, char *buffer, Tcl_Interp *interp)
 {
   /* unscale forces ! */
-  Tcl_PrintDouble(interp, part.f.f[0]/(0.5*time_step*time_step), buffer);
+  Tcl_PrintDouble(interp, part->f.f[0]/(0.5*time_step*time_step), buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.f.f[1]/(0.5*time_step*time_step), buffer);
+  Tcl_PrintDouble(interp, part->f.f[1]/(0.5*time_step*time_step), buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.f.f[2]/(0.5*time_step*time_step), buffer);
+  Tcl_PrintDouble(interp, part->f.f[2]/(0.5*time_step*time_step), buffer);
   Tcl_AppendResult(interp, buffer, (char *)NULL);
-  return;
 }
 
-void part_print_position(Particle part, char *buffer, Tcl_Interp *interp)
+void part_print_position(Particle *part, char *buffer, Tcl_Interp *interp)
 {
   double ppos[3];
   int img[3];
-  memcpy(ppos, part.r.p, 3*sizeof(double));
-  memcpy(img, part.l.i, 3*sizeof(int));
+  memcpy(ppos, part->r.p, 3*sizeof(double));
+  memcpy(img, part->l.i, 3*sizeof(int));
   unfold_position(ppos, img);
   Tcl_PrintDouble(interp, ppos[0], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
@@ -459,66 +455,73 @@ void part_print_position(Particle part, char *buffer, Tcl_Interp *interp)
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
   Tcl_PrintDouble(interp, ppos[2], buffer);
   Tcl_AppendResult(interp, buffer, (char *)NULL);
-  return;
 }
 
-void part_print_folded_position(Particle part, char *buffer, Tcl_Interp *interp)
+void part_print_folded_position(Particle *part, char *buffer, Tcl_Interp *interp)
 {
-  Tcl_PrintDouble(interp, part.r.p[0], buffer);
+  Tcl_PrintDouble(interp, part->r.p[0], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.r.p[1], buffer);
+  Tcl_PrintDouble(interp, part->r.p[1], buffer);
   Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-  Tcl_PrintDouble(interp, part.r.p[2], buffer);
+  Tcl_PrintDouble(interp, part->r.p[2], buffer);
   Tcl_AppendResult(interp, buffer, (char *)NULL);
-  return;
 }
 
-void part_print_bonding_structure(Particle part, char *buffer, Tcl_Interp *interp, IntList *bl)
+void part_print_bonding_structure(Particle *part, char *buffer, Tcl_Interp *interp)
 {
-  int i=0,j,size;
+  int i,j,size;
   Tcl_AppendResult(interp, " { ", (char *)NULL);
-  while(i<bl->n) {
-    size = bonded_ia_params[bl->e[i]].num;
-    sprintf(buffer, "{%d ", bl->e[i]); i++;
+  i = 0;
+  while(i<part->bl.n) {
+    size = bonded_ia_params[part->bl.e[i]].num;
+    sprintf(buffer, "{%d ", part->bl.e[i]); i++;
     Tcl_AppendResult(interp, buffer, (char *)NULL);
     for(j=0;j<size-1;j++) {
-      sprintf(buffer, "%d ", bl->e[i]); i++;
+      sprintf(buffer, "%d ", part->bl.e[i]); i++;
       Tcl_AppendResult(interp, buffer, (char *)NULL);
     }
-    sprintf(buffer, "%d} ", bl->e[i]); i++;
+    sprintf(buffer, "%d} ", part->bl.e[i]); i++;
     Tcl_AppendResult(interp, buffer, (char *)NULL);
   }
   Tcl_AppendResult(interp, "} ", (char *)NULL);
-  return;
 }
 
+#ifdef EXCLUSIONS
+void part_print_exclusions(Particle *part, char *buffer, Tcl_Interp *interp)
+{
+  int i;
+  for (i = 0; i < part->el.n; i++) {
+    sprintf(buffer, "%d ", part->el.e[i]);
+    Tcl_AppendResult(interp, buffer, (char *)NULL);
+  }
+}
+#endif
+
 #ifdef EXTERNAL_FORCES
-void part_print_fix(Particle part, char *buffer, Tcl_Interp *interp)
+void part_print_fix(Particle *part, char *buffer, Tcl_Interp *interp)
 {
   int i;
   for (i = 0; i < 3; i++) {
-    if (part.l.ext_flag & COORD_FIXED(i))
+    if (part->l.ext_flag & COORD_FIXED(i))
       Tcl_AppendResult(interp, "1 ", (char *)NULL);
     else
 	    Tcl_AppendResult(interp, "0 ", (char *)NULL);
   }
-  return;
 }
 
-void part_print_ext_force(Particle part, char *buffer, Tcl_Interp *interp)
+void part_print_ext_force(Particle *part, char *buffer, Tcl_Interp *interp)
 {
-  if(part.l.ext_flag & PARTICLE_EXT_FORCE) {
-    Tcl_PrintDouble(interp, part.l.ext_force[0], buffer);
+  if(part->l.ext_flag & PARTICLE_EXT_FORCE) {
+    Tcl_PrintDouble(interp, part->l.ext_force[0], buffer);
 	  Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-	  Tcl_PrintDouble(interp, part.l.ext_force[1], buffer);
+	  Tcl_PrintDouble(interp, part->l.ext_force[1], buffer);
 	  Tcl_AppendResult(interp, buffer, " ", (char *)NULL);
-	  Tcl_PrintDouble(interp, part.l.ext_force[2], buffer);
+	  Tcl_PrintDouble(interp, part->l.ext_force[2], buffer);
 	  Tcl_AppendResult(interp, buffer, (char *)NULL);
   }
   else {
     Tcl_AppendResult(interp, "0.0 0.0 0.0 ", (char *)NULL);
   }
-  return;
 }
 #endif
 
@@ -529,7 +532,6 @@ int printParticleToResult(Tcl_Interp *interp, int part_num)
 {
   char buffer[50 + TCL_DOUBLE_SPACE + TCL_INTEGER_SPACE];
   Particle part;
-  IntList *bl = &(part.bl);
 
   if (get_particle_data(part_num, &part) == TCL_ERROR)
     return (TCL_ERROR);
@@ -537,7 +539,7 @@ int printParticleToResult(Tcl_Interp *interp, int part_num)
   sprintf(buffer, "%d", part.p.identity);
   Tcl_AppendResult(interp, buffer, (char *)NULL);
   Tcl_AppendResult(interp, " pos ", (char *)NULL);
-  part_print_position(part, buffer, interp);
+  part_print_position(&part, buffer, interp);
   Tcl_AppendResult(interp, " type ", (char *)NULL);
   sprintf(buffer, "%d", part.p.type);
   Tcl_AppendResult(interp, buffer, " molecule ", (char *)NULL);
@@ -551,40 +553,47 @@ int printParticleToResult(Tcl_Interp *interp, int part_num)
   Tcl_PrintDouble(interp, part.p.q, buffer);
 #endif
   Tcl_AppendResult(interp, buffer, " v ", (char *)NULL);
-  part_print_v(part, buffer, interp);
+  part_print_v(&part, buffer, interp);
   Tcl_AppendResult(interp, " f ", (char *)NULL);
-  part_print_f(part, buffer, interp);
+  part_print_f(&part, buffer, interp);
 
 #ifdef ROTATION
   /* print information about rotation */
   Tcl_AppendResult(interp, " quat ", (char *)NULL);
-  part_print_quat(part, buffer, interp);
+  part_print_quat(&part, buffer, interp);
 
   Tcl_AppendResult(interp, " omega ", (char *)NULL);
-  part_print_omega(part, buffer, interp);
+  part_print_omega(&part, buffer, interp);
 
   Tcl_AppendResult(interp, " torque ", (char *)NULL);
-  part_print_torque(part, buffer, interp);
+  part_print_torque(&part, buffer, interp);
+#endif
+
+#ifdef EXCLUSIONS
+  if (part.el.n > 0) {
+    Tcl_AppendResult(interp, buffer, " exclude ", (char *)NULL);
+    part_print_exclusions(&part, buffer, interp);
+  }
 #endif
 
 #ifdef EXTERNAL_FORCES
   /* print external force information. */
   if (part.l.ext_flag & PARTICLE_EXT_FORCE) {
     Tcl_AppendResult(interp, " ext_force ", (char *)NULL);
-    part_print_ext_force(part, buffer, interp);
+    part_print_ext_force(&part, buffer, interp);
   }
 
   /* print fix information. */
   if (part.l.ext_flag & COORDS_FIX_MASK) {
     Tcl_AppendResult(interp, " fix ", (char *)NULL);
-    part_print_fix( part, buffer, interp);
+    part_print_fix(&part, buffer, interp);
   }
 #endif
 
   /* print bonding structure */
-  if(bl->n > 0) {
-    part_print_bonding_structure(part, buffer, interp, bl);
-  }
+  if (part.bl.n > 0)
+    part_print_bonding_structure(&part, buffer, interp);
+
   free_particle(&part);
   return (TCL_OK);
 }
@@ -624,7 +633,6 @@ int part_parse_print(Tcl_Interp *interp, int argc, char **argv,
 
   char buffer[TCL_DOUBLE_SPACE + TCL_INTEGER_SPACE];
   Particle part;
-  IntList *bl = &(part.bl);
 
   if (part_num > max_seen_particle) {
     Tcl_AppendResult(interp, "na", (char *)NULL);
@@ -643,11 +651,11 @@ int part_parse_print(Tcl_Interp *interp, int argc, char **argv,
       Tcl_AppendResult(interp, buffer, (char *)NULL);
     }
     else if (ARG0_IS_S("position"))
-      part_print_position(part, buffer, interp);
+      part_print_position(&part, buffer, interp);
     else if (ARG0_IS_S("force"))
-      part_print_f(part, buffer, interp);
+      part_print_f(&part, buffer, interp);
     else if (ARG0_IS_S("folded_position"))
-      part_print_folded_position(part, buffer, interp);
+      part_print_folded_position(&part, buffer, interp);
     else if (ARG0_IS_S("type")) {
       sprintf(buffer, "%d", part.p.type);
       Tcl_AppendResult(interp, buffer, (char *)NULL);
@@ -670,25 +678,31 @@ int part_parse_print(Tcl_Interp *interp, int argc, char **argv,
     }
 #endif
     else if (ARG0_IS_S("v"))
-      part_print_v(part, buffer, interp);
+      part_print_v(&part, buffer, interp);
 
 #ifdef ROTATION
     else if (ARG0_IS_S("quat"))
-      part_print_quat(part, buffer, interp);
+      part_print_quat(&part, buffer, interp);
     else if (ARG0_IS_S("omega"))
-      part_print_omega(part, buffer, interp);
+      part_print_omega(&part, buffer, interp);
     else if (ARG0_IS_S("torque"))
-      part_print_torque(part, buffer, interp);
+      part_print_torque(&part, buffer, interp);
 #endif
 
 #ifdef EXTERNAL_FORCES
     else if (ARG0_IS_S("ext_force"))
-      part_print_ext_force(part,buffer,interp);
+      part_print_ext_force(&part,buffer,interp);
     else if (ARG0_IS_S("fix"))
-      part_print_fix(part, buffer, interp);
+      part_print_fix(&part, buffer, interp);
+#endif
+
+#ifdef EXCLUSIONS
+    else if (ARG0_IS_S("exclusions"))
+      part_print_exclusions(&part, buffer, interp);
 #endif
+
     else if (ARG0_IS_S("bonds"))
-      part_print_bonding_structure(part, buffer, interp, bl);
+      part_print_bonding_structure(&part, buffer, interp);
     else {
       Tcl_ResetResult(interp);
       Tcl_AppendResult(interp, "unknown particle data \"", argv[0], "\" requested", (char *)NULL);
@@ -1171,7 +1185,7 @@ int part_parse_bond(Tcl_Interp *interp, int argc, char **argv,
   } else {
     if(type_num < 0 || type_num >= n_bonded_ia) {
       Tcl_AppendResult(interp, "invalid bonded interaction type_num"
-		       "(Set bonded interaction parameters first)", (char *) NULL);
+		       " (set bonded interaction parameters first)", (char *) NULL);
       return TCL_ERROR;
     }
     /* check partners */
@@ -1218,6 +1232,62 @@ int part_parse_bond(Tcl_Interp *interp, int argc, char **argv,
   return TCL_OK;
 }
 
+#ifdef EXCLUSIONS
+int part_parse_exclusion(Tcl_Interp *interp, int argc, char **argv,
+			 int part_num, int * change)
+{
+  int delete = 0, partner;
+
+  *change = 0;
+
+  /* check number of arguments */
+  if (argc < 1) {
+    Tcl_AppendResult(interp, "exclusion requires at least 1 arguments: "
+		     "[delete] { <partner 1> <partner 2> ... }", (char *) NULL);
+    return TCL_ERROR;
+  }
+
+  /* parse away delete eventually */
+  if (ARG0_IS_S("delete")) {
+    delete = 1;
+    argc--;
+    argv++;
+    (*change)++;
+
+    if (argc < 1) {
+      Tcl_AppendResult(interp, "exclusion requires at least 1 arguments: "
+		       "[delete] { <partner 1> <partner 2> ... }", (char *) NULL);
+      return TCL_ERROR;
+    }
+  }
+
+  /* parse partners */
+  if (!particle_node)
+    build_particle_node();
+
+  while (argc > 0) {
+    if (!ARG0_IS_I(partner)) {
+      /* seems to be the next property */
+      Tcl_ResetResult(interp);
+      break;
+    }
+
+    /* set/delete exclusion */
+    if (change_exclusion(part_num, partner, delete) != TCL_OK) {
+      if (delete)
+	Tcl_AppendResult(interp, "exclusion to delete did not exist", (char *)NULL);
+      else
+	Tcl_AppendResult(interp, "particle to exclude from interaction does not exist or is the same", (char *)NULL);	
+      return TCL_ERROR;
+    }
+
+    argc--; argv++;
+    (*change)++;
+  }
+  return TCL_OK;
+}
+#endif
+
 int part_parse_cmd(Tcl_Interp *interp, int argc, char **argv,
 		   int part_num)
 {
@@ -1285,6 +1355,11 @@ int part_parse_cmd(Tcl_Interp *interp, int argc, char **argv,
     else if (ARG0_IS_S("bond"))
       err = part_parse_bond(interp, argc-1, argv+1, part_num, &change);
 
+#ifdef EXCLUSIONS
+    else if (ARG0_IS_S("exclude"))
+      err = part_parse_exclusion(interp, argc-1, argv+1, part_num, &change);
+#endif
+
     else {
       Tcl_AppendResult(interp, "unknown particle parameter \"",
 		       argv[0],"\"", (char *)NULL);
@@ -1313,20 +1388,44 @@ int part(ClientData data, Tcl_Interp *interp,
   /* if no further arguments are given, print out all stored particles */
   if (argc == 1)
     return part_print_all(interp);
+  
+  /* parse away all non particle number arguments */
+  if (ARG1_IS_S("deleteall")) {
+    if (argc != 2) {
+      Tcl_AppendResult(interp, "part deleteall takes no arguments", (char *)NULL);
+      return TCL_ERROR;
+    }
+    remove_all_particles();
+    return TCL_OK;
+  }
+#ifdef EXCLUSIONS
+  else if (ARG1_IS_S("delete_exclusions")) {
+    if (argc != 2) {
+      Tcl_AppendResult(interp, "part delete_exclusions takes no arguments", (char *)NULL);
+      return TCL_ERROR;
+    }
+    remove_all_exclusions();
+    return TCL_OK;
+  }
+  else if (ARG1_IS_S("auto_exclusions")) {
+    int bond_neighbors = 2;
+    if (argc > 3) {
+      Tcl_AppendResult(interp, "usage: part auto_exclusions {distance (bond neighbors to exclude, defaults to 2)}", (char *)NULL);
+      return TCL_ERROR;
+    }
+    if (argc == 3 && !ARG_IS_I(2, bond_neighbors))
+      return TCL_ERROR;
+    auto_exclusion(bond_neighbors);
+    return TCL_OK;
+  }
+#endif
 
   /* only one argument is given */
   if (argc == 2) {
 
-    if (ARG1_IS_S("deleteall")) {
-      remove_all_particles();
-
-      return TCL_OK;
-    }
-
     if (ARG1_IS_I(part_num)) {
       if (printParticleToResult(interp, part_num) == TCL_ERROR)
 	Tcl_AppendResult(interp, "na", (char *)NULL);
-
       return TCL_OK;
     }
 
@@ -1597,16 +1696,13 @@ int change_particle_bond(int part, int *bond, int delete)
 
   if (pnode == -1)
     return TCL_ERROR;
-  if(delete != 0 || bond==NULL) {
-    delete = 1; bond=NULL; }
-  else {
-    if(bond[0] < 0 ) {
-      fprintf(stderr, "ERROR: Failed upon changing bonds of particle %d due to invalid bonded interaction type_num %d (must be >0)!\n", part,bond[0]);
-      return TCL_ERROR; }
-    if(bond[0] >= n_bonded_ia) {
-      fprintf(stderr, "ERROR: Failed upon changing bonds of particle %d due to currently unknown bonded interaction %d!\n", part,bond[0]);
-      fprintf(stderr, "       Specify its parameters first (using 'inter <bond_type_number> { FENE | Angle | ... } <parameters> ...') "
-	      "before adding bonds of that type!\n");
+  if(delete != 0 || bond == NULL)
+    delete = 1;
+
+  if (bond != NULL) {
+    if (bond[0] < 0 || bond[0] >= n_bonded_ia) {
+      char *errtxt = runtime_error(128 + TCL_INTEGER_SPACE);
+      sprintf(errtxt, "invalid/unknown bonded interaction type %d", bond[0]);
       return TCL_ERROR;
     }
   }
@@ -1804,6 +1900,37 @@ int local_change_bond(int part, int *bond, int delete)
   return TCL_OK;
 }
 
+int try_delete_bond(Particle *part, int *bond)
+{
+  IntList *bl = &part->bl;
+  int i, j, type, partners;
+
+  if (!bond) {
+    realloc_intlist(bl, bl->n = 0);
+    return TCL_OK;
+  }
+
+  for (i = 0; i < bl->n;) {
+    type = bl->e[i];
+    partners = bonded_ia_params[type].num;
+    if (type != bond[0])
+      i += 1 + partners;
+    else {
+      for(j = 1; j <= partners; j++)
+	if (bond[j] != bl->e[i + j])
+	  break;
+      if (j > partners) {
+	bl->n -= 1 + partners;
+	memcpy(bl->e + i, bl->e + i + 1 + partners, bl->n - i);
+	realloc_intlist(bl, bl->n);
+	return TCL_OK;
+      }
+      i += 1 + partners;
+    }
+  }
+  return TCL_ERROR;
+}
+
 void remove_all_bonds_to(int identity)
 {
   Cell *cell;
@@ -1841,10 +1968,76 @@ void remove_all_bonds_to(int identity)
   }
 }
 
+#ifdef EXCLUSIONS
+void local_change_exclusion(int part1, int part2, int delete)
+{
+  Cell *cell;
+  int p, np, c;
+  Particle *part;
+
+  if (part1 == -1 && part2 == -1) {
+    /* delete all exclusions */
+    for (c = 0; c < local_cells.n; c++) {
+      cell = local_cells.cell[c];
+      np = cell->n;
+      part = cell->part;
+      for (p = 0; p < np; p++)
+	realloc_intlist(&part[p].el, part[p].el.n = 0);
+    }
+    return;
+  }
+
+  /* part1, if here */
+  part = local_particles[part1];
+  if (part) {
+    if (delete)
+      try_delete_exclusion(part, part2);
+    else
+      try_add_exclusion(part, part2);
+  }
+
+  /* part2, if here */
+  part = local_particles[part2];
+  if (part) {
+    if (delete)
+      try_delete_exclusion(part, part1);
+    else
+      try_add_exclusion(part, part1);
+  }
+}
+
+void try_add_exclusion(Particle *part, int part2)
+{
+  int i;
+  for (i = 0; i < part->el.n; i++)
+    if (part->el.e[i] == part2)
+      return;
+  
+  realloc_intlist(&part->el, part->el.n + 1);
+  part->el.e[part->el.n++] = part2;
+}
+
+void try_delete_exclusion(Particle *part, int part2)
+{
+  IntList *el = &part->el;
+  int i;
+
+  for (i = 0; i < el->n;) {
+    if (el->e[i] == part2) {
+      el->n--;
+      memcpy(el->e + i, el->e + i + 1, el->n - i);
+      realloc_intlist(el, el->n);
+      break;
+    }
+  }
+}
+#endif
+
 void send_particles(ParticleList *particles, int node)
 {
   int pc;
-  IntList local_bi;
+  /* Dynamic data, bonds and exclusions */
+  IntList local_dyn;
 
   PART_TRACE(fprintf(stderr, "%d: send_particles %d to %d\n", this_node, particles->n, node));
 
@@ -1852,22 +2045,30 @@ void send_particles(ParticleList *particles, int node)
   MPI_Send(particles->part, particles->n*sizeof(Particle),
 	   MPI_BYTE, node, REQ_SNDRCV_PART, MPI_COMM_WORLD);
 
-  init_intlist(&local_bi);
+  init_intlist(&local_dyn);
   for (pc = 0; pc < particles->n; pc++) {
     Particle *p = &particles->part[pc];
-    realloc_intlist(&local_bi, local_bi.n + p->bl.n);
-    memcpy(&local_bi.e[local_bi.n], p->bl.e, p->bl.n*sizeof(int));
-    local_bi.n += p->bl.n;
+    int size =  local_dyn.n + p->bl.n;
+#ifdef EXCLUSIONS
+    size += p->el.n;
+#endif
+    realloc_intlist(&local_dyn, size);
+    memcpy(local_dyn.e + local_dyn.n, p->bl.e, p->bl.n*sizeof(int));
+    local_dyn.n += p->bl.n;
+#ifdef EXCLUSIONS
+    memcpy(local_dyn.e + local_dyn.n, p->el.e, p->el.n*sizeof(int));
+    local_dyn.n += p->el.n;
+#endif
   }
 
-  PART_TRACE(fprintf(stderr, "%d: send_particles sending %d bond ints\n", this_node, local_bi.n));
-  if (local_bi.n > 0) {
-    MPI_Send(local_bi.e, local_bi.n*sizeof(int),
+  PART_TRACE(fprintf(stderr, "%d: send_particles sending %d bond ints\n", this_node, local_dyn.n));
+  if (local_dyn.n > 0) {
+    MPI_Send(local_dyn.e, local_dyn.n*sizeof(int),
 	     MPI_BYTE, node, REQ_SNDRCV_PART, MPI_COMM_WORLD);
-    realloc_intlist(&local_bi, 0);
+    realloc_intlist(&local_dyn, 0);
   }
 
-  /* remove particles from this nodes local list */
+  /* remove particles from this nodes local list and free data */
   for (pc = 0; pc < particles->n; pc++) {
     local_particles[particles->part[pc].p.identity] = NULL;
     free_particle(&particles->part[pc]);
@@ -1880,7 +2081,7 @@ void recv_particles(ParticleList *particles, int node)
 {
   int transfer, read, pc;
   MPI_Status status;
-  IntList local_bi;
+  IntList local_dyn;
 
   PART_TRACE(fprintf(stderr, "%d: recv_particles from %d\n", this_node, node));
 
@@ -1894,10 +2095,13 @@ void recv_particles(ParticleList *particles, int node)
 	   REQ_SNDRCV_PART, MPI_COMM_WORLD, &status);
   particles->n += transfer;
 
-  local_bi.n = 0;
+  local_dyn.n = 0;
   for (pc = particles->n - transfer; pc < particles->n; pc++) {
     Particle *p = &particles->part[pc];
-    local_bi.n += p->bl.n;
+    local_dyn.n += p->bl.n;
+#ifdef EXCLUSIONS
+    local_dyn.n += p->el.n;
+#endif
 
     PART_TRACE(fprintf(stderr, "%d: recv_particles got particle %d\n", this_node, p->p.identity));
 #ifdef ADDITIONAL_CHECKS
@@ -1910,10 +2114,10 @@ void recv_particles(ParticleList *particles, int node)
 
   update_local_particles(particles);
 
-  PART_TRACE(fprintf(stderr, "%d: recv_particles expecting %d bond ints\n", this_node, local_bi.n));
-  if (local_bi.n > 0) {
-    alloc_intlist(&local_bi, local_bi.n);
-    MPI_Recv(local_bi.e, local_bi.n*sizeof(int), MPI_BYTE, node,
+  PART_TRACE(fprintf(stderr, "%d: recv_particles expecting %d bond ints\n", this_node, local_dyn.n));
+  if (local_dyn.n > 0) {
+    alloc_intlist(&local_dyn, local_dyn.n);
+    MPI_Recv(local_dyn.e, local_dyn.n*sizeof(int), MPI_BYTE, node,
 	     REQ_SNDRCV_PART, MPI_COMM_WORLD, &status);
   }
   read = 0;
@@ -1921,181 +2125,143 @@ void recv_particles(ParticleList *particles, int node)
     Particle *p = &particles->part[pc];
     if (p->bl.n > 0) {
       alloc_intlist(&p->bl, p->bl.n);
-      memcpy(p->bl.e, &local_bi.e[read], p->bl.n*sizeof(int));
+      memcpy(p->bl.e, &local_dyn.e[read], p->bl.n*sizeof(int));
       read += p->bl.n;
     }
     else
       p->bl.e = NULL;
+#ifdef EXCLUSIONS
+    if (p->el.n > 0) {
+      alloc_intlist(&p->el, p->el.n);
+      memcpy(p->el.e, &local_dyn.e[read], p->el.n*sizeof(int));
+      read += p->el.n;
+    }
+    else
+      p->el.e = NULL;
+#endif
   }
-  if (local_bi.n > 0)
-    realloc_intlist(&local_bi, 0);
+  if (local_dyn.n > 0)
+    realloc_intlist(&local_dyn, 0);
 }
 
 #ifdef EXCLUSIONS
-/**truncates the size (=number of bond partners) for each particle's bond list to w*/
-void recreate_bond_list(int w)
-{
-
-  int *tmp_partBondPartners;
-  int  i, j=0;
-  tmp_partBondPartners = (int *) malloc(n_total_particles*w*sizeof(int));
-
-  for(i=0;i<n_total_particles*w;i++)
-    *(tmp_partBondPartners+j)=-1;
-
-
-  for(i=0;i<n_total_particles;i++)
-  {
-    for(j=0;j<w;j++)
-       *(tmp_partBondPartners+i*w+j)= *(partBondPartners+i*(w+1)+j);
-  }
-
-  mpi_bcast_bond_partners(w, REALLOC_BOND_PARTNERS);
-
-  for(i=0;i<n_total_particles*w;i++)
-    *(partBondPartners+i) = *(tmp_partBondPartners+i);
-
-  free(tmp_partBondPartners);
-  tmp_partBondPartners=NULL;
-}
-
-/**Resorts the *partBondPartners array to minimise memory usage. Input parameter: current size of bondlist for each particle
-returns 0 if resort is not possible and 1 if possible. */
-int resort_bond_list(int w)
-{
-  int i, j, check=0;
-  if(w!=1)
-  {
-    for(i=0;i<n_total_particles;i++)
-    {
-      if (*(partBondPartners+i*w+w-1) != -1)
-      {
-        int flag=0;
-	for(j=0;j<w;j++)
-        {
-	   if(*(partBondPartners + *(partBondPartners+i*w+w-1)*w + j) == i)
-	   {
-	     *(partBondPartners + i*w + w - 1) = -1;
-	     check = 1;
-	     flag=1;
-	     break;
-	   }
-        }
-	if(!flag)
-	{
-	  for(j=0;j<w-1;j++)
-	  {
-	    if(*(partBondPartners + *(partBondPartners+i*w+w-1)*w + j) == -1)
-	    {
-	       *(partBondPartners + *(partBondPartners + i*w + w - 1)*w + j) = i;
-	       *(partBondPartners + i*w + w - 1) = -1;
-	       check = 1;
-	       break;
-
-	    }
-
-	  }
-	}
 
-      }
-    }
-  }
-  return(check);
+int change_exclusion(int part1, int part2, int delete)
+{
+  if (!particle_node)
+    build_particle_node();
+
+  if (part1 < 0 || part1 > max_seen_particle ||
+      part2 < 0 || part2 > max_seen_particle ||
+      part1 == part2 ||
+      particle_node[part1] == -1 ||
+      particle_node[part2] == -1)
+    return TCL_ERROR;
+
+  mpi_send_exclusion(part1, part2, delete);
+  return TCL_OK;
 }
 
-/**returns the number of integers(or size) allocated for each particle*/
-int create_bond_list(int w)
+void remove_all_exclusions()
 {
-  Particle *p;
-  int i=0, j, k, l, itr, loc;
-  IntList bl_tmp;
-  Particle *pcfg;
-  pcfg = (Particle *) malloc(n_total_particles*sizeof(Particle));
+  mpi_send_exclusion(-1, -1, 1);
+}
 
-  mpi_get_particles(pcfg, &bl_tmp);
+/* keep a unique list for particle i. Particle j is only added if it is not i
+   and not already in the list. */
+MDINLINE void add_partner(IntList *il, int i, int j, int distance)
+{
+  int k;
+  if (j == i) return;
+  for (k = 0; k < il->n; k += 2)
+    if (il->e[k] == j)
+      return;
+  realloc_intlist(il, il->n + 2);
+  il->e[il->n++] = j;
+  il->e[il->n++] = distance;
+}
 
-  if (n_total_particles != max_seen_particle + 1)
-    return 0;
+void auto_exclusion(int distance)
+{
+  int count, p, i, j, p1, p2, p3, dist1, dist2;
+  Bonded_ia_parameters *ia_params;
+  Particle *part1, *part2;
+  /* partners is a list containing the currently found excluded particles for each particle,
+     and their distance, as a interleaved list */
+  IntList *partners, *ptns1, *ptns2;
+
+  updatePartCfg(WITH_BONDS);
+
+  /* setup bond partners and distance list. Since we need to identify particles via their identity,
+     we use a full sized array */
+  partners    = malloc((max_seen_particle + 1)*sizeof(IntList));
+  for (p = 0; p <= max_seen_particle; p++)
+    init_intlist(&partners[p]);
 
-  mpi_bcast_bond_partners(w, REALLOC_BOND_PARTNERS);
-  for(i=0;i<w*n_total_particles;i++)
-    *(partBondPartners+i)=-1;
-
-  for (i=n_total_particles-1;i>=0;i--)
-  {
-    p = &pcfg[i];
-    loc = w*p->p.identity;
-    k=0;
-    j=0;
-      while (k<p->bl.n)
-      {
-        l=0;
-	while(l<bonded_ia_params[p->bl.e[k]].num && bonded_ia_params[p->bl.e[k]].num!=3)
-        {
-	   int condition = 0,jj ;
-	   for(jj=0;jj<j;jj++)
-	   {
-	      condition = condition || (*(partBondPartners + loc+jj) == p->bl.e[k+l+1]);
-	      if(condition) break;
-	   }
-	   if(!condition)
-	   {
-	     *(partBondPartners+loc+j) = p->bl.e[k+l+1];
-	     j++;
-	   }
-
-	   l++;
-        }
-        k+=(bonded_ia_params[p->bl.e[k]].num+1);
+  /* determine initial connectivity */
+  for (p = 0; p < n_total_particles; p++) {
+    part1 = &partCfg[p];
+    ptns1 = &partners[part1->p.identity];
+    for (i = 0; i < part1->bl.n;) {
+      ia_params = &bonded_ia_params[part1->bl.e[i++]];
+      if (ia_params->num == 1) {
+	part2 = &partCfg[part1->bl.e[i++]];
+	ptns2 = &partners[part2->p.identity];
 
+	/* you never know whether the user does, may bond a particle to itself...? */
+	if (part2->p.identity != part1->p.identity) {
+	  add_partner(ptns1, part1->p.identity, part2->p.identity, 1);
+	  add_partner(ptns2, part2->p.identity, part1->p.identity, 1);
+	}
       }
-  }
-  free(pcfg);
-  pcfg=NULL;
-  itr=1;
-  while(itr)
-  {
-    int check = 1;
-    for(i=w-1;i<n_total_particles;i=i+w)
-    {
-      check = check && (*(partBondPartners+i)==-1);
-      if(!check) break;
+      else
+	i += ia_params->num;
     }
-    if(check)
-    {
-     recreate_bond_list(--w);
-     itr=1;
+  }
+
+  /* calculate transient connectivity. For each of the current neighbors,
+     also exclude their close enough neighbors.
+  */
+  for (count = 1; count < distance; count++) {
+    for (p1 = 0; p1 <= max_seen_particle; p1++) {
+      for (i = 0; i < partners[p1].n; i += 2) {
+	p2 = partners[p1].e[i];
+	dist1 = partners[p1].e[i + 1];
+	if (dist1 > distance) continue;
+	/* loop over all partners of the partner */
+	for (j = 0; j < partners[p2].n; j += 2) {
+	  p3 = partners[p2].e[j];
+	  dist2 = dist1 + partners[p2].e[j + 1];
+	  if (dist2 > distance) continue;
+	  add_partner(&partners[p1], p1, p3, dist2);
+	  add_partner(&partners[p3], p3, p1, dist2);
+	}
+      }
     }
-    else
-     itr=resort_bond_list(w);
+  }
 
+  /* setup the exclusions and clear the arrays. We do not setup the exclusions up there,
+     since on_part_change clears the partCfg, so that we would have to restore it
+     continously. Of course this could be optimized by bundling the exclusions, but this
+     is only done once and the overhead is as much as for setting the bonds, which
+     the user apparently accepted.
+  */
+  for (p = 0; p <= max_seen_particle; p++) {
+    for (j = 0; j < partners[p].n; j++)
+      if (p < partners[p].e[j]) change_exclusion(p, partners[p].e[j], 0);
+    realloc_intlist(&partners[p], 0);
   }
-  return(w);
+  free(partners);
 }
 
-/**returns 1 if i and j are bond partners, otherwise 0*/
-int is_bond_partner(int i, int j)
+int do_nonbonded(Particle *p1, Particle *p2)
 {
-  int check=0, k;
-  if(ia_excl!=-1)
-  {
-    //checking if j is a bond partner of i
-    for(k=0;k<ia_excl;k++)
-    {
-      check = check || (*(partBondPartners + (i*ia_excl) +k) == j);
-      if(check) break;
-    }
-    //if j is not a bond partner of i then check if i is a bond partner of j
-    if(!check)
-    {
-      for(k=0;k<ia_excl;k++)
-      {
-        check = check || (*(partBondPartners + (j*ia_excl) +k) == i);
-        if(check) break;
-      }
-    }
-  }
-  VERLET_TRACE(if(check) fprintf(stderr,"%d: %d-%d pair excluded from verlet list \n", this_node, i, j));
-  return(check);
+  int i, i2;
+  /* check for particle 2 in particle 1's exclusion list. The exclusion list is
+     symmetric, so this is sufficient. */
+  i2  = p2->p.identity;
+  for (i = 0; i < p1->el.n; i++)
+    if (i2 == p1->el.e[i]) return 0;
+  return 1;
 }
 #endif
diff --git a/particle_data.h b/particle_data.h
index 09722f5245f..980cfd4d9e2 100644
--- a/particle_data.h
+++ b/particle_data.h
@@ -148,8 +148,15 @@ typedef struct {
   ///
   ParticleLocal l;
 
-  /** bonded interactions list. */
+  /** bonded interactions list. The format is pretty simple:
+      Just the bond type, and then the particle ids. The number of particle ids can be determined
+      easily from the bonded_ia_params entry for the type. */
   IntList bl;
+
+#ifdef EXCLUSIONS
+  /** list of particles, with which this particle has no nonbonded interactions */
+  IntList el;
+#endif
 } Particle;
 
 /** List of particles. The particle array is resized using a sophisticated
@@ -228,9 +235,6 @@ void init_particle(Particle *part);
 /** Deallocate the dynamic storage of a particle. */
 void free_particle(Particle *part);
 
-/** Remove bond from particle if possible */
-int try_delete_bond(Particle *part, int *bond);
-
 /*    Functions acting on Particle Lists        */
 /************************************************/
 
@@ -416,13 +420,27 @@ int set_particle_fix(int part,  int flag);
 /** Call only on the master node: change particle bond.
     @param part     identity of principal atom of the bond.
     @param bond     field containing the bond type number and the
-    identity of all bond partners (secundary atoms of the bond).
+    identity of all bond partners (secundary atoms of the bond). If NULL, delete all bonds.
     @param delete   if true, do not add the bond, rather delete it if found
     @return TCL_OK on success or TCL_ERROR if no success
     (e. g. particle or bond to delete does not exist)
 */
 int change_particle_bond(int part, int *bond, int delete);
 
+#ifdef EXCLUSIONS
+/** Call only on the master node: change particle constraints.
+    @param part1    identity of particle for which the exclusion is set.
+    @param part2    identity of particle for which the exclusion is set. If -1, delete all exclusions.
+    @param delete   if true, do not add the exclusion, rather delete it if found
+    @return TCL_OK on success or TCL_ERROR if no success
+    (e. g. particles do not exist / did not have exclusion set)
+*/
+int change_exclusion(int part, int part2, int delete);
+
+/** remove all exclusions. */
+void remove_all_exclusions();
+#endif
+
 /** remove particle with a given identity. Also removes all bonds to the particle.
     @param part     identity of the particle to remove
     @return TCL_OK on success or TCL_ERROR if particle does not exist
@@ -449,6 +467,11 @@ void remove_all_bonds_to(int part);
 */
 void updatePartCfg(int bonds_flag );
 
+/** release the partCfg array. Use this function, since it also frees the
+    bonds, if they are used.
+*/
+void freePartCfg();
+
 /** sorts the \ref partCfg array. This is indicated by setting
     \ref partCfgSorted to 1. Note that for this to work the particles
     have to be stored consecutively starting with 0.
@@ -483,6 +506,14 @@ void added_particle(int part);
 */
 int local_change_bond(int part, int *bond, int delete);
 
+/** Used for example by \ref mpi_send_exclusion.
+    Locally add a exclusion to a particle.
+    @param part1 the identity of the first exclusion partner
+    @param part2 the identity of the second exclusion partner
+    @param delete if true, delete the exclusion instead of add
+*/
+void local_change_exclusion(int part1, int part2, int delete);
+
 /** Used by \ref mpi_remove_particle, should not be used elsewhere.
     Remove a particle on this node.
     @param part the identity of the particle to remove
@@ -512,26 +543,10 @@ void send_particles(ParticleList *particles, int node);
 void recv_particles(ParticleList *particles, int node);
 
 #ifdef EXCLUSIONS
-/** Determines if particle ids i and j are connected by bonds
-   (e.g FENE, HAMONIC, RIGID_BOND etc ) and angles (e.g BOND_ANGLE_HARMONIC).
-    Note: Can be used only if non-bonded interaction exclusions are set, and is
-    probably pretty slow
-*/
-int is_bond_partner(int i, int j);
+/** Determines if the non bonded interactions between p1 and p2 should be calculated */
+int do_nonbonded(Particle *p1, Particle *p2);
 #endif
 
-/** Invoked from \ref "create_bond_list(int w)" **/
-void recreate_bond_list(int w);
-
-/** Invoked from \ref "create_bond_list(int w)" */
-int  resort_bond_list(int w);
-
-/** Creates 1D array \ref partBondPartners of minimum size to store the bond
-    partners of each particle. For example; If particle j and k are bond partners
-    of i then partBondPartners[i]=j and partBondPartners[i+1]=k. But,
-    partBondPartners[j] and partBondPartners[k] do not contain i. */
-int  create_bond_list(int w);
-
 /** Complain about a missing bond partner. Just for convenience, replaces the old checked_particle_ptr.
     @param id particle identity.
  */
diff --git a/pressure.c b/pressure.c
index 2907fafa738..9bab4fcfe6f 100644
--- a/pressure.c
+++ b/pressure.c
@@ -596,7 +596,7 @@ int calc_p_tensor(double volume, IntList *p_list, int flag)
 	p_tensor.data.e[k*3 + l] += (p1.m.v[k])*(p1.m.v[l])*PMASS(p1);
 
     /* bonded interactions */
-    i=0;
+    i = 0;
     while(i < p1.bl.n) {
       if((flag==1) || intlist_contains(p_list,p1.bl.e[i+1])) {
 	type_num = p1.bl.e[i++];
diff --git a/pressure.h b/pressure.h
index 6c3ca9f78d6..caeb701d38e 100644
--- a/pressure.h
+++ b/pressure.h
@@ -201,7 +201,7 @@ MDINLINE void add_bonded_virials(Particle *p1)
 
   int i, type_num, type;
   
-  i=0;
+  i = 0;
   while(i<p1->bl.n) {
     type_num = p1->bl.e[i++];
     iaparams = &bonded_ia_params[type_num];
diff --git a/rattle.c b/rattle.c
index b687660abd9..0a267566f48 100644
--- a/rattle.c
+++ b/rattle.c
@@ -15,6 +15,8 @@
 #include "domain_decomposition.h"
 #include "rattle.h"
 
+int n_rigidbonds = 0;
+
 #ifdef BOND_CONSTRAINT
 
 /** \name Private functions */
@@ -401,28 +403,28 @@ int check_tol_vel()
     p  = cell->part;
     np = cell->n;
     for(i = 0; i < np; i++) {
-       k=0;
-       while(k<p[i].bl.n)
-       {
-	 b_ia = &bonded_ia_params[p[i].bl.e[k++]];
-	 if(b_ia->type == BONDED_IA_RIGID_BOND)
-	   {
-	     Particle *p2 = local_particles[p[i].bl.e[k++]];
-	     if (!p2) {
-	       char *errtxt = runtime_error(128 + 2*TCL_INTEGER_SPACE);
-	       sprintf(errtxt,"{ rigid bond broken between particles %d and %d (particles not stored on the same node)} ",
-		       p[i].p.identity, p[i].bl.e[k-1]);
-	       return 0;
-	     }
-
-	     get_mi_vector(r_ij, p[i].r.p , p2->r.p);
-	     vector_subt(v_ij, p[i].m.v , p2->m.v);
-	     tol = fabs(scalar(r_ij, v_ij));
-	     repeat = repeat || (tol > b_ia->p.rigid_bond.v_tol);
-	   }
-	 else
-	   k += b_ia->num;
-       } //while k loop
+      k=0;
+      while(k<p[i].bl.n)
+	{
+	  b_ia = &bonded_ia_params[p[i].bl.e[k++]];
+	  if(b_ia->type == BONDED_IA_RIGID_BOND)
+	    {
+	      Particle *p2 = local_particles[p[i].bl.e[k++]];
+	      if (!p2) {
+		char *errtxt = runtime_error(128 + 2*TCL_INTEGER_SPACE);
+		sprintf(errtxt,"{ rigid bond broken between particles %d and %d (particles not stored on the same node)} ",
+			p[i].p.identity, p[i].bl.e[k-1]);
+		return 0;
+	      }
+	      
+	      get_mi_vector(r_ij, p[i].r.p , p2->r.p);
+	      vector_subt(v_ij, p[i].m.v , p2->m.v);
+	      tol = fabs(scalar(r_ij, v_ij));
+	      repeat = repeat || (tol > b_ia->p.rigid_bond.v_tol);
+	    }
+	  else
+	    k += b_ia->num;
+	} //while k loop
     } //for i loop
  }// for c loop
  return(repeat);
@@ -495,13 +497,13 @@ void correct_vel_shake()
 
 void print_bond_len()
 {
- int i, k, c, np;
- Cell *cell;
- Particle *p;
- Bonded_ia_parameters *b_ia;
- double r_ij[3];
- printf("%d: ", this_node);
- for (c = 0; c < local_cells.n; c++) {
+  int i, k, c, np;
+  Cell *cell;
+  Particle *p;
+  Bonded_ia_parameters *b_ia;
+  double r_ij[3];
+  printf("%d: ", this_node);
+  for (c = 0; c < local_cells.n; c++) {
     cell = local_cells.cell[c];
     p  = cell->part;
     np = cell->n;
@@ -527,8 +529,8 @@ void print_bond_len()
 	       k += b_ia->num;
        } //while k loop
     } //for i loop
- }// for c loop
- printf("\n");
+  }// for c loop
+  printf("\n");
 }
 
 #endif
diff --git a/rattle.h b/rattle.h
index 4e7c94851b0..dd345729a0a 100644
--- a/rattle.h
+++ b/rattle.h
@@ -20,6 +20,9 @@
 #include "global.h"
 #include "particle_data.h"
 
+/** number of rigid bonds */
+extern int n_rigidbonds;
+
 /** Transfers the current particle positions from r.p[3] to r.p_pold[3]
     of the \ref Particle structure. Invoked from \ref correct_pos_shake() */
 void save_old_pos();
diff --git a/verlet.c b/verlet.c
index 38d656a7457..8f181d66ad1 100644
--- a/verlet.c
+++ b/verlet.c
@@ -126,7 +126,7 @@ void build_verlet_lists()
 	/* Loop neighbor cell particles */
 	for(j = j_start; j < np2; j++) {
 #ifdef EXCLUSIONS
-          if(!is_bond_partner(p1[i].p.identity, p2[j].p.identity))
+          if(do_nonbonded(&p1[i], &p2[j]))
 #endif
 	  {
 	    dist2 = distance2(p1[i].r.p, p2[j].r.p);
@@ -227,7 +227,7 @@ void build_verlet_lists_and_calc_verlet_ia()
 	/* Loop neighbor cell particles */
 	for(j = j_start; j < np2; j++) {
 #ifdef EXCLUSIONS
-          if(!is_bond_partner(p1[i].p.identity, p2[j].p.identity))
+          if(do_nonbonded(&p1[i], &p2[j]))
 #endif
 	  {
 	  dist2 = distance2vec(p1[i].r.p, p2[j].r.p, vec21);