DIM/QM: Discrete Interaction Model/Quantum Mechanics¶

DRF¶

The Discrete Solvent Reaction Field (DRF) model is a hybrid Quantum mechanical and Molecular Mechanics (QM/MM) model for studying solvation effects on (time-dependent) molecular properties such as dipole moments, excitation energies and (hyper)polarizabilities[141-145]. The classical solvent molecules are represented using distributed atomic charges and polarizabilities.

With the ADF-GUI it is simple to setup a DRF calculation. Another possibility is to use and/or modify a python script that is explained in the example DIMQM_DRF_plams. Both the ADF-GUI and the DRF python script have default settings which can be easily modified if this is required. In the simplest case the user just needs to provide the QM and DRF regions. By default, atomic charges that are used, in case of the ADF-GUI or the DRF python script, for the DRF are computed with LDA functional, DZP basis set and normal numerical quality. Also by default, atomic polarizabilities (computed with the Thole’s model) are taken from a inner database including H, C, N, O, F, S, Cl, Br, I atoms [146]. If one does not use the ADF-GUI or the DRF python script, the user has to specify this by hand.

Within the Discrete Solvent Reaction Field model the QM/MM operator is

where the first term, v^el , is the electrostatic operator and describes the Coulombic interaction between the QM system and the permanent charge distribution of the solvent molecules. The second term, v^pol , is the polarization operator and describes the many-body polarization of the solvent molecules, i.e. the change in the charge distribution of the solvent molecules due to interaction with the QM part and other solvent molecules. The charge distribution of the solvent is represented by atomic point charges and the many-body polarization by induced atomic dipoles at the solvent molecules. The induced atomic dipole at site s is found by solving a set of linear equations

where \(\alpha_{s,\alpha \beta}\) is a component of the atomic polarizability tensor at site s. The screened dipole interaction tensor is given by

where the damping functions f^T _st and f^E _st have been introduced, see also [146]. A smeared-out point charge model [147] is used for short-range damping of the QM/MM operator

The scaled distance, S_st , then replaces the normal distance, R_st , in the QM/MM operator.

In order to perform a DRF calculation two types of parameters (model atomic charges and atomic polarizabilities) for each type of atom in the MM part are required. The point charges should represent at least the permanent molecular dipole moment, and the distributed atomic polarizabilities the full molecular polarizability tensor. The atomic charges can straightforward be obtained using e.g. Multipole Derived Charges (MDC) [See section on MDC] and the distributed polarizabilities by adopting standard parameters or refitting them to match the calculated polarizability tensor [146,147]. This allows for a simple procedure to obtain the solvent model parameters which subsequently can be used in the DRF calculation.

Surface-enhanced response properties¶

In addition to the calculation of solvation effects on response properties, DIM/QM can also be used for simulating nonlinear optical properties of molecules near metal surfaces of a nanoparticle. A TDDFT description of the molecule is coupled to an electrodynamical treatment of the nanoparticle, including electrostatics and polarization of the atoms. In this way DIM/QM can be used to study surface-enhance Raman scattering (SERS) [362], surface-enhanced Raman optical activity (SEROA) [382], and plasmon-enhanced photochemistry.

Starting from ADF2018 the Polarizability Interaction Model (PIM) can be used to calculate damped first hyperpolarizabilities \(\beta\) and damped second hyperpolarizabilities \(\gamma\). Different types of \(\beta\)-tensors can be calculated, such as STATIC, OPTICALR, EOPE, and SHG. This allows for simulation of surface-enhanced hyper-Raman scattering (SEHRS) [487]. All the (sub)keywords for regular damped \(\beta\), such as BETA and QUADRATIC, can also be used. Also different types of \(\gamma\)-tensors can be calculated, such as STATIC, EFIOR, OKE, IDRI, EFISHG, THG and TPA. This allows for simulation of plasmon-enhanced two-photon absorption (PETPA). Except CUBIC, all the other (sub)keywords for regular damped \(\gamma\) can also be used.

Input options¶

