November 13th, 2011TutorialsJuan Sebastian Totero Gongora 0 Comments

In place, Fourier analysis of fields and energy patterns

At the basic level, Nanocpp provides a comprehensive FDTD simulator that can be used by properly setting the *.in files. At a higher level, Nanocpp can be used as an FDTD library to perform more advanced electrodynamics parallel calculations. In these example, we will study how to enable a in place, Fourier analysis of electromagnetic fields. This is extremely useful as it allows with a single parallel simulation to investigate the spectral properties of a system at an arbitrary number of different frequencies.

In order to use Nanocpp as a Library, we need to modify the file located in your simulation directory (default to <nanocpp root>/newsim).

The standard file is as follows:

using namespace std;

int main(int argc, char* argv[]){

////////////////////// INIT THE CODE
nanocpp n(argc,argv);
cout<<n<<endl; //to screen

///////////////////// CYCLE

double sta=MPI_Wtime();


  n.time=n.tcurr*n.dt;//update time;

  n.single_step();//single step


//////////////////// SAVE OUTPUT AND LOGOUT /////////////

n.post_process(); //process data
n.goodbye(sta);   //get time and print goodbye



The file begins by including the nanocpp class, and then init the code:

nanocpp n(argc,argv);
cout<<n<<endl; //to screen

in this phase the files *.in are read, and the necessary FDTD data is set up. After the initialization, the FDTD main cycle starts:

double sta=MPI_Wtime();
  n.time=n.tcurr*n.dt;//update time;
  n.single_step();//single step

The for loop execute the FDTD single_step() routine several times, up to the reaching of the timestep chosen in the file. After the FDTD simulation, data is processed and the program terminates:

n.post_process(); //process data
n.goodbye(sta); //get time and print goodbye

In order to enable the DFT analysis in Nanocpp, we start by including the corresponding class:


This module allows to define a volume 3D (or a surface in 2D), where the Fourier Transform (FT) of the electromagnetic fields if calculated during the simulation, without the need of saving probes in the time domain. This requires to save a considerably amount of time, as well as of memory. The FT is computed on the fly by a discrete trapezoidal integration method:

\tilde{\psi}(\omega)=\int e^{i\omega t}\psi(t)dt=\sum_n e^{i\omega ndt}\psi(ndt)

We begin by specifying the volume where the FT analysis needs to be computed. To this extent, we define a volume object vv an initialize it:

double cds[2*FD_N]={0.5,1.5,0.5,1.5};
int fc[FD_N]={CE,CE};
volume<FD_N> vv;


FD_N is a macro and contains the FDTD dimensions, in this 2D example FD_N=2. The volume class init constructor requires:

  1. an array of boundaries of the type xsta,xend,ysta,yend,... In this example we choose a surface limited by 0.5 to 1.5 in both x and z.
  2. an array of faces where the surface is centered. The DFT class works by definition at the YEE cell center, so we need to specify CE (=center) along every direction.

For more info on the volume class, look at the file “surfaces.h“, which contains several parallel primitives, including boxes and planes.

Once the volume has been defined, we initialize a DFT object:

double fw[20];
for(int i=0;i<20;i++)
int wsiz=20;
double dt=n.get_dt();
dft_vol df(vv,fw,wsiz,dt);

This operation requires:

  1. the definition of an array of frequencies (fw) and its relative length (wsiz).
  2. the FDTD dt integration step, available via the public nanocpp function get_dt().
  3. a volume (vv).

We then need to tell Nanocpp to accumulate FT data after each time step. To this extent, we modify the standard loop as follows:

n.single_step(); //single step
n.erg_exch(); //exchange boundaries·
df.accum_dft_all(n,n.tcurr);   // accumulate dft results

After the standard FDTD step, we call:

1) the routine erg_exch(), which exchanges the necessary boundaries for the computation of the field at the YEE cell center

2) accum_dft_all(), which computes the in place FT transform for every field. The routine is optimized and takes very little time.

When the simulation is completed, we can calculate and save the energy:

string fn="res/data.bin";

as well as the real part of the electromagnetic fields (for the imaginay part, just add the string “imag” as a final argument to the function save_fields):


in this case we need to only specify a directory (“res“) where to save data. Each field will be saved as <fn>/Xiw.bin, with X=E or H and i=x,y or z. For more info on the class DFT, please refer to the file “dft_in_place.h“.

In the following example, generated by the files, and, we perform a parallel FT analysis, with 20 different frequencies, of a TFSF source in air. In order to compile the code, type:

make clean;make

and then execute it via the standard mpiexec command. This is the output analyzed through Matlab:


Fig.1: Energy(x=1,z=1,λ)


Figure 2. Ey field at λ=0.25


Figure 3. Energy (x,z) distribution at·λ=0.25 in a surface delimited by (0.5,1,0.5,1) with a TFSF region 1 3 1 3

For a plane wave of the type

\psi\propto \cos(\omega_0 t -k z),

we have

\hat{\psi}(\omega>0)\propto\delta(\omega-\omega_0)\exp^{i k z}

as perfectly verified in the figures above.

Optical Forces Calculation


