TURBOMOLE manual
v8.0

21.5 Polarizable embedding calculations

Realizable embedding (PE) calculations are a based on a hybrid model of quantum mechanics and molecular mechanics (QM/MM) in which the classical region is represented by an electrostatic potential with up to octupole moments and induced point dipole moments. The main improvement over the more common QM/MM approaches without polarizable MM sites can be found for the description of electronic excitations but also for any other process which causes a significant change in the QM density and which is accompanied by a fast response of the environment.

In TURBOMOLE, ground state energies computed with the dscf, ridft, and ricc2 module and electronic excitation properties based on RI-CC2 and RI-ADC(2) are implemented. The excited state analytic gradients are also available at the RI-ADC(2) level. The general theory is presented in ref. [380] and [285, 286], the PERI-CC2 model and the TURBOMOLE implementation is described in ref. [283].

21.5.1 Theory

In the following, only the most important ideas are presented and discussed with a focus on the PERI-CC2 model. The essential concept is the introduction of an environment coupling operator \(\hat G(\mathbf{D}^{\text{CC}})\)

\[\begin{align} \hat G(\mathbf{D}^{\text{CC}}) & = \hat G^{\text{es}} + \hat G^{\text{pol}}(\mathbf{D}^{\text{CC}}) \end{align}\](21.14)

with the electrostatic contribution

\[\begin{align} \hat G^{\text{es}} &= \sum_{m = 1}^M \sum_{k = 0}^K \sum_{pq} \Theta^{(k)}_{m,pq} \mathbf{Q}_m^{(k)} \hat{E}_{pq} \end{align}\](21.15)

and the polarization contribution

\[\begin{align} \hat G^{\text{pol}}(\mathbf{D}^{\text{CC}}) &= \sum_{u = 1}^U \sum_{pq} \Theta^{(1)}_{u,pq} \mathbf{\mu}_u^{\text{ind}}(\mathbf{D}^{\text{CC}}) \hat{E}_{pq} \, . \end{align}\](21.16)

Here, \(\Theta^{(k)}_{m,pq}\) are multipole interaction integrals of order \(k\) and \(\mathbf{\mu}_u^{\text{ind}}\) are the induced dipoles which can be obtained from the electric field \(\mathbf{F}_u\) and the polarizability \(\mbox{\boldmath$\alpha$}_u\) at a site \(u\):

\[\begin{equation} \mathbf{\mu}_u^{\text{ind}} = \mathbf{F}_u \mbox{\boldmath$\alpha$}_u \end{equation}\](21.17)

Because the induced dipoles depend on the electron density and vice versa, their computations enter the self-consistent part of the HF cycle. Introducing \(\hat G(\mathbf{D}^{\text{CC}})\) into standard equations for the HF reference state and the CC2 equations leads to a general PE-CC2 formulation. To maintain efficiency, a further approximation has been introduced which makes the operator only dependent on a CCS-like density term. These general ideas define the PERI-CC2 model and allow to formulate the corresponding Lagrangian expression

\[\begin{align} L_{\text{PERI-CC2}}(\mathbf{t},\bar{\mathbf{t}}) &= E_{\text{PE-HF}} + \ensuremath{\langle \text{HF}|} \hat{W}(\hat{T}_1 + \hat{T}_2 + \frac{1}{2} \hat{T}_1^2) \ensuremath{|\text{HF}\rangle} + \\ & \sum_{\mu_1} \bar{t}_{\mu_1} \ensuremath{\langle \mu_1|} \tilde{W} + [ \hat{F}^{\text{PE}}, \hat{T}_1 ] + [\tilde{W}, \hat{T}_2] \ensuremath{|\text{HF}\rangle} + \\ & \sum_{\mu_2} \bar{t}_{\mu_2} \ensuremath{\langle \mu_2|} \tilde{W} + [ \hat{F}^{\text{PE}}, \hat{T}_2 ] \ensuremath{|\text{HF}\rangle} \\ & - \frac{1}{2} \sum_{uv} F_u^{\text{elec}}(\mathbf{D}^{\Delta\prime}) R_{uv} F_v^{\text{elec}}(\mathbf{D}^{\Delta\prime}) \end{align}\](21.18)