DIMQM
  <DRF|PIM|CPIM>
  NOPOL
  NOCHAR
  NOCROSS
  DDA
  FULLGRID
  LOCALFIELD
  EFIELD x y z
  SCREEN <ERF|EXP|ESP|GAU|NONE> {length}
  :: FREQUENCY-DEPENDENT PARAMETERS
  FREQUENCY
  OM_H value
  OM_C value
  OM_O value
  OM_N value
  :: CONTROL OVER SOLVER
  ALGORITHM <BEST|DIRECT|BRUTE|SINGLE|MULTI>
  TOLERANCE tol
  NITER iterations
  VOLUME vol_in_nm^3
  MULTIPLIER aterm bterm cterm
  GRID <COARSE|MEDUIM|FINE>
  :: GRADIENT OPTIONS
  FORCEFIELD
  CHEMBOND qmindex dimindex
  COORDDEPEND
  CHEMISORPTION
  COORDPAR atomtype e0 e1 r0 r1 CNmax Rmax Rmin
  CHEMPAR atomtype e0 e1 r0 r1 cutoff
  PROJECTIONMATRIXPOINTS <ALL|CUTOFF radius|OFF>
  :: OUTPUT CONTROL
  PRINTATOMICDIPOLES
  PRINTLJPAR
  DEBUG
END

<DRF|CPIM|PIM>

DIM/QM relies on one of three descriptions of the system: Discrete Reaction Field (DRF), where the atoms interact via induced dipoles and static charge, Capacitance Polarizability Interaction Model (CPIM), where the atoms interact via induced dipoles and induced charges, and Polarizability Interaction Model (PIM), where the atoms interact via induced dipoles only. DRF is best for solvents, CPIM is best for small metal nanoparticles, and PIM is best for large metal nanoparticles. One and only one of these three keys must be included in every DIM/QM calculation.

NOPOL

The NOPOL key turns off the polarization terms, and thus all induced dipoles are zero. This key is only valid for DRF or CPIM calculations.

NOCHAR

The NOCHAR key turns off the charge terms, and thus all induced or static charges are zero. This key is only valid for DRF or CPIM calculations.

NOCHOSS

The NOCROSS key turns off the charge-dipole interactions. This key is only valid for CPIM calculations.

DDA

By default, the dipole-dipole, charge-dipole and charge-charge interactions are screened to take into account that atoms are not point charges. The DDA key will turn off this screening so that the results can be compared directly to the discrete dipole approximation (DDA).

FULLGRID

This is used in conjunction with the frozen density approximation.

LOCALFIELD

When the molecule interacts with a (for example) metal nanoparticle, there are two types of interactions: the image field and the local field. The image field is caused by the dipoles induced into the nanoparticle by the molecule’s electron density. This is always taken into account in a DIM/QM calculation. The local field arises by direct interactions of the nanoparticle with an external field. Addition of the LOCALFIELD key causes the DIM/QM calculation to include this effect, but by default this is not included in a DIM/QM calculation. In ADF2016 the correct behavior with local fields for the electric dipole-dipole, electric dipole-quadrupole, and electric dipole-magnetic dipole response tensors in the velocity gauge is implemented. This allows for correct calculation of surface-enhanced Raman optical activity (SEROA) (Ref. [382]) and plasmonic circular dichroism (Ref. [383]).

EFIELD x y z

The EFIELD key is used to include an external static electric field in the vector x y z. Internally, the static charges used in a DRF calculation use ADF’s EFIELD block, and therefore use of the EFIELD block is not allowed with the DIMQM block. This key is included so the user may include an electric field that would normally be included by the EFIELD block.

SCREEN <ERF|EXP|ESP|GAU|NONE> {length}

The SCREEN key indicates what functional form is used to screen the interactions between each DIM atom and the QM density. The choices are ERF (error function), EXP (exponential), ESP (error function for potential operator only), GAU (Gaussian), or NONE. For CPIM and PIM the default is GAU; for DRF, the default is ESP. In all cases, the default screening length is 1.0, but this may be changed with the optional length parameter.

FREQUENCY

The FREQUENCY key turns on frequency-dependent parameters.

OM_[HCON] value

The OM_H, OM_C, OM_O, and OM_N keys provide the resonance frequency (in atomic units) for the elements H, C, O and N, respectively. These keys are only for use with DRF and only when the older EXTERNALS block is used.

ALGORITHM <BEST|DIRECT|BRUTE|SINGLE|MULTI>

