Input file creation
How to make your own input files
FreePATHS is used by providing a config file to the program that controls the simulation. This config file is in the form of a python script. In this python script, you can either define the parameters directly like this NUMBER_OF_PARTICLES = 100, or you can also use python code to set the parameters like this, for example TIMESTEP = 300e-9 / 10000.
Each parameter has a default value, which is defined in the default_config.py file. If you provide no config file, FreePATHS will simply run the default config. And each parameter that you do not set will use the default value. If you input a non-existent parameter, it will be ignored and no error will be raised, so double-check the parameters.
All the parameters are explained on this site, grouped logically. For some parameters, not only what they do but also some theoretical explanations and assumptions behind them are detailed.
The first section explains the most important parameters, the second section presents more in-depth aspects of the simulation, and the last section shows the parameters mainly used for adjusting the simulation outputs.
Although the default values shown on this page are set to realistic values and should work fine if not provided, I still recommend reading this entire page to learn about all of FreePATHS's features and to avoid overlooking anything when setting up your simulation.
Please be advised that the information on this page might not be quite up-to-date.
After learning about the parameters, please look at the example files provided here, some of which have further explanations on the wiki.
Basic parameters
Most basic parameters
These parameters control the most fundamental elements of the simulation.
OUTPUT_FOLDER_NAME = 'Si nanowire at 300 K'
NUMBER_OF_PARTICLES = 10000
TIMESTEP = 2e-12
NUMBER_OF_TIMESTEPS = 200000
T = 300➡️ OUTPUT_FOLDER_NAME : string
The outputs of the simulation will be saved in the Results folder. In this folder, another folder with this name will be created, which will contain the output files. So, in this case, the result files will be in Results/Si nanowire at 300 K. Keep in mind that the the Results folder will be created in the folder you executed the FreePATHS command. Also, pay attention that if the simulation is run again, the results will be overwritten without a warning.
A useful trick is to use f-strings to automatically name the output folders. For example, if the same simulation is to be run at multiple temperatures, using OUTPUT_FOLDER_NAME = f'Simulation at {T}K' will automatically put the simulation temperature in the output folder name. Just make sure to define the parameters (T in this case) before.
➡️ NUMBER_OF_PARTICLES : int
This will define how many particles are simulated. Since the Monte Carlo simulation approach is inherently statistical, more phonons should result in more stable results with less standard variation between the results of different simulations at the cost of more calculation time.
➡️ TIMESTEP : float
The phonons are not simulated in a continuous fashion but only every timestep. If the timestep is small, the phonon behavior will be more realistic, but the simulation time will increase. And vice versa for a large timestep. Because the phonons are only simulated every timestep the time between two scattering events of a particular phonon cannot be smaller than the timestep so take this into account, especially when simulating at high temperatures where scattering events are more frequent. If you experience wrong or unexpected phonon behavior, reducing the timestep can also help with this.
For electron simulations, the default phonon timestep of 2e-12 s is far too large. Electrons travel much faster than phonons, so their mean free path is covered in a single timestep, making scattering statistics meaningless. Use TIMESTEP = 1e-14 s (or smaller) for electron transport simulations.
➡️ NUMBER_OF_TIMESTEPS : int
The phonon is simulated until it reaches a cold side or until the number of timesteps is reached. This is to prevent infinite calculation times if a phonon gets stuck somewhere. Thus, this parameter should usually be high enough so that most phonons reach the cold side. Otherwise, you will see a warning at the end of the simulation. This is especially important for the thermal conductivity calculation: phonons that are removed because they ran out of timesteps are missing from the later timeframes, which biases the temperature and heat flux profiles towards fast phonons.
➡️ T : float
The temperature of the simulation in Kelvin. Since increasing the temperature increases the number of scattering events and phonons take longer to traverse the structure, simulations at high temperatures take significantly longer than at low temperatures.
Simulation domain
The simulation domain consists of a box. These parameters control the size of the box. The unit is meters.
➡️ THICKNESS : float
Defines the thickness of the simulation domain in meters. This corresponds to the z coordinate and the box will span from -THICKNESS/2 to THICKNESS/2 on the z axis.
➡️ WIDTH : float
Defines the width of the simulation domain in meters. This corresponds to the x coordinate and the box will span from -WIDTH/2 to WIDTH/2 on the x axis.
➡️ LENGTH : float
Defines the length of the simulation domain in meters. This corresponds to the y coordinate and the box will span from 0 to LENGTH on the y axis.
Simulation boundaries
Each side of the simulation domain can either be a wall that phonons scatter on, a cold side that phonons disappear on, a hot side where phonons re-thermalize or empty so that phonons can fly through it. The floor and ceiling of the box are always physical walls. By default, the bottom is assumed to be hot, the top is assumed cold, and the left and right walls are just physical walls.
There will not be a section for every parameter in this section because it would be redundant.
To set a wall to be solid, set the corresponding INCLUDE_..._SIDEWALL = True and set COLD_SIDE_POSITION_... = False and HOT_SIDE_POSITION_... = False. The same goes if you want to set a wall to a hot or cold side. If a wall is assigned multiple functions, an error will occur. And if you want a wall to be open, set all three parameters to False.
Particle sources
The phonons are not emitted by the hot sides but by phonon sources. A source of phonons is an area where phonons are generated in a given direction. This area can be placed anywhere in the structure.
➡️ PARTICLE_SOURCES : list
The list needs to contain one or multiple Source objects. Often this will be one source on the hot side. Here is an example (note that all values of a Source that are not set will be zero):
It is convenient to use WIDTH, LENGTH and THICKNESS like in this example:
Note that the source angle distribution is also adjusted. The angle distribution can be chosen from among one of those shown in the image below. In the case of multiple sources, the phonons will be emitted from them with equal probability.
Holes and pillars
While holes are not necessarily required for a simulation, which you can see in the empty default value, they are a key aspect of FreePATHS which is why they are in the necessary parameters section.
➡️ HOLES : list
To build any structure in FreePATHS holes are used. A hole has a certain shape and cuts through the simulation domain in the z direction. A selection of holes and their parameters are shown in the image below. If you want to look at the holes and their parameters in more detail, take a look at the holes.py file. I also recommend taking a look at the all_shapes.py example file.
Note that ParabolaBottom and ParabolaTop are special because you cannot place them anywhere in the simulation domain. They will always appear at the top or bottom side of the structure. See the parabolic_lens_focusing.py example.
To add holes to the simulation, simply put them into the list. It is often very useful to generate these holes using some simple python code. For example, the following code will create a 5x6 square lattice of circular holes:
There are multiple ways to add arbitrary shapes into the simulation. The simplest one is based on mathematical equations and requires practically no programming knowledge. Check out this tutorial.
➡️ PILLARS : list
Pillars are a more experimental feature, and the only pillar available at the time is CircularPillar. Pillars work the same way as holes but instead of preventing phonons from entering a certain area of the simulation domain they extend the simulation domain in z direction locally.
➡️ INTERFACES : list
Interfaces represent vertical planes on which phonon can either pass, or be scattered according to usual rules of scattering on walls. The interface can be set up as follows:
The transmission will be calculated using the equations from this section.
➡️ BULKS : list
Bulks represent rectangular inclusions of another material embedded in the simulation domain (for example a SiGe block inside a Si matrix). Phonon transmission and reflection at the inclusion boundary are calculated with the same spectral hybrid transmission model as interfaces, explained here. See the RectangularBulk class in scatterers.py for the available parameters.
Multiprocessing and resources
➡️ NUMBER_OF_PROCESSES : int
Every phonon is simulated independently, one after the other. To speed up the calculation, the phonons should be distributed across multiple processes, which will each simulate phonons independently. This value should be set to a value close to the number of threads your processor has. Please take note that the progress percentage displayed in the terminal is the progress of a single process, and that some processes will take longer than others to finish.
➡️ LOW_MEMORY_USAGE : bool
When set to true, it will use less memory, which may help with heavy calculations, but will not save the mean free path data.
Advanced simulation parameters
Material
➡️ MEDIA : str
This parameter describes what material the simulation domain is made of. Phonons speed and internal scattering behavior are examples of what is affected by this. Current choices are: Si, SiGe, SiC, and Graphite. Diamond and AlN have tabulated dispersions but no relaxation-time model yet, so they cannot be used for full simulations.
➡️ IS_TWO_DIMENSIONAL_MATERIAL : bool
If this is set to True the z dimension will be ignored and the simulation will take place only in the x-y plane. This is usually used for Graphene sheet simulation.
Roughness
When phonons scatter on a surface, the surface roughness influences the probability of specular or diffuse scattering on the surface.
The roughness values are expressed in meters. Most of the variable names are pretty self-explanatory. To clarify, TOP_ROUGHNESS references the simulation domain boundary in the positive z direction and BOTTOM_ROUGHNESS in the negative z direction. SIDE_WALL_ROUGHNESS are all other simulation domain boundaries.
If the roughness is set to a very low value, the scattering will be mostly specular, which can be used on the side walls, for example, as a symmetric/periodic boundary condition. It can also be used for testing purposes to check if the scattering behavior is programmed correctly.
Internal scattering
➡️ INCLUDE_INTERNAL_SCATTERING : bool
If this is set to False phonons will not experience internal scattering. This means that they will only be scattered on Holes and simulation boundaries. This is mainly used for debugging and testing purposes.
Phonon frequencies and branches are always sampled from the real phonon dispersion of the material, using the density of states weighted by the mode heat capacity and group velocity — the physically correct emission spectrum for a hot reservoir (a flux source emits proportionally to how fast each mode carries energy away). At each inelastic internal scattering event (Umklapp, 4-phonon) the phonon's branch and frequency are redrawn from the same dispersion-weighted distribution, modelling the anharmonic coupling that redistributes phonon energy among modes. Elastic events (impurity/mass-disorder) do not rethermalize — they only redirect the phonon, conserving its frequency and branch.
➡️ USE_DISPERSION_HEAT_CAPACITY : bool
If True (default), the deposited phonon energy is converted to temperature using the heat capacity computed directly from the same dispersion branches that are sampled in the simulation (acoustic branches only). This makes the temperature profile and the resulting thermal conductivity self-consistent with the RTA integral. If False, the full experimental heat capacity (including optical branches) is used, which is physically more accurate for real materials but makes the simulated thermal conductivity incomparable to the model's own RTA prediction.
Grain boundary scattering
For polycrystalline materials, FreePATHS can add an extra scattering channel representing grain boundaries, on top of the usual internal scattering.
➡️ GRAIN_SIZE : float or None
The mean grain diameter, in meters. Leave as None (default) to disable grain boundary scattering entirely.
➡️ GRAIN_SIZE_STD : float
The standard deviation of the grain size, in meters. Each phonon is assigned a grain size drawn from a lognormal distribution with this mean and standard deviation. Set to 0 for monodisperse grains (all the same size).
➡️ GRAIN_ROUGHNESS : float
The RMS disorder width of the grain boundary, in meters, used in a Soffer-type specularity factor: at low frequencies (long wavelengths) the boundary appears smooth and phonons pass through, while at high frequencies it acts as a diffuse scatterer. Typical values are 100 nm–10 µm for GRAIN_SIZE and 0.5–2 nm for GRAIN_ROUGHNESS.
Time
Do not confuse the "virtual" timesteps discussed in this section with the timesteps discussed in the Most basic parameters section. The basic NUMBER_OF_TIMESTEPS parameter defines the maximum time a phonon has to travel through the structure, while the parameters of this section are used to make sure the thermal simulation can reach the steady state.
Considering the Thermal map.pdf and the resulting Temperature profile.pdf please consider that the physics of the entire simulation behaves with the temperature of the parameter T even if Temperature profile.pdf shows a significantly higher temperature. This is because the temperature in Temperature profile.pdf results from the amount of heat that enters the structure, which is dependent on NUMBER_OF_PARTICLES. Thus, the temperatures in Temperature profile.pdf should not be taken at face value. For the thermal conductivity calculation, the gradient of this profile is used, as explained here.
➡️ NUMBER_OF_VIRTUAL_TIMESTEPS : int
The phonons do not all enter the structure at the same time, but a virtual start time is assigned to each phonon randomly, and the range of these start times is controlled with this parameter. In other words, this parameter defines the total virtual time span of the thermal simulation: phonon emission is spread uniformly over it, and it is this time span that is divided into the timeframes described below. Because no phonons are generated before the simulation starts, the first moments of the simulation are not useful, as all phonons are at the beginning of the structure and none are towards the end — this initial period is discarded via NUMBER_OF_STABILIZATION_TIMEFRAMES. So this parameter should be at least a couple of times larger than the time it takes phonons to traverse the structure, so that the steady state can establish before the measurement timeframes begin. The time it takes phonons to traverse the structure can be determined with Distribution of travel times.pdf (determining the 95% or 99% quantile by eye should be sufficient). Note also that each phonon is only simulated for at most NUMBER_OF_TIMESTEPS of its own flight time, so if the virtual time span is much longer than that, make sure that almost all phonons reach the cold side within their flight time — otherwise the later timeframes will be missing the slow phonons and the profiles will be biased.
➡️ NUMBER_OF_TIMEFRAMES : int
The simulation time determined by NUMBER_OF_VIRTUAL_TIMESTEPS is divided into several timeframes to observe the evolution of the system with time. The number of the timeframes should be around 8 for reasonable output.
➡️ NUMBER_OF_STABILIZATION_TIMEFRAMES : int
The system requires some time to reach a steady state, in which the thermal conductivity should be calculated. Thus, this parameter controls how many first timeframes are skipped before the measurement begins. For example, this can be at least half of the NUMBER_OF_TIMEFRAMES parameter.
Electron parameters
➡️ IS_CARRIER_ELECTRON : bool
Choose whether the charge carriers are electrons or holes.
➡️ ELECTRON_MFP : float
Set the mean free path of the carriers. Typically it is a few nanometers. This is the crystal (acoustic-phonon-limited) mean free path, taken energy-independent. When DOPING_CONCENTRATION > 0 it is combined with the ionized-impurity mean free path (see below) through Matthiessen's rule.
➡️ DOPING_CONCENTRATION : float
Ionized-dopant concentration N_I (in m⁻³, i.e. cm⁻³ × 10⁶). When greater than zero, the carrier mean free path becomes energy- and doping-dependent through the Brooks–Herring ionized-impurity model, combined with ELECTRON_MFP via Matthiessen's rule: 1/Λ(E) = 1/ELECTRON_MFP + 1/Λ_ii(E). Ionized impurities scatter low-energy carriers most strongly (τ_ii ∝ E^{3/2}), which shortens the mean free path with doping and raises the Seebeck coefficient. This requires the material to define a static dielectric constant (currently silicon). The default 0.0 recovers the constant ELECTRON_MFP.
➡️ SURFACE_POTENTIAL : float
Surface band bending φ_s (in eV) used to model surface depletion of carriers. When greater than zero, a dead layer of width W = sqrt(2·ε_s·φ_s / (e·N)) (the abrupt-junction depletion approximation) is excluded around every free surface — side walls, top/bottom, and hole edges — for carriers only (phonons are unaffected). This represents Fermi-level pinning by surface states or the native oxide, which depletes mobile carriers near surfaces and reduces the conductivity of nanostructures. φ_s is the total band bending set by the pinning position (typically a few tenths of an eV). It needs a doping (DEPLETION_DOPING or DOPING_CONCENTRATION) and a material dielectric constant. The default 0.0 disables depletion. Note that a neck between two holes is depleted from both sides, so the conducting slit is neck − 2W; if 2W exceeds the neck it pinches off entirely (σ → 0).
➡️ DEPLETION_DOPING : float
Doping N (in m⁻³) used only for the surface-depletion width formula above. When 0.0/None it falls back to DOPING_CONCENTRATION. Set this explicitly (and leave DOPING_CONCENTRATION = 0) to add a surface dead layer driven by the real doping without engaging Brooks–Herring ionized-impurity scattering — i.e. to isolate the depletion effect on top of a constant-ELECTRON_MFP transport calibration (so that a fixed MEAN_MAPPING_CONSTANT obtained without depletion stays valid).
➡️ MEDIA_FERMI_LEVEL : float or None
The Fermi level of the material in joules, measured from the conduction band minimum (i.e. negative values mean the Fermi level is below the band edge). If set to None, the material's default Fermi level is used. Setting this explicitly is important for doped semiconductors, where the Fermi level depends on the carrier concentration. For example, for n-doped silicon at 3×10¹⁸ cm⁻³, a value of -50e-3 * electron_volt places the Fermi level 50 meV below the conduction band. This value is used to mark the operating point on the output plots of conductivity, Seebeck coefficient, and power factor.
➡️ MEAN_MAPPING_CONSTANT : float or None
This is the calibration constant C that maps the raw MC travel-time results to physically correct units. It must be obtained from a two-step workflow:
Step 1 — calibration run: Set MEAN_MAPPING_CONSTANT = None and run the simulation on a pristine structure with no holes or other scatterers. FreePATHS will automatically compute C and save it to Data/Mapping constant.csv. Note the value at your material's Fermi level.
Step 2 — production runs: Set MEAN_MAPPING_CONSTANT to the value obtained in step 1 (e.g. MEAN_MAPPING_CONSTANT = 5e-6) and run the simulation on your nanostructured geometry. FreePATHS will use this fixed C for all subsequent calculations.
If you skip step 1 and leave MEAN_MAPPING_CONSTANT = None while running on a nanostructured sample (with holes, pores, etc.), FreePATHS will recompute C from the nanostructured travel times, which conflates the geometry effect with the calibration and gives incorrect results. The value 5e-6 shown in the example is only a placeholder — you must determine it for your specific material and temperature. See this page for more details on C.
➡️ FERMI_LEVEL_LOWER_BOUND / FERMI_LEVEL_UPPER_BOUND : float
These define the range of Fermi levels (in eV, measured from the conduction band minimum) over which the post-processing integrals for σ, S, and PF are evaluated. The default range of −0.2 to +0.1 eV covers the intrinsic and lightly-doped regimes at 300 K. Extend the lower bound (e.g. to −0.4 eV) if you need to cover near-intrinsic or p-type conditions, or raise the upper bound for heavily n-doped samples. Note that this is a post-processing sweep only — it does not affect which electrons are simulated.
Output parameters
Thermal maps
The heat flux maps, the thermal map and the pixel volumes plot are all based on the same pixel grid. The thermal conductivity calculation is heavily dependent on these plots.
Concerning the heat flux maps, it is important to consider that the Heat flux map.pdf file displays the absolute magnitude of heat flux. So if a phonon travels through a place twice in opposite directions, the heat flux will be added. In the Heat flux map x.pdf and Heat flux map y.pdf the directional heat flux is calculated which means that if a phonon travels through a place twice in opposite directions, the heat flux will cancel out.
➡️ NUMBER_OF_PIXELS_X NUMBER_OF_PIXELS_Y : int
These parameters define the pixel grid which is used for the map creation, so a higher number will generally result in higher quality maps. Note that this is not purely for looks, but the thermal conductivity calculations rely on the values in these maps. I recommend setting these using this code snippet, where the pixel size can be adjusted:
➡️ IGNORE_FAULTY_PARTICLES : bool
Sometimes, particles may escape the structure and get trapped outside the structure or travel outside the simulation domain. This can cause the maps to not look nice because the holes are not empty. If this parameter is set to True particles that are outside the structure are simply ignored in the map generation, which improves both the look and the calculations. Try keeping this on False so that you notice if particles leave the structure, which can be an indicator of some errors or bugs. If just a few particles leave the structure, this can be turned on to still get clean data.
➡️ GRADIENT_FIT_RANGE : tuple
This parameter defines the portion of the structure length, as a pair of fractions between 0 and 1, over which the temperature gradient is fitted and the heat flux is averaged in the thermal conductivity calculation. Within a few phonon mean free paths of the hot and cold contacts the transport is quasi-ballistic, so the temperature profile deviates from linear near the contacts (temperature jumps, similar to those next to thermostats in molecular dynamics simulations). Restricting the fit to the interior of the structure excludes these contact regions, which yields a more accurate thermal conductivity, especially when the structure length is not much larger than the phonon mean free paths. By default, (0.1, 0.9), the outer 10% of the length at each contact is excluded; set it to (0.0, 1.0) to use the whole length, or to something like (0.2, 0.8) when the contact regions are more pronounced. The linear fit shown in Temperature profile.pdf uses the same range.
Structure plots
➡️ OUTPUT_SCATTERING_MAP : bool
If this is set to True an additional output file Scattering map.pdf will be generated in which the position and type of all scattering events is shown. This can be very useful for debugging and testing.
➡️ OUTPUT_TRAJECTORIES_OF_FIRST : int
This parameter defines how many trajectories of phonons are saved into the Particle paths.csv output file and how many phonon trajectories are plotted in the Particle paths XY.pdf and Particle paths YZ.pdf output files.
➡️ OUTPUT_STRUCTURE_COLOR : str
This variable defines the background color in the Particle paths XY.pdf and Particle paths YZ.pdf output plots. You can use a hexadecimal color or a color name like 'blue' or any color that the matplotlib library accepts.
Line plots
➡️ NUMBER_OF_NODES : int
This value affects the number of bins for the histogram output plots, like for example Distribution of angles.pdf.
➡️ NUMBER_OF_LENGTH_SEGMENTS : int
A few plots display information in segments along the y-axis like Scattering rate profile.pdf and Time spent in segments.pdf. The number of segments for these plots can be adjusted with this parameter.
Animation
Animations are currently not working. See this issue.
Note that NUMBER_OF_TIMESTEPS should not be too large, otherwise the generation of animation may take a very long time because one frame for each time step will be created. A few hundred time steps is a reasonable value.
➡️ OUTPUT_PATH_ANIMATION : bool
Set this to True to generate an animation.
➡️ OUTPUT_ANIMATION_FPS : int
Each timestep corresponds to one frame. This parameter determines the playback speed of the frames in the generated video. Please note that all particles will start flying at the beginning of the video.
Last updated