private void breakAndMakeBonds(int i, int j) {

		rBond = null;
		RadialBond rb = null;
		synchronized (bonds.getSynchronizationLock()) {
			for (Iterator it = bonds.iterator(); it.hasNext();) {
				rb = (RadialBond) it.next();
				if ((rb.contains(atom[i]) && !rb.contains(atom[j])) || (rb.contains(atom[j]) && !rb.contains(atom[i]))) {
					if (rBond != null)
						return; // both i and j are bonded, their reaction is not permitted.
					rBond = rb;
				}
			}
		}
		if (rBond == null)
			return; // no required bond found

		dpot = 0.0;
		int i1 = rBond.atom1.getIndex();
		int i2 = rBond.atom2.getIndex();
		int k = -1; // the bonded atom not colliding with the radical
		int m = -1; // the atom bonded to k
		int n = -1; // the free radical: m-n the colliding pair
		if (i1 != i && i1 != j) {
			k = i1;
			m = i2;
		}
		else if (i2 != i && i2 != j) {
			k = i2;
			m = i1;
		}
		if (m == i) {
			n = j;
		}
		else if (m == j) {
			n = i;
		}

		if (k != -1 && m != -1 && n != -1) {

			if (atom[k].getID() == atom[n].getID()) {// useless collision

				// if(allowUselessCollision){}

			}

			else {

				// figure out the activation energies for different scenarios of collision.
				ecol = Double.MAX_VALUE;

				if (type instanceof Reaction.A2_B2__2AB) {

					double vAA = type.getParameter(Reaction.A2_B2__2AB.VAA).doubleValue();
					double vBB = type.getParameter(Reaction.A2_B2__2AB.VBB).doubleValue();
					double vAB = type.getParameter(Reaction.A2_B2__2AB.VAB).doubleValue();
					double vA2B = type.getParameter(Reaction.A2_B2__2AB.VA2B).doubleValue();
					double vAB2 = type.getParameter(Reaction.A2_B2__2AB.VAB2).doubleValue();
					int idA = ((Reaction.A2_B2__2AB) type).getIDA();
					int idB = ((Reaction.A2_B2__2AB) type).getIDB();

					if (atom[k].getID() == idA) { // i.e. atom n must be an B
						if (atom[m].getID() == idA) { // A2+B*-->A*+AB
							if (vAB > vAA) {
								ecol = vA2B;
							}
							else {
								ecol = vA2B + vAA - vAB;
							}
							/*
							 * (imagine vA2B=0) If the final state is more stable than the initial one (vAB>vAA), the
							 * transition should be easy. Conversely, if the final state is less stable than the initial
							 * one (vAB<vAA), the transition should be harder. Same applied to the following
							 * conditional branches.
							 */
						}
						else if (atom[m].getID() == idB) { // AB+B*-->A*+B2
							if (vAB > vBB) {
								ecol = vAB2 + vAB - vBB;
							}
							else {
								ecol = vAB2;
							}
						}
					}
					else if (atom[k].getID() == idB) { // i.e. atom n must be an A
						if (atom[m].getID() == idB) { // A*+B2-->AB+B*
							if (vAB > vBB) {
								ecol = vAB2;
							}
							else {
								ecol = vAB2 + vBB - vAB;
							}
						}
						else if (atom[m].getID() == idA) { // A*+AB-->A2+B*
							if (vAB > vAA) {
								ecol = vA2B + vAB - vAA;
							}
							else {
								ecol = vA2B;
							}
						}
					}

				}

				else if (type instanceof Reaction.A2_B2_C__2AB_C) {

					double vAA = type.getParameter(Reaction.A2_B2_C__2AB_C.VAA).doubleValue();
					double vBB = type.getParameter(Reaction.A2_B2_C__2AB_C.VBB).doubleValue();
					double vAB = type.getParameter(Reaction.A2_B2_C__2AB_C.VAB).doubleValue();
					double vAC = type.getParameter(Reaction.A2_B2_C__2AB_C.VAC).doubleValue();
					double vBC = type.getParameter(Reaction.A2_B2_C__2AB_C.VBC).doubleValue();
					double vBA2 = type.getParameter(Reaction.A2_B2_C__2AB_C.VBA2).doubleValue();
					double vAB2 = type.getParameter(Reaction.A2_B2_C__2AB_C.VAB2).doubleValue();
					double vCA2 = type.getParameter(Reaction.A2_B2_C__2AB_C.VCA2).doubleValue();
					double vCB2 = type.getParameter(Reaction.A2_B2_C__2AB_C.VCB2).doubleValue();
					double vABC = type.getParameter(Reaction.A2_B2_C__2AB_C.VABC).doubleValue();
					double vBAC = type.getParameter(Reaction.A2_B2_C__2AB_C.VBAC).doubleValue();
					int idA = ((Reaction.A2_B2_C__2AB_C) type).getIDA();
					int idB = ((Reaction.A2_B2_C__2AB_C) type).getIDB();
					int idC = ((Reaction.A2_B2_C__2AB_C) type).getIDC();

					if (atom[k].getID() == idA) { // i.e. atom n must be an B or C
						if (atom[n].getID() == idB) {
							if (atom[m].getID() == idA) { // B*+A2-->AB+A*
								if (vAB > vAA) {
									ecol = vBA2;
								}
								else {
									ecol = vBA2 + vAA - vAB;
								}
							}
							else if (atom[m].getID() == idB) { // B*+AB-->B2+A*
								if (vBB > vAB) {
									ecol = vAB2;
								}
								else {
									ecol = vAB2 + vAB - vBB;
								}
							}
							else if (atom[m].getID() == idC) { // AC+B*-->A*+BC
								if (vBC > vAC) {
									ecol = (vABC + vBAC) * 0.5;
								}
								else {
									ecol = 0.5 * (vABC + vBAC) + vAC - vBC;
								}
							}
						}
						else if (atom[n].getID() == idC) {
							if (atom[m].getID() == idA) { // A2+C*-->A*+AC
								if (vAC > vAA) {
									ecol = vCA2;
								}
								else {
									ecol = vCA2 + vAA - vAC;
								}
							}
							else if (atom[m].getID() == idB) { // AB+C*-->A*+BC
								if (vBC > vAB) {
									ecol = vABC;
								}
								else {
									ecol = vABC + vAB - vBC;
								}
							}
						}
					}
					else if (atom[k].getID() == idB) { // i.e. atom n = A or C
						if (atom[n].getID() == idA) {
							if (atom[m].getID() == idA) { // A*+AB-->A2+B*
								if (vAA > vAB) {
									ecol = vBA2;
								}
								else {
									ecol = vBA2 + vAB - vAA;
								}
							}
							else if (atom[m].getID() == idB) { // A*+B2-->AB+B*
								if (vAB > vBB) {
									ecol = vAB2;
								}
								else {
									ecol = vAB2 + vBB - vAB;
								}
							}
							else if (atom[m].getID() == idC) { // A*+BC-->AC+B*
								if (vAC > vBC) {
									ecol = 0.5 * (vABC + vBAC);
								}
								else {
									ecol = 0.5 * (vABC + vBAC) + vBC - vAC;
								}
							}
						}
						else if (atom[n].getID() == idC) {
							if (atom[m].getID() == idA) { // C*+AB-->AC+B*
								if (vAC > vAB) {
									ecol = vBAC;
								}
								else {
									ecol = vBAC + vAC - vAB;
								}
							}
							else if (atom[m].getID() == idB) { // C*+B2-->BC+B*
								if (vBC > vBB) {
									ecol = vCB2;
								}
								else {
									ecol = vCB2 + vBB - vBC;
								}
							}
						}
					}
					else if (atom[k].getID() == idC) { // i.e. atom n = A or B
						if (atom[n].getID() == idA) {
							if (atom[m].getID() == idA) { // A*+AC-->A2+C*
								if (vAA > vAC) {
									ecol = vCA2;
								}
								else {
									ecol = vCA2 + vAC - vAA;
								}
							}
							else if (atom[m].getID() == idB) { // A*+BC-->AB+C*
								if (vAB > vBC) {
									ecol = vABC;
								}
								else {
									ecol = vABC + vBC - vAB;
								}
							}
						}
						else if (atom[n].getID() == idB) {
							if (atom[m].getID() == idA) { // B*+AC-->AB+C*
								if (vAB > vAC) {
									ecol = vBAC;
								}
								else {
									ecol = vBAC + vAC - vAB;
								}
							}
							else if (atom[m].getID() == idB) { // B*+BC-->B2+C*
								if (vBB > vBC) {
									ecol = vCB2;
								}
								else {
									ecol = vCB2 + vBC - vBB;
								}
							}
						}
					}

				}

				else if (type instanceof Reaction.O2_2H2__2H2O) {

					double vHH = type.getParameter(Reaction.O2_2H2__2H2O.VHH).doubleValue();
					double vOO = type.getParameter(Reaction.O2_2H2__2H2O.VOO).doubleValue();
					double vHO = type.getParameter(Reaction.O2_2H2__2H2O.VHO).doubleValue();
					double vHO2 = type.getParameter(Reaction.O2_2H2__2H2O.VHO2).doubleValue();
					double vOH2 = type.getParameter(Reaction.O2_2H2__2H2O.VOH2).doubleValue();
					int idH = ((Reaction.O2_2H2__2H2O) type).getIDH();
					int idO = ((Reaction.O2_2H2__2H2O) type).getIDO();

					if (atom[k].getID() == idH) { // i.e. atom n must be an O
						if (atom[m].getID() == idH) { // H2+O*-->H*+HO
							if (vHO > vHH) {
								ecol = vOH2;
							}
							else {
								ecol = vOH2 + vHH - vHO;
							}
							/*
							 * (imagine vOH2=0) If the final state is more stable than the initial one (vHO>vHH), the
							 * transition should be easy. Conversely, if the final state is less stable than the initial
							 * one (vHO<vHH), the transition should be harder. Same applied to the following
							 * conditional branches.
							 */
						}
						else if (atom[m].getID() == idO) { // HO+O*-->H*+O2
							if (vHO > vOO) {
								ecol = vHO2 + vHO - vOO;
							}
							else {
								ecol = vHO2;
							}
						}
					}
					else if (atom[k].getID() == idO) { // i.e. atom n must be
						// an H
						if (atom[m].getID() == idO) { // H*+O2-->HO+O*
							if (vHO > vOO) {
								ecol = vHO2;
							}
							else {
								ecol = vHO2 + vOO - vHO;
							}
						}
						else if (atom[m].getID() == idH) { // H*+HO-->H2+O*
							if (vHO > vHH) {
								ecol = vOH2 + vHO - vHH;
							}
							else {
								ecol = vOH2;
							}
						}
					}

				}

				if (ecol >= Double.MAX_VALUE - 1)
					return;

				/* calculate the line-of-centers energy */
				rijsr = 1.0 / Math.sqrt(rij);
				vxij = atom[i].vx - atom[j].vx;
				vyij = atom[i].vy - atom[j].vy;
				vxy = vxij * xij * rijsr + vyij * yij * rijsr;

				/* use reduced mass to compute head-on kinetic energy */
				if (atom[i].mass * atom[j].mass / (atom[i].mass + atom[j].mass) * vxy * vxy * EV_CONVERTER > ecol) {

					double bx2 = atom[i1].rx - atom[i2].rx;
					bx2 *= bx2;
					double by2 = atom[i1].ry - atom[i2].ry;
					by2 *= by2;
					double br2 = bx2 + by2;
					br2 = Math.sqrt(br2) - rBond.getBondLength();
					dpot += 0.5 * rBond.getBondStrength() * br2 * br2;
					dpot -= rBond.getChemicalEnergy();

					bondLength = getBondLength(atom[i], atom[j]);
					bondStrength = getBondStrength(atom[i], atom[j]);
					bondEnergy = getBondEnergy(i, j);
					br2 = Math.sqrt(rij) - bondLength;
					dpot -= 0.5 * bondStrength * br2 * br2;
					dpot += bondEnergy;

					atom[i].storeCurrentVelocity();
					atom[j].storeCurrentVelocity();
					atom[k].storeCurrentVelocity();
					if (conserve(i, j, k, dpot)) {

						bonds.remove(rBond);
						rBond.destroy();
						removeAssociatedAngularBonds();
						bonds.add(new RadialBond.Builder(atom[i], atom[j]).bondLength(bondLength).bondStrength(
								bondStrength).chemicalEnergy(bondEnergy).build());
						/*
						 * temporarily indicate that the atoms should not be treated as a free radical. CHECK!!! In H2O
						 * reaction, this is not true any more, but in order for the OH fragments not to react with
						 * other fragments, this had better be set.
						 */
						atom[i].setRadical(false);
						atom[j].setRadical(false);
						bondChanged = true;

					}
					else {

						atom[i].restoreVelocity();
						atom[j].restoreVelocity();
						atom[k].restoreVelocity();

					}

				}

			}

		}
		else {
			System.err.println("Error in <bonds>...");
		}

	}