Solution: The Generic Magic Square

The code...

declare
fun {MagicSquare N}
   NN=N*N
   L1N={List.number 1 N 1} % [1..N]
in
   proc {$ Square}
      fun {Field I J}
         Square.((I-1)*N + J)
      end
      proc {Assert F}
         {FD.sum {Map L1N F} '=:' Sum}
      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}
      %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 trees for the first three solutions at N=5 with split (left) and naive (right) distribution strategy.

The first three magic squares yopu get with domain splitting are

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

well, they seem to come in packs somehow, don't they?