as suggested in a previous answer, numerical integration is a solution in this case. Below you will find 2 ways of achieving this using R. Hope this helps.
set.seed(1)
x1 = rnorm(250)
#Computing of bandwidth
h1 = bw.ucv(x1)
#Estimate
est = density(x1,kernel = "gaussian",bw = h1)
# Create approx func obj and integrate()
splxy = splinefun(est$x, (est$y - dnorm(est$x))^2)
integrate(splxy, lower = min(x1), upper = max(x1)) # res 0.0006681713
# Simple numerical integration in pkg sfsminsc
sfsmisc::integrate.xy(x = est$x, (est$y - dnorm(est$x))^2) # res 0.0006682952