10.13 Finite-Field Calculation of (Hyper)Polarizabilities

10.13.2 Romberg Finite-Field Procedure

Whereas the FF procedure described in Section 10.13.1 is a straightforward, finite-difference implementation of the derivatives suggested in Eq. (10.61), in the Romberg procedure210 one combines energy values obtained for a succession of k external electric fields with amplitudes that form a geometric progression:

F(k)=akF0. (10.62)

The FF expressions are obtained by combining truncated Taylor expansions of the energy with different amplitudes and/or external field directions. For example, in the case of the diagonal 𝜷-tensor components the Romberg FF expression is

βiii(k,0)=3([E(F-i(k+1)-E(Fi(k+1))]-a[E(F-i(k)-E(Fi(k))]a(a2-1)(akF0)3). (10.63)

The field index ±i refers to the possible field directions, i.e., ±x, ±y or ±z. Truncation of the Taylor expansions means that the results are contaminated by higher-order hyperpolarizabilities, and to remove this contamination, successive “Romberg iterations” are performed using a recursive expression. For a component of a (hyper)polarizability tensor 𝜻, the recursive expression is

ζ(k,n)=a2nζ(k,n-1)-ζ(k+1,n-1)a2n-1. (10.64)

This expression leads to a triangular Romberg table enabling monitoring of the convergence of the numerical derivative.210

As with any finite-difference procedure, the FF method for computing (hyper)polarizabilities is sensitive to the details of numerical differentiation. The Romberg procedure allows one to find a field window, defined by its upper and lower bounds, where the finite-difference procedure is stable. Energy values for field amplitudes below that window suffer from too-large round-off errors, which are proportional to the energy convergence threshold. The upper bound is imposed by the critical field amplitude corresponding to the intersection between the ground and excited-state energies. In the Romberg procedure, this stability window defines a sub-triangle, determination of which is the primary goal in analyzing Romberg tables.

The automatic procedure based on the analysis of field amplitude errors is implemented in scripts provided with Q-Chem’s distribution. The field amplitude error is defined as the difference between ζ-values obtained for consecutive field amplitudes at the same Romberg iteration:

ϵk(n)=ζ(k+1,n)-ζ(k,n). (10.65)

By virtue of Romberg’s recursive expression, the field error is expected to decrease with each iteration. Convergence of the Romberg procedure can be probed using the iteration (order) error, defined as

ϵn(k)=ζ(k,n+1)-ζ(k,n). (10.66)

Automatic analysis of these quantities is described in detail in Ref. 210.

Note:  The automatic procedure can fail if the field window is not chosen wisely. The Romberg procedure can be performed either for one specific diagonal direction or for all Cartesian components. In the case of the second hyperpolarizability, only the iiii, iijj components are available.

10.13.2.1 How to execute Romberg’s differentiation procedure

To perform Romberg calculations of (hyper)polarizabilities, Q-Chem provides the following scripts in the $QC/bin/Romberg directory: input-Q-Chem-t-rex.sh parse-t-rex.sh input-Q-Chem-t-rex-3.0.f90 tddft_read.f90 eom_read.f90 T-REX-3.0.3.f90 To set up a calculation, copy these files into your home directory and compile Fortran files (*.f90). The executables should be named input-Q-Chem-t-rex-3.0, tddft_read, eom_read, and T-REX. The compilation command is gfortran file.f90 -o file After compilation, put the binaries in the install directory (which can be the same as the directory with *.sh and inputs), and add this directory to the $PATH variable:

bash syntax:
export PATH="$PATH":full_path_to_the_fortran_dir
csh/tsch syntax:
set path=($path full_path_to_the_fortran_dir)

(This line can be added to your .bashrc/.cshrc file for future runs.)

To run the calculation, create an input that specifies your molecule, an electronic structure method, and several additional $rem variables that (i) turn off symmetry (SYMMETRY and CC_SYMMETRY for CC/EOM calculations), (ii) request higher-precision printing (CC_PRINT_PREC), and (iii) set up very tight convergence (SCF_CONVERGENCE, CC_CONVERGENCE, EOM_DAVIDSON_CONVERGENCE, etc.). An example of an input file is given below.

Example 10.34  Input for Romberg calculations of ethylene molecule using B3LYP.

$molecule
0 1
 C        0.00000000   0.00000000 -0.66880000
 H        0.94859916   0.00000000 -1.19917145
 H       -0.94859916   0.00000000 -1.19917145
 C       -0.00000000   0.00000000  0.66880000
 H       -0.54409413   0.77704694  1.19917145
 H        0.54409413  -0.77704694  1.19917145
$end

$rem
   BASIS            =  sto-3g
   EXCHANGE         =  B3lyp
   SCF_CONVERGENCE  =  13   Need tight convergence for finite field calculations
   SCF_MAX_CYCLES   =  200
   SYMMETRY         =  false    All symmetries need to be turned off
   CC_PRINT_PREC    =  16  16 decimal points of total energies will be printed
$end

Run the script input-Q-Chem-t-rex.sh and answer the questions regarding the parameters of the geometrical progression of field amplitudes (see example below). The script will create multiple input files for the FF calculations based on the basic input file that you provided. After running the calculations, parse the output files using the script parse-t-rex.sh. Then run the T-REX program. Answer the the questions to compute the dipole moment and (hyper)polarizabilities (see example below).

Romberg differentiation is only available for methods where printing the total energy to high precision has been enabled. In the current version of Q-Chem, CC_PRINT_PREC is implemented for the following methods: HF, DFT, MP2, RI-MP2, MP3, CCD, CCSD, CCSD(T), QCISD, QCISD(T), TDDFT, and EOM-CCSD.

Note:  When using excited-state methods such as TDDFT, CIS, and EOM-CC, state ordering may switch when the external field is large.

With the T-REX program, you can compute the static dipole moment, polarizability, and first and second hyperpolarizabilities for one specific diagonal direction or for all Cartesian components, except that second hyperpolarizabilities are limited to iiii and iijj components. When computing all the components, you can obtain the norm of the dipole moment and polarizability, the Hyper-Rayleigh scattering first hyperpolarizability, the mean of the second hyperpolarizability, and other information

10.13.2.2 Step-by-step Example of a Finite-Field Calculation

Put the sample input file (given above) for a DFT calculation in the new directory; the name of the input file should be input. Run the input-Q-Chem-t-rex.sh script. The following questions are asked:

Components: x=1 y=2 z=3 all_beta=4 all_alpha=5
4
File name
eth-
Number of field amplitudes
5
Smallest field amplitude F_0
0.0004
a   (F_k=a^k F_0)
2.0
Number of files created         101

In this example, the answers correspond to FF calculations for F(k)=2.0k×0.0004 and for a geometric progression (k=0,1,2,3,4,5) of external field amplitudes. The calculation is set up for all the components.

After the FF calculations (i.e., after executing Q-Chem jobs for all generated inputs), parse the outputs using the parse-t-rex.sh script. The energies are written in a file called prelogfile.

Run the T-REX program and answer the questions:

Components: x=1 y=2 z=3 all_beta=4
4
Number of methods
1
Number of field amplitudes k_max+1
5
Smallest field amplitude   F_0
0.0004
Step-size a
2.0

The energies are ordered in a file called logfile and the results are printed in the results file.