CPDDA: A Python Package for Discrete Dipole Approximation Accelerated by CuPy

Abstract

Discrete Dipole Approximation (DDA) is a rapidly developing numerical method in recent years. DDA has found wide application in many research fields including plasmonics and atmospheric optics. Currently, few DDA packages based on Python have been reported. In this work, a Python package called CPDDA is developed. It can be used to simulate the light-scattering and -absorption properties of arbitrarily shaped particles. CPDDA uses object-oriented programming, offers high flexibility and extensibility, and provides a comprehensive database of refractive indices. The package uses the biconjugate gradient method and fast Fourier transform for program acceleration and memory optimization, and it uses parallel computation with graphics processing units to enhance program performance. The accuracy and performance of CPDDA were demonstrated by comparison with Mie theory, the MATLAB package MPDDA, and the Python package pyGDM2. Finally, CPDDA was used to simulate the variations in light-absorption and -scattering properties of ZnO@Au core-shell nanorods based on the particle size. CPDDA is useful for calculating light-scattering and -absorption properties of small particles and selecting materials with excellent optical properties.

1. Introduction

The interaction of light with matter is a hot topic in physics. In recent years, the interaction of light with particles has emerged in applications in various fields such as biomedicine [1], atmospheric science [2], and sensing [3]. Although Mie theory can accurately and efficiently analyze such interactions, its applicability is limited to simple geometries like spheres. There are various numerical methods for calculating the light-scattering and -absorption properties of particles with arbitrary shapes. These include Finite Difference Time Domain (FDTD) [4], Finite Element Method (FEM) [5], Boundary Element Method (BEM) [6], the null-field method [7,8], and Discrete Dipole Approximation (DDA) [9]. Compared to FDTD and FEM, DDA is faster to compute. For the detailed characteristics of these numerical methods, see the published paper written by Amirjani and Sadrnezhaad [10].

DDA was originally proposed by Purcell and Pennypacker to calculate light scattering from grains in the interstellar medium [11]. Over the years, DDA has been applied to calculate the optical properties of arbitrarily shaped particles. In DDA, the particle is discretized into an array of N dipole points, and these polarizable points respond to the local electric field, thereby acquiring a dipole moment. For a finite array of point dipoles, the scattering problem can be solved exactly. Although DDA is simple to implement and feature-rich, it requires solving a 3N × 3N linear system of equations. For particles larger than the wavelength of the incident light, N becomes extremely large (close to a million), and the matrix becomes dense. Clearly, directly solving such a linear system is impractical. This problem can be solved using iterative methods and fast Fourier transform (FFT). In addition to reducing the computational complexity from n² to nlog(n), FFT also makes the interaction matrices sparse [12], which better satisfies the conditions for using iterative methods like Lanczos-Type Product Methods (LTPMs) [13].

Currently, several DDA open-source programs are available on the web. Firstly, DDSCAT [14], developed by Draine and Flatau in Fortran, has undergone significant growth and iteration over time. ADDA [15], written in C by Yurkin and Hoekstra, is a computational package designed for DDA calculations in cases involving high refractive indices. It supports OpenCL hybrid programming, enabling the use of graphics processing units (GPUs) to accelerate the computations. Donald et al. developed OpenDDA [12] using C, and it provides two parallel computation modes: OpenMP and MPI. Loke et al. wrote DDA-SI [16] using MATLAB, which can calculate the case involving surface interactions. However, FFT and iterative methods are not used in it to speed up and reduce program memory. Shabaninezhad et al. developed a MATLAB package called MPDDA [17], which enables parallel acceleration through GPU. Chaumet et al. implemented IF-DDA [18] using Fortran with a graphical user interface (GUI), allowing the computation of interactions between inhomogeneous, anisotropic objects within multilayer structures. Additionally, Chaumet compared the performance of various iterative methods in DDA [19] and tested two methods called Induced Dimension Reduction and Generalized Product-Bi-conjugate Gradient Stabilized, both of which had never been applied to DDA before.

It is a clear trend in computational sciences over the last two decades that researchers tend to move from low-level to high-level languages like MATLAB and Python. These languages allow for fast algorithm implementation and are easy to read, attracting more researchers and simplifying code iteration. High-performance computing libraries in C/C++, such as BLAS and LAPACK, can be directly called from Python. With just a few lines of code, Python can achieve computation speed comparable to multi-line compiled C/C++ code. For example, NumPy [20] is a library widely used for linear algebra and array operations. Additionally, GPUs can bypass the limitations of Python’s Global Interpreter Lock (GIL) to effectively accelerate the CPU-bound part of code. If DDA code can be easily implemented in Python without significant loss of computational efficiency, it could open up a new perspective for developing DDA packages in Python. It is worth noting that another volume integral method similar to DDA, known as the Green Dyadic Method (GMD) [21,22], has also been implemented in Python. However, unlike DDA, the GMD method relies on the environment for its self-term, while the self-term of DDA depends on the discrete object itself. Furthermore, DDA incurs less overhead when computing anisotropic media. It is worth mentioning that the GDM python package pyGDM2 [22] and the null-field method python package SMUTHI [23] both have GPU acceleration.

