Does the matrix square root have directional derivatives at semipositive points? $\newcommand{\psym}{\operatorname{P}_{\ge 0}}$
Let $\psym$ be the set of symmetric positive semidefinite (real) matrices.
Let $\sqrt \cdot :\psym  \to \psym $ be the unique positive semidefinite square root.
Let $A \in \psym$ be a matrix on the boundary, i.e $\det A=0$. Let $B$ be a symmetric matrix such that $A+tB$ is positive semidefinite for $t>0$ small enough**. 
Quesion: Characterize all such pairs $A,B$ where the one sided directional derivative
$$  (d\sqrt{\cdot})_A(B)=\left. \frac{d}{dt}\right|_{t=0}  \sqrt{A+tB}:=\lim_{t \to 0^+} \frac{ \sqrt{A+tB}-\sqrt{A}}{t} \, \, \text{ exist.}$$

Partial results:
(1) For $A=0$ the limit exists only for $B=0$.
(2) For every $A=B$ the limit exists (and equals $\frac{1}{2}\sqrt{A}$).
(3) For every $A$, there are a "lot" of matrices $B$ for which the limit does not exist (see the answer of loup blanc):  
By orthogonal diagonalization, after noting that 
$(d\sqrt{\cdot})_A(B)$ exists $\iff (d\sqrt{\cdot})_{V^TAV}(V^TBV)$ exists for every $V \in O_n$, 
we reduce to the case $A=diag(0,a_2,\cdots,a_n)$ where $a_i\geq 0$: 
For every matrix $B$ of the form $B=diag(\lambda,C)$ where $\lambda > 0,C\in M_{n-1}$ is symmetric s.t. $diag(a_2,\cdots,a_n)+C\geq 0$, 
$(d\sqrt{\cdot})_A(B)$ does not exist.

In this question, it is shown that $\sqrt{\cdot} $ is not differentiable in the standard sense. In fact, this answer shows $\psym$ is not a manifold (with boundary) at all, so "standard differentiability" does not really make sense here.

**
Actually, as noted by loup blanc, if $B$ is a symmetric matrix such that $A+B \in \psym$, then $A+tB \in \psym$ for small enough $t$. Indeed,
$A,A+B \in \psym$ and $\psym$ is convex, hence $(1-t)A+t(A+B)=A+tB \in \psym$.
Note that the reverse implication is false in general; It can happen that $A+tB \in \psym$, but $A+B \notin \psym$, see here.
 A: ** is true because $P_{\geq 0}$ is convex (as the closure of a convex).
Now, if $A$ is in the boundary, then, up to an orthonormal change of basis, $A=diag(0,a_2,\cdots,a_n)$ where $a_i\geq 0$. For non-differentiability, it suffices to choose $B=diag(1,C)$ where $C\in M_{n-1}$ is symmetric s.t. $diag(a_2,\cdots,a_n)+C\geq 0$.
A: This is not really answer, but rather a long comment:$\newcommand{\psym}{\operatorname{P}_{\ge 0}}$
When:
(1) $B$ is invertible
(2) $A,B$ commute
(3) $B^{-1}A \in \psym$
it is possible to reduce the problem to the case of $B=I$:
$$ \sqrt{A+tB}-\sqrt{A}=\sqrt{B(B^{-1}A+tI)} -\sqrt{BB^{-1}A} \stackrel{(*)}{=} \sqrt{B}\sqrt{B^{-1}A+tI}-\sqrt{B}\sqrt{B^{-1}A}=$$
$$ \sqrt{B} \cdot (\sqrt{B^{-1}A+tI}-\sqrt{B^{-1}A}) $$
where equality (*) comes from lemma 2 below. Note that $\sqrt{B^{-1}A+tI},\sqrt{B^{-1}A}$ are well defined by assumption (3).
To conclude, $$ (d\sqrt{\cdot})_A(B)=\sqrt{B}\cdot (d\sqrt{\cdot})_{B^{-1}A}(I)$$
when this equality is interpreted in the sense that the left hand exists $\iff$ the  right hand exists, and in that case they are equal.

Lemma(1):
Let $A,B \in \psym$ be commuting matrices. Then $\sqrt{A},\sqrt{B}$ commute.
Proof:
See here.
Lemma(2):
Let $A,B \in \psym$ be commuting matrices. Then $\sqrt{AB}=\sqrt{A}\sqrt{B}$.
Proof:
By lemma(1) $\sqrt{A}\sqrt{B}$ is symmetric. In fact $\sqrt{A}\sqrt{B} \in \psym$ (since the product is symmetric, see here), and its square is $AB$:
$$(\sqrt{A}\sqrt{B})^2=\sqrt{A}\sqrt{B}\sqrt{A}\sqrt{B}=\sqrt{A}\sqrt{A}\sqrt{B}\sqrt{B} =AB$$ (where the last equality is by lemma(1))
Finally, by the uniqueness of the symmetric semipositive square root, we get $$\sqrt{AB}=\sqrt{A}\sqrt{B} $$.
