X

Search Results

Searching....

13.5 Nuclear–Electronic Orbital Method

13.5.2 Overview of Available NEO Models

(May 21, 2025)

13.5.2.1 NEO-Hartree-Fock

The simplest method within the NEO framework is the Hartree-Fock (NEO-HF) method, where the total nuclear-electronic wavefunction is approximated as a product of electronic (𝚽0e(𝐱e)) and nuclear (𝚽0p(𝐱p)) Slater determinants composed of electronic and protonic spin orbitals, respectively:

𝚿NEO-HF(𝐱e,𝐱p)=𝚽0e(𝐱e)𝚽0p(𝐱p)=|0e0p. (13.33)

Here, 𝐱e and 𝐱p are collective spatial and spin coordinates of the quantum electrons and protons. The NEO-HF energy for a restricted Hartree-Fock (RHF) treatment of the electrons and a high-spin open-shell treatment of the quantum protons is

ENEO-HF=2iNe/2hiie+iNe/2jNe/2(2(ii|jj)-(ij|ij))+INphIIp+12INpJNp((II|JJ)-(IJ|IJ))-2iNe/2INp(ii|II). (13.34)

The i,j,, indices denote occupied spatial electronic orbitals, and the I,J,, indices correspond to occupied spatial protonic orbitals. In Eq. (13.34), hije and (ij|kl) are conventional electronic core Hamiltonian and two-electron integrals, respectively, and the corresponding terms for quantum protons are defined analogously. The last term in Eq. (13.34) is the Coulomb interaction between the electrons and the quantum protons. The spatial electronic and protonic orbitals [ψie(𝐫e) and ψIp(𝐫p)] are expanded as linear combinations of electronic or protonic Gaussian basis functions [ϕμe(𝐫e) or ϕμp(𝐫p)]:

ψie(𝐫e)= μNebfCμieϕμe(𝐫e) (13.35a)
ψIp(𝐫p)= μNpbfCμIpϕμp(𝐫p). (13.35b)

The lower-case Greek letters without and with primes denote basis functions for electrons and protons, respectively, and Cμie and CμIp are electronic and protonic MO expansion coefficients, respectively.

Analogous to the conventional electronic Hartree-Fock method, the electronic and protonic coefficients are determined by variationally minimizing the energy in Eq. (13.34) via the self-consistent field (SCF) procedure. This procedure leads to a set of coupled electronic and protonic Roothaan equations:

𝐅e𝐂e =𝐒e𝐂e𝐄e (13.36a)
𝐅p𝐂p =𝐒p𝐂p𝐄p, (13.36b)

where 𝐒e and 𝐒p are electronic and protonic overlap matrices, respectively. The electronic and protonic Fock elements in Eqs. (13.36a) and (13.36b) are given by

Fμνe =hμνe+ρλPλρe((μν|ρλ)-12(μλ|ρν))-μνPνμp(μν|μν) (13.37a)
Fμνp =hμνp+ρλPλρp((μν|ρλ)-(μλ|ρν))-μνPνμe(μν|μν). (13.37b)

The electronic and protonic density matrix elements in Eqs. (13.37a) and (13.37b) are defined as

Pνμe =2iNe/2CνieCμie* (13.38a)
Pνμp =INpCνIpCμIp*. (13.38b)

The generalization to the unrestricted Hartree-Fock (NEO-UHF) treatment of electrons is accomplished by introducing separate spatial orbitals for α and β electron spins.

The analytical gradients of the NEO-HF energy 1370 Webb S. P., Iordanov T., Hammes-Schiffer S.
J. Chem. Phys.
(2002), 117, pp. 4106.
Link
with respect to the classical nuclear coordinates (or coordinates of the centers of the quantum proton basis functions) are available. These gradients allow geometry optimizations within the NEO framework. The analytical Hessians of the NEO-HF energy with respect to the classical nuclear coordinates are also available. 1163 Schneider P. E. et al.
J. Chem. Phys.
(2021), 154, pp. 054108.
Link
The Hessians can identify whether the optimized geometries are minima or transition states on the ground state vibronic potential energy surface.

The essential $rem variables for toggling a NEO-HF calculation are NEO = TRUE and METHOD = HF. Supported JOBTYPE’s include SP, FORCE, OPT, TS, and FREQ. Consult Section 13.5.4 for the full list of variables as well as see the examples in Section 13.5.6.

13.5.2.2 NEO-DFT

NEO density functional theory (NEO-DFT) 991 Pak M. V., Chakraborty A., Hammes-Schiffer S.
J. Phys. Chem. A
(2007), 111, pp. 4522.
Link
, 223 Chakraborty A., Pak M. V., Hammes-Schiffer S.
Phys. Rev. Lett.
(2008), 101, pp. 153001.
Link
, 224 Chakraborty A., Pak M. V., Hammes-Schiffer S.
J. Chem. Phys.
(2009), 131, pp. 124115.
Link
is an extension of DFT to multicomponent systems within the NEO framework. The Hohenberg-Kohn theorems have been extended to multicomponent systems, where the reference is expressed as the product of electronic and nuclear Slater determinants composed of Kohn-Sham orbitals. The NEO-DFT total energy is

E[ρe,ρp]=Eext[ρe,ρp]+Eref[ρe,ρp]+Eexc[ρe]+Epxc[ρp]+Eepc[ρe,ρp]. (13.39)

In this equation, Eext[ρe,ρp] is the interaction of the electronic and protonic densities with the external potential created by the classical nuclei, and Eref[ρe,ρp] contains the electron–electron, proton–proton, and electron-proton classical Coulomb energies, as well as the noninteracting kinetic energies of the quantum particles. The terms Eexc[ρe], Epxc[ρp], and Eepc[ρe,ρp] are the electron-electron exchange-correlation functional, the proton-proton exchange-correlation functional, and the electron-proton correlation functional, respectively. The quantities

ρe(𝐫1e) =2i=1Ne/2|ψie(𝐫1e)|2 (13.40a)
ρp(𝐫1p) =I=1Np|ψIp(𝐫1p)|2 (13.40b)

are the electron and proton densities, respectively, and ψie(𝐫1e) and ψIp(𝐫1p) are the electronic and protonic Kohn-Sham spatial orbitals, respectively. These orbitals are obtained by solving two sets of coupled Kohn-Sham equations for the electrons and quantum protons:

(-122+veffe(𝐫1e))ψie =ϵieψie (13.41a)
(-12mp2+veffp(𝐫1p))ψIp =ϵIpψIp. (13.41b)

The effective potentials veff and veff are obtained by taking the derivative of the total energy expression in Eq. (13.39) with respect to electron density and proton density, respectively. Analogous to NEO-HF, these electronic and protonic Kohn-Sham orbitals are expanded as linear combinations of electronic or protonic Gaussian basis functions (ϕμe(𝐫e) and ϕμp(𝐫p)). The extension to open-shell electron systems is analogous to the NEO-UHF method.

The practical implementation of the NEO-DFT method requires an electron-electron exchange-correlation functional, a proton-proton exchange-correlation functional, and an electron-proton correlation functional. Any conventional electron-electron exchange-correlation functional can be used within the NEO-DFT framework. 160 Brorsen K. R., Schneider P. E., Hammes-Schiffer S.
J. Chem. Phys.
(2018), 149, pp. 044110.
Link
Because the proton-proton exchange and correlation are negligible in molecular systems, only the exchange at the NEO-Hartree-Fock level is included to eliminate self-interaction error in the NEO-DFT method. A suitable electron-proton correlation functional is essential for obtaining accurate proton densities and energies, and the epc17-2 1448 Yang Y. et al.
J. Chem. Phys.
(2017), 147, pp. 114113.
Link
, 161 Brorsen K. R., Yang Y., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2017), 8, pp. 3488.
Link
and epc19 1289 Tao Z., Yang Y., Hammes-Schiffer S.
J. Chem. Phys.
(2019), 151, pp. 124102.
Link
functionals are designed to achieve this goal. These two functionals are based on the multicomponent extension of the Colle-Salvetti formalism. The epc17-2 functional is of the local density approximation (LDA) type with the functional form:

Eepc[ρe,ρp]=-𝑑𝐫ρe(𝐫)ρp(𝐫)a-b[ρe(𝐫)ρp(𝐫)]1/2+cρe(𝐫)ρp(𝐫). (13.42)

The epc19 functional is its multicomponent generalized gradient approximation (GGA) extension that depends on the electron and proton density gradients and is of the form:

Eepc[ρe,ρp,^ρe,^ρp]=-d𝐫ρe(𝐫)ρp(𝐫)a-b[ρe(𝐫)ρp(𝐫)]1/2+cρe(𝐫)ρp(𝐫)×{1-d([ρe(𝐫)ρp(𝐫)]-1/3(1+mp)2[mp2^2ρe(𝐫)ρe(𝐫)-2mp^ρe(𝐫)^ρp(𝐫)ρe(𝐫)ρp(𝐫)+^2ρp(𝐫)ρp(𝐫)])exp[-k[ρe(𝐫)ρp(𝐫)]1/6]}. (13.43)

In addition to the parameters a, b, and c in the epc17-2 functional, 161 Brorsen K. R., Yang Y., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2017), 8, pp. 3488.
Link
the epc19 functional 1289 Tao Z., Yang Y., Hammes-Schiffer S.
J. Chem. Phys.
(2019), 151, pp. 124102.
Link
has the d and k parameters and also depends on the proton mass mp. Analogous to the NEO-HF analytical energy gradients, the NEO-DFT analytical gradients 1288 Tao Z. et al.
J. Chem. Theory Comput.
(2021), 17, pp. 5110.
Link
are available for these two functionals, allowing geometry optimizations on the ground state vibronic potential energy surface. The NEO-DFT analytical Hessians 1288 Tao Z. et al.
J. Chem. Theory Comput.
(2021), 17, pp. 5110.
Link
are available for the epc17-2 functional or when no electron-proton correlation functional is used and allow characterization of the stationary points.

