(Hyper-)Polarizabilities, ORD, magnetizabilities, Verdet constants¶

RESPONSE: Polarizabilities¶

The calculation of frequency-dependent (hyper)polarizabilities and related properties (Raman, ORD) is activated with the block key RESPONSE

RESPONSE
END

In this example only the zz component of the dipole polarizability tensor is calculated, at zero frequency. The orientation of the molecule is therefore crucial. Be aware that the program may modify the orientation of the molecule if the input coordinates do not agree with the symmetry conventions in ADF! (This calculation could equivalently be done through a finite field method).

See also the alternative implementation with the AORESPONSE key that offers some unique features like magnetizability, and lifetime options.

The impact of various approximations on the quality of computed polarizabilities has been studied in, for instance, Refs. [74, 82, 89]. If you are new to this application field, we strongly recommend that you study a few general references first, in particular when you consider hyperpolarizability calculations. These have many pitfalls, technically (which basis sets to use, application of the DEPENDENCY key) and theoretically (how do theoretical tensor components relate to experimental quantities, different conventions used). Please, take a good look both at ADF-specific references [75-77, 90] and at general references related to this subject: Refs. [91-93], the entire issues of Chem.Rev.94, the ACS Symposium Series #628, and further references in the ADF-specific references.

RESPONSE
  ALLCOMPONENTS
  Nfreq Nfreq
  FrqBeg FirstFreq
  FrqEnd LastFreq
  ALLTENSOR
  Quadrupole
  Octupole
  Traceless
END

Entire tensor or only one component

You specify the ALLCOMPONENTS subkey to get the entire polarizability tensor, instead of just the zz component.

Frequencies or wavelengths

Instead of performing the calculation at zero frequency (which results in the static polarizability), one can specify an even-spaced sequence of frequencies, using the subkeys Nfreq, FrqBeg, and FrqEnd with obvious meaning. The unit of FrqBeg, and FrqEnd is eV.

Higher multipole polarizabilities

Instead of just calculating the dipole-dipole polarizability, one may address the dipole-quadrupole, quadrupole-quadrupole, dipole-octupole, quadrupole-octupole, and octupole-octupole polarizability tensors. These can all be calculated in a single run, using the subkey ALLTENSOR. If only quadrupole-quadrupole or octupole-octupole tensors are needed, the subkey quadrupole or octupole should be used.

Traceless

If the subkey TRACELESS is included this ensures the calculation of tensors using the traceless form of the quadrupole operator.

RESPONSE: Accuracy and convergence¶

RESPONSE
  erralf 1e-6
  erabsx 1e-6
  errtmx 1e-6
  ncycmx 30
END

erralf, erabsx, errtmx

The subkeys erralf, erabsx, errtmx determine the convergence criteria for a polarizability calculation. The strict defaults are shown. It is rarely necessary to change the defaults, as these are rather strict but do not lead to many iterations.

ncycmx

The maximum number of attempts within which the algorithm has to converge. The default appears to be adequate in most cases.

RESPONSE: Hyperpolarizabilities¶

Hyperpolarizabilities

RESPONSE
  HYPERPOL LaserFreq
END

The first hyperpolarizability tensor \(\beta\) is calculated (in atomic units in the ‘theoreticians convention’, i.e. convention T=AB in Ref. [92]) if the subkey HYPERPOL is present with a specification of the laser frequency (in Hartree units). If also the subkey ALLCOMPONENTS is specified, all components of the hyperpolarizability tensor will be obtained.

As mentioned before, by default only the static dipole hyperpolarizability tensor is computed. If one is interested in the frequency-dependent hyperpolarizability, the input could look like:

RESPONSE
  ALLCOMPONENTS
  HYPERPOL 0.01
  DYNAHYP
END

The subkey DYNAHYP has to be added and the main frequency \(\omega\) has to be specified in Hartrees after the subkey hyperpol. In the output all nonzero components of the tensors governing the static first hyperpolarizability, second harmonic generation, electro-optic Pockels effect, and optical rectification are printed.

