(*Annotated Biplane vs. Octa example. Allan Stewart
FCC lattice models for protein folding (Thesis 2010)*)

(*Get FCC *)
b=Normal@LatticeData["FaceCenteredCubic","Basis"]
(*One layer of biplane*)
b1 =Flatten[Table[i*b[[1]]+j*b[[2]],{i,0,2},{j,0,2}],1]
(*3-unit triangle for octa*)
p3 = {{1,1,0}-{3,0,0},{1,1,0}+-1*b[[1]]-{3,0,0},{1,1,0}+b[[2]]-{3,0,0}}
(*6-unit triangle for octa*)
p6=Join[Table[i*b[[2]],{i,-1,1}],Table[b[[1]]+i*b[[2]],{i,0,1}], {2*b[[1]]+b[[2]]}]
(*Reflection transformation function *)
rt=ReflectionTransform[{2,0}];
(*Biplane equals two "b1" planes*)
proBiplane=Join[b1,Map[(#+b[[3]])&,b1]]
(*I've modded the biplane to be 20 AAs: 3x3x2 + 2 -- for optimality purposes*)
proBiplane={{0,0,0},{1,-1,0},{2,-2,0},{-1,-1,0},{0,-2,0},{1,-3,0},{-2,-2,0},{-1,-3,0},{0,-4,0},{0,1,-1},{1,0,-1},{2,-1,-1},{-1,0,-1},{0,-1,-1},{1,-2,-1},{-2,-1,-1},{-1,-2,-1},{0,-3,-1},{1,0,0},{1,0,0}+b[[3]]};
(*Octahedron is constructed by layering 1,p3, p6, p6, p3, 1*)
proOcta=Join[
(*top layer  3 unit*)Map[(Join[Flatten[rt[{#[[1]],#[[2]]}]],{0}]+{0,-2,0})&,p3],
(*mid-bottom layer  6 unit*)
Map[(Join[Flatten[rt[{#[[1]],#[[2]]}]],{0}]+b[[3]]+{1,-1,-1})&,p6],(*mid-top layer  6 unit *)Map[(#+2*b[[3]]+{2,-2,1})&,p6],
(*bottom layer 3 unit *)
Map[(Join[Flatten[rt[{#[[1]],#[[2]]}]],{0}]+3*b[[3]]+{0,-6,0})&,p3],
(*add singletons in most optimal manner*)
{{3,0,0},{0,-2,-3}}]
(*END "octahedron*)

(*Plot function*)
{Graphics3D[Map[Sphere[#,.2]&,proBiplane],Boxed->True,Axes->True],Graphics3D[Map[Sphere[#,.2]&,proOcta],Boxed->True,Axes->True]}
(*Contact counting function*)
ContactsH[p_,rv_]:=(
contacts:=0;
For[i=1,i≤Length[p],i++,For[j=i+1,j≤Length[p],j++,If[And[rv[[i]]==0,rv[[j]]==0,Total[Abs[Flatten[p[[i]]-p[[j]]]]]≤2],contacts++;]]]
Return[contacts];
)

(*Count contacts for biplane model*)
rv=ConstantArray[0,Length[proBiplane]];
biplaneContacts=0; p=proBiplane;
For[i=1,i≤Length[p],i++,For[j=i+1,j≤Length[p],j++,If[And[rv[[i]]==0,rv[[j]]==0,Total[Abs[Flatten[p[[i]]-p[[j]]]]]≤2],biplaneContacts++;]]];
(*Count contacts for octahedral model*)
rv=ConstantArray[0,Length[proOcta]];
octaContacts=0; p=proOcta;
For[i=1,i≤Length[p],i++,For[j=i+1,j≤Length[p],j++,If[And[rv[[i]]==0,rv[[j]]==0,Total[Abs[Flatten[p[[i]]-p[[j]]]]]≤2],octaContacts++;]]];
(*Print contacts number*)
biplaneContacts
octaContacts