Regarding glm.nb() and my parameter I have been doing a negative binomial regression model using the following code

My my estimate here comes out as 3.48. (the exponential of the intercept).
The data was taken randomly (with set seed) from a distribution with my 3.2 so yeah it makes sense for it to be this close. BUT
When i try and find the maximum likelihood function it appears to be peaking at a my value of around 3.13, which to me doesn't make much sense. Why aren't these to my's the same. I would assume the GLM.NB is using the maximum likelihood estimator on a log-scale? 
Here is the script for reproducing the likelihood function:

 A: You are forgetting that glm.nb has to estimate the size parameter of the negative binomial distribution (called theta in glm.nb) as well as the intercept coefficient.
glm.nb is maximizing the likelihood with respect to two parameters (theta and the intercept coefficient) whereas you are maximizing with respect to only one (the intercept). So naturally glm.nb will find a higher value of the likelihood than you do and a different coefficient value.
When you ran summary(fit.what) the output would have shown you that glm.nb has estimated theta to be 0.3913. If you evaluated the negative binomial likelihood with size = 0.3913 instead of your value then your MLE for the intercept would agree with that from glm.nb.
I wonder whether you may have misunderstood what the size parameter represents for the negative binomial distribution. The variance of the distribution is given by
$${\rm var}(y) = \mu + \mu^2 / \theta$$
where $\mu$ is the mean and $\theta$ is the size parameter. glm.nb assumes a constant value for $\theta$, not one that is proportion to $\mu$, as you have made it in your LL code.
If you try to make $\theta$ depend on $\mu$, then you will destroy the quadratic mean-variance relationship that is a hallmark of the negative binomial distribution.
There is a careful discussion of negative binomial GLMs in Chapter 10 of my book with Peter Dunn (Dunn and Smyth, 2018), including some code examples using the glm.nb function.
Reference
Dunn, PK, and Smyth, GK (2018). Generalized linear models with examples in R. Springer, New York, NY. https://doi.org/10.1007/978-1-4419-0118-7
A: thanks for the answer (Gordon Smyth). Super useful. I can now reproduce the LL using a size parameter of dnbinom set to constant 0.3913, then i get the peak at aronud 3.48. However our assignment was through this expression.

This yields me the likelihood function directly through this script
LL2 <- function(my)
{
R2 = gamma(y+my/7)/(factorial(y)*gamma(my/7)7^(my/7))((7/(1+7))^(y+my/7))
sum(log(R2))
}
Vectorize(LL2)
LL2(3.13425)
par(mfrow=c(1,1))
curve(logLike2,from=2,to=4,xname="my")
I was wondering if there is a way of implementing the constant theta into this equation so it will reproduce the likelihoodfunction that peaks at 3.48 instead of my previous one where the size parameter is proportional to the mu that  peaks at 3.12. I'm not sure how to interpret the size parameter of p(y), since if i set the mu = 3.2 and w = 7, the function will no longer have any mu for the loglikelihoodfunction. 
We were told to write the loglikelihoodfunction of the p(y) given above as a function of mu, given w = 7. Im not sure which way is better to combat this - which LL fucntion would be best to answer this? The one peaking at 3.12 or the one peaking at 3.48?
