  Dense stellar systems such as globular clusters and dense nuclear clusters
are the breeding ground of sources of gravitational waves for the advanced
detectors LIGO and Virgo. These systems deserve a close study to estimate rates
and parameter distribution. This is not an easy task, since the evolution of a
dense stellar cluster involves the integration of $N$ bodies with high
resolution in time and space and including hard binaries and their encounters
and, in the case of gravitational waves (GWs), one needs to take into account
important relativistic corrections. In this work we present the first
implementation of the effect of spin in mergers in a direct-summation code,
NBODY6. We employ non-spinning post-Newtonian corrections to the Newtonian
accelerations up to 3.5 post-Newtonian (PN) order as well as the spin-orbit
coupling up to next-to-lowest order and the lowest order spin-spin coupling. We
integrate spin precession and add a consistent treatment of mergers. We analyse
the implementation by running a set of two-body experiments and then we run a
set of 500 simulations of a relativistic stellar cluster. In spite of the large
number of mergers in our tests, the application of the algorithm is robust. We
find in particular the formation of a runaway star whose spin decays with the
mass it wins, independently of the initial value of the spins of the stars.
More remarkably, the subset of compact objects that do not undergo many
mergers, and hence represent a more realistic system, has a correlation between
the final absolute spin and the initial choice for the initial distribution,
which could provide us with information about the evolution of spins in dense
clusters once the first detections have started.
