DIM/QM: Discrete Interaction Model/Quantum Mechanics
****************************************************

The Discrete Interaction Model/Quantum Mechanics (DIM/QM) method
facilitates the calculation of the optical properties of molecules
under the influence of a discrete solvent or a metal nanoparticle;
see, for example, Ref. [1]. DIM/QM relies on one of three descriptions
of the system: Discrete Reaction Field (DRF), where the atoms interact
via induced dipoles and static charges, 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.

The DRF module is now part of DIM/QM, so DRF is a submodule of DIM/QM.

To perform a DIM/QM calculation, two block keys are required. The
first is the DIMQM block, which activates the DIM/QM module. The
parameters for each DIM atom must also be given, either with the
DIMPAR block, which is most convenient for metal nanoparticles, or
with the EXTERNALS block, which is designed for DRF.

In ADF2019, the bond energies calculated with DIM/QM in ADF have been
corrected.


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 [7]
[8] [9] [10] [11]. The classical solvent molecules are represented
using distributed atomic charges and polarizabilities.

With the ADF-GUI, it is simple to set up 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 that 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, the atomic charges used (in the ADF-GUI
or the DRF Python script) for the DRF are computed with the LDA
functional, DZP basis set, and normal numerical quality. Also by
default, atomic polarizabilities (computed with the Thole model) are
taken from an internal database including H, C, N, O, F, S, Cl, Br,
and I atoms [2]. If one does not use the ADF-GUI or the DRF Python
script, one has to specify this by hand.

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

   H_\text{QM/MM} = \sum_i v_\text{DRF} (r_i, \omega) = \sum_i
   V^\text{el} (r_i) + V^\text{pol} (r_i, \omega)

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

   \mu^\text{ind}_{s,\alpha} (\omega) = \alpha_{s,\alpha \beta} \left(
   F^\text{init}_{s,\beta} (\omega) + \sum_{t \neq s}
   T^{(2)}_{st,\beta \gamma} \mu^\text{ind}_{t,\gamma} (\omega)
   \right)

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

   T^{(2)}_{st,\beta \gamma} = 3f^T_{st} R_{st,\alpha} R_{st,\beta} /
   R^5_{st} - f^E_{st} \delta_{\alpha \beta} / R^3_{st}

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

   1/R_{st}  \rightarrow 1/S_{st}  = erf(R_{st})/R_{st}

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 atom type in the
MM part are required. The point charges should represent at least the
permanent molecular dipole moment, and the distributed atomic
polarizabilities should represent the full molecular polarizability
tensor. The atomic charges can be straightforwardly 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 [2] [3].
This allows for a simple procedure to obtain solvent model parameters
that can subsequently be used in the DRF calculation.

DRF gradients can be calculated by using the DRF2 keyword in the DIMQM
block. This is an expert option and should be used with great caution.
The DRF2 keyword uses DRF, but with Gaussian screening like that used
in the PIM method. It also allows input of parameters by atom type,
like the PIM method. The user must specify a radius and charge for the
atom type, as shown below.

   DIMPAR
     O
       rad=0.9028
       char=-0.6690
       STATIC
     SUBEND
     H
       rad=0.2479
       char=0.3345
       STATIC
     SUBEND
     Ag
       rad=1.445
       char=0.0000
     SUBEND
   END

Here, the STATIC keyword is used to tell DIM/QM to ignore the
frequency-dependent response for those atom types.


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-enhanced Raman scattering (SERS) [1], surface-enhanced Raman
optical activity (SEROA) [4], 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) [5]. All (sub)keywords for regular damped \beta, such as BETA
and QUADRATIC, can also be used. Different types of \gamma tensors can
also 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 <Resp|All|Static>
     NOIMAGE
     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|MEDIUM|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
     LOCALDIM
     LOCALDIRECT
     DIMBUFFER
     QMBUFFER
     CUTOFFDIST
   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 charges, 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. However, in ADF2019,
   there is also the expert option DRF2, which allows for DRF
   gradients but should be used cautiously.

"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.

"NOCROSS"
   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 account for the fact 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 <Resp|All|Static>"
   Default is All, when no option is specified. When the molecule
   interacts with, for example, a 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 from 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 was implemented. This allows for correct calculation of
   surface-enhanced Raman optical activity (SEROA) (Ref. [4]) and
   plasmonic circular dichroism (Ref. [6]).

   The default option should be used in most situations, as the other
   options are expert options. In an effort to speed up calculation
   times, one can request calculation of the local fields only during
   the response section of the calculation by specifying "LOCALFIELD
   Resp". All will calculate the local field for both ground-state
   calculations and response calculations, while Static will calculate
   the local field during the ground state only.

