Can Agresti-Coull binomial confidence intervals be negative? According to http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm
Agresti-Coull intervals cannot be negative; however using the formula from Wikipedia as well as the binom.confint function in R with x = 0 and n = 5 produces an interval of (-0.05457239, 0.4890549).
Are there conflicting definitions of Agresti-Coull?
Edit - 
Using binom.confint:
binom.confint(0, 5, .95, 'agresti-coull')
         method x n mean       lower     upper
1 agresti-coull 0 5    0 -0.05457239 0.4890549

Manually calculating the lower bound (using the formula from Wikipedia):
> n <- 5
> x <- 0
> n2 <- n+qnorm(0.025)^2
> p <- 1/n2 * (0 + .5 * zn^2) 
> zn <- qnorm(0.025)
1/n2 * (0 + .5 * zn^2) + zn * sqrt(1/n2 * p * (1 - p))
[1] -0.05457239

 A: The lower limit of the formula from your link cannot be negative. But the interval from your link is not the Agresti-Coull interval, it is the Wilson interval. The formulas from your link are for the so called Wilson interval and not the Agresti-Coull interval. Agresti and Coull list the formulas from your link in their paper and call it the score confidence interval (page 120). In the superb paper from Brown et al. (2001) Interval estimation for a binomial proportion, it is called the Wilson interval. It is more commonly know as Wilson interval. They show in their article that the Wilson interval performs well even with small n. In the output of binom.confint the Wilson interval is denoted as wilson and can be calculated by setting the method methods="wilson" in binom.confint. Here is the R code for the (modified) Wilson confidence interval:
n <- 5
x <- 0

alpha <- 0.05

p.hat <- x/n

upper.lim <- (p.hat + (qnorm(1-(alpha/2))^2/(2*n)) + qnorm(1-(alpha/2)) * sqrt(((p.hat*(1-p.hat))/n) + (qnorm(1-(alpha/2))^2/(4*n^2))))/(1 + (qnorm(1-(alpha/2))^2/(n)))

lower.lim <- (p.hat + (qnorm(alpha/2)^2/(2*n)) + qnorm(alpha/2) * sqrt(((p.hat*(1-p.hat))/n) + (qnorm(alpha/2)^2/(4*n^2))))/(1 + (qnorm(alpha/2)^2/(n)))

#==============================================================================
# Modification for probabilities close to boundaries
#==============================================================================

if ((n <= 50 & x %in% c(1, 2)) | (n >= 51 & n <= 100 & x %in% c(1:3))) {
  lower.lim <- 0.5 * qchisq(alpha, 2 * x)/n
}

if ((n <= 50 & x %in% c(n - 1, n - 2)) | (n >= 51 & n <= 100 & x %in% c(n - (1:3)))) {
  upper.lim <- 1 - 0.5 * qchisq(alpha, 2 * (n - x))/n
}

upper.lim
[1] 0.4344825

lower.lim
[1] 3.139253e-17

Here, the lower limit is clearly 0 (the remaining is a numerical error). The Wilson interval in the output of binom.confint can be calculated by setting the option methods="wilson" and this is the same as the one we've calculated above:
library(binom)

binom.confint(x=0, n=5, methods="wilson")

  method x n mean        lower     upper
1 wilson 0 5    0 3.139253e-17 0.4344825

The function binom.confint implements the formulas given on the Wikipedia page for the Agresti-Coull interval:
n.hat <- n + qnorm(1-(alpha/2))^2

p.hat <- (1/n.hat) * (x + (1/2)*qnorm(1-(alpha/2))^2)

upper.lim2 <- p.hat + qnorm(1-(alpha/2))*sqrt((1/n.hat)*p.hat*(1-p.hat))

lower.lim2 <- p.hat - qnorm(1-(alpha/2))*sqrt((1/n.hat)*p.hat*(1-p.hat))

upper.lim2
[1] 0.4890549

lower.lim2
[1] -0.05457239

They are the same as the ones by binom.confint with the option methods="ac":
library(binom)

binom.confint(x=0, n=5, methods="ac")

         method x n mean       lower     upper
1 agresti-coull 0 5    0 -0.05457239 0.4890549