NANOCPP includes a set of routines designed to compute the electromagnetic force acting on an object in a 2D system. Starting from the Maxwell Stress Tensor formulation, the average force vector acting on an object is given by

\vec{F}=\frac{1}{T} \int_0^T dt \int_S dS\,\bar{T}\cdot\hat{n}


 T_{ij} \equiv \epsilon_0 \left( E_i E_j - \frac{1}{2} \delta_{ij} E^2\right) + \frac{1}{\mu_0} \left(H_i H_j - \frac{1}{2} \delta_{ij} H^2\right)

Here S is a generic surface surrounding the object of interest and T_{ij} is the Maxwell Stress Tensor. This means that to compute the electromagnetic force, we must compute the flux of \Bar{T} across the surface S for every t\in[0,T], and then compute its time average.

In NANOCPP, the flux of the Maxwell Tensor across a surface S is given, at fixed time, by the function


that takes as arguments a box, representing the integration surface, and an integer (here \hat{x} corresponds to 0, and \hat{z} to 1) that selects the force component to be computed.

Consider a system in which a 40nm silica sphere (n=1.5) is placed in a quasi-gaussian field, generated using an opportune TE-polarized tfsf-single source(see attached and files). The particle position is (2.15e^{-6}, 2e^{-6}) and we want to enclose it in a 0.1\ \mu m \times 0.1\ \mu m box.


Total energy after 1000 timesteps

In our we have to create the integration surface as

//create a volume

double vol_coords[2*FD_N]={2.1e-6,2.2e-6,1.95e-6,2.05e-6};
int vol_faces[FD_N]={CE,CE};
volume<FD_N> vol;

//create a box

double mycoords[2*FD_N];
int myfaces[FD_N*FD_N*FD_N];
for(int i=0;i<FD_N,i++){

for(int i=0;i<FD_N*FD_N*FD_N,i++)

box<FD_N> mybox;


This way, mybox is a parallel surface, i.e. it is automatically distributed on all the available processors.To store the flux of the Maxwell Tensor time by time, we create a dynamic array using the class tensor

tensor<1,double> fx(&n.total_steps,0.)

and procede to insert inside the time cycle

/////////////////// MY CODE ////////////////////
fx[i]=n.get_steady_optical_force(mybox,0); // force along x

At the end of the cycle, every processor will store its part of the force, so we need to perform an opportune MPI_Reduce on the processor 0


where gsum is a tensor belonging to the processor 0, that can perform the average over time simply as


The variable force, at this point, contains the force expressed in Newton (in the case under exam we obtained Fx=-0.33pN).

Repeating the simulations changing the position along x of the scatterer, one obtains the following


where the value of the force at x=0, and the negative sign of the derivative of F(x) in x=0 suggests the presence of an optical trap in the x direction (as expected in the case of a gaussian beam).

Maxwell-Bloch Equations

Nanocpp encompasses a set of fast functions for solving the full vector set of Maxwell-Bloch equations ins 2D/3D. Theory in 3D is based on A. Fratalocchi et al., PRA 78, 013806 (2008), while 2D implementations are described in details in G. Slavcheva et al., PRA 66, 063418 (2002).

In order to define a material supporting gain, we need to specify in the file the parameters


(please see the papers above for a definition of the quantities).

@ 1 2.62 1e26 6.2788e15 1e-10 1e-13 1e-14

which correspond to a nondispersive gain material with

\varepsilon_r=2.62 and  N_a=10^{26}\omega_0=6.2788 \cdot 10^{15}q_0=0.1 nmT_0=100 fsT_1=10 fs.

Dispersive gain can be introduced by expressing the poles first and then the parameters for the gain medium.

In the file, we then need to enable the Maxwell-Bloch functions:


By default Nanocpp enables pumping fields on all polarizations (Ex, Ey and Ez in 3D, Ex and Ez in 2D). In order to specify a single polarization the following function can be used in the

//init nanocpp object
nanocpp n(argc,argv);
// field type is "Ex", "Ey" or "Ez".
n.set_mb_field(string <field_type>,int <material id>);

Quantum noise is provided by an extremely fast random number generator based on the Mersenne-Twister algorithm, which is characterized by a practically infinite period of 2^{19937}-1. The noise variance of the electric field is specified in the code by the noise_std value in the file:


the following and simulate a cylindrical resonator laser for TM polarized waves:


Fig. 1. Energy snapshot at t=6ps


Fig. 2. Dynamics of Ex versus time

 Advanced Absorption Calculations

In this example we are going to analyze the power absorption in a metal slab, using different methods. Files used in this tutorial:

To calculate the Pointing flux along a closed box, NANOCPP has the function

double erg_flux(const box<FD_N>& vol); 

This function calculates the Pointing flux over all the sides of the box vol. Here is an example of initialization of a box in 3D for energy computations.

double cds[6]={1.5e-6,2.5e-6,1.5e-6,2.5e-6,1.5e-6,2.5e-6};

for(int i=0;i<2*FD_N*3;i++)

box<FD_N> kkk;

For more information, please refer to the class box defined in surfaces.h.

Casimir Forces computation

(In preparation).