for (i=0; i<Nv; i++) {
	
	// Do this for X direction
	E_nxt_x = 0; E_prv_x = 0; E_nxt_y = 0; E_prv_y = 0; E_nxt_z = 0; E_prv_z = 0;
	
	// Loop over all those cells that are affected by the vertex i
	
	for (j=0; j<4; j++) {
		Cell_sel = V_cells[i][j];
		
		if (Cell_sel < 0) continue;
		
		// Assing the parameters of this cell
		
		K_par[0] = A_K_cell[Cell_sel]; K_par[1] = B_K_cell[Cell_sel]; K_par[2] = V_K_cell[Cell_sel];
		K_par[3] = L_K_cell[Cell_sel]*C_side_factor[Cell_sel][0]; K_par[4] = L_K_cell[Cell_sel]*C_side_factor[Cell_sel][1];
		K_par[5] = L_K_cell[Cell_sel]*C_side_factor[Cell_sel][2]; K_par[6] = L_K_cell[Cell_sel]*C_side_factor[Cell_sel][3];
		K_par[7] = D_K_cell[Cell_sel];
		
		O_par[0] = A_0_cell[Cell_sel]; O_par[1] = B_0_cell[Cell_sel]; O_par[2] = V_0_cell[Cell_sel];
		O_par[3] = L_0_cell[Cell_sel]; O_par[4] = L_0_cell[Cell_sel]; O_par[5] = L_0_cell[Cell_sel];
		O_par[6] = L_0_cell[Cell_sel]; O_par[7] = D_0_cell[Cell_sel];
		
		
		//**********************************************
		// Do for X
		//cout << "For cell " << Cell_sel << " K0 = [" << Par[4][0] << ", " << Par[4][1] << ", " << Par[5][0] << ", " << Par[5][1]<< "]" << endl;
		// Perturb the vertex in the X direction
		Vx_nxt[i] += dx; Vx_prv[i] -= dx;
		
		// Assing apical and basal verticies for the +X pertubation
		for (k=0; k<4; k++) {
			AP[k][0] = Vx_nxt[C_apver[Cell_sel][k]]; 
			AP[k][1] = Vy_nxt[C_apver[Cell_sel][k]]; 
			AP[k][2] = Vz_nxt[C_apver[Cell_sel][k]];
			
			BS[k][0] = Vx_nxt[C_bsver[Cell_sel][k]]; 
			BS[k][1] = Vy_nxt[C_bsver[Cell_sel][k]]; 
			BS[k][2] = Vz_nxt[C_bsver[Cell_sel][k]];
		}
		//E_nxt_x += E_from_cell(AP, BS, Par);
		//E_nxt_x += E_from_cell_circ(AP, BS, Par);
		//E_nxt_x += E_from_cell_center(AP, BS, Par);
		E_nxt_x += E_from_cell_V2_diag(AP, BS, K_par, O_par);
		
		// Assing apical and basal verticies for the -X pertubation
		for (k=0; k<4; k++) {
			AP[k][0] = Vx_prv[C_apver[Cell_sel][k]]; 
			AP[k][1] = Vy_prv[C_apver[Cell_sel][k]]; 
			AP[k][2] = Vz_prv[C_apver[Cell_sel][k]];
			
			BS[k][0] = Vx_prv[C_bsver[Cell_sel][k]]; 
			BS[k][1] = Vy_prv[C_bsver[Cell_sel][k]]; 
			BS[k][2] = Vz_prv[C_bsver[Cell_sel][k]];
		}
		//E_prv_x += E_from_cell(AP, BS, Par);
		//E_prv_x += E_from_cell_circ(AP, BS, Par);
		//E_prv_x += E_from_cell_center(AP, BS, Par);
		E_prv_x += E_from_cell_V2_diag(AP, BS, K_par, O_par);
		
		// Remove the pertubation
		Vx_nxt[i] = Vx[i]; Vx_prv[i] = Vx[i];
		
		//**********************************************
		// Repeat for  for Y
		Vy_nxt[i] += dx; Vy_prv[i] -= dx;
		
		// Assing apical and basal verticies
		for (k=0; k<4; k++) {
			AP[k][0] = Vx_nxt[C_apver[Cell_sel][k]];
			AP[k][1] = Vy_nxt[C_apver[Cell_sel][k]];
			AP[k][2] = Vz_nxt[C_apver[Cell_sel][k]];
			
			BS[k][0] = Vx_nxt[C_bsver[Cell_sel][k]]; 
			BS[k][1] = Vy_nxt[C_bsver[Cell_sel][k]]; 
			BS[k][2] = Vz_nxt[C_bsver[Cell_sel][k]];
		}
		//E_nxt_y += E_from_cell(AP, BS, Par);
		//E_nxt_y += E_from_cell_circ(AP, BS, Par);
		//E_nxt_y += E_from_cell_center(AP, BS, Par);
		E_nxt_y += E_from_cell_V2_diag(AP, BS, K_par, O_par);
		
		// Assing apical and basal verticies
		for (k=0; k<4; k++) {
			AP[k][0] = Vx_prv[C_apver[Cell_sel][k]]; 
			AP[k][1] = Vy_prv[C_apver[Cell_sel][k]]; 
			AP[k][2] = Vz_prv[C_apver[Cell_sel][k]];
			
			BS[k][0] = Vx_prv[C_bsver[Cell_sel][k]]; 
			BS[k][1] = Vy_prv[C_bsver[Cell_sel][k]]; 
			BS[k][2] = Vz_prv[C_bsver[Cell_sel][k]];
		}				
		//E_prv_y += E_from_cell(AP, BS, Par);
		//E_prv_y += E_from_cell_circ(AP, BS, Par);
		//E_prv_y += E_from_cell_center(AP, BS, Par);
		E_prv_y += E_from_cell_V2_diag(AP, BS, K_par, O_par);
		
		Vy_nxt[i] = Vy[i]; Vy_prv[i] = Vy[i];
		
		//**********************************************
		// Repeat for Z
		Vz_nxt[i] += dx; Vz_prv[i] -= dx;
		
		// Assing apical and basal verticies
		for (k=0; k<4; k++) {
			AP[k][0] = Vx_nxt[C_apver[Cell_sel][k]]; 
			AP[k][1] = Vy_nxt[C_apver[Cell_sel][k]]; 
			AP[k][2] = Vz_nxt[C_apver[Cell_sel][k]];
			
			BS[k][0] = Vx_nxt[C_bsver[Cell_sel][k]]; 
			BS[k][1] = Vy_nxt[C_bsver[Cell_sel][k]]; 
			BS[k][2] = Vz_nxt[C_bsver[Cell_sel][k]];
		}
		//E_nxt_z += E_from_cell(AP, BS, Par);
		//E_nxt_z += E_from_cell_circ(AP, BS, Par);
		//E_nxt_z += E_from_cell_center(AP, BS, Par);
		E_nxt_z += E_from_cell_V2_diag(AP, BS, K_par, O_par);
		
		// Assing apical and basal verticies
		for (k=0; k<4; k++) {
			AP[k][0] = Vx_prv[C_apver[Cell_sel][k]]; 
			AP[k][1] = Vy_prv[C_apver[Cell_sel][k]];
			AP[k][2] = Vz_prv[C_apver[Cell_sel][k]];
			
			BS[k][0] = Vx_prv[C_bsver[Cell_sel][k]]; 
			BS[k][1] = Vy_prv[C_bsver[Cell_sel][k]]; 
			BS[k][2] = Vz_prv[C_bsver[Cell_sel][k]];
		}
		//E_prv_z += E_from_cell(AP, BS, Par);
		//E_prv_z += E_from_cell_circ(AP, BS, Par);
		//E_prv_z += E_from_cell_center(AP, BS, Par);
		E_prv_z += E_from_cell_V2_diag(AP, BS, K_par, O_par);
		
		Vz_nxt[i] = Vz[i]; Vz_prv[i] = Vz[i];
		//**********************************************
	}
	
	//**********************************************
	// Now if it is a basal vertex also take into account the energy contribution from the yolk
	
	if (yolk_present == 1 && i < Nv/2){
		
		// For X 
		Yolk_V_tmp = Yolk_V;
		for (j=0; j<4; j++) {
			Cell_sel = V_cells[i][j];
			
			if (Cell_sel < 0) continue;
			
			// Perturb the vertex in the X direction
			Vx_nxt[i] += dx;
			
			// Assing apical and basal verticies for the +X pertubation
			for (k=0; k<4; k++) {
				AP[k][0] = Vx_nxt[C_apver[Cell_sel][k]]; 
				AP[k][1] = Vy_nxt[C_apver[Cell_sel][k]]; 
				AP[k][2] = Vz_nxt[C_apver[Cell_sel][k]];
				
				BS[k][0] = r_out*AP[k][0]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]); 
				BS[k][1] = AP[k][1]; 
				BS[k][2] = r_out*AP[k][2]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]);
			}
			
			Yolk_V_tmp += Volume_Calc(AP, BS) - Yolk_volumes[Cell_sel];
			Vx_nxt[i] = Vx[i];
		}
		
		E_nxt_x += K_yolk*pow(Yolk_V_tmp - Yolk_O, 2);
		
		Yolk_V_tmp = Yolk_V;
		for (j=0; j<4; j++) {
			Cell_sel = V_cells[i][j];
			
			if (Cell_sel < 0) continue;
			
			// Perturb the vertex in the X direction
			Vx_prv[i] -= dx;
			
			// Assing apical and basal verticies for the +X pertubation
			for (k=0; k<4; k++) {
				AP[k][0] = Vx_prv[C_apver[Cell_sel][k]]; 
				AP[k][1] = Vy_prv[C_apver[Cell_sel][k]]; 
				AP[k][2] = Vz_prv[C_apver[Cell_sel][k]];
				
				BS[k][0] = r_out*AP[k][0]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]); 
				BS[k][1] = AP[k][1]; 
				BS[k][2] = r_out*AP[k][2]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]);
			}
			
			Yolk_V_tmp += Volume_Calc(AP, BS) - Yolk_volumes[Cell_sel];
			Vx_prv[i] = Vx[i];
		}
		
		E_prv_x += K_yolk*pow(Yolk_V_tmp - Yolk_O, 2);
		
		// For Y 
		Yolk_V_tmp = Yolk_V;
		for (j=0; j<4; j++) {
			Cell_sel = V_cells[i][j];
			
			if (Cell_sel < 0) continue;
			
			// Perturb the vertex in the X direction
			Vy_nxt[i] += dx;
			
			// Assing apical and basal verticies for the +X pertubation
			for (k=0; k<4; k++) {
				AP[k][0] = Vx_nxt[C_apver[Cell_sel][k]]; 
				AP[k][1] = Vy_nxt[C_apver[Cell_sel][k]]; 
				AP[k][2] = Vz_nxt[C_apver[Cell_sel][k]];
				
				BS[k][0] = r_out*AP[k][0]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]); 
				BS[k][1] = AP[k][1]; 
				BS[k][2] = r_out*AP[k][2]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]);
			}
			
			Yolk_V_tmp += Volume_Calc(AP, BS) - Yolk_volumes[Cell_sel];
			Vy_nxt[i] = Vy[i];
		}
		
		E_nxt_y += K_yolk*pow(Yolk_V_tmp - Yolk_O, 2);
		
		Yolk_V_tmp = Yolk_V;
		for (j=0; j<4; j++) {
			Cell_sel = V_cells[i][j];
			
			if (Cell_sel < 0) continue;
			
			// Perturb the vertex in the X direction
			Vy_prv[i] -= dx;
			
			// Assing apical and basal verticies for the +X pertubation
			for (k=0; k<4; k++) {
				AP[k][0] = Vx_prv[C_apver[Cell_sel][k]]; 
				AP[k][1] = Vy_prv[C_apver[Cell_sel][k]]; 
				AP[k][2] = Vz_prv[C_apver[Cell_sel][k]];
				
				BS[k][0] = r_out*AP[k][0]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]); 
				BS[k][1] = AP[k][1]; 
				BS[k][2] = r_out*AP[k][2]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]);
			}
			
			Yolk_V_tmp += Volume_Calc(AP, BS) - Yolk_volumes[Cell_sel];
			Vy_prv[i] = Vy[i];
		}
		
		E_prv_y += K_yolk*pow(Yolk_V_tmp - Yolk_O, 2);
		
		// For Z 
		Yolk_V_tmp = Yolk_V;
		for (j=0; j<4; j++) {
			Cell_sel = V_cells[i][j];
			
			if (Cell_sel < 0) continue;
			
			// Perturb the vertex in the X direction
			Vz_nxt[i] += dx;
			
			// Assing apical and basal verticies for the +X pertubation
			for (k=0; k<4; k++) {
				AP[k][0] = Vx_nxt[C_apver[Cell_sel][k]]; 
				AP[k][1] = Vy_nxt[C_apver[Cell_sel][k]]; 
				AP[k][2] = Vz_nxt[C_apver[Cell_sel][k]];
				
				BS[k][0] = r_out*AP[k][0]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]); 
				BS[k][1] = AP[k][1]; 
				BS[k][2] = r_out*AP[k][2]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]);
			}
			
			Yolk_V_tmp += Volume_Calc(AP, BS) - Yolk_volumes[Cell_sel];
			Vz_nxt[i] = Vz[i];
		}
		
		E_nxt_z += K_yolk*pow(Yolk_V_tmp - Yolk_O, 2);
		
		Yolk_V_tmp = Yolk_V;
		for (j=0; j<4; j++) {
			Cell_sel = V_cells[i][j];
			
			if (Cell_sel < 0) continue;
			
			// Perturb the vertex in the X direction
			Vz_prv[i] -= dx;
			
			// Assing apical and basal verticies for the +X pertubation
			for (k=0; k<4; k++) {
				AP[k][0] = Vx_prv[C_apver[Cell_sel][k]]; 
				AP[k][1] = Vy_prv[C_apver[Cell_sel][k]]; 
				AP[k][2] = Vz_prv[C_apver[Cell_sel][k]];
				
				BS[k][0] = r_out*AP[k][0]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]); 
				BS[k][1] = AP[k][1]; 
				BS[k][2] = r_out*AP[k][2]/sqrt(AP[k][0]*AP[k][0] + AP[k][2]*AP[k][2]);
			}
			
			Yolk_V_tmp += Volume_Calc(AP, BS) - Yolk_volumes[Cell_sel];
			Vz_prv[i] = Vz[i];
		}
		
		E_prv_z += K_yolk*pow(Yolk_V_tmp - Yolk_O, 2);
		
	}
	
	//**********************************************
	// Now take into account the energy penalty due to the plane constrains
	
	if (plane_constraint_present == 1){
		//**********************************************
		// Do for X
		
		for (j=0; j<24; j++) {
			
			if ((Plane_test_p1[i][j] < 0) || (Plane_test_p2[i][j] < 0) || (Plane_test_p3[i][j] < 0))
				continue;
			
			Point1[0] = Vx[Plane_test_p1[i][j]];Point1[1] = Vy[Plane_test_p1[i][j]];Point1[2] = Vz[Plane_test_p1[i][j]];
			Point2[0] = Vx[Plane_test_p2[i][j]];Point2[1] = Vy[Plane_test_p2[i][j]];Point2[2] = Vz[Plane_test_p2[i][j]];
			Point3[0] = Vx[Plane_test_p3[i][j]];Point3[1] = Vy[Plane_test_p3[i][j]];Point3[2] = Vz[Plane_test_p3[i][j]];
			
			// Pertubation for x
			Point0[0] = Vx[i] + dx;Point0[1] = Vy[i];Point0[2] = Vz[i];				
			Dis_tmp = Dis_to_plane(Point1, Point2, Point3, Point0);
			if (Dis_tmp < R_max) {
				E_nxt_x += Plane_potential[int(Dis_tmp/dx_plane)];
			}
			
			Point0[0] = Vx[i] - dx;Point0[1] = Vy[i];Point0[2] = Vz[i];				
			Dis_tmp = Dis_to_plane(Point1, Point2, Point3, Point0);
			if (Dis_tmp < R_max) {
				E_prv_x += Plane_potential[int(Dis_tmp/dx_plane)];
			}
			
			// Pertubation for y
			Point0[0] = Vx[i];Point0[1] = Vy[i] + dx;Point0[2] = Vz[i];
			Dis_tmp = Dis_to_plane(Point1, Point2, Point3, Point0);
			if (Dis_tmp < R_max) {
				E_nxt_y += Plane_potential[int(Dis_tmp/dx_plane)];
			}
			
			Point0[0] = Vx[i];Point0[1] = Vy[i] - dx;Point0[2] = Vz[i];
			Dis_tmp = Dis_to_plane(Point1, Point2, Point3, Point0);
			if (Dis_tmp < R_max) {
				E_prv_y += Plane_potential[int(Dis_tmp/dx_plane)];
			}
			
			// Pertubation for z
			Point0[0] = Vx[i];Point0[1] = Vy[i];Point0[2] = Vz[i] + dx;
			Dis_tmp = Dis_to_plane(Point1, Point2, Point3, Point0);
			if (Dis_tmp < R_max) {
				E_nxt_z += Plane_potential[int(Dis_tmp/dx_plane)];
			}
			
			Point0[0] = Vx[i];Point0[1] = Vy[i];Point0[2] = Vz[i] - dx;
			Dis_tmp = Dis_to_plane(Point1, Point2, Point3, Point0);
			if (Dis_tmp < R_max) {
				E_prv_z += Plane_potential[int(Dis_tmp/dx_plane)];
			}
			
		}
		
		
	}
	
	//**********************************************
	
	V_grad[i][0] = -(E_nxt_x - E_prv_x);
	V_grad[i][1] = -(E_nxt_y - E_prv_y);
	V_grad[i][2] = -(E_nxt_z - E_prv_z);
	
}