from which all PERI-CC2 equations including the linear response terms may be derived. Note that the dependency on the density couples the CC amplitude and multiplier equations for the ground state solution vector.
This coupling is avoided by the simplified polarizable embedding method (sPE) described in ref. [287]. \[\begin{equation} \begin{split} L_{\text{sPERI-CC2}}(\mathbf{t},\bar{\mathbf{t}}) =& E^{\text{PE-HF}} \\ &+ \langle\Lambda|\hat{g}_{\text{N}} + \hat{F}^{\text{PE}}_{\text{N}} + \hat{G}_{\text{N}}^{\text{pol}}({\bf D}^{\Delta\text{CC}})|\text{CC}\rangle \; . \end{split} \end{equation}\](21.19) The subscript N indicates that the operator is normal ordered with respect to the Hartree–Fock state. Here, a polarization operator \(\hat{G}^{\text{pol}}({\bf D})\) was introduced, \[\begin{equation} \hat{G}^{\text{pol}}({\bf D}^{\Delta\text{CC}}) = -\frac{1}{2} \sum_{uv}^M \left(\hat{{\bf F}}_u\right)^\top {\bf R}_{uv} {\bf F}^{\text{el}}_v({\bf D}^{\Delta\text{CC}}) \; . \end{equation}\](21.20) which depends on the elements of the difference density matrix \({\bf D}^{\Delta\text{CC}}\). These are defined as \[\begin{equation} D_{pq}^{\Delta\text{CC}}=\langle\text{HF}| \hat{a}_p^\dagger \hat{a}^{\phantom{\dagger}}_q|\text{CC}\rangle - D_{pq}^{\text{HF}} \; , \end{equation}\](21.21) hence, they do not depend on the Lagrangian multipliers.

21.5.2 Computational details: SCF calculations

To carry out a PE-SCF calculation with the DSCF or RIDFT module, you have to specify the following in the control file:

$point_charges pe [options]
<length unit>
<no. MM sites> <order k> <order pol> <length exclude list>
<list of MM sites: exclude list, xyz coords, multipole mom., pol. tensor>
length unit

specifies the unit for the MM site coordinates (use AA or AU)

no. MM sites

the amount of MM sites (length of the list)

order k

the order of multipoles used (0: point charges, 1: dipole moments, 2: quadrupole moments, 3: octupole moments)

order pol

the treatment of polarizabilities (0: none, 1: isotropic, 2: anisotropic)

length exclude list

number of elements in the exclude list

list of MM sites

each MM sites is described on one line, entries separated by blanks; first entry is the exclude list of with as much elements as defined in the head line (If the first element in the exclusion list of one site occurs in the exclude list of another site, they do not contribute to each others polarization); next follows the MM site coordinates in (x,y,z positions), the point charge, the dipole moment (for \(k \ge 1\), x,y,z component), the quadrupole moment (for \(k \ge 2\), xx, xy, xz, yy, yz, zz component), the octupole moment (for \(k = 3\), xxx, xxy, xxz, xyy, xyz, xzz, yyy, yyz, yzz, zzz component), the polarizability ( one component for pol-order 1, xx, xy, xz, yy, yz, zz component for pol-order 2)

An example for a polarizable embedding with coordinates given in Å, point charges and isotropic polarizabilities:

$point_charges pe
AA
6 0  1  1
  39   -0.2765102481    2.5745845304    3.5776314866    0.038060  15.217717
  39    1.3215071687    2.3519378014    2.8130403183   -0.009525  14.094642
  39   -0.5595582934    1.2645007691    4.7571719292   -0.009509  14.096775
  39   -1.5471918244    2.5316479230    2.3240961995   -0.009519  14.096312
  39   -0.3207417883    4.1501938400    4.4162313889   -0.009507  14.096476
  41   -1.1080691595    4.9228723099   -1.6753825535    0.038060  15.217717
  41   -0.9775910525    6.5274614891   -2.4474576239   -0.009525  14.094642
  41   -2.5360480539    4.8923046027   -0.6040781123   -0.009509  14.096775
  41    0.3630448878    4.6028736791   -0.7155647205   -0.009519  14.096312
  41   -1.2817317422    3.6689143712   -2.9344225518   -0.009507  14.096476

