Context. The star formation process produces hierarchical clustered stellar distributions through gravoturbulent fragmentation of molecular clouds. Simulating stellar dynamics in such an environment is numerically challenging due to the strong coupling between young stars and their surrounding and the large range of length and time scales. Aims. This paper is the first of a suite aimed at investigating the complex early stellar dynamics in star-forming regions, from the initial collapse of the molecular cloud to the phases of complete gas removal. We present a new simulation framework. This advanced framework is the key to generating a larger set of simulations, enabling statistical analysis, which is mandatory to address the stochastic nature of dynamical interactions. Methods. Methods originating from the stellar dynamics community, including regularisation and slow-down methods (SDAR), have been added to the hydrodynamical code Phantom to produce simulations of embedded cluster early dynamics. This is completed by a novel prescription of star formation to initialise stars with a low numerical cost, but in a way that is consistent with the gas distribution during the cloud collapse. Finally, a prescription for H II region expansion has been added to model the gas removal. Results. We have run test-case simulations following the dynamical evolution of stellar clusters from the cloud collapse to a few million years. Our new numerical methods fulfil their function by speeding up the calculation. The N-body dynamics with our novel implementation never appear as a bottleneck that stalls the simulation before its completion. Our new star formation prescription avoids the need to sample individual star formations within the simulated molecular clouds with high resolution. Overall, these new developments allow accurate hybrid simulations in minimal calculation time. Our first simulations show that massive stars largely impact the star formation process and shape the dynamics of the resulting cluster. Depending on the position of these massive stars and the strength of their H II regions, they can prematurely dismantle part of the cloud or trigger a second event of cloud collapse, preferentially forming low-mass stars. This leads to different stellar distributions for numerical simulations with similar initial conditions and confirms the need for statistical studies. Quantitatively, and despite the implementation of feedback effects, the final star formation efficiencies are too high compared with those measured in molecular clouds of the Milky Way. This is probably due to the lack of feedback mechanisms other than H II regions, in particular jets, non-ionising radiation, or Galactic shear. Conclusions. Our new Phantom N-Body framework, coupled with the novel prescription of star formation, enables the efficient simulation of the formation and evolution of star clusters. It enables the statistical analysis needed to establish a solid theoretical framework for the dynamical evolution of embedded star clusters, continuing the work done in the stellar dynamics community.