12 Molecules in Complex Environments: Solvent Models, QM/MM and QM/EFP Features, Density Embedding

12.8 Polarizable Embedding Model

The polarizable embedding (PE) model is a fragment-based quantum-classical explicit embedding scheme to model molecular properties in complex heterogeneous environments 688, 689. The theory is explained thorougly in literature 688, 689, 578. In essence, the environment is represented by a multi-center multipole expansion to model electrostatic interactions, whereas polarization is taken into account by dipole-dipole polarizabilities placed at the expansion points. Polarization effects can thus be treated fully self-consistently by mutual polarization of the environment and the quantum region.

A recent tutorial review on how to prepare PE calculations in general (creating embedding potentials) is also available 879. For automated generation of embedding potentials, please refer to the PyFraME tool 11 1 https://gitlab.com/FraME-projects/PyFraME which is also explained in the aforementioned review.

PE can be used for Hartree-Fock and density-functional theory ground-state SCF methods. In addition, PE has been combined with the algebraic-diagrammatic construction for the polarization propagator (ADC) 816, explained in the subsequent section.

12.8.0.1 PE-ADC

The combined scheme of the PE model and ADC (PE-ADC) 816 is built on top of a PE-HF ground-state calculation and takes into account perturbative corrections of the excitation energies in a density-driven manner. That is, after the Hartree-Fock ground-state calculations, the induced dipole moments in the environment are kept frozen and an ADC calculation is performed as usual. Thereafter, perturbative corrections of the electronic excitation energies are calculated based on i) the transition density (perturbative linear-response-type correction, ptLR), and ii) the difference density (perturbative state-specific correction, ptSS) for each excited state.

12.8.0.2 PE Job Control

The PE job control is accomplished in two sections, $rem and $pe. To enabling PE-ADC, specification of the ADC method and other ADC job control parameters (thresholds, max. iterations etc.) should be set in the $rem section. PE-ADC supports the excited state analysis (STATE_ANALYSIS) carried out by the libwfa module.

PE
       Turns PE on.
TYPE:
       BOOLEAN
DEFAULT:
       False
OPTIONS:
       True Perform a PE calculation. False Don’t perform a PE calculation.
RECOMMENDATION:
       Set the $rem variable PE to TRUE to start a PE calculation.

Note:  Turning PE on disables symmetry by setting SYM_IGNORE to TRUE.

Note:  Setting the REM variables USE_LIBQINTS and GEN_SCFMAN to TRUE is required to run PE.

The PE-specific options can be set in the $pe input section. The format of the $pe section requires key and value pairs separated by a space character:

$pe
   <keyword>  <parameter>
$end

Note:  The following job control variables belong only in the $pe section. Do not place them in the $rem section.

POTFILE
       Path of the potential file.
INPUT SECTION: $pe
TYPE:
       STRING
DEFAULT:
       potfile.pot
OPTIONS:
       Provide the path/name of the potential file.
RECOMMENDATION:
       None

DIIS
       Use DIIS acceleration to obtain induced moments.
INPUT SECTION: $pe
TYPE:
       BOOLEAN
DEFAULT:
       TRUE
OPTIONS:
       TRUE Turns DIIS acceleration on. FALSE Turns DIIS acceleration off (normal Jacobi solver is used).
RECOMMENDATION:
       TRUE

CONVERGENCE_INDUCED
       Threshold for induced moments convergence. Converge induced moments to a residual norm of 10-CONVERGENCE_INDUCED.
INPUT SECTION: $pe
TYPE:
       INTEGER
DEFAULT:
       8 Corresponding to 10-8
OPTIONS:
       n12 Corresponding to 10-n
RECOMMENDATION:
       Use the default unless higher accuracy is desired.

MAXITER
       Maximum number of iterations for induced moments.
INPUT SECTION: $pe
TYPE:
       INTEGER
DEFAULT:
       50
OPTIONS:
       n1
RECOMMENDATION:
       Use the default. If more iterations are required to converge the induced moments, there might be an error in the system setup.

BORDER
       Activate border redistribution/removal options for sites in proximity to the QM/MM border.
INPUT SECTION: $pe
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       TRUE Enable border options. FALSE Disable border options.
RECOMMENDATION:
       None

BORDER_TYPE
       Remove or redistribute multipole moments/polarizabilities.
INPUT SECTION: $pe
TYPE:
       STRING
