Solving the integral of the equation of a mixed model in R

I downloaded your dataset as Data.csv. I had to do some formatting to get it to work on my local machine:

library(ggplot2)
library(nlme)
library(data.table)


##################
# Format data   ##
##################


dat <- read.table("Data.csv", 
                  sep=";", 
                  dec=",", 
                  colClasses=c("character",
                                rep("numeric",4)), 
                  skip=1)
setDT(dat)
format(dat,decimal.mark=".")

dat[, Var2 := V1]
dat[, Var3 := as.numeric(V2)]
dat[, Var4 := as.numeric(V3)]
dat[, Var1 := as.numeric(V4)]
dat[, Var5 := as.numeric(V5)]
dat

## this is name used in OP code
dados <- copy(dat[,c("Var2","Var3","Var4","Var1","Var5")])

I rewrote the code a little bit so I could reproduce your graphic -- Here's your code with some minor formatting changes:

################### BEGIN OP CODE ####################

fitmixedmodel <-  lme( log(Var1) ~  I(exp(Var3/Var4))+ I((Var5/Var4)^3),
                      random = ~1|Var2, 
                      data = dados, 
                      method="REML")

summary(fitmixedmodel)

volume <- dados[dados$Var5 == 0.1,]

fmixedmodel <- function(Var3, Var5, Var4){
  (pi/40000)*
  (Var3^2)*
  (coefficients(summary(fitmixedmodel))[1] + 
   coefficients(summary(fitmixedmodel))[2]*I(exp(Var3/Var4)) + 
   coefficients(summary(fitmixedmodel))[3]*(I((Var5/Var4)^3)))
}

vmixedmodel <- function(Var3, Var5, Var4){ 
  integrate(Vectorize(fmixedmodel), 
            lower = 0.1, 
            upper = Var4, 
            Var3 = Var3,  
            Var4 = Var4)$value  
}

mixed.vol <- mapply(FUN = vmixedmodel,  
                    Var5 = as.list(volume$Var5),
                    Var3 = as.list(volume$Var3),
                    Var4 = as.list(volume$Var4))

################# END OP CODE ##################

## now verify the graph.  looks good.
ggplot() +
  geom_point(aes(y=mixed.vol, x=volume$Var3, color=volume$Var3))

So at this point I was able to reproduce your graphic.

I initially thought there were two options for incorporating random intercepts. One is to "integrate them out" which would involve the variance of the random intercepts and a double integral. But it turns out for linear regression this type of marginalization doesn't change the outcome. To prove this to ourselves, look at the following code that goes through the trouble of fitting a double integral to integrate out the random intercept b that follows a Normal(0, 0.1691067^2) distribution. Because the integral with respect to b can just isolate b by itself and the E[b] = 0, this approach is not materially different than the OP approach.

# Option 1: integrate over the random intercept distribution
# this will require the random intercept variance as well as 
# double integration.

## to be able to accommodate a random intercept, we need to integrate
## over the random intercepts, which are distributed as N(0, sig2)
## where sig2 is 0.1691067^2 as seen from the fitmixedmodel output:

# 
# Random effects:
#   Formula: ~1 | Var2
#          (Intercept)  Residual
# StdDev:   0.1691067 0.2559742

## add "b" random intercept, multiply whole thing by normal density dnorm
integrand <- function(x,  Var3, Var4){
  
  Var5 <- x[1]
  b <- x[2]
  
  (pi/40000)*(Var3^2)*
  (coefficients(summary(fitmixedmodel))[1] + b  + 
   coefficients(summary(fitmixedmodel))[2]*I(exp(Var3/Var4)) + 
   coefficients(summary(fitmixedmodel))[3]*(I((Var5/Var4)^3))) *
    dnorm(b, sd = 0.1691067)
}

vmixedmodel.option1 <- function(Var5, Var3, Var4){ 
  pcubature(integrand, 
            lower = c(0.1,-Inf), 
             upper = c(Var4,Inf), 
             Var3 = Var3, 
             Var4 = Var4)$integral  
}

## this is slow. And unnecessary.  Because the E[b] = 0
mixed.vol.option1 <- mapply(FUN = vmixedmodel.option1,  
                            Var5 = as.list(volume$Var5),
                            Var3 = as.list(volume$Var3),
                            Var4 = as.list(volume$Var4))

max(abs(mixed.vol - mixed.vol.option1))

ggplot() +
  geom_point(aes(y=mixed.vol.option1, x=volume$Var3, color=volume$Var3))

The second approach is plugging in the estimated random intercept value, much like how the OP approach plugs in values for Var4 and Var3. To pursue this avenue, we first create volume_ri which is the same as the volume dataset but has the estimate values of b:

## Option 2: plug in the random intercept value.

rand_int <- data.table(Var2 = rownames(fitmixedmodel$coeff$random$Var2), 
                          b = fitmixedmodel$coeff$random$Var2 )
setnames(rand_int, names(rand_int), c("Var2","b"))
rand_int

## merge into `volume` (or `dados` and then re-subset)
volume_ri <- merge(volume,
                   rand_int)

And then basically we adust the OP code to accommodate this b as an argument or value where appropriate:

## throw in a b argument
fmixedmodel_ri <- function(Var3, Var5, Var4, b){
  (pi/40000)*
    (Var3^2)*
    (coefficients(summary(fitmixedmodel))[1] + b + 
       coefficients(summary(fitmixedmodel))[2]*I(exp(Var3/Var4)) + 
       coefficients(summary(fitmixedmodel))[3]*(I((Var5/Var4)^3)))
}

## throw in a b argument
vmixedmodel_ri <- function(Var3, Var5, Var4, b){ 
  integrate(Vectorize(fmixedmodel_ri), 
            lower = 0.1, 
            upper = Var4, 
            Var3 = Var3,  
            Var4 = Var4,
               b = b)$value  
}

## plug in the b values
mixed.vol_ri <- mapply(FUN = vmixedmodel_ri,  
                      Var5 = as.list(volume_ri$Var5),
                      Var3 = as.list(volume_ri$Var3),
                      Var4 = as.list(volume_ri$Var4),
                         b = as.list(volume_ri$b))

## now verify the graph. only 8 levels of Var2, so use color
ggplot() +
  geom_point(aes(y=mixed.vol_ri, x=volume_ri$Var3, color=volume_ri$Var2))

Below are old questions that are answered in comments

Old questions:

My concern/confusion is the line where I comment DID YOU MEAN Var5 = Var5 -- would you mind double checking that. and leaving me a comment with the answer?

Also, a separate question regarding accounting for the random intercept:

do you want to integrate over all random intercept values

or

do you want to plug in the random intercept estimate for each unique Var2 from the mixed model fit?