About building a list until it meets conditions

step(Ls) --> [add5(Ls, L)], step(L).

This doesn't do what you want. It describes a list element of the form add5(Ls, L). Presumably Ls is bound to some value when you get here, but L is not bound. L would become bound if Ls were a non-empty list of the correct form, and you executed the goal add5(Ls, L). But you are not executing this goal. You are storing a term in a list. And then, with L completely unbound, some part of the code that expects it to be bound to a list will throw this error. Presumably that sort/2 call is inside all_different/1.

Edit: There are some surprisingly complex or inefficient solutions posted here. I think both DCGs and CLP are overkill here. So here's a relatively simple and fast one. For enforcing the correct 2/10/14 order this uses a state argument to keep track of which ones we have seen in the correct order:

puzzle(Solution) :-
    run([0], seen_nothing, ReverseSolution),
    reverse(ReverseSolution, Solution).
    
run(FinalList, seen_14, FinalList).
run([Head | Tail], State, Solution) :-
    dif(State, seen_14),
    step(Head, Next),
    \+ member(Next, Tail),
    state_next(State, Next, NewState),
    run([Next, Head | Tail], NewState, Solution).
    
step(Number, Next) :-
    (   Next is Number + 5
    ;   Next is Number + 7
    ;   nth_integer_root_and_remainder(2, Number, Next, 0) ),
    Next =< 60,
    dif(Next, Number).  % not strictly necessary, added by request

    
state_next(State, Next, NewState) :-
    (   State = seen_nothing,
        Next = 2
    ->  NewState = seen_2
    ;   State = seen_2,
        Next = 10
    ->  NewState = seen_10
    ;   State = seen_10,
        Next = 14
    ->  NewState = seen_14
    ;   NewState = State ).

Timing on SWI-Prolog:

?- time(puzzle(Solution)), writeln(Solution).
% 13,660,415 inferences, 0.628 CPU in 0.629 seconds (100% CPU, 21735435 Lips)
[0,5,12,17,22,29,36,6,11,16,4,2,9,3,10,15,20,25,30,35,42,49,7,14]
Solution = [0, 5, 12, 17, 22, 29, 36, 6, 11|...] .

The repeated member calls to ensure no duplicates make up the bulk of the execution time. Using a "visited" table (not shown) takes this down to about 0.25 seconds.

Edit: Pared down a bit further and made 100x faster:

prev_next(X, Y) :-
    between(0, 60, X),
    (   Y is X + 5
    ;   Y is X + 7
    ;   X > 0,
        nth_integer_root_and_remainder(2, X, Y, 0) ),
    Y =< 60.

moves(Xs) :-
    moves([0], ReversedMoves),
    reverse(ReversedMoves, Xs).
    
moves([14 | Moves], [14 | Moves]) :-
    member(10, Moves).
moves([Prev | Moves], FinalMoves) :-
    Prev \= 14,
    prev_next(Prev, Next),
    (   Next = 10
    ->  member(2, Moves)
    ;   true ),
    \+ member(Next, Moves),
    moves([Next, Prev | Moves], FinalMoves).

?- time(moves(Solution)), writeln(Solution).
% 53,207 inferences, 0.006 CPU in 0.006 seconds (100% CPU, 8260575 Lips)
[0,5,12,17,22,29,36,6,11,16,4,2,9,3,10,15,20,25,30,35,42,49,7,14]
Solution = [0, 5, 12, 17, 22, 29, 36, 6, 11|...] .

The table of moves can be precomputed (enumerate all solutions of prev_next/2, assert them in a dynamic predicate, and call that) to gain another millisecond or two. Using a CLP(FD) instead of "direct" arithmetic makes this considerably slower on SWI-Prolog. In particular, Y in 0..60, X #= Y * Y instead of the nth_integer_root_and_remainder/4 goal takes this up to about 0.027 seconds.

Given that the question seems to have shifted from using DCGs to solving the puzzle, I thought I might post a more efficient approach. I am using clp(fd) on SICStus, but I included a modified version that should work with clpz on Scryer (replacing table/2 with my_simple_table/2).

:- use_module(library(clpfd)).
:- use_module(library(lists)).

move(X,Y):-
    (
      X+5#=Y
    ;
      X+7#=Y
    ;
      X#=Y*Y
    ).

