x = Transpose[{{x1,y1,z1,x2,y2,z2,x3,y3,z3}}];

   G = {	{ 4, 0, 0,-2, 0, 0,-2, 0, 0},
		{ 0, 4, 0, 0,-2, 0, 0,-2, 0},
		{ 0, 0, 4, 0, 0,-2, 0, 0,-2},
		{-2, 0, 0, 4, 0, 0,-2, 0, 0},
		{ 0,-2, 0, 0, 4, 0, 0,-2, 0},
		{ 0, 0,-2, 0, 0, 4, 0, 0,-2},
		{-2, 0, 0,-2, 0, 0, 4, 0, 0},
		{ 0,-2, 0, 0,-2, 0, 0, 4, 0},
		{ 0, 0,-2, 0, 0,-2, 0, 0, 4}
	};

p = 2 ( (x1-x2)^2 + (x1-x3)^2 + (x2-x3)^2 
	+	(y1-y2)^2 + (y1-y3)^2 + (y2-y3)^2
	+	(z1-z2)^2 + (z1-z3)^2 + (z2-z3)^2 );

f = Function[x, If[Mod[x,3]==0,1,0]  ]
n=3;
MatrixForm[ G = 
		2 n IdentityMatrix[3 n]
	-	2   Table[ f[i-j], {i,3 n},{j, 3 n}] 
	]