In the OP's particular case, the situation is somehwat simpler than the general case that José discusses. That's because the family of left-invariant metrics on $\mathrm{SU}(4)$ that the OP wants to consider has special properties, although just how special does not become apparent until one looks at the problem from a rather different viewpoint, using the fact that $\mathrm{SU}(4)$ is $\mathrm{Spin}(6)$. (In fact, one has $\mathrm{SU}(4)/\{\pm I_4\}=\mathrm{SO}(6)$, and the problem is much easier to describe and treat as a problem on $\mathrm{SO}(6)$, as will be seen.)

First, though, a quick review of the geodesic equations for a left-invariant metric on a compact, semi-simple Lie group $G$: If $\kappa:{\frak{g}}\times{\frak{g}}\to\mathbb{R}$ is the Killing form on ${\frak{g}} = T_eG$, and $\omega:TG\to{\frak{g}}$ is the canonical left-invariant form, then the standard bi-invariant metric on $G$ is given by $\mathrm{d}s^2 = -\kappa(\omega,\omega)$. Any other left-invariant metric on $G$ can be written uniquely in the form $\mathrm{d}\bar s^2 = -\kappa(B\omega,\omega)$, where $B:{\frak{g}}\to{\frak{g}}$ is a positive definite $\kappa$-symmetric linear isomorphism. To find the $\mathrm{d}\bar s^2$-geodesic passing through $g_0\in G$ with initial velocity $L'_{g_0}(v_0)\in T_{g_0}G = L'_{g_0}\bigl({\frak{g}}\bigr)$, one has a $2$-step procedure: First, one finds the curve $v:\mathbb{R}\to{\frak{g}}$ that satisfies the *Euler equation* (a nonlinear ODE initial value problem)
$$
v'(t) = B^{-1}\bigl[v(t),Bv(t)\bigr],\qquad v(0) = v_0
$$
and then the curve $g:\mathbb{R}\to G$ satisfying the *Lie equation*
$$
\omega\bigl(g'(t)\bigr) = v(t),\qquad g(0) = g_0\,.
$$
(When $G$ is a matrix group, this latter equation is just $g'(t) = g(t) v(t)$, with initial value $g(0) = g_0$.)

Note that, when $v_0$ is an eigenvector of $B$, the solution of the Euler equation is $v(t) = v_0$, and so the geodesic is just $g(t) = g_0 \exp(tv_0)$ (i.e., the left-translation of a $1$-parameter subgroup). More generally, if $B$ preserves a subalgebra ${\frak{s}}\subset {\frak{g}}$ that contains $v_0$, then the problem reduces to finding the geodesic in the corresponding subgroup $S\subset G$ (which is totally geodesic in $G$ with respect to the metric $\mathrm{d}\bar s^2$).

Next, in the OP's specific case, one has ${\frak{g}} = {\frak{su}}(4) = {\frak{so}}(6)$ and the OP has prescribed an orthogonal basis $\mathbf{b}$ consisting of 15 elements in ${\frak{su}}(4)$ and wants to consider, all together, the $15$-dimensional cone of metrics determined by the set of positive definite symmetric transformations $B:{\frak{su}}(4)\to {\frak{su}}(4)$ that preserve the $15$ lines spanned by the elements of $\mathbf{b}$. What is not apparent in the OP's description is the great deal of symmetry that the basis $\mathbf{b}$ possesses.

This is much more apparent when one, instead, uses the alternative form ${\frak{so}}(6)$, i.e., the skew-symmetric linear transformations of $\mathbb{R}^6$ with its standard inner product. In this form, one can describe the OP's basis $\mathbf{b}$ as follows: Let $e_1,\dots,e_6$ be an orthonormal basis of $\mathbb{R}^6$ and let $E_{ij}\in {\frak{so}}(6)$ for $1\le i<j\le 6$ be the rank $2$ linear transformation that satisfies $E_{ij}(e_i) = e_j$ and $E_{ij}(e_j) = - e_i$. Then the basis $\mathbf{b} = \bigl(E_{ij}\bigr)_{i<j}$ is orthonormal with respect to the Killing form of ${\frak{so}}(6)$, and it corresponds, under an appropriate isomorphism, to the OP's prescribed basis of ${\frak{su}}(4)$, at least up to signs (which are immaterial to the problem). (Verifying this is an interesting exercise for the reader.)

Now, a subalgebra ${\frak{s}}\subset {\frak{so}}(6)$ is invariant under all of the positive definite linear transformations $B:{\frak{so}}(6)\to {\frak{so}}(6)$ that preserve $\mathbf{b}$ up to multiples if and only if it has a basis that is a subset of $\mathbf{b}$. There are many such subspaces, and this makes it easy to compute the geodesics for many initial values $v_0$:

There are $15$ such maximal tori ${\frak{t}}\subset {\frak{so}}(6)$: For any permutation $\pi = \bigl(\pi(1),\ldots,\pi(6)\bigr)$ let ${\frak{t}}_\pi$ be spanned by the three elements $E_{\pi(1)\pi(2)}$, $E_{\pi(3)\pi(4)}$, and $E_{\pi(5)\pi(6)}$. Then, for $v_0\in {\frak{t}}_\pi$, the solution to the Euler equation is $v(t) = v_0$, so the corresponding geodesics for all of the $15$-parameter family of left-invariant metrics are left-translates of $1$-parameter subgroups.

There are $20$ such copies of ${\frak{so}}(3)\subset {\frak{so}}(6)$: For any triple $(i,j,k)$ with $1\le i<j<k\le 6$, let ${\frak{so}}(3)_{ijk}$ be spanned by the elements
$E_{ij}$, $E_{ik}$, and $E_{jk}$. Then this defines a subgroup of $\mathrm{SO}(6)$ that is totally geodesic for all of the metrics in the OP's class, and each of these metrics restricts to be a left-invariant metric on each such $\mathrm{SO}(3)$. (Unfortunately, these include the general left-invariant metrics on $\mathrm{SO}(3)$, and, as is well-known, the geodesic equations for the generic such metric on $\mathrm{SO}(3)$ can only be integrated using the Jacobian elliptic functions. [See any good book on mechanics for this integration, where it is described as solving the rigid body problem. Also, note the onset of chaos already in this simple case.] As a result, it follows that it is hopeless to expect a general solution in any explicit form, even for the Euler equation.) Note, by the way, that these $20$ copies of totally geodesic $\mathrm{SO}(3)$s in $\mathrm{SO}(6)$ can be grouped into $10$ pairs that commute with each other, which generates $10$ totally geodesic copies of $\mathrm{SO}(3)\times\mathrm{SO}(3)$ on which the geodesic equations for all the metrics in the family reduce to solving independent pairs of $3$-dimensional rigid body problems.

There are, of course, other subgroups that are totally geodesic for the entire $15$-dimensional cone of metrics, such as $15$ copies of $\mathrm{SO}(2)\times\mathrm{SO}(4)$, and $6$ copies of $\mathrm{SO}(5)$. But the Euler equations become progressively harder to solve, and, as far as I know, there is no general solution known for this family of left-invariant metrics on $\mathrm{SO}(5)$ and maybe not even for $\mathrm{SO}(4)$. (Even the $\mathrm{SO}(2)\times\mathrm{SO}(4)$ case is not easy: Even though the Lie algebra of $\mathrm{SO}(4)$ splits as the direct sum of two subalgebras, this splitting is not preserved by the generic linear transformation $B$ in the $15$-dimensional family, and, as a result, the Euler equations do not usually uncouple to simpler equations.)

My conclusion is that, while one can compute the geodesics for this family of left-invariant metrics on $\mathrm{SO}(6)$ for special subspaces of initial conditions for the Euler equations, to get the general solution in any explicit form is probably not possible.