Question❓
Please bear with me I am not an experienced programmer and I might be doing something silly.

I am trying to integrate
p = 200.0
function psi(du,u,p,t)
du[1] = -im * p * u[1]
du[2] = -im * p * u[2]
end

function g(u,resid)
resid[1] = abs(sol.u[1]) - 1
resid[2] = abs(sol.u[2]) - 1
end

cb = ManifoldProjection(g)

sol = solve(prob,callback = cb)
and I get

ArgumentError: Cannot create a dual over scalar type ComplexF64. If the type behaves as a scalar, define ForwardDiff.can_dual(::Type{ComplexF64}) = true.

Stacktrace:
[1] throw_cannot_dual(V::Type)
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/dual.jl:41
[2] ForwardDiff.Dual{ForwardDiff.Tag{DiffEqCallbacks.NonAutonomousFunction{typeof(g), false}, ComplexF64}, ComplexF64, 2}(value::ComplexF64, partials::ForwardDiff.Partials{2, ComplexF64})

After many web-searches it seems to me that this is an old problem since 2018. My question is motivated by the stochastic Schroedinger equation where one needs an efficient integrator capable of preserving the norm of a complex vector for ensemble simulations.
Many thanks!