move_table(Table):-
    findall([X,Y],(
            X in 0..60,
            Y in 0..60,
            move(X,Y),
            labeling([], [X,Y])
         ),Table).
      

% Naive version
%%post_move(X,Y):- move(X,Y).
%%
% SICSTUS clp(fd)
%%post_move(X,Y):-
%%  move_table(Table),
%%  table([[X,Y]],Table).
%%
% clpz is mising table/2
post_move(X,Y):-
    move_table(Table),
    my_simple_table([[X,Y]],Table).

my_simple_table([[X,Y]],Table):-
      transpose(Table, [ListX,ListY]),
      element(N, ListX, X),
      element(N, ListY, Y).


post_moves([_]):-!.
post_moves([X,Y|Xs]):-
    post_move(X,Y),
    post_moves([Y|Xs]).

state(N,Xs):-
    length(Xs,N),
    domain(Xs, 0, 60),
    all_different(Xs),
    post_moves(Xs),
    % ordering: 0 is first, 2 comes before 10, and 14 is last.
    Xs=[0|_],
    element(I2, Xs, 2),
    element(I10, Xs, 10),
    I2#<I10,
    last(Xs, 14).

try_solve(N,Xs):-
    state(N, Xs),
    labeling([ffc], Xs).
try_solve(N,Xs):-
    N1 is N+1,
    try_solve(N1,Xs).


solve(Xs):-
    try_solve(1,Xs).

Two notes of interest:

• It is much more efficient to create a table of the possible moves and use the table/2 constraint rather than posting a disjunction of constraints. Note that we are recreating the table every time we post it, but we might as well create it once and pass it along.
• This is using the element/3 constraint to find and constraint the position of the numbers of interest (in this case just 2 and 10, because we can fix 14 to be last). Again, this is more efficient than checking the order as filtering after solving the constraint problem.

Edit:

Here is an updated version to conform to the bounty constraints (predicate names, -hopefully- SWI-compatible, create the table only once):

:- use_module(library(clpfd)).
:- use_module(library(lists)).

generate_move_table(Table):-
    X in 0..60,
    Y in 0..60,
    (    X+5#=Y 
    #\/  X+7#=Y 
    #\/  X#=Y*Y 
    ),
    findall([X,Y],labeling([], [X,Y]),Table).
      
%post_move(X,Y,Table):- table([[X,Y]],Table). %SICStus
post_move(X,Y,Table):- tuples_in([[X,Y]],Table). %swi-prolog
%post_move(X,Y,Table):- my_simple_table([[X,Y]],Table). %scryer

my_simple_table([[X,Y]],Table):- % Only used as a fall back for Scryer prolog
    transpose(Table, [ListX,ListY]),
    element(N, ListX, X),
    element(N, ListY, Y).

post_moves([_],_):-!.
post_moves([X,Y|Xs],Table):-
    post_move(X,Y,Table),
    post_moves([Y|Xs],Table).

puzzle_(Xs):-
    generate_move_table(Table),
    
    N in 4..61, 
    indomain(N),
    length(Xs,N),
    
    %domain(Xs, 0, 60), %SICStus
    Xs ins 0..60, %swi-prolog, scryer
    
    all_different(Xs),
    post_moves(Xs,Table),
    
    % ordering: 0 is first, 2 comes before 10, 14 is last.
    Xs=[0|_],
    element(I2, Xs, 2),
    element(I10, Xs, 10),
    I2#<I10,
    last(Xs, 14).

label_puzzle(Xs):-
    labeling([ffc], Xs).

solve(Xs):-
    puzzle_(Xs),
    label_puzzle(Xs).