DIM/QM can choose between several solver algorithms. The DIRECT method solves the linear system of equations directly with a LAPACK routine; this should be considered the most robust method, but scales poorly with the number of atoms (O(N³ ) where N is the number of atoms in the system). The other three methods use an iterative technique, The BRUTE method (brute force) takes into account all atoms in the matrix-vector multiply step, and scales as O(N² ). The SINGLE method uses the single-level cell-multipole-method (CMM), wherein dipoles that are spatially similar are collected into a multipole moment which effectively reduces the system size. This also scales as O(N² ) but with a lower prefactor than BRUTE. The MULTI method uses the multi-level cell-multipole-method, which uses larger and large multipole the farther apart the dipoles are. This is the fastest method and scales as O(N log N).

Due to technical limitations, CPIM can only use DIRECT. Further, depending on the system size DIRECT or SINGLE may be more efficient than MULTI. To simplify choosing the solving algorithm, there is a BEST option that chooses the best algorithm for the particular system. BEST is the default option for algorithm.

TOLERANCE tol

The TOLERANCE key allows the user to specify a tolerance for the iterative solver. By default the tolerance is based on ADF’s INTEGRATION key. This has no effect with the DIRECT solver.

NITER iterations

The NITER key allows the user to specify the maximum number of iterations for the iterative solver. By default this is MAX(N/100, 200) where N is the number of DIM atoms.

VOLUME vol_in_nm3

The VOLUME key is used to specify the DIM system volume. The volume is used to determine how to partition the system for the cell multipole method (ALGORITHM options SINGLE or MULTI), and is also use to determine the scattering efficiencies for frequency-dependent polarizability calculations. The volume does not need to be supplied; if it is missing, it will be calculated based on each atom’s radius and the MULTIPLIER key.

MULTIPLIER aterm bterm cterm

An efficient way to get a approximation for the volume of the system is to sum the volume of each atom in the system, modified to account for the space between the atoms. This is done by modifying the atomic radius by a formula that takes into account the number of DIM atoms so that the effective radius changes with surface-to-bulk ratio. This formula is given by

r_eff = -ar/N^b + c

where r is the atomic radius, r_eff is the effective radius, N is the number of atoms, and the a, b, and c terms are the three parameters defined by the MULTIPLIER key. If the MULTIPLIER key is missing, the default values are 0.7, 0.5, and 1.13, respectively.

GRID <COARSE|MEDIUM|FINE>

In the cell multipole method (ALGORITHM options SINGLE or MUTLI), a certain number of the closest atom interactions must be calculated explicitly. The GRID key controls how many atoms must be calculated this way, with COARSE being the least and FINE being the most. COARSE will be the fastest to calculate but may be numerically unstable. FINE is slowest to calculate but is the most stable. If the GRID key is missing, the default is MEDIUM.

FORCEFIELD

The FORCEFIELD key indicates that the DIM/QM calculation will include the DIM/QM force field. Currently the only maintained potential is the Lennard-Jones 12-6 potential (see Ref. [362]). This key is required to perform a DIM/QM geometry optimization and vibrational frequencies. By default, the DIM/QM force field is not included into the calculation. DRF gradients are not implemented.

Currently, DIM/QM geometry optimizations must be done in Cartesian coordinates which is specified in the GEOMETRY block. The user should be aware that ADF’s default convergence criterion for a geometry optimization are relatively low, thus it is strongly suggested for a DIM/QM calculation to set the numerical integration quality (BeckeGrid) to good and change the convergence criterion of the max gradient to 1E-4. It is not possible to do geometry optimization with DRF.

COORDDEPEND

The COORDDEPEND key indicates that the DIM/QM force field will be coordination dependent. This only effects the Lennard-Jones parameters for DIM atoms (see Ref. [362]). By default, the DIM/QM force field is not coordination dependent.

CHEMISORPTION

The CHEMISORPTION key will include chemisorption corrections for all atoms that have chemisorption parameters within a given cutoff radius. By default, the DIM/QM force field does not included chemisorption corrections.

COORDPAR atomtype e0 e1 r0 r1 CNmax Rmax Rmin

The COORDPAR key allows the user to add additional coordination dependent parameter for a selected element type. atomtype specifies the element type (i.e., Ag for silver) for the given parameter set. e0, e1, r0, and r1 are the coordination dependent Lennard-Jones parameters; see Ref. [362]) for more details. The coordination numbers of the DIM atoms are computed as an effective coordination number. This scheme requires a maximum and minimum cutoff distances, Rmax and Rmin respectively, and a maximum coordination number, CNmax. All parameters need to be in atomic units.

