Thanks to Kit Baum a new package is available on SSC. simulate2 and psimulate2 enhance and parallalise Stata's build-in simulate command.

simulate2 follows the syntax of simulate and eases Monte-Carlo Type simulations. The novelties of simulate2 are that it allows the programs called by simulate2 to return macros (strings) to r() and e() which are then saved into the dataset in use. In addition it has advanced options to assign random number states (seed, seed stream and type) to a specific repetition of the simulation. Seeds can be saved to a Stata dataset (dta) or a frame to allow for a seamless replication of simulation results. Simulation results can be directly saved into frames.

psimulate2 parallels simulate2. The number of replications are split into equal blocks and each block is run on a separate Stata instance. This speeds up simulations and allows for a more time efficient use of computational power if the computer in use has two or more cores. psimulate2 creates for each block a do file and a batch file. The batch file is then called in the separate Stata instance.

Currently, psimulate2 only works in Windows machines. Both psimulate2 and simulate2 require Stata 16.

Examples
Make a dataset containing the OLS coefficient, standard error, the current time and save the seeds in a frame called seed_frame. Perform the experiment 1000 times:

Code:
program define testsimul, rclass
                version 16
                syntax anything
                clear
                set obs `anything'
                gen x = rnormal(1,4)
                gen y = 2 + 3*x + rnormal()
                reg y x
                matrix b = e(b)
                matrix se = e(V)
                ereturn clear
                return scalar b = b[1,1]
                return scalar V = se[1,1]
                return local time "`r(current_time)'"
end
Run the simulation 1000 times and save the seed in a frame:
Code:
simulate2 time = r(time) b = r(b) V = r(V), reps(1000) saveseed(seed_frame,frame): testsimul 100
Run the simulation 1000 times but with two parallel instance and use the seed:
Code:
 psimulate2 time = r(time) b = r(b) V = r(V),  reps(500) parallel(2) seed(seed_frame seed, frame): testsimul 100
For more info, please see the help file.

simulate2 was inspired by comments from Alan Riley (StataCorp) and Tim Morris to a discussion during the Stata User Group meeting 2019 in London.