 	float rho = 0;
 	float scorr = 0;
 
-	FOR_EACH_NEIGHBOUR
+	vec3 grad_pi_Ci = vec3 (0, 0, 0);
+	FOR_EACH_NEIGHBOUR(j)
 	{			
 		// compute rho_i (equation 2)
-		float len = length (particle.position - particles[neighbourid].position);
+		float len = length (particle.position - particles[j].position);
 		float tmp = Wpoly6 (len);
 		rho += tmp;
 	

	
	vec3 deltap = vec3 (0, 0, 0);
	
	FOR_EACH_NEIGHBOUR
	{
		// accumulate position corrections (part of equation 12)
		deltap += (lambda + lambdas[neighbourid]) * gradWspiky (particle.position - particles[neighbourid].position);
	}
	END_FOR_EACH_NEIGHBOUR


		scorr += -0.1 * tmp * tmp;
			
		// sum gradients of Ci (equation 8 and parts of equation 9)
		// use j as k so that we can stay in the same loop
		uint k = j;
		// for k=i it holds that particle.position==particles[k].position, so
		// the summand does no harm
		vec3 grad_pk_Ci = vec3 (0, 0, 0);
		grad_pk_Ci = gradWspiky (particle.position - particles[k].position);
		grad_pk_Ci /= rho_0;
		sum_k_grad_Ci += dot (grad_pk_Ci, grad_pk_Ci);
		
		// now use j as j again and accumulate grad_pi_Ci for the case k=i
		// from equation 8
		grad_pi_Ci += gradWspiky (particle.position - particles[j].position); 
	}
	END_FOR_EACH_NEIGHBOUR
	// add grad_pi_Ci to the sum
	sum_k_grad_Ci += dot (grad_pi_Ci, grad_pi_Ci) / rho_0;
	
	// compute lambda_i (equations 1 and 9)
	float C_i = rho / rho_0 - 1;