DEFAULT:
       REMOVE
OPTIONS:
       REMOVE remove multipole moments/polarizabilities. REDIST redistribute multipole moments/polarizabilities.
RECOMMENDATION:
       None

BORDER_RMIN
       Minimum distance from QM atoms to MM sites to be taken into account for removal/redistribution
INPUT SECTION: $pe
TYPE:
       FLOAT
DEFAULT:
       2.2 (AU)
OPTIONS:
       r>0 (Unit depends on BORDER_RMIN_UNIT)
RECOMMENDATION:
       None

BORDER_RMIN_UNIT
       Unit of BORDER_RMIN, default is atomic units (AU)
INPUT SECTION: $pe
TYPE:
       STRING
DEFAULT:
       AU
OPTIONS:
       AU Use atomic units. AA Use Angstrom.
RECOMMENDATION:
       None

BORDER_REDIST_ORDER
       Order from which on moments are removed. For example, if set to 1 (default), only charges are redistributed and all higher order moments are removed.
INPUT SECTION: $pe
TYPE:
       INTEGER
DEFAULT:
       1
OPTIONS:
       n=0,1,2,3,
RECOMMENDATION:
       None

BORDER_N_REDIST
       Number of neighbor sites to redistribute multipole moments/polarizabilities to. The default (-1) redistributes to all sites which are not in the border region.
INPUT SECTION: $pe
TYPE:
       INTEGER
DEFAULT:
       -1
OPTIONS:
       n=-1,1,2,3,,number of MM sites
RECOMMENDATION:
       Use the default value.

BORDER_REDIST_POL
       Redistribute polarizabilities? If set to FALSE, polarizabilities are removed.
INPUT SECTION: $pe
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       TRUE Redistribute polarizabilities. FALSE Remove polarizabilities.
RECOMMENDATION:
       None

Example 12.30  Input for a PE-HF calculation of 4-Nitroaniline in presence of six water molecules

$molecule
0 1
C          8.64800        1.07500       -1.71100
C          9.48200        0.43000       -0.80800
C          9.39600        0.75000        0.53800
C          8.48200        1.71200        0.99500
C          7.65300        2.34500        0.05500
C          7.73200        2.03100       -1.29200
H         10.18300       -0.30900       -1.16400
H         10.04400        0.25200        1.24700
H          6.94200        3.08900        0.38900
H          7.09700        2.51500       -2.01800
N          8.40100        2.02500        2.32500
N          8.73400        0.74100       -3.12900
O          7.98000        1.33100       -3.90100
O          9.55600       -0.11000       -3.46600
H          7.74900        2.71100        2.65200
H          8.99100        1.57500        2.99500
$end

$pe
 potfile  gen_scfman_pe_potfile.pot
$end

$comment
The corresponding potential file \texttt{gen\_scfman\_pe\_potfile.pot} can be
found in the samples folder.
$end

$rem
   METHOD               HF
   BASIS                STO-3G
   PE                   TRUE
   SYM_IGNORE           TRUE
   USE_LIBQINTS         TRUE
$end

12.8.0.3 PE output

After SCF convergence, the PE module prints a summary of PE energy contributions, for example:

----------------------------------------------------------------------
Polarizable Embedding Summary:

   Electrostatics:
      Electronic:                   0.30901227399981
      Nuclear:                      -0.32134940137969
      Multipole:                    0.00000000000000
      Total:                        -0.01233712737988

   Polarization:
      Electronic:                   -0.01817189734325
      Nuclear:                      0.01717961961137
      Multipole:                    -0.02091890381649
      Total:                        -0.02191118154837

   Total Energy:                    -0.03424830892825
----------------------------------------------------------------------

If a PE-ADC calculation is carried out, the perturbative corrections are printed together with the excitation energies:

  Excited state 1 (singlet, A)                                     [converged]
  ----------------------------------------------------------------------------
  Term symbol:  2 (1) A                                     R^2 =  1.84764e-13

  Total energy:                                           -483.3704138865 a.u.
  Excitation energy:                                               3.906651 eV
  -------------------------------------------
  PE ptSS energy correction:                                      -0.001804 eV
  Corrected Excitation Energy (ptSS):                              3.904847 eV
  -------------------------------------------
  PE ptLR energy correction:                                      -0.000096 eV
  Corrected Excitation Energy (ptLR):                              3.906554 eV
  -------------------------------------------