The essential $rem variables for toggling a NEO-DFT calculation are NEO = TRUE and one of the two options: (1) specify METHOD or (2) specify EXCHANGE and CORRELATION. We inherit Q-Chem’s conventional, electronic DFT library, and therefore METHOD, EXCHANGE, and CORRELATION should be set to a functional listed in Section 5.3. Although it has not been tested across all possible combinations and permutations, NEO-DFT should be able to support the vast majority of the electronic XC functionals listed. Supported JOBTYPE’s include SP, FORCE, OPT, TS, and FREQ. Consult Section 13.5.4 for the full list of variables as well as see the examples in Section 13.5.6. Lastly, the electron-proton correlation functional should be set to one of the supported types:

NEO_EPC

NEO_EPC
       Specifies the electron-proton correlation functional.
TYPE:
       STRING
DEFAULT:
       No default
OPTIONS:
       NAME Use NEO_EPC = NAME, where NAME can be either epc172 or epc19.
RECOMMENDATION:
       Consult the NEO literature to guide your selection.

13.5.2.3 NEO-MSDFT

The NEO multistate DFT (NEO-MSDFT) method was developed to describe hydrogen transfer and hydrogen tunneling systems within the NEO framework. 1461 Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 10106.
Link
, 317 Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170.
Link
Analogous to the conventional electronic MSDFT method, 919 Mo Y., Bao P., Gao J.
Phys. Chem. Chem. Phys.
(2011), 13, pp. 6760.
Link
, 410 Gao J. et al.
J. Phys. Chem. Lett.
(2016), 7, pp. 5143.
Link
, 477 Grofe A. et al.
J. Chem. Theory Comput.
(2017), 13, pp. 1176.
Link
, 833 Lu Y., Gao J.
J. Phys. Chem. Lett.
(2022), 13, pp. 7762.
Link
the NEO-MSDFT method linearly combines localized NEO-DFT states in a nonorthogonal configuration interaction (NOCI) scheme. 1207 Skone J. H., Pak M. V., Hammes-Schiffer S.
J. Chem. Phys.
(2005), 123, pp. 134108.
Link
This approach captures the delocalized, bilobal proton densities needed to describe hydrogen tunneling and avoids complications of local minima in orbital space.

Consider a system with N transferring protons, where each proton moves in a double-well potential in the Born-Oppenheimer picture. In NEO-DFT calculations, the protonic density tends to localize in one well of a symmetric or nearly symmetric double-well potential instead of delocalizing over both wells. Quantizing each transferring proton with the NEO approach leads to 2N diabatic NEO-DFT states. Each diabatic state has the protonic density of each transferring proton localized in one of the two wells of its corresponding double-well potential. In practice, higher-energy diabatic states can be excluded from the NOCI expansion. The set of all diabatic NEO-DFT states is {|Ψ~0,|Ψ~1,|Ψ~n}, where n=2N-1. Each diabatic state is the product of a Kohn-Sham electronic and protonic determinant, as discussed in Section 13.5.2.2.

The adiabatic NEO-MSDFT states {|Ψ0,|Ψ1,|Ψn} are linear combinations of all diabatic NEO-DFT states:

|Ψ0 =D00|Ψ~0+D10|Ψ~1++Dn0|Ψ~n (13.44)
|Ψ1 =D01|Ψ~0+D11|Ψ~1++Dn1|Ψ~n
              
|Ψn =D0n|Ψ~0+D1n|Ψ~1++Dnn|Ψ~n

The coefficients in Eq. (13.44) are determined by solving the 2N×2N matrix eigenvalue problem

𝐇𝐃=𝐒𝐃𝐄. (13.45)

The overlap matrix 𝐒 and effective Hamiltonian matrix 𝐇 contain the overlap and couplings, respectively, between pairs of localized diabatic states. Note that the diagonal elements of 𝐒 are unity, and the diagonal elements of 𝐇 are the NEO-DFT energies of the diabatic states. For the sake of brevity, the analytical forms of the off-diagonal terms of these matrices are excluded here but are given in previous work. 317 Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170.
Link

The limitations of the epc functionals and the resulting inaccuracies of the overlap between two localized NEO-DFT states, as well as the approximate form of the off-diagonal Hamiltonian matrix elements, can be accounted for by applying a simple correction function to the off-diagonal elements of the 𝐒 matrix:

Sij=α(Sij)β. (13.46)

Users can control the values of α and β used in NEO-MSDFT calculations, and their default values were determined though a fitting process discussed in our previous work. 1461 Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 10106.
Link
Setting both these parameters to unity produces the original overlap matrix. However, in general the alpha and beta parameters should be kept at their default values unless another parameterization procedure is performed to determine alternative values.

A NEO-MSDFT calculation can be enabled by setting the $rem variable NEO_MSDFT = 1 in addition to the variables in the NEO-DFT subsection (see  13.5.2.2):

NEO_MSDFT

NEO_MSDFT
       Enable a NEO-MSDFT calculation
TYPE:
       INTEGER
DEFAULT:
       0 No NEO-MSDFT calculation.
OPTIONS:
       1 Enable a NEO-MSDFT calculation. 0 Disable a NEO-MSDFT calculation.
RECOMMENDATION:
       See Section 13.5.2.3 for details on customizing a NEO-MSDFT calculation.

Analytical gradients have also been implemented, 1464 Yu Q., Schneider P. E., Hammes-Schiffer S.
J. Chem. Phys.
(2022), 156, pp. 114115.
Link
allowing users to optimize geometries on specified NEO-MSDFT surfaces by setting JOBTYPE = OPT. By default, analytical derivative couplings between NEO-MSDFT states are computed directly following an analytical gradient calculation, 1463 Yu Q., Roy S., Hammes-Schiffer S.
J. Chem. Theory Comput.
(2022), 18, pp. 7132.
Link
but this calculation can be disabled. Note that the only gradients and derivative couplings printed out are those pertaining to the ground and first-excited NEO-MSDFT states, as these are the only states of physical relevance (as explained elsewhere 317 Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170.
Link
). Users can also request a semi-numerical Hessian calculation by setting IDERIV = 2. Further control of a NEO-MSDFT calculation is given by the $neo_msdft section, whose job-control options are given below:

OPT_STATE
       Controls which NEO-MSDFT state is geometry-optimized when JOBTYPE = OPT or which NEO-MSDFT surface the semi-numerical Hessian will be calculated on when IDERIV = 2.
INPUT SECTION: $neo_msdft
TYPE:
       INTEGER
DEFAULT:
       0 Ground-State
OPTIONS:
       n Indicates optimization/semi-numerical Hessian calculation will occur on the nth NEO-MSDFT surface.
RECOMMENDATION:
       Ensure that n is strictly less than the number of diabatic states included in the adiabatic state expansion of Eq. (13.44).

ALPHA
       Sets the α parameter used for the overlap correction of Eq. (13.46).
INPUT SECTION: $neo_msdft
TYPE:
       FLOAT
DEFAULT:
       0.0604
OPTIONS:
       User-defined.
RECOMMENDATION:
       Keep default value unless another parameterization procedure is performed.

BETA
       Sets the β parameter used for the overlap correction of Eq. (13.46).
INPUT SECTION: $neo_msdft
TYPE:
       FLOAT
DEFAULT:
       0.492
OPTIONS:
       User-defined.
RECOMMENDATION:
       Keep default value unless another parameterization procedure is performed.

NACV
       Controls if analytical derivative couplings are calculated following an analytical gradient calculation.
INPUT SECTION: $neo_msdft
TYPE:
       INTEGER
DEFAULT:
       1 Enables derivative coupling calculation.
OPTIONS:
       1 Enables derivative coupling calculation. 0 Disables derivative coupling calculation.
RECOMMENDATION:
       None.

DENPLT
       Controls the generation of proton density cube files for NEO-MSDFT states.
INPUT SECTION: $neo_msdft
TYPE:
       INTEGER
DEFAULT:
       -1 Disables generation of proton density cube files.
OPTIONS:
       n Enables generation of proton density cube files for NEO-MSDFT states {0,,n}.
RECOMMENDATION:
       Users can also generate electron density cube files for each of the specified NEO-MSDFT states via the $plots section (refer to Section 10.5.4.1 for details). Also note that if n>0, a directory neo_msdft_denplt will be created in the current working directory where all the generated cube files will be written.

CPSCF_THRESH
       Controls the convergence criteria for the NEO-CPSCF routine. Solving the NEO-CPSCF equations is required for calculating the analytical gradients of the NEO-MSDFT energies. 1464 Yu Q., Schneider P. E., Hammes-Schiffer S.
J. Chem. Phys.
(2022), 156, pp. 114115.
Link
The NEO-CPSCF routine is considered converged when the error is less than 10-CPSCF_THRESH.

INPUT SECTION: $neo_msdft
TYPE:
       INTEGER
DEFAULT:
       8
OPTIONS:
       User-defined.
RECOMMENDATION:
       Tightening the NEO-CPSCF convergence will improve the reliability of the analytical gradients.

CPSCF_NSTEPS
       Controls the maximum number of NEO-CPSCF iterations permitted. Solving the NEO-CPSCF equations is required for calculating the analytical gradients of the NEO-MSDFT energies. 1464 Yu Q., Schneider P. E., Hammes-Schiffer S.
J. Chem. Phys.
(2022), 156, pp. 114115.
Link
The NEO-CPSCF routine will fail if it does not converge within a CPSCF_NSTEPS number of steps.