Note: Second hyperpolarizabilities are currently not available analytically in case of RESPONSE. Some can however be obtained by calculating the first hyperpolarizability in a finite field.

The effect of using different DFT functionals (LDA, LB94, BLYP) on hyperpolarizabilities in small molecules has been investigated in [77].

RESPONSE: Optical rotation dispersion (ORD)¶

RESPONSE
  OPTICALROTATION
END

OPTICALROTATION

With the subkey OPTICALROTATION the (frequency dependent) optical rotation [80, 94] will be calculated. For correct calculations one should calculate the entire tensor (see also the subkey ALLCOMPONENTS), which is done automatically.

An alternative implementation uses the AORESPONSE key, in which life time effects can be included.

AORESPONSE: Polarizabilities¶

If the block key AORESPONSE is used, by default, the polarizability is calculated.

AORESPONSE
  FREQUENCY      Nfreq freq1 freq2 ... freqN units
  FREQRANGE      freq1 freqN TotFreq units
  LIFETIME width
END

FREQUENCY      Nfreq freq1 freq2 ... freqN units

To calculate time-dependent properties, one needs to specify frequency of perturbation field. Here Nfreq specifies the number of frequencies that follow. The last item on the line specifies the units and is one of EV, HARTREE, ANGSTROM.

FREQRANGE freq1 freqN TotFreq units

This key is useful when it is necessary to specify more than 20 equally spaced frequencies for the response calculations. The first frequency is freq1 and the last one is freqN. The total number of frequencies including the first and the last one is TotFreq. The last item specifies the units: EV, HARTREE or ANGSTROM.

LIFETIME width

Specify the resonance peak width (damping) in Hartree units. Typically the lifetime of the excited states is approximated with a common phenomenological damping parameter. Values are best obtained by fitting absorption data for the molecule, however, the values do not vary a lot between similar molecules, so it is not hard to estimate values. A value of 0.004 Hartree was used in Ref. [266].

The spin-orbit ZORA polarizability code (Ref. [311]) is automatically selected if the AORESPONSE keyword is given in a spin-orbit coupled calculation. In this case a spin-restricted calculation is required, but, unlike the rest of AORESPONSE, also SYMMETRY NOSYM. Spin-polarization terms in the XC response kernel are neglected. In Ref. [311] the imaginary polarizability dispersion curves (spin-restricted) match well the broadened spin-orbit TDDFT data from Ref. [182]. Thus the corrections from the spin-polarization terms appear to be rather minor. No picture change corrections were applied in the ZORA formalism. AORESPONSE icw hybrids icw spin-orbit is not implemented.

AORESPONSE: Technical parameters and expert options¶

AORESPONSE
 ...
  SCF        {NOCYC} {NOACCEL} {CONV=conv} {ITER=niter}
  GIAO
  FITAODERIV
  COMPONENTS {XX} {XY} {XZ} {YX} {YY} {YZ} {ZX} {ZY} {ZZ}
  ALDA|XALPHA
  ALPHA
END

SCF {NOCYC} {NOACCEL} {CONV=conv} {ITER=niter}

Specify CPKS parameters such as the degree of convergence and the maximum number of iterations:

NOCYC - disable self-consistence altogether

NOACCEL - disable convergence acceleration

CONV - convergence criterion for CPKS. The default value is 10^-6 . The value is relative to the uncoupled result (i.e. to the value without self-consistence).

ITER - maximum number of CPKS iterations, 50 by default.

Specifying ITER=0 has the same effect as specifying NOCYC.

GIAO

Include the Gauge-Independent Atomic Orbitals (GIAO). This option should not be used with damping (LIFETIME keyword) and the VELOCITYORD option should be used instead.

FITAODERIV

Use fitted AO Derivatives. This will improve the density fitting, can only be used in cae of STO fitting. In case of ZlmFit one can improve the fitting with the ZLMFIT block key.

