psi[x_]= Exp[-q x^2]
ke=2 Integrate[psi'[x]^2,{x,0,Infinity}]
pe=Integrate[psi[x]^2 (1-Exp[-x])^2,{x,-Infinity,Infinity}]
n=2 Integrate[psi[x]^2,{x,0,Infinity}]
Simplify[(ke/alpha^2+pe)/n,Re[q] > 0]
% /. alpha->10
Plot[%,{q,0,10}]
FindMinimum[%%,{q,5}]
 {0.104178, {q -> 5.40905}}

compare exact: 1/alpha-.25/alpha^2 /. alpha->10
Out[10]= 0.0975

-----

u[n_,l_,r_]=1/(n Sqrt[Pochhammer[(n-l), 2 l +1]]) (2 r/n)^(l+1) LaguerreL[n-l-1,2 l +1,2 r/n] Exp[-r/n]

psi1s=u[1,0,r] SphericalHarmonicY[0,0,theta,phi]
psi2p=u[2,1,r] SphericalHarmonicY[1,0,theta,phi]
Integrate[Integrate[Integrate[psi1s r Sin[theta]Cos[phi] psi2p Sin[theta],{r,0,Infinity}],{theta,0,Pi}],{phi,0,2 Pi}]
Integrate[Integrate[Integrate[psi1s r Sin[theta]Sin[phi] psi2p Sin[theta],{r,0,Infinity}],{theta,0,Pi}],{phi,0,2 Pi}]
Integrate[Integrate[Integrate[psi1s r Cos[theta] psi2p Sin[theta],{r,0,Infinity}],{theta,0,Pi}],{phi,0,2 Pi}]
        128 Sqrt[2]
Out[6]= -----------
           243

Alternatively do all three angular integrals first:
Integrate[Integrate[{Sin[theta]Cos[phi],Sin[theta]Sin[phi],Cos[theta]}
SphericalHarmonicY[0,0,theta,phi] SphericalHarmonicY[1,0,theta,phi] Sin[theta],{theta,0,Pi}],{phi,0,2 Pi}]
                  1
Out[7]= {0, 0, -------}
               Sqrt[3]

Then do the radial integral: (Note: no r^2 dr as these are u rather than R wavefunctions)
Integrate[u[1,0,r] r u[2,1,r] ,{r,0,Infinity}]
                 2
        128 Sqrt[-]
                 3
Out[8]= -----------
            81

and combine the results:
Last[%7]*%8
        128 Sqrt[2]
Out[9]= -----------
            243
N[%]
Out[10]= 0.744936


what is omega?
Ry*(1-1/4)/hbar
1.5503e16

A=1.5503e16*(1.5503e16*e*0.744936*a)^2/(3*pi*epsilon0*c^3*hbar)
6.2683e8
1/A=1.5953e-9 sec


-------------
Note: Mathematica can learn values of physical constants:

<<PhysicalConstants`
<<Units`

w=SI[Rydberg](1-1/4)/PlanckConstantReduced
                  16
        1.55032 10
Out[3]= ------------
           Second

A=w^3(0.744936*BohrRadius*ElectronCharge)^2/(3 Pi VacuumPermittivity SpeedOfLight^3 PlanckConstantReduced)

                  8        2
        6.26856 10  Coulomb  Volt
Out[4]= -------------------------
                             2
          Ampere Joule Second

however mathematica apparently does not know Coulomb*Volt=Joule, Ampere=Coulomb/Second

%4 /. {Volt->Joule/Coulomb, Ampere->Coulomb/Second}
                   8
         6.26856 10
Out[5]= -----------
           Second