INPUT SECTION: $neo_msdft
TYPE:
       INTEGER
DEFAULT:
       300
OPTIONS:
       User-defined.
RECOMMENDATION:
       If CPSCF_THRESH is higher than the default, increasing the number of iterations permitted in the NEO-CPSCF routine may be needed.

FDIFF_STEP
       Sets the distance each center is perturbed if a semi-numerical Hessian calculation is requested.
INPUT SECTION: $neo_msdft
TYPE:
       FLOAT
DEFAULT:
       10-4 Å
OPTIONS:
       User-defined.
RECOMMENDATION:
       If INPUT_BOHR = FALSE, then the step size input by the user will be measured in Å. If INPUT_BOHR = TRUE, then the step size input by the user will be measured in bohr.

As discussed in previous work, 1461 Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 10106.
Link
, 317 Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170.
Link
each quantum transferring proton in NEO-MSDFT is given two basis function centers: one center localized near the donor and another center localized near the acceptor. When requesting a NEO-MSDFT calculation, in the $molecule  section, each quantum transferring proton must have one center input as a standard H atom center and its other center input as a ghost H atom center. These centers must be listed consecutively in the $molecule section. Additionally, it is assumed that all ghost H atom centers in the $molecule section of a NEO-MSDFT calculation belong to a quantum transferring proton.

By default, all possible 2N NEO-DFT diabatic states are included in the adiabatic state expansion of Eq. (13.44). If this is not desired, users can provide an input to the $neo_msdft_diabat_control section to specify which diabatic states to include. This input section takes in an array of Boolean values (input as ones or zeros), where each row is a unique string corresponding to a particular diabatic state. These Booleans refer to whether or not a quantum transferring proton is localized on each proton basis function center in the order in which they appear in the $molecule section. Thus, if the number of diabatic states the user wants to include is m and the number of quantum transferring protons is N, the $neo_msdft_diabat_control section should contain an array with m rows and 2N columns, where the i-th row specifies which proton basis function centers have a localized quantum proton in the i-th diabatic state. Note that other protons can be quantized with NEO that do not correspond to a transferring proton and therefore are represented by a single proton basis function center; these protons should not be included in the $neo_msdft_diabat_control section. Examples showing how to set up NEO-MSDFT calculations/input files with both one and multiple transferring protons can be found in Section 13.5.6.

The overlap correction of Eq. (13.46) is parameterized to produce accurate proton densities and tunneling splittings when the proton basis function centers are optimized for separate NEO-DFT diabatic states (see the procedures explained in our previous work). 1461 Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 10106.
Link
, 317 Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170.
Link
However, users may wish to perform geometry optimizations on NEO-MSDFT surfaces instead. One can either optimize the positions of just the proton basis function centers for fixed geometries of classical nuclei or optimize the positions of both the classical nuclei and proton basis function centers simultaneously. Both cases require the user to set JOBTYPE = OPT, but in the former case, the $opt section must be used to freeze the classical nuclei during the optimization (see Section 9.4.3 for details). In either case, the input of the $molecule section serves as an initial guess for the optimization, so users must make an informed guess as to where to initially place each of the two proton basis function centers for each quantum transferring proton.

13.5.2.4 CNEO

In NEO-SCF calculations, the nuclear solution corresponds to the variational nuclear wavefunction (density) for the specific method and basis sets. The nuclear orbitals in molecular systems are localized, and therefore the quantized nuclei can be treated as distinguishable when analyzing their position expectation values. These expectation values are generally not equal to the input coordinates of the NEO nuclear basis function centers. Control over the position expectation values of the nuclei can be achieved with the constrained (C)NEO approach, 1438 Xu X., Yang Y.
J. Chem. Phys.
(2020), 152, pp. 084107.
Link
based on the original NEO method 1370 Webb S. P., Iordanov T., Hammes-Schiffer S.
J. Chem. Phys.
(2002), 117, pp. 4106.
Link
and in analogy to the conventional constrained DFT (Section 5.11). This is accomplished by imposing the following constraint:

ψIp|𝐫1p|ψIp=𝐑I (13.47)

where 𝐑I is the chosen value for the position expectation value for nucleus I. For geometry optimizations and vibrational frequency calculations, typically 𝐑I is the position of the NEO nuclear basis function center, as defined by the input coordinates. The constraint is imposed by using Lagrange multipliers. In this case, the protonic part of the coupled Kohn-Sham equations of NEO-DFT (Eq. (13.41b)) becomes:

(-12mp2+veffp(𝐫1p)+𝐟I𝐫1p)ψIp=ϵIpψIp (13.48)

where 𝐟I is the Lagrange multiplier. In practice, the value of 𝐟I is optimized numerically to satisfy the constraint in Eq. (13.47).

CNEO calculations can be invoked by setting CNEO = TRUE in addition to the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2):

CNEO

CNEO
       Enable a CNEO calculation. Using CNEO with multiple quantum protons automatically activates the nuclear Hartree product approximation.
TYPE:
       LOGICAL/INTEGER
DEFAULT:
       FALSE No CNEO calculation.
OPTIONS:
       TRUE (or 1) Enable a CNEO calculation. FALSE (or 0) Disable a CNEO calculation.
RECOMMENDATION:
       None.

Analytical CNEO gradients 1439 Xu X., Yang Y.
J. Chem. Phys.
(2020), 153, pp. 074106.
Link
and Hessians 1440 Xu X., Yang Y.
J. Chem. Phys.
(2021), 154, pp. 244110.
Link
are available, allowing geometry optimizations with JOBTYPE = OPT and vibrational frequency calculations with JOBTYPE = FREQ. CNEO can be used with Hartree-Fock or with DFT in conjunction with the epc17-2 or epc19 electron-proton correlation functional and the electron exchange-correlation functionals available in Q-Chem. Note that using CNEO with multiple quantum protons invokes the nuclear Hartree product approximation (Section 13.5.2.5).

13.5.2.5 Nuclear Hartree Product Approximation

In NEO-HF, the total wavefunction is a product of electronic and nuclear Slater determinants (Section 13.5.2.1); relatedly, in NEO-DFT, the reference is a product of Slater determinants (Section 13.5.2.2). Because the nuclear orbitals in molecular systems are highly localized and the overlap between them is effectively zero, the nuclear exchange interactions are numerically insignificant, 1003 Pavošević F., Culpitt T., Hammes-Schiffer S.
Chem. Rev.
(2020), 120, pp. 4222.
Link
and we can treat the quantum nuclei as distinguishable. This means that we can approximate the nuclear wavefunction as a Hartree product instead of a Slater determinant: 56 Auer B., Hammes-Schiffer S.
J. Chem. Phys.
(2010), 132, pp. 084110.
Link

𝚽Hartreep(𝐱p)=I=1NpχIp(𝐱Ip). (13.49)

Here, the product is over Np quantum protons, 𝐱p are the collective spatial and spin coordinates for the quantum protons, and χIp(𝐱Ip) is the spin orbital for the Ith quantum proton. The total wavefunction is then (see also Eq. (13.33)):

𝚿NEO-HF(𝐱e,𝐱p)=𝚽0e(𝐱e)I=1NpχIp(𝐱Ip). (13.50)

Note that a high-spin configuration is assumed for the protonic wavefunction, and therefore the protonic spin orbitals can be replaced with spatial orbitals multiplied by the same spin function. For NEO-HF, Eq. (13.50) leads to Np coupled one-nucleus Roothaan equations instead of one Np-nucleus Roothaan equation (Eq. (13.36b)). For NEO-DFT, it leads to Np coupled one-nucleus Kohn-Sham equations (Eq. (13.41b)). Making use of the localized nature of the quantum nuclei, the spatial nuclear orbitals for each nucleus (Eq. (13.35)) are expanded in the local basis set for that nucleus (i.e., different nuclei do not share basis functions).

The nuclear Hartree product approximation can be activated by setting NEO_N_HARTREE_PROD = TRUE in addition to the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2):

NEO_N_HARTREE_PROD

NEO_N_HARTREE_PROD
       Enable the nuclear Hartree product approximation.
TYPE:
       LOGICAL/INTEGER
DEFAULT:
       FALSE
OPTIONS:
       TRUE (or 1) Enable a nuclear Hartree product calculation. FALSE (or 0) Disable a nuclear Hartree product calculation.
RECOMMENDATION:
       Use simultaneous NEO-SCF for faster convergence.

It is recommended to use the simultaneous NEO-SCF algorithms for faster convergence (Section 13.5.4.3). By default, the Coulomb and exchange self-interaction terms for each nucleus are computed. Similarly to NEO-SCF calculations with one quantum proton, the evaluation of these terms can be turned off with NEO_VPP = FALSE to further accelerate the calculations. Note that this does not change the NEO-SCF energy and the occupied orbitals, but it does affect the nuclear virtual orbitals. 1275 Swalina C., Pak M. V., Hammes-Schiffer S.
Chem. Phys. Lett.
(2005), 404, pp. 394.
Link
Note that the supported JOBTYPE’s are the same as for CNEO 13.5.2.4 and include SP, FORCE, OPT, TS, and FREQ.

13.5.2.6 LR-NEO-TDDFT

The LR-NEO-TDDFT method 1449 Yang Y., Culpitt T., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2018), 9, pp. 1765.
Link
is a multicomponent extension of the linear response TDDFT method within the NEO framework. It allows the simultaneous calculation of the electronic and proton vibrational excitation energies. In the (LR)-NEO-TDDFT method, the linear response of the NEO Kohn-Sham system to perturbative external fields is computed. The NEO-TDDFT working equation is