COMPONENTS {XX} {XY} {XZ} {YX} {YY} {YZ} {ZX} {ZY} {ZZ}

Limit the tensor components to the specified ones. Using this option may save the computation time.

ALDA|XALPHA

If ALDA is specified the VWN kernel is used. This option is the default. If ALPHA is specified the X\(\alpha\) kernel is used instead of the default VWN one.

ALPHA

Writes perturbed density matrix to TAPE16.

AORESPONSE: Damped First Hyperpolarizabilities¶

In Ref.[451] an implementation of finite lifetimes into TDDFT for the calculation of quadratic response properties in ADF is described. All \(\beta\) tensor components (27 in total) will be printed in the output.

AORESPONSE
  BETA|QUADRATIC
  FREQUENCY  2 freq1 freq2 units
  LIFETIME width
  STATIC|OPTICALR|EOPE|SHG
END

BETA

Option to calculate the damped first hyperpolarizability (\(\beta\)) using quadratic response theory based on the 2n+1 rule. Two input frequencies are required for this calculation and the property \(\beta (-(\omega_1+\omega_2); \omega_1, \omega_2)\) will be the output. Note that one can choose certain values of the two frequencies to calculate different types of \(\beta\), such as static case \(\beta(0;0,0)\), optical rectification \(\beta(0;\omega,-\omega)\), electro-optical Pockels effect \(\beta(-\omega;\omega,0)\), and second harmonic generation \(\beta(-2\omega;\omega,\omega)\). Alternatively, these can be efficiently calculated using the (sub-)keywords STATIC, OPTICALR, EOPE, and SHG, respectively. Note that the needed input frequencies all rely on \(\omega_1\) (freq1) when using the (sub)keywords above.

QUADRATIC

This option possess the same functionality as BETA, i.e., calculate the damped \(\beta\), except not adapting the 2n+1 rule. The (sub)keywords STATIC, OPTICALR, EOPE,and SHG are also applicable here. Note that this approach facilitates the direct partitioning of the response into contributions from localized orbitals and is important for natural bond analysis.

Note: Please only use HARTREE or EV as the unit for the input frequencies. The unit option ANGSTROM does not work correctly due to the current implementation structure.

AORESPONSE: Damped Second Hyperpolarizabilities¶

In Ref.[452] a general implementation for damped cubic response properties is presented using time-dependent density functional theory (TDDFT) and Slater-type orbital basis sets. To directly calculate two-photon absorption (TPA) cross sections, an implementation of a reduced damped cubic response approach is described in Ref.[452]. All \(\gamma\) tensor components (81 in total) will be printed in the output.

AORESPONSE
  GAMMA|CUBIC
  FREQUENCY  3 freq1 freq2 freq3 units
  LIFETIME width
  STATIC|EFIOR|OKE|IDRI|EFISHG|THG|TPA
END

GAMMA

Option to calculate the damped second hyperpolarizability (\(\gamma\)) using cubic response theory based on the 2n+1 rule. Three input frequencies are required for this calculation and the property \(\gamma(-(\omega_1+\omega_2+\omega_3); \omega_1, \omega_2, \omega_3)\) will be the output. Note that one can choose certain values of the three frequencies to calculate different types of \(\gamma\), such as static case \(\gamma(0;0,0,0)\), electric field induced optical rectification \(\gamma(0;\omega,-\omega,0)\), optical Kerr effect \(\gamma(-\omega;\omega,0,0)\), intensity dependent refractive index \(\gamma(-\omega;\omega,\omega,-\omega)\), electric field induced second harmonic generation \(\gamma(-2 \omega;\omega,\omega,0)\), and third harmonic generation \(\gamma(-3 \omega;\omega,\omega,\omega)\). Alternatively, these can be efficiently calculated using the (sub)keywords STATIC, EFIOR, OKE, IDRI, EFISHG, and THG, respectively. The (sub)keyword TPA can be used to calculate the \(\gamma\) corresponding to the two photon absorption process (i.e., the reduced form of \(\gamma(-\omega;\omega,\omega,-\omega)\)), however, it can ONLY be used with keyword GAMMA. Note that the needed input frequencies all rely on \(\omega_1\) (freq1) when using the (sub)keywords above.

