* @GuardedBy("this")
	 */

	public synchronized double computeForce(final int time) {

		double vsum = 0.0;
		if (time <= 0) {
			updateList = true; // the neighbor list has to be updated in order for the force vector to be plotted
		}

		for (int i = 0; i < numberOfAtoms; i++) {
			atom[i].fx = 0;
			atom[i].fy = 0;
		}
		if (obstacles != null) {
			RectangularObstacle obs;
			synchronized (obstacles.getSynchronizationLock()) {
				if (!obstacles.isEmpty()) {
					for (int jobs = 0, nobs = obstacles.size(); jobs < nobs; jobs++) {
						obs = obstacles.get(jobs);
						obs.ax = 0;
						obs.ay = 0;
					}
				}
			}
		}

		vsum += computeForceForElectrons(time);
		// computeForceForElectrons(time);

		if (numberOfAtoms == 1)
			return computeForceForSingleAtom(time);

		if (updateParArray) {
			resetParArray();
			updateParArray = false;
		}

		xbox = boundary.width;
		ybox = boundary.height;
		double coul = 0.0;
		double rCD = 1.0;
		if (hasCoulomb && interCoulomb)
			rCD = universe.getCoulombConstant() / universe.getDielectricConstant();

		if (updateList) {

			for (int i = 0; i < numberOfAtoms; i++) {
				rx0[i] = atom[i].rx;
				ry0[i] = atom[i].ry;
			}

			nlist = 0;

			for (int i = 0, imax1 = numberOfAtoms - 1; i < imax1; i++) {

				pointer[i] = nlist;
				rxi = atom[i].rx;
				ryi = atom[i].ry;
				fxi = atom[i].fx;
				fyi = atom[i].fy;

				for (int j = i + 1; j < numberOfAtoms; j++) {

					if (!ljBetweenBondPairs && bondTable[i][j])
						continue;

					rxij = rxi - atom[j].rx;
					ryij = ryi - atom[j].ry;
					minimumImageConvention();
					rijsq = rxij * rxij + ryij * ryij;

					if (rijsq < listSquareMatrix[atom[i].id][atom[j].id]) {
						neighborList[nlist++] = j;
					}

					if (rijsq < cutOffSquareMatrix[atom[i].id][atom[j].id]) {

						sr2 = sigab[i][j] / rijsq;
						/* check if this pair gets too close */
						if (sr2 > 2.0) {
							sr2 = 2.0;
							rijsq = 0.5 * sigab[i][j];
						}
						sr6 = sr2 * sr2 * sr2;
						sr12 = sr6 * sr6;

						if (isRepulsive(atom[i].id, atom[j].id)) {
							vij = sr6 * epsab[i][j] * crossRepulsionIntensity;
							wij = vij;
						} else {
							vij = (sr12 - sr6) * epsab[i][j];
							wij = vij + sr12 * epsab[i][j];
						}
						if (cutOffShift) {
							vij -= poten_LJ[atom[i].id][atom[j].id];
							wij -= slope_LJ[atom[i].id][atom[j].id];
						}

						vsum += vij;
						fij = wij / rijsq * SIX_TIMES_UNIT_FORCE;
						fxij = fij * rxij;
						fyij = fij * ryij;
						fxi += fxij;
						fyi += fyij;
						atom[j].fx -= fxij;
						atom[j].fy -= fyij;

					}

					if (hasCoulomb && interCoulomb) {
						if (Math.abs(atom[i].charge) > ZERO && Math.abs(atom[j].charge) > ZERO) {
							coul = atom[i].charge * atom[j].charge / Math.sqrt(rijsq) * rCD;
							vsum += coul;
							fij = coul / rijsq * GF_CONVERSION_CONSTANT;
							fxij = fij * rxij;
							fyij = fij * ryij;
							fxi += fxij;
							fyi += fyij;
							atom[j].fx -= fxij;
							atom[j].fy -= fyij;
						}
					}

				}

				atom[i].fx = fxi;
				atom[i].fy = fyi;

			}

			if (numberOfAtoms > 0)
				pointer[numberOfAtoms - 1] = nlist;

		} else {

			for (int i = 0, imax1 = numberOfAtoms - 1; i < imax1; i++) {

				rxi = atom[i].rx;
				ryi = atom[i].ry;
				fxi = atom[i].fx;
				fyi = atom[i].fy;

				jbeg = pointer[i];
				jend = pointer[i + 1];

				if (jbeg < jend) {

					for (int jnab = jbeg; jnab < jend; jnab++) {

						int j = neighborList[jnab];
						if (!ljBetweenBondPairs && bondTable[i][j])
							continue;
						rxij = rxi - atom[j].rx;
						ryij = ryi - atom[j].ry;
						minimumImageConvention();
						rijsq = rxij * rxij + ryij * ryij;

						if (rijsq < cutOffSquareMatrix[atom[i].id][atom[j].id]) {

							sr2 = sigab[i][j] / rijsq;
							/* check if this pair gets too close */
							if (sr2 > 2.0) {
								sr2 = 2.0;
								rijsq = 0.5 * sigab[i][j];
							}
							sr6 = sr2 * sr2 * sr2;
							sr12 = sr6 * sr6;

							if (isRepulsive(atom[i].id, atom[j].id)) {
								vij = sr6 * epsab[i][j] * crossRepulsionIntensity;
								wij = vij;
							} else {
								vij = (sr12 - sr6) * epsab[i][j];
								wij = vij + sr12 * epsab[i][j];
							}
							if (cutOffShift) {
								vij -= poten_LJ[atom[i].id][atom[j].id];
								wij -= slope_LJ[atom[i].id][atom[j].id];
							}

							vsum += vij;

							fij = wij / rijsq * SIX_TIMES_UNIT_FORCE;
							fxij = fij * rxij;
							fyij = fij * ryij;
							fxi += fxij;
							fyi += fyij;
							atom[j].fx -= fxij;
							atom[j].fy -= fyij;

						}
					}
				}

				if (hasCoulomb && interCoulomb) { // no Verlet neighbor list for long-range electrostatic potentials

					for (int j = i + 1; j < numberOfAtoms; j++) {
						if (!ljBetweenBondPairs && bondTable[i][j])
							continue;
						if (Math.abs(atom[i].charge) > ZERO && Math.abs(atom[j].charge) > ZERO) {
							rxij = rxi - atom[j].rx;
							ryij = ryi - atom[j].ry;
							minimumImageConvention();
							rijsq = rxij * rxij + ryij * ryij;
							coul = atom[i].charge * atom[j].charge / Math.sqrt(rijsq) * rCD;
							vsum += coul;
							fij = coul / rijsq * GF_CONVERSION_CONSTANT;
							fxij = fij * rxij;
							fyij = fij * ryij;
							fxi += fxij;
							fyi += fyij;
							atom[j].fx -= fxij;
							atom[j].fy -= fyij;
						}
					}

				}

				atom[i].fx = fxi;
				atom[i].fy = fyi;

			}

		}

		double inverseMass;
		for (int i = 0; i < numberOfAtoms; i++) {
			atom[i].fx += atom[i].hx * GF_CONVERSION_CONSTANT;
			atom[i].fy += atom[i].hy * GF_CONVERSION_CONSTANT;
			inverseMass = 1.0 / atom[i].mass;
			atom[i].fx *= inverseMass;
			atom[i].fy *= inverseMass;
		}

		// pointwise space restraints do not contribute to the internal pressure
		// but they maintains the energy conservation law.
		for (int i = 0; i < numberOfAtoms; i++) {
			double etemp;
			Atom a = atom[i];
			if (a.restraint != null) {
				a.restraint.dyn(a);
				etemp = a.restraint.getEnergy(a);
				vsum += etemp;
			}
			if (a.friction > ZERO) {
				inverseMass = GF_CONVERSION_CONSTANT * a.friction / a.mass;
				switch (a.dampType) {
				case AIR_DRAG:
					a.fx -= inverseMass * a.vx * universe.getViscosity();
					a.fy -= inverseMass * a.vy * universe.getViscosity();
					break;
				case DRY_FRICTION:
					a.fx -= inverseMass * Math.signum(a.vx) * DRY_FRICTION_CONVERTOR;
					a.fy -= inverseMass * Math.signum(a.vy) * DRY_FRICTION_CONVERTOR;
					break;
				}
			}
			if (a.getUserField() != null) {
				a.getUserField().dyn(a);
			}
		}

		double etemp = 0.0;
		synchronized (fields) {
			for (VectorField f : fields) {
				if (f instanceof GravitationalField) {
					GravitationalField gf = (GravitationalField) f;
					if (obstacles != null && !obstacles.isEmpty()) {
						RectangularObstacle obs = null;
						synchronized (obstacles.getSynchronizationLock()) {
							for (int jobs = 0, nobs = obstacles.size(); jobs < nobs; jobs++) {
								obs = obstacles.get(jobs);
								gf.dyn(obs);
								etemp = gf.getPotential(obs, time);
								vsum += etemp;
							}
						}
					}
					for (int i = 0; i < numberOfAtoms; i++) {
						gf.dyn(atom[i]);
						etemp = gf.getPotential(atom[i], time);
						vsum += etemp;
					}
				} else if (f instanceof ElectricField) {
					ElectricField ef = (ElectricField) f;
					for (int i = 0; i < numberOfAtoms; i++) {
						if (Math.abs(atom[i].charge) > ZERO) {
							ef.dyn(universe.getDielectricConstant(), atom[i], time);
							etemp = ef.getPotential(atom[i], time);
							vsum += etemp;
						}
					}
				} else if (f instanceof MagneticField) {
					MagneticField mf = (MagneticField) f;
					for (int i = 0; i < numberOfAtoms; i++) {
						if (Math.abs(atom[i].charge) > ZERO)
							mf.dyn(atom[i]);
					}
				} else if (f instanceof AccelerationalField) {
					AccelerationalField af = (AccelerationalField) f;
					if (obstacles != null && !obstacles.isEmpty()) {
						RectangularObstacle obs = null;
						synchronized (obstacles.getSynchronizationLock()) {
							for (int jobs = 0, nobs = obstacles.size(); jobs < nobs; jobs++) {
								obs = obstacles.get(jobs);
								af.dyn(obs);
								etemp = af.getPotential(obs, time);
								vsum += etemp;
							}
						}
					}
					for (int i = 0; i < numberOfAtoms; i++) {
						af.dyn(atom[i]);
						etemp = af.getPotential(atom[i], time);
						vsum += etemp;
					}
				}
			}
		}

		FieldArea[] fa = view.getFieldAreas();
		if (fa != null) {
			for (FieldArea x : fa) {
				for (int i = 0; i < numberOfAtoms; i++)
					x.interact(atom[i]);
				if (photons != null && !photons.isEmpty()) {
					List<Photon> absorbedPhotons = new ArrayList<Photon>();
					synchronized (photons) {
						for (Photon p : photons)
							if (x.absorb(p))
								absorbedPhotons.add(p);
					}
					if (!absorbedPhotons.isEmpty())
						photons.removeAll(absorbedPhotons);
				}
				if (freeElectrons != null && !freeElectrons.isEmpty()) {
					List<Electron> absorbedElectrons = new ArrayList<Electron>();
					synchronized (freeElectrons) {
						for (Electron e : freeElectrons)
							if (x.absorb(e))
								absorbedElectrons.add(e);
					}
					if (!absorbedElectrons.isEmpty())
						freeElectrons.removeAll(absorbedElectrons);
				}
			}
		}

		if (time < 0) {
			for (int i = 0; i < numberOfAtoms; i++) {
				atom[i].fx *= atom[i].mass;
				atom[i].fy *= atom[i].mass;
			}
		}

		return numberOfAtoms > 0 ? vsum / numberOfAtoms : vsum;

	}