[𝐀e𝐁e𝐂𝐂𝐁e𝐀e𝐂𝐂𝐂T𝐂T𝐀p𝐁p𝐂T𝐂T𝐁p𝐀p][𝐗e𝐘e𝐗p𝐘p]=ω[𝐈0000-𝐈0000𝐈0000-𝐈][𝐗e𝐘e𝐗p𝐘p] (13.51)

where

Aia,jbe =(ϵa-ϵi)δabδij+aj|ib+δ2EexcδPjbeδPaie+δ2EepcδPjbeδPaie (13.52)
Bia,jbe =ab|ij+δ2EexcδPjbeδPiae+δ2EepcδPjbeδPiae (13.53)
AIA,JBp =(ϵA-ϵI)δABδIJ+AJ|IB+δ2EpxcδPJBpδPAIp+δ2EepcδPJBpδPAIp (13.54)
BIA,JBp =AB|IJ+δ2EpxcδPJBpδPIAp+δ2EepcδPJBpδPIAp (13.55)
Cia,JB =-aB|iJ+δ2EepcδPJBpδPaie (13.56)

Here, the occupied electronic orbitals are denoted with indices i and j, whereas the unoccupied electronic orbitals are denoted with indices a and b. The analogous upper case indices denote protonic orbitals. The solution of Eq. (13.51) provides the electronic and proton vibrational excitation energies ω, as well as the transition excitation and de-excitation amplitudes, 𝐗 and 𝐘, respectively. Analogous to the TDDFT method, the Tamm-Dancoff approximation (TDA) can be imposed within the NEO framework, defining the NEO-TDDFT-TDA method that is represented by

[𝐀e𝐂𝐂T𝐀p][𝐗e𝐗p]=ω[𝐗e𝐗p]. (13.57)

The extension of the NEO-TDDFT and NEO-TDDFT-TDA approaches to open-shell electron systems is straightforward. 281 Culpitt T. et al.
J. Chem. Phys.
(2019), 150, pp. 201101.
Link
NEO-TDHF and NEO-CIS have similar forms as NEO-TDDFT and NEO-TDA without electron-proton, electron-electron, or proton-proton correlation. The analytical gradients for NEO-CIS/NEO-TDA/NEO-TDHF/NEO-TDDFT are available, 1288 Tao Z. et al.
J. Chem. Theory Comput.
(2021), 17, pp. 5110.
Link
enabling geometry optimizations on the excited state vibronic potential energy surfaces. For NEO-TDA and NEO-TDDFT, analytical gradients are available for the epc17-2 functional or when no electron-proton correlation functional is used. The transition densities can be analyzed to determine the percentages of electronic and protonic character for each vibronic excited state.

The LR-NEO-TDHF (or TDDFT) method is activated by the following keyword:

RPA

RPA
       Do a NEO-TDDFT or NEO-TDHF calculation.
TYPE:
       LOGICAL
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do a NEO-TDA or NEO-CIS calculation. TRUE Do a NEO-TDDFT or NEO-TDHF calculation.
RECOMMENDATION:
       Consult the NEO literature to guide your selection.

Furthermore, the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2) should be included. Supported JOBTYPE’s include SP and OPT. Additional customization can be controlled by the following $rem variables:

CIS_N_ROOTS

CIS_N_ROOTS
       Sets the number of NEO excited state roots to find by Davidson or display the number of roots obtained by direct diagonalization.
TYPE:
       INTEGER
DEFAULT:
       0 Do not look for any excited states.
OPTIONS:
       n n>0 Looks for n NEO excited states.
RECOMMENDATION:
       None

RPA

RPA
       Do a NEO-TDDFT or NEO-TDHF calculation.
TYPE:
       LOGICAL
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do a NEO-TDA or NEO-CIS calculation. TRUE Do a NEO-TDDFT or NEO-TDHF calculation.
RECOMMENDATION:
       Consult the NEO literature to guide your selection.

DIRECT_DIAG

DIRECT_DIAG
       Perform direct diagonalization to obtain all the NEO excitation energies.
TYPE:
       INTEGER
DEFAULT:
       0 Use Davidson algorithm.
OPTIONS:
       1 Do the direct diagonalization. 0 Use Davidson algorithm.
RECOMMENDATION:
       Only use this option when Davidson solutions are not stable.

CIS_STATE_DERIV

CIS_STATE_DERIV
       This keyword is used to specify for which NEO excited state the gradient or geometry optimization is needed.
TYPE:
       INTEGER
DEFAULT:
       No default.
OPTIONS:
       n n>0 Looks to calculate gradient or conduct geometry optimization for the nth NEO excited state.
RECOMMENDATION:
       Consult the keyword NEO_SET_ESTATE if gradient is desired for a vibronic excited state with dominant electronic character.

NEO_SET_ESTATE

NEO_SET_ESTATE
       This keyword is used to specify for which vibronic excited state with dominant electronic character the gradient or geometry optimization is needed.
TYPE:
       INTEGER
DEFAULT:
       No default.
OPTIONS:
       n n>0 Looks to calculate gradient or conduct geometry optimization for the nth NEO vibronic excited state with dominant electronic character.
RECOMMENDATION:
       Make sure enough roots are requested by the CIS_N_ROOTS keyword because the vibronic excited states with dominant protonic character usually come before.

NEO_SET_OPT

NEO_SET_OPT
       Enable a NEO excited state geometry optimization.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       1 Enable a NEO excited state geometry optimization. 0 Disable a NEO excited state geometry optimization.
RECOMMENDATION:
       Need to use with CIS_STATE_DERIV. Consult the keyword NEO_SET_ESTATE if geometry optimization is desired for a vibronic excited state with dominant electronic character.

NEO_ZVEC_LINEAR

NEO_ZVEC_LINEAR
       Use linear solver for Z-vector equations for NEO excited state gradient.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       1 Use linear solver 0 Use iterative conjugate gradient solver
RECOMMENDATION:
       Use the default iterative conjugate gradient solver because it is more memory efficient.

NEO_ZVEC_CG_MAXITER

NEO_ZVEC_CG_MAXITER
       Controls the maximum number of iterative gradient solver iterations permitted.
TYPE:
       INTEGER
DEFAULT:
       300
OPTIONS:
       n Use n>0 iterations.
RECOMMENDATION:
       None.

NEO_ZVEC_CG_CONV

NEO_ZVEC_CG_CONV
       The convergence threshold (10-NEO_ZVEC_CG_CONV) for the iterative gradient solver for NEO Z-vector equations.
TYPE:
       INTEGER
DEFAULT:
       8
OPTIONS:
       n Use n>0 iterations.
RECOMMENDATION:
       None.

SET_SUBSPACE

SET_SUBSPACE
       Specify the number of protonic guess vectors for NEO-TDDFT
TYPE:
       INTEGER
DEFAULT:
       Number of states desired (as set by CIS_N_ROOTS) if the number is smaller than the size of the protonic subspace (number of protonic occupied orbitals × number of protonic virtual orbitals) or the size of the protonic subspace
OPTIONS:
       n Use n>0 vectors.
RECOMMENDATION:
       None.

13.5.2.7 RT-NEO-TDDFT

Real-time NEO time-dependent density functional theory (RT-NEO-TDDFT or RT-NEO for brevity) is a multicomponent extension of conventional electronic RT-TDDFT within the NEO framework 1479 Zhao L. et al.
J. Phys. Chem. Lett.
(2020), 11, pp. 4052.
Link
and is an approach that enables the modeling of nonequilibrium, non-Born-Oppenheimer nuclear-electronic quantum dynamics. This approach assumes the form of the single Slater determinant product ansatz of the NEO-DFT reference system, where 𝐱e and 𝐱p are collective spatial and spin coordinates of the quantum electrons and protons, 𝐏e (or 𝐏p) is the single particle density matrix and 𝐅e (or 𝐅p) is the Kohn-Sham matrix for the electrons (or protons), and t is the time. The time-dependent Schrödinger equation

itΨNEO(𝐱e , 𝐱p ; t)=H(𝐱e , 𝐱p ; t)ΨNEO(𝐱e , 𝐱p ; t) (13.58)

can be propagated according to the set of multicomponent von Neumann equations:

it𝐏e(t)=[𝐅e(t,𝐏e(t),𝐏p(t)),𝐏e(t)] (13.59a)
it𝐏p(t)=[𝐅p(t,𝐏p(t),𝐏e(t)),𝐏p(t)] (13.59b)

Note that these equations are not propagated independently, but are actually coupled through Fock terms that depend on both particle densities. The electron and proton densities are evolved together in time. The resulting electronic and protonic dipole moments can be obtained at every time step and analyzed via post-processing to provide relevant spectral features. Fourier transformation of these signals from the corresponding time-domain to the frequency-domain can yield both electronic and vibrational spectra that have been shown to agree well with LR-NEO-TDDFT. In addition, this approach provides the nonequilbrium, real-time dynamics of the electronic and protonic densities.

