Many features of GADGET-4 are controlled by a parameterfile that has to be specified whenever the code is started. Each parameter value is set by specifying a keyword, followed by a numerical value or a character string, separated by whitespace (spaces, tabs). For each keyword, a separate line needs to be used, with the value and keyword appearing on the same line. Between keywords, arbitrary amounts of whitespace (including empty lines) may be used. The order in which the keywords are specified is arbitrary, but each keyword needs to appear exactly once, even if its value is not relevant for the particular simulation (otherwise you'll get an error message). Note that the keywords are type-sensitive.

Lines with a leading '%' or '#' are ignored. In lines with keywords, comments may also be added after the specification of the value for the corresponding keyword, and in this case do not need to begin with any special character.

In the following, each keyword and the meaning of its value are discussed in detail, and typical example values are provided as an illustration. Some keywords may be changed during a run, i.e. changes of the corresponding values will be taken into account upon resuming the code from restartfiles whenever this is reasonable, but changes of certain other parameter values will be ignored. For example, while changing the memory- or cpu-time limit specified for the code is always possible, changing the cosmological parameters in the middle of a run will be prevented. If a change in the middle of a run is accepted by the code, this will also be reflected in the log messages when starting the code. Note, however, that you normally do not need to make any changes in the parameterfile when you restart a run from restart-files.

Filenames and file formats

OutputDir /home/volker/galaxy_collision

This is the pathname of the directory that holds all the output generated by the simulation (snapshot files, restart files, diagnostic files). The code will try to create this directory if it does not yet exist, but the directory's parent directory needs to exist otherwise an error message will be produced.

SnapshotFileBase snapshot

From this string, the name of the snapshot file is derived by adding an underscore, and the number of the snapshot in a 3-digits format. If NumFilesPerSnapshot > 1 , each snapshot is distributed into several files, with a group of processors writing their data to one of the files (these files can be written concurrently). In this case, the filenames are supplemented with a tailing .n , where n designates the rank of the file in the group of files representing the snapshot. If the HDF5 file format is used, an identifier .hdf5 is also appended automatically.

SnapFormat 3

A flag that specifies the file-format to be used for writing snapshot files. A value of 1 selects the simple legacy binary file-format of GADGET-1/2/3, while a value of 2 selects a more convenient variant of this simple binary format (which has been available from GADGET-2 onwards). A value of 3 selects the use of HDF5 instead, which is the strongly recommended format. This is because this data format allows a simple browsing of its contents, access to individual data items is easily possible through the name of the data set, conversions between endianness and single/double precision are automatically done if needed, and pesky I/O errors (like reading too many items from a given data set) are detected reliably. HDF5 output has been introduced for snapshots already in GADGET-2/3, but in GADGET-4, it is now also available for group catalogues, merger trees, light-cone data, etc., and this parameter also regulates the file format of these outputs. If HDF5 is selected, the filenames will be appended automatically with a .hdf5 suffix. The structure of the snapshot files and of other outputs is discussed in a separate part of the documentation.

ICFormat 1

This flag selects the file format of the initial conditions read in by the code upon start-up. The possible format choices are the same as for SnapFormat. It is therefore possible to use different formats for the initial conditions and the produced snapshot files. In case your initial conditions is in a different file format, we recommend that you convert your IC files to one of the three formats supported by GADGET-4 (preferably HDF5). Alternatively, you could incorporate a customised reading routine directly into the GADGET-4 code, but this requires an intimate understanding of the internal workings of the code.

InitCondFile /home/volker/ICs/galaxy.dat

This sets the filename of the initial conditions to be read in at start-up. Note that the initial conditions file (or files) does not have to reside in the output directory. The initial conditions can be distributed into several files, in the same way as snapshot files. In this case, only the basename without the tailing .n number should be specified as initial conditions filename. Likewise the .hdf5 file-name suffix should be omitted if HDF5 is used. The code will recognise the number of files that make up the initial conditions from the file header entries, and load all of these files accordingly. There is no constraint on the number of these files in relation to the processor number used.

NumFilesPerSnapshot 2

The code can distribute each snapshot onto several files. This leads to files that are easier to handle in case the simulation is very large, and also speeds up I/O, because these files can then be written or read in parallel. The number of processors should be equal or larger than NumFilesPerSnapshot , because each snapshot file will hold the data of a group of processors (otherwise the code reduces the value to the number of MPI ranks used). Optimum I/O throughput is reached if the number of processors is equal to, or a multiple of NumFilesPerSnapshot , and if MaxFilesWithConcurrentIO is reasonably large. With the setting NumFilesPerSnapshot=1 it is possible to write all particle data in just one snapshot file but then no parallel I/O is used.

MaxFilesWithConcurrentIO 8

This sets the number of concurrent I/O operations the code is allowed to carry out. If 0 is specified, the code adopts a value equal to the number of MPI ranks, and the same is done if the specified value is larger than the number of MPI rank. However, it can sometimes be sensible to limit this to a smaller number, especially on very large MPI partitions in order to prevent that the I/O subsystem is overloaded. Whether or not this is an issue can be found out by starting the code with restartflag 9, which carries out a special I/O bandwidth test where MaxFilesWithConcurrentIO is systematically varied from the number of MPI-ranks down to 1 by factors of 2, and in each case a write-test is performed using parallel I/O from all ranks. This can be used to determine a reasonable setting for MaxFilesWithConcurrentIO that is not causing choking of the filesystem.

CPU-time limit and restarts

TimeLimitCPU 40000.0

This is the wallclock time limit for the current execution of the code, in seconds. Often the code will be run through a submission to a computing queue, and hence this value should be matched to the corresponding time limit of the computing queue or job submission script, as appropriate. The run will automatically interrupt itself and write restart files if 85% of this time has elapsed. The extra 15% is introduced to make sure that there is always enough time left to safely finish the current time step (or FOF/SUBFIND group finding) and for writing the restart files before the time limit is reached. Note that this time refers to the wall-clock time on one processor only. The total CPU time consumed by the code is obtained by multiplying with the total number of cores that are used/occupied by the run.

CpuTimeBetRestartFile 7200

This is the maximum amount of wall-clock time, in seconds, that may elapse before the code writes a new set of restart file for regular checkpointing. With this parameter the code can hence be asked to write a restart file every once in a while. This is meant to provide some protection against hardware or software failures, in which case one can resume a simulation from the last set of restart files. In the above example, a restart file would be written automatically every 2 hours, so that the lost time in case of such an issue would be at most 2 hours of computing. The old set of restart files is renamed into a set of bak-restart.X files before the new files are written, so that there is some protection against a crash during the writing of the restart files themselves. The latter could happen, for example, because of a disk-full error. It is not possible to resume a simulation from a corrupted set of restart files, or with restart files that correspond to a mix of different output times. In case the code should really crash while writing restart files, it is best to discard all restart.X files, and then to rename the bak-restart.X files into restart.X files (where X runs from 0 to the number of MPI ranks minus 1), for example with the command:

rename bak-restart restart bak-restart.*

Note that depending on your system, the convenient rename command may have a slightly different syntax or may not be available at all (consult the man pages).

FlushCpuTimeDiff 120

The GADGET-4 code provides quite verbose output in its log-files, with each timestep producing some entries. If you carry out a simulation with a very large number of very short timesteps and you have a slow or busy filesystem, this I/O can slow down the code if the filesystem buffers are flushed to disk after every output. With this parameter, the flush operations are only carried out with a reduced frequency, in the above example every 120 seconds. This should avoid any significant cost of this I/O, but it also means that you may not see right away what the code has been doing when inspecting the log-files, because the information in these files will typically only be updated when a new flush operation is triggered by the code.

Memory allocation

MaxMemSize 2000

This value gives the maximum amount of memory (in MByte) the code is allowed to use per MPI process. The code will strictly enforce this limit, and terminate if a higher use is attempted by the code. If this should occur, a table with the memory allocated for different code parts is output to the log file, together with information in which line of the code each allocation has happened. Note that the code will automatically try to make good use (in communication phases) of any extra memory you specify here. It is therefore a good idea to set MaxMemSize to something close to the amount of physical memory that can be used per MPI rank on the target compute nodes. The code will check whether the memory on the target machines is actually sufficient for the setting chosen, and otherwise terminate (this check only works on Linux). This check should safely prevent that you accidentally run the code with a too high memory use that could make a compute-node start swapping. (Driving a node into swapping is generally a really bad thing and would lead to dismal performance and/or node crashes. Most HPC systems prevent this nowadays at the system level for their compute nodes.) If you want to find out the smallest amount of memory a given simulation would have needed to complete, you can grep the log-file memory.txt for the values reported behind the phrase "Largest Allocation Without Generic". This gives the maximum that was needed by any of the MPI ranks over the course of the simulation and corresponds to the lower possible limit for MaxMemSize for the given run minus a small amount of communication buffer space that you have to allow the code to use. So in practice, the setting for MaxMemSize needs to be at least slightly larger than the "Largest Allocation Without Generic" value reported in memory.txt.

Simulated time span and spatial extent

TimeBegin 0

This initialises the time variable of the simulation when a run is started from initial conditions (in internal units). If comoving integration is selected (ComovingIntegrationOn=1), the time variable is the dimensionless expansion factor a itself, i.e. TimeBegin = a = 1 / (1 + z_start) , otherwise it is simply physical time in the internal system of units of the code.

TimeMax 3.0

This marks the end of the simulation. The simulation will run up to this point, then write a restart-file and a snapshot file corresponding to this time (even if the time TimeMax is not in the normal sequence of snapshot files). If TimeMax is increased later on, the simulation can be simply continued from the last restart-file. Note that this last snapshot file will then be overwritten in case this was a special dump out of the normally expected output sequence. For comoving integrations, the time variable is the expansion factor, e.g. TimeMax=1.0 will stop the simulation at redshift z=0 . Otherwise the value of TimeMax refers to physical time.

ComovingIntegrationOn 0

This flag enables or disables comoving integration in an expanding universe. For ComovingIntegrationOn=0 , the code uses plain Newtonian physics with vacuum or periodic boundary conditions. Time, positions, velocities, and masses are measured in the internal system of units, as specified by the selected system of units. For ComovingIntegrationOn=1 , the integration is carried out in an expanding universe, using a cosmological model as specified by Omega0 , OmegaLambda , etc. In this cosmological mode, coordinates are comoving, and the time variable is the natural logarithm of the expansion factor itself. If the code has not been compiled with the PERIODIC makefile option, the underlying model makes use of vacuum boundary conditions, i.e. density fluctuations outside the particle distribution are assumed to be zero. This requires that your particle distribution represents a spherical region of space around the origin. If PERIODIC is enabled, the code expects the particle coordinates to lie in the interval [0, BoxSize].

BoxSize 10000.0

The size of the periodic box (in code units) encompassing the simulation volume. This parameter is only relevant if the PERIODIC option is activated.

Cosmological parameters

Omega0 0.3

Cosmological matter density parameter in units of the critical density at z=0. Relevant only for comoving integration.

OmegaLambda 0.7

Cosmological vacuum energy density (cosmological constant) in units of the critical density at z=0. For a geometrically flat universe, one has Omega0 + OmegaLambda = 1. Important: For simulations in Newtonian space that do not account for cosmological expansion, this parameter has to be set to zero.

OmegaBaryon 0.04

Baryon density in units of the critical density at z=0. This is not explicitly used in the time integration of GADGET-4, but the parameter is relevant when initial conditions are created, or when dark matter-only initial conditions are outfitted with gas particles upon code start-up with the GENERATE_GAS_IN_ICS option.

HubbleParam 0.7

This dimensionless parameter enters the definition of GADGET's system of units, and can be used to eliminate an explicit dependence on the value of the Hubble constant, like it has been traditionally done in cosmology. Most often, the value of HubbleParam is chosen to express the value of the Hubble constant in units of 100 km/s/Mpc. While the value of HubbleParam is not needed (in fact, it does not enter the computations at all) in purely collisionless simulations, the value of HubbleParam is still relevant when conversions to physical cgs units are required, for example to compute rate equations in radiative cooling physics.

Hubble 100.0

Value of the Hubble constant in internal units. Since the internal units contain a factor HubbleParam, one can basically choose whether one wants to set the Hubble constant via HubbleParam (then Hubble has the same value in all simulations, even if the cosmological factors and the Hubble constant change), or one sets HubbleParam to unity and uses Hubble to directly set the Hubble constant. Both is possible, and intermediate forms in principle as well.

System of units

UnitLength_in_cm 3.085678e21

This sets the internal length unit in cm/h, where H_0 = 100 h km/s/Mpc. The above choice is convenient for cosmology, as it sets the length unit to 1.0 kpc/h.

UnitMass_in_g 1.989e43

This sets the internal mass unit in g/h, where H_0 = 100 h km/s/Mpc. The above choice is convenient for cosmology, as it sets the mass unit to 10^10 M_sun/h.

UnitVelocity_in_cm_per_s 1.0e5

This sets the internal velocity unit in cm/sec. The above choice corresponds to a velocity unit of km/sec, which is the commonly used and most convenient unit in cosmology. Note that the specification of UnitLength_in_cm , UnitMass_in_g and UnitVelocity_in_cm_per_s also determines the internal unit of time. The definitions made above imply that in internal units the Hubble constant has a numerical value independent of h (where h is given by HubbleParam ). For the numerical examples above, the Hubble constant has always the value 0.1 in internal units, independent of h, and the Hubble time is always 10.0 in internal units, with one internal time unit corresponding to 9.8 x 10^8 yr/h. However, of course, you are free to choose a different system of units if you like. Note again that this implies that for purely gravitational dynamics, the code will not need to know the value of h at all. HubbleParam is nevertheless kept in the parameterfile because additional physics in the hydrodynamical sector may require it.

GravityConstantInternal 0

The numerical value of the gravitational constant G in internal units depends on the system of units you choose. For example, for the numerical choices made above, the physical value of G corresponds to G=43007.1 in internal units. For GravityConstantInternal=0 (the normal choice for cosmological simulations), the code calculates the internal value corresponding to the physical value of G automatically. But sometimes, you might want to set G yourself. For example, in scale-free test simulations, specifying GravityConstantInternal=1, UnitLength_in_cm=1, UnitMass_in_g=1, and UnitVelocity_in_cm_per_s=1, yields a natural system of units in which one may also want to adopt G=1 as well, which can then be achieved by specifying a non-zero value for GravityConstantInternal, in this example GravityConstantInternal=1.

Gravitational force accuracy

TypeOfOpeningCriterion 0

This selects the type of cell-opening criterion used in the tree walks for computing gravitational forces. A value of 0 results in a geometric opening criterion which is primarily governed by the opening angle theta, while 1 selects a relative criterion that tries to limit the absolute truncation error of the multipole expansion for every particle-cell or cell-cell interaction. The latter scheme usually gives slightly higher accuracy at a comparable level of computational cost compared with the geometric criterion. When it becomes more demanding to calculate accurate forces due to strong cancellation effects (such as at high redshift with nearly uniform mass distribution), the relative criterion automatically invests more effort because of the small residual forces there. This adaptivity makes it advantageous especially for cosmological simulations, where a single value of theta for the geometric criterion is not ideal at all redshifts.

ErrTolTheta 0.7

This is the accuracy criterion parameter (the opening angle theta) of the tree algorithm if the geometric opening criterion (i.e. TypeOfOpeningCriterion=0 is used). If TypeOfOpeningCriterion=1 is adopted, then theta and the geometric opening criterion are only used for a first force computation whose purpose is only to provide an estimate for the current acceleration of each particle, which in turn is needed to compute forces with the relative cell opening criterion. Hence, theta needs to be set to a sensible value even if the relative criterion is used and only forces computed with the relative criterion enter the dynamics.

ErrTolForceAcc 0.005

This controls the accuracy of the relative cell-opening criterion (if enabled). Here, alpha is given by ErrTolForceAcc. Note that independent of this relative criterion, the code will always open nodes if the point of reference lies within a geometric boundary box around the cubical cell (unless TREE_NO_SAFETY_BOX is enabled). This protects against the possibility of an occurrence of unusually large force errors for very particular particle configurations.

ErrTolThetaMax 1.0

When the relative opening criterion is used, the effective opening angle allowed for a node of little mass may grow very large, possibly approaching the convergence radius of the multipole expansion. To protect against the possibility to get unexpectedly large errors from this, the maximum allowed geometric opening angle can be limited with this parameter.

ActivePartFracForPMinsteadOfEwald 0.1

This parameter is only needed when the TreePM scheme is used in combination with the TREEPM_NOTIMESPLIT option. Then the time integration does not distinguish between a long range and a short range force, instead the total force is integrated with a single (variable) timestep. The TreePM method here only functions as a method for accelerating the computation of the total force, with the alternative being to do it with a pure tree calculation (if needed with Ewald correction for periodic boundaries). If only a small number of particles is active, doing the force calculation as a pure tree can be faster than doing it with the TreePM approach, because for the PM part always full FFTs have to be computed independent of the number of active particles. This parameter is used to specify a threshold above which the active fraction needs to lie before TreePM is applied for a given timestep, otherwise the force calculation is done with a pure tree. The same applies to FMM-PM and pure FMM. In principle, this should only affect the performance of the calculation, not its accuracy in any significant way (this is true in the limit when the TreePM and pure Tree force errors are comparable in magnitude, and are both negligible).

Time integration accuracy

MaxSizeTimestep 0.01

This parameter sets the maximum timestep a particle may take. This should be set to a sensible value in order to protect against too large timesteps for particles with very small acceleration. Usually, a few percent of the dynamical time of the system gives sufficient accuracy. For cosmological simulations, the parameter specifies the maximum allowed step in ln(a), because the natural logarithm of the scale factor is discretized for the time integration in this case. Hence, specifying the maximum allowed timestep for cosmological simulations is equivalent to specifying it as a fraction of the current Hubble time. A value of ~0.01 is usually accurate enough for most cosmological runs.

MinSizeTimestep 0

If a particle requests a timestep smaller than the specified value of this parameter, the simulation terminates with an error message. This is meant to prevent simulations from continuing when the timestep has dropped to an unreasonably small value, because such behaviour typically indicates a problem of some sort. Setting the parameter to zero disables this safety check.

ErrTolIntAccuracy 0.025

This dimensionless parameter controls the accuracy of the simple kinematical timestep criterion commonly employed in cosmological simulations, and which is also used in GADGET-4. The timestep constraint is given by dt = sqrt(2 eta epsilon/|a|), where eta=ErrTolIntAccuracy, epsilon is the gravitational softening length, and a the acceleration experienced by the particle. The actual timestep taken by the particle will always be shorter than dt, as the particle will be forced onto the power-of-two hierarchy of allowed timestep sizes by reducing the step to the next available shorter step.

CourantFac 0.15

This sets the value of the Courant coefficient used in the determination of the hydrodynamical timestep of SPH particles. Note that GADGET-4's definition of the SPH smoothing length differs by a factor of 2 from that found in some part of the SPH literature. As a consequence, comparable settings of CourantFac may be a factor of 2 smaller in GADGET4 when compared with codes using a different convention.

Domain decomposition

ActivePartFracForNewDomainDecomp 0.01

A new domain decomposition is not necessarily determined for every single timestep. A value of ActivePartFracForNewDomainDecomp=0.01, for example, means that the domain decomposition is reconstructed whenever there are at least 0.01 N particles active at the current synchronization time, where N is the total particle number. Note that the gravitational tree is always reconstructed in every step, whereas the neighbor search tree is only reconstructed in case a domain decomposition is done for the current step. Otherwise it expands its nodes as needed to accommodate all the SPH particles that were grouped into each node.

TopNodeFactor 2.5

The domain decomposition involves the construction of a coarse oct-tree whose leaf nodes tessellate the simulation volume. The TopNodeFactor regulates how fine this top-level tree gets (it is designated as f_top in the code paper). The code will roughly produce TopNodeFactor * NTask * f_mult leaf-nodes in the top-level tree, focusing always on refining the most loaded one first until the desired fineness is reached. Here f_mult refers to the number of different cost categories that are balanced simultaneously. Since a tree node can only be assigned in full to individual domains, this parameter influences the level of discreteness fluctuations present in the load among the set of multiple domains that are mapped to individual MPI ranks. A value of a few for this parameter is usually good enough. If a very large value is adopted, the top-level tree (which is identically stored on all nodes) may get very large, making its memory use and construction time costly.

Output frequency

OutputListOn 0

A value of 1 signals that the output times are given in the file specified by OutputListFilename Otherwise, the output times are generated automatically in the way described below. We note that the code will only generate snapshot files if full timesteps have been finished (this is different from GADGET-2) and thus the full system is synchronized in time. This means that in GADGET-4, each particle has finished an integer number of full KDK timesteps when stored in snapshots. Desired output times are given either in the file with output times or are created in a regularly spaced way (as described below). The corresponding desired output times will always be mapped to the closest available output time. The set of these available output times is basically given by the simulation timespan divided by maximum used time step size. The mapping then means that the actual output time of a snapshot can deviate from the desired output at most by 0.5 times the maximum timestep actually used in the simulation when the output occurs. If you want to have many outputs with very fine spacing, it makes sense to set MaxSizeTiStep sufficiently small, in particular smaller than the desired output spacing otherwise the number of created snapshots could be lower than desired in case the simulation takes timesteps for some particles that are larger than the desired output spacing.

OutputListFilename output_times.txt

This specifies the name of a file that contains a list of desired output times. If OutputListOn is set to 1, this list will determine the times when snapshot-files are desired. The file given by OutputListFilename should just contain the floating point values of the desired output times in plain ASCII format. The times do not have to be ordered in time, but there may be at most 1100 values (this is the default, but it can be enlarged if desired by setting the MAXLEN_OUTPUTLIST contant in Output times that are in the past relative to the current simulation time will always be ignored.

TimeOfFirstSnapshot 0.047619048

This variable selects the time for the first snapshot (relevant only if OutputListOn=0). For comoving integration, the above choice would therefore produce the first dump at redshift z=20.

TimeBetSnapshot 1.0627825

If OutputListOn=1 this parameter is ignored. Otherwise, after a snapshot has been written, the time for the next snapshot is determined by either adding TimeBetSnapshot to TimeOfFirstSnapshot, or by multiplying TimeOfFirstSnapshot with TimeBetSnapshot. The latter is done for comoving integration, and will hence lead to a series of outputs that are equally spaced in ln(a). The above example steps down to redshift z=0 in 50 logarithmically spaced steps.

TimeBetStatistics 0.1

This determines the interval of time between two subsequent computations of the total potential energy of the system. This information is then written to the file energy.txt, together with information about the kinetic energies of the different particle types. A first energy statistics is always produced at the start of the simulation at time TimeBegin.

SPH parameters

DesNumNgb 64

This is the desired number of SPH smoothing neighbours. Normally, the effective number of neighbours (defined as the mass inside the kernel divided by the particle mass) is kept constant very close to this value. Should it ever try to get outside a range +/- MaxNumNgbDeviation from DesNumNgb, the code will readjust the smoothing length such that the number of neighbours is again in this range.

MaxNumNgbDeviation 2

This sets the allowed variation of the number of neighbours around the target value DesNumNgb.

InitGasTemp 10000

This sets the initial gas temperature in Kelvin when initial conditions are read. However, the gas temperature is only set to a certain temperature if InitGasTemp > 0 and if at the same time the temperature of the gas particles in the initial conditions file was found to be zero, otherwise the initial gas temperature is left at the value stored in the IC file. If the temperature is set through this parameter, and if it is below 10^4 K, a mean molecular weight corresponding to neutral gas of primordial abundance is assumed, otherwise complete ionisation is assumed.

ArtBulkViscConst 1.0

This sets the value of the artificial viscosity parameter alpha_visc used by GADGET-4. See code paper for details.


This sets the minimum value of the artificial viscosity parameter when a time-dependent viscosity is enabled.

Gravitational softening

The code distinguishes between different particle types. As far as gravity is concerned, all the types are treated equally by the code. The particles of the first type (type=0) are treated as SPH particles and receive an additional hydrodynamic acceleration from pressure gradients. The concept of particle types is primarily introduced to simplify analysis and to give certain particles an easily identifiable role.

The default number of types is NTYPES=6, which was also used as a fixed setting in GADGET-1/2/3. There, the six particle types were referred to with the symbolic tags "Gas", "Halo", "Disk", "Bulge", "Stars", and "Bndry", in this order, but these names are now dropped in favour of just numerical type specifiers. The number of available types can be enlarged or reduced if needed, but a value equal to 6 needs to be used if backwards compatibility to the format of older versions of GADGET is desired.

Normally, each particle type is mapped to a certain gravitational softening length. The number of available different softening lengths is given by NSOFTCLASSES, and does not necessarily have to be equal to NTYPES (this is however the default).

SofteningClassOfPartType0 0

Specifies the softening class that should be assigned to particle type=0. Depending on the setting of NTYPES, additional such parameters are needed, one for each particle type. One hence needs to specify SofteningClassOfPartType0, SofteningClassOfPartType1, ..., SofteningClassOfPartTypeX, where X = NTYPES - 1. The values that are assigned to these parameters need to be in the range [0, NSOFTCLASSES - 1]. It is allowed to map several particle types to the same softening class.

SofteningComovingClass0 0.5

This specifies the (comoving) softening length of the first softening class. Depending on the setting of NSOFTCLASSES, additional such parameters are needed, one for each softening class. One hence needs to specify SofteningComovingType0, SofteningComovingType1, ..., SofteningComovingTypeX, where X = NSOFTCLASSES - 1.

Gravity is softened with a spline kernel in GADGET-4, as outlined in the code paper. The softenings quoted here all refer to epsilon, the equivalent Plummer softening length. Note that for the spline that is used, the force will be exactly Newtonian beyond r = 2.8 epsilon, and the potential of a point mass m at zero lag is phi(0) = -G*m/epsilon. The softening lengths are given in internal length units. For comoving integration, the softening refers to the one employed in comoving coordinates, which usually stays fixed during the simulation.

SofteningMaxPhysClass0 0.5

In cosmological simulations, one sometimes wants to start a simulation with a softening epsilon_com that is fixed in comoving coordinates (where the physical softening, epsilon_phys = a * epsilon_com , then grows proportional to the scale factor a ), but at a certain redshift one may want to freeze the resulting growth of the physical softening epsilon_phys at a certain maximum value. These maximum softening lengths are specified by the SofteningMaxPhysClassX parameters. In the actual implementation, the code uses epsilon_com = min(epsilon_com, epsilon_phys^max / a) as comoving softening. Note that this feature is only enabled for ComovingIntegrationOn=1, otherwise the SofteningMaxPhysClassX values are ignored. The specific parameter SofteningMaxPhysClass0 specifies the maximum physical softening of the first softening class. Depending on the setting of NSOFTCLASSES, additional such parameters are needed, one for each softening class. One hence needs to specify SofteningMaxPhysClass0, SofteningMaxPhysClass1, ..., SofteningMaxPhysClassX, where X = NSOFTCLASSES - 1.

GasSoftFactor 1.5

This parameter is only needed if ADAPTIVE_HYDRO_SOFTENING is activated. In this case, the gravitational softening for the gas particles is individually selected based on their smoothing length from a logarithmic table of available softening classes. The organization of this table is described by the parameters below. In practice, the smoothing length of the gas particle is multiplied by GasSoftFactor and then the closest softening (in terms of smallest difference in the log of the softenings) from the table is assigned as the softening class of the gas particle.

MinimumComovingHydroSoftening 0.001

Specifies the smallest allowed gaseous softening value. Together with the multiplicative factor AdaptiveHydroSofteningSpacing that describes the increase from step to step, this defines the discrete table of available softenings for SPH particles when ADAPTIVE_HYDRO_SOFTENING is used.

AdaptiveHydroSofteningSpacing 1.05

This parameter defines the spacing between two adjacent available SPH softening classes when ADAPTIVE_HYDRO_SOFTENING is enabled. The total number of the set of available softenings classes is given by the constant NSOFTCLASSES_HYDRO, which is normally set to 64 as default. (This can be changed, if desired, by overriding the value of this constant in The largest available gaseous softening length is then given by MinimumComovingHydroSoftening * AdaptiveHydroSofteningSpacing ^ (NSOFTCLASSES_HYDRO - 1).

Subfind parameters

DesLinkNgb 20

This sets the number of neighboring particles that are examined in the SUBFIND algorithm to identify locally isolated peaks in the density field, and to link overdensity candidates across saddle points. The typical value for this quantity is 20, and results of SUBFIND should be quite insensitive to the exact choice. The parameter needs only be specified if SUBFIND is actually enabled in

Initial conditions generation

The following initial conditions parameters are only used if NGENIC is activated in the code configuration file. The value of NGENIC is interpreted as the size of the FFT grid used to compute the displacement field. One should have NGENIC >= Nsample. The redshift of the initial conditions is the same as the defined starting redshift of the simulation, hence is given by TimeBegin.

NSample 128

This sets the maximum wave number k that the code uses, i.e. this effectively determines the Nyquist frequency that the code assumes, k_Nyquist = 2*PI / BoxSize * Nsample/2 Normally, one chooses Nsample such that Ntot = Nsample^3, where Ntot is the total number of particles.

GridSize 128

This parameter is only needed if CREATE_GRID is activated. In this case, the code will create the initial particle load itself in terms of a uniform Cartesian grid with particles of constant mass. The total number of particles that is used is then GridSize^3 . In cold dark matter, perturbations should then be imprinted all the way to the Nyquist frequency of this particle grid, i.e. one should pick GridSize=Nsample. If CREATE_GRID is not active, then this parameter is not present. In this case, the initial unperturbed particle load is read in as the specified IC file.

PowerSpectrumType 2

This can be used to select the parameterization of the linear theory input spectrum. For a value of 1, an analytic fitting function by Eisenstein & Hu is selected, while 2 uses a tabulated power spectrum in the file specified with PowerSpectrumFile A value of 3 (or any other value) will use an analytic parameterization from Efstathiou.

PowerSpectrumFile cmb_code_wmap7_spectrum.txt

This file gives a tabulated linear theory input power spectrum, which can be computed by a Boltzmann code, like CAMB, for example. The file format is ASCII, and should contain two columns, with a pair of values on every line:

log(k) log(Delta^2)

Here log is the base10 logarithm, and k is given in units of h / cm / InputSpectrum_UnitLength_in_cm (see below). Delta^2 refers to the dimensionless power spectrum at wavenumber k, which is related to the ordinary power spectrum P(k) through Delta^2 = 4 PI k^3 P(k).

InputSpectrum_UnitLength_in_cm 3.085678e24

Defines the length unit used in the tabulated input spectrum in cm/h. If desired, this can be chosen differently from UnitLength_in_cm, so that one can for example have an input spectrum table based on Mpc/h while the simulation one carries out uses kpc/h.

ShapeGamma 0.21

This parameter is only relevant when PowerSpectrumType=3 is used (Efstathiou parameterization). In this case, ShapeGamma should usually be set to something close to Omega0 * HubbelParam.

PrimordialIndex 1.0

This may be used to specify a tilt in the primordial index, provided this has not already been taken care of in a tabulated input spectrum. Effectively, this option multiplies the spectrum (whatever the source) with an additional factor k^(PrimordialIndex-1.0), i.e. if one does not want to do this, one has to set this parameter to a value of 1. In particular, if one uses a tabulated input spectrum and this already accounts for a primordial tilt, this parameter nevertheless still has to be set to 1.0.

Sigma8 0.86

Normalization of the linear theory input power spectrum when extrapolated to z=0. As is commonly done, sigma8 gives the rms density contrast fluctuations in top-hat spheres of radius 8 Mpc/h.

ReNormalizeInputSpectrum 1

If set to zero, the tabulated input spectrum is assumed to be correctly normalized already in its amplitude to the adopted starting redshift, otherwise the normalization is recomputed based on the specified sigma8 value, and the linear theory growth factor for the specified cosmology. Normally, this renormalization is always recommended.

SphereMode 0

If this is activated by setting it to 1, only modes with |k| < k_Nyquist are used (i.e. a sphere in k-space), otherwise modes with |k_x|,|k_y|,|k_z| < k_Nyquist are used (i.e. a cube in k-space).

Seed 123456

An integer number that serves as seed for the random number generator used by the IC code. Because the k-modes are filled systematically "inside-out" in k-space, changing the resolution but keeping the seed the same will yield the same large-scale modes. This means that for a fixed seed, one can easily carry out resolution studies on an object-by-object basis, if desired.

Lightcone output

The lightcone output will first collect particles that cross the backwards lightcone in an intermediate storage buffer. The particle positions and velocities are extrapolated to the point of the light-cone crossing. Once the intermediate storage buffer has accumulated a certain particle number, the particles are written to disk in a light-cone file similarly to a snapshot file. After that the intermediate buffer is emptied, and gradually filled again by the continuing lightcone until the next output occurs. The individual lightcone files are thus spherical shells (or segments thereof) around the observer position. When concatenated they yield the full lightcone.

LightConeDefinitionFile lightcones.txt

When continuous light-cone outputting is activated via the switch LIGHTCONE, the geometry of the lightcone(s) is specified via this file. In particular, the start and end redshifts of the lightcone(s) are specified in this file. Note that the lightcone feature only makes sense for cosmological simulations with ComovingIntegrationOn=1. The constant LIGHTCONE_ALLOC_FAC (default value 0.1) determines how much storage space the code sets aside for the intermediate buffer before lightcone files are flushed to disk. This is done by prescribing the particle number that can fit into the intermediate light-cone buffer as a fraction of the total particle number. This indirectly determines on how many output files the lightcone is distributed. By overriding the constant LIGHTCONE_ALLOC_FAC in the, this can be influenced if desired.

The file defining the lightcones has the following format. Each line defines a separate lightcone, and is defined at least by four numbers:

Here can be either 0, 1, 2, 3, or 4, and defines the geometric selection of the specific lightcone, according to: 0 = full sky, 1 = one octant, 2 = a pencil beam cone with circular aperture, 3 = a disk like region, 4 = a pencil beam with a square-shaped aperture.

sets the far edge of the lightcone to redshift zmax = 1/astart - 1 sets the near edge of the lightcone to redshift zend = 1/aend - 1

is normally 0. A value of 1 only makes sense combined with the SUBFIND_ORPHAN_TREATMENT option and then restricts the output to formerly most bound particles.

For type=1 (octant), an additional number with value 0, 1, ..., or 7 is needed to select the specific octant.

For type=2 (cone), three additional numbers are expected to define the principal direction vector of the cone (this does not need to be normalized), followed by its half opening angle in degrees.

For type=3 (disk), three additional number are expected that define the normal of the disk region, followed by a further number defining its comoving thickness.

For type=4 (pyramid with square base), first three additional number are expected to define the direction vector of the pencil beam, then a further vector is expected that is used to set the orientation of the x-direction of the patch of sky that is mapped by the lightcone, with the y-direction being orthogonal to that. Finally, a last number gives the half opening angle of the pyramid-shaped pencil beam in degrees.

Multiple lightcones, also of the same type, can be specified if desired. Note that the output will get extremely large if you select even a moderate redshift depth, because the code will automatically periodically replicate the simulation box as needed to cover the specified lightcone geometry.

An example for a lightcone definition file could look like this:

0 0   0.5    1.0
1 0   0.4    1.0    0

This would define a full-sky light cone from z=1 to z=0, and an octant covering positive x>0,y>0,z>0 from redshift z=1.5 to z=0.

LightConeOriginsFile lightcone_origins.txt

Only when LIGHTCONE_MULTIPLE_ORIGINS is activated, this option is required. One can then supply a list of coordinate triples, each of which is a possible origin of a lightcone as defined above. The invividual lightcone defintions from above then require one additional number for each lightcone at the end. This number is an index into the listed origins, and thus selected the corresponding origin.

LightConeMassMapsNside 12

The healpix Nside parameter defining the angular resolution used for the mass maps projections of lightcone particles. To enable this, the LIGHTCONE_MASSMAPS option needs to be set.

LightConeMassMapThickness 25

Comoving thickness of one lightcone massmap shell.

LightConeMassMapMaxRedshift 5.0

Redshift out to which the massmaps should extend.

Cooling and star formation

The following parameters refer to the simple cooling and star formation model described in Springel & Hernquist (2003,

MaxSfrTimescale 1.5

Gas consumption time scale in internal time units at the threshold density for star formation. This sets the parameter t_0^star in the above paper.

TempClouds 1000.0

This is given in Kelvin and corresponds to the T_c parameter (or correspondingly to u_c when expressed as thermal energy per unit mass) in the above paper.

TempSupernova 1.0e8

This is given in Kelvin and corresponds to the T_SN parameter in the above paper, or correspondingly to u_SN when expressed as thermal energy per unit mass. The relation between these two quantities is T_SN = (2/3) m u_SN / k_B, where m is the mean molecular mass and k_B the Boltzmann constant.

FactorEVP 1000.0

This is the cloud evaporation parameter A_0 in the above paper.

FactorSN 0.1

This parameter (which corresponds to the symbol beta in the above paper) is the mass fraction of short-lived massive stars (>8 Msun) formed for each initial population of stars. This depends on the adopted stellar initial mass function.

CritPhysDensity 0

The critical physical hydrogen number density in cm^(-3) above which star formation is allowed. In the model of Springel & Hernquist (2003), this is computed via eqn (23) of the this paper if the parameter is set to zero (which is the recommended setting for this model).

CritOverDensity 57.7

If a cosmological simulation would be started at very high redshift, then the physical baryon density could exceed the prescribed physical star formation threshold computed above. To prevent this, a second criterion is imposed, namely to require a minimum comoving overdensity as well, which is given in dimensionless form by this parameter. CritOverDensity=57.7 is the extrapolated overdensity for an NFW halo at the R200 radius, i.e. this choice corresponds to requiring that star-forming gas should be contained inside the virial radius of halos. This cures the high-z problem without restricting low-redshift star formation in any way.

TreecoolFile data/TREECOOL_fg_dec11

This file is used in cosmological simulations to tabulate the time evolution of an externally imposed UV background that is responsible for cosmic reionization. Typically something like a Haard & Madau or Faucher-Giguere model is used. The file tabulates the base10 logarithm of (1+z), followed by the photoionization rates of HI, HeI and HeII, and the associated heating rates. The simulation code will then linearly extrapolate from this table to set the UV background parameters at any given redshift. Only in cosmological simulations with comoving integration the UV background will be used, although the file is always expected if COOLING is enabled.

MinEgySpec 0

This parameter can be used to effectively impose a minimum allowed temperature onto the gas. This is however done in terms of a minimum energy per unit mass u , i.e. if MinEgySpec is set to some finite value, u will not be allowed to drop below this value.

Special features

A_StaticHQHalo 5.0

In case the EXTERNALGRAVITY_STATICHQ option is activated, this specifies the scale length of a static Hernquist halo whose gravitational potential is added to the force calculation. The halo is centered at the origin, and the scale length is given in internal length units.

Mass_StaticHQHalo 100.0

This parameter is only active when EXTERNALGRAVITY_STATICHQ is enabled, and then gives the total mass (in internal units) of the halos that is added as a static potential to the force computation.