I do not have SWI-prolog installed so I can't test the efficiency requirement (or that it actually runs at all) but on my machine and with SICStus, the new version of the solve/1 predicate takes 16 to 31 ms, while the puzzle/1 predicate in Isabelle's answer (https://stackoverflow.com/a/65513470/12100620) takes 78 to 94 ms.

As for elegance, I guess this is in the eye of the beholder. I like this formulation, it is relatively clear and is showcasing some very versatile constraints (element/3, table/2, all_different/1), but one drawback of it is that in the problem description the size of the sequence (and hence the number of FD variables) is not fixed, so we need to generate all sizes until one matches. Interestingly, it appears that all the solutions have the very same length and that the first solution of puzzle_/1 produces a list of the right length.

An alternative that uses dcg only to build the list. The 2,10,14 constraint is checked after building the list, so this is not optimal.

num(X) :- between(0, 60, X).

isqrt(X, Y) :- nth_integer_root_and_remainder(2, X, Y, 0). %SWI-Prolog

% list that ends with an element.
list([0], 0) --> [0].
list(YX, X) --> list(YL, Y), [X], { append(YL, [X], YX), num(X), \+member(X, YL),
                                    (isqrt(Y, X); plus(Y, 5, X); plus(Y, 7, X)) }.
soln(X) :-
    list(X, _, _, _),
    nth0(I2, X, 2), nth0(I10, X, 10), nth0(I14, X, 14),
    I2 < I10, I10 < I14.

?- time(soln(X)).
% 539,187,719 inferences, 53.346 CPU in 53.565 seconds (100% CPU, 10107452 Lips)
X = [0, 5, 12, 17, 22, 29, 36, 6, 11, 16, 4, 2, 9, 3, 10, 15, 20, 25, 30, 35, 42, 49, 7, 14] 

I tried a little magic set. The predicate path/2 does search a path without giving us a path. We can therefore use commutativity of +5 and +7, doing less search:

step1(X, Y) :- N is (60-X)//5, between(0, N, K), H is X+K*5,
         M is (60-H)//7, between(0, M, J), Y is H+J*7.
step2(X, Y) :- nth_integer_root_and_remainder(2, X, Y, 0).

:- table path/2.
path(X, Y) :- step1(X, H), (Y = H; step2(H, J), path(J, Y)).

We then use path/2 as a magic set for path/4:

step(X, Y) :- Y is X+5, Y =< 60.
step(X, Y) :- Y is X+7, Y =< 60.
step(X, Y) :- nth_integer_root_and_remainder(2, X, Y, 0).

/* without magic set */
path0(X, L, X, L).
path0(X, L, Y, R) :- step(X, H), \+ member(H, L), 
   path0(H, [H|L], Y, R).

/* with magic set */
path(X, L, X, L).
path(X, L, Y, R) :- step(X, H), \+ member(H, L), 
   path(H, Y), path(H, [H|L], Y, R).

Here is a time comparison:

SWI-Prolog (threaded, 64 bits, version 8.3.16)

/* without magic set */
?- time((path0(0, [0], 2, H), path0(2, H, 10, J), path0(10, J, 14, L))), 
   reverse(L, R), write(R), nl.
% 13,068,776 inferences, 0.832 CPU in 0.839 seconds (99% CPU, 15715087 Lips)
[0,5,12,17,22,29,36,6,11,16,4,2,9,3,10,15,20,25,30,35,42,49,7,14]

/* with magic set */
?- abolish_all_tables.
true.

?- time((path(0, [0], 2, H), path(2, H, 10, J), path(10, J, 14, L))), 
   reverse(L, R), write(R), nl.
% 2,368,325 inferences, 0.150 CPU in 0.152 seconds (99% CPU, 15747365 Lips)
[0,5,12,17,22,29,36,6,11,16,4,2,9,3,10,15,20,25,30,35,42,49,7,14]

Noice!

I managed to solve it without DCG's, it takes about 50 minutes on my machine to solve for the length N=24. I suspect this is because the order check is done for every list from scratch.

:- use_module(library(lists)).
:- use_module(library(clpz)).
:- use_module(library(dcgs)).
:- use_module(library(time)).

%% integer sqrt
isqrt(X, Y) :- Y #>= 0, X #= Y*Y.

before(X, Y, Z, L) :-
        %% L has a suffix [X|T], and T has a suffix of [Y|_].
        append(_, [X|T], L),
        append(_, [Y|TT], T),
        append(_, [Z|_], TT).

order(L) :- before(2,10,14, L).

game([X],X).
game([H|T], H) :- ((X #= H+5); (X #= H+7); (isqrt(H, X))), X #\= H, H #=< 60, X #=< 60,  game(T, X). % H -> X.

searchN(N, L) :- length(L, N), order(L), game(L, 0).