Does it make any statistical sense to subtract / add RMSE from the predicted results to make the model more accurate? If RMSE represents the standard deviation of the residuals (prediction errors), while  residuals are a measure of how far from the regression line data points are, does it make any sense to subtract / add RMSE to the predicted results as some sort of adjustment to make the model more accurate?
 A: Your question indicates that you are thinking about the sampling variability for the parameter estimates in the model that you would obtain from a sample data set.  And, while we cannot make a single prediction "more accurate", your idea of improving the model may be restated as ¿can we add more confidence to our predictions with the rmse?  (If I have completely misrepresented your query, please let me know.)
The answer is:  yes.  And your intuition that it would involve the rmse is essentially spot on.  But, the manner in which we add/subtract this value is a bit convoluted.
I will provide here the formula you would use for a bivariate regression model (one scalar predictor and one scalar response variable).  If you imagine starting with your bivariate regression model, you can make a point-estimate prediction:
$$\hat{y} = \beta_0 + \beta_1 · x_0$$
Now, if you want to add confidence to this prediction, we will add/subtract the rmse...but we also have to account for some of the variability in our estimates for the slope and intercept.  To update the estimate, we have the following formula:
$$\hat{y} = \left(\beta_0 + \beta_1 · x_0\right) \pm t_\text{c.v.}·\text{rmse}·\sqrt{1 + \frac{1}{n} + \frac{(x_0 - \bar{x})^2}{s_x^2(n-1)}}$$
where $t_\text{c.v.}$ is the critical value associated with the degree of confidence you wish to have in your final estimate, $n$ is the sample size, $\bar{x}$ is the mean of the independent variable, and $s_x$ is the standard deviation of the independent variable.
Note, the 1 under the square root is multiplied by the rmse, and this captures the idea that you were suggesting in your post.
I hope this was useful.
A: Let's simulate this and see what happens.
set.seed(2022)

# Define sample size
#
N <- 100 

# Define predictor variable
#
x <- sort(runif(N, -2, 2))

# Define y
#
y <- 3 - x + rnorm(N, 0, 1) 

# Fit linear model
#
L <- lm(y ~ x) # Fit linear model

# Calculate RMSE
#
rmse <- sqrt(mean((predict(L) - y)^2)) # Calculate RMSE

# # Add/subtract the RMSE to/from the predicted values
#
adjusted_predictions <- predict(L) + sample(c(-1, 1), N, replace = T)*rmse 

# Calculate the RMSE of the adjusted predictions
#
adjusted_rmse <- sqrt(mean((adjusted_predictions - y)^2)) 

rmse # I get ~1.066
adjusted_rmse # I get ~ 1.57, which is greater than the original RMSE

# Let's do it many times

# Define the number of times to try the prediction adjustment
#
R <- 1000 
adjusted_rmses <- rep(NA, R) # Blank vector to hold RMSE values

# Loop
for (i in 1:R){ 
    
    # Add/subtract the RMSE to/from the predicted values
    #
    adjusted_predictions <- predict(L) + sample(c(-1, 1), N, replace = T)*rmse 
    
    # Calculate the RMSE of the adjusted predictions
    #
    adjusted_rmses[i] <- sqrt(mean((adjusted_predictions - y)^2)) 
    
}

mean(adjusted_rmses) # I get ~ 1.50, which is greater than the original RMSE