CUBIC

This option possess the same functionality as GAMMA, i.e., calculate the damped \(\gamma\), except not adapting the 2n+1 rule. The (sub)keywords STATIC, EFIOR, OKE, IDRI, EFISHG and THG are also applicable here. Note that this approach facilitates the direct partitioning of the response into contributions from localized orbitals and is important for natural bond analysis.

Note: Please only use HARTREE or EV as the unit for the input frequencies. The unit option ANGSTROM does not work correctly due to the current implementation structure.

AORESPONSE: ORD¶

AORESPONSE
  OPTICALROTATION
  VELOCITYORD
  FREQUENCY      Nfreq freq1 freq2 ... freqN units
  FREQRANGE      freq1 freqN TotFreq units
  LIFETIME width
END

OPTICALROTION

Specify OPTICALROTION to calculate optical rotatory dispersion spectrum instead of polarizabilities.

VELOCITYORD

This option should be used instead of OPTICALROT with GIAO if the finite lifetime effects need to be taken into account (LIFETIME option).

AORESPONSE: magnetizabilities, Verdet constants, Faraday B term¶

AORESPONSE
  MAGNETICPERT
  MAGOPTROT
  FREQUENCY      Nfreq freq1 freq2 ... freqN units
  FREQRANGE      freq1 freqN TotFreq units
  LIFETIME width
END

MAGNETICPERT

Calculate static or time-dependent magnetizability, see also Ref. [230].

MAGOPTROT

Specify MAGOPTROT to calculate the Verdet constant instead of polarizability, see for the details of the implementation Ref. [307]. When it is specified together with the LIFETIME key the real and imaginary part of the damped Verdet constant will be calculated. Combination of three keys MAGOPTROT, LIFETIME and FREQRANGE yields the magnetic optical rotatory dispersion and magnetic circular dichroism spectrum (Faraday A and B terms) calculated simultaneously in the range from freq1 to freqN. It is also possible to combine MAGOPTROT, LIFETIME and FREQUENCY. In order to obtain the Faraday B terms from the Verdet constant calculations it is necessary to perform several steps, involving a fit of the imaginary Verdet data to the MCD spectrum. You can request SCM for details on the fitting procedure. For details of the method, see Ref. [272].

AORESPONSE: Raman¶

AORESPONSE
  RAMAN
  FREQUENCY  1  freq1 units
  LIFETIME width
END

RAMAN

Calculates the Raman scattering factors. The AOREPONSE-Raman only works with one frequency. If one frequency is specified the Raman scattering factors are calculated at that frequency. The Raman option is compatible with the lifetime option so that resonance Raman scattering can be calculated. For details of this method, see Ref. [266]. To get Raman intensities with AORESPONSE, numerical frequencies need to be calculated using a FREQUENCIES key in the GEOMETRY input block. Non-resonance Raman intensities can also be obtained using the RESPONSE key or, alternatively, using RAMANRANGE in combination with analytically or numerically pre-calculated frequencies.

Applications of AORESPONSE¶

It may be useful to consult the following applications of the AORESPONSE key in ADF:

• Calculation of static and dynamic linear magnetic response in approximate time-dependent density functional theory [230]
• Calculation of CD spectra from optical rotatory dispersion, and vice versa, as complementary tools for theoretical studies of optical activity using time-dependent density functional theory [231]
• Calculation of origin independent optical rotation tensor components for chiral oriented systems in approximate time-dependent density functional theory [232]
• Time-dependent density functional calculations of optical rotatory dispersion including resonance wavelengths as a potentially useful tool for determining absolute configurations of chiral molecules [233]
• Calculation of optical rotation with time-periodic magnetic field-dependent basis functions in approximate time-dependent density functional theory [234]
• A Quantum Chemical Approach to the Design of Chiral Negative Index Materials [235]
• Calculation of Verdet constants with time-dependent density functional theory. Implementation and results for small molecules [236]
• Calculations of resonance Raman [266,267]
• Calculations of surface-enhanced Raman scattering (SERS) [268,269]
• Calculation of magnetic circular dichroism spectra from damped Verdet constants [272]
• Calculation of the polarizability in case of spin-orbit coupling [311]

