I didn’t try running the model, but some general ideas:
This looks very suspicious - ignoring the sum-to zero constraint on mu it should be completely equivalent to something like
target += normal_lpdf(theta | X*gamma, sqrt(2));
i.e. mu probably can be integrated out and having it in the model could lead to problematic posterior geometry (I would expect to see strong correlations between corresponding elements of mu and theta). Couldn’t you transform the sum-to-zero constraint on mu into a constraint on theta and don’t treat mu as a parameter at all?
This form of sum-to-zero is potentially problematic. There has been quite an extensive discussion on this - probably more than you need to know: Test: Soft vs Hard sum-to-zero constrain + choosing the right prior for soft constrain
TLDR: try using sum(x) ~ normal(0, <small_number> * <num_lements>) or the code in this answer in the aformentioned thread + the one below (which is in my experience slightly superior). For more details read the whole topic.
It appears that the problem is that beta (as a parameter to partial_sum_lpmf) is defined as real, but you are indexing it as if it was array/vector. To define it as array, you would have real[] beta in the function signature.
Best of luck with your model!