Is there a better way to build a matrix with position coordinates conditions?

Question

Suppose that I have a list of coordinates {x,y} following some geometry in space, for example

TRIANGLE = Triangle[{{0, 0}, {1, 1}, {2, 0}}];
pts = RandomPoint[TRIANGLE, 3000];

And, I want to construct a matrix in which follow some rule

First we create the matrix with zeros (I would like to know if its really necessary to create it)

MATRIX = Table[0, {i, 1, Length[pts]}, {j, 1, Length[pts]}];

Or using Array

MATRIX = Array[0 &, {Length[pts], Length[pts]}];

And then, fill the matrix with the following condition (Is there a better way to construct this part?)

Table[
  If[
   EuclideanDistance[pts[[i]], pts[[j]]] == 
    EuclideanDistance[pts[[1]], pts[[2]]](*or equal to other condition*),

   MATRIX[[i, j]] = r
   ,
   Nothing
   ]
  , {i, 1, Length[pts]}, {j, 1, Length[pts]}
  ];

Edit1: (Timing check for 1000x1000 matrices)

The method above gives: 4.26563

And

MAT = SparseArray[{{i_, j_} /; 
       EuclideanDistance[pts[[i]], pts[[j]]] == 
        EuclideanDistance[pts[[1]], pts[[2]]] -> r}, {Length[pts], 
     Length[pts]}]; 

Gives: 4.40625

I know this example is not so good because pts follows a random distribution, but I'm using it because I want to create big matrices like that.

I would like to know if there is another way to do it, in view of, for larger matrices it takes a lot of computational time.

Edit2: For a realistic geometrical system (The real problem!)

I want to construct a bilayered system in which consist of two layers interacting with each other. In the following code, I proceed with only one layer because the other will be analogous with its respectives lists.

So, we are going to check in both layers separated when the points are near by a distance of 1.

m = 1.;
n = m + 1;
t=1;

(*Basis vectors*)
a1 = {Sqrt[3]/2, -1/2}*Sqrt[3];
a2 = {Sqrt[3]/2, 1/2}*Sqrt[3];
(*Unit cell vectors*)
t1 = m*a1 + n*a2;
t2 = (n + m)*a1 - m*a2;
(*k vector (Eigenvalues will deppend upon it)*)
k = {kx, ky};

(*This function is responsible to distribute the hexagon points in x,y plane*)
pts[x_, y_] := 
  Block[{j, k}, 
   Flatten[Table[{{Sqrt[3] j, 1 k} Sqrt[3]}, {j, 0, x}, {k, 0, y}], 
    2]];

hexagon = {{0, 1/2} Sqrt[3], {Sqrt[3]/6, 1} Sqrt[
    3], {Sqrt[3]/2, 1} Sqrt[3], {(2 Sqrt[3])/3, 1/2} Sqrt[
    3], {Sqrt[3]/2, 0} Sqrt[3], {Sqrt[3]/6, 0} Sqrt[3]};

(*Angle of rotation between the layers*)
\[Theta] = ArcCos[1/2*(n^2 + 4 n*m + m^2)/(n^2 + n*m + m^2)];




rot = RotationTransform[N[\[Theta]], {0, 0}]; 

