private double computeForceForElectrons(final int time) {
		if (freeElectrons.isEmpty())
			return 0;

		synchronized (freeElectrons) {
			List<Electron> toRemove = null;
			// List<Electron> toAdd = null;
			for (Electron e : freeElectrons) {
				e.fx = 0;
				e.fy = 0;
				for (int i = 0; i < numberOfAtoms; i++) {
					if (atom[i].near(e, atom[i].sigma * 0.5)) {
						if (atom[i].hasElectrons()) { // if the atom contains electrons
							// Electron x = atom[i].collideWithElectron(e);
							// if (x != null) {
							// if (toAdd == null)
							// toAdd = new ArrayList<Electron>();
							// toAdd.add(x);
							// }
						} else { // if the atom is a positive ion, recombine
							if (threeBodyRecombination == null)
								threeBodyRecombination = new ThreeBodyRecombination(this);
							if (threeBodyRecombination.recombine(atom[i], e)) {
								if (toRemove == null)
									toRemove = new ArrayList<Electron>();
								toRemove.add(e);
							}
						}
						break;
					}
				}
			}
			if (toRemove != null)
				freeElectrons.removeAll(toRemove);
			// if (toAdd != null) freeElectrons.addAll(toAdd);
		}

		double vsum = 0;
		double vij = 0;
		double rCD = universe.getCoulombConstant() / universe.getDielectricConstant();
		// interactions between free electrons
		int n = freeElectrons.size();
		if (n >= 2) {
			Electron ei, ej;
			synchronized (freeElectrons) {
				for (int i = 0; i < n - 1; i++) {
					ei = freeElectrons.get(i);
					rxi = ei.rx;
					ryi = ei.ry;
					fxi = ei.fx;
					fyi = ei.fy;
					for (int j = i + 1; j < n; j++) {
						ej = freeElectrons.get(j);
						rxij = rxi - ej.rx;
						ryij = ryi - ej.ry;
						minimumImageConvention();
						rijsq = rxij * rxij + ryij * ryij;
						vij = rCD / Math.sqrt(rijsq);
						vsum += vij;
						fij = vij / rijsq * GF_CONVERSION_CONSTANT;
						fxij = fij * rxij;
						fyij = fij * ryij;
						fxi += fxij;
						fyi += fyij;
						ej.fx -= fxij;
						ej.fy -= fyij;
					}
					ei.fx = fxi;
					ei.fy = fyi;
				}
			}
		}
		// interactions between atoms and free electrons: v(r)=alpha*sigma^11/r^12-Q/r
		double sig11, r12, coul;
		for (int i = 0; i < numberOfAtoms; i++) {
			// if an atom is charged, compute the electrostatic force
			if (Math.abs(atom[i].charge) > ZERO) {
				synchronized (freeElectrons) {
					for (Electron e : freeElectrons) {
						rxij = e.rx - atom[i].rx;
						ryij = e.ry - atom[i].ry;
						minimumImageConvention();
						rijsq = rxij * rxij + ryij * ryij;
						coul = -atom[i].charge / Math.sqrt(rijsq) * rCD;
						sig11 = getPower(atom[i].sigma, 11);
						r12 = getPower(rijsq, 6);
						vij = alpha * sig11 / r12 + coul;
						vsum += vij;
						fij = GF_CONVERSION_CONSTANT * (12 * alpha * sig11 / r12 + coul) / rijsq;
						fxij = fij * rxij;
						fyij = fij * ryij;
						e.fx += fxij;
						e.fy += fyij;
						atom[i].fx -= fxij;
						atom[i].fy -= fyij;
					}
				}
			}
		}
		synchronized (freeElectrons) {
			double imass = 1.0f / Electron.mass;
			for (Electron e : freeElectrons) {
				e.fx *= imass;
				e.fy *= imass;
				for (VectorField f : fields) {
					if (f instanceof ElectricField) {
						ElectricField ef = (ElectricField) f;
						ef.dyn(universe.getDielectricConstant(), e, time);
						vsum += ef.getPotential(e, time);
					} else if (f instanceof MagneticField) {
						((MagneticField) f).dyn(e);
					}
				}
			}
		}
		return vsum;
	}