CHEMPAR atomtype e0 e1 r0 r1 CUTOFF

The CHEMPAR key allows the user to add additional chemisorption dependent parameter for a selected element type. atomtype specifies the element type (i.e., N for nitrogen) for the given parameter set. e0, e1, r0, and r1 are the coordination dependent Lennard-Jones parameters; see Ref. [362]) for more details. The code determines if there is a chemical bond by a cutoff distance parameter, CUTOFF. If the QM-DIM bond is within the cutoff, the code uses the chemisorption parameter; otherwise, the code uses the standard parameter set. All parameters need to be in atomic units.

CHEMBOND qmindex dimindex

The CHEMBOND key indicates that there will a chemisorption correction used for the bond between the specified QM and DIM atoms. The user may repeat the CHEMBOD key up to 50 times to specify up to 50 different chemical bonds for the force field. The qmindex is an integer based on the order of atoms in the ATOMS block; i.e. the fifth QM atom in the ATOMS block would have qmindex = 5. The dimindex is the same but corresponds to the DIM atom involved in the bond. This key should only be used when the CHEMISORPTION key is also specified. When using CHEMBOND, the cutoff distance parameter for chemisorption correction parameter key will be ignored. It is suggested to use CHEMBOND if the user is generating a potential energy surface with a chemisorbed QM system.

PROJECTIONMATRIXPOINTS <ALL|CUTOFF radius|OFF>

The PROJECTIONMATRIXPOINTS key specifies what DIM atoms to include for the projection matrix when removing rigid motions out of the gradient. The methods available are ALL, CUTOFF, or OFF. The ALL option causes PROJECTIONMATRIXPOINTS to include all DIM atoms. OFF will turn off the removal of rigid motions. CUTOFF includes any DIM atom points within a cutoff radius from the center of mass of the QM system to the DIM atom points and requires a cutoff radius to be given in Angstrom. This key only applies to a geometry optimization. If the PROJECTIONMATRIXPOINTS key is not given, the option CUTOFF with a cutoff radius of 25.4 Angstrom is assumed.

PRINTATOMICDIPOLES

The PRINTATOMICDIPOLES key causes all the induced dipole moments of each DIM atom to be printed at the conclusion of each SCF cycle and each RESPONSE or AORESPONSE polarizability calculation. Because DIM/QM is typically used with many thousands of atoms, this can result in a large output file, but they may be useful for debugging purposes or to calculate electric fields. By default these are not printed.

PRINTLJPAR

The PRINTLJPAR key specifies that all Lennard-Jones parameters used for the calculation will be printed in the output file. The QM atoms’ Lennard-Jones parameters are also printed with the DEBUG key.

DEBUG

The DEBUG key will print out extra information in the process of the calculation.

EXTERNALS key for DRF¶

The EXTERNALS block key controls the input data for the MM atoms. The EXTERNALS block is designed for DRF calculations. For each MM atom the following data are required:

EXTERNALS
 atm num grp-nam grp-num, char, x, y, x, pol
 ...
 GROUP
 {...}
end

atm

Type of atom, i.e., H, O, ...

num

number of atoms (optional)

grp-nam

Name of the group to which the atom belongs

grp-num

Number of the group to which the atom belong

char

atomic charge (in atomic units)

x

x-coordinate

y

y-coordinate

z

z-coordinate

pol

atomic polarizability (in atomic units)

GROUP

Indicates the end of group

The separation of molecules into GROUP’s are important. Since in the many-body polarization operator only inter-molecular interactions, i.e. only interaction between sites which do not belong to the some group, are included. Therefore, it is important that the combined string (grp-nam + grp-num) is unique for each GROUP.

An example of a EXTERNALS block for two water molecules:

EXTERNALS
 O 4 water 2, -0.6690, -11.380487, -11.810553, -4.515226, 9.3005
 H 5 water 2,  0.3345, -13.104751, -11.837669, -3.969549, 0.0690
 H 6 water 2,  0.3345, -10.510898, -12.853311, -3.320199, 0.0690
 GROUP
 O 7 water 3, -0.6690,  -1.116350,   9.119186, -3.230948, 9.3005
 H 8 water 3,  0.3345,  -2.822714,   9.717033, -3.180632, 0.0690
 H 9 water 3,  0.3345,  -0.123788,  10.538199, -2.708607, 0.0690
 GROUP
 {...}
