The discrete dipole approximation for simulation of light scattering by
  particles much larger than the wavelength
The discrete dipole approximation for simulation of light scattering by particles much larger than the wavelength
The discrete dipole approximation for simulation of light scattering 
by particles much larger than the wavelength 
M.A. Yurkina,b∗, V.P. Maltsevb,c, and A.G. Hoekstraa∗
a Section Computational Science, Faculty of Science, University of Amsterdam, Kruislaan 
403, 1098 SJ, Amsterdam, The Netherlands 
b Institute of Chemical Kinetics and Combustion, Siberian Branch of the Russian Academy of 
Sciences, Institutskaya Str. 3, 630090, Novosibirsk, Russia 
c Novosibirsk State University, Pirogova Str. 2, 630090, Novosibirsk, Russia 
Abstract 
In this manuscript we investigate the capabilities of the Discrete Dipole Approximation 
(DDA) to simulate scattering from particles that are much larger than the wavelength of the 
incident light, and describe an optimized publicly available DDA computer program that 
processes the large number of dipoles required for such simulations. Numerical simulations of 
light scattering by spheres with size parameters x up to 160 and 40 for refractive index 
 and 2 respectively are presented and compared with exact results of the Mie theory. 
Errors of both integral and angle-resolved scattering quantities generally increase with m and 
show no systematic dependence on x. Computational times increase steeply with both x and 
m, reaching values of more than 2 weeks on a cluster of 64 processors. The main distinctive 
feature of the computer program is the ability to parallelize a single DDA simulation over a 
cluster of computers, which allows it to simulate light scattering by very large particles, like 
the ones that are considered in this manuscript. Current limitations and possible ways for 
improvement are discussed. 
05.1=m
Keywords: discrete dipole approximation, light scattering simulation, computer program 
                                                 
∗ Corresponding authors. Tel.: +31-20-525-7462; fax: +31-20-525-7490. 
  e-mail addresses: myurkin@science.uva.nl, alfons@science.uva.nl 
1 Introduction 
The discrete dipole approximation (DDA) is a general method to calculate scattering and 
absorption of electromagnetic waves by particles of arbitrary geometry and composition. The 
DDA was first proposed by Purcell and Pennypacker [1] and was reviewed by Draine and 
Flatau in 1994 [2]. A recent review [3] describes the current state of the DDA and its 
historical development. It also explains the equivalence of the DDA and methods based on the 
volume integral equation formulation. The reader is referred to this review for an in-depth 
discussion of the DDA. 
There are a number of computer programs based on the DDA, some of which were 
recently compared by Penttila et al. [4]. The most popular among them is DDSCAT [5], 
which has been widely used by many researchers for more than 10 years. In this paper we 
present a new program, Amsterdam DDA (ADDA), which recently has been put in the public 
domain.1 Its main distinctive feature is the ability to parallelize a single DDA simulation over 
a cluster of computers, which allows simulation of light scattering by very large particles. 
This is demonstrated for a number of test cases in this manuscript. Validation of ADDA by 
simulating light scattering by wavelength-sized particles and comparing it with other DDA 
programs was reported elsewhere [4]. 
Section 2 describes in detail the ADDA computer code, showing its advantages 
compared to other codes. A number of numerical tests are shown in Section 3, demonstrating 
that DDA is actually capable processing large particles, and showing the current capabilities 
of ADDA. Results of these simulations are discussed in Section 4; the errors are compared 
with previous results for much smaller particles. Section 5 concludes the manuscript and 
discusses possible future work. 
2 ADDA computer code 
ADDA has been developed over a period of more than 10 years at the University of 
Amsterdam [6-8]. Its main feature (distinctive from other DDA codes) has always been the 
capability of running on a cluster of computers, parallelizing a single DDA computation, in 
contrast with e.g. DDSCAT [5] that allows farming several instantiations of a DDA 
simulation to different processors. This allows using a practically unlimited number of 
dipoles, since ADDA is not limited by the memory of a single computer [8,9]. Recently the 
overall performance of the code has been improved significantly, together with some 
optimizations specifically for single-processor mode. ADDA's source code and 
documentation is freely available.
Most of ADDA is written in ANSI C, which ensures wide portability on the source-code 
level. The code is fully operational under Linux and, in sequential mode, on Windows based 
systems. The parallelization over multiple processors is based on a geometric decomposition 
of the particle and the single-program-multiple-data paradigm of parallel computing. The 
code is written for distributed memory systems using the message passing interface (MPI).2 
Note that ADDA should in principle also run on shared memory computers, but so far this 
was not explicitly tested. The fast Fourier transform (FFT) used for the matrix-vector products 
in the iterative solver is performed either using routines by Temperton [10] or the more 
advanced package “Fastest Fourier transform in the West” (FFTW) [11]. The latter is 
generally considerably faster but requires a separate package installation. 
ADDA has four options implemented for dipole polarizabilities: Clausius-Mossotti [1], 
radiative reaction correction [12], lattice dispersion relation (LDR) [13], and corrected LDR 
[14]. It includes four iterative methods: conjugate gradient applied to normalized equation 
with minimization of residual norm (CGNR) [15], Bi-CG stabilized (Bi-CGSTAB) [15], Bi-
                                                 