Up to this point, the equations of motion presented above [Eq. (13.59a) and Eq. (13.59b) assume that the classical nuclei are held at a fixed geometry throughout the entire simulation. To combine the real-time dynamics of the electrons and protons together with the mean-field motion of the classical nuclei, the RT-NEO-TDDFT-Ehrenfest approach has been formulated. 1480 Zhao L. et al.
J. Chem. Phys.
(2020), 153, pp. 224111.
Link
In this approach, the classical nuclei are evolved according to Newtonian dynamics, while the nonequilibrium electrons and protons are still being propagated according to Eq. (13.59a) and Eq. (13.59b), respectively:

MA𝐑¨A=-AΨNEO|H^NEO|ΨNEO (13.60)

Here we denote MA and 𝐑A to be the mass and the position, respectively, of the Ath classical nucleus, and the abbreviated form of the total nuclear gradient expression is AΨNEO|H^NEO|ΨNEO. We support both of these approaches, i.e., RT-NEO with fixed classical nuclei and RT-NEO with moving classical nuclei, and we have designated those methods as METHOD = REALTIME and METHOD = EHRENFEST, respectively.

In addition to electrons, it is important that the proton dynamics is adequately described throughout the course of a given trajectory. This is especially important in cases where proton delocalization may start to occur, such as during a proton transfer reaction or when the geometry of the classical nuclei is no longer held fixed. To provide the flexibility for describing proton delocalization, one may use a larger protonic basis set and/or utilize multiple fixed proton basis (FPB) function centers placed at pre-designated positions. Alternatively, a traveling proton basis (TPB) approach has been formulated and will be made available in a future release. 1480 Zhao L. et al.
J. Chem. Phys.
(2020), 153, pp. 224111.
Link
Examples of using either one or multiple FPB function centers can be found in Section 13.5.6) of this manual .

In certain applications, the requirement of using a small time step due to the fast electronic dynamics can prohibit long simulation times. For systems that are electronically adiabatic, the electronic Born-Oppenheimer (BO) approximation can be invoked by quenching the electrons to an instantaneous electronic ground state via the standard self-consistent-field (SCF) procedure instead of propagating Eq. (13.59a) at every time step: 788 Li T. E., Hammes-Schiffer S.
J. Chem. Phys.
(2023), 158, pp. 114118.
Link

𝐏e(t)=SCF[𝐏p(t),{𝐑A(t)}] (13.61)

This approach enables the use of a time step that is several orders of magnitude larger than that required for integrating Eq. (13.59a), therefore rendering longer simulation times more tractable. The electronic BO approximation can be invoked for RT-NEO with fixed classical nuclei or for RT-NEO with moving classical nuclei by toggling the methods METHOD = BO-REALTIME or METHOD = BO-EHRENFEST, respectively.

A RT-NEO-HF or RT-NEO-TDDFT calculation can be invoked by setting NEO_TDKS = TRUE in the $rem section. In addition, the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2) should be included, and note that the JOBTYPE $rem variable is not used as it’s not applicable:

NEO_TDKS

NEO_TDKS
       Enable a real-time NEO-TDHF or NEO-TDDFT calculation.
TYPE:
       LOGICAL
DEFAULT:
       FALSE No RT-NEO calculation.
OPTIONS:
       TRUE Enable a RT-NEO calculation. FALSE Disable a RT-NEO calculation.
RECOMMENDATION:
       None.

Additional job control is provided in its own input section. The input file for a RT-NEO propagation is illustrated in several examples in Section 13.5.6, however, the available controls are also discussed below. Note that the keywords and field types defined for $tdks (Section 7.4) are not the same parameters as those defined for the $neo_tdks, and vice versa. Only parameters/keywords documented below are supported. The parameters of the jobs are set using options specified in the $neo_tdks input section. The format of the $neo_tdks section is analogous to the $rem section:

$neo_tdks
  <Keyword>  <parameter/option>
$end

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

METHOD
       Specifies which RT-NEO dynamics approach will be used to propagate the system.
INPUT SECTION: $neo_tdks
TYPE:
       STRING
DEFAULT:
       NONE
OPTIONS:
       Realtime Fixed classical nuclei 1479 Zhao L. et al.
J. Phys. Chem. Lett.
(2020), 11, pp. 4052.
Link
Ehrenfest Moving classical nuclei 1480 Zhao L. et al.
J. Chem. Phys.
(2020), 153, pp. 224111.
Link
BO-Realtime Fixed classical nuclei with electronic BO approximation 788 Li T. E., Hammes-Schiffer S.
J. Chem. Phys.
(2023), 158, pp. 114118.
Link
BO-Ehrenfest Moving classical nuclei with electronic BO approximation 788 Li T. E., Hammes-Schiffer S.
J. Chem. Phys.
(2023), 158, pp. 114118.
Link

RECOMMENDATION:
       None.

DT
       Specifies the time step Δt, in atomic units.
INPUT SECTION: $neo_tdks
TYPE:
       DOUBLE
DEFAULT:
       0.04
OPTIONS:
       Δt Δt>0 User-selected.
RECOMMENDATION:
       None.

MAXITER
       Specifies the maximum number of time steps to simulate.
INPUT SECTION: $neo_tdks
TYPE:
       INTEGER
DEFAULT:
       NONE
OPTIONS:
       n n>0 User-selected.
RECOMMENDATION:
       The total propagation length is Δt×𝐌𝐀𝐗𝐈𝐓𝐄𝐑.

ELECTRONIC_HOMO_TO_LUMO
       Performs a HOMO LUMO excitation for a single electron prior to time propagation. An unrestricted calculation is performed when this keyword is invoked.
INPUT SECTION: $neo_tdks
TYPE:
       LOGICAL
DEFAULT:
       FALSE
OPTIONS:
       TRUE FALSE
RECOMMENDATION:
       Should not use in conjunction with the electronically adiabatic approximation.

To perturb the system, a time-dependent electric field pulse 𝓔(t) can be applied via a dipole coupling term added to the electronic and/or protonic Fock matrices,

𝐅e =𝐅e,0-𝝁e𝓔(t) (13.62a)
𝐅p =𝐅p,0-𝝁p𝓔(t), (13.62b)

where 𝐅e,0 (or 𝐅p,0) is the unperturbed Kohn-Sham matrix and 𝝁e (or 𝝁p) is the dipole moment matrix for the electrons (or protons). Note that the dipole moment (i.e., 𝝁e and 𝝁p) have the opposite sign for electrons and protons. Tunable parameters for FIELD_TYPE are provided in the additional keywords found below. For the supported field types below, note that the electric field vector 𝓔=(x,y,z) is a quantity whose magnitude is controlled by FIELD_AMP and whose direction is controlled by FIELD_DIRECTION.

  1. 1.

    Delta simulates a Dirac δ-function kick with a field that is turned on only at time t0=0:

    𝓔(t)=𝐀0δ(t-t0) (13.63)
    • The amplitude of 𝐀0 and its direction are set using FIELD_AMP and FIELD_DIRECTION, respectively.

  2. 2.

    Gaussian simulates the following impulse field that is turned on at all times t0:

    𝓔(t)=𝐀0exp((t-tpeak)2τ2)sin(ωt) (13.64)
    • The amplitude of 𝐀0 and its direction are set using FIELD_AMP and FIELD_DIRECTION, respectively.

    • The center of the pulse tpeak is set using the FIELD_PEAK keyword.

    • The pulse half-width τ is set using the FIELD_TAU keyword.

    • The pulse frequency ω is set using the FIELD_FREQUENCY keyword.

FIELD_TYPE
       The external applied field
INPUT SECTION: $tdks
TYPE:
       STRING
DEFAULT:
       NONE
OPTIONS:
       DELTA δ-function kick GAUSSIAN Impulse field (Gaussian envelope) NONE No field
RECOMMENDATION:
       No recommendation.

FIELD_AMP
       The amplitude of the external field, in atomic units.
INPUT SECTION: $tdks
TYPE:
       DOUBLE
DEFAULT:
       0.01
OPTIONS:
       NONE
RECOMMENDATION:
       No recommendation.

FIELD_PEAK
       The peak position tpeak (in atomic units of time) for the Gaussian impulse field.
INPUT SECTION: $tdks
TYPE:
       DOUBLE
DEFAULT:
       0.0
OPTIONS:
       NONE
RECOMMENDATION:
       No recommendation.

FIELD_TAU
       The value of τ (in atomic units of time) for the Gaussian impulse field.
INPUT SECTION: $tdks
TYPE:
       DOUBLE
DEFAULT:
       800.0
OPTIONS:
       NONE
RECOMMENDATION:
       No recommendation.

FIELD_FREQUENCY
       The frequency ω of the external field, in eV units.
INPUT SECTION: $tdks
TYPE:
       DOUBLE
DEFAULT:
       6.0
OPTIONS:
       NONE
RECOMMENDATION:
       No recommendation.

FIELD_DIRECTION
       The direction of the external applied field vector.
INPUT SECTION: $tdks
TYPE:
       STRING
DEFAULT:
       NONE
OPTIONS:
       XYZ Pulse along the direction {1, 1, 1} X Pulse along the direction {1, 0, 0} Y Pulse along the direction {0, 1, 0} Z Pulse along the direction {0, 0, 1} NONE
RECOMMENDATION:
       No recommendation.

FIELD_PARTICLE
       The subsystem on which to apply the external field perturbation.
INPUT SECTION: $tdks
TYPE:
       STRING
DEFAULT:
       NONE
OPTIONS:
       BOTH Electrons and quantized protons ELECTRONS PROTONS
RECOMMENDATION:
       No recommendation.

13.5.2.8 NEO-DFT(V)

By quantizing select nuclei, NEO gives rise to a potential energy surface with fewer nuclear degrees of freedom, which prevents a direct calculation of the vibrational frequencies of the entire molecule. Consequently, diagonalization of a coordinate Hessian in the NEO framework yields vibrational frequencies and accompanying normal modes of only the classical nuclei, with the quantum nuclei responding instantaneously to the motion of the classical nuclei. 1163 Schneider P. E. et al.
J. Chem. Phys.
(2021), 154, pp. 054108.
Link
Although the fundamental anharmonic vibrational frequencies of the quantum nuclei can be accurately obtained through LR-NEO-TDDFT, 281 Culpitt T. et al.
J. Chem. Phys.
(2019), 150, pp. 201101.
Link
the couplings between the vibrations of the classical and quantum nuclei are missing. To obtain the fully coupled molecular vibrations, an effective strategy denoted NEO-DFT(V) was developed. 1450 Yang Y. et al.
J. Phys. Chem. Lett.
(2019), 10, pp. 1167.
Link
, 282 Culpitt T. et al.
J. Chem. Theory Comput.
(2019), 15, pp. 6840.
Link
The NEO-DFT(V) method has been shown to incorporate key anharmonic effects in full molecular vibrational analyses and to produce accurate molecular vibrational frequencies compared to experiments. 1450 Yang Y. et al.
J. Phys. Chem. Lett.
(2019), 10, pp. 1167.
Link
, 282 Culpitt T. et al.
J. Chem. Theory Comput.
(2019), 15, pp. 6840.
Link