Your proposed methodology results is weaker performance, not stronger.
I think I see why you would expect it to be better, however. If you know that the average amount by which the model misses the truth is the RMSE (not quite what RMSE means, but it's close), then adjust the prediction, right? The trouble is that, since we don't know the truth when we make predictions for real (if we did, there would not be anything to predict), we have to guess about adding or subtracting the value by which we are adjusting the prediction.
If we guess right, then we improve performance, but if we guess wrong, we lower performance. Interestingly, having equal numbers of right and wrong guesses seems not to cancel out to result in equal performance.
Further, we can overshoot the true value, even if we move in the right direction. For instance, if the prediction is lower than the observation and we guess that we should increase the prediction, we might wind up with a larger residual. As an example, consider an $RMSE$ adjustment of $10$, for a true value of $5$ and a prediction of $4$. We guess that our prediction is low, so we add $10$ to get an adjusted prediction of $14$, but now our residual has a magnitude of $9$ instead of $1$.
NOW LET'S PROVE IT
My simulation assumes a model of $y_i = \beta_0 + \beta_1x +\epsilon_i$ for $iid$ $\epsilon_i$. $$RMSE = \sqrt{\mathbb E\left[\epsilon_i^2\right]} = \sqrt{\mathbb E\left[\left(Y_i - \hat Y_i\right)^2\right]}$$
I am going to drop the square root so I don't have to keep writing it, but if we show that $MSE_{original}<MSE_{adjusted}$, then $RMSE_{original}<RMSE_{adjusted}$.
Let $\tilde Y_i = \hat Y_i + Z_i$ be an adjusted prediction, where $Z_i$ takes $\pm RMSE$ with equal probability. Note that $Z_i$ has positive variance and zero expected value, so its second moment is equal to its variance. $\hat Y_i$ and $Z_i$ are independent, so when we get to the step in the proof where we calculate $\mathbb E\left[\hat Y_iZ\right]$, we can take that as $\mathbb E\left[\hat Y_i\right]\mathbb E\left[Z\right]$.
The $MSE$ of the adjusted predictions is $\mathbb E\left[
\left(
Y_i - \tilde{Y_i}
\right)^2
\right]$.
$$
=
\mathbb E\left[
Y_i^2 - 2Y_i\tilde Y_i + \tilde Y_i^2
\right]\\=
\mathbb E\left[
Y_i^2 -2Y_i(\hat Y_i + Z_i) + (\hat Y_i + Z_i)^2
\right]\\=
\mathbb E\left[
Y_i^2 - 2Y_i\hat Y_i -2Y_iZ_i + \hat Y_i^2 + 2\hat Y_i Z_i + Z_i^2
\right]\\=
\mathbb E\left[
Y_i^2 - 2Y_i\hat Y_i+\hat Y_i^2 -2Y_iZ_i+2\hat Y_iZ_i + Z_i^2
\right]\\=
\mathbb E\left[
(Y_i -\hat Y_i)^2-2Y_iZ_i+2\hat Y_iZ_i + Z_i^2
\right]\\=
\mathbb E\left[
\left(
Y_i - \hat{Y_i}
\right)^2
\right] - 2\mathbb E\left[Y_i Z_i\right] + \mathbb E\left[\hat Y_i Z_i\right] + \mathbb E\left[Z_i^2\right]
$$
Since $Y_i$ is a constant, we can pull it out of the expectation.
$$
=
\mathbb E\left[
\left(
Y_i - \hat{Y_i}
\right)^2
\right] - 2Y_i \mathbb E\left[Z_i\right] + \mathbb E\left[\hat Y_i Z_i\right] + \mathbb E\left[Z_i^2\right]
\\=
\mathbb E\left[
\left(
Y_i - \hat{Y_i}
\right)^2
\right] - 2Y_i \mathbb \times 0 + \mathbb E\left[\hat Y_i\right]\mathbb E\left[Z\right] + \mathbb E\left[Z_i^2\right]
\\=
\mathbb E\left[
\left(
Y_i - \hat{Y_i}
\right)^2
\right] - 2Y_i \mathbb \times 0 + \mathbb E\left[\hat Y_i\right]\times 0 + \mathbb E\left[Z_i^2\right]
\\=
\mathbb E\left[
\left(
Y_i - \hat{Y_i}
\right)^2
\right] + \mathbb E\left[Z_i^2\right] 
\\>\mathbb E\left[
\left(
Y_i - \hat{Y_i}
\right)^2
\right] = MSE_{original}\\\square
$$
