Table of Contents

Nuclear-electronic orbital (NEO) method

The Nuclear-electron orbital (NEO) method pioneered by Hammes-Schiffer and coworkers is available in Molpro for density fitted spin-restricted NEO-Hartree-Fock as well as a local-density fitting variant. It allows to handle a selected number of hydrogen nuclei as quantum particles by building a second Fock-matrix for the latter, coupling both subsystems (electrons and quantum protons) by a Coulomb operator. Further information about the method can be found here.

Currently, both options require the gdirect option and are not available with symmetry.

Density fitting NEO-Hartree-Fock

Using

DF-NEO-RHF, options

enables the density fitting NEO-RHF program. Through the density fitting approximation in the electronic subsystem as well as the Coulomb coupling the computational scaling for small to mid-size systems is drastically reduced. The calculation parameters can be fine tuned with the options described in the SCF program section and density fitting section. However, NEO calculations require some additional parameters explained in the following.

Local density fitting NEO-Hartree-Fock

Using

LDF-NEO-RHF, options

enables the local density fitting NEO-RHF program. The local density fitting of the electronic subsystem leads to further speed-ups in particular for large molecular systems. The specific parameters for local density fitting can be adjusted with the options given in the local density fitting Hartree-Fock section.

NEO specific options

Basis sets

The basis definition for NEO calculations must be given accordingly to the following basis block layout

 basis={
 default=minao         #Basis definition for the electronic subsystem

 set,nucbas
 default=neo-basis
 H1=pb4-f2             #Basis definition for the nuclear subsystem

 set,nucfit
 default=neo-basis
 H1=10s10p10d10f       #Basis definition for the nuclear density fitting
 }

The electronic basis set can be freely chosen from the Molpro basis set library. At the current stage no user defined mixed basis sets are possible within the NEO programs.

The nuclear basis set is defined via the nucbas keyword. The default basis for nuclear basis sets must be defined in every case as the neo-basis. Afterwards, the selected NEO centers can be assigned with the desired basis set. It is highly recommended to use the specifically tailored PB basis sets for multicomponent methods developed by Hammes-Schiffer and coworkers. Note that all NEO centers need to be assigned individually with the same basis set.

The density fitting basis for the nuclear subsystem is defined via the nucfit keyword. In order to avoid issues in basis set assignments for the classical nuclei, the default basis must be assigned as the neo-basis. Afterwards all NEO centers must be assigned the same fitting basis set (two have been included in the basis library), or a new set must be defined. For the fitting of the PB basis sets the even tempered 10s10p10d10f to 12s12p12d12f12g basis sets are recommended.

NEO centers

The desired NEO centers must be declared immediately before the NEO computation explicitly via

qnuc, H1, ...

Additionally, the chosen quantum mechanical nuclei must be given as first atoms in the geometry definition as shown for a water molecule below

3
Water molecule with one NEO center
H1  -3.5008791    1.2736107    0.7596000
O   -3.9840791    1.3301107   -0.0574000
H   -4.9109791    1.2967107    0.1521000

Starting options

In order to provide suitable starting orbitals for the NEO computation three options can be chosen.

Thresholds

The thresholds for the NEO computation can be adjusted with the following keywords

For robust convergence it is recommended to use higher thresholds for the SCF computations of both subsystems than the overall NEO energy.

Additional options

Adaptive NEO

Optimization of quantum nuclei positions with the adaptive NEO approach, where the nuclear centroids are computed on-the-fly during the SCF iterations. This procedure is available by using the

ADAPTIVE

keyword in the NEO program input card.

Threshold

The thresholds for the convergence criteria of the nuclear centers during an adaptive NEO computation can be adjusted with the following keyword

Damping

The shift of the nuclear basis function center towards the charge centroid can be damped with the following keyword

NEO examples

The first example shows the input of a DF-NEO-RHF calculation for a water molecule with two NEO centers starting with the NEOATDEN option and individual thresholds.

memory,50,m
gdirect
nosym

geometry={
3

H1  -3.5008791    1.2736107    0.7596000
H2  -4.9109791    1.2967107    0.1521000
O   -3.9840791    1.3301107   -0.0574000
}