1 http://www.science.uva.nl/research/scs/Software/adda/
2 http://www.mpi-forum.org
http://www.science.uva.nl/research/scs/Software/adda/
http://www.mpi-forum.org/
CG [16], and quasi minimal residual (QMR) [16]. The last two iterative methods employ the 
complex-symmetric property of the DDA interaction matrix to halve the calculation time [16]. 
The default stopping criterion of the iterative method in ADDA is the relative norm of the 
residual ε, which must be . 510−<
The usual formulation of DDA can be written as [2,3]:  
jijii EPGP =− ∑
−α , (1)
where iα  is the tensor of dipole polarizability,  is incident electric field, 
iE ijG  is the free-
space Green’s tensor (complex symmetric), and Pi is the unknown dipole polarization. If the 
polarizability tensor is diagonal for all dipoles then there always exists a iβ  such that 
iii αββ = , i.e. ii αβ = . Moreover, iβ  is then complex symmetric, and so is the matrix with 
elements 
A , (2)
where I  is an identity tensor. A  is the interaction matrix that is used in ADDA, i.e. the 
following system of linear equations is solved: 
jjijii
jij ExGxxA βββ =−= ∑∑
, (3)
where iii Px
1−= β  is a new unknown vector. Eq. (3) is equivalent to the use of Jacobi-
preconditioning [15] together with keeping the interaction matrix complex-symmetric (for any 
distribution of refractive index inside the scatterer and for any of the supported polarization 
prescriptions). We have not studied, however, whether this Jacobi-preconditioning improves 
the convergence of the iterative solver. Flatau showed [17] that in some test cases it helps, 
while in others there is no improvement. It is important to note also that DDA is not limited to 
diagonal or symmetric polarizabilities. Any other tensor may be used, but then the interaction 
matrix is not complex-symmetric; hence, QMR and Bi-CG are less efficient. 
ADDA can perform orientation averaging of the scattering quantities over three Euler 
angles (α, β, γ) of the particle orientation. Averaging over the angle α is done with a single 
computation of internal fields by computing scattering in different scattering planes, which is 
comparably fast. Averaging over the other two Euler angles is done by independent DDA 
simulations. The averaging itself is performed using a Romberg integration [18], which may 
be used adaptively (i.e. automatically simulating the required number of different orientations 
to reach a prescribed accuracy) but limits the possible number of values for each orientation 
angle to be , where n is an integer. Moreover, symmetries of the scatterer may be used 
to decrease the intervals of Euler angles, over which to average, and hence accelerate the 
calculation. This feature of ADDA was tested in a recent benchmark study [4]. 
12 +n
Other features of ADDA include computation of scattering by a tightly focused 
Gaussian beams [6], a checkpoint system to allow for long runs on queuing systems that 
enforce upper limits on wall clock time for execution as is usually the case on massively 
parallel supercomputers, calculation of radiation forces on each of the dipoles [19], use of 
rotational symmetry of the scatterer to halve the simulation time, and an extended command 
line interface. Some other features, such as applicability to anisotropic scatterers and a large 
set of predefined shapes, are planned to be implemented in the near future. 
There are several factors that allow ADDA's performance to compare favorably with 
other codes, which was shown in a benchmark study by Penttila et al. [4]. First of all, the 
FFTW 3 package that is used automatically adapts itself to optimally perform on any 
particular hardware. Moreover, ADDA does not perform complete 3D FFT transforms in one 
run, but decomposes them into a set of 1D transforms with data transposition in between. This 
allows employing the fact that input data for the forward transform contains many zeros, and 
0 20 40 60 80 100 120 140 160
ε ∈(10−5,10−3)
Size parameter x
ε =10−5
70 GB
Fig. 1. Current capabilities of the ADDA for spheres with different x and m. The striped region 
corresponds to full convergence and densely hatched region to incomplete convergence. The dashed 
lines show two levels of memory requirements for the simulation, according to the “rule of thumb” 
(see main text for explanation). 
only part of the output data of the backward transform is used [8]. Second, we have 
implemented four different Krylov-space-based iterative solvers, allowing us to choose the 
most suitable one for a particular application. As is known from the literature [17,20,21] and 
demonstrated in Section 3, there is not a best iterative solver for DDA. Depending on all 
details of the scattering problem, any of the methods may outperform the others. Third, 
dynamic memory allocation and optimized data structures allow all computations, except the 
FFT, to be performed only for the real (non-void) dipoles and not for the whole computational 
box. This also decreases ADDA's memory consumption. Moreover, symmetry of the 
interaction matrix is used to decrease memory required for its Fourier transform. Finally, all 
float variables in ADDA are represented in double precision. This accelerates convergence in 
cases when machine precision becomes important. Moreover, basic operations with double-
precision numbers can be faster than with single-precision ones on modern processors. This 
acceleration comes at a cost of increased memory consumption, which is, however, still lower 
than for other computer codes [4]. 
More information on ADDA can be found in an extensive manual included in the 
distribution package.
3 Numerical simulations 
3.1 Simulation parameters 
In our tests we used ADDA v.0.75, compiled with the Intel C compiler v.9.0 with maximum 
possible optimizations (default options in ADDA’s makefile). All the tests were run on the 
Dutch compute cluster LISA,3 using 32 nodes (each dual Intel Xeon 3.4 GHz processor with 
4 GB RAM). LDR was used as the most common polarization formulation. We have tried 
three different iterative solvers: QMR, Bi-CG, and Bi-CGSTAB. For all of them a default 
stopping criterion  was used. 510−=ε
                                                 
3 http://www.sara.nl/userinfo/lisa/description/
http://www.sara.nl/userinfo/lisa/description/
Table 1. Parameters of the numerical simulations. 
m x λ/md 
Number of 
dipolesa Iterative method 
Number of 
iterations 
20 9.6 2.6×105 Bi-CGSTAB 6 
30 9.6 8.8×105 Bi-CGSTAB 7 
40 9.6 2.1×106 Bi-CGSTAB 9 
60 9.6 7.1×106 Bi-CGSTAB 14 
80 9.6 1.7×107 Bi-CGSTAB 20 
100 9.6 3.3×107 Bi-CGSTAB 27 
130 10.3 9.0×107 Bi-CGSTAB 40 
1.05 
160 9.6 1.3×108 Bi-CGSTAB 65 
20 10.5 5.1×105 QMR 86 
30 11.2 2.1×106 QMR 223 
40 10.5 4.1×106 QMR 598 
60 9.8 1.1×107 QMR 2120 
80 10.5 3.3×107 Bi-CGSTAB 21748 
100 10.1 5.7×107 Bi-CGSTAB 6169 
130 10.3 1.3×108 Bi-CGSTAB 29200 
20 10.8 8.8×105 QMR 1344 
30 10.8 3.0×106 QMR 16930 
40 10.8 7.1×106 QMR 8164 
60 9.6 1.7×107 Bi-CG 127588 
20 11.0 1.4×106 QMR 8496 1.6 
30 10.5 4.1×106 Bi-CG 69748 
20 11.2 2.1×106 QMR 28171 1.8 
30 10.2 5.5×106 Bi-CG 118383 
2 20 10.1 2.1×106 QMR 58546 
a This is the total number of dipoles in the rectangular computational grid, which is the main factor determining 
the computation time of one iteration. For spheres the number of dipoles occupied by the scatterer itself is almost 
two times smaller. 
0 5000 10000 15000 20000 25000
Number of iterations  
Fig. 2. Convergence of the QMR iterative solver for the sphere with x = 20 and m = 1.8. The residual 
as a function of the iteration number is shown. The system of linear equations contains 3×106 
unknowns. 
Spheres were used as test objects. Their size parameter x was varied from 20 to 160 and 
their refractive index m was varied from 1.05 to 2. We limited ourselves to the case of real m. 
The current capabilities of ADDA are shown as a region of the (x,m)-plane in Fig. 1. The 
striped region corresponds to full convergence, the densely hatched region corresponds to 
those cases where ADDA could not fully converge to the required residual norm, but only to 
20 40 60 80 100 120 140 160
106 1 week
1 day
Size parameter x
          m =
 1.05
