Description

Calculation of nonlinear spectroscopy from non-perturbative propagation of density matrix defined in spectral representation.

Applications

For example, 2D third-order spectra can be obtained for photon-echo experiments.

Reference

1. T. Mančal, A. V. Pisliakov, and G. R. Fleming, The Journal of Chemical Physics 124, 234504 (2006).

2. B. Brüggemann, P. Kjellberg, and T. Pullerits, Chemical Physics Letters 444, 192 (2007).

Quick Startup

Before starting up, you need to set Makefile.in to make sure it is compilable. The options to check are CC, CFLAGS and LFLAGS. The program requires following dependencies,

• BLAS and LAPACK - Matrix inverse is used for Seidner's method
• GSL - ODE solver
• MPICH2
• libconfig - used to read parameters.cfg with nice format. But probably I will try JSON soon.

Several templates have already been there so just choose the suitable one and comment out all the others. Then you can check it by follow the next explanations step by step.

1D Nonlinear spectroscopy for single molecule (Seidner's method):

This is the option you can play on a destop machine.

• In file src/Makefile, set SRC = seidner_main_single.cc
• In file parameters.cfg, set esmb.n_esmb = 1, mvar.ny = 1, pols.ppar_calc_method = "seidner"
• Be sure you create the folder ./res to store results
• Run ./test
• After the main program exits, create folder ./fig to store image files. Then run python plot_ppar.py to see nonlinear signals in different directions. The index in file name specified in the script plot_ppar.py suggest the direction of the signal, e.g., 37 indicates the wave vector direction of [-1, 1, 1], whose correspondence can be found in file parameters.cfg. coo2 in the file name is for z-direction of the emitting field.

2D Nonlinear spectroscopy for single molecule (Seidner's method):

This option should be run on nodes with MPI avaiable. Something like 20 nodes are good. You can configure the required number of nodes by just modifying nohup mpirun -n NUMBER ./test & in file Makefile

• In file src/Makefile, set SRC = main_mpi.cc
• In file parameters.cfg, set esmb.n_esmb = 1, mvar.ny = 20, pols.ppar_calc_method = "seidner"
• Be sure you create the folder ./res to store results before running.
• Run make run
• After the main program exists, go to ana/seidner, make, then run ./test. Be sure res is also present in this folder.
• Check result by running python plot_ppar_2d.py. The naming convention is simlar like that in 1D case. Enjoy!

Code Specification

Number of equations of motion

For an n-level system (n_lvl = n), there are n terms for population ($rho_{ii}$, which has only real part), and n*(n-1)/2 terms contributing to a half coherence part ($rho_{ij}, i \neq j$, which has real and imaginary parts), the other half are just the complex conjugates of the former. Therefore, if the equations are all for real variables, the number of equation of motions is given by n_eom = n + 2 * (n*(n-1)/2) = n * n, as is defined in para.cc.

Array of polarization:

The array for generic polarization is prepared or cleaned in src/polar.cc. The polarization for Seidner's method has an extra rank in order to distinguish results from different spatial phases:

• Generic: pol[nx][n_dpl][n_dim]

• Seidner: pol[n_phase][nx][n_dpl][n_dim]

The variable nx may include all varying variables. In the simplest case, nx = nt, which is just a 1D grid along the detection time. For 2D case, the coherence time $\tau$ is also varied for typical 2D spectroscopy, nx = nt * n_tau.

Directions of field vectors

The k-unit-vector in lab frame if you define the euler angle (phi, theta, psi):

$$ kuvL_x = + sin(\phi) sin(\theta) $$

$$ kuvL_y = - cos(\phi) sin(\theta) $$

$$ kuvL_z = cos(\theta) $$

The E-unit-vector in lab frame if you define angle chi relative to the k-unit-vector:

$$ EuvL_x = + sin(\phi) cos(\theta) cos(\psi+\chi) + cos(\phi) sin(\psi+\chi) $$

$$ EuvL_y = - cos(\phi) cos(\theta) cos(\psi+\chi) + sin(\phi) sin(\psi+\chi) $$

$$ EuvL_z = - sin(\theta) cos(\psi+\chi) $$

So in order to control field polarization, If you want EuvL_x = 0: $$ \chi = -\psi + atan( - cos(\theta) * tan(\phi) ) $$

If you want EuvL_y = 0: $$ \chi = -\psi + atan( + cos(\theta) / tan(\phi) ) $$

If you want EuvL_z = 0: $$ \chi = -\psi + \pi / 2 $$

Config Data Format (CDF)

The config file may contain a chunk of data that is not easy to be expressed with key-value pairs as shown in parameters.cfg. So we introduced another config file also with the suffix name .cfg. It allows defining data directly in the file. It is useful if you have input parameters from output directly generated by other programs. In order to read config files with different formats, one needs to define their formats and transfer data to para structure by invoking sscanf function. The example/template can be found in file repr_generic.cc. Two steps are needed:

• Call cdf_read( file_name, n_lines, sscanf_structure, ps )
• Define the callback funciton sscanf_structure

The callback functions sscanf_st_energy and sscanf_st_dipole show how energy and dipole are defined. In the routine, it is also a good place to check if config data has any problem.

The current format in each line for energy and dipole data are as following,

• energy.cfg (2 columns):

1	2
level_index	energy_value

• dipole.cfg (5 columns):

1	2	3	4	5
level_index_a	level_index_b	dipole_x	dipole_y	dipole_z

• Each line in configuration file should contain characters less than 1024 as defined in cdf.cc.
• One can include blank or comment lines (// like in C++) in config files which are neglected.
• The number of data lines (not including blank or comment lines) in file energy.cfg should agree with n_lvl defined in parameters.cfg. Similarly, the number of data lines in file dipole.cfg should agree with n_dpl defined in parameters.cfg. If not, function cdf_read will complain.

Redfield Tensor

the following restrictions are imposed, but can be extended for further work:

• real part of half-sided Fourier transform of bath correlation function (neglect Lamb shift) that the tensors are all real.

• secular approximation; for R[m][n][k][l], only for population transfer (m=n, k=l) and coherence dephasing (m=k, n=l). Both the EOMs and redfield tensor generator only consider these terms. Going beyond the secular approximation, the exp( i (E_mn - E_kl) ) should be multiplied that the tensor becomes complex variable, the population terms formally are not real (although physically they are).

• the system operators, which requires representation transform, is not yet fully implemented. Here only consider the system coordinate directly defined in each electronic state, there is no overlapping between electronic states for Redfield relaxation: delta function is used, only consider terms r==l==n && s==m==k. Indices (l, n) and (m, k) are for spectral representation, r and s are indices for basis of representation where the system operator for coupling is defined. E.g., if the system operator S is defined in local site representation labeled by r, after transformation, it becomes S_ln in spectral representation, S_ln = sum_{ll,nn}( S_{ll,nn} * c_ll^l* c_nn^n ), where c_ll^l and c_nn^n are coefficients for representation transform between local site (basis labeled by ll) and exciton representation (basis labeled by l). If the basis between site and exciton ones are not orthogonal, overlapping leads to population transfer between states in the spectral representation.

TODO List

• let main_mpi.cc work for repr.type = GENERIC, and try population time loop
  • temporarily change to 2-level system with 1 dipole transition: energy.cfg, dipole.cfg and parameters.cfg
  • 2D calcuation: loop over tau is changed to loop over T in mvar.cc
• add dacay rate for new ODE
• CDF may use absolute path, so I cannot use it in subdir, adjust this.
• index with JSON format
• parameters with JSON format

FIXME

(2014/10/29)

• mvar.cc: mvar_update remember to update the time intervals properly.
• parameters.cfg:
  • 5-order -> 3-order, remember to use 5th order for calculation.
  • n_esmb = 40000
• src/Makefile: flag DUSE_MPI