In fact, EVD does allow users to simulate nested logit deviates. The function is rmvevd (random multivariate EVD). The documentation provides an example of how to simulate deviates from the multivariate EVD, but it's a little cumbersome, so I'll provide an example for a simple nested logit problem.
D = 4 alternatives, alternative 1 is in nest 1; alternatives 2,3,4 are in nest 2.
lambda_2 = 0.5 is the nesting parameter.
B = the set of all possible sets that can be formed from the 4 alternatives: (B = {1}, {2}, {3}, {4}, {1,2} {1,3}, {1,4}, {2,3}, {2,4} {3,4}, {1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}, {1,2,3,4})
mvevd requires the user to specify the arguments dep and asy. dep is a vector of nesting parameters on the sets in B of length 2 or greater. asy (in the nested logit example) is a list that specifies the nesting structure. In the example here asy will be vector of 1s in the 1st and 14th entries (corresponding to alternative 1 in its own group, and alternatives 2,3,4 in group 2). dep will have a 1 in all entries except for the 10th entry, which will equal lambda_2 (as {2,3,4} is the 10th element of B that is not of unit length). A code example to simulate 10 4-d nested logit deviates according to this example is as follows:
library('evd')
vdep <- c(rep(1,9), lambda, 1)
asy_list <- list(1, 0, 0, 0, c(0,0), c(0,0), c(0,0), c(0,0), c(0,0), c(0,0),
c(0,0,0), c(0,0,0), c(0,0,0), c(1,1,1), c(0,0,0,0))
deviates <- rmvevd(10, dep = vdep, asy = asy_list, model = "alog", d = 4)