How to get rid of the perspective effect in a 3D graphics

Maybe this will help a little (adapting documentation exaple for Slider2D):

 DynamicModule[{p = {2 π, 0}},
   Row @ {Slider2D[Dynamic[p], {{2 Pi, 0}, {0, Pi}}],
          Plot3D[Exp[-(x^2 + y^2)], {x, -3, 3}, {y, -3, 3}, 
            ImageSize -> {700, 700}, PlotRange -> All, ViewAngle -> .0015,
            ViewPoint -> Dynamic[1200 {Cos[p[[1]]] Sin[p[[2]]], 
                                       Sin[p[[1]]] Sin[p[[2]]], 
                                       Cos[p[[2]]]}
                         ]
          ]   
   }]

Notice that for parallel projection (like here) box edges do not seem to be parallel but they are and for default Mathematica display they are not but they look parallel

I think the only way to do this is by dynamically reseting the ViewMatrix to be an orthographic projection. It was beyond my ability, patience, or inclination to figure how to decompose the ViewMatrix that is created when the graphic is moved into the components ViewPoint, ViewVertical, etc. It seemed to me that the front end usually make a discontinuous jump from the view point in the orthographic projection to the initial view point in rotated graphic. Discouraged by this apparent behavior, I opted for a hybrid solution.

I used an ordinary Graphics3D with Dynamic view properties and used Inset to insert the graphics tom2 into it. The view properties of the outer Graphics3D are used to compute a corresponding orthographic projection ViewMatrix to be used to display the graphics tom2. So the mouse rotates the outer graphics and the code rotates the inset tom2 in the same way. The only difference is that tom2 is projected orthographically instead of perspectively.

To do this, I used and adapted code from Heike in this answer and Alexey Popkov in this answer. I needed Alexey Popkov's completePlotRange because AbsoluteOptions would not return the actual plot range.

(* Heike *)
theta[v1_] := ArcTan[v1[[3]], Norm[v1[[;; 2]]]];
phi[v1_] := If[Norm[v1[[;; 2]]] > .0001, ArcTan[v1[[1]], v1[[2]]], 0];
alpha[vert_, v1_] := ArcTan[{-Sin[phi[v1]], Cos[phi[v1]], 0}.vert, 
                             Cross[v1/Norm[v1], {-Sin[phi[v1]], Cos[phi[v1]], 0}].vert];
tt[v1_, vert_, center_, r_, scale_] := TransformationMatrix[
   RotationTransform[-alpha[vert/scale, v1], {0, 0, 1}] .
   RotationTransform[-theta[v1], {0, 1, 0}] .
   RotationTransform[-phi[v1], {0, 0, 1}] .
   ScalingTransform[r {1, 1, 1}] .
   TranslationTransform[-center]];

(* orthographic projection *)
pp = N@{{1, 0, 0, 1}, {0, 1, 0, 1}, {0, 0, -1, 0}, {0, 0, 0, 2}};

(* Alexey Popkov *)
completePlotRange[plot : (_Graphics | _Graphics3D)] := 
 Quiet@Last@
   Last@Reap[
     Rasterize[
      Show[plot, Axes -> True, Ticks -> (Sow[{##}] &), 
       DisplayFunction -> Identity], ImageResolution -> 1]]

(* OP *)
a0 = 1; ceng = 1; znn = 3*ceng; ynn = 5; xnn = 5; a1 = 2; cznn = 3; hh = 2;
p1[θ_] := RotationTransform[θ, {0, 0, 1}];
posx1 = (nx1 - 1) a0 + ((ny1 - 1) a0)/2 + (Mod[nz1 - 1, 3] a0)/2;
posy1 = Sqrt[3]/2 (ny1 - 1) a0 + (Sqrt[3] Mod[nz1 - 1, 3] a0)/6;
posz1 = Sqrt[3]/2 (nz1 - 1) a0 + Floor[(nz1 - 1)/3] a0;
sinpq2 = Transpose[Table[{posx1, posy1, posz1},
                         {ny1, ynn}, {nz1, znn}, {nx1, -xnn, 0}]];
rr = 0.3 a0;
(* updated - switched Graphics3D and Table, adjusted options *)
tom2 = Graphics3D[
   Table[{If[Mod[cc, 3] == 2, Pink, Yellow], Sphere[N@sinpq2[[cc, bb, aa]], rr]},
         {aa, xnn}, {bb, ynn}, {cc, znn}],
   Axes -> False, BoxStyle -> Directive[Orange, Opacity[1]], 
   BoxRatios -> Automatic, Background -> Black];

(* adapted from Heike *)
bb = completePlotRange@tom2;
scale = 1/Abs[#1 - #2] & @@@ bb;
center = Mean /@ bb;
vv = {Flatten[Differences /@ bb] + center, center};
v1 = (vv[[1]] - center);
vert = {0, 0, 1} - {0, 0, 1}.v1 v1;
viewAngle = 2 ArcCot[2.];

(* final graphic *)
Graphics3D[{},
 Epilog -> Inset[
   Show[tom2, 
        ViewMatrix ->
          Dynamic@{tt[v1, vert, center, Cot[viewAngle/2]/Norm[v1], scale], pp},
        PlotRange -> bb],
   Center, Center, 1],
 ViewAngle    -> Dynamic[viewAngle], 
 ViewVector   -> Dynamic[vv, (vv = #; center = vv[[2]]; v1 = vv[[1]] - center) &], 
 ViewVertical -> Dynamic[vert],
 SphericalRegion -> True, PlotRange -> bb, Boxed -> False]

Output after manual rotation:

I just found the simplest solution.

Don't set ViewPoint to infinity!!! Just set it any number large enough, the effect will equivalent to setting to infinity.

Try for example the following

ListPointPlot3D[Tuples[Range[10], 3], BoxRatios -> Automatic, 
 ViewPoint -> {0, 0, 1000}]

for your example just change ViewPoint -> {0, 0, ∞} to ViewPoint -> {0, 0, 1000} or any number large than 1000 as you wish.

You'll be amazed at how simple it is to get rid of perspective. Rotate freely! : )