How to solve this trigonometric complex ODE system?

It is convenient and also provides additional insight to rescale t by s, for which the ODE system becomes,

ss = NDSolveValue[{(p''[t] - 1/2 Sin[2 p[t]] (q'[t])^2) == c Sin[p[t]] q'[t] - 
    a s Sin[p[t]]/I, 
    (Sin[p[t]]^2 q''[t] + Sin[2 p[t]] p'[t] q'[t]) == -c Sin[p[t]] p'[t], 
    p[0] == 2 π/3, p[max/s] == 2 π/3, q[0] == π/4, q[max/s] == π/4}, {p, q}, {t, 0, max/s}]
    // Flatten;

Thus, s now can be absorbed into a in the first ODE, if desired, and does not appear at all in the second ODE. Note, though, that the range of integration now is {0, max/s}, so s is not completely eliminated. (b has been dropped here, because it is set to zero in the question. However, it can be reintroduced easily, if desired.) It also is important to observe that the coefficient of q''[t] is proportional to Sin[p[t]]^2. Thus, the ODE system becomes singular whenever p becomes equal to π anywhere in the domain of integration.

User64494 used Maple to obtain a solution for

max = 1/2; s = 1/100; a = 1; b = 0; c = 1/2;

It can be obtained with Mathematica by

ss = NDSolveValue[{(p''[t] - 1/2 Sin[2 p[t]] (q'[t])^2) == c Sin[p[t]] q'[t] - 
    a s Sin[p[t]]/I, 
    (Sin[p[t]]^2 q''[t] + Sin[2 p[t]] p'[t] q'[t]) == -c Sin[p[t]] p'[t], 
    p[0] == 2 π/3, p[max/s] == 2 π/3, q[0] == π/4, q[max/s] == π/4}, {p, q}, 
    {t, 0, max/s}, Method -> {"Shooting", "ImplicitSolver" -> {"Newton", 
    "StepControl" -> "LineSearch"}, "StartingInitialConditions" -> 
    p'[0] == 1/100 - 20/100 I, q'[0] == 1/100}}] // Flatten;
Plot[Evaluate@ReIm@Through[ss[t]], {t, 0, max/s}]

Here, the real and imaginary parts of p are blue and orange, and of q are green and red. (See the answer to 124457 for a discussion of "ImplicitSolver" -> {"Newton", "StepControl" -> "LineSearch"}, without which NDSolve crashes.) With modest effort max can be increased to 55/100 by means of "StartingInitialConditions" -> {p'[0] == 4/100 - 29/100 I, q'[0] == 1/100 - 8/100 I}.

Visibly, the ODE system almost is singular, and I have been unable to obtain a solution for max == 56/100. There may, of course, be solutions for stable islands of larger max, but I have not tried to find them. In that regard, it is easy to obtain solutions for max == 40/100 and max == 42/100 but not for max == 41/100, again because the ODE system is singular there.

An alternative approach to finding a solution for the specific parameters in the question is to begin with small a and gradually increase it. Unfortunately, I was unable to increase a beyond 4 10^-2 (and max == 10), for which with "StartingInitialConditions" -> {p'[0] == 0.00026069739639169307 - 0.1478258792803156 I, q'[0] == -0.1646363401819224 + 0.01211739772818876 I}, the solution is

Larger values of a result in

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

but increasing WorkingPrecision does not help and is painfully slow besides. In closing, I have been able to find solutions for many parameters but not for those in the question.

Again, let me give a solution based on finite difference method (FDM). I'll use pdetoae for the generation of difference equation.

Probably you've already noticed that, finding proper initial guess for this problem is no longer as easy as in e.g. your previous problem. We need to guess it in a cleverer way. Inspired by the observations of user64494 and bbgodfrey i.e. the system can be solved more easily when max isn't that large, I found that using solutions at smaller max as initial guess seems to work wonders.

Clear@max
{s = 1/100, c = 1/2, points = 900, difforder = 4};
initialguess = ConstantArray[1/2, points];
{a, b, {{p0, pmax}, {q0, qmax}}} = {1, 0, {{2 Pi/3, 2 Pi/3}, {Pi/4, Pi/4}}};
{eq, bc} = {{I s (p''[t] - 1/2 Sin[2 p[t]] (q'[t])^2) == 
     I c Sin[p[t]] q'[t] - a Sin[p[t]] + b Cos[p[t]] Cos[q[t]], 
    I s (Sin[p[t]]^2 q''[t] + Sin[2 p[t]] p'[t] q'[t]) == -I c Sin[p[t]] p'[t] - 
      b Sin[p[t]] Sin[q[t]]}, {p[0] == p0, p[max] == pmax, q[0] == q0, q[max] == qmax}};

grid = Array[# &, points, {0, max}];
ptoafunc = pdetoae[{p, q}[t], grid, difforder];
delete = #[[2 ;; -2]] &;
ae = delete /@ ptoafunc[eq];
solandmaxlst = NestList[Function[guessandmax, Block[{guess, max, guessfunc, sollst},
      {guess, max} = guessandmax;
      guessfunc = ListInterpolation[#, {0, max}] & /@ guess;
      sollst = 
       ArrayReshape[
        FindRoot[{ae, bc}, 
          Flatten[#, 1] &@
           Table[{{p, q}[[var]][index], guessfunc[[var]][index]}, {var, 2}, {index, 
             grid}], MaxIterations -> 1000][[All, -1]], {2, points}];
      {sollst, max + 1}]], {{initialguess, initialguess}, 1}, 10]; // AbsoluteTiming
(* {21.842398, Null} *)

funclst = Table[
   ListInterpolation[#, {0, lst[[2]] - 1}] & /@ lst[[1]], {lst, solandmaxlst // Rest}];

We can observe how the solution evolutes by:

plotlst = Table[
  Plot[{Re@#[t], Im@#[t]}, {t, 0, #["Domain"][[1, -1]]}, PlotRange -> All] & /@ pq, {pq, 
   funclst}];
ListAnimate@plotlst

The final solution is

plotlst[[-1]] // GraphicsColumn

Further check by adjusting points etc. shows the solution seems to be reliable.

Finally the following is the solution for {a, b, {{p0, pmax}, {q0, qmax}}} = {1/2, Sqrt[3]/2, {{3 Pi/4, 3 Pi/4}, {Pi/4, 3 Pi/4}}}: