What (discrete) distribution might fit best here? I have some count data I need to fit to a distribution.  I have tried several distributions but none of them start out steep enough.
Here is the fit of a geometric distribution:

My question is: Which distribution should I fit here?  (And (how) is that possible in R?)
In case you want to play:
Due to its size I put the data on pastebin.  You can read it into R using the pastebin package:
## library("devtools")
## devtools::install_github("hrbrmstr/pastebin")
library("pastebin")
ttt <- get_paste("hkuUg0m3")$text
ttt <- as.numeric(strsplit(ttt, "\r\n")[[1]])

I fit the geometric distribution using the fitdistrplus package:
library("fitdistrplus")
fit <- fitdist(ttt, distr = "geom")

And this is what I use for visualization:
viz_fitdist(fit, distr = "geom")

using this function:
viz_fitdist <- function(fit, distr="geom", annotation_colour="#326a97")
{
  library("ggplot2")
  library("grid")
  library("gridExtra")

  data <- data.frame(dat=fit$data)

  plotdata <- as.data.frame(table(data))
  plotdata <- as.data.frame(table(data)/nrow(data))
  plotdata$data <- as.numeric(as.character(plotdata$data))
  p <- plotdata %>% ggplot(aes(x=data, y=Freq)) + geom_bar(stat="identity")

  if (distr=="geom") {
    xsupport <- 1:max(data$dat)
    } else {
      xsupport <- unique(data$dat)
    }
  if (distr!="none") {
    dfun <- paste0("d", distr)
    y <- do.call(dfun, c(list(x = xsupport), as.list(fit$estimate)))
  }

  if (distr!="none") {
    fitdata <- data.frame(x=xsupport, y=y)

    p <- p +
      geom_point(aes(x=x, y=y), data=fitdata, colour=annotation_colour)

    p <- p + geom_line(aes(x=x, y=y), data=fitdata, colour=annotation_colour, size=1)

    p <- p +
      annotation_custom(grob=textGrob(sprintf("Parameter: %.2f", fit$estimate),
                                      x=1-0.07,
                                      y=1-0.07,
                                      just=c("right", "top")))
  }

  p + scale_x_continuous("Data", breaks=xsupport) + ylab("Frequency")
}

EDIT: more distributions
Here is the fit of a negative binomial:
ttt %>%
  fitdist(distr = "nbinom") %>%
  viz_fitdist(distr = "nbinom")


And a zero inflated poisson:
library("extraDistr")
ttt %>%
  fitdist(distr = "zip", start = list(lambda = 1.6, pi = 0.5), lower = c(-Inf, 0), upper = c(Inf, 1)) %>%
  viz_fitdist(distr = "zip")


And a Weibull:
ttt %>%
  fitdist(distr = "weibull") %>%
  viz_fitdist(distr = "weibull")


 A: Let me open by saying that searching for and identifying a model that appears to describe the data is not where I would normally begin my modelling; it begins with the information whuber and I were seeking on how the model will be used, with what is being measured - and further, to consider any associated theoretical or expert practical understanding of its properties, the performance of models on similar data sets, and so on. Considerations like these are essential to the choice of an appropriate model. 
Nevertheless, as requested, I will attempt to describe some steps one could take in examining data like this and seeing if it could suggest a model.
If you're going to use the resulting distributional model for inference (hypothesis tests, or confidence intervals, or prediction intervals for example), there's a problem that you chose the model based on your data so the nominal properties of the inference will no longer hold. (If you really can't use the characteristics of the variable and subject area knowledge, having no other basis than the sample itself, then it's better to pull off a random subset and use that for model identification.)
Given the heavy tail, one option is to look at the way the log of the empirical survivor function relates to $x$. My first plots were to plot $\log(\hat{S}(x))$ against $\log(x)$ (not shown; linear at the top but then somewhat curved near the middle) and then against $\sqrt x$. The appearance of these then suggested to me that a fair first model might be a zeta distribution. However, they can also suggest a number of other possibilities (such as perhaps directly modelling the seeming linear relationship in the first plot below; it might provide a somewhat better fit, at the expense of being less familiar)
So I then looked at the empirical pmf, plotting $\log(\hat{p}(x))$ against $\log{x}$ as a check on the potential for fitting a zeta function. 

That indeed suggests a fairly straight line relationship, and it looks like a parameter in the ballpark of 2 or perhaps a bit more would be a reasonable first model -- note that the most precisely estimated proportions are the first few, and fitting those accurately will matter more for a lot of applications -- but not all. However, it's not a perfect fit - with a large sample size like you have, I believe you would easily reject a pure zeta model (or indeed any other simple model).

(NB the parameter value here has not been optimized; it's just a reasonable guess at a suitable value)
So if I had to suggest a simple model on the basis of only looking at the data, the zeta distribution would be an obvious candidate; if an initial guess at the parameter value is needed, I'd start at 2.
On fitting power laws, such as the zeta distribution, also see Cosma Shalizi's essay So you think you have a power law? which debunks some of the enthusiasm for power laws (a lot of things will look more or less power-law like) and the paper it discusses, which has a lot of useful information.
A: Have you tried a Negative Binomial distribution or a zero-inflated Poisson? 
 These distributions are often used to take care of overdispersion.
You have been suggested below a Weibull distribution, which may fit the data but it is not appropriate for count data.  A Weibull is a distribution for continuous measurements, and not for discrete data.  Negative Binomial distributions are routinely used for RNAseq and Microbiome data, among other types of data. So are the zero-inflated Poisson distribution. They are standard tools and are easily implemented in R packages  See Cameron-Trivedi (2013, Regression Analysis of Count Data). Now, there are also Discrete Weibull distributions, but they are not so much used. See the R package "DiscreteWeibull" and related manuscripts. 
The negative binomial fit above is similar to the Weibull and it is more appropriate, since it takes into account the nature of the data. 
If you really wanna be flexible, mixtures of Poisson-Gamma may help you to capture overdispersion in a Bayesian fit. They are of course not trivial. But if the question is about fitting the data, then all the options should be proposed. 
Now, your data show that histogram because there are probably bumps at the right end of the x-axis. No parametric model can capture those bumps correctly. That's why I was suggesting mixtures of parametric models. And in such case, the Bayesian fit is probably the way to go.
The right model depends of course on the objective of your inference, which is not clear from your question. 
