Motion described by $m \frac{\mathrm{d}^2 x}{\mathrm{d}t^2}=-k\frac{\mathrm{d}^{\frac12 }x}{\mathrm{d}t^{\frac12}}$
If $D^n$ denotes the $n$th derivative and $D^{-n}$ the $n$th integral, then we have that,
$$D^n f(t) = D^m[D^{-(m-n)}f(t)]$$
providing $m \geq \lceil{n}\rceil$. For our half derivative, we choose $n=1/2$, and $m=2$, in which case we have,
$$D^{1/2}f(t) = D^2[D^{-(3/2)}f(t)]$$
There is a general formula for the $n$th integral of a function, one of my favorite results of Cauchy:
$$f^{-(n)}(t) = \frac{1}{\Gamma(n)}\int_{0}^t (t-u)^{n-1}f(u) \, du$$
which is essentially a convolution $f(t) \ast t^{n-1}$. Applying it, we find,
$$D^{1/2}f(t) = \frac{d^2}{dt^2} \left[ \frac{2}{\sqrt{\pi}}\int_0^t (t-u)^{1/2}f(u) \, du\right]$$
Given the differential equation,
$$\frac{d^2 x(t)}{dt^2} = -\frac{k}{m} \frac{d^{1/2} x(t)}{dt^{1/2}}$$
we can substitute in our definition of $D^{1/2}x(t)$, and conclude,
$$x(t) = -\frac{2k}{m\sqrt{\pi}}\int_{0}^t (t-u)^{1/2}x(u)\, du + c_1t +c_2$$
for $c_1,c_2 \in \mathbb{R}$ which is an integral equation. If we can assume $x(t)$ is supported on $[0,\infty)$ only, then the integral is a convolution $x(t) \ast \sqrt{t}$ and taking the Laplace transform, we find,
$$X(s) = \left( 1+ \frac{k}{ms^{3/2}}\right)^{-1} \left( \frac{c_1}{s^2} + \frac{c_2}{s} \right) = \frac{m(c_1 + c_2 s)}{k\sqrt{s}+ms^2}$$
The solution $x(t)$ is then the inverse Laplace transform of $X(s)$. Formaly, this is given by,
$$x(t) = \frac{1}{2\pi i} \int_{\Gamma} e^{st} \frac{m(c_1 + c_2 s)}{k\sqrt{s}+ms^2} \, ds$$
where the contour $\Gamma$ is in the complex plane; it is a vertical line of infinite length with all poles of $F(s)$ to its left. In practice, we close the contour with an additional contour, ensure the second integral tends to zero (e.g. by the estimation lemma), and use the residue theorem.
The integrand, which we denote $F(s)$, has three poles located at $s^3 = k^2/m^2$, or equivalently,
$$s_1 = \omega^{4/3}_0, \quad s_2 = \frac{1}{2}(1+i\sqrt{3})\omega^{4/3}_0, \quad s_3 = \frac{1}{2}(i\sqrt{3}-1)\omega^{4/3}_0$$
as well as at $s_0= 0$, where we define $\omega^2_0 := k/m$. The vertical contour should begin after $s_1$ so all poles are to the left. However, doing so analytically is somewhat tedious. I chose to use a numerical method for the evaluation of inverse Laplace transforms due to H.E. Salzer which uses aquadrature formula. With Mathematica, I managed to reconstruct $x(t)$ partially:
in the simplified case when $c_1 = c_2 = k/m = 1$. It seems, by visual inspection, the solution resembles that of damped harmonic motion, such as when one introduces a damping term $\gamma \dot{x}$ in the equations of motion of a standard harmonic oscillator.
I am no mathematician and am a little afraid that my answer is too simple to be true, but here goes:
I use Fourier transforms to define the fractional derivative. $x(\omega)$ is defined such that
$$ x(t) = \int_{-\infty}^\infty \, \frac{\text{d}\omega}{2\pi} \text{e}^{i \omega t} \, x(\omega) \, .$$
Then any integer derivatives is
$$ \frac{\text{d}^n}{\text{d} t^n} x(t) = \int_{-\infty}^\infty \, \frac{\text{d}\omega}{2\pi} \text{e}^{i \omega t} \, (i \omega )^n \, x(\omega) \, ,$$
so that
$$ \left(\frac{\text{d}^n x}{\text{d} t^n}\right)(\omega) = (i \omega )^n \, x(\omega) \, .$$
Then this can be generalised to any real number $n$. In particular,
$$ \frac{\text{d}^{1/2}}{\text{d} t^{1/2}} x(t) = \int_{-\infty}^\infty \, \frac{\text{d}\omega}{2\pi} \text{e}^{i \omega t} \, \sqrt{i \omega } \, x(\omega) \, .$$
Since your equation is linear it can be solved separately for each Fourier modes. In Fourier space it becomes
$$ \left(- m \omega^2 + k \sqrt{i \omega}\right) x(\omega) = 0 \, .$$
This tells us that either $x(\omega) = 0$ or $ \omega^2 = \frac{k}{m} \sqrt{i \omega}$. The first solution is trivial. It is equivalent to $x(t)= 0$. The second however has four solutions:
\begin{align} & \omega_1 = \left(\frac{k}{m}\right)^{2/3} \text{e}^{i\,\pi/6} \, ,\\ & \omega_2 = \left(\frac{k}{m}\right)^{2/3} \text{e}^{i\,\pi/2} \, , \\ & \omega_3 = \left(\frac{k}{m}\right)^{2/3} \text{e}^{i \, 5\pi/6} \, , \\ & \omega_4 = 0 \, . \end{align}
Then the most general solution to your problem is
$$x(t) = C_1 \text{e}^{i\omega_1 t} + C_2 \text{e}^{i\omega_2 t} + C_3 \text{e}^{i\omega_3 t} + C_4\, .$$
$C_i$ are integration constants. Explicitly one finds,
$$x(t) = C_1 \, \text{e}^{i \,t \, (k/m)^{2/3} \sqrt{3}/2} \, \text{e}^{-(k/m)^{2/3} t/2} + C_2 \, \text{e}^{- (k/m)^{2/3} t} + C_3 \, \text{e}^{-i \, t \, (k/m)^{2/3} \sqrt{3}/2} \, \text{e}^{-(k/m)^{2/3} t/2} + C_4 \, .$$
We find three solutions that decay exponentially and one constant. Two of the decaying solutions oscillate as well.
In order to make a plot, I set C_4 = 0 because this amounts to a simple vertical shift of $x(t)$, I choose $C_1 = C_2^* = (c_1 + i c_2)/2$ and $C_3 = C_3^*$ so that $x(t)$ is real, I rescale the time according to $\tau = t \, (k/m)^{2/3}/2$ and rescale $x(t)$ according to $y = x \, C_3$. Then my solution becomes,
$$y(\tau) = \text{e}^{-\tau} \left[ c_1 \, \cos(\sqrt{3} \, \tau) + c_2 \sin(\sqrt{3} \, \tau) \right] + \text{e}^{-2 \tau} \, .$$
Here is a plot of $y(\tau)$ for different values of $c_{1,2}$:
You can simply take the semi-derivative of your equation again, which yields
$$\begin{align} m\frac{d^2}{dt^2}\underbrace{\frac{d^{\tfrac12}x}{dt^{\tfrac12}}}_{=-\frac mk\frac{d^2x}{dt^2}} &= -k\frac{dx}{dt} \\\Rightarrow m^2\frac{d^4x}{dt^4} &= k^2\frac{dx}{dt} \tag{*} \end{align}$$
and then solve that ODE. But, similarly to squaring an algebraic equation to eliminate roots, you (probably) have to discard halve of the solutions to satisfy the original equation.
As hinted at by user121330's comment, I assumed commutativity of $D^2$ (where $D:=d/dt$) and $D^{\tfrac12}$ which is not granted in general. I suspect that basically correlates with my above statement that some solutions of $(*)$ must actually be discarded. I'll try and think about this some more, but for now please be aware that there may be flaws or even uselessness in this answer...