private float calculateTBonds() {

		int n = model.tBonds.size();
		if (n <= 0)
			return 0;
		TBond tBond;
		float energy = 0;

		synchronized (model.tBonds) {

			for (int i = 0; i < n; i++) {

				tBond = model.tBonds.get(i);
				atom1 = tBond.getAtom1();
				atom2 = tBond.getAtom2();
				atom3 = tBond.getAtom3();
				atom4 = tBond.getAtom4();
				if (!atom1.isMovable() && !atom2.isMovable() && !atom3.isMovable() && !atom4.isMovable())
					continue;
				angle = tBond.getAngle();
				strength = tBond.getStrength();
				rxij = atom1.rx - atom2.rx;
				ryij = atom1.ry - atom2.ry;
				rzij = atom1.rz - atom2.rz;
				vector1.set(rxij, ryij, rzij);
				rxlk = atom4.rx - atom3.rx;
				rylk = atom4.ry - atom3.ry;
				rzlk = atom4.rz - atom3.rz;
				vector2.set(rxlk, rylk, rzlk);
				rijsq = rxij * rxij + ryij * ryij + rzij * rzij;
				rlksq = rxlk * rxlk + rylk * rylk + rzlk * rzlk;
				rij = (float) Math.sqrt(rijsq);
				rlk = (float) Math.sqrt(rlksq);
				sr2 = rxij * rxlk + ryij * rylk + rzij * rzlk;
				theta = vector1.angle(vector2);
				sintheta = (float) Math.sin(theta);
				if (Math.abs(sintheta) < MIN_SINTHETA) {// zero or 180 degree disaster
					sintheta = sintheta > 0 ? MIN_SINTHETA : -MIN_SINTHETA;
				}
				sr12 = (float) Math.sin(tBond.getPeriodicity() * theta - tBond.getAngle());
				sr6 = 0.5f * tBond.getPeriodicity() * strength * sr12 / (rij * rlk * sintheta);

				rijsq = 1.0f / rijsq;
				rlksq = 1.0f / rlksq;
				fxi = sr6 * (rxlk - sr2 * rxij * rijsq);
				fyi = sr6 * (rylk - sr2 * ryij * rijsq);
				fzi = sr6 * (rzlk - sr2 * rzij * rijsq);
				fxl = sr6 * (rxij - sr2 * rxlk * rlksq);
				fyl = sr6 * (ryij - sr2 * rylk * rlksq);
				fzl = sr6 * (rzij - sr2 * rzlk * rlksq);

				inverseMass1 = GF_CONVERSION_CONSTANT / atom1.mass;
				inverseMass2 = GF_CONVERSION_CONSTANT / atom2.mass;
				inverseMass3 = GF_CONVERSION_CONSTANT / atom3.mass;
				inverseMass4 = GF_CONVERSION_CONSTANT / atom4.mass;
				atom1.fx += fxi * inverseMass1;
				atom1.fy += fyi * inverseMass1;
				atom1.fz += fzi * inverseMass1;
				atom4.fx += fxl * inverseMass4;
				atom4.fy += fyl * inverseMass4;
				atom4.fz += fzl * inverseMass4;
				atom2.fx -= fxi * inverseMass2;
				atom2.fy -= fyi * inverseMass2;
				atom2.fz -= fzi * inverseMass2;
				atom3.fx -= fxl * inverseMass3;
				atom3.fy -= fyl * inverseMass3;
				atom3.fz -= fzl * inverseMass3;

				// note that we use 1-cos(...) instead of 1+cos(...) as used on the following page:
				// http://en.wikipedia.org/wiki/AMBER
				// This reduced the equilibrium energy to zero, as in the case of radial and angular bonds
				energy += strength * (1.0f - Math.cos(tBond.getPeriodicity() * theta - tBond.getAngle()));

			}

		}

		return 0.5f * energy;

	}