"NOIMAGE"
   Expert option. By default, the image field is always included in a
   DIM/QM calculation. With this keyword one can turn it off.

"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 AMS's System%ElectrostaticEmbedding, and therefore
   the use of the AMS System%ElectrostaticEmbedding is not allowed in
   combination with the DIMQM block. This key is included so the user
   may include an external static electric field.

"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. 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 it
   scales poorly with the number of atoms (O(N^3 ) where N is the
   number of atoms in the system). The other three methods use
   iterative techniques. The BRUTE method (brute force) takes all
   atoms into account in the matrix-vector multiplication step and
   scales as O(N^2 ). The SINGLE method uses the single-level cell-
   multipole method (CMM), in which dipoles that are spatially similar
   are collected into a multipole moment, which effectively reduces
   the system size. This also scales as O(N^2 ), but with a lower
   prefactor than BRUTE. The MULTI method uses the multi-level cell-
   multipole method, which uses larger and larger multipoles 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
   system at hand. BEST is the default option for the 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" setting. 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
   used 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 an 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 the 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 MULTI), a
   certain number of the closest atom interactions must be calculated
   explicitly. The GRID key controls how many atom-atom interactions
   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 the 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. [1]). This
   key is required to perform a DIM/QM geometry optimization and
   vibrational frequencies. By default, the DIM/QM force field is not
   included in 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 criteria for a
   geometry optimization are relatively loose; thus, for a DIM/QM
   calculation it is strongly recommended 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 affects the Lennard-Jones
   parameters for DIM atoms (see Ref. [1]). 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 include
   chemisorption corrections.

"COORDPAR atomtype e0 e1 r0 r1 CNmax Rmax Rmin"
   The COORDPAR key allows the user to add additional coordination-
   dependent parameters 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. [1] for more details. The
   coordination numbers of the DIM atoms are computed as effective
   coordination numbers. This scheme requires 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 parameters 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. [1] for more details. The code
   determines whether 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 a chemisorption correction will be
   used for the bond between the specified QM and DIM atoms. The user
   may repeat the CHEMBOND 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 the 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 which DIM atoms to include
   in the projection matrix when removing rigid motions from 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 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 the dipoles may be useful for
   debugging purposes or for calculating 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 extra information during the calculation.

"LocalDIM"
   LocalDIM performs the calculation fully for a specified CutOffDist
   (in Bohr) but ignores the image dipoles from atoms outside this
   distance. This option should converge to the full simulation, and
   the distance should be chosen with care. The LocalDirect keyword
   must also be used.

"LocalDirect"
   Tells the LocalDIM calculation to use the direct solver for the DIM
   dipoles. This solver is generally much faster than the other
   solvers available when the system is small. Currently, LocalDIM
   only works with the intention of using this faster solver, so this
   option must be used with "LocalDIM".

"DIMBUFFER"
   Default value = 0.85. This keyword tells DIM how to decide, for
   DRF2 calculations, which atoms are in the same molecular grouping.
   This grouping determines whether the charges polarize other atomic
   dipoles. For example, charges for atoms in the same molecule will
   not polarize each other, but they will polarize atomic dipoles in
   other molecules. This value should be adjusted based on the XYZ
   coordinates you supply and checked in the output file to ensure
   proper grouping behavior.

"QMBUFFER"
   Default value = 0.85. This keyword tells DIM how to decide if the
   QM atoms are "bonded" to one another in order to determine the
   force-field parameters needed for calculations requiring DIM/QM
   gradient terms. As with the DIMBUFFER keyword, this value should be
   adjusted for each system and checked in the output file to ensure
   proper behavior.

"CutOffDist"
   Distance in Bohr that determines whether the image dipoles of a DIM
   atom are included in a LocalDIM calculation. Any DIM atom that is
   inside the supplied distance from any QM atom will be fully
   included in the calculation.


DRF atomic parameters
=====================

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 belongs

"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 the group

The separation of molecules into GROUPs is important, since in the
many-body polarization operator only intermolecular interactions, i.e.
only interactions between sites that do not belong to the same group,
are included. Therefore, it is important that the combined string
(grp-nam + grp-num) is unique for each GROUP.

An example of an 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