end

In this block, the parameters for the DIM atoms are defined.

DIMPAR
  Element
    RAD val
    POL val
    CAP val
    CHAR val
    OM val
    OM1 val
    OM2 val
    GM1 val
    GM2 val
    SIZE val
    BOUND val
    EXP /path/to/experimental/dielectric/file
    DRUDE plasma damping {EV}
    FERMI val
    <LRTZ|LRTZ1> osc res damp {EV}
    LRTZ2 osc pls res damp {EV}
    LRTZ3 pls res damp {EV}
  SUBEND
  XYZ
    {/absolute/path/to/coordinates.xyz}
    {natoms
     elem x.xxx y.yyy z.zzz
     elem x.xxx y.yyy z.zzz
     elem x.xxx y.yyy z.zzz
     ...}
  SUBEND
  GROUP groupname
    natoms
     elem x.xxx y.yyy z.zzz {pol=pol} (cap=cap} {char=char} {..}
     elem x.xxx y.yyy z.zzz {pol=pol} (cap=cap} {char=char} {..}
     elem x.xxx y.yyy z.zzz {pol=pol} (cap=cap} {char=char} {..}
     ...
  SUBEND
END

Element

Within the DIMPAR block, you will need a sub-block that defines the parameters for each element that is in your DIM system. You will need to replace ‘Element’ with the element you are assigning parameters to, as in:

Ag
 ...
SUBEND

if you are assigning parameters to Ag. Note that the first letter MUST be capitalized and the second MUST be lowercase.

RAD val

RAD specifies the atomic radius in the unit defined by the input file. RAD is required for all PIM calculations, all calculations with ALOGORITHM options SINGLE or DIRECT, and all frequency-dependent calculations where the AORESPONSE LIFETIME key is given.

POL val

POL specifies the polarizability parameter (in a.u.) used in DRF or CPIM.

CAP val

CAP specifies the capacitance parameter (in a.u.) used in CPIM.

CHAR val

CHAR specifies the atomic charge (in a.u) used in DRF.

OM val

OM specifies the resonance frequency (in a.u) used in DRF. This replaces the OM_[HCON] key in the DIMQM block.

OM1 val

OM1 specifies the \(\omega\)₁ parameter (in a.u) used in CPIM.

OM2 val

OM2 specifies the \(\omega\)₂ parameter (in a.u) used in CPIM.

GM1 val

GM1 specifies the \(\gamma\)₁ parameter (in a.u) used in CPIM.

GM2 val

GM2 specifies the \(\gamma\)₂ parameter (in a.u) used in CPIM.

SIZE val

SIZE specifies the size-dependent parameter used in CPIM.

EXP /absolute/path/to/experimental/dielectric/file

