  Introducing the DRAGON simulation project, we present direct $N$-body
simulations of four massive globular clusters (GCs) with $10^6$ stars and 5$\%$
primordial binaries at a high level of accuracy and realism. The GC evolution
is computed with NBODY6++GPU and follows the dynamical and stellar evolution of
individual stars and binaries, kicks of neutron stars and black holes, and the
effect of a tidal field. We investigate the evolution of the luminous (stellar)
and dark (faint stars and stellar remnants) GC components and create mock
observations of the simulations (i.e. photometry, color-magnitude diagrams,
surface brightness and velocity dispersion profiles). By connecting internal
processes to observable features we highlight the formation of a long-lived
'dark' nuclear subsystem made of black holes (BHs), which results in a
two-component structure. The inner core is dominated by the BH subsystem and
experiences a core collapse phase within the first Gyr. It can be detected in
the stellar (luminous) line-of-sight velocity dispersion profiles. The outer
extended core - commonly observed in the (luminous) surface brightness profiles
- shows no collapse features and is continuously expanding. We demonstrate how
a King (1966) model fit to observed clusters might help identify the presence
of post core-collapse BH subsystems. For global observables like core and
half-mass radii the direct simulations agree well with Monte-Carlo models.
Variations in the initial mass function can result in significantly different
GC properties (e.g. density distributions) driven by varying amounts of early
mass loss and the number of forming BHs.