The NEO-DFT(V) method involves diagonalization of an extended NEO Hessian composed of partial second derivatives of the coordinates of the classical nuclei (𝐫c) and the expectation values of the quantum nuclei (𝐫q). This extended Hessian matrix is composed of three sub-matrices: 𝐇0=(2E/𝐫c2)|𝐫c=𝐫q, 𝐇1=2E/𝐫q𝐫c, and 𝐇2=2E/𝐫q2, where in each case, all other coordinates of the classical nuclei and expectation values of the quantum nuclei are fixed. The extended Hessian has the following structure:

𝐇DFT(V)=[𝐇0𝐇1𝐇1𝐇2] (13.65)

where

𝐇2 =𝐔𝛀𝐌𝐔 (13.66a)
𝐇1 =-𝐇2𝐑 (13.66b)
𝐇0 =𝐇NEO+𝐑𝐇2𝐑. (13.66c)

The quantity 𝐑=d𝐫q/d𝐫c and the NEO Hessian matrix is 𝐇NEO=d2E/d𝐫c2 (without the constraint that the expectation values of the quantum nuclei are fixed). In the expression for the 𝐇2 matrix, 𝐌 is the diagonal mass matrix, and 𝛀 is the diagonal matrix with elements corresponding to the squares of the LR-NEO-TDDFT fundamental vibrational frequencies. 1450 Yang Y. et al.
J. Phys. Chem. Lett.
(2019), 10, pp. 1167.
Link
𝐔 is a unitary matrix that transforms 𝛀 to the target coordinate system and is approximated with the transition dipole moment vectors afforded by a LR-NEO-TDDFT calculation. 282 Culpitt T. et al.
J. Chem. Theory Comput.
(2019), 15, pp. 6840.
Link

Diagonalization of 𝐇DFT(V) produces the fully coupled molecular vibrational frequencies including anharmonic effects associated with the quantum protons. NEO-SCF(V) calculations can be invoked by setting NEO_SCFV = TRUE:

NEO_SCFV

NEO_SCFV
       Enable a NEO-SCFV calculation
TYPE:
       INTEGER
DEFAULT:
       0 No NEO-SCFV calculation.
OPTIONS:
       1 Enable a NEO-SCFV calculation. 0 Disable a NEO-SCFV calculation.
RECOMMENDATION:
       None.

In addition, the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2) should be included, and JOBTYPE should be set to FREQ. Furthermore, the NEO-HF(V) method involves building the extended NEO-Hessian based on the NEO-HF Hessian and inputs from NEO-TDHF, and therefore the variables specific to LR-NEO-TDHF (TDDFT) ( 13.5.2.6) can be also be used to tune the calculation.

13.5.2.9 NEO-CC

An alternative route for inclusion of correlation effects between quantum particles (i.e., electrons and protons) is with wave functions methods that are systematically improvable and parameter-free. 1370 Webb S. P., Iordanov T., Hammes-Schiffer S.
J. Chem. Phys.
(2002), 117, pp. 4106.
Link
, 1003 Pavošević F., Culpitt T., Hammes-Schiffer S.
Chem. Rev.
(2020), 120, pp. 4222.
Link
Among the various developed multicomponent wave function methods, the NEO coupled cluster (NEO-CC) methods have been particularly successful. 1002 Pavošević F., Culpitt T., Hammes-Schiffer S.
J. Chem. Theory Comput.
(2018), 15, pp. 338.
Link
, 1004 Pavošević F., Hammes-Schiffer S.
J. Chem. Phys.
(2019), 151, pp. 074104.
Link
, 1005 Pavošević F., Rousseau B. J. G., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 1578.
Link
, 1006 Pavošević F., Tao Z., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2021), 12, pp. 1631.
Link
The NEO-CC wave function is given by

|ΨNEO-CC=eT^|0e0p, (13.67)

where T^ is the cluster operator that incorporates the correlation effects between quantum particles, and |0e0p is the NEO-HF reference wave function. In the NEO-CCSD method, 1002 Pavošević F., Culpitt T., Hammes-Schiffer S.
J. Chem. Theory Comput.
(2018), 15, pp. 338.
Link
the cluster operator is given by

T^=T^1+T^2=taiaia+tAIaIA+14tabijaijab+14tABIJaIJAB+taAiIaiIaA=αtαaα, (13.68)

where aα=aα={aia,aIA,aijab,aIJAB,aiIaA} are the excitation operators expressed in terms of creation/annihilation (ap/ap) fermionic operators, and α is the excitation rank. Here, the i,j, indices denote occupied electronic orbitals, the a,b, indices denote unoccupied electronic orbitals, and the p,q, indices denote general electronic orbitals. The protonic orbitals are denoted analogously using the capitalized indices. The unknown tα wave function parameters (amplitudes) are determined by solving the set of nonlinear equations for each α: 1002 Pavošević F., Culpitt T., Hammes-Schiffer S.
J. Chem. Theory Comput.
(2018), 15, pp. 338.
Link

0e0p|aαe-T^1-T^2H^NEOeT^1+T^2|0e0p=0. (13.69)

In this equation, HNEO=hqpapq+12grspqapqrs+hQPapq+12gRSPQaPQRS-gqQpPapPqQ is the second-quantized NEO Hamiltonian, where hqp=q|h^e|p and grspq=rs|pq are conventional electronic core Hamiltonian and two-electron integrals, respectively. The remaining protonic (hQP and gRSPQ) and electron-proton (gqQpP) integrals are defined analogously. Lastly, the NEO-CCSD energy is calculated from

ENEO-CCSD=0e0p|e-T^1-T^2H^NEOeT^1+T^2|0e0p. (13.70)

To increase the computational efficiency and reduce the memory requirements for the NEO-CCSD method, the two-particle integrals can be approximated with the density fitting (DF) approximation, 1006 Pavošević F., Tao Z., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2021), 12, pp. 1631.
Link
in which the two-particle four-center integrals are factorized into a sum of products of three-center and two-center two-particle integrals. In particular, the four-center two-electron integrals are approximated by

(μν|ρσ)=μρ|νσXY(μν|X)(X|Y)-1(Y|ρσ), (13.71)

where (μν|X) and (X|Y) are three-center and two-center two-electron integrals, respectively. In this equation, μ,ν, and X,Y, indices denote electronic and auxiliary electronic basis functions, respectively. The four-center two-proton integrals are approximated analogously by

(μν|ρσ)=μρ|νσXY(μν|X)(X|Y)-1(Y|ρσ), (13.72)

where primed indices denote protonic basis functions and (μν|X) and (X|Y) are three-center and two-center two-proton integrals, respectively. Finally, the four-center electron-proton integrals are approximated as

(μν|μν)=μμ|ννXY(μν|X)(X|Y)-1(Y|μν). (13.73)

By employing the DF approximation, the memory requirements for storing four-center two-particle integrals are reduced from Nbf4 to Nbf2×Naux, where Nbf and Naux are the number of electronic or protonic basis functions and auxiliary basis functions, respectively. 1006 Pavošević F., Tao Z., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2021), 12, pp. 1631.
Link

Currently only JOBTYPE = SP is supported. NEO-CC is activated by the following keyword in addition to the variables for NEO-HF (NEO = TRUE and METHOD = HF):

NEO_RICCSD

NEO_RICCSD
       Enable a NEO-RICCSD calculation.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       1 Enable this option. 0 Disable this option.
RECOMMENDATION:
       Both electronic and protonic auxiliary basis sets must be specified.

The following additional $rem variables can also be used to customize the NEO-RICCSD calculation (also see the examples in Section 13.5.6):

NEO_CCSD_MAX_CYCLES

NEO_CCSD_MAX_CYCLES
       Controls the maximum number of CC iterations permitted.
TYPE:
       INTEGER
DEFAULT:
       5000
OPTIONS:
       n Set the maximum number of iterations to n>0.
RECOMMENDATION:
       None

NEO_CCSD_CONVERGENCE

NEO_CCSD_CONVERGENCE
       NEO-RICCSD is considered converged when the energy error is less than 10-NEO_CCSD_CONVERGENCE.
TYPE:
       INTEGER
DEFAULT:
       8
OPTIONS:
       User-defined
RECOMMENDATION:
       None

13.5.2.10 NEO-MP2

Second-order Møller-Plesset perturbation theory provides a useful framework for wave function based correlation effects with lower cost than NEO-CC 1275 Swalina C., Pak M. V., Hammes-Schiffer S.
Chem. Phys. Lett.
(2005), 404, pp. 394.
Link
. The NEO-MP2 correlation energy is composed of the electron-electron, proton-proton and electron-proton contributions:

ENEO(2)=Eee(2)+Epp(2)+Eep(2), (13.74)

where

Eee(2) =14ijab|ij|ab-ij|ba|2εi+εj-εa-εb (13.75a)
Epp(2) =14IJAB|IJ|AB-IJ|BA|2εI+εJ-εA-εB (13.75b)
Eep(2) =iIaA|iI|aA|2εi+εI-εa-εA. (13.75c)