In PIM, the atomic polarizabilities are calculated from the dielectric constant. If you have access to the experimental dielectric constant, this may be supplied directly to DIM/QM. The values will be splined, so it is not necessary to ensure that each frequency at which you are calculating be in the file. DIM/QM expects the file to be formatted with the wavelength (in nm) in the first column, the real part of the dielectric in the second column, and the imaginary part of the dielectric in the third column. All other columns that may exist will be ignored, as well as lines beginning with the hash (#) symbol.

BOUND val

A Drude function is typically written as

\[\epsilon_\infty - \frac{\omega_p^2}{\omega (\omega + i \gamma)}\]

with the second term being the Drude function, and the first term accounting for bound electrons. For a conductor with no bound electrons, \(\epsilon\)_∞ = 1 which is the default value for BOUND. To account for bound electrons you may set BOUND to a value greater than 1. This key only affects PIM.

DRUDE plasma damping {EV}

The formula for a Drude function is

\[\epsilon_\infty - \frac{\omega_p^2}{\omega (\omega + i \gamma)}\]

where \(\epsilon\)_∞ represents the bound electrons (as discussed for BOUND), \(\omega\)_p (plasma) is the plasma frequency, and \(\gamma\) (damping) is the damping parameter (decay rate). Optionally, EV may be added to specify the values be read in units of electron volts, otherwise they are read in units of a.u. This key only affects PIM.

FERMI val

The FERMI key is used to specify a Fermi velocity (in m/s) so that the Drude function may be size-corrected using a modified Drude function:

\[\epsilon_\infty - \frac{\omega_p^2}{\omega (\omega + i (\gamma + v_\text{fermi}/R_\text{eff}))}\]

where +v_fermi is the Fermi velocity and R_eff is the effective nanoparticle radius. This can also be used in conjunction with EXP and DRUDE to size-correct experimental dielectric parameters. This key only affects PIM:

<LRTZ|LRTZ1> osc res damp {EV}
LRTZ2 osc pls res damp {EV}
LRTZ3 pls res damp {EV}

There are three forms of the Lorentzian function seen in the literature:

\[\begin{split}& \sum_n \frac{f_n \Omega_{0,n}^2}{\Omega_{0,n}^2 - \omega^2 - i \Gamma_n \omega} \\ & \sum_n \frac{f_n \omega_{p}^2}{\Omega_{0,n}^2 - \omega^2 - i \Gamma_n \omega} \\ & \sum_n \frac{\Omega_{p,n}^2}{\Omega_{0,n}^2 - \omega^2 - i \Gamma_n \omega}\end{split}\]

where \(\Omega\)_0,n (res) is a bound electron resonance frequency, f_n (osc) is a bound electron oscillator strength, \(\Gamma\)_n (damp) is a bound electron excited state decay rate (or damping parameter), \(\omega\)_p (pls) is the free electron plasma frequency, and \(\Omega\)_p,n (pls) is the bound electron plasma frequency. You may choose the Lorentzian for against which your parameters were parametrized. The top form is LRTZ1, the middle is LRTZ2, and the bottom is LRTZ3. Because LRTZ1 is the most common, it is also aliased as LRTZ. Optionally, EV may be added to specify the values be read in units of electron volts, otherwise they are read in units of a.u. This key only affects PIM. You may give any of the form of the LRTZ key up to 50 times supply up to 50 Lortenzian functions to make a Drude-Lorentz function.

XYZ

The XYZ sub-block is where the DIM atom coordinates are given. Two methods of supplying coordinates are allowed.

In-File Coordinates

As an example of how to supply coordinates in-file, imagine you wish to calculate a Au dimer system on the Z-axis. You might define your coordinates as:

XYZ
 2
 Au 0.0 0.0 0.0
 Au 0.0 0.0 3.0
SUBEND

The first line gives the number of atoms to follow. Every line after that contains the element in the first column (first letter MUST be capitalized, second MUST be lowercase), then the x-component, then the y-component, then the z-component. You may not number the atoms.

External File Coordinates

When the DIM system size becomes large, it is often more convenient to keep the DIM coordinates in a separate file. The XYZ block would then look like:

XYZ
 /absolute/path/to/coordinates.xyz
SUBEND

Note that you MUST include the absolute path to your file.

The .xyz file is set up identically to the in-file table, except that there is a space between the number of atoms and the first coordinate in case a comment need be added. The .xyz file for our dimer system would be:

2
A gold dimer (this line will be ignored)
Au 0.0 0.0 0.0
Au 0.0 0.0 3.0

GROUP groupname

The GROUP sub-block is where the DIM atom coordinates and parameters are given. Combines the XYZ subblock and Element subblock. A (unique) groupname is required (maximum 10 characters).

Example for a water molecule:

GROUP water1
3
O  0.00000  0.00000  0.59372  pol=5.7494  char=-0.6690
H  0.00000  0.76544 -0.00836  pol=2.7927  char=0.3345
H  0.00000 -0.76544 -0.00836  pol=2.7927  char=0.3345
SUBEND

The first line gives the number of atoms to follow. Every line after that contains the element in the first column (first letter MUST be capitalized, second MUST be lowercase), then the x-component, then the y-component, then the z-component. You may not number the atoms. At the end of the line for each element you can specify the DIM/QM parameters, as in the format in the example. The DIM/QM parameters can be ‘POL’, ‘CAP’, ‘CHAR’, ‘OM1’, ‘OM2’, ‘GM1’, ‘GM2’, ‘SIZE’, ‘RAD’, ‘SPILLOUT’, ‘BOUND’, ‘DRUDE’, ‘FERMI’, ‘LRTZ’, or ‘EXP’.