In this paper, we have developed a Python package called CPDDA, which is available for free download on GitLab (https://gitlab.com/dmcxu1/cpdda (accessed on 21 March 2025)). CPDDA is object-oriented and provides two modes of operation: CPU and GPU. The GPU mode offers significantly higher computational efficiency compared to the CPU mode. In addition, the package includes a comprehensive database of refractive indices (3107 experimental data for 398 materials), which can offer powerful assistance for researching and optimizing the optical properties of various particles. Finally, the light-absorption and -scattering properties of ZnO@Au core-shell nanorods were investigated using CPDDA as an example.

2. Fundamentals of DDA

DDA is a numerical method based on volume-integral discretization derived directly from Maxwell’s equations [14]. It calculates optical properties like absorption, scattering, and extinction cross sections for particles of any shape and composition. The method approximates a continuous target by a finite array of polarizable points. The target is discretized into N points, each represented as a dipole and periodically arranged on a cubic lattice with a side length of d. The points acquire dipole moments in response to the local electric field. Therefore, the only approximation in the DDA is replacing the continuum target with an array of N-point dipoles.

Assuming the center coordinates of each dipole are {r_p, p = 1, 2, …, N}, and the polarizability is {α_p, p = 1, 2, …, N}, each dipole generates a dipole moment under the influence of the local electric field:

(1)$P_p={\alpha }_p{E}_p,$

where E_p is the electric field at r_p due to the incident wave E_inc,p plus the contribution of each of the other N − 1 dipoles:

(2)$E_p=E_{\mathrm{inc},p}-\sum _{q\ne p}{A}_{pq}{P}_q,$

where A_pq is the 3 × 3 interaction matrix containing the Green tensor, which can be expressed as:

(3)$A_{pq}={\displaystyle \frac{\exp (ik{r}_{pq})}{r_{pq}}}\times \left[\left.k^2\left({\widehat{r}}_{pq}{\widehat{r}}_{pq}-I_3\right)+{\displaystyle \frac{ik{r}_{pq}-1}{r_{pq}^2}}\left(3{\widehat{r}}_{pq}{\widehat{r}}_{pq}-I_3\right)\right]\right.,\hspace{1em}\hspace{0.33em}p\ne q,$

where $k=\omega c$ is the wave number, $r_{pq}=\hspace{0.33em}|r_p-r_q|$ is the distance from dipole q to p, ${\widehat{r}}_{pq}=(r_p-r_q)/r_{pq}$ is a unit vector in the direction from dipole q to p, and I₃ represents the 3 × 3 identity matrix. The diagonal elements of the interaction matrix are defined as $A_{pp}={{\alpha }_p}^{-1}{I}_3$, which leads to the simplified form of Equation (2) as:

(4)${\displaystyle \sum _{q=1}^N{A}_{pq}}{P}_{pq}=E_{\mathrm{inc},p}.$

Equation (4) is the final linear equation system to be solved and the primary objective of DDA. Once Equation (4) has been solved for the unknown polarization P_p, the extinction, absorption, and scattering cross sections can be evaluated:

(5)$C_{\mathrm{ext}}={\displaystyle \frac{4\pi k}{E_0^2}}\sum _{p=1}^N\mathrm{Im}(E_{\mathrm{inc},p}^{\ast }\cdot P_p),$

(6)$C_{\mathrm{abs}}={\displaystyle \frac{4\pi k}{E_0^2}}\sum _{p=1}^N\left\{\mathrm{Im}[P_p\cdot {({\alpha }_p^{-1})}^{\ast }{P}_p^{\ast }]-{\displaystyle \frac{2}{3}}{k}^3|P_p{|}^2\right\},$

(7)$C_{\mathrm{sca}}=C_{\mathrm{ext}}-C_{\mathrm{abs}}.$

After discretizing the target object into dipoles, each dipole can be considered to represent a specific polarizability of the target’s material and can interact with other N − 1 dipoles. The dielectric properties of a material are directly linked to the polarizability. When dipoles are arranged on an infinite cubic lattice, this relationship is described by the simple and precise Clausius–Mossotti relation:

(8)${{\alpha }_p}^{CM}={\displaystyle \frac{3d^3}{4\pi }}{\displaystyle \frac{{ϵ}_p-1}{{ϵ}_p+2}},$

where ϵ_p is the dielectric function of the target’s material at location r_p. However, the derivation of this correction was based on the assumption that the electric field remains uniform across cubic regions with a volume of d³. Thus, we consider the Lattice Dispersion Relation (LDR) polarizability. Since the dipoles are arranged on a periodic cubic lattice and coupling occurs between the electromagnetic waves and the lattice, consider that for the LDR polarizability α, the infinite lattice of polarizable points has the same dispersion relation as the continuous refractive index m. In the limit kd << 1, the polarizability α is expressed as a power series expansion in terms of kd and m² = $ϵ$:

(9)${\alpha }^{\mathrm{LDR}}\approx {\displaystyle \frac{{\alpha }^{\mathrm{CM}}}{1+({\alpha }^{\mathrm{CM}}/d^3)[(b_1+b_2{m}^2+b_3{m}^{2}S){(kd)}^2-(2/3)i{(kd)}^3]}},$

where:

(10)$$\begin{gathered} b_1=-1.891531,\hspace{1em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}{b}_2=0.1648469, \\ {b}_3=-1.7700004,\hspace{1em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}S={\displaystyle \sum _{j=1}^3{({\widehat{a}}_p{\widehat{e}}_p)}^2.} \end{gathered}$$

In order to ensure the accuracy of the DDA calculation results, the rule of the following equation is followed:

(11)$kd\left|m\right|<<1.$

For irregular particles, the criteria may be appropriately adjusted to account for their specific characteristics [24]. The subsequent section describes the implementation of DDA in Python, including acceleration techniques, and provides a quantitative assessment of the code’s accuracy, reliability, and performance.

3. Implementation and Validation of CPDDA

3.1. Implementation of CPDDA

Our primary algorithms are based on MPDDA, and we have extended this framework to develop CPDDA using Python. CPDDA uses an object-oriented form to program the DDA. The structure of the CPDDA package and the main steps to set up and run a simulation are schematically depicted in Figure 1. The core of CPDDA is the simulation object, which contains information about the structure and the incident field used in the simulation. In Appendix A, a simple Python script for CPDDA to create a simulation object is supplied. The geometry of the particle and the refractive indices of the particle and the environment are stored in an instance of structures.struct. The geometry module is used to establish the grid coordinates of the particle, and the materials object contains the dipole material properties, including the refractive index of the particles and the corresponding dielectric constant. It is worth noting that CPDDA contains a rich database of refractive indices of materials, which provides 398 kinds of materials (including 3107 kinds of experimental data). Table 1 presents a summary of the database. The diverse range of refractive indices in this package’s database enables researchers to explore novel optical properties of materials and conduct comparative studies within material classes to identify high-performance candidates for specific applications. After setting up the first two objects, the polarizability of each dipole is determined by placeholders.

The second key ingredient of CPDDA simulation is the incident field. The part about the incident field is set in the field.efield object, which contains the propagation direction and polarization of the incident field. The current CPDDA implementation contains two incident fields: plane wave and Gaussian beam.

The above two parts together set up and define the simulation object, thus setting E_inc,p and the diagonal elements in A_pq in Equation (4). The main simulation of CPDDA is divided into constructing matrix–vector multiplications and solving them. For Equation (4), it can be written as follows:

(12)$\left[\begin{array}{ccc}{A}_{xx}& A_{xy}& A_{xz}\\ A_{yx}& A_{yy}& A_{yz}\\ A_{zx}& A_{zy}& A_{zz}\end{array}\right]\left[\begin{array}{c}{P}_x\\ P_y\\ P_z\end{array}\right]=\left[\begin{array}{c}{E}_{\mathrm{inc},x}\\ E_{\mathrm{inc},y}\\ E_{\mathrm{inc},z}\end{array}\right],$

where $A_{\mu \nu }$ ($\mu$ and $\nu$ is x, y, or z) is an N×N matrix that represents the interaction matrix, and $P_{\mu }$ and $E_{\mathrm{inc},\mu }$ ($\mu$ and $\nu$ is x, y, or z) are denoted as the three polarization vectors of N×1 and the three components of the incident electric field, respectively. Since the nine blocks of the interaction matrix are symmetric regarding changing the indices, only six independent blocks of A (A_xx, A_xy, A_xz, A_yy, A_yx, A_zx) are required.

After constructing the matrix, it is time to solve for this matrix–vector multiplication. Due to the size of the matrix, it takes a significant amount of computing resources to generate the interaction matrix’s inverse, and solving Equation (12) in this case is not possible. Therefore, the solution is to use an iterative method to solve the linear system. The principle of an iterative method for solving the linear system is to create a sequence E_s such that:

(13)$r_q={\displaystyle \frac{\parallel \overline{A}{E}_q-E_{\mathrm{inc},p}\parallel }{\parallel E_{\mathrm{inc},p}\parallel }},$

with the residue r_q tending towards zero when q increases, in which the q-th approximation is derived from the previous ones. The iterative process is stopped when r_q < η, where the value of η is set by the user and depends on the desired precision of the field [13,19]. Bi-conjugate gradient (Bi-CG) has been employed in several DDA implementations [12,17], and CPDDA also adopts this iterative method.

The iterative method involves several matrix–vector multiplications, which are further accelerated by 3D FFT. The interaction matrix is constructed as a Toeplitz matrix, where all diagonal elements are equal, and the off-diagonal elements along each parallel diagonal are equal. Thus, the Toeplitz matrix can be entirely defined by its first row and first column as [25]:

(14)$T_n=[t_{q,p};q,p=0, 1, \cdots ,\hspace{0.33em}n-1].$

By setting t_q,p = t_q-p and appending the first column (excluding the initial term) in reverse order to the first row, then completing the diagonal, the Toeplitz matrix is transformed into the circulant matrix $C^{n^{\prime }}$, where $n^{\prime }=2n-1$. Arbitrary circulant matrices can be diagonalized using Fourier matrices:

(15)$F^{n^{\prime }}={\displaystyle \frac{1}{\sqrt{n^{\prime }}}}{\omega }^{pq},\hspace{1em}p,\hspace{0.33em}q=0,\hspace{0.33em}\cdots ,\hspace{0.33em}{n}^{\prime }-1,\hspace{1em}\omega =e^{-{\displaystyle \frac{2\pi i}{n^{\prime }}}},$

(16)$F^{n^{\prime }}{C}^{n^{\prime }}{F}^{\mathrm{H}}=diag(\sqrt{n^{\prime }}{F}^{n^{\prime }}{C}^{n^{\prime }})=diag({\lambda }_c).$

Equation (15) is the expression for the Fourier matrix using the fact that the Fourier matrix is unitary, ${}^{[n^{\prime }]}F^{-1}={}^{[n^{\prime }]}F^{\mathrm{H}}$. Where the H superscript is the transposed conjugate, the multiplication of an arbitrary vector with a circulant matrix can be performed as an element-wise multiplication in the Fourier domain via the Convolution Theorem. For a comprehensive discussion of specific algorithms for iterative methods and FFT, readers may refer to the studies by Shabaninezhad et al. and Donald et al. [12,17]. All of the above operations for matrix–vector multiplication can be found in the core.DDA module. core.DDA is divided into CPU and GPU compute modes. The next section details how CPU-bound prat blocks are accelerated in GPU mode.

After calculating the polarization vector P_p, CPDDA can compute the near-field and far-field parameters and visualize them through the visu module. Specifically, the three main cross-sections, efficiency and volume coefficients, and the electric field distribution around the particles can be calculated; the corresponding code is available in the post-processing module.

3.2. GPU Mode of CPDDA

For CPU-bound tasks like DDA, parallel computing is an efficient computational method. Due to the GIL, explicit multi-threading cannot be implemented in Python to achieve parallelism. Using multi-threading through tools like Cython limits programming to basic atomic operations, which prevents using the high-performance capabilities of libraries like NumPy and increases programming complexity [26,27]. Therefore, GPU mode is added to core.DDA for parallel computation.

In general, GPUs have thousands of cores and are capable of running several thousands of threads simultaneously. Although modern CPUs can have an even larger number of cores, they are still very small compared to GPUs. Since the GPU is a foreign device compared to the host CPU, any data that need to be processed on a GPU must be offloaded to it first. The data transfer between CPU and GPU is performed via a PCIe cable. PCIe data communication can be bottlenecked by large data transfers, and in extreme cases, data transfer costs may exceed computation times. Therefore, considering the hardware communication overhead, only the CPU-bound part of the code is offloaded to the GPU, as shown in Figure 2. CPU-bound iterative methods and matrix–vector multiplications were offloaded to the GPU, with only two communications between the CPU and GPU during the entire computation. In the Python implementation, we used CuPy and NumPy for high-performance computing and transfered data between the CPU and GPU using the cp.array and cp.asnumpy. Additionally, the memory pooling technique in CuPy manually replaces some of Python’s garbage collection techniques for more efficient memory management. It is worth mentioning that the CuPy library offers a highly compatible interface with NumPy, enabling GPU-accelerated data processing and scientific computing, particularly for large-scale data. GPU acceleration is available when the user selects GPU mode in the core.DDA module.

3.3. Numerical Validation of CPDDA

To verify the accuracy of CPDDA, the calculations were performed using a LAPTOP-3EE0017T with Windows 10 operating system. The laptop was equipped with an Intel(R) i5-7300HQ processor with 2.50 GHz and 8192 MB of total RAM as well as a NVIDIA GeForce RTX 1050ti graphics card with 4096 MB of RAM, and Python 3.10 with 13.1.0 cupy-cuda12x was installed on the laptop. We calculated the extinction efficiency of the Au nanosphere with a diameter of 30 nm in an aqueous solution for two lattice spacings, d = 0.5 nm and 0.25 nm, and compared the results with the Mie theory using PyMieLab_V1.0 [28], as shown in Figure 3a. The CPDDA results closely matched Mie theory with a relative error of less than 2% for d = 0.25 nm, and the spectra converged to the exact solution as the lattice spacing was reduced. Additionally, scattering efficiency of ZnO@Au nanorods, characterized by a core length of 60 nm, a shell thickness of 5 nm, and an aspect ratio of 12, were calculated and compared with the results from MPDDA [17], as shown in Figure 3b. The results are in excellent agreement, demonstrating the accuracy and reliability of CPDDA.

To demonstrate the GPU acceleration of CPDDA, we compare the computation time for the CPU and GPU under the same configurations. Figure 4a shows the variations between the CPU and GPU computation times with lattice spacing d for an Au nanocube with an isomeric diameter of 40 nm in a medium with a refractive index 1.33 and incident light wavelength of 450 nm. As d decreases, the GPU acceleration becomes more pronounced. As the lattice pitch decreases, CPU computation time increases significantly, while GPU computation time remains relatively constant. At d = 0.3 nm, CPU time is approximately 26 times that of the GPU. To further evaluate performance, we compared the GPU mode of CPDDA with that of MPDDA on the same configured system, both using Bi-CG as the iterative solver. The computation time ratio between MPDDA and CPDDA demonstrates the acceleration advantage of CPDDA, as shown in Figure 4b. The results indicate that CPDDA is generally faster, but MPDDA outperforms CPDDA when d reaches 1.9 nm. CPDDA shows a computational speed advantage at smaller lattice spacings. To further evaluate the computational efficiency, we compared the performance of pyGDM2 and CPDDA in GPU mode. As Figure 4c shows, we simulated Au nanospheres (50 nm diameter in a vacuum) using both pyGDM2 and CPDDA under GPU acceleration based on the example provided in the pyGDM2 documentation. Due to the high memory demand of pyGDM2 compared to DDA for the same number of discrete cells, the analysis was limited to approximately 1000 discrete units. The results demonstrate that while CPDDA exhibits superior performance at shorter wavelengths, pyGDM2 outperforms at longer wavelengths.

4. Optical Properties of ZnO@Au Nanorods

ZnO@Au core-shell nanoparticles, consisting of ZnO core and Au shell, show significant potential in applications such as surface-enhanced Raman scattering, photocatalysis, and biosensing [29,30,31,32]. ZnO, a semiconductor with a bandgap of 3.37 eV at room temperature, exhibits good electrical conductivity, making it valuable in electrochemical devices [33]. Moreover, these nanoparticles are low in toxicity, biodegradable, and easy to synthesize [34]. When coated with Au, ZnO@Au nanoparticles exhibit enhanced biocompatibility. ZnO@Au core-shell nanorods have been synthesized and widely used in gas sensing [35], and other fields [36]. In this section, ZnO@Au core-shell nanorods are used to quantitatively analyze the variation in their light-absorption and -scattering properties with particle size using CPDDA. The refractive index data of Au and ZnO were obtained from ref. [37] and ref. [38], respectively, and the refractive index of the surrounding medium was set to 1.44 (representing subcutaneous fat) [39]. The light-absorption and -scattering characteristics were quantified by the volume-absorption coefficient A_abs and volume-scattering coefficient A_sca, defined as follows:

(17)$A_{\mathrm{abs}}=C_{\mathrm{abs}}/V_{\mathrm{p}},$

(18)$A_{\mathrm{sca}}=C_{\mathrm{sca}}/V_{\mathrm{p}},$

where V_p represents the overall volume of the particle (including the shell).

The schematic diagram of interaction of light with the ZnO@Au core-shell nanorod is shown in Figure 5. The ZnO@Au core-shell nanorod is represented by three geometrical parameters, where l_c denotes the core length, d_c denotes the core diameter, and t_s denotes the shell thickness. The ratio of core length to diameter is defined as the aspect ratio AR, and the particle geometry can also be expressed in terms of AR, t_s, and l_c. The incident light propagates along the short axis of the nanorod, and the electric field is polarized along the long axis of the nanorod.

4.1. Effect of Aspect Ratio

To quantitatively analyze the effect of the core aspect ratio AR on the absorption and scattering spectra, we calculated the spectra for ZnO@Au core-shell nanorods with AR values of 3, 6, 9, and 12, fixing the core length l_c at 80 nm and the shell thickness t_s at 5 nm. The results are shown in Figure 6a,c. As the AR increased from 3 to 12, the resonance wavelength redshifted from 965 nm to 1055 nm (see Figure 6b,d), the volume-absorption coefficient increased from 0.5256 nm⁻¹ to 1.6854 nm⁻¹ (see Figure 6b), and the volume-scattering coefficient first increased from 0.2024 nm⁻¹ to 0.2340 nm⁻¹, then decreased to 0.2160 nm⁻¹ (see Figure 6d). The imaginary part of the dielectric function is expected to exhibit an increasing trend with wavelength, as the collision frequency of electrons in the metal becomes comparable to the optical frequency in the lower energy (longer wavelength) regime. At these wavelengths, the absorption contribution is anticipated to rise due to the enhancement of the metal’s intrinsic dielectric function. This effect is particularly pronounced for Au, where the increase in the imaginary dielectric function is more significant. This could be the reason for the decrease in the maximum value of the bulk-scattering coefficient at longer wavelengths for high-aspect-ratio nanorods [40]. The optical properties and electronic structure of the nanoparticles changed with AR. As the AR increased, the charge accumulation on the nanoparticle surface decreased, weakening the electron-restoring force and lowering the Local Surface Plasmon Resonance (LSPR) frequency, which resulted in a redshift of the resonance wavelength [41].

4.2. Effect of Shell Thickness

To further analyze the effect of the shell thickness t_s on the absorption and scattering spectra, we calculated the spectra for ZnO@Au core-shell nanorods with t_s values of 2.5, 5, 7.5, and 10 nm, fixing the core length l_c at 80 nm and the core aspect ratio AR at 12. The results are shown in Figure 7a,c. As the t_s increased from 2.5 to 10 nm, The resonance wavelength of the ZnO@Au core-shell nanorods underwent a significant blueshift (see Figure 7b,d), the volume-absorption coefficient decreased from 1.8834 nm⁻¹ to 0.8885 nm⁻¹ (see Figure 7b), and the volume-scattering coefficient increased from 0.0521 nm⁻¹ to 0.5204 nm⁻¹ (see Figure 7d). As t_s increased, the accumulated surface charge on the nanoparticle enhanced the electron-restoring force, raising the LSPR frequency and inducing a blueshift in the resonance wavelength.

4.3. Effect of Core Length

To further analyze the effect of the core length l_c on the absorption and scattering spectra, we calculated the spectra for ZnO@Au core-shell nanorods with l_c values of 40, 60, 80, and 100 nm, fixing the shell thickness t_s at 5 nm and the core aspect ratio AR at 12. The results are shown in Figure 8a,c. As the l_c increased from 40 to 100 nm, the resonance wavelength redshifted from 875 nm to 1300 nm (see Figure 8b,d), the volume-absorption coefficient decreased from 1.8589 nm⁻¹ to 1.6166 nm⁻¹ (see Figure 8b), and the volume-scattering coefficient increased from 0.1444 nm⁻¹ to 0.2440 nm⁻¹ (see Figure 8d). As l_c increased, the vibrational distance of free electrons within the Au shell increased, decreasing their vibrational frequency and lowering the LSPR frequency, which induced a redshift in the resonance wavelength.

5. Conclusions

This paper presents CPDDA, a Python package for simulating the light-scattering and -absorption properties of arbitrarily shaped particles. CPDDA uses Bi-CG and FFT algorithms to accelerate the CPU-bound part of DDA and reduce memory usage. To improve performance and overcome GIL limitations, CuPy was used to implement parallel computation on the GPU, significantly enhancing computation speed compared to the CPU. Having a rich database of refractive indices of materials (3107 refractive index experimental data of 398 materials) is one of the features of this package, which can provide powerful help for studying and optimizing electromagnetic properties of particles with multiple materials. Compared to Mie theory, MPDDA, and pyGDM2, the reliability and performance of CPDDA were validated. Finally, ZnO@Au core-shell nanorods were used to study how their light-absorption and -scattering properties varied with aspect ratio, shell thickness, and core length. The high-level Python API, along with various tools for rapid data analysis and visualization, renders standard optical simulations of particles relatively straightforward. In the future, more optical parameters will be added to CPDDA to expand the functionality of the package, and more high-performance algorithms will be integrated to support broader applications.

Author Contributions

Conceptualization: D.X., P.T. and S.L.; methodology: D.X. and P.T.; software: D.X., P.T. and S.L.; validation: D.X., P.T. and L.J.; formal analysis: D.X., P.T. and S.L.; investigation: D.X., P.T., S.L., M.W. and L.J.; resources: P.T. and S.L.; data curation: D.X. and M.W.; writing—original draft preparation: D.X., P.T. and S.L.; writing—review and editing: D.X., P.T. and S.L.; visualization: D.X., M.W. and L.J.; supervision: P.T. and S.L.; project administration: P.T.; funding acquisition: P.T. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in this study are included in the article; further inquiries can be directed to the corresponding authors.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Structure of the CPDDA package and workflow of the typical simulation: (1) Setup of the geometry, material, and incident electric field. They are bundled in an instance of the simulation object. (2) Main simulation. (3) Post-processing. (4) Visualization of the results.

Figure 2. Schematic diagram of CPDDA data transfer from CPU to GPU.

Figure 3. (a) Calculated extinction spectra of Au nanosphere using Mie theory and CPDDA (lattice spacing is 0.5 nm and 0.25 nm, respectively). (b) Calculated scattering spectra of ZnO@Au nanorod using MPDDA and CPDDA. In the simulation, the diameter of Au nanosphere is 30 nm, and core length, shell thickness, and aspect ratio of core-shell nanorod are 60 nm, 5 nm, and 12, respectively. The refractive indices of the environment are (a) 1.33 and (b) 1.44, respectively.

Figure 4. (a) Variation of CPU and GPU computation time and their ratio with the lattice spacing d. (b) Variation of the ratio of MPDDA to CPDDA computation time with the lattice spacing d. In the simulation, the particle is an Au nanocube with an isomeric diameter of 40 nm, the refractive index of the environment is 1.33, and the incident light is a plane wave with a wavelength of 450 nm. (c) Variation of the ratio of pyGDM2 to CPDDA computation time with the wavelengths λ, the diameter of Au nanosphere is 50 nm, the number of discrete units for both is around 1000.

Figure 5. Schematic of the interaction of ZnO@Au core-shell nanorod with light.

Figure 6. Calculated (a) absorption and (c) scattering spectra of ZnO@Au core-shell nanorods with different aspect ratios AR. Variation of (b) absorption and (d) scattering resonance wavelength and peak intensity with AR. The core length lc and shell thickness ts were 80 nm and 5 nm, respectively.

Figure 7. Calculated (a) absorption and (c) scattering spectra of ZnO@Au core-shell nanorods with different shell thicknesses ts. Variation of (b) absorption and (d) scattering resonance wavelength and peak intensity with ts. The core length lc and aspect ratio AR are 80 nm and 12, respectively.

Figure 8. Calculated (a) absorption and (c) scattering spectra of ZnO@Au core-shell nanorods with different core lengths lc. Variation of (b) absorption and (d) scattering resonance wavelength and peak intensity with lc. The shell thickness ts and aspect ratio AR were 5 nm and 12, respectively.

Table 1

Description of material refractive index database.

Material Classification	Description
MAIN	simple inorganic materials
GLASS	glasses
ORGANIC	organic materials
OTHER	miscellaneous materials

Appendix A. A CPDDA Code Example for Setting the Simulation Object

import numpy as npfrom CPDDA import materialsfrom CPDDA import structuresFrom CPDDA import fieldsfrom CPDDA import corefrom CPDDA import post_processingd = 2r_eff = 20wavelength = np.arange(450, 570, 5, dtype=np.float32)material1=materials.FromDatabase(shelf_name="main",book_name="Au",page_name="Johnson.yml",wl=wavelength, nb=1.33)material = [material1]geometry = structures.sphere(r_eff, d) occupied = structures.INDEX_in("spherical", geometry)struct = structures.struct(d, material, occupied, geometry)print(struct)field_generator = fields.plan_wave()K_kwargs = dict(K0=np.array([1, 0, 0]))E_kwargs = dict(E0=np.array([0, 0, 1]))efield = fields.efield(field_generator, wavelength, K_kwargs, E_kwargs)print(efield)sim = core.simulation(struct, efield)

References

1. Kumar, A.; Kim, S.; Nam, J.-M. Plasmonically Engineered Nanoprobes for Biomedical Applications. J. Am. Chem. Soc.; 2016; 138, pp. 14509-14525. [DOI: https://dx.doi.org/10.1021/jacs.6b09451] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27723324]

2. Li, D.; Yue, W.; Gao, P.; Gong, T.; Wang, C.; Luo, X. Surface-enhanced Raman spectroscopy (SERS) for the characterization of atmospheric aerosols: Current status and challenges. TrAC Trends Anal. Chem.; 2024; 170, 117426. [DOI: https://dx.doi.org/10.1016/j.trac.2023.117426]

3. Khurana, K.; Jaggi, N. Localized Surface Plasmonic Properties of Au and Ag Nanoparticles for Sensors: A Review. Plasmonics; 2021; 16, pp. 981-999. [DOI: https://dx.doi.org/10.1007/s11468-021-01381-1]

4. Taflove, A.; Hagness, S.C. Computational Electrodynamics: The Finite-Difference Time-Domain Method; 3rd ed. Artech House: Boston, MA, USA, 2005.

5. Jin, J.-M. The Finite Element Method in Electromagnetics; 3rd ed. IEEE Press: Piscataway, NJ, USA, 2014.

6. Liu, Y.J.; Mukherjee, S.; Nishimura, N.; Schanz, M.; Ye, W.; Sutradhar, A.; Pan, E.; Dumont, N.A.; Frangi, A.; Saez, A. Recent Advances and Emerging Applications of the Boundary Element Method. Appl. Mech. Rev.; 2011; 64, 30802. [DOI: https://dx.doi.org/10.1115/1.4005491]

7. Doicu, A.; Mishchenko, M.I. An overview of the null-field method. I: Formulation and basic results. Phys. Open; 2020; 5, 100020. [DOI: https://dx.doi.org/10.1016/j.physo.2020.100020]

8. Doicu, A.; Mishchenko, M.I. An overview of the null-field method. II: Convergence and numerical stability. Phys. Open; 2020; 3, 100019. [DOI: https://dx.doi.org/10.1016/j.physo.2020.100019]

9. Chaumet, P.C. The Discrete Dipole Approximation: A Review. Mathematics; 2022; 10, 3049. [DOI: https://dx.doi.org/10.3390/math10173049]

10. Amirjani, A.; Sadrnezhaad, S.K. Computational electromagnetics in plasmonic nanostructures. J. Mater. Chem. C; 2021; 9, pp. 9791-9819. [DOI: https://dx.doi.org/10.1039/d1tc01742j]

11. Purcell, E.M.; Pennypacker, C.R. Scattering and Absorption of Light by Nonspherical Dielectric Grains. Astrophys. J.; 1973; 186, pp. 705-714. [DOI: https://dx.doi.org/10.1086/152538]

12. Donald, J.M.; Golden, A.; Jennings, S.G. Opendda: A Novel High-Performance Computational Framework for the Discrete Dipole Approximation. Int. J. High Perform. Comput. Appl.; 2009; 23, pp. 42-61. [DOI: https://dx.doi.org/10.1177/1094342008097914]

13. Horiuchi, I.; Aihara, K.; Suzuki, T.; Ishiwata, E. Global GPBiCGstab(L) method for solving linear matrix equations. Numer. Algorithms; 2023; 93, pp. 295-319. [DOI: https://dx.doi.org/10.1007/s11075-022-01415-7]

14. Draine, B.T.; Flatau, P.J. Discrete-Dipole Approximation For Scattering Calculations. J. Opt. Soc. Am. A; 1994; 11, pp. 1491-1499. [DOI: https://dx.doi.org/10.1364/JOSAA.11.001491]

15. Yurkin, M.A.; Hoekstra, A.G. The discrete-dipole-approximation code ADDA: Capabilities and known limitations. J. Quant. Spectrosc. Radiat. Transf.; 2011; 112, pp. 2234-2247. [DOI: https://dx.doi.org/10.1016/j.jqsrt.2011.01.031]

16. Loke, V.L.Y.; Pinar Mengüç, M.; Nieminen, T.A. Discrete-dipole approximation with surface interaction: Computational toolbox for MATLAB. J. Quant. Spectrosc. Radiat. Transf.; 2011; 112, pp. 1711-1725. [DOI: https://dx.doi.org/10.1016/j.jqsrt.2011.03.012]

17. Shabaninezhad, M.; Awan, M.G.; Ramakrishna, G. MATLAB package for discrete dipole approximation by graphics processing unit: Fast Fourier Transform and Biconjugate Gradient. J. Quant. Spectrosc. Radiat. Transf.; 2021; 262, 107501. [DOI: https://dx.doi.org/10.1016/j.jqsrt.2020.107501]

18. Chaumet, P.C.; Sentenac, D.; Maire, G.; Rasedujjaman, M.; Zhang, T.; Sentenac, A. IFDDA, an easy-to-use code for simulating the field scattered by 3D inhomogeneous objects in a stratified medium: Tutorial. J. Opt. Soc. Am. A; 2021; 38, pp. 1841-1852. [DOI: https://dx.doi.org/10.1364/JOSAA.432685]

19. Chaumet, P.C. A comparative study of efficient iterative solvers for the discrete dipole approximation. J. Quant. Spectrosc. Radiat. Transf.; 2023; 312, 108816. [DOI: https://dx.doi.org/10.1016/j.jqsrt.2023.108816]

20. Van Der Walt, S.; Colbert, S.C.; Varoquaux, G. The NumPy Array: A Structure for Efficient Numerical Computation. Comput. Sci. Eng.; 2011; 13, pp. 22-30. [DOI: https://dx.doi.org/10.1109/MCSE.2011.37]

21. Wiecha, P.R. pyGDM—A python toolkit for full-field electro-dynamical simulations and evolutionary optimization of nanostructures. Comput. Phys. Commun.; 2018; 233, pp. 167-192. [DOI: https://dx.doi.org/10.1016/j.cpc.2018.06.017]

22. Wiecha, P.R.; Majorel, C.; Arbouet, A.; Patoux, A.; Brûlé, Y.; Francs, G.C.D.; Girard, C. “pyGDM”—New functionalities and major improvements to the python toolkit for nano-optics full-field simulations. Comput. Phys. Commun.; 2022; 270, 108142. [DOI: https://dx.doi.org/10.1016/j.cpc.2021.108142]

23. Egel, A.; Czajkowski, K.M.; Theobald, D.; Ladutenko, K.; Kuznetsov, A.S.; Pattelli, L. SMUTHI: A python package for the simulation of light scattering by multiple particles near or between planar interfaces. J. Quant. Spectrosc. Radiat. Transf.; 2021; 273, 107846. [DOI: https://dx.doi.org/10.1016/j.jqsrt.2021.107846]

24. Zubko, E.; Kochergin, A.; Videen, G. Convergence of the DDA for ensembles of objects of irregular shape. J. Quant. Spectrosc. Radiat. Transf.; 2024; 314, 108854. [DOI: https://dx.doi.org/10.1016/j.jqsrt.2023.108854]

25. Zhu, Z.; Wakin, M.B. On the Asymptotic Equivalence of Circulant and Toeplitz Matrices. IEEE Trans. Inf. Theory; 2017; 63, pp. 2975-2992. [DOI: https://dx.doi.org/10.1109/TIT.2017.2676808]

26. Behnel, S.; Bradshaw, R.; Citro, C.; Dalcin, L.; Sverre Seljebotn, D.; Smith, K. Cython: The best of both worlds. Comput. Sci. Eng.; 2011; 13, pp. 31-39. [DOI: https://dx.doi.org/10.1109/MCSE.2010.118]

27. Smith, K.W. Cython: A Guide for Python Programmers; 1st ed. O’Reilly Media: Tokyo, Japan, 2015.

28. Ma, D.; Tuersun, P.; Cheng, L.; Zheng, Y.; Abulaiti, R. PyMieLab_V1.0: A software for calculating the light scattering and absorption of spherical particles. Heliyon; 2022; 8, e11469. [DOI: https://dx.doi.org/10.1016/j.heliyon.2022.e11469]

29. Zhou, M.; Diao, K.; Zhang, J.; Wu, W. Controllable synthesis of plasmonic ZnO/Au core/shell nanocable arrays on ITO glass. Phys. E Low-Dimens. Syst. Nanostructures; 2014; 56, pp. 59-63. [DOI: https://dx.doi.org/10.1016/j.physe.2013.08.012]

30. Zamiri, R.; Zakaria, A.; Jorfi, R.; Zamiri, G.; Mojdehi, M.S.; Ahangar, H.A.; Zak, A.K. Laser assisted fabrication of ZnO/Ag and ZnO/Au core/shell nanocomposites. Appl. Phys. A; 2013; 111, pp. 487-493. [DOI: https://dx.doi.org/10.1007/s00339-012-7517-y]

31. Senthilkumar, N.; Ganapathy, M.; Arulraj, A.; Meena, M.; Vimalan, M.; Potheher, I.V. Two step synthesis of ZnO/Ag and ZnO/Au core/shell nanocomposites: Structural, optical and electrical property analysis. J. Alloys Compd.; 2018; 750, pp. 171-181. [DOI: https://dx.doi.org/10.1016/j.jallcom.2018.03.348]

32. Beatto, T.G.; Gomes, W.E.; Etchegaray, A.; Gupta, R.; Mendes, R.K. Dopamine levels determined in synthetic urine using an electrochemical tyrosinase biosensor based on ZnO@Au core–shell. RSC Adv.; 2023; 13, pp. 33424-33429. [DOI: https://dx.doi.org/10.1039/d3ra06277e]

33. Srivastava, V.; Gusain, D.; Sharma, Y.C. Synthesis, characterization and application of zinc oxide nanoparticles (n-ZnO). Ceram. Int.; 2013; 39, pp. 9803-9808. [DOI: https://dx.doi.org/10.1016/j.ceramint.2013.04.110]

34. Mendes, R.; Arruda, B.; de Souza, E.; Nogueira, A.; Teschke, O.; Bonugli, L.; Etchegaray, A. Determination of Chlorophenol in Environmental Samples Using a Voltammetric Biosensor Based on Hybrid Nanocomposite. J. Braz. Chem. Soc.; 2016; 28, pp. 1212-1219. [DOI: https://dx.doi.org/10.21577/0103-5053.20160282]

35. Ponnuvelu, D.V.; Pullithadathil, B.; Prasad, A.K.; Dhara, S.; Ashok, A.; Mohamed, K.; Tyagi, A.K.; Raj, B. Rapid synthesis and characterization of hybrid ZnO@Au core–shell nanorods for high performance, low temperature NO₂ gas sensor applications. Appl. Surf. Sci.; 2015; 355, pp. 726-735. [DOI: https://dx.doi.org/10.1016/j.apsusc.2015.07.143]

36. Tran, V.T.; Tran, T.H.; Le, M.P.; Pham, N.H.; Nguyen, V.T.; Do, D.B.; Nguyen, X.T.; Trinh, B.N.Q.; Van Nguyen, T.T.; Pham, V.T. et al. Highly efficient photo-induced surface enhanced Raman spectroscopy from ZnO/Au nanorods. Opt. Mater.; 2022; 134, 113069. [DOI: https://dx.doi.org/10.1016/j.optmat.2022.113069]

37. Johnson, P.B.; Christy, R.W. Optical Constants of the Noble Metals. Phys. Rev. B; 1972; 6, pp. 4370-4379. [DOI: https://dx.doi.org/10.1103/physrevb.6.4370]

38. Querry, M.R. Optical Constants. 1985; Available online: https://apps.dtic.mil/sti/citations/ADA158623 (accessed on 21 March 2025).

39. Duck, F.A. Physical Properties of Tissue: A Comprehensive Reference Book; Academic Press: London, UK, 1990.

40. Lee, K.-S.; El-Sayed, M.A. Dependence of the Enhanced Optical Scattering Efficiency Relative to That of Absorption for Gold Metal Nanorods on Aspect Ratio, Size, End-Cap Shape, and Medium Refractive Index. J. Phys. Chem. B; 2005; 109, pp. 20331-20338. [DOI: https://dx.doi.org/10.1021/jp054385p]

41. A Garcia, M. Surface plasmons in metallic nanoparticles: Fundamentals and applications. J. Phys. D Appl. Phys.; 2012; 45, 389501. [DOI: https://dx.doi.org/10.1088/0022-3727/45/38/389501]