Hello,
As far as I know in MoveIt! there are no functions to get the acceleration of the end-effector. However, since the API is quite extensive, it reasonable to expect that I am wrong! If you - or anyone - find such function, I would be glad to know it :)
That said, to solve your problem you might consider using Orocos KDL. Here there is a list of available solvers that can compute forward/inverse kinematics at position and velocity level. Regarding acceleration, the solver is not available yet (the class is abstract). However that could be a starting point for a web search, and maybe someone actually implemented the solver somewhere.
As a 'quick and dirty' workaround, you can consider using this solver, which allows you to get the derivative of the jacobian. Note that in the documentation this is discouraged, as the solver is meant for internal use. In addition, note that using the formula a=J * qdd + Jd * qd is rather inefficient compared to a recursive implementation (e.g., the one used in the Newton-Euler algorithm for inverse dynamics).
Hope it helps!
EDIT
The recursive implementation I am referring to is the one that allows to recursively compute the acceleration of a kinematic chain's end-effector knowing joint positions, velocities and accelerations. It is usually exploited when dealing with Newton-Euler equations of dynamics. See, eg, this paper. The equations of insterest are (13) and (14), where V is a spatial velocity vector (ie, it contains both the linear and angular velocities). If you go through the document you will have all the theoretical information to implement the algorithm!
Regarding the statement 'formula a=Jadd+Jdqd is rather inefficient': KDL implementation of J and Jdot solvers is actually using a recursive formulation similar to the one I mentioned. If you get the acceleration using a=J * qdd+Jd * qd (I had mis-typed it before!) you will actually run two times the recursive equations to get the Jacobian and its derivative. Then, you will multiply the joint acceleration vector by J and the joint velocity vector by Jdot. In total you will perform 2x6xdof multiplications (dof = degrees of freedom of the robot, ie, number of joints) and almost the same number of additions to get end-effector's acceleration. Thereofre, if having J and Jdot is not strictly necessary (which seems to be your case, isn't it?), implementing directly a recursive acceleration computation would be faster as you 'traverse' the kinematic chain just once. This is obviously just a rough explanation, but I hope it clarify a little what I said :)
Concluding note: If you do not have particular efficiency needs, using J and Jdot is probably ok ;)
Can you explain to me how (d/dtJ(q)) is not just zero? When we have a function x = f(q). x_dot = df/dqq_dot. x_dd = df/dq/dtq_dot + df/dq*q_dd. But isnt df/dq/dt = 0?