DIM atomic parameters
=====================

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
       STATIC
       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
     LCLFLD /absolute/path/to/localfields.txt
       {freq1 freq2 freq3 ..
        .....
       }
     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 ALGORITHM 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_1  parameter (in a.u) used in CPIM.

   "OM2 val"
      OM2 specifies the \omega_2  parameter (in a.u) used in CPIM.

   "GM1 val"
      GM1 specifies the \gamma_1  parameter (in a.u) used in CPIM.

   "GM2 val"
      GM2 specifies the \gamma_2  parameter (in a.u) used in CPIM.

   "SIZE val"
      SIZE specifies the size-dependent parameter used in CPIM.

   "STATIC"
      Calculate the response properties for this atom type as static
      during a frequency-dependent calculation.

   "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 every frequency at which you are calculating is 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 will 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 that the values are 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:

         & \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}

      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 form 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 that the values are read in units of electron volts;
      otherwise, they are read in units of a.u. This key only affects
      PIM. You may give any form of the LRTZ key up to 50 times to
      supply up to 50 Lorentzian functions that make up 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 an 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, and 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 and
      the file name must end with ".xyz".

      The .xyz file is set up identically to the in-file table, except
      that there is a comment line between the number of atoms and the
      first atom. 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

"LCLFLD /absolute/path/to/localfields.txt"
   This is a text file containing the local-field-induced dipoles that
   can be read in to reduce calculation time. The file must be in the
   following format:

      freq1 freq2 freq3  ...
      1
      X
      Atom1R-DipoleX    Atom1R-DipoleY    Atom1R-DipoleZ
      Atom1I-DipoleX    Atom1I-DipoleY    Atom1I-DipoleZ
      Atom2R-DipoleX    Atom2R-DipoleY    Atom2R-DipoleZ
      Atom2I-DipoleX    Atom2I-DipoleY    Atom2I-DipoleZ
      ...
      Y
      Atom1R-DipoleX    Atom1R-DipoleY    Atom1R-DipoleZ
      Atom1I-DipoleX    Atom1I-DipoleY    Atom1I-DipoleZ
      Atom2R-DipoleX    Atom2R-DipoleY    Atom2R-DipoleZ
      Atom2I-DipoleX    Atom2I-DipoleY    Atom2I-DipoleZ
      ...

   Here freq1 freq2 freq3 ... are the frequencies that are being used
   in the response calculation, and 1 tells the program that the
   following dipoles are to be used with the first frequency. X and Y
   tell the program that these are the induced dipoles from
   perturbations in the X and Y Cartesian directions, respectively.

"GROUP groupname"
   The GROUP sub-block is where the DIM atom coordinates and
   parameters are given. This combines the "XYZ" sub-block and the
   "Element" sub-block. 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, and 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, in the format shown 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'.

[1] J.L. Payton, S.M. Morton, Justin E. Moore, and Lasse Jensen, *A
    discrete interaction model/quantum mechanical method for
    simulating surface-enhanced Raman spectroscopy*, Journal of
    Chemical Physics 136, 214103 (2012)

[2] P.T. van Duijnen and M. Swart, *Molecular and Atomic
    Polarizabilities: Thole's Model Revisited*, Journal of Physical
    Chemistry A 102, 2399 (1998)

[3] L. Jensen, P.-O. Astrand, A. Osted, J. Kongsted and K.V.
    Mikkelsen, *Polarizability of molecular clusters as calculated by
    a dipole interaction model*, Journal of Chemical Physics 116, 4001
    (2002)

[4] D.V. Chulhai and L. Jensen, *Simulating Surface-Enhanced Raman
    Optical Activity Using Atomistic Electrodynamics-Quantum
    Mechanical Models*, Journal of Physical Chemistry A 118, 9069
    (2014)

[5] Z. Hu, D.V. Chulhai, and L. Jensen, *Simulating Surface-Enhanced
    Hyper-Raman Scattering Using Atomistic Electrodynamics-Quantum
    Mechanical Models*, Journal of Chemical Theory and Computation 12,
    5968 (2016)

[6] D.V. Chulhai and L. Jensen, *Plasmonic Circular Dichroism of 310-
    and α-Helix Using a Discrete Interaction Model/Quantum Mechanics
    Method*, Journal of Physical Chemistry A 119, 5218 (2015)

[7] L. Jensen, P.T. van Duijnen and J.G. Snijders, *A discrete solvent
    reaction field model within density functional theory*, Journal of
    Chemical Physics 118, 514 (2003)

[8] L. Jensen, P.T. van Duijnen and J.G. Snijders, *A discrete solvent
    reaction field model for calculating molecular linear response
    properties in solution*, Journal of Chemical Physics 119, 3800
    (2003)

[9] L. Jensen, P.T. van Duijnen and J.G. Snijders, *A discrete solvent
    reaction field model for calculating frequency-dependent
    hyperpolarizabilities of molecules in solution*, Journal of
    Chemical Physics 119, 12998 (2003)

[10] L. Jensen, M. Swart and P.T. van Duijnen, *Microscopic and
     macroscopic polarization within a combined quantum mechanics and
     molecular mechanics model*, Journal of Chemical Physics 122,
     34103 (2005)

[11] L. Jensen, Modelling of optical response properties: Application
     to nanostructures, PhD thesis, Rijksuniversiteit Groningen, 2004.