(*----Constructing the layers----*)
LAYER1 = TranslationTransform[# - hexagon[[4]]][hexagon] & /@ 
   pts[(t1 + t2)[[1]]/2, -(t1 + t2).RotationMatrix[-60 Degree][[2]]/2];

LAYER2 = rot /@ TranslationTransform[# - hexagon[[3]]][hexagon] & /@ 
   pts[(t1 + t2)[[1]]/2, -(t1 + t2).RotationMatrix[-60 Degree][[2]]/2];

(*--NOW I'LL PROCEED ONLY WITH THE ROTATED PART--*)

L2 = Flatten[LAYER2, 1];


(*---ROTATED---*)
ROT = DeleteDuplicates[Table[
    If[
     RegionMember[
      ConvexHullMesh[({t1, t2, t1 + t2, {0, 0}} + 
          0.1).RotationMatrix[-60 Degree]], L2[[i]]]
     
     , L2[[i]], Nothing]
    
    , {i, 1, Length[L2]}]];


(*--OUTSIDE THE UNIT CELL TWISTED--*)
OUTREGIONTWISTED = 
  DeleteDuplicates[
   Table[If[
     RegionMember[
      ConvexHullMesh[({{-2, 2}, t1 + {2, 2}, (t1 + t2) + {2, -2}, 
           t2 + {-2, -2}, {-2, 2}} + 0.1).RotationMatrix[-60 Degree]],
       L2[[u]]], L2[[u]], Nothing], {u, 1, Length[L2]}]];

(*--FUNCTION TO REMOVE THE DUPLICATED POINTS INSIDE OUTREGIONTWISTED--*)
SIZOR[a_List, b_List] := Module[{c, o, x}, c = Join[b, a];
  o = Ordering[c];
  x = 1 - 2 UnitStep[-1 - Length[b] + o];
  x = FoldList[Max[#, 0] + #2 &, x];
  x[[o]] = x;
  Pick[c, x, -1]]

(*Points to construct the periodic system*)
PERIODICTWISTED = SIZOR[OUTREGIONTWISTED, ROT];


(*TWISTED*)

TWISTEDMATRIX= Table[0, {i, 1, Length[ROT]}, {j, 1, Length[ROT]}];

(*-----------ELEMENTS RELATED TO THE UNITCELL AND OUTSIDE-------------*)

Table[
  If[

(*----CHECK THE FIRST DISTANCE BETWEEN THE POINTS INSIDE THE UNIT \
CELL AND OUTSITE WITH THAT WE GOT THE "j" POSITION IN MATRIX---------*)

   
   EuclideanDistance[PERIODICTWISTED[[h]], ROT[[j]]] == 1 && 

(*----------ANOTHER PERIODIC CONTITION TO OBTAIN THE "i" POSITION \
IN MATRIX---------*)
    
    PERIODICTWISTED[[h]] - (nn t1 + mm t2) == ROT[[i]],
   
(*--------THEN THE ELEMENT POSITION----------*)
   
   TWISTEDMATRIX[[i, j]] = t*Exp[I*k.(ROT[[j]] - PERIODICTWISTED[[h]])];
   ,
   Nothing]
  
  , {nn, -1, 1}, {mm, -1, 1}, {i, 1, Length[ROT]}, {j, 1, 
   Length[ROT]}, {h, 1, Length[PERIODICTWISTED]}];

(*------------ELEMENTS RELATED WITH INSIDE THE UNIT CELL------------*)

(*HERE I CAN USE WHAT kglr SUGGESTED*)
(*MAT1=t*(1-Unitize[DistanceMatrix[ROT]-1])*)
(*AND THEN DO MAT1+TWISTEDMATRIX*)

Flatten[Table[If[EuclideanDistance[ROT[[i]], ROT[[j]]] == 1,
    
    TWISTEDMATRIX[[i, j]] = t
    
    , Nothing], {i, 1, Length[ROT]}, {j, 1, Length[ROT]}], 1];

Then, the idea is to combine all the matrices doing ArrayFlatten[{{TWISTEDMATRIX,SOME_INTERACTION_MATRIX},{SOME_INTERACTION_MATRIX,NOTTWISTEDMATRIX}}], and solve the Eigenvalues for kx and ky, for example from 0 to 2Pi to make a surfaceplot (ListPlot3D) or the other option is doing as follows

MAT[kx_,ky_]=Table[TWISTEDMATRIX[[i,j]],{i,1,Length[TWISTEDMATRIX]},{j,1,Length[TWISTEDMATRIX]}];

ClearAll[\[CapitalGamma], K, M, b1, b2, T1, T2, B1, B2];

T2 = Flatten[ArrayReshape[t2, {1, 3}]];
T1 = Flatten[ArrayReshape[t1, {1, 3}]];

B2 = 2 Pi (Cross[{0, 0, 1}, T1]/T2.Cross[{0, 0, 1}, T1]);
B1 = 2 Pi (Cross[T2, {0, 0, 1}]/T1.Cross[T2, {0, 0, 1}]);

(*--Reciprocal Vectors--*)
b2 = {B2[[1]], B2[[2]]};
b1 = {B1[[1]], B1[[2]]};


\[CapitalGamma] = {0, 0};
M = b1/2;
K = (b1 - b2)/3;
\[CapitalGamma]x = \[CapitalGamma][[1]];
\[CapitalGamma]y = \[CapitalGamma][[2]];
Kx = K[[1]];
Ky = K[[2]];
Mx = M[[1]];
My = M[[2]];
N1 = 500;


datK\[CapitalGamma] = Module[{k, kx, ky},
   Transpose[Table[
     k = -nn*
       Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2]/N1;
      kx = k*(Kx - \[CapitalGamma]x)/
       Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2]; 
     ky = k*(Ky - \[CapitalGamma]y)/
       Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2]; \
({k + 0.92, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 
      1, N1}]]];

dat\[CapitalGamma]M = Module[{k, kx, ky}, Transpose[Table[
     k = nn*
       Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2]/N1;
      kx = \[CapitalGamma]x + 
       k*(Mx - \[CapitalGamma]x)/
        Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2]; 
     ky = \[CapitalGamma]y + 
       k*(My - \[CapitalGamma]y)/
        Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2]; \
({Sqrt[(\[CapitalGamma]x - Kx)^2 + (\[CapitalGamma]y - Ky)^2] + 
          k, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 
      1, N1}]]];

datMK = Module[{k, kx, ky},
   Transpose[Table[
     kx = Mx + k*(Kx - Mx)/Sqrt[(Kx - Mx)^2 + (Ky - My)^2]; 
     ky = My + k*(Ky - My)/Sqrt[(Kx - Mx)^2 + (Ky - My)^2]; 
     k = nn*Sqrt[(Kx - Mx)^2 + (Ky - My)^2]/
       N1; ({Sqrt[(\[CapitalGamma]x - Kx)^2 + (\[CapitalGamma]y - 
             Ky)^2] + 
          Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2] \
+ k, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 1, 
      N1}]]];

ListLinePlot[Union[datK\[CapitalGamma], dat\[CapitalGamma]M, datMK]]

Edit3: Implementing the Henrik suggestion

The following code is what I got above together with what Henrik have developed, for both cases, for example for m=5 I get imaginary eigenvalues, as can be seen in Listplot corresponding to blank spaces. First of all, our matrix is supposed to be Hermitian, so, only real eigenvalues are expected. My guess is that some elements in upper triangular part are numerically different from the botom triangular part of MAT[kx,ky].

m = 5.; (*for example*)
n = m + 1;
a1 = {Sqrt[3]/2, -1/2}*Sqrt[3];
a2 = {Sqrt[3]/2, 1/2}*Sqrt[3];
r = n*a1 + m*a2;
t1 = m*a1 + n*a2;
t2 = (n + m)*a1 - m*a2;


pts[x_, y_] := 
  Block[{j, k}, 
   Flatten[Table[{{Sqrt[3] j, 1 k} Sqrt[3]}, {j, 0, x}, {k, 0, y}], 
    2]];
\[Theta] = ArcCos[1/2*(n^2 + 4 n*m + m^2)/(n^2 + n*m + m^2)];
hexagon = {{0, 1/2} Sqrt[3], {Sqrt[3]/6, 1} Sqrt[
    3], {Sqrt[3]/2, 1} Sqrt[3], {(2 Sqrt[3])/3, 1/2} Sqrt[
    3], {Sqrt[3]/2, 0} Sqrt[3], {Sqrt[3]/6, 0} Sqrt[3]};

rot = RotationTransform[N[\[Theta]], {0, 0}]; 

LAYER1 = TranslationTransform[# - hexagon[[4]]][hexagon] & /@ 
   pts[(t1 + t2)[[1]]/2, -(t1 + t2).RotationMatrix[-60 Degree][[2]]/2];

LAYER2 = rot /@ TranslationTransform[# - hexagon[[3]]][hexagon] & /@ 
   pts[(t1 + t2)[[1]]/2, -(t1 + t2).RotationMatrix[-60 Degree][[2]]/2];



L2 = Flatten[LAYER2, 1];

ROT = DeleteDuplicates[Table[
    If[
     RegionMember[
      ConvexHullMesh[({t1, t2, t1 + t2, {0, 0}} + 
          0.1).RotationMatrix[-60 Degree]], L2[[i]]]
     
     , L2[[i]], Nothing]
    
    , {i, 1, Length[L2]}]];


OUTREGIONTWISTED = 
  DeleteDuplicates[
   Table[If[
     RegionMember[
      ConvexHullMesh[({{-2, 2}, t1 + {2, 2}, (t1 + t2) + {2, -2}, 
           t2 + {-2, -2}, {-2, 2}} + 0.1).RotationMatrix[-60 Degree]],
       L2[[u]]], L2[[u]], Nothing], {u, 1, Length[L2]}]];

SIZOR[a_List, b_List] := Module[{c, o, x}, c = Join[b, a];
  o = Ordering[c];
  x = 1 - 2 UnitStep[-1 - Length[b] + o];
  x = FoldList[Max[#, 0] + #2 &, x];
  x[[o]] = x;
  Pick[c, x, -1]]

PERIODICTWISTED = SIZOR[OUTREGIONTWISTED, ROT];



R = Developer`ToPackedArray[ROT];

P = Developer`ToPackedArray[PERIODICTWISTED];

\[Epsilon] = 0.00001;
A = ConstantArray[0., {Length[R], Length[R]}];

Do[If[Abs[EuclideanDistance[R[[i]], R[[j]]] - 1] <= \[Epsilon], 
   A[[i, j]] += t Exp[I k.(R[[j]] - R[[i]])];], {i, 1, Length[R]}, {j,
    1, Length[R]}];

Do[If[Abs[EuclideanDistance[P[[h]], R[[j]]] - 1] <= \[Epsilon], 
    Do[If[Norm[P[[h]] - (R[[i]] + (nn t1 + mm t2))] <= \[Epsilon], 
       A[[i, j]] += t Exp[-I k.(P[[h]] - R[[j]])];];, {nn, -1, 
      1}, {mm, -1, 1}, {i, 1, Length[R]}]];, {j, 1, Length[R]}, {h, 1,
    Length[P]}];

Or in alternative way

(*Make sure that `SparseArray` assemble is additive.*)
SetSystemOptions[
  "SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}];
(*Create NearestFunction for R.*)
Rnf = Nearest[R -> "Index"];
(*Find all pairs {i,j} satisfying \
Abs[EuclideanDistance[R[[i]],R[[j]]]-1]\[LessEqual]\[Epsilon].*)
{ilist1, jlist1} = 
  Transpose[
   Join @@ Map[Thread[{First[#], Rest[#]}] &, 
     Rnf[R, {\[Infinity], 1 + \[Epsilon]}]]];
(*Use the {i,j} pairs to assemble the matrix.*)
A1 = SparseArray[
   Transpose[{ilist1, jlist1}] -> 
    t Exp[I (R[[jlist1]] - R[[ilist1]]).k], {Length[R], Length[R]}];
stencil = Flatten[Table[(nn t1 + mm t2), {nn, -1, 1}, {mm, -1, 1}], 1];
(*Displacing the points in R by each vector in stencil.*)

T = Flatten[Outer[Plus, R, stencil, 1], 1];
(*Create NearestFunction for T.*)
Tnf = Nearest[T -> "Index"];
(*Find all pairs {j,h} satisfying \
Abs[EuclideanDistance[P[[h]],R[[j]]]-1]\[LessEqual]\[Epsilon].*)
{jcandidates, hcandidates} = 
  Transpose[
   Join @@ MapIndexed[Thread[{#1, First[#2]}] &, 
     Rnf[P, {\[Infinity], 1 + \[Epsilon]}]]];
(*For each pair {j,h} above,search for the corresponding i satisfying \
Norm[P[[h]]-(R[[i]]+(nn t1+mm t2))]\[LessEqual]\[Epsilon] by using \
the NearestFunction of T.*)

icandidates = 
  Quotient[Tnf[P[[hcandidates]], {\[Infinity], \[Epsilon]}] - 1, 
    Length[stencil]] + 1;
(*Create all {i,j,h} triplets.*)
{ilist2, jlist2, hlist2} = 
  Transpose[
   Join @@ MapThread[
     Thread[{#1, #2, #3}] &, {icandidates, jcandidates, hcandidates}]];
(*Use these triplets to assemble the matrix.*)
A2 = SparseArray[
   Transpose[{ilist2, jlist2}] -> 
    t Exp[-I (P[[hlist2]] - R[[jlist2]]).k], {Length[R], Length[R]}];

A = A1 + A2;

Now, we proceed to the calculation of the Eigenvalues as showed above

MAT[kx_, ky_] = Table[A[[i, j]], {i, 1, Length[A]}, {j, 1, Length[A]}];
ClearAll[\[CapitalGamma], K, M, b1, b2, T1, T2, B1, B2];

T2 = Flatten[ArrayReshape[t2, {1, 3}]];
T1 = Flatten[ArrayReshape[t1, {1, 3}]];

B2 = 2 Pi (Cross[{0, 0, 1}, T1]/T2.Cross[{0, 0, 1}, T1]);
B1 = 2 Pi (Cross[T2, {0, 0, 1}]/T1.Cross[T2, {0, 0, 1}]);

(*--Reciprocal Vectors--*)
b2 = {B2[[1]], B2[[2]]};
b1 = {B1[[1]], B1[[2]]};


\[CapitalGamma] = {0, 0};
M = b1/2;
K = (b1 - b2)/3;
\[CapitalGamma]x = \[CapitalGamma][[1]];
\[CapitalGamma]y = \[CapitalGamma][[2]];
Kx = K[[1]];
Ky = K[[2]];
Mx = M[[1]];
My = M[[2]];
N1 = 500;


datK\[CapitalGamma] = 
  Module[{k, kx, ky}, 
   Transpose[
    ParallelTable[
     k = -nn*Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \
\[CapitalGamma]y)^2]/N1;
     kx = 
      k*(Kx - \[CapitalGamma]x)/
        Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2];
     ky = 
      k*(Ky - \[CapitalGamma]y)/
        Sqrt[(Kx - \[CapitalGamma]x)^2 + (Ky - \[CapitalGamma]y)^2]; \
({k + 0.92, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 1, N1}]]];

dat\[CapitalGamma]M = 
  Module[{k, kx, ky}, 
   Transpose[
    ParallelTable[
     k = nn*Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \
\[CapitalGamma]y)^2]/N1;
     kx = \[CapitalGamma]x + 
       k*(Mx - \[CapitalGamma]x)/
         Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2];
     ky = \[CapitalGamma]y + 
       k*(My - \[CapitalGamma]y)/
         Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2]; \
({Sqrt[(\[CapitalGamma]x - Kx)^2 + (\[CapitalGamma]y - Ky)^2] + 
          k, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 1, N1}]]];

datMK = Module[{k, kx, ky}, 
   Transpose[
    ParallelTable[
     kx = Mx + k*(Kx - Mx)/Sqrt[(Kx - Mx)^2 + (Ky - My)^2];
     ky = My + k*(Ky - My)/Sqrt[(Kx - Mx)^2 + (Ky - My)^2];
     k = nn*
       Sqrt[(Kx - Mx)^2 + (Ky - My)^2]/
        N1; ({Sqrt[(\[CapitalGamma]x - Kx)^2 + (\[CapitalGamma]y - 
               Ky)^2] + 
          
          Sqrt[(Mx - \[CapitalGamma]x)^2 + (My - \[CapitalGamma]y)^2] \
+ k, #} &) /@ Sort[Eigenvalues[MAT[kx, ky]]], {nn, 1, N1}]]];

ListLinePlot[Union[datK\[CapitalGamma], dat\[CapitalGamma]M, datMK]]

Answer 1

Allow me to first post a refactored version of OP's code that is more expressive towards the intentions. Note that I replaced Table by Do where its return value was not used. Moreover I introduced shorter and more eye-friendly notation for the variables. It is also not a good idea to compared floating point numbers such as distances with the == operator. Instead, I introduced a tolerance ϵ and replaced equality checks by inequalities involving this tolerance. Finally, I switched the order of the two main loops because the (originally) second loop had easier logic.

m = 1.;
n = m + 1;
a1 = {Sqrt[3]/2, -1/2}*Sqrt[3];
a2 = {Sqrt[3]/2, 1/2}*Sqrt[3];
t1 = m*a1 + n*a2;
t2 = (n + m)*a1 - m*a2;
k = {kx, ky};

t = 2. Pi;
kx = 2;
ky = 3;

R = Developer`ToPackedArray[{{1.85714, 0.742307}, {2.78571, 1.11346}, {1.07143, 1.3609}, {1.21429, 2.35064}, {2.14286, 2.72179}, {2.92857, 2.1032}, {1.5, 4.33013}, {2.28571, 3.71154}, {4.64286, 1.85577}, {3.85714, 2.47436}, {4., 3.4641}, {4.92857, 3.83526}, {3.21429, 4.08269}, {5.07143, 4.825}}];

P = Developer`ToPackedArray[{{-0.928571, -0.371154}, {0., 0.}, {0.785714, -0.61859}, {-1.71429, 0.247436}, {-1.57143, 1.23718}, {-0.642857, 1.60833}, {0.142857, 0.989743}, {-0.5, 2.59808}, {1.71429, -0.247436}, {3.57143, 0.494872}, {3.42857, -0.494872}, {2.5, -0.866025}, {0.428571, 2.96923}, {0.571429, 3.95897}, {-0.214286, 4.57756}, {0.857143, 5.93846}, {1.64286, 5.31987}, {4.5, 0.866025}, {5.57143, 2.22692}, {6.35714, 1.60833}, {6.21429, 0.61859}, {5.28571, 0.247436}, {5.71429, 3.21667}, {3.35714, 5.07243}, {4.28571, 5.44359}, {2.57143, 5.69102}, {2.71429, 6.68077}, {4.42857, 6.43333}, {6.64286, 3.58782}, {6.78571, 4.57756}, {7.71429, 4.94872}, {6., 5.19615}, {6.14286, 6.1859}, {7.85714, 5.93846}}];

ϵ = 0.00001;
A = ConstantArray[0., {Length[R], Length[R]}];

Do[
  If[
   Abs[EuclideanDistance[R[[i]], R[[j]]] - 1] <= ϵ
   ,
   A[[i, j]] += t Exp[I k.(R[[j]] - R[[i]])];
   ],
  {i, 1, Length[R]}, {j, 1, Length[R]}];

Do[
  If[
    Abs[EuclideanDistance[P[[h]], R[[j]]] - 1] <= ϵ
    ,
    Do[
     If[
       Norm[P[[h]] - (R[[i]] + (nn t1 + mm t2))] <= ϵ
       ,
       A[[i, j]] += t Exp[-I k.(P[[h]] - R[[j]])];
       ];
     , {nn, -1, 1}, {mm, -1, 1}, {i, 1, Length[R]}]
    ];
  , {j, 1, Length[R]}, {h, 1, Length[P]}];

This is my reimplementation of the first loop: Instead of building a Graph first (what I suggested in the comments), I simply generate a NearestFunction for the set R and use it to generate the sparsitit pattern of the matrix to assemble (i.e., the list ilist1, jlist1). The nonzero values of that matrix can then be computed easily. Finally, I use SparseArray to assemble the matrix. Most entries in the matrix are equal to 0., so a SparseArray is also a good way to safe memory.

(*Make sure that `SparseArray` assemble is additive.*)
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}];
(*Create NearestFunction for R.*)
Rnf = Nearest[R -> "Index"];
(*Find all pairs {i,j} satisfying Abs[EuclideanDistance[R[[i]],R[[j]]]-1]≤ϵ.*)
{ilist1, jlist1} = Transpose[Join @@ Map[Thread[{First[#], Rest[#]}] &, Rnf[R, {∞, 1 + ϵ}]]];
(*Use the {i,j} pairs to assemble the matrix.*)
A1 = SparseArray[
    Transpose[{ilist1, jlist1}] -> t Exp[I (R[[jlist1]] -R[[ilist1]]).k], 
    {Length[R], Length[R]}
    ];

The second loop has a more complicated structure because two checks are involved. But the idea is basically the same. I filter first with repect to one check to obtain all potential pairs {i,h} and then see which j I find by a second Nearest query.

stencil = Flatten[Table[(nn t1 + mm t2), {nn, -1, 1}, {mm, -1, 1}], 1];
(*Displacing the points in R by each vector in stencil.*)

T = Flatten[Outer[Plus, R, stencil, 1], 1];
(*Create NearestFunction for T.*)
Tnf = Nearest[T -> "Index"];
(*Find all pairs {j,h} satisfying Abs[EuclideanDistance[P[[h]],R[[j]]]-1]≤ϵ.*)
{jcandidates, hcandidates} =  Transpose[Join @@ MapIndexed[Thread[{#1, First[#2]}] &, Rnf[P, {∞, 1 + ϵ}]]];
(*For each pair {j,h} above, search for the corresponding i satisfying Norm[P[[h]]-(R[[i]]+(nn t1+mm t2))]≤ϵ by using the NearestFunction of T.*)

icandidates = Quotient[Tnf[P[[hcandidates]], {∞, ϵ}] - 1, Length[stencil]] + 1;
(*Create all {i,j,h} triplets.*)
{ilist2, jlist2, hlist2} = Transpose[Join @@ MapThread[ Thread[{#1, #2, #3}] &, {icandidates, jcandidates, hcandidates}]];
(*Use these triplets to assemble the matrix.*)
A2 = SparseArray[
    Transpose[{ilist2, jlist2}] -> t Exp[-I (P[[hlist2]] - R[[jlist2]]).k],
    {Length[R], Length[R]}
    ];

The final result is the sum of A1 and A2:

Max[Abs[SparseArray[A1 + A2 - A]]]

For large inputs on can probably safe some time by assembling the some in one go as follows:

B = SparseArray[
   Rule[
    Transpose[{Join[ilist1, ilist2], Join[jlist1, jlist2]}],
    Join[t Exp[I (R[[jlist1]] - R[[ilist1]]).k], 
     t Exp[-I (P[[hlist2]] - R[[jlist2]]).k]]
    ],
   {Length[R], Length[R]}
   ];

Max[Abs[SparseArray[B - A]]]

Answer 2

Here is a MMA like implementation using "Outer"

n = 5;
TRIANGLE = Triangle[{{0, 0}, {1, 1}, {2, 0}}];
pts = RandomPoint[TRIANGLE, n];
dist = EuclideanDistance[pts[[1]], pts[[2]]];
Outer[If[EuclideanDistance[#1, #2] == dist, r, 0] &, pts, pts, 1]

(* {{0, r, 0, 0, 0}, {r, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 
  0}, {0, 0, 0, 0, 0}} *)

However, because the probability of getting 2 points exactly a given distance apart is zero, you will always only get 2 "r" in your matrix