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
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.
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
the rank of the file in the group of files representing the
snapshot. If the HDF5 file format is used, an identifier
also appended automatically.
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.
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
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
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
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.
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
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
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.
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
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.
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
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
restart.X files, and then to rename the
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).
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.
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
Simulated time span and spatial extent
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.
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
z=0 . Otherwise the value of
TimeMax refers to
This flag enables or disables comoving integration in an expanding
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
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
PERIODIC is enabled, the code expects the particle
coordinates to lie in the interval
The size of the periodic box (in code units) encompassing the
simulation volume. This parameter is only relevant if the
option is activated.
Cosmological matter density parameter in units of the critical density
z=0. Relevant only for comoving integration.
Cosmological vacuum energy density (cosmological constant) in units of
the critical density at
z=0. For a geometrically flat universe, one
Omega0 + OmegaLambda = 1. Important: For simulations in
Newtonian space that do not account for cosmological expansion, this
parameter has to be set to zero.
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
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.
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
has the same value in all simulations, even if the cosmological
factors and the Hubble constant change), or one sets
unity and uses
Hubble to directly set the Hubble constant. Both is
possible, and intermediate forms in principle as well.
System of units
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.
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.
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
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.
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
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
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
in this example
Gravitational force accuracy
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.
This is the accuracy criterion parameter (the opening angle theta) of
the tree algorithm if the geometric opening criterion
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.
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.
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.
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
Time integration accuracy
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.
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.
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
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
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
A new domain decomposition is not necessarily determined for every
single timestep. A value of
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.
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.
A value of 1 signals that the output times are given in the file
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.
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 Config.sh). Output times that are in
the past relative to the current simulation time will always be
This variable selects the time for the first snapshot (relevant only
OutputListOn=0). For comoving integration, the above choice would
therefore produce the first dump at redshift
OutputListOn=1 this parameter is ignored. Otherwise, after a
snapshot has been written, the time for the next snapshot is
determined by either adding
TimeOfFirstSnapshot, or by multiplying
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.
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
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 +/-
DesNumNgb, the code will readjust the
smoothing length such that the number of neighbours is again in this
This sets the allowed variation of the number of neighbours around the
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.
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.
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
NTYPES (this is however the default).
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
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.
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
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.
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
epsilon_phys at a certain maximum value. These maximum
softening lengths are specified by the
parameters. In the actual implementation, the code uses
min(epsilon_com, epsilon_phys^max / a) as comoving softening. Note
that this feature is only enabled for
SofteningMaxPhysClassX values are ignored. The
SofteningMaxPhysClass0 specifies the maximum
physical softening of the first softening class. Depending on the
NSOFTCLASSES, additional such parameters are needed, one
for each softening class. One hence needs to specify
X = NSOFTCLASSES - 1.
This parameter is only needed if
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.
Specifies the smallest allowed gaseous softening value. Together with
the multiplicative factor
describes the increase from step to step, this defines the discrete
table of available softenings for SPH particles when
ADAPTIVE_HYDRO_SOFTENING is used.
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
NSOFTCLASSES_HYDRO, which is normally set to 64 as
default. (This can be changed, if desired, by overriding the value of
this constant in
Config.sh). The largest available gaseous softening
length is then given by
AdaptiveHydroSofteningSpacing ^ (NSOFTCLASSES_HYDRO - 1).
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
be quite insensitive to the exact choice. The parameter needs only be
SUBFIND is actually enabled in
Initial conditions generation
The following initial conditions parameters are only used if
is activated in the code configuration file. The value of
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
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
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
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.
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.
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:
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).
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.
This parameter is only relevant when
PowerSpectrumType=3 is used
(Efstathiou parameterization). In this case,
usually be set to something close to
Omega0 * HubbelParam.
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.
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.
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
sigma8 value, and the linear theory growth factor for the
specified cosmology. Normally, this renormalization is always
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).
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.
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.
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
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
Config.sh-file, 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
For type=1 (octant), an additional number
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.
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.
Comoving thickness of one lightcone massmap shell.
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, http://adsabs.harvard.edu/abs/2003MNRAS.339..289S).
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.
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.
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.
This is the cloud evaporation parameter A_0 in the above paper.
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.
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).
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
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.
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.
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
u will not be allowed to drop below this value.
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.
This parameter is only active when
enabled, and then gives the total mass (in internal units) of the
halos that is added as a static potential to the force computation.