Custom Fractal Aggregate Aerosols
In this tutorial, we will learn:
How to generate your own .mieff files for custom paticle sizes, numbers of monomers, monomer radii, or fractal dimension values.
The VIRGA v2.0 aerosol database makes the assumption that \(N_{mon}\)=1000 for all species. To explore other options, you need to calculate your own .mieff grids by calculating \(Q_{ext}\), \(Q_{sca}\) and the asymmetry parameter \(g\) for each wavelength, particle size, shape and condensible species that forms in the clouds. We do this using Modified Mean Field theory, which is built into the code optool : https://github.com/cdominik/optool . For more details on MMF calculations in general, see the following papers:
https://doi.org/10.3847/1538-4357/aac32d (Tazaki R. and Tanaka H., 2018)
https://doi.org/10.1093/mnras/stad3743 (Lodge M. et al., 2024)
Optool should be compiled before starting this tutorial (with “make multi=True”, if possible, to massively speed up calculations). See optool documentation for further install instructions.
Once compiled, add the filepaths into the cell below.
[ ]:
# --------------------------------------------- UPDATE THESE FILEPATHS FOR YOUR OWN SETUP ------------------------------------------------------------------
# set the required filepaths
virga_dir = '/PICASO/virga_refractive_indices' # set filepath for VIRGA refractive index files (.refrind). These are usually located wherever you placed the .mieff database
output_dir = '/optical_db' # set filepath for the output .mieff files that you want to create
optool_dir = '/Programs/optool' # set root filepath for main OPTOOL code (code should be compiled already). IMPORTANT: It is assumed that .lnk files are stored within this directory in a "/lnk_data" subdirectory, as standard for installations of optool.
# --------------------------------------------------------------------------------------------------------------------------------------------
[ ]:
# standard import statements
import numpy as np
import virga.justdoit as jdi
Choose gas and define radius grid
You can choose any gas or gasses that you would like to consider as aggregates – they just need to be from the VIRGA database to be able to make clouds out of them. For this example, we will choose KCl.
[ ]:
aggregate_list = ['KCl'] # choose gas or list of gasses (submit as an array of strings)
# optool needs to read the refractive index for each gas, so this function will convert the refractive index file from the virga refractive index database (.refrind format) into a format
# that optool can read (.lnk format), and then save it in the "lnk_data" folder within the optool directory. This only needs to be done once per gas, but it's a very quick function anyway.
jdi.convert_refrind_to_lnk(aggregate_list, virga_dir, optool_dir)
# SET RADIUS GRID -- (in cm!)
radius_min = 1e-8 # in cm, so this is 0.0001 um - set minimum radius of particles to calculate optical properties for
radius_max = 1e-2 # in cm, so this is 100 um - set maximum radius of particles to calculate optical properties for
num_radii = 40 # number of radii that we want to calculate properties for
logspace=True # make the radius grid in log-space (this is the default option)
# the wavelength grid will simply be all of the wavelengths that are in the .refrind file.
Spheres
First, calculate optical properties for spheres on the radius grid made above. VIRGA should print the number of wavelengths and details about the radius bins that are analysed (note: it calculates an average of the optical properties from 6 sub-bins centered around each of these mean radii to smooth out resonance features).
If this process takes a really long time, consider the following:
Check that your radii are in cm. Larger particles take longer computation times, so check your units.
Do you need all of the wavelengths in your .refrind file? If not, you may wish to delete some of them, especially if you are using an experimental dataset that has UV wavelengths in there but not actually analysing them in the clouds (short wavelengths take longer).
Consider decreasing the resolution of your wavelength grid. Fewer wavelengths = Fewer calculations
[ ]:
# ------------------------------------------ SPHERES ------------------------------------------
# First, calculate optical properties for spheres
# Calculate extinction, absorption and scattering efficiencies for the radius grid specified within the function below
optical_properties = jdi.calc_mie_db(aggregate_list, virga_dir, output_dir, optool_dir,
rmin = radius_min, # 0.0001 um
rmax = radius_max, # 100 um
nradii = num_radii, # number of radii that we want to calculate properties for
aggregates=False, # consider condensed gas as spheres
logspace=True)
Aggregates
The code for calculating .mieff files for aggregates is very similar, but we have to provide it with a few more details about the particle shape:
a) The fractal dimension (general shape) of the aerosols
b) EITHER the number of monomers (N_mon) OR the radius of the monomer (r_mon).
Below we create files for 5 different fractal dimensions: Df = [1.2, 1.6, 2.0, 2.4, 2.8] and \(N_{mon}\) = 2,000. To guide you through the creation of the .mieff grid, VIRGA will print each sub-bin that it is currently calculating optical properties for, as well as the exact aggregate properties in each sub-bin. It does calculations for 6 sub-bins, finds the average, and then moves onto the next radius bin (separated by dashed lines). The calculation order is:
Pick a sub-bin radius, beginning with the smallest particles.
Calculate optical properties for all wavelengths at this radius.
Additional guidance
Tip: If optool is compiled with the “make multi=True” command, it will calculate properties for multiple wavelengths on different cores, giving a signficant speed boost (especially if combined with HPC architecture).
If you want to calculate optical properties for fixed \(r_{mon}\) instead of \(N_{mon}\), a second example codeblock is commented out below (see Figure 2 and discussion in Moran & Lodge (2025) for details on both methods).
You may notice that you set \(N_{mon}\) = 2,000 but for the smallest particle in the grid, the N_mon value printed is much smaller. This is because \(r_{mon}\) is not allowed to be smaller than an atom, but VIRGA may still require small particles, so the smallest particles in the grid are allowed to deviate from the fixed N_mon value to remain physical. See Figure 2 and discussion in Moran & Lodge (2025) for more details on this. In each sub-bin, the values of \(N_{mon}\), \(R_{gyr}\), \(k_{0}\) and \(r_{mon}\) are printed for clarity on the exact parameters that MMF is using to calculate optical properties.
If your calculations are taking a long time, check the guidance on reducing calculation time for spheres above. However, note that MMF calculations for large numbers of radii and wavelengths can just take a while.
Finally, check your wavelength ranges in the .refrind files. If the wavelength range spans > 3 orders of magnitude, due to computational memory issues multithreading may be disabled and calculations can be limited to one core at a time, which can slow things down (a warning will be printed if this happens). Avoid this by only submitting wavelengths that you NEED in your .refrind file (and subsequent optool .lnk file).
[ ]:
d_f = np.linspace(1.2,2.8,5) # test 5 different values of fractal dimension
N_monomers = 2000 # Number of monomers in an aggregate
# k0 will be prescribed using Tazeki Eq. 2, if not set by the user
# R_gyr will be calculated within the code
for i_fractal in range(len(d_f)):
# Calculate extinction, absorption and scattering efficiencies for the radius grid specified within the function below
optical_properties = jdi.calc_mie_db(aggregate_list, virga_dir, output_dir, optool_dir,
rmin = radius_min, # 0.001 um
rmax = radius_max, # 0.5 um
nradii = num_radii, # number of radii that we want to calculate properties for
aggregates=True, # consider condensed gas as aggregates
Df=d_f[i_fractal], # loop through each fractal dimension in the array
N_mon=N_monomers, # number of monomers
logspace=True)
'''
# OR if you want to fix r_mon instead of N_mon, use the following code:
d_f = np.linspace(1.2,2.8,5) # test 5 different values of fractal dimension
r_monomers = 1e-6 # Here we provide r_mon in cm, so this would be r_mon = 0.01 um
# k0 will be prescribed using Tazeki Eq. 2, if not set by the user
# R_gyr will be calculated within the code
for i_fractal in range(len(d_f)):
# Calculate extinction, absorption and scattering efficiencies for the radius grid specified within the function below
optical_properties = jdi.calc_mie_db(aggregate_list, virga_dir, output_dir, optool_dir,
rmin = radius_min, # 0.001 um
rmax = radius_max, # 0.5 um
nradii = num_radii, # number of radii that we want to calculate properties for
aggregates=True, # consider condensed gas as aggregates
Df=d_f[i_fractal], # loop through each fractal dimension in the array
r_mon=r_monomers, # monomer radius
logspace=True)
'''