UV/Vis spectra, CD spectra¶

If one includes the POLTDDFT keyword the (real and imaginary part of the) diagonal of the polarizability tensor and rotatory strengths will be calculated, which can be used to calculate the photoabsorption and circular dichroism (CD) spectra.

POLTDDFT
  IRREP
     Irrep1
     Irrep2
     ...
  SUBEND
  KGRID eVgrid ngrid
  FREQRANGE eVi eVf nfreq
  LIFETIME eVwidth
  CUTOFF eVcutoff
  LAMBDA lambda
  VELOCITY
END

IRREP

Subblock key for selecting which symetry irreps of the excitations to calculate (all excitations by default). In the subkey data block list the symmetry irrep labels, like B1, for example.

KGRID eVgrid ngrid

Keyword KGRID is used to discretize the energy scale for calculating the complex dynamical polarizability. Only pairs of an occupied and virtual orbital are included, for which the orbital energy difference is lower than eVgrid (9 eV by default). Ngrid is the number of points within the energy grid (180 by default).

FREQRANGE eVi eVf nfreq

Keyword FREQRANGE is used to specify the equally spaced points in the spectrum for which one would like to calculate the complex dynamical polarizability. The first point is eVi (0 eV by default) and the last one is eVf (5 eV by default). The total number of points in the spectrum is nfreq (100 by default).

LIFETIME eVwidth

Specify the resonance peak width (damping) in eV units. Typically the lifetime of the excited states is approximated with a common phenomenological damping parameter. Values are best obtained by fitting absorption data for the molecule, however, the values do not vary a lot between similar molecules, so it is not hard to estimate values. Default value is 0.1 eV.

CUTOFF eVcutoff

For a given point in the spectrum, only include pairs of an occupied and virtual orbital, for which the orbital energy difference is lower than the energy of the point in the spectrum plus eVcutoff. The default value for eVcutoff is 4 eV.

LAMBDA lamba

Jacob’s scaling factor [482] for the study of plasmonic resonances. This factor, 0<lambda<1 (default lambda=1), turns on the coupling matrix K.

VELOCITY

If the subkey VELOCITY is included ADF calculates the dipole moment in velocity gauge. Default the dipole-length representation is used.

Reduced fit set¶

For POLTDDFT it very important is to use specially made auxiliary fit sets. These are available only for a very limited amount of elements, and for a limited amount of basis sets. These special basis sets can be found in the atomic database directories, with a directory name REDFIT, for example, the directory $ADFHOME/atomicdata/ZORA/TZP/REDFIT. An example:

Basis
  Au ZORA/DZ/REDFIT/Au.4f
End

DISPER program: Dispersion Coefficients¶

The DISPER program was originally written by V.Osinga [79]. The original documentation was written by S.J.A. van Gisbergen.

Van der Waals dispersion coefficients

The program DISPER computes Van der Waals dispersion coefficients up to C₁₀ for two arbitrary closed-shell molecules. ADF itself can already compute some C₆ and C₈ coefficients between two identical closed-shell molecules. These coefficients describe the long-range dispersion interaction between two molecules. It requires previous ADF-TDDFT calculations for the polarizability tensors at imaginary frequencies for the two interacting molecules. Each such ADF calculation produces a file TENSOR (if suitable input for ADF is given). The TENSOR files must be renamed tensorA and tensorB, respectively and must be present as local files for DISPER. The DISPER program takes no other input and prints a list of dispersion coefficients.