Here, i and j are occupied electronic spin orbitals with energies εi and εj, a and b are virtual (unoccupied) electronic spin orbitals with energies εa and εb, and uppercase indices indicate the analogous protonic spin orbitals and energies.

To provide accuracy that is competitive with NEO-CCSD, NEO-MP2 has been extended to include orbital optimization and empirical spin-component and electron-proton correlation scaling. 1005 Pavošević F., Rousseau B. J. G., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 1578.
Link
, 1006 Pavošević F., Tao Z., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2021), 12, pp. 1631.
Link
This leads to the NEO equivalents of SOS-OOMP2 (O2) 823 Lochan R. C., Head-Gordon M.
J. Chem. Phys.
(2007), 126, pp. 164101.
Link
and SOS-MP2 471 Grimme S.
J. Chem. Phys.
(2003), 118, pp. 9095.
Link
, 637 Jung Y. et al.
J. Chem. Phys.
(2004), 121, pp. 9793.
Link
(Sections 6.6.5 and 6.6.6). The NEO equivalents are referred to as NEO-SOS-OOMP2 and NEO-SOS-MP2 to indicate the inclusion of an electron-proton correlation scaling factor:

ENEO(2)=cssEeess(2)+cosEeeos(2)+cepEep(2)+Epp(2), (13.76)

where Eeess(2) and Eeeos(2) are the same- and opposite-spin parts of the electron-electron correlation energy, which are scaled by css and cos, respectively, and cep is the electron-proton scaling factor.

Orbital optimization is performed by minimizing the NEO-MP2 energy with respect to the unitary operator eX^-X^ formed from the rotation operator X^=X^e+X^p=xaiaia+xAIaIA, with 𝐱={xai,xAI} being the set of unknown electronic and protonic orbital rotation parameters. Optimal orbital rotation parameters are solved for by finding the stationary points of electronic and protonic orbital gradients 𝐰e and 𝐰p, which have elements

(𝐰e)ai=ENEO-MP2(𝐱)xia|𝐱e=0 (13.77a)
and
(𝐰p)AI=ENEO-MP2(𝐱)xIA|𝐱p=0, (13.77b)

and are solved in an alternating fashion until self-consistency is achieved. For computational efficiency, all NEO-MP2 methods are implemented with the density fitting approximation for two-particle integrals, 371 Fetherolf J. H. et al.
J. Phys. Chem. Lett.
(2022), 13, pp. 5563.
Link
as described in more detail in Section 13.5.2.9.

Currently only JOBTYPE = SP is supported. NEO-MP2 methods are invoked with NEO_RIMP2 = 1 to run without orbital optimization and NEO_RIMP2 = 2 to run with orbital optimization. Additionally, the variables for NEO-HF (NEO = TRUE and METHOD = HF) should be included:

NEO_RIMP2

NEO_RIMP2
       Enable a NEO-MP2 or NEO-OOMP2 calculation.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       2 Perform a NEO-OOMP2 calclulation. 1 Perform a NEO-MP2 calculation. 0 Disable this option.
RECOMMENDATION:
       Both electronic and protonic auxiliary basis sets must be specified.

Spin-component scaling settings are controlled with the SCS variable. Custom css and cos can be set with SSS_FACTOR and SOS_FACTOR, respectively, when SCS = 3 is set (arbitrary scaling). By default, cep is set to the same value as cos, but a custom electron-proton scaling factor can be input with the EP_FACTOR variable, which will always override the default set by SCS. These settings are as follows (also see the examples in Section 13.5.6):

SCS

SCS
       Set the type of spin-component scaling.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       1 Turns on spin-component scaling with SCS (css=0.33, cos=1.2). 2 Turns on spin-component scaling with SOS (css=0.0, cos=1.2 for MP2, cos=1.3 for OOMP2). 3 arbitrary SCS (set with SSS_FACTOR and SOS_FACTOR). 0 no spin-component scaling.
RECOMMENDATION:
       NONE

SSS_FACTOR

SSS_FACTOR
       Controls the strength of the same-spin component of the MP2 electron-electron correlation energy.
TYPE:
       INTEGER
DEFAULT:
       1000000
OPTIONS:
       n Corresponding to css=n/106.
RECOMMENDATION:
       NONE

SOS_FACTOR

SOS_FACTOR
       Controls the strength of the opposite-spin component of the MP2 electron-electron correlation energy.
TYPE:
       INTEGER
DEFAULT:
       1000000
OPTIONS:
       n Corresponding to cos=n/106.
RECOMMENDATION:
       NONE

EP_FACTOR

EP_FACTOR
       Controls the strength of the electron-proton component of the NEO-MP2 correlation energy.
TYPE:
       INTEGER
DEFAULT:
       1000000
OPTIONS:
       n Corresponding to cep=n/106.
RECOMMENDATION:
       NONE

13.5.2.11 NEO Configuration Interaction

The NEO Configuration Interaction (NEO-CI) wavefunction is a linear combination of NEO configurations:

ΨNEO(𝐑c)=iIciI|Φie|ΦIn (13.78)

where |Φie is electronic determinant i and |ΦIn is nuclear determinant I. The electronic determinant is composed of electronic orbitals and the nuclear determinant is composed of protonic orbitals. 842 Malbon C., Hammes-Schiffer S.
J. Chem. Theory Comput.
(2025), 21, pp. 3968.
Link

The NEO-CI wavefunction is almost always truncated to include only single, or single and double excitations into the virtual orbital space. In the single-reference case, the electronic and nuclear determinants use the NEO-HF electronic and nuclear orbitals, respectively. In the multi-reference case, the electronic and nuclear orbitals are optimized via NEO-MCSCF. The orbitals can be optimized for multiple states with state-averaging.

Available models include:

  • NEO-MCSCF(CASSCF,RASSCF)

  • NEO-CI

  • NEO-MR-CI

  • NEO-OOMRCI

In addition to setting NEO = TRUE, METHOD = HF, and JOBTYPE = SP, the job settings for these available models are controlled by the following $rem variables:

Note:  Wavefunction Keywords

NEO_WF_N_ELECS

NEO_WF_N_ELECS
       Total number of electrons in wavefunction.
TYPE:
       INTEGER
DEFAULT:
       (Read from $molecule input.)
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_E_FROZEN_CORE

NEO_WF_E_FROZEN_CORE
       Number of electronic frozen core orbitals. Orbitals are doubly occupied in all configurations.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       Generally, 1s orbitals can be treated as frozen core orbitals.

NEO_WF_E_FROZEN_VIRT

NEO_WF_E_FROZEN_VIRT
       Number of electronic frozen virtual orbitals. These orbitals are unoccupied in all configurations.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       Virtual orbitals contribute to the dynamic correlation recovered by CI calculations, so it is not recommended that this value is set, except in special cases.

NEO_WF_E_DOCC

NEO_WF_E_DOCC
       Number of electronic reference doubly-occupied orbitals minus the CI frozen core orbitals. There are NEO_WF_E_XLEVEL excitations out of these orbitals.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_E_RAS

NEO_WF_E_RAS
       Number of electronic reference restricted active space orbitals.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_E_RAS_NHOLES

NEO_WF_E_RAS_NHOLES
       Maximum Number of electronic reference excitations out of the RAS orbitals and into the CAS and AUX orbitals.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_E_CAS

NEO_WF_E_CAS
       Number of electronic reference complete active space orbitals.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_E_AUX

NEO_WF_E_AUX
       Number of electronic reference auxilliary space orbitals.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_E_AUX_NELECS

NEO_WF_E_AUX_NELECS
       Maximum number of electronic reference excitations into the AUX orbitals from the RAS and AUX orbitals.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_E_XLEVEL

NEO_WF_E_XLEVEL
       Electronic excitation level of NEO-CI expansion: 1 = single, 2 = single and double, etc.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_N_FROZEN_VIRT

NEO_WF_N_FROZEN_VIRT
       Number of nuclear frozen virtual orbitals. These orbitals are unoccupied in all configurations.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       Virtual orbitals contribute to the dynamic correlation recovered by CI calculations, so it is not recommended that this value is set, except in special cases.

NEO_WF_N_DOCC

NEO_WF_N_DOCC
       Number of nuclear reference doubly-occupied orbitals. There are NEO_WF_N_XLEVEL excitations out of these orbitals.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_N_CAS

NEO_WF_N_CAS
       Number of nuclear reference complete active space orbitals.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_WF_N_XLEVEL

NEO_WF_N_XLEVEL
       Nuclear excitation level of NEO-CI expansion: 1 = single, 2 = single and double, etc.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

Note:  NEO-MRCI Keywords

NEO_MRCI

NEO_MRCI
       Perform NEO Multireference Configuration Interaction
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       TRUE Run a NEO-MRCI Calculation. FALSE Do not run a NEO-MRCI Calculation.
RECOMMENDATION:
       NONE

NEO_MRCI_MAX_ITER

NEO_MRCI_MAX_ITER
       Maximum NEO CI iterations.
TYPE:
       INTEGER
DEFAULT:
       100
OPTIONS:
       User-defined
RECOMMENDATION:
       NONE

NEO_MRCI_NUM_ROOTS

NEO_MRCI_NUM_ROOTS
       Number of NEO CI states to converge.
TYPE:
       INTEGER
DEFAULT:
       1
OPTIONS:
       User-defined
RECOMMENDATION:
       Due to the spectrum and structure of the NEO CI Hamiltonian, it is advisable to set NEO_MRCI_NUM_ROOTS to be greater than the number of desired roots.

NEO_MRCI_KSPACE_MIN

NEO_MRCI_KSPACE_MIN
       Minimum size of Krylov subspace.