1 min
1 hour
Fig. 3. Total simulation wall clock time (on 64 processors) for spheres with different x and m. Time is 
shown in logarithmic scale. Horizontal dotted lines corresponding to a minute, an hour, a day, and a 
week are shown for convenience. 
20 40 60 80 100 120 140 160
Size parameter x
          m =
 1.05
Fig. 4. Relative errors of the extinction efficiency in logarithmic scale for spheres with different x and 
)10,10( 35 −−∈ε . Although this incomplete convergence probably affects the final accuracy of 
the scattering quantities only slightly, we remove such results from further consideration 
because a separate study is required to quantify this effect (see Section 4). For fully converged 
results, the errors of scattering quantities due to the numerical convergence are much smaller 
than the total errors (data not shown). 
A complete set of (x,m) pairs, for which ADDA converged, is shown in Table 1. It also 
shows the number of dipoles per wavelength in the medium ( md/λ  where d is the size of the 
dipole). We tried to keep it equal to 10 according to the “rule of thumb” as formulated by 
Draine and Flatau [2]; however, it was slightly different because we varied the size of the 
dipole grid to optimize the parallel efficiency of ADDA.4 The total number of dipoles in a 
rectangular computational grid, shown in Table 1, was varied from 2.6×105 to 1.3×108, it can 
                                                 
4 The best parallel performance is obtained when grid size divides the number of processors. However, ADDA 
works with any grid size. 
20 40 60 80 100 120 140 160
Size parameter x
          m =
 1.05
Fig. 5. Same as Fig. 4 but now for the asymmetry parameter. 
20 40 60 80 100 120 140 160
Size parameter x
          m =
 1.05
Fig. 6. Maximum relative errors of S11(θ ) in logarithmic scale for spheres with different x and m. 
be approximately determined as . Both memory requirements and computation time 
of one iteration are proportional to this number. Two dashed lines are shown in 
3)18.3( xm
Fig. 1 to 
indicate the memory requirements for different x and m. They correspond to typical memory 
of a modern desktop computer (2 GB) and the maximum total memory used in our 
simulations (70 GB), respectively. 
For each sphere we computed the extinction efficiency, the asymmetry parameter, and 
all Mueller matrix elements in one scattering plane, which is a symmetry plane of the cubical 
discretization of the sphere. Exact results for the same spheres were obtained using the Mie 
theory [22]. Spherical symmetry was used by ADDA to get all results from calculations for 
only one polarization state of the incident field. Therefore computation time is a factor of two 
smaller than for non-symmetric scatterers with the same x and m. We employed a volume 
correction to ensure equal volumes of sphere and its dipole representation [2]. Note, however, 
that for the very large spheres this correction is extremely small. 
3.2 Results 
Table 1 shows the iterative solver that provided the best performance for each particular case 
and the number of iterations to achieve convergence. Fig. 2 illustrates one specific example of 
20 40 60 80 100 120 140 160
Size parameter x
          m =
 1.05
Fig. 7. Same as Fig. 6 but now for RMS relative errors. 
0 30 60 90 120 150 180
170 175 180
Scattering angle θ, deg
Fig. 8. DDA results (dotted line) of S11(θ ) in logarithmic scale for a sphere with x = 160 and m = 1.05, 
compared with the results of the Mie theory (solid line). 
convergence of the DDA iterative solver. This is QMR applied to the system of 3⋅106 linear 
equations obtained for the sphere with 20=x  and 8.1=m . The total simulation wall clock 
time t for all particles is shown in Fig. 3. Fig. 4 and Fig. 5 show the relative errors of the 
extinction efficiency Qext and the asymmetry parameter >< θcos  respectively. Maximum - 
and root-mean-squared (RMS) relative errors of S11 over the whole range of scattering angle 
are shown in Fig. 6 and Fig. 7 respectively. Errors of other non-trivial Mueller matrix 
elements behave in a similar way (data not shown). 
DDA results of S11(θ) for a sphere with 160=x  and 05.1=m  are compared with the 
Mie theory in Fig. 8. The inset shows a magnification of the backscattering region. This is, to 
the best of our knowledge, the largest particle ever simulated with DDA. Fig. 9 and Fig. 10 
show the same comparisons but for 60=x , 4.1=m  and 20=x , 2=m  respectively. 
0 30 60 90 120 150 180
Scattering angle θ, deg
Fig. 9. Same as Fig. 8 but now for x = 60 and m = 1.4. 
0 30 60 90 120 150 180
Scattering angle θ, deg
Fig. 10. Same as Fig. 8 but now for x = 20 and m = 2. 
4 Discussion 
The convergence of the QMR iterative solver shown in Fig. 2, featuring plateaus and steep 
descents, is in agreement both with its behavior in general [16] and with particular examples 
of its application to DDA [20,23]. A distinctive feature of this graph compared to the 
literature data is that the convergence slows down with iteration number, i.e. the logarithm of 
the residual norm decreases slower than linearly. This is probably due to the large size of the 
scatterer and loss of numerical precision (see discussion below). 
The total computation times t increase steeply both with x and m (Fig. 3). The time is 
displayed in a logarithmic scale covering a range from 4 seconds to more than 2 weeks. For 
, the increase of t with x is mostly due to the increasing number of dipoles to model 
the scatterer, since the number of iterations increase at a slower pace (
05.1=m
Table 1). For larger m 
these two effects are comparable, combining into a very unfavorable scaling, which can be 
approximately described by a power law , where )()( mxmCt α≈ 6>α  for . It should 
be noted that both the number of iterations and t do not always increase monotonically with x. 
For example for ,  and 
2.1≥m
80=x 2.1=m 30=x , 4.1=m  the execution times are unusually high. 
This may be caused by a large condition number of DDA interaction matrices for these two 
particular particles. Moreover, when the convergence is slow it may suffer from machine 
precision, the latter determining the limit of x and m, for which ADDA will converge at all. 
Therefore, current size limitations of the DDA for  are due to the practically 
unbearable computation times, and not due to memory requirements.
2.1≥m
5 Simulations for larger 
m are far from the memory limit shown in Fig. 1. Moreover, simply using more processors 
does not solve the problem. Improving numerical performance is required, e.g. dedicated 
preconditioning of the iterative solver [15]. On the other hand, extension to larger sizes for 
 is feasible if more computer resources are available. This facilitates, for example, 
simulating scattering of visible light by almost all biological cells in suspension. 
2.1<m
The increase of the number of iterations with m is a well-known fact [12,17,21,24]; 
however, there is still no theoretical foundation to describe it in details. Rahola [24] provided 
theoretical predictions of the dependence of the number of iterations on m, valid for scatterers 
smaller than the wavelength. However, these conclusions are not applicable to the scatterers 
studied in this manuscript. The general reason for the slowing down of the convergence with 
increasing m is increased interaction between dipoles and, hence, an increased condition 
number of the interaction matrix. Absorption, if present, should decrease the overall 
interaction between dipoles in a large scatterer. Therefore, it is expected that convergence for 
complex refractive indices should be better than for the purely real ones that we consider here. 
The same was suggested by Budko and Samokhin [25] based on the analysis of the spectrum 
of the integral scattering operator. However, this proposition is still to be verified by 
numerical tests. 
Another parameter that may greatly affect the computation time is the convergence 
threshold ε. In this paper it is set to a de-facto default value of 10-5 [2], which ensures 
negligibly small numerical errors compared to the model errors. However, in many cases 
numerical errors are small enough already for , i.e. the difference of the scattering 
quantities between simulations with  and  is significantly smaller than the 
difference between the latter and the exact values (data not shown). 
310−=ε
310−=ε 510−=ε
Fig. 2 shows that QMR 
for a particular case converges to  and  three and six times faster 
respectively than to . Results for other simulated particles and iterative solvers show 
similar trends and even larger acceleration with increasing ε in some cases (data not shown). 
Therefore, if one can determine an optimum ε for a particular case, it can decrease the 
computation time significantly. However, we do not pursue this issue further in this 
manuscript. 
310−=ε 3102 −×=ε
510−=ε
Fig. 4 shows the deterioration of the accuracy of Qext with increasing m, while there is 
no clear dependence on x (the only exception is a single result for ). Results for 2=m
>< θcos  (Fig. 5) behave in a similar way. These results are in good agreement with results of 
other researchers for smaller size parameters [2,13,26], both in terms of the errors themselves 
and their dependence on m. To express errors on the angular dependencies of S11 we use two 
integral parameters: the maximum - and RMS relative errors (Fig. 6 and Fig. 7 respectively). 
Although these parameters are not completely objective, as they are significantly influenced 
by the values of S11 in deep minima, which are completely irrelevant to most real 
experiments, they do provide a consistent measure of the DDA accuracy. To relate these 
integral parameters to some other criteria, e.g. visual agreement, three examples are presented 
in Fig. 8 – Fig. 10. Errors of S11(θ ) show the same tendencies as the integral scattering 
quantities, except that errors for  are relatively large (larger than those for  in 
the range ) and generally decrease with x. This is due to the relative nature of the 
measured errors and the huge dynamical range of S
05.1=m 2.1=m
11(θ ) for small refractive indices (see Fig. 
8). Results for smaller size parameters found in the literature [2,26] show a similar increase of 
                                                 
5 The boundary value of m is not well-defined, as it depends on particular hardware and restrictions on 
computation time; 1.2 is just a convenient value to guide the reader. 
errors with m: however, the errors themselves are considerably smaller. For instance, 
maximum relative errors of S11(θ ) for 10<x  and m up to i4.15.2 +  are smaller than 0.4. This 
is due to the general differences between functions S11(θ ) for particles comparable to and 
much larger than the wavelength. The latter has deeper minima and a larger overall dynamic 
range. It is important to note that refractive indices as small as 1.05 are rarely used in DDA 
simulations [26], therefore it is hard to make any definite conclusions concerning the behavior 
of errors in this case. 
In what follows, the traditional “rule of thumb” [2] is discussed, which states that for 
10/ =mdλ  errors of cross sections and asymmetry parameter are expected to be a few 
percents, and maximum errors in the angular dependence of S11 on the order of 20-30 %. 
Results for both Qext and >< θcos  do satisfy the “rule of thumb,” however this rule does not 
describe the decrease of errors by two orders of magnitude with decreasing m. The latter can 
be used to cut down the number of dipoles and hence computation time in cases when only 
integral scattering quantities need to be calculated for small m. Relative errors of S11(θ ) are 
much larger than that predicted by the “rule of thumb,” which is due to the fact that the latter 
was derived based on test simulations for x smaller than 10 [2]. See, however, the discussion 
below on possible changes for complex refractive index and non-spherical shapes. To 
conclude, the “rule of thumb” has very limited application for the range of x and m here. More 
elaborate empirical functions are required to estimate the number of dipoles needed to reach a 
prescribed accuracy. They will also allow a more realistic estimate of DDA computational 
complexity, i.e. the computation time needed to reach a certain accuracy of some scattering 
quantities for particular x and m. This topic is left for the future study. 
The test results shown in this paper are limited to real refractive indices and spherically 
shaped scatterers. In the following we try to generalize our conclusions to complex refractive 
index and non-spherical shapes. However, we want to stress that this generalization is 
speculative, and more numerical tests are clearly needed to verify them. It is expected that 
accuracies of integral scattering quantities should not change significantly for more general 
cases. Their accuracy should deteriorate both with increasing real and imaginary parts of the 
refractive index. The situation for angle-resolved scattering quantities is expected to be 
different. Large relative errors observed in this paper are due to deep minima that are a 
consequence of both spherical symmetry and purely real refractive index. It is expected that 
visual agreement between the DDA results and the exact solution (as shown in Fig. 8 – Fig. 
10) should not change significantly for more general cases, however it will result in smaller 
relative errors, especially for larger x and smaller m. 
5 Conclusion 
In this paper we present the ADDA, a computer program to simulate light scattering by 
arbitrarily shaped particles. ADDA can parallelize a single DDA simulation, which allows it 
not to be limited by the memory of a single computer. Moreover, ADDA is heavily optimized, 
which allows it to compare favorably with other programs based on DDA when running on a 
single processor. We showed its capabilities for simulating light scattering by spheres with x 
up to 160 and m up to 2. The maximum reachable x on a cluster of 64 modern processors 
decrease rapidly with increasing m: it is 160 for 05.1=m  and only 20-40 (depending on the 
convergence threshold) for . This is mostly due to the slow convergence of the iterative 
solver leading to practically unbearable computation times. It is expected that larger particle 
sizes can be reached if m has a significant imaginary part. 
Errors of both integral and angle-resolved scattering quantities show no systematic 
dependence on x, but generally increase with m. Errors of Qext and >< θcos  range from less 
than 0.01 % to 6 %. Maximum - and RMS relative errors of S11(θ ) are in the ranges 0.2–18 
and 0.04–1 respectively. Error predictions of the traditional “rule of thumb” have very limited 
application in this range of x and m: it describes the upper limit of errors of Qext and >< θcos , 
however it does not account for the decrease of the errors with m. 
Currently, the ADDA is capable of simulating light scattering by almost all biological 
cells in suspension; however, its performance for other cases can be improved. These 
improvements, left for future work, may include improving the convergence of the iterative 
solver by preconditioning. It also is desirable to conduct a detailed study of the dependence of 
the accuracy of the final results on the size of the dipole and convergence thresholds of the 
iterative solver for different scatterers. Such a study should result in a reduction of the 
computation time and provide a realistic estimate of DDA complexity over a wide range of x 
and m. 
Acknowledgements 
We thank Gorden Videen for critically reading the manuscript and anonymous reviewer for 
valuable comments. Our research is supported by Siberian Branch of the Russian Academy of 
Sciences through the grant 2006-03. 
References 
 [1]  Purcell EM, Pennypacker CR. Scattering and adsorption of light by nonspherical dielectric grains. 
Astrophys J 1973;186:705-714. 
 [2]  Draine BT, Flatau PJ. Discrete-dipole approximation for scattering calculations. J Opt Soc Am A 
1994;11:1491-1499. 
 [3]  Yurkin MA, Hoekstra AG. The discrete dipole approximation: an overview and recent developments. J 
Quant Spectrosc Radiat Transf 2007, doi:10.1016/j.jqsrt.2007.01.034. 
 [4]  Penttila A, Zubko E, Lumme K, Muinonen K, Yurkin MA, Draine BT, Rahola J, Hoekstra AG, Shkuratov 
Y. Comparison between discrete dipole implementations and exact techniques. J Quant Spectrosc Radiat 
Transf 2007, doi:10.1016/j.jqsrt.2007.01.026. 
 [5]  Draine BT, Flatau PJ. User guide for the discrete dipole approximation code DDSCAT 6.1. 
http://xxx.arxiv.org/abs/astro-ph/0409262, 2004. 
 [6]  Hoekstra AG. Computer simulations of elastic light scattering. PhD thesis. University of Amsterdam, 
Amsterdam, 1994. 
 [7]  Hoekstra AG, Sloot PMA. Coupled dipole simulations of elastic light scattering on parallel systems. Int J 
Mod Phys C 1995;6:663-679. 
 [8]  Hoekstra AG, Grimminck MD, Sloot PMA. Large scale simulations of elastic light scattering by a fast 
discrete dipole approximation. Int J Mod Phys C 1998;9:87-102. 
 [9]  Yurkin MA, Semyanov KA, Tarasov PA, Chernyshev AV, Hoekstra AG, Maltsev VP. Experimental and 
theoretical study of light scattering by individual mature red blood cells with scanning flow cytometry and 
discrete dipole approximation. Appl Opt 2005;44:5249-5256. 
 [10]  Temperton C. Self-sorting mixed-radix fast Fourier transforms. J Comp Phys 1983;52:1-23. 
 [11]  Frigo M, Johnson SG. FFTW: an adaptive software architecture for the FFT. Proc ICASSP 1998;3:1381-
1384. 
 [12]  Draine BT. The discrete-dipole approximation and its application to interstellar graphite grains. Astrophys 
J 1988;333:848-872. 
 [13]  Draine BT, Goodman JJ. Beyond clausius-mossotti - wave-propagation on a polarizable point lattice and 
the discrete dipole approximation. Astrophys J 1993;405:685-697. 
 [14]  Gutkowicz-Krusin D, Draine BT. Propagation of electromagnetic waves on a rectangular lattice of 
polarizable points. http://xxx.arxiv.org/abs/astro-ph/0403082, 2004. 
 [15]  Barrett R, Berry M, Chan TF, Demmel J, Donato J, Dongarra J, Eijkhout V, Pozo R, Romine C, van der 
Vorst HA. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd ed. 
SIAM, 1994. 
 [16]  Freund RW. Conjugate gradient-type methods for linear-systems with complex symmetrical coefficient 
matrices. SIAM J Sci Stat Comp 1992;13:425-448. 
 [17]  Flatau PJ. Improvements in the discrete-dipole approximation method of computing scattering and 
absorption. Opt Lett 1997;22:1205-1207. 
 [18]  Davis PJ, Rabinowitz P. Methods of Numerical Integration. New York: Academic Press, 1975. 
 [19]  Hoekstra AG, Frijlink M, Waters LBFM, Sloot PMA. Radiation forces in the discrete-dipole 
approximation. J Opt Soc Am A 2001;18:1944-1953. 
 [20]  Rahola J. Solution of dense systems of linear equations in the discrete-dipole approximation. SIAM J Sci 
Comp 1996;17:78-89. 
http://xxx.arxiv.org/abs/astro-ph/0409262,
http://xxx.arxiv.org/abs/astro-ph/0403082,
 [21]  Fan ZH, Wang DX, Chen RS, Yung EKN. The application of iterative solvers in discrete dipole 
approximation method for computing electromagnetic scattering. Microwave Opt Tech Lett 
2006;48:1741-1746. 
 [22]  Bohren CF, Huffman DR. Absorption and scattering of Light by Small Particles. New York: Wiley, 1983. 
 [23]  Rahola J. Iterative solution of dense linear systems arising from integral equations. Appl Parall Comput , 
Lect Not Comp Sci 1998;1541:460-467. 
 [24]  Rahola J. On the eigenvalues of the volume integral operator of electromagnetic scattering. SIAM J Sci 
Comp 2000;21:1740-1754. 
 [25]  Budko NV, Samokhin AB. Spectrum of the volume integral operator of electromagnetic scattering. SIAM 
J Sci Comp 2006;28:682-700. 
 [26]  Hoekstra AG, Rahola J, Sloot PMA. Accuracy of internal fields in volume integral equation simulations 
of light scattering. Appl Opt 1998;37:8482-8497. 
	1  Introduction
	2 ADDA computer code
	3 Numerical simulations
	3.1 Simulation parameters
	3.2 Results
	4 Discussion
	5 Conclusion
	Acknowledgements
	References
  /ASCII85EncodePages false
  /AllowTransparency false
  /AutoPositionEPSFiles true
  /AutoRotatePages /None
  /Binding /Left
  /CalGrayProfile (Dot Gain 20%)
  /CalRGBProfile (sRGB IEC61966-2.1)
  /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2)
  /sRGBProfile (sRGB IEC61966-2.1)
  /CannotEmbedFontPolicy /Error
  /CompatibilityLevel 1.4
  /CompressObjects /Tags
  /CompressPages true
  /ConvertImagesToIndexed true
  /PassThroughJPEGImages true
  /CreateJDFFile false
  /CreateJobTicket false
  /DefaultRenderingIntent /Default
  /DetectBlends true
  /DetectCurves 0.0000
  /ColorConversionStrategy /CMYK
  /DoThumbnails false
  /EmbedAllFonts true
  /EmbedOpenType false
  /ParseICCProfilesInComments true
  /EmbedJobOptions true
  /DSCReportingLevel 0
  /EmitDSCWarnings false
  /EndPage -1
  /ImageMemory 1048576
  /LockDistillerParams false
  /MaxSubsetPct 100
  /Optimize true
  /OPM 1
  /ParseDSCComments true
  /ParseDSCCommentsForDocInfo true
  /PreserveCopyPage true
  /PreserveDICMYKValues true
  /PreserveEPSInfo true
  /PreserveFlatness true
  /PreserveHalftoneInfo false
  /PreserveOPIComments true
  /PreserveOverprintSettings true
  /StartPage 1
  /SubsetFonts true
  /TransferFunctionInfo /Apply
  /UCRandBGInfo /Preserve
  /UsePrologue false
  /ColorSettingsFile ()
  /AlwaysEmbed [ true
  /NeverEmbed [ true
  /AntiAliasColorImages false
  /CropColorImages true
  /ColorImageMinResolution 300
  /ColorImageMinResolutionPolicy /OK
  /DownsampleColorImages true
  /ColorImageDownsampleType /Bicubic
  /ColorImageResolution 300
  /ColorImageDepth -1
  /ColorImageMinDownsampleDepth 1
  /ColorImageDownsampleThreshold 1.50000
  /EncodeColorImages true
  /ColorImageFilter /DCTEncode
  /AutoFilterColorImages true
  /ColorImageAutoFilterStrategy /JPEG
  /ColorACSImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /ColorImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /JPEG2000ColorACSImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /JPEG2000ColorImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /AntiAliasGrayImages false
  /CropGrayImages true
  /GrayImageMinResolution 300
  /GrayImageMinResolutionPolicy /OK
  /DownsampleGrayImages true
  /GrayImageDownsampleType /Bicubic
  /GrayImageResolution 300
  /GrayImageDepth -1
  /GrayImageMinDownsampleDepth 2
  /GrayImageDownsampleThreshold 1.50000
  /EncodeGrayImages true
  /GrayImageFilter /DCTEncode
  /AutoFilterGrayImages true
  /GrayImageAutoFilterStrategy /JPEG
  /GrayACSImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /GrayImageDict <<
    /QFactor 0.15
    /HSamples [1 1 1 1] /VSamples [1 1 1 1]
  /JPEG2000GrayACSImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /JPEG2000GrayImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 30
  /AntiAliasMonoImages false
  /CropMonoImages true
  /MonoImageMinResolution 1200
  /MonoImageMinResolutionPolicy /OK
  /DownsampleMonoImages true
  /MonoImageDownsampleType /Bicubic
  /MonoImageResolution 1200
  /MonoImageDepth -1
  /MonoImageDownsampleThreshold 1.50000
  /EncodeMonoImages true
  /MonoImageFilter /CCITTFaxEncode
  /MonoImageDict <<
    /K -1
  /AllowPSXObjects false
  /CheckCompliance [
    /None
  /PDFX1aCheck false
  /PDFX3Check false
  /PDFXCompliantPDFOnly false
  /PDFXNoTrimBoxError true
  /PDFXTrimBoxToMediaBoxOffset [
    0.00000
    0.00000
    0.00000
    0.00000
  /PDFXSetBleedBoxToMediaBox true
  /PDFXBleedBoxToTrimBoxOffset [
    0.00000
    0.00000
    0.00000
    0.00000
  /PDFXOutputIntentProfile ()
  /PDFXOutputConditionIdentifier ()
  /PDFXOutputCondition ()
  /PDFXRegistryName ()
  /PDFXTrapped /False
  /Description <<
    /CHS <FEFF4f7f75288fd94e9b8bbe5b9a521b5efa7684002000410064006f006200650020005000440046002065876863900275284e8e9ad88d2891cf76845370524d53705237300260a853ef4ee54f7f75280020004100630072006f0062006100740020548c002000410064006f00620065002000520065006100640065007200200035002e003000204ee553ca66f49ad87248672c676562535f00521b5efa768400200050004400460020658768633002>
    /CHT <FEFF4f7f752890194e9b8a2d7f6e5efa7acb7684002000410064006f006200650020005000440046002065874ef69069752865bc9ad854c18cea76845370524d5370523786557406300260a853ef4ee54f7f75280020004100630072006f0062006100740020548c002000410064006f00620065002000520065006100640065007200200035002e003000204ee553ca66f49ad87248672c4f86958b555f5df25efa7acb76840020005000440046002065874ef63002>
    /DAN <FEFF004200720075006700200069006e0064007300740069006c006c0069006e006700650072006e0065002000740069006c0020006100740020006f007000720065007400740065002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e007400650072002c0020006400650072002000620065006400730074002000650067006e006500720020007300690067002000740069006c002000700072006500700072006500730073002d007500640073006b007200690076006e0069006e00670020006100660020006800f8006a0020006b00760061006c0069007400650074002e0020004400650020006f007000720065007400740065006400650020005000440046002d0064006f006b0075006d0065006e0074006500720020006b0061006e002000e50062006e00650073002000690020004100630072006f00620061007400200065006c006c006500720020004100630072006f006200610074002000520065006100640065007200200035002e00300020006f00670020006e0079006500720065002e>
    /DEU <FEFF00560065007200770065006e00640065006e0020005300690065002000640069006500730065002000450069006e007300740065006c006c0075006e00670065006e0020007a0075006d002000450072007300740065006c006c0065006e00200076006f006e002000410064006f006200650020005000440046002d0044006f006b0075006d0065006e00740065006e002c00200076006f006e002000640065006e0065006e002000530069006500200068006f006300680077006500720074006900670065002000500072006500700072006500730073002d0044007200750063006b0065002000650072007a0065007500670065006e0020006d00f60063006800740065006e002e002000450072007300740065006c006c007400650020005000440046002d0044006f006b0075006d0065006e007400650020006b00f6006e006e0065006e0020006d006900740020004100630072006f00620061007400200075006e0064002000410064006f00620065002000520065006100640065007200200035002e00300020006f0064006500720020006800f600680065007200200067006500f600660066006e00650074002000770065007200640065006e002e>
    /ESP <FEFF005500740069006c0069006300650020006500730074006100200063006f006e0066006900670075007200610063006900f3006e0020007000610072006100200063007200650061007200200064006f00630075006d0065006e0074006f00730020005000440046002000640065002000410064006f0062006500200061006400650063007500610064006f00730020007000610072006100200069006d0070007200650073006900f3006e0020007000720065002d0065006400690074006f007200690061006c00200064006500200061006c00740061002000630061006c0069006400610064002e002000530065002000700075006500640065006e00200061006200720069007200200064006f00630075006d0065006e0074006f00730020005000440046002000630072006500610064006f007300200063006f006e0020004100630072006f006200610074002c002000410064006f00620065002000520065006100640065007200200035002e003000200079002000760065007200730069006f006e0065007300200070006f00730074006500720069006f007200650073002e>
    /FRA <FEFF005500740069006c006900730065007a00200063006500730020006f007000740069006f006e00730020006100660069006e00200064006500200063007200e900650072002000640065007300200064006f00630075006d0065006e00740073002000410064006f00620065002000500044004600200070006f0075007200200075006e00650020007100750061006c0069007400e90020006400270069006d007000720065007300730069006f006e00200070007200e9007000720065007300730065002e0020004c0065007300200064006f00630075006d0065006e00740073002000500044004600200063007200e900e90073002000700065007500760065006e0074002000ea0074007200650020006f007500760065007200740073002000640061006e00730020004100630072006f006200610074002c002000610069006e00730069002000710075002700410064006f00620065002000520065006100640065007200200035002e0030002000650074002000760065007200730069006f006e007300200075006c007400e90072006900650075007200650073002e>
    /ITA <FEFF005500740069006c0069007a007a006100720065002000710075006500730074006500200069006d0070006f007300740061007a0069006f006e00690020007000650072002000630072006500610072006500200064006f00630075006d0065006e00740069002000410064006f00620065002000500044004600200070006900f900200061006400610074007400690020006100200075006e00610020007000720065007300740061006d0070006100200064006900200061006c007400610020007100750061006c0069007400e0002e0020004900200064006f00630075006d0065006e007400690020005000440046002000630072006500610074006900200070006f00730073006f006e006f0020006500730073006500720065002000610070006500720074006900200063006f006e0020004100630072006f00620061007400200065002000410064006f00620065002000520065006100640065007200200035002e003000200065002000760065007200730069006f006e006900200073007500630063006500730073006900760065002e>
    /JPN <FEFF9ad854c18cea306a30d730ea30d730ec30b951fa529b7528002000410064006f0062006500200050004400460020658766f8306e4f5c6210306b4f7f75283057307e305930023053306e8a2d5b9a30674f5c62103055308c305f0020005000440046002030d530a130a430eb306f3001004100630072006f0062006100740020304a30883073002000410064006f00620065002000520065006100640065007200200035002e003000204ee5964d3067958b304f30533068304c3067304d307e305930023053306e8a2d5b9a306b306f30d530a930f330c8306e57cb30818fbc307f304c5fc59808306730593002>
    /KOR <FEFFc7740020c124c815c7440020c0acc6a9d558c5ec0020ace0d488c9c80020c2dcd5d80020c778c1c4c5d00020ac00c7a50020c801d569d55c002000410064006f0062006500200050004400460020bb38c11cb97c0020c791c131d569b2c8b2e4002e0020c774b807ac8c0020c791c131b41c00200050004400460020bb38c11cb2940020004100630072006f0062006100740020bc0f002000410064006f00620065002000520065006100640065007200200035002e00300020c774c0c1c5d0c11c0020c5f40020c2180020c788c2b5b2c8b2e4002e>
    /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken die zijn geoptimaliseerd voor prepress-afdrukken van hoge kwaliteit. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.)
    /NOR <FEFF004200720075006b00200064006900730073006500200069006e006e007300740069006c006c0069006e00670065006e0065002000740069006c002000e50020006f0070007000720065007400740065002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e00740065007200200073006f006d00200065007200200062006500730074002000650067006e0065007400200066006f00720020006600f80072007400720079006b006b0073007500740073006b00720069006600740020006100760020006800f800790020006b00760061006c0069007400650074002e0020005000440046002d0064006f006b0075006d0065006e00740065006e00650020006b0061006e002000e50070006e00650073002000690020004100630072006f00620061007400200065006c006c00650072002000410064006f00620065002000520065006100640065007200200035002e003000200065006c006c00650072002000730065006e006500720065002e>
    /PTB <FEFF005500740069006c0069007a006500200065007300730061007300200063006f006e00660069006700750072006100e700f50065007300200064006500200066006f0072006d00610020006100200063007200690061007200200064006f00630075006d0065006e0074006f0073002000410064006f0062006500200050004400460020006d00610069007300200061006400650071007500610064006f00730020007000610072006100200070007200e9002d0069006d0070007200650073007300f50065007300200064006500200061006c007400610020007100750061006c00690064006100640065002e0020004f007300200064006f00630075006d0065006e0074006f00730020005000440046002000630072006900610064006f007300200070006f00640065006d0020007300650072002000610062006500720074006f007300200063006f006d0020006f0020004100630072006f006200610074002000650020006f002000410064006f00620065002000520065006100640065007200200035002e0030002000650020007600650072007300f50065007300200070006f00730074006500720069006f007200650073002e>
    /SUO <FEFF004b00e40079007400e40020006e00e40069007400e4002000610073006500740075006b007300690061002c0020006b0075006e0020006c0075006f00740020006c00e400680069006e006e00e4002000760061006100740069007600610061006e0020007000610069006e006100740075006b00730065006e002000760061006c006d0069007300740065006c00750074007900f6006800f6006e00200073006f00700069007600690061002000410064006f0062006500200050004400460020002d0064006f006b0075006d0065006e007400740065006a0061002e0020004c0075006f0064007500740020005000440046002d0064006f006b0075006d0065006e00740069007400200076006f0069006400610061006e0020006100760061007400610020004100630072006f0062006100740069006c006c00610020006a0061002000410064006f00620065002000520065006100640065007200200035002e0030003a006c006c00610020006a006100200075007500640065006d006d0069006c006c0061002e>
    /SVE <FEFF0041006e007600e4006e00640020006400650020006800e4007200200069006e0073007400e4006c006c006e0069006e006700610072006e00610020006f006d002000640075002000760069006c006c00200073006b006100700061002000410064006f006200650020005000440046002d0064006f006b0075006d0065006e007400200073006f006d002000e400720020006c00e4006d0070006c0069006700610020006600f60072002000700072006500700072006500730073002d007500740073006b00720069006600740020006d006500640020006800f600670020006b00760061006c0069007400650074002e002000200053006b006100700061006400650020005000440046002d0064006f006b0075006d0065006e00740020006b0061006e002000f600700070006e00610073002000690020004100630072006f0062006100740020006f00630068002000410064006f00620065002000520065006100640065007200200035002e00300020006f00630068002000730065006e006100720065002e>
    /ENU (Use these settings to create Adobe PDF documents best suited for high-quality prepress printing.  Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.)
  /Namespace [
    (Adobe)
    (Common)
    (1.0)
  /OtherNamespaces [
    <<
      /AsReaderSpreads false
      /CropImagesToFrames true
      /ErrorControl /WarnAndContinue
      /FlattenerIgnoreSpreadOverrides false
      /IncludeGuidesGrids false
      /IncludeNonPrinting false
      /IncludeSlug false
      /Namespace [
        (Adobe)
        (InDesign)
        (4.0)
      ]
      /OmitPlacedBitmaps false
      /OmitPlacedEPS false
      /OmitPlacedPDF false
      /SimulateOverprint /Legacy
    >>
    <<
      /AddBleedMarks false
      /AddColorBars false
      /AddCropMarks false
      /AddPageInfo false
      /AddRegMarks false
      /ConvertColors /ConvertToCMYK
      /DestinationProfileName ()
      /DestinationProfileSelector /DocumentCMYK
      /Downsample16BitImages true
      /FlattenerPreset <<
        /PresetSelector /MediumResolution
      >>
      /FormElements false
      /GenerateStructure false
      /IncludeBookmarks false
      /IncludeHyperlinks false
      /IncludeInteractive false
      /IncludeLayers false
      /IncludeProfiles false
      /MultimediaHandling /UseObjectSettings
      /Namespace [
        (Adobe)
        (CreativeSuite)
        (2.0)
      ]
      /PDFXOutputIntentProfileSelector /DocumentCMYK
      /PreserveEditing true
      /UntaggedCMYKHandling /LeaveUntagged
      /UntaggedRGBHandling /UseDocumentProfile
      /UseDocumentBleed false
    >>
>> setdistillerparams
  /HWResolution [2400 2400]
  /PageSize [612.000 792.000]
>> setpagedevice
