Tutorial:Optical Spectra from TD
From OctopusWiki
In this tutorial we will learn how to obtain the absorption spectrum of a molecule from the explicit solution of the time-dependent Kohn-Sham equations. We chose as a test case methane (CH4).
Contents |
Ground state
Before starting out time-dependent simulations, we need to obtain the ground state of the system. For this we use basically the same inp file as in a previous tutorial:
CalculationMode= gsUnits= eV_angstromradius= 4spacing= 0.175 CH = 1.2 %Coordinates"C" | 0 | 0 | 0 "H" | CH/sqrt(3) | CH/sqrt(3) | CH/sqrt(3) "H" | -CH/sqrt(3) |-CH/sqrt(3) | CH/sqrt(3) "H" | CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3) "H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3) %
After running Octopus, we will have the Kohn-Sham wave-functions of the ground-state in the directory restart/gs. As we are going to propagate these wave-functions, they have to be well converged. It is important not to converge the energy (that is relatively easy to converge), but the density. The default ConvRelDens = 1e-03 is usually enough, though.
Time-dependent run
To calculate absorption, we excite the system with a very short electric pulse, and then propagate the time-dependent Kohn-Sham equations for a certain time T. The spectrum can then be evaluated from the time-dependent dipole moment.
The input file then reads:
CalculationMode= tdUnits= eV_angstromfromScratch= yesradius= 4spacing= 0.175 CH = 1.2 %Coordinates"C" | 0 | 0 | 0 "H" | CH/sqrt(3) | CH/sqrt(3) | CH/sqrt(3) "H" | -CH/sqrt(3) |-CH/sqrt(3) | CH/sqrt(3) "H" | CH/sqrt(3) |-CH/sqrt(3) | -CH/sqrt(3) "H" | -CH/sqrt(3) | CH/sqrt(3) | -CH/sqrt(3) %TDDeltaStrength= 0.01TDPolarizationDirection= 1 tmax = 10TDEvolutionMethod= aetrsTDTimeStep= 0.002TDMaximumIter= tmax/TDTimeStepTDPolarizationEquivAxes= 3
The first part of the input file is basically just a copy of the one used for the ground-state. The interesting part comes afterwards. First of all we define our perturbation. Its strength is given by TDDeltaStrength. This number should be small to keep the response linear, but should be sufficiently large to avoid numerical problems. The variable TDPolarizationDirection sets the polarization of our perturbation to be on the first axis, and we state that all three axes of the molecule are equivalent (see below for more details).
We then set up the time-evolution. Octopus has a large number of possible propagators that you can use (see J. Chem. Phys 121 3425-3433 (2004) for an overview). Here we use the simple Approximately Enforced Time-Reversal Symmetry (aetrs). Then we need a time step (variable TDTimeStep). The choice of this number depends on the propagation method, and it is a quite crucial parameter. You should have learned how to set it up in the previous tutorial. Finally, we set the number of time steps with the variable TDMaximumIter. To have a maximum propagation time of 10 
we will need 5000 iterations.
Note that you will be calculating the singlet dipole spectrum. You can also obtain the triplet by using TDDeltaStrengthMode, and other multipole responses by using TDKickFunction.
You can now start Octopus and go for a coffee. Propagations are slow, but the good news is that they scale very well with the size of the system. This means that even if methane is very slow, a molecule with 200 atoms can still be calculated without big problems.
Symmetries
The dynamic polarizability (related to optical absorption cross-section via 
) is, in its most general form, a 3x3 tensor. The reason is that we can shine light on the system polarized in any of the three Cartesian axes, and for each of these three cases measure how the dipole of the molecule oscillates along the three Cartesian axes. This usually means that to obtain the full dynamic polarizability of the molecule we usually need to apply 3 different perturbations along 
, by setting TDPolarizationDirection to 1, 2, or 3.
However, if the molecule has some symmetries, it is in general possible to reduce the total number of calculations from 3 to 2, or even 1. This is explained in detail here. To use this formalism in Octopus, you can use the variables TDPolarization, TDPolarizationDirection, TDPolarizationEquivAxes, and TDPolarizationWprime. The block TDPolarization defines a basis (not necessarily orthogonal) chosen to maximize the number of equivalent axes, and TDPolarizationDirection, and TDPolarizationWprime are vectors specified in that basis.
Let us give a couple of examples.
The methane molecule has Td symmetry, which means that the response is identical for all directions. This means that we only need one propagation to obtain the whole tensor. This propagation can be performed in any direction we wish. So we could use the input
%TDPolarization1 | 0 | 0 0 | 1 | 0 0 | 0 | 1 %TDPolarizationDirection= 1TDPolarizationEquivAxes= 3 %TDPolarizationWprime0 | 0 | 1 %
Note that we had omitted the blocks TDPolarization and TDPolarizationWprime in the previous input file, as these are their default values.
Now let us look at a linear molecule. In this case, you might think that we need two calculations to obtain the whole tensor, one for the direction along the axis of the molecule, and another for the axis perpendicular to the molecule. The fact is that we need only one, in a specially chosen direction, so that our field has components both along the axis of the molecule and perpendicular to it. Let us assume that the axis of the molecule is oriented along the 
axis. Then we can use
%TDPolarization1/sqrt(2) | -1/sqrt(2) | 0 1/sqrt(2) | 1/sqrt(2) | 0 1/sqrt(2) | 0 | 1/sqrt(2) %TDPolarizationDirection= 1TDPolarizationEquivAxes= 3 %TDPolarizationWprime0 | 0 | 1 %
You should try to convince yourself that the three axes are indeed equivalent and linearly independent.
Finally, let us look at a general planar molecule (in the xy plane). In principle we need only two calculations (that is reduced to one if more symmetries are present like, e.g., in benzene). In this case we chose one of the polarization axes on the plane and the other two rotated 45 degrees:
%TDPolarizationDirection1/sqrt(2) | 0 | 1/sqrt(2) 1/sqrt(2) | 0 |-1/sqrt(2) 0 | 1 | 0 %TDPolarizationEquivAxes= 2
In this case, we need two runs, one for TDPolarizationDirection equal to 1, and another for it equal to 3. Note that if there are less than 3 equivalent axes, TDPolarizationWprime is irrelevant.
Optical spectra
When the run is finished, you should find the file td.general/multipoles that contains all the information needed to compute the spectrum of methane. Now, let us rename this file to td.general/multipoles.1. Yes, you guessed, the suffix .1 has the meaning that it was the run for TDPolarizationDirection equal to 1. If we had a completely asymmetric molecule, for which we would need three time-dependent runs, the multipoles files should be called td.general/multipoles.1, td.general/multipoles.2, and td.general/multipoles.3.
To obtain the spectrum, you simply run in your working directory (not in td.general) the utility oct-cross_section (renamed to oct-propagation_spectrum in the development version). You will notice that one of two files are created: cross_section_vector.1 or cross_section_tensor. If you have all the information (either a symmetric molecule and one multipole file, or a less symmetric molecule and multiple multipole files), oct-cross-section will generate the whole polarizability tensor, cross_section_tensor. If not enough multipole files are available, it can only generate the response to the particular polarization you chose in you time-dependent run, cross_section_vector.1. Let us look first at cross_section_vector.1:
# nspin 1 # kick mode 0 # kick strength 0.010000000000 # pol(1) 1.000000000000 0.000000000000 0.000000000000 # pol(2) 0.000000000000 1.000000000000 0.000000000000 # pol(3) 0.000000000000 0.000000000000 1.000000000000 # direction 1 # Equiv. axis 3 # wprime 0.000000000000 0.000000000000 1.000000000000 #%
The beginning of the files just repeats some information concerning the run.
# Number of time steps = 5000 #SpecDampMode= 2 #SpecDampFactor= 0.1500 #SpecStartTime= 0.0000 #SpecEndTime= 4.0000 #SpecMaxEnergy= 20.0000 #SpecEnergyStep= 0.0100 #%
Now comes a summary of the variables used to calculate the spectra. Note that you can change all these settings just by adding these variables to the input file before running oct-cross-section. Of special importance are perhaps SpecMaxEnergy that sets the maximum energy that will be calculated, and SpecEnergyStep that determines how many points your spectrum will contain. To have smoother curves, you should reduce this last variable.
# Electronic sum rule = 4.820887 # Polariz. (sum rule) = 2.841141 #%
Now comes some information from the sum rules. The first is just the f-sum rule, which should yield the number of active electrons in the calculations. We have 8 electrons (4 from carbon and 4 from hydrogen), but the sum rule gives a number that is much smaller than 8! The reason is that we are just summing our spectrum up to 20 eV (see SpecMaxEnergy), but for the sum rule to be fulfilled, we should go to infinity. Of course infinity is a bit too large, but increasing 20 to a somewhat larger number will improve dramatically the f-sum rule. The second number is the static polarizability also calculated from the sum rule. Again, do not forget to converge this number with SpecMaxEnergy.
# Energy sigma(1, nspin=1) sigma(2, nspin=1) sigma(3, nspin=1) StrengthFunction(1)
# [eV] [A^2] [A^2] [A^2] [1/eV]
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.10000000E-01 0.90639407E-08 -0.86214443E-13 -0.29295335E-12 0.82578310E-08
...
Finally comes the spectrum. The first column is the energy (frequency), the next three columns are a row of the polarizability tensor, and the last one is the strength function for this run.
Let us look now at the second file, cross_section_tensor:
# nspin 1
# kick mode 0
# kick strength 0.010000000000
# pol(1) 1.000000000000 0.000000000000 0.000000000000
# pol(2) 0.000000000000 1.000000000000 0.000000000000
# pol(3) 0.000000000000 0.000000000000 1.000000000000
# direction 1
# Equiv. axis 3
# wprime 0.000000000000 0.000000000000 1.000000000000
# Energy (1/3)*Tr[sigma] Anisotropy[sigma] sigma(1,1,1) sigma(1,2,1) sigma(1,3,1) ...
# [eV] [A^2] [A^2] [A^2] [A^2] [A^2] ...
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.0000000 ...
0.10000000E-01 0.90639407E-08 0.59845845E-12 0.90639407E-08 -0.86214443E-13 -0.2929533 ...
...
The columns are now the energy, the average absorption coefficient (Tr σ / 3), the anisotropy, and then all the 9 components of the tensor. The anisotropy is defined as
![(\Delta \sigma)^2 = \frac{1}{3} * \{ 3*{\rm Tr}(\sigma^2) - [{\rm Tr}(\sigma)]^2 \}](/programs/octopus/wiki/images/math/a/2/f/a2f25b59bdead0382c01c55c0b74d8ec.png)
The reason for this definition is that it is identically equal to:
![(\Delta \sigma)^2 = (1/3) * [ (\sigma_1-\sigma_2)^2 + (\sigma_1-\sigma_3)^2 + (\sigma_2-\sigma_3)^2 ]\,](/programs/octopus/wiki/images/math/8/1/8/81833ebb1ecb8887181231a4e3b7309c.png)
where 
are the eigenvalues of 
. An "isotropic" tensor is characterized by having three equal eigenvalues, which leads to zero anisotropy. The more different that the eigenvalues are, the larger the anisotropy is.
If you now plot the average absorption coefficient, you should obtain the plot shown on the right. Of course, you should now try to converge this spectrum with respect to the calculation parameters. In particular:
- Increasing the total propagation time will reduce the width of the peaks. In fact, the width of the peaks in this methods is absolutely artificial, and is inversely proportional to the total propagation time. Do not forget, however, that the area under the peaks has a physical meaning: it is the oscillator strength of the transition.
- As you are running in a box with zero boundary conditions, you will see the box states in the spectrum. This can be improved by increasing the radius of the box. However, as one is usually interested in bound-bound transitions in the optical or near-ultraviolet regimes, a fairly small box (like the one used in this tutorial) is often enough.
.
