star formation - Simulation packages or theory to work with gravitational collapse of massive molecular clouds?

This sounds very ambitious! I think your way in to this topic might be to have a look at the works of two authors - Matthew Bate and Mark Krumholz.

The reason I suggest these two names is that they are representative of the two basic approaches to simulating the collapse of molecular clouds (which has little to do with n-body simulations). These approaches are known as Smoothed-Particle Hydrodynamics (SPH) and Adaptive Mesh Refinement (AMR).

Crudely put, SPH treats the fluid as discrete particles whose properties are smoothed over the size of a kernel - the resolution is controlled by the number of particles in the simulation and the smoothing length. AMR is a more "traditional" grid-based approach where the size of the grid is varied adaptively with time and space to get the requisite resolution.

A quick google reveals a number of public codes -e.g. for SPH there is SPHYSICS or GRLAB; for AMR there is ENZO or PARAMESH.

One of the main activities in model simulations these days is trying to incorporate "feedback" into star formation. i.e. the winds and radiation from newly born high-mass stars and the outflows from protostars. Such things do not form part of standard codes.