TYPE:
       INTEGER
DEFAULT:
       NEO_MRCI_NUM_ROOTS + 1
OPTIONS:
       User-defined
RECOMMENDATION:
       Increasing the number of guess vectors improves CI convergence but has increased memory costs.

NEO_MRCI_KSPACE_MAX

NEO_MRCI_KSPACE_MAX
       Maximum size of Krylov subspace.
TYPE:
       INTEGER
DEFAULT:
       NEO_MRCI_KSPACE_MIN + 5
OPTIONS:
       User-defined
RECOMMENDATION:
       Increasing the number of guess vectors before contraction improves CI convergence but has increased memory costs.

NEO_MRCI_PREDIAG_ROUTINE

NEO_MRCI_PREDIAG_ROUTINE
       Choose how initial guess vectors are created.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       0 Linear combination of lowest NEO_MRCI_PREDIAG_DIM NEO-CI diagonal elements. 1 Eigenvectors of NEO_MRCI_PREDIAG_DIM x NEO_MRCI_PREDIAG_DIM block of NEO-CI Hamiltonian. 2 MCSCF vectors are used for initial guess.
RECOMMENDATION:
       0

NEO_MRCI_PREDIAG_DIM

NEO_MRCI_PREDIAG_DIM
       Dimension for NEO_MRCI_PREDIAG_ROUTINE
TYPE:
       INTEGER
DEFAULT:
       1000
OPTIONS:
       User-defined
RECOMMENDATION:
       The larger the size, the better the initial guess.

NEO_MRCI_CONV_THRESH

NEO_MRCI_CONV_THRESH
       Convergence tolerance for roots of NEO-CI. It is defined as 10-NEO_MRCI_CONV_THRESH
TYPE:
       INTEGER
DEFAULT:
       5
OPTIONS:
       4 Sufficient for energies. 5 Gradients and couplings.
RECOMMENDATION:
       4

NEO_MRCI_PRINT_CIVECS

NEO_MRCI_PRINT_CIVECS
       Print the final NEO CI vectors to civec.dat.
TYPE:
       BOOLEAN
DEFAULT:
       TRUE
OPTIONS:
       User-defined
RECOMMENDATION:
       None

NEO_MRCI_READ_CIVECS

NEO_MRCI_READ_CIVECS
       Read in NEO CI vectors from civec.dat.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       User-defined
RECOMMENDATION:
       This can be used to restart from a previously unconverged CI calculation.

NEO_MRCI_READ_MOS

NEO_MRCI_READ_MOS
       Read in MO coefficients from mocoefs.electronic and mocoefs.protonic.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       User-defined
RECOMMENDATION:
       None

NEO_MRCI_PRINT_LEVEL

NEO_MRCI_PRINT_LEVEL
       Determines print level for NEO CI calculation.
TYPE:
       INTEGER
DEFAULT:
       1
OPTIONS:
       0 No convergence information is printed. Only iteration number and timings. 1 Convergence information is printed. >1 Detailed printing for debugging.
RECOMMENDATION:
       1

NEO_MRCI_CIVEC_ANALYSIS

NEO_MRCI_CIVEC_ANALYSIS
       Perform brief analysis for CI vectors, largest contributing configurations.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE No analysis is performed. TRUE Analysis is performed.
RECOMMENDATION:
       FALSE

NEO_MRCI_CIDEN_PRINT

NEO_MRCI_CIDEN_PRINT
       Print the NEO CI protonic densities of the converged roots to ⁢.cube files
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do not print protonic densities. TRUE Print protonic densities.
RECOMMENDATION:
       None

Note:  NEO-MCSCF Keywords

NEO_MCSCF

NEO_MCSCF
       Run a NEO-MCSCF calculation.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do not run a NEO-MCSCF. TRUE Run a NEO-MCSCF.
RECOMMENDATION:
       None

NEO_MCSCF_DEBUG

NEO_MCSCF_DEBUG
       Run NEO-MCSCF in debug mode.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       0 Do not debug NEO-MCSCF. 1 Debug NEO-MCSCF
RECOMMENDATION:
       0

NEO_MCSCF_MACRO_MAX_ITER

NEO_MCSCF_MACRO_MAX_ITER
       Maximum number of NEO-MCSCF macroiterations.
TYPE:
       INTEGER
DEFAULT:
       50
OPTIONS:
       User-defined
RECOMMENDATION:
       Most CAS calculations should require few macroiterations. More complex RASSCF wavefunctions may require more iterations as convergence can be slow.

NEO_MCSCF_MACRO_CONV_THRESH

NEO_MCSCF_MACRO_CONV_THRESH
       Orbital gradient convergence tolerance for NEO-MCSCF. Defined as 10-NEO_MCSCF_MACRO_CONV_THRESH.
TYPE:
       INTEGER
DEFAULT:
       5
OPTIONS:
       User-defined
RECOMMENDATION:
       5

NEO_MCSCF_MICRO_MAX_ITER

NEO_MCSCF_MICRO_MAX_ITER
       Maximum NEO-CI iterations.
TYPE:
       INTEGER
DEFAULT:
       1000
OPTIONS:
       User-defined
RECOMMENDATION:
       Set as high as necessary to converge all roots required.

1000

NEO_MCSCF_MICRO_OPT_MAX_ITER

NEO_MCSCF_MICRO_OPT_MAX_ITER
       Maximum number of orbital optimizations per macroiteration.
TYPE:
       INTEGER
DEFAULT:
       5
OPTIONS:
       User-defined
RECOMMENDATION:
       5

NEO_MCSCF_EMC_CONV_THRESH

NEO_MCSCF_EMC_CONV_THRESH
       Convergence threshold for NEO-MCSCF energy. Defined as 10-NEO_MCSCF_EMC_CONV_THRESH.
TYPE:
       INTEGER
DEFAULT:
       8
OPTIONS:
       User-defined
RECOMMENDATION:
       8

NEO_MCSCF_NSTATES

NEO_MCSCF_NSTATES
       Number of NEO-CI roots to compute during microiterations.
TYPE:
       INTEGER
DEFAULT:
       1
OPTIONS:
       User-defined
RECOMMENDATION:
       None

NEO_MCSCF_SA

NEO_MCSCF_SA
       Perform state-averaged NEO-MCSCF.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE Perform CASSCF. TRUE Perform SA-MCSCF
RECOMMENDATION:
       None

NEO_MCSCF_SA_NSTATES

NEO_MCSCF_SA_NSTATES
       Number of states to include in state-averaging.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       User-defined
RECOMMENDATION:
       None

NEO_MCSCF_PRINT_CIVECS

NEO_MCSCF_PRINT_CIVECS
       Print the NEO-MCSCF CI vectors to mcvec.dat.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE No printing. TRUE MCSCF CI vectors are printed to mcvec.dat and wavefunction is printed to mcwf.dat.
RECOMMENDATION:
       None

NEO_MCSCF_READ_MOS

NEO_MCSCF_READ_MOS
       Read MCSCF MO coefficients from mocoefs.electronic and mocoefs.protonic.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE NEO-HF MO coefficients used for initial guess. TRUE MO coefficients from files used as initial guess.
RECOMMENDATION:
       None

NEO_MCSCF_PRINT_MOS

NEO_MCSCF_PRINT_MOS
       Print MCSCF MO coefficients to mocoefs.electronic.new and mocoefs.protonic.new.
TYPE:
       BOOLEAN
DEFAULT:
       TRUE
OPTIONS:
       TRUE Print MO coefficients. FALSE MO coefficients are not printed.
RECOMMENDATION:
       TRUE

NEO_MCSCF_PRINT_LEVEL

NEO_MCSCF_PRINT_LEVEL
       Print level for NEO-MCSCF.
TYPE:
       INTEGER
DEFAULT:
       1
OPTIONS:
       0 Only convergence information is printed. 1 Macro and micro iteration information is printed. 2 Detailed printing.
RECOMMENDATION:
       1

NEO_MCSCF_OPT_METHOD

NEO_MCSCF_OPT_METHOD
       Orbital optimization method for NEO-MCSCF.
TYPE:
       INTEGER
DEFAULT:
       2
OPTIONS:
       0 GDM macroiterations. No microiterations. 1 Augmented Hessian macroiterations. Augmented Hessian microiterations. 2 Augmented Hessian macroiterations. GDM microiterations.
RECOMMENDATION:
       2

NEO_MCSCF_KSPACE_MIN

NEO_MCSCF_KSPACE_MIN
       Minimum size of the Krylov subspace for NEO-CI microiterations.
TYPE:
       INTEGER
DEFAULT:
       NEO_MCSCF_NSTATES + 2
OPTIONS:
       User-defined
RECOMMENDATION:
       None

NEO_MCSCF_KSPACE_MAX

NEO_MCSCF_KSPACE_MAX
       Maximum size of the Krylov subspace for NEO-CI microiterations.
TYPE:
       INTEGER
DEFAULT:
       NEO_MCSCF_KSPACE_MIN + 10
OPTIONS:
       User-defined
RECOMMENDATION:
       None

NEO_MCSCF_CIDEN_PRINT

NEO_MCSCF_CIDEN_PRINT
       Print the CI protonic density for MCSCF states to *.cube files.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do not print density. TRUE Print density.
RECOMMENDATION:
       None

NEO_MCSCF_REORDER_MOS

NEO_MCSCF_REORDER_MOS
       Reorder initial guess molecular orbitals. Requires reorder_mo.txt and reorder_nuc_mo.txt.
TYPE:
       BOOLEAN
DEFAULT:
       FALSE
OPTIONS:
       FALSE Orbitals are not reordered. TRUE Orbitals are reordered.
RECOMMENDATION:
       None