A schematic example, taken from the set of sample runs, for the usage of DISPER is the following:

Step1: run ADF for, say, the HF molecule. In the input file you specify the RESPONSE data block:

RESPONSE
  MaxWaals 8    ! Compute dispersion coefficients up to C8
  ALLTENSOR     ! This option must be specified in the ADF calc for a
                ! subsequent DISPER run
  ALLCOMPONENTS ! Must also be specified for DISPER
End

At the end of the run, copy the local file ‘TENSOR’ to a file ‘tensorA’. For simplicity, we will now compute the dispersion coefficients between two HF molecules. Therefore, copy ‘tensorA’ to ‘tensorB’.

Now run DISPER (without any other input). It will look for the local files ‘tensorA’ and ‘tensorB’ and compute corresponding dispersion coefficients to print them on standard output.

$ADFBIN/disper -n1 << eor
eor

The output might look something like this:

DISPER 2000.02 RunTime: Apr04-2001 14:14:13
  C-COEFFICIENTS
  n   LA  KA  LB  KB  L  coefficient(Y)    coefficient(P)
  6   0   0   0   0   0   28.29432373         28.29432373
  6   2   0   0   0   2    7.487547697         3.348533127
  8   0   0   0   0   0  416.1888455         416.1888455
  8   0   0   2   0   2    0.4323024202E-05    0.1933315197E-05
  8   2   0   0   0   2  402.3556946         179.9389368
  8   2   0   2   0   4    0.4238960180E-05
  8   4   0   0   0   4  -36.67895539        -12.22631846
  8   4   0   2   0   6   -0.2000286301E-05

The n-value in the first column refers to the long-range radial interaction. The case n=6 refers to the usual dipole-dipole type interaction related to a 1/R⁶ dependence in the dispersion energy. The n=7 case relates to a dipole-quadrupole polarizability on one system and a dipole-dipole polarizability on the other (this is not symmetric!). The n=8 term may contain contributions from a quadrupole-quadrupole polarizability on one system in combination with a dipole-dipole polarizability on the other as well as contributions from two dipole-quadrupole polarizabilities.

Terms which are zero by symmetry are not printed. In the example above, this is the case for all n=7 terms, because the systems (apparently) are too symmetric to have a nonzero dipole-quadrupole polarizability. The best known and most important coefficients are the isotropic ones, determining the purely radial dependence of the dispersion energy. They are characterized by the quantum numbers: 6 0 0 0 0 0 (or 8 0 0 0 0 0 etc.) Other combinations of quantum numbers refer to different types of angular dependence. The complete set determines the dispersion energy for arbitrary orientations between the two subsystems A and B.

The complete expressions are rather involved and lengthy. We refer the interested reader to the paper [79] which contains a complete description of the meaning of the various parts of the output, as well as references to the earlier literature which contain the mathematical derivations. In particular, a useful review, which was at the basis of the ADF implementation, is given in [414]. Of particular significance is Eq.(8) of the JCP paper mentioned above, as it defines the meaning of the calculated coefficients \(C_n^{L_A, K_A, L_B, K_B, L}\) as printed above.

For highly symmetric systems, a different convention is sometimes employed. It is based on Legendre polynomials (hence the ‘P’ in the final column) instead of on the spherical harmonics (the ‘Y’ in the column before the last). The ‘P’ coefficients are defined only for those coefficients that are nonzero in highly symmetric systems and never contain additional information with respect to the ‘Y’ coefficients. They are defined [Eq. (14) in the mentioned J. Chem. Phys. paper] in terms of the ‘Y’ coefficients by:

Because the quality of the dispersion coefficients is determined by the quality of the polarizabilities that are the input for DISPER, it is important to get good polarizabilities from ADF. For that it is important, in the case of small systems, to use an asymptotically correct XC potential (several choices are available in ADF, such as SAOP or GRAC) and a basis set containing diffuse functions. We refer to the ADF User’s Guide for details.