Solution: The Prime Magic Square

The code...

declare
fun {MagicSquare N}
   NN=N*N
   L1N={List.number 1 N 1} % [1..N]
in
   proc {$ Square}
      fun {Field A B} % modified to align bounds
         I=(A mod N) + N* (A div N)
         J=(B mod N) + N *(B div N)
      in
         Square.((I-1)*N + J)
      end
      proc {Assert F}
         {FD.sum {Map L1N F} '=:' Sum}
      end
      proc {AssertPrime F}
         F :: [ 1 2 3 5 7 11 13 17 19 23] % Primes between 1 and 25
      end
      
      Sum={FD.decl}
   in
      {FD.tuple square NN 1#NN Square}
      {FD.distinct Square}
      %Diagonals
    
      {Assert fun {$ I} {Field I I} end}

      {Assert fun {$ I} {Field I N+1-I} end}
      %Columns
      {For 1 N 1
       proc {$ I} {Assert fun {$ J} {Field I J} end} end}
      %Rows
      {For 1 N 1
       proc {$ J} {Assert fun {$ I} {Field I J} end} end}
      %Assert first Row to consist of primes
      {For 1 N 1 proc {$ I} {AssertPrime {Field 1 I}} end}
      %Assert other 'T' constituents to be primes
      {ForAll [[2 1] [2 3] [2 5] [3 3]]
       proc {$ I}
          [A B]=I
       in
          {AssertPrime {Field A B}}
       end
      }
      
      %Eliminate Symmetries
      {Field 1 1} <: {Field N N}
      {Field N 1} <: {Field 1 N}
      {Field 1 1} <: {Field N 1}
      %Redundant: sum of all fields = (number rows) * Sum
      NN*(NN+1) div 2 =: N*Sum

      {FD.distribute split Square}
   end
end

{ExploreOne {MagicSquare 5}}

...and the search tree

This first solution is

square( 3  7 13 19 23
       11 22 17 10  5
       12 21  2 14 16
       15  6 25 18  1
       24  9  8  4 20)