Indeed, this problem should be solved with Fourier cosine transform rather than Fourier transform. If you use Fourier transform, you won't be able to take the b.c. at $x=0$ into consideration, and the equation will then be solved in $-\infty<x<\infty$ with implicit b.c.s at $\pm\infty$. (For this part check this post for more information. )
The key point is, when $f(\infty)=0$, Fourier cosine transform
$$
\mathcal{F}_t^{(c)}[f(t)](\omega)=\sqrt{\frac{2}{\pi }}\int _0^{\infty } f(t) \cos (\omega t) d t
$$
has the following property:
$$
\mathcal{F}_t^{(c)}\left[f''(t)\right](\omega)=-\omega^2 \mathcal{F}_t^{(c)}[f(t)](\omega)-\sqrt{\frac{2}{\pi }} f'(0)
$$
You can easily verify this property using integration by parts.
So, by applying this property on your equation, we have
$$
\mathcal{F}_x^{(c)}\left[u^{(1,0)}(t,x)\right](\omega )=k \left(-\omega ^2 \mathcal{F}_x^{(c)}[u(t,x)](\omega )-\sqrt{\frac{2}{\pi }} u^{(0,1)}(t,0)\right)-a k \mathcal{F}_x^{(c)}[u(t,x)](\omega )
$$
Substitute the b.c. into the equation, we obtain a simple initial value problem (IVP) of linear ODE:
$$U'(t)=-a k U(t)-k \omega ^2 U(t)$$
$$U(0)=F$$
where $U(t)=\mathcal{F}_x^{(c)}[u(t,x)](\omega )$, $F=\mathcal{F}_x^{(c)}[f(x)](\omega)$.
If you have difficulty in solving the IVP, check the wikipedia page). Anyway, we can easily find its solution:
$$U(t)=F e^{-k t \left(a + \omega ^2\right)}$$
The last step is to transform back with inverse Fourier Cosine transform
$${\mathcal{F}_\omega^{(c)}}^{-1}[F(\omega)](t)=\sqrt{\frac{2}{\pi }} \int_0^{\infty } F(\omega ) \cos (t \omega) \, d\omega $$
and the solution is
$$
u(t,x)=\sqrt{\frac{2}{\pi }} \int_0^{\infty } e^{-k t \left(a + \omega ^2\right)} \mathcal{F}_x^{(c)}[f(x)](\omega ) \cos (\omega t) \, d\omega
$$
Notice this solution is probably the same as the one given by doraemonpaul. I guess he has just chosen a different convention for Fourier parameters.
Long solution
Being the domain infinite in one variable you have to use the Fourier transform or, if you prefer, the Laplace transform. Let me show you how to do it, using Fourier transform.
First off we take the Fourier transform of both sides of the PDE and get
$$
\mathcal F\{u_t\} = \mathcal F\{u_{xx}\}\implies \frac{\partial}{\partial t} \hat{u}(k,t) = -k^2\hat{u}(k,t)
$$
This was done by using the simple property of derivation under Fourier transform (all properties are listed on the linked wikipedia page). The function $\hat u$ is the Fourier transform of $u$. Continuing with the solution, we now have a simpler PDE to which the solution is well known
$$
\frac{\partial}{\partial t} \hat{u}(k,t) = -k^2\hat{u}(k,t) \\
\hat{u}(k,t) = c(k)e^{-k^2t} \tag 1
$$
where $c(k)$ is a function to be found by the initial conditions. If we take te Fourier transform of the initial condition we can find $c(k)$. In fact
$$
\mathcal F\{u(x,0)\} = \mathcal F\{e^{-x^2}\} = \frac{1}{\sqrt{2\pi}}\int_{\mathbb R} e^{-ikx}e^{-x^2}\,\mathbb d x = e^{-\frac{k^2}{4}}
$$
A very well known Fourier transform is the one of the Gaussian function: indeed the transform of a Gaussian is itself a Gaussian (with some numeric factors depending on your definition of Fourier transform). You can evaluate it yourself, but there's a table in the wikipedia page.
Plugging this into the function found in $(1)$ what we get is
$$
\hat{u}(k,0) = c = e^{-\frac{k^2}{4}}
$$
So our solution to the differential equation $(1)$ is
$$
\hat{u}(k,t) = e^{-\frac{k^2}{4}}e^{-k^2 t} = e^{-\left(t+\frac{1}{4}\right)k^2}
$$
To go back to the function we just take the inverse Fourier transform of both sides
$$
\mathcal F^-1\{\hat{u}(k,t)\} = \mathcal F^{-1}\left\{e^{-\left(t+\frac{1}{4}\right)k^2}\right\}\\
u(x,t) = \frac{1}{\sqrt{2\pi}}\int_{\mathbb R} e^{ikx}e^{-\left(t+\frac{1}{4}\right)k^2}\,\mathbb dk
$$
In our case we have that
$$
\underbrace{\frac{c}{\sqrt{2\pi}}\int_{\mathbb R} e^{ikt}e^{-\left(t+\frac{1}{4}\right)k^2}\,\mathbb d k}_{\text{Inverse transform of }e^{-\left(t+\frac{1}{4}\right)k^2}} = \sqrt{\frac{2}{4t+1}}e^{-\frac{x^2}{4t+1}} = u(x,t)
$$
It follows clearly that
$$
u_x(x,t) = -\frac{2x}{4t+1}\sqrt{\frac{2}{4t+1}}e^{-\frac{x^2}{4t+1}}\\[20pt]
\implies u_x(0,t) = 0
$$
So at the end, our solution to the problem is
$$
u(x,t) = \sqrt{\frac{2}{4t+1}}e^{-\frac{x^2}{4t+1}} = \frac{2}{\sqrt{8t+2}}e^{-\frac{x^2}{4t+1}}
$$
the last step is totally not mandatory.
Faster solution
This solution requires you to know the Green's function on the standard heat equation
$$
u_t = Du_{xx}
$$
which is
$$
\mathcal G(x,x') = \frac{1}{\sqrt{4\pi D t}}e^{-\frac{(x-x')^2}{4Dt}} \tag 2
$$
From this you can easily find the solution by integration
$$
u(x,t) = \int_{\mathbb R}\mathcal G(x,x')u(x',0)\,\mathbb d x' \tag 3
$$
Your initial condition is pretty convenient! If we plug it in $(3)$ with $(2)$ what we get is
$$
u(x,t) = \frac{1}{\sqrt{4\pi D t}} \int_{\mathbb R}e^{-\frac{(x-x')^2}{4Dt}}e^{-x'^2}\,\mathbb d x'
$$
which upon closer inspection it is clearly a convolution between two Gaussians! And again, by a very cool property of the Gaussian, the convolution of two Gaussians is itself a Gaussian.
$$
\left( e^{-\frac{x^2}{4Dt}} * e^{-x^2}\right) = \color{red}{\frac{1}{\sqrt{2\pi}}}\frac{\sqrt{\pi}}{\sqrt{\frac{1}{4Dt}+1}} e^{-\frac{x^2}{4Dt +1 }}
$$
the red factor comes from the convolution theorem an the definition of Fourier transform that I'm using. Plugging $D=1$ you would get the solution.
Best Answer
Separation of variables definitely works. $$f(x)g'(t)=\alpha f''(x)g(t)+f(x)g(t)$$ Divide through by $f(x)g(t)$: $$\frac{g'(t)}{g(t)}=\alpha \frac{f''(x)}{f(x)}+1=\text{constant}=-K$$ So you have some differential equations $$g'(t)=-K~g(t)$$ $$ f''(x)=-\frac{K+1}{\alpha} f(x)$$ More work... $$g(t)=b e^{-K t}$$ $$f(x)=a_1 \sin\left(\sqrt{\frac{K+1}{\alpha}}x\right)+a_2 \cos\left(\sqrt{\frac{K+1}{\alpha}}x\right)$$ Meaning our $x$ partial derivative is $$\partial_x u(x,t)=be^{-Kt}\sqrt{\frac{K+1}{\alpha}}\left(a_1 \cos\left(\sqrt{\frac{K+1}{\alpha}}x\right)-a_2 \sin\left(\sqrt{\frac{K+1}{\alpha}}x\right)\right)$$ So our boundary conditions are $$\partial_x u(0,t)=0\implies a_1 \cos\left(\sqrt{\frac{K+1}{\alpha}}\cdot 0\right)-a_2 \sin\left(\sqrt{\frac{K+1}{\alpha}}\cdot 0\right)=0$$ $$\implies a_1=0$$ Let's rename $b\cdot a_2 \to C$ $$u(x,t)=C\cdot e^{-Kt} \cos\left(\sqrt{\frac{K+1}{\alpha}}x\right)$$ Can you perhaps apply the second boundary condition now? EDIT: The second boundary condition yields $$\sin\left(\sqrt{\frac{K+1}{\alpha}}\right)=0\implies \sqrt{\frac{K+1}{\alpha}}=n\pi$$ So we have separation constants $$K_n = \alpha n^2\pi^2 -1$$ Meaning our full solution is an arbitrary linear combination (for choices of $\{C_n\}$ that converge) $$u(x,t)=\sum_{n=0}^\infty C_n e^{-(n^2\pi^2-1)t}\cos\left(n\pi x\right)$$ Meaning $$u(x,0)=\sum_{n=0}^\infty C_n \cos\left(n\pi x\right)$$ But also $$u(x,0)=\sin^2(\pi x)=\frac{1}{2}-\frac{1}{2}\cos(2\pi x)$$ So $C_0=1/2$, $C_2=-1/2$, and all others are $0$.