Additional customization to the NEO Methods in the previous section can also be invoked. Table 13.1 shows the current customization options we have made available, and Table 13.2 shows the combination of NEO Methods upon which these features have been tested.
Feature | Single-point | Analytical | Analytical |
---|---|---|---|
Energies | Gradients | Hessian | |
Implicit Solvation | ✓ | ✓ | ✓ |
Electrostatic Embedding | ✓ | ✓ | – |
Empirical Dispersion Correction | ✓ | ✓ | ✓ |
Pseudopotentials | ✓ | ✓ | ✓ |
Level of Theory | Implicit | Electrostatic | Empirical | Pseudopotentials |
---|---|---|---|---|
(Algorithm) | Solvation | Embedding | Dispersion Correction | |
NEO-HF/DFT | ✓ | ✓ | ✓ | ✓ |
NEO-MSDFT | – | – | – | – |
CNEO | ✓ | ✓ | ✓ | ✓ |
NEO Hartree Product | ✓ | ✓ | ✓ | ✓ |
LR-NEO-TDHF/TDDFT | – | – | – | – |
RT-NEO-TDHF/TDDFT | – | ✓ | ✓ | ✓ |
NEO-SCF(V) | – | – | – | – |
NEO-CC | – | – | – | – |
NEO-MP2 | – | – | – | – |
NEO-CI | – | – | – | – |
Bulk solvent effects can be directly incorporated into NEO calculations
through the application of various implicit solvation models (Section 11.2) within the NEO framework.
1403
J. Chem. Theory Comput.
(2022),
18,
pp. 1340.
Link
The polarizable continuum model (PCM) constitutes one family of implicit solvation models and itself encompasses
several different formulations:
734
Chem. Phys. Lett.
(2011),
509,
pp. 77.
Link
,
548
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021),
11,
pp. e1519.
Link
C-PCM,
1316
Chem. Phys. Lett.
(1995),
240,
pp. 253.
Link
,
77
J. Phys. Chem. A
(1998),
102,
pp. 1995.
Link
IEF-PCM,
238
J. Chem. Phys.
(2000),
112,
pp. 5558.
Link
,
181
J. Chem. Phys.
(2001),
114,
pp. 4744.
Link
SS(V)PE,
238
J. Chem. Phys.
(2000),
112,
pp. 5558.
Link
etc.
548
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021),
11,
pp. e1519.
Link
In the PCM approach, the solute molecule is placed in a cavity that is embedded in dielectric continuum solvent, and the cavity surface
is discretized into tesserae grid points. The solvent response is represented by a partial charge centered at
each tesserae grid point .
1304
Chem. Rev.
(2005),
106,
pp. 2999.
Link
,
548
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021),
11,
pp. e1519.
Link
For each SCF iteration, the current electronic and protonic densities, together with the fixed classical nuclei, define the solute’s charge distribution. This charge distribution gives rise to the solute’s electrostatic potential at each tesserae grid point:
(13.79) |
The solute electrostatic potential is then used to compute using standard PCM methods. Once obtained, the set of tesserae charges is included as an additional one-electron (one-proton) contribution to the electronic (protonic) Fock or analogous Kohn-Sham matrix:
(13.80a) | ||||
(13.80b) |
where and refer to the gas-phase, electronic [Eq. (13.37a)] and protonic [Eq. (13.37b)] Fock or analogous Kohn-Sham matrix elements, respectively.
NEO-PCM calculations involve iterative, self-consistent convergence of the nuclear-electronic wavefunction in the presence of the dielectric continuum solvent. NEO-PCM analytic gradients and analytic Hessians are also implemented. The calculation can be invoked by setting SOLVENT_METHOD = PCM in the $rem input section. Non-electrostatic contributions to the solvation free energy can also be included by setting SOLVENT_METHOD = SMD:
SOLVENT_METHOD
The SMD approach builds upon
PCM by including non-electrostatic interactions, namely the cavitation, dispersion, and solvent structure energy terms. A caveat
is that SMD was originally parameterized to reproduce experimental free energies of solvation for conventional, electronic DFT calculations (i.e., solute nuclei represented as point charges).
274
Acc. Chem. Res.
(2008),
41,
pp. 760.
Link
In the simplest approach, the cavity surface is discretized into point charges. However, a more sophisticated approach utilizing Gaussian-smeared
charges is also supported.
731
J. Phys. Chem. Lett.
(2010),
1,
pp. 556.
Link
,
732
J. Chem. Phys.
(2010),
133,
pp. 244111.
Link
,
729
Mol. Phys.
(2020),
118,
pp. e1644384.
Link
,
548
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021),
11,
pp. e1519.
Link
Selection of these various PCM schemes and related variables can be set in the following input sections. For NEO-PCM
calculations, use the $pcm and $solvent input sections. For NEO-SMD calculations, use the $pcm and $smx input sections. The capabilities from Q-Chem’s solvation library
have been inherited to work with NEO methods, and the full list of variables that have been tested to work with NEO are listed in the pages that follow. See Section 11.2 for more comprehensive
information on these solvation models as well as consult John Herbert’s review article
548
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021),
11,
pp. e1519.
Link
to guide your selection. Lastly, examples on setting up a solvated NEO calculation can be found in Section 13.5.6.
The format of the $pcm section is analogous to the $rem section:
$pcm <Keyword> <parameter/option> $end
Note: The following job control variables belong only in the $pcm section. Do not place them in the $rem section.
Theory
Specifies the which polarizable continuum model will be used.
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
CPCM
OPTIONS:
CPCM (or CPCM1)
Conductor-like PCM of Ref.
268
J. Comput. Chem.
(2003),
24,
pp. 669.
Link
,
with .
CPCM2
Original conductor-like screening model of Ref.
, with
.
IEFPCM
IEF-PCM with an asymmetric matrix.
SSVPE
SS(V)PE model, equivalent to IEF-PCM with a symmetric matrix.
RECOMMENDATION:
We recommend C-PCM with NEO methods. The IEF-PCM model is more exact and better suited for low-dielectric solvents, but analytic Hessians are not supported.
Method
Specifies which surface discretization method will be used.
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
SwiG
OPTIONS:
SwiG
Switching/Gaussian method
Spherical
Use a single, fixed sphere for the cavity surface.
Fixed
Use discretization point charges instead of smooth Gaussians.
RECOMMENDATION:
Use the default.
Solver
Specifies the algorithm used to solve the PCM equations.
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
INVERSION
OPTIONS:
INVERSION
Direct matrix inversion
CG
Iterative conjugate gradient
RECOMMENDATION:
Use the default.
SwitchThresh
Threshold for discarding grid points on the cavity surface.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
8
OPTIONS:
Discard grid points when the switching function is less than .
RECOMMENDATION:
Use the default.
Radii
Specifies which set of atomic van der Waals radii will be used to define the solute cavity.
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
Bondi
OPTIONS:
Bondi
Use the (extended) set of Bondi radii (Ref.
1136
J. Phys. Chem.
(1996),
100,
pp. 7384.
Link
).
FF
Use Lennard-Jones radii from a molecular mechanics force field.
UFF
Use radii from the Universal Force Field (Ref.
1097
J. Am. Chem. Soc.
(1992),
114,
pp. 10024.
Link
).
Read
Read the atomic radii from a $van_der_waals input section.
RECOMMENDATION:
Do not alter this section if doing SMD calculations.
vdwScale
Scaling factor for the atomic van der Waals radii used to define the solute cavity.
INPUT SECTION: $pcm
TYPE:
FLOAT
DEFAULT:
1.2
OPTIONS:
Use a scaling factor of .
RECOMMENDATION:
Use the default.
SASradius
Form a “solvent accessible” surface with the given solvent probe radius.
INPUT SECTION: $pcm
TYPE:
FLOAT
DEFAULT:
0.0
OPTIONS:
Use a solvent probe radius of , in Å.
RECOMMENDATION:
Use the default.
SurfaceType
Selects the solute cavity surface construction.
INPUT SECTION: $pcm
TYPE:
STRING
DEFAULT:
VDW_SAS
OPTIONS:
VDW_SAS
van der Waals or solvent-accessible surface
SES
solvent-excluded surface
RECOMMENDATION:
Use the default.
Note: The acceptable values for the number of Lebedev points per sphere are , 26, 38, 50, 86, 110, 146, 170, 194, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294.
HeavyPoints
The number of Lebedev grid points to be placed non-hydrogen atoms in the QM system.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
194
OPTIONS:
Acceptable values are listed above.
RECOMMENDATION:
Use the default for geometry optimizations. For absolute solvation energies,
the user may want to examine convergence with respect to .
HPoints
The number of Lebedev grid points to be placed on H atoms in the QM system.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
110
OPTIONS:
Acceptable values are listed above.
RECOMMENDATION:
Use the default for geometry optimizations. For absolute solvation energies,
the user may want to examine convergence with respect to .
PrintLevel
Controls the printing level during PCM calculations.
INPUT SECTION: $pcm
TYPE:
INTEGER
DEFAULT:
0
OPTIONS:
0
Prints PCM energy and basic surface grid information. Minimal additional printing.
1
Level 0 plus PCM solute-solvent interaction energy components and Gauss’ Law error.
2
Level 1 plus surface grid switching parameters and a .PQR file for
visualization of the cavity surface apparent surface charges.
3
Level 2 plus a .PQR file for visualization of the electrostatic
potential at the surface grid created by the converged solute.
4
Level 3 plus additional surface grid information, electrostatic potential
and apparent surface charges on each SCF cycle.
5
Level 4 plus extensive debugging information.
RECOMMENDATION:
Use the default unless further information is desired.
CavityRadius
Specifies the solute cavity radius.
INPUT SECTION: $pcm
TYPE:
FLOAT
DEFAULT:
None
OPTIONS:
Use a radius of , in Ångstroms.
RECOMMENDATION:
None.
CavityCenter
Specifies the center of the spherical solute cavity.
INPUT SECTION: $pcm
TYPE:
FLOAT
DEFAULT:
0.0 0.0 0.0
OPTIONS:
Coordinates of the cavity center, in Ångstroms.
RECOMMENDATION:
The format is CavityCenter followed by three floating-point values,
delineated by spaces. Uses the same coordinate system as the $molecule
section.
The format of the $solvent section is analogous to the $rem section:
$solvent <Keyword> <parameter/option> $end
Note: The following job control variables belong only in the $solvent section. Do not place them in the $rem section.
Dielectric
Specifies the static dielectric constant for the PCM solvent.
INPUT SECTION: $solvent
TYPE:
FLOAT
DEFAULT:
78.39
OPTIONS:
Use a dielectric constant .
RECOMMENDATION:
The static (i.e., zero-frequency) dielectric constant is what is usually called “the”
dielectric constant. The default corresponds to water at 25C.
The format of the $smx section is analogous to the $rem section:
$smx <Keyword> <parameter/option> $end
Note: The following job control variables belong only in the $smx section. Do not place them in the $rem section.
Solvent
Sets the SM solvent
INPUT SECTION: $smx
TYPE:
STRING
DEFAULT:
water
OPTIONS:
Any name from the list of solvents given in Table 11.7.
RECOMMENDATION:
NONE
A set of external charges can be incorporated into a NEO calculation by specifying the $external_charges input section along with setting QM_MM = TRUE in the $rem input section. The format is shown below and consists of Cartesian coordinates and the value of the point charge, with one charge per line. The charge is in atomic units and the coordinates are in Ångstroms, unless bohrs are selected by setting the $rem keyword INPUT_BOHR to TRUE. The external charges are rotated with the molecule into the standard nuclear orientation and are specified in the following format:
$external_charges x-coord1 y-coord1 z-coord1 charge1 x-coord2 y-coord2 z-coord2 charge2 x-coord3 y-coord3 z-coord3 charge3 $end
Note that the electric field and analytic gradients are automatically saved to disk when the $external_charges keyword is present. This therefore enables NEO-QM/MM calculations to be performed using Q-Chem paired an external, molecular mechanics package. When using an external, molecular mechanics driver, the following $rem keyword should be set to FALSE to avoid issues with double counting:
SKIP_CHARGE_SELF_INTERACT
SKIP_CHARGE_SELF_INTERACT
Ignores the electrostatic interactions among external charges in a QM/MM calculation.
TYPE:
LOGICAL
DEFAULT:
FALSE
OPTIONS:
TRUE
No electrostatic interactions among external charges.
FALSE
Computes the electrostatic interactions among external charges.
RECOMMENDATION:
None
All of the available DFT-D empirical dispersion methods for modeling non-covalent interactions introduced in Section 5.7.3 have been extended to NEO (including analytic gradients and analytic Hessians) and can be requested via the $rem variable DFT_D, which is discussed below.
Two caveats are worth pointing out: First, note that the functional form of the DFT-D empirical dispersion potential, , is a sum over all
pairwise atomic contributions. For convenience, our implementation of the DFT-D correction uses the positions of the quantum proton(s)
basis function center(s) for computing the interatomic distances .
788
J. Chem. Phys.
(2023),
158,
pp. 114118.
Link
Second, the DFT-D dispersion potentials have been parameterized for classical nuclei represented as point charges.
DFT_D
DFT_D
Controls the empirical dispersion correction to be added.
TYPE:
LOGICAL
DEFAULT:
None
OPTIONS:
FALSE
(or 0) Do not apply the DFT-D2, DFT-CHG, or DFT-D3 scheme
EMPIRICAL_GRIMME
DFT-D2 dispersion correction from Grimme
474
J. Comput. Chem.
(2006),
27,
pp. 1787.
Link
EMPIRICAL_CHG
DFT-CHG dispersion correction from Chai and Head-Gordon
217
Phys. Chem. Chem. Phys.
(2008),
10,
pp. 6615.
Link
EMPIRICAL_GRIMME3
DFT-D3(0) dispersion correction from Grimme (deprecated as
of Q-Chem 5.0)
D3_ZERO
DFT-D3(0) dispersion correction from Grimme et al.
466
J. Chem. Phys.
(2010),
132,
pp. 154104.
Link
D3S_ZERO
DFT-D3S(0) dispersion correction from Tkachenko and Head-Gordon
1303
J. Chem. Theory Comput.
(2024),
20,
pp. 9741.
Link
D3_BJ
DFT-D3(BJ) dispersion correction from Grimme et al.
468
J. Comput. Chem.
(2011),
32,
pp. 1456.
Link
D3S_BJ
DFT-D3S(BJ) dispersion correction from Tkachenko and Head-Gordon
1303
J. Chem. Theory Comput.
(2024),
20,
pp. 9741.
Link
D3_CSO
DFT-D3(CSO) dispersion correction from Schröder et al.
1165
J. Chem. Theory Comput.
(2015),
11,
pp. 3163.
Link
D3_ZEROM
DFT-D3M(0) dispersion correction from Smith et al.
1220
J. Phys. Chem. Lett.
(2016),
7,
pp. 2197.
Link
D3_BJM
DFT-D3M(BJ) dispersion correction from Smith et al.
1220
J. Phys. Chem. Lett.
(2016),
7,
pp. 2197.
Link
D3_OP
DFT-D3(op) dispersion correction from Witte et al.
1414
J. Chem. Theory Comput.
(2017),
13,
pp. 2043.
Link
D3
Automatically select the “best” available D3 dispersion correction
D4
DFT-D4 dispersion correction from Caldeweyher et al.
172
J. Chem. Phys.
(2017),
147,
pp. 034112.
Link
,
173
J. Chem. Phys.
(2019),
150,
pp. 154122.
Link
,
174
Phys. Chem. Chem. Phys.
(2020),
22,
pp. 8499.
Link
RECOMMENDATION:
Use D4 if the specified functional is avialable. Currently, only a subset of functionals in DFT-D4 is supported.
It includes B3LYP, B97, B1LYP, PBE0, PW6B95, M06L, M06, WB97, WB97X, CAMB3LYP, PBE02, PBE0DH, MPW1K, MPWB1K, B1B95, B1PW91, B2GPPLYP, B2PLYP, B3P86, B3PW91, O3LYP, REVPBE,
REVPBE0, REVTPSS, REVTPSSH, SCAN, TPSS0, TPSSH, X3LYP, TPSS, BP86, BLYP, BPBE, MPW1PW91, MPW1LYP, PBE, RPBE, and PW91.
Note: The following variable is applicable only to the DFT-CHG dispersion model.
DFT_D_A
DFT_D_A
Controls the strength of dispersion corrections in the Chai–Head-Gordon DFT-D scheme, Eq. (5.26).
TYPE:
INTEGER
DEFAULT:
600
OPTIONS:
Corresponding to .
RECOMMENDATION:
Use the default.
Note: The following input section is only relevant to the DFT-D2 dispersion model. This input section is optional and is only used if the user wants to change the default values recommended by Grimme.
¯$empirical_dispersion ¯¯S6 S6_value ¯¯D D_value ¯¯C6 element_1 C6_value_for_element_1 element_2 C6_value_for_element_2 ¯¯VDW_RADII element_1 radii_for_element_1 element_2 radii_for_element_2 ¯$end ¯
Note: The following variables are applicable only to the DFT-D3 dispersion models.
DFT_D3_S6
DFT_D3_S6
The linear parameter in eq. (5.27). Used in all forms of DFT-D3.
TYPE:
INTEGER
DEFAULT:
100000
OPTIONS:
Corresponding to .
RECOMMENDATION:
NONE
DFT_D3_RS6
DFT_D3_S8
DFT_D3_S8
The linear parameter in Eq. (5.27). Used in DFT-D3(0),
DFT-D3(BJ), DFT-D3M(0), DFT-D3M(BJ), and DFT-D3(op).
TYPE:
INTEGER
DEFAULT:
100000
OPTIONS:
Corresponding to .
RECOMMENDATION:
NONE
DFT_D3_RS8
DFT_D3_A1
DFT_D3_A2
DFT_D3_POWER
DFT_D3_POWER
The nonlinear parameter in Eq. (5.32). Used in
DFT-D3(op). Must be greater than or equal to 6 to avoid divergence.
TYPE:
INTEGER
DEFAULT:
600000
OPTIONS:
Corresponding to .
RECOMMENDATION:
NONE
DFT_D3_3BODY
DFT_D3_3BODY
Controls whether the three-body interaction in Grimme’s DFT-D3 method should be applied
(see Eq. (14) in Ref.
466
J. Chem. Phys.
(2010),
132,
pp. 154104.
Link
).
TYPE:
LOGICAL
DEFAULT:
FALSE
OPTIONS:
FALSE (or 0)
Do not apply the three-body interaction term
TRUE
Apply the three-body interaction term
RECOMMENDATION:
NONE
Note: The following variables are applicable only to the DFT-D4 dispersion model.
DFT_D4_S6
DFT_D4_S6
The linear parameter . Used in DFT-D4.
TYPE:
INTEGER
DEFAULT:
Optimized number for the specified functional
OPTIONS:
Corresponding to .
RECOMMENDATION:
NONE
DFT_D4_S8
DFT_D4_S8
The linear parameter . Used in DFT-D4.
TYPE:
INTEGER
DEFAULT:
Optimized number for the specified functional
OPTIONS:
Corresponding to .
RECOMMENDATION:
NONE
DFT_D4_S10
DFT_D4_S10
The linear parameter . Used in DFT-D4.
TYPE:
INTEGER
DEFAULT:
Optimized number for the specified functional
OPTIONS:
Corresponding to .
RECOMMENDATION:
NONE
DFT_D4_A1
DFT_D4_A1
The nonlinear parameter . Used in DFT-D4.
TYPE:
INTEGER
DEFAULT:
Optimized number for the specified functional
OPTIONS:
Corresponding to .
RECOMMENDATION:
NONE
DFT_D4_A2
DFT_D4_A2
The nonlinear parameter . Used in DFT-D4.
TYPE:
INTEGER
DEFAULT:
Optimized number for the specified functional
OPTIONS:
Corresponding to .
RECOMMENDATION:
NONE
DFT_D4_S9
DFT_D4_S9
The linear parameter . Used in DFT-D4.
TYPE:
INTEGER
DEFAULT:
Optimized number for the specified functional
OPTIONS:
Corresponding to .
RECOMMENDATION:
NONE
DFT_D4_WF
DFT_D4_WF
Weighting factor for Gaussian weighting.
TYPE:
INTEGER
DEFAULT:
600000000
OPTIONS:
Corresponding to .
RECOMMENDATION:
Use default
DFT_D4_GA
DFT_D4_GA
Charge scaling
TYPE:
INTEGER
DEFAULT:
300000000
OPTIONS:
Corresponding to .
RECOMMENDATION:
Use default
DFT_D4_GC
DFT_D4_GC
Charge scaling
TYPE:
INTEGER
DEFAULT:
200000000
OPTIONS:
Corresponding to .
RECOMMENDATION:
Use default
Q-Chem’s Effective Core Potential (ECP) package has been integrated with the NEO method
to enable the description of relativistic and core electronic effects for systems in which some of the atoms may bear pseudopotentials.
This is done by adding an additional one-electron potential contribution, which serves to model the effects of the core electrons, into the electronic Fock or analagous
Kohn-Sham matrix (also see Section 8.10). In both the electronic and protonic Fock or analogous Kohn-Sham matrices, the one-electron term corresponding to the
interaction of electrons or quantum protons with the classical nucleus with an ECP utilizes the effective nuclear charge shielded by the core electrons.
789
J. Phys. Chem. Lett.
(2024),
15,
pp. 751.
Link
The following $rem variable controls which ECP is used:
ECP
ECP
Defines the effective core potential and associated basis set to be used
TYPE:
STRING
DEFAULT:
No ECP
OPTIONS:
General, Gen
User defined. ($ecp keyword required)
Symbol
Use standard ECPs discussed above.
RECOMMENDATION:
ECPs are recommended for first row transition metals and heavier
elements. Also consult Section 8.10 for more details.