All values are given in atomic units (except coordinates if stated otherwise). These data are mandatory. An alternative input format can also be used by specifying daltoninp as option on the $point_charges line. The format is completely compatible with the current Dalton 2015 input format (The definition can be found in the manual of the Dalton program at https://daltonprogram.org/documentation/. See there for more information).

In addition, you can specify further options on the same line as the $point_charges flag. These are:

  • rmin=<float>: minimum distance between an active MM site and any QM center (in a.u.), treatment is handle by option iskip, (DEFAULT: 0.00 a.u.)

  • iskip=(1,2) : treatment of too close MM sites

    • (1) zeroing all contributions

    • (2) distribute values to nearest non-skipped MM site (DEFAULT)

  • rmax=<float>: maximum distance between an active MM site and QM center of coordinates (in a.u.), sites too far away are skipped (zeroed) (DEFAULT: 1000.00 a.u.)

  • nomb: no treatment of many body effects between induced dipoles (all interaction tensors on the off-diagonal of the response matrix are set to Zero); works best with isotropic polarizabilities, speeds up calculations (especially for large response matrices), has reduced accuracy, not well tested so far

  • longprint=(1,2,3): sets a flag for additional output

    • (1) print all MM site input information

    • (2) additionally: print all induced dipoles due to nulcei/multipole/electron electric filed

    • (3) additionally: print response matrix

  • file=<input file>: specifies a file from which the data group $point_charges is read. Note that all options which are following on the line in the control file are then ignored because reading continues in the input file (But here, further options can be specified after the $point_charges flag). The file has to start with $point_charges as top line and should be finished with $end

  • ccdens: activates the polarizable embedding method with an approximate non-hermitian density described in ref. [287] (only implemented for excitation energies and one-photon transition moments).

Limitations with respect to standard SCF computations:

  • In PE-SCF computations, symmetry cannot be exploited.

  • PE-SCF computations do not work in parallel (MPI parallelization).

  • For two-component all-electron calculations, the decoupling of the embedding potential is neglected, i.e. it is affected by a picture-change error just like the two-electron integrals.

The energy of a PE-SCF calculation printed in the output contains the following terms:

\[\begin{equation} E_{\text{PE-SCF}} = E_{\text{QM}} + E_{\text{QM/MM,es}} + E_{\text{pol}} \end{equation}\](21.22)

Here, \(E_{\text{QM}}\) is the energy of the quantum mechanical method of your choice, \(E_{\text{QM/MM,es}}\) the electrostatic interaction energy between the QM and the MM region, and \(E_{\text{pol}}\) the energy gain due to the total of induced dipole moments. If necessary, missing terms can be computed without knowledge of the electron distribution.

At the moment, TURBOMOLE does not offer the possibility to generate the necessary potentials or to create a potential file from a set of coordinates. Embedding potentials can be obtained from literature or generated by approaches like the LoProp method.[381] Atom centered polarizabilities are also available from other methods or from experiment. Finally, there are some polarizable force fields which, in principle, can be used for the PE method (for example, the AMOEBA force field).

21.5.3 Computational details for post-SCF methods

PERI-CC2 calculations:

Apart from the definition of the embedding described above, the input for PERI-CC2 calculations is the same as without polarizable embedding.

There are several limitations for the use of PERI-CC2:

  • only ground state energies, excitation energies and transition moments are supported (no other properties or gradients and so on)

  • no use of symmetry

  • no MPI parallelization is available (but SMP binaries work)

  • open-shell systems are not covered (exception: two-component references)

PE-MP2 within PTED Reaction-Field Scheme:

In addition to the PERI-CC2 method, the calculation of the ground-state energy is available at PE-MP2 level within the framework of PTE and PTED reaction-field schemes. The implementation of PE-MP2 using the PTED approach is inspired by the work published by Lunkenheimer and Köhn.[281] For the details about the PTED reaction-field scheme read section 21.2.5. For the calculation of the ground-state energy the following data groups must be included in the control in addition to the PE data groups given in section 21.5.2 and the $ricc2 data group:

$ricc2
  mp2
$response
  fop relaxed
$reaction_field
  scrf state=(x)
  PTED
  cycle

The cycle flag in the $reaction_field data group controls the PTED macro-iterations between the SCF and post-SCF calculations. In order to invoke the PTED-PE-MP2 calculations, the pecc2 script must be used. The combination of PTED-PE-MP2 with SCS and SOC is also available. Furthermore, the calculation of ground-state energy is possible for closed-shell and open-shell systems with the SMP (OpenMP) and MPI parallelizations.

PE-ADC(2) within post-SCF Reaction-Field Scheme

:The new PE-ADC(2) method[382] implemented in the framework of the post-SCF reaction field scheme makes the following calculations available in the ricc2 module for both closed-shell and open-shell systems:

  • Vertical excitation energy and transition moments

  • Excited-state analytic gradients

  • Excited-state properties

For these calculations, in addition to the typical data groups $ricc2, $excitations and $point_charges, the following keywords should be added to the control file:

 $reaction_field
   post-SCF
   ccs-like