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

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 = DeveloperToPackedArray[ROT];

P = DeveloperToPackedArray[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[
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[
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[
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]]


• Look up SparseArray Commented Jan 25, 2022 at 21:34
• I'have already tried using SparseArray, but the gain is not that big. I think the "problem" is when its reading the list, and checking every element to follow the condition. Commented Jan 25, 2022 at 21:48
• try MAT1 =r (1 - Unitize[DistanceMatrix[pts] - EuclideanDistance @@ pts[[;; 2]]])?
– kglr
Commented Jan 25, 2022 at 22:09
• yes, replace EuclideanDistance @@ pts[[;; 2]]] with 0.5
– kglr
Commented Jan 25, 2022 at 22:22
• Please, correct the posted code. There are undefined symbols and the result matrix does not make any sense. Commented Jan 27, 2022 at 21:04

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 = DeveloperToPackedArray[{{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 = DeveloperToPackedArray[{{-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]]]

• First of all, thank you for your contribution. I have edited the post and explained I little better what I'm calculating. The problem is for m>5 we got imaginary values to datMK. Commented Jan 29, 2022 at 20:41
• In summary, you suggestion works pretty well and is faster than two loops to constructing the matrices. So, for the proposite of this Question it is well done. But I got other problem with eigenvalues in which maybe I should open another question... Commented Jan 29, 2022 at 20:48
• You're welcome. And indeed, a new thread would be better. Commented Jan 29, 2022 at 20:51
• Should I call Needs["Developer"] to use this recourse? Commented Jan 29, 2022 at 23:55
• Not necessary. The "Developer" context should be already present right from the start of the kernel. Commented Jan 30, 2022 at 0:55

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

• I made an edit to include a more realistic system. Commented Jan 27, 2022 at 1:11