charge=0

basis={
default=cc-pvdz

set,nucbas
default=neo-basis
H1=pb4-f2
H2=pb4-f2

set,nucfit
default=neo-basis
H1=10s10p10d10f
H2=10s10p10d10f
}

qnuc,H1,H2

{df-neo-rhf,maxdis=10,maxit=200,df_basis=cc-pvdz
neothre,1.d-6
neothrie,1.d-7
neothrin,1.d-7
neothrg,1.d-7
neothrd,1.d-7
neoatden}

The second example shows the input of a LDF-NEO-RHF computation of the same molecule starting from a prior RHF calculation. In this example a cube file is requested. This will output the quantum nuclei density.

memory,50,m
gdirect
nosym

geometry={
3

H1  -3.5008791    1.2736107    0.7596000
H2  -4.9109791    1.2967107    0.1521000
O   -3.9840791    1.3301107   -0.0574000
}

charge=0

basis={
default=cc-pvdz

set,nucbas
default=neo-basis
H1=pb4-f2
H2=pb4-f2

set,nucfit
default=neo-basis
H1=10s10p10d10f
H2=10s10p10d10f
}

{rhf}

qnuc,H1,H2

{ldf-neo-rhf,maxdis=10,maxit=200,df_basis=cc-pvdz}

{cube,nuclear.cube;density,2102.2}

The following example shows a NEO calculation, where a user-defined mixed basis set is used. Thereby, the electronic basis set at the quantum nuclei is larger than for regular hydrogen atoms. The use of the NEOMIXBAS requires the additional definition of the elebas and elefit basis sets as shown below.

memory,50,m
gdirect
nosym

geometry={
3

H1  -3.5008791    1.2736107    0.7596000
H2  -4.9109791    1.2967107    0.1521000
O   -3.9840791    1.3301107   -0.0574000
}

charge=0

basis={
default=cc-pvtz
H1=cc-pv5z

set,nucbas
default=neo-basis
H1=pb4-f2

set,nucfit
default=neo-basis
H1=10s10p10d10f

set,elebas
default=cc-pvtz
H1=cc-pv5z

set,elefit,context=jkfit
default=cc-pvtz
H1=cc-pv5z
}

qnuc,H1

{df-neo-rhf,maxdis=10,maxit=1000,df_basis=elefit
neoatden
neomixbas
}

The example below shows the input for an adaptive NEO calculation, where the nuclear basis function centers convergence is set below 1E-5 bohr and a damping factor of 0.5 is applied.

memory,50,m
gdirect
nosym

geometry={
3

H1  -3.5008791    1.2736107    0.7596000
H2  -4.9109791    1.2967107    0.1521000
O   -3.9840791    1.3301107   -0.0574000
}

charge=0

basis={
default=cc-pvdz

set,nucbas
default=neo-basis
H1=pb4-f2

set,nucfit
default=neo-basis
H1=10s10p10d10f
}

qnuc,H1

{df-neo-rhf,maxdis=10,maxit=500,df_basis=cc-pvdz
adaptive
adthres,1.d-5
addump,0.5
}

Bibliography

NEO methodology

Simon P. Webb, Tzvetelin Iordanov, and Sharon Hammes-Schiffer Multiconfigurational nuclear-electronic orbital approach: Incorporation of nuclear quantum effects in electronic structure calculations J. Chem. Phys. 2002 117 (9), 4106–4118.

Fabijan Pavošević, Tanner Culpitt, and Sharon Hammes-Schiffer Multicomponent Quantum Chemistry: Integrating Electronic and Nuclear Quantum Effects via the Nuclear–Electronic Orbital Method Chemical Reviews 2020 120 (9), 4222-4253.

PB basis sets

Qi Yu, Fabijan Pavošević, and Sharon Hammes-Schiffer Development of nuclear basis sets for multicomponent quantum chemistry methods J. Chem. Phys. 2020 152 (24), 244123.

(L)DF-NEO-RHF

Lukas Hasecke, and Ricardo A. Mata Nuclear Quantum Effects Made Accessible: Local Density Fitting in Multicomponent Methods J. Chem. Theory Comput. 2023 19 (22), 8223–8233.