3.31 NEO

This functionality was tested to our best knowledge, but there is a possibility for remaining bugs or inconsistencies. If you have any questions regarding usage and/or functionality, if you encouter problems/bugs, or if you have suggestions - please contact jxu@unc.edu or ryzhou@live.unc.edu. Please cite our publication of NEO-DFT[317], RT-NEO-TDDFT[316], CNEO-DFT[202] when you publish results obtained with these functionalities.

The current implementation incorporates a large set of control keywords of which only the most important are noted here. Several keywords controlling debugging/development features are not captured here but can be found in the code (e.g. neo/src/neo_variable.f90). Please contact us if you want to learn more about this.

Currently, one can do NEO-DFT simulations for finite and periodic systems. All modes of parallelism are supported.

Theory

Multicomponent DFT via NEO method for Kohn-Sham system

The NEO formalism in the framework of multicomponent DFT leads to a set of coupled electron-proton Kohn-Sham (KS) equations,

	$\displaystyle\hat{H}^{e}_{\bf{k}}\psi^{e}_{i,\bf{k}}(\bm{r}^{e})$	$\displaystyle=\left[-\frac{1}{2}\nabla^{2}+v_{\text{eff}}^{e}(\bm{r}^{e})% \right]\psi^{e}_{i,\bf{k}}=\epsilon^{e}_{i,\bf{k}}\psi^{e}_{i,\bf{k}}(\bm{r}^{% e}),$		(3.114)
	$\displaystyle\hat{H}^{p}\psi^{p}_{i}(\bm{r}^{p})$	$\displaystyle=\left[-\frac{1}{2M^{p}}\nabla^{2}+v_{\text{eff}}^{p}(\bm{r}^{p})% \right]\psi^{p}_{i}=\epsilon^{p}_{i}\psi^{p}_{i}(\bm{r}^{p}),$		(3.115)

where $M^{p}$ is the proton mass, $\psi^{e}_{i,\bf{k}}$ and $\epsilon^{e}_{i,\bf{k}}$ are the KS orbitals and eigenvalues of the electrons, and $\psi^{p}_{i}$ and $\epsilon^{p}_{i}$ are the KS orbitals and eigenvalues of the protons. The effective potentials are defined as

	$\displaystyle v_{\text{eff}}^{e}(\bm{r}^{e})$	$\displaystyle=v_{\text{ext}}(\bm{r}^{e})+v_{\text{es}}^{e}(\bm{r}^{e})-v_{% \text{es}}^{p}(\bm{r}^{e})+\frac{\delta E_{\text{xc}}^{e}[\rho^{e}]}{\delta% \rho^{e}}+\frac{\delta E_{\text{epc}}[\rho^{e},\rho^{p}]}{\delta\rho^{e}},$		(3.116)
	$\displaystyle v_{\text{eff}}^{p}(\bm{r}^{p})$	$\displaystyle=-v_{\text{ext}}(\bm{r}^{p})-v_{\text{es}}^{e}(\bm{r}^{p})+v_{% \text{es}}^{p}(\bm{r}^{p})+\frac{\delta E_{\text{xc}}^{p}[\rho^{p}]}{\delta% \rho^{p}}+\frac{\delta E_{\text{epc}}[\rho^{e},\rho^{p}]}{\delta\rho^{p}}.$		(3.117)

where $v_{\text{ext}}$ , $v_{\text{es}}^{e}$ , and $v_{\text{es}}^{p}$ are electrostatic potentials for the classical nuclei, electrons, and quantum protons. $E_{\text{xc}}^{e}$ and $E_{\text{xc}}^{p}$ are the exchange-correlation energies of the electrons and the quantum protons, respectively. $E_{\text{epc}}$ is the correlation energy between the electrons and quantum protons, which is given as a functional of both the electron density and the proton density in the multicomponent DFT formalism. In most cases, the KS wave functions for the quantum protons are extremely localized in real space such that the overlap of proton wave functions from different unit cells is close to zero. Thus, the Brillouin zone (BZ) integration can be neglected for the protons even for extended systems, whereas the electron KS wave functions require the BZ integration. With Equations 3.114 and 3.115, the ground state of the multicomponent system of electrons and quantum protons needs to be solved self-consistently in the presence of the electrostatic potential from the classical nuclei. Further discussions can be found in our publication [317].

Constrained NEO-DFT method for Kohn-Sham system

Following by the NEO formalism, the CNEO method is to impose a constraint such that the expectation value of the position operator for quantum protons gives a specified point $\bm{R}^{p}_{n}$ such that

\langle\psi_{n}^{p}|\hat{\bm{r}}|\psi_{n}^{p}\rangle=\bm{R}^{p}_{n}

(3.118)

The Lagrangian function can be written as

	$\displaystyle L=$	$\displaystyle E[\rho^{e},\rho^{p}]+\sum_{n}\bm{f}_{n}\cdot(\langle\psi_{n}^{p}% \|\hat{\bm{r}}\|\psi_{n}^{p}\rangle-\bm{R}^{p}_{n}\langle\psi_{n}^{p}\|\psi_{n}^{% p}\rangle)$		(3.119)
		$\displaystyle+\int d\mathbf{k}\sum_{n}\varepsilon_{n\mathbf{k}}^{e}(\langle% \psi^{e}_{n\mathbf{k}}\|\psi^{e}_{n\mathbf{k}}\rangle-1)+\sum_{n}\varepsilon_{n% }^{p}(\langle\psi^{p}_{n}\|\psi^{p}_{n}\rangle-1)$		(3.119)

where $\bm{f}_{n}$ is the Lagrange multiplier of the constraint on the $n$ -th proton. We then obtain a set of new eigenvalue equations for proton orbitals. The proton KS equation is

\displaystyle\left[-\frac{1}{2M^{p}}\nabla_{p}^{2}+v_{\mathrm{eff}}^{p}(\bm{r}% ^{p})+\bm{f}_{n}\cdot(\hat{\bm{r}}-\bm{R}^{p}_{n})\right]\psi_{n}^{p}(\bm{r}^{% p})=\epsilon_{n}^{p}\psi_{n}^{p}(\bm{r}^{p}).

(3.120)

Since the position constraints only act on quantum protons, the form of the electronic KS equation remains unchanged. The first and second-order derivative of $L$ with respect to $\mathbf{f}$ could be expressed as

	$\displaystyle\frac{dL}{d\bm{f}_{n}}=$	$\displaystyle\frac{\partial L}{\partial\bm{f}_{n}}+\sum_{m}\frac{\partial L}{% \partial\psi^{p}_{m}}\frac{d\psi^{p}_{m}}{d\bm{f}_{n}}+\sum_{m}\frac{\partial L% }{\partial\psi^{p}_{m}}\frac{d\psi^{p}_{m}}{d\bm{f}_{n}}$		(3.121)
	$\displaystyle=$	$\displaystyle\langle\psi^{p}_{n}\|\hat{\bm{r}}\|\psi^{p}_{n}\rangle-\bm{R}^{p}_{% n}\langle\psi^{p}_{n}\|\psi^{p}_{n}\rangle.$		(3.121)

\frac{\delta^{2}L}{\delta\bm{f}_{n}^{2}}{\approx}2\sum_{m{\neq n}}^{\text{% unocc}}\frac{\langle\psi_{m}^{p}|\hat{\bm{r}}|\psi_{n}^{p}\rangle\langle\psi_{% n}^{p}|\hat{\bm{r}}|\psi_{m}^{p}\rangle}{\epsilon^{p}_{m}-\epsilon_{n}^{p}},% \quad\frac{\delta L}{\delta\bm{f}_{m}\delta\bm{f}_{n}}{\approx}0

(3.122)

So the Lagrange multipliers could be determined iteratively using Newton-Raphson method,

\bm{f}^{n+1}=\bm{f}^{n}-[\frac{d^{2}L}{d\bm{f}^{2}}]^{-1}\frac{dL}{d\bm{f}}.

(3.123)

When the Lagrange multipliers are converged (i.e. the Lagrangian is stationary with respect to $\bm{f}$ ), the expectation values of the position operator for protons coincide with the constrained positions. We also implemented a Wannier process here to handle occasions that proton density itinerant or delocalized with multiple geometric centers. Further discussions can be found in our publication [202].

Coupled quantum dynamics of electrons and protons via RT-NEO-TDDFT methods

The Lagrangian for the coupled quantum electron-proton system with classical nuclei is

$\displaystyle L(t)$	$\displaystyle=\int\text{d}\bm{r}^{e}\int\text{d}\mathbf{k}\sum_{j}{\psi_{j% \mathbf{k}}^{e}}^{*}(\bm{r}^{e},t)\left[i\frac{\partial}{\partial t}+\frac{1}{% 2m^{e}}\nabla_{\bm{r}^{e}}^{2}\right]\psi_{j\mathbf{k}}^{e}(\bm{r}^{e},t)$	(3.124)
	$\displaystyle-\frac{1}{2}\int\text{d}\bm{r}^{e}\text{d}\bm{r}^{\prime e}\frac{% e^{2}}{\|\bm{r}^{e}-\bm{r}^{\prime e}\|}\rho^{e}(\bm{r}^{e},t)\rho^{e}(\bm{r}^{% \prime e},t)-E_{XC}^{e}[\rho^{e}]$
	$\displaystyle+\int\text{d}\bm{r}^{p}\sum_{j}{\psi_{j}^{p}}^{*}(\bm{r}^{p},t)% \left[i\frac{\partial}{\partial t}+\frac{1}{2M^{p}}\nabla_{\bm{r}^{p}}^{2}% \right]\psi_{j}^{p}(\bm{r}^{p},t)$
	$\displaystyle-\frac{1}{2}\int\text{d}\bm{r}^{p}\text{d}\bm{r}^{\prime p}\frac{% e^{2}}{\|\bm{r}^{p}-\bm{r}^{\prime p}\|}\rho^{p}(\bm{r}^{p},t)\rho^{p}(\bm{r}^{% \prime p},t)-E_{XC}^{p}[\rho^{p}]$
	$\displaystyle+\frac{1}{2}\int\text{d}\bm{r}^{e}\text{d}\bm{r}^{p}\frac{e^{2}}{% \|\bm{r}^{e}-\bm{r}^{p}\|}\rho^{e}(\bm{r}^{e},t)\rho^{p}(\bm{r}^{p},t)-E_{EPC}[% \rho^{e},\rho^{p}]$
	$\displaystyle+\sum_{I}^{N}\frac{1}{2}M_{I}\left[\frac{d}{dt}\bm{R}_{I}(t)% \right]^{2}-\sum_{I<J}\frac{Z_{I}Z_{J}e^{2}}{\|\bm{R}_{I}(t)-\bm{R}_{J}(t)\|}$
	$\displaystyle-\int\text{d}\bm{r}^{p}\rho^{p}(\bm{r}^{p},t)\sum_{I}^{N}\frac{Z_% {I}e^{2}}{\|\bm{r}^{p}-\bm{R}_{I}(t)\|}+\int\text{d}\bm{r}^{e}\rho^{e}(\bm{r}^{e% },t)\sum_{I}^{N}\frac{Z_{I}e^{2}}{\|\bm{r}^{e}-\bm{R}_{I}(t)\|}$

where $\psi_{j\mathbf{k}}^{e}(\bm{r}^{e},t)$ , $\psi_{j}^{p}(\bm{r}^{p},t)$ , and $\rho^{e/p}(\bm{r}^{e/p},t)$ are the time-dependent Kohn-Sham (KS) orbitals and the quantum-mechanical probability density for either electron or proton degrees of freedom. $m^{e}$ and $M^{p}$ are the electron and proton mass. $E_{XC}^{e/p}[\rho^{e/p}]$ is the electron or proton exchange-correlation energy functional, and $E_{EPC}[\rho^{e},\rho^{p}]$ is the correlation energy functional between quantum electrons and protons. $M_{I}$ , $\bm{R}_{I}(t)$ , and $Z_{I}$ are the mass, the coordinates, and the charge, respectively, of classical nucleus $I$ . In the RT-NEO-TDDFT formulation, the equations of motion are obtained only for the quantum-mechanical degrees of freedom, namely the electrons and the quantum protons by applying the Dirac-Frenkel variational principle to the action. The resulting time-dependent Kohn-Sham (TDKS) equations govern the coupled quantum dynamics of the electrons and protons. Within our periodic NEO-DFT framework[317], the electron and proton TDKS wave functions are propagated on an equal footing,

	$\displaystyle i\frac{\partial}{\partial t}\psi^{e}_{i,\bf{k}}(\bm{r}^{e},t)=$	$\displaystyle\hat{H}^{e}\psi^{e}_{i,\bf{k}}(\bm{r}^{e},t)=\left[-\frac{1}{2m^{% e}}\nabla^{2}+\hat{v}_{\text{eff}}^{e}(\bm{r}^{e},t)+\hat{v}_{\text{ext}}(t)% \right]\psi^{e}_{i,\bf{k}}(\bm{r}^{e},t)$		(3.125)
	$\displaystyle i\frac{\partial}{\partial t}\psi^{p}_{i}(\bm{r}^{p},t)=$	$\displaystyle\hat{H}^{p}\psi^{p}_{i}(\bm{r}^{p},t)=\left[-\frac{1}{2M^{p}}% \nabla^{2}+\hat{v}_{\text{eff}}^{p}(\bm{r}^{p},t)-\hat{v}_{\text{ext}}(t)% \right]\psi^{p}_{i}(\bm{r}^{p},t).$		(3.126)

The external potential, $\hat{v}_{\text{ext}}(t)$ , includes not only the terms that derive from the electrostatic interaction from the classical atomic nuclei but also the (time-dependent) external electromagnetic field in the length gauge. $\hat{v}_{\text{eff}}$ is the effective potential that includes the Hartree potential (Coulomb interaction terms), the exchange-correlation (XC) potential, and the electron-proton correlation (EPC) potential. Generally speaking, the XC potential could have a dependence on all former time, but we adopt the adiabatic approximation such that all the history dependence is ignored, removing the explicit time dependence of $v_{\text{XC}}$ . Further discussions can be found in our publication [316].

Tags for general section of control.in

Core Settings

Tag: NEO_type(control.in)

Usage: NEO_type type
Purpose: Performing a NEO-DFT simulation
type: String that specifies the type of NEO calculation to be performed.

•

0/none: Do not include NEO-DFT. (Default)
•

1/scf_post: NEO-DFT as a post-SCF procdure. (Use DFT electron as starting point)
•

2/rt-tddft: RT-NEO-TDDFT.
•

3/scf: NEO-DFT together with the conventional SCF loop.

The overall switch to turn on a NEO-DFT calculation in FHI-aims. Must set manually to enable NEO simulations.

Tag: NEO_basis_type(control.in)

Usage: NEO_basis_type type
type: String that specifies the type of basis used for the quantization of protons.

•

1/primitive: Primitive spherical gaussian basis whose exponent and weight are provided later by NEO_basis keyword for each individual protons. (Default)
•

2/even_tempered/ET: Even-tempered spherical gaussian basis whose radial part is defined as $\phi_{i,l}(r)=C_{i}*r^{l}*e^{-\alpha*\beta^{i-1}*r^{2}}$ , where $i$ runs from 1 to $N$ (NEO_basis_n), $l$ runs from 0 to $l_{max}$ (NEO_basis_l_max), $C_{i}$ is the automatically generated normalize constant, $\alpha$ and $\beta$ is defined by NEO_basis_alpha and NEO_basis_beta keyword.

Tag: NEO_basis_preset_type(control.in)

Usage: NEO_basis_preset_type type
type: String that specifies the basis set used for the quantization of protons if NEO_basis_type is set to 3/’preset’.

•

1/pb4d: The PB4-D basis set.
•

2/pb4f2: The PB4-F2 basis set. (Default)
•

3/pb5g: The PB5-G basis set.

The preset basis set are optimized Gaussian-type nuclear basis sets defined in [323].

Tag: NEO_density_fitting_type(control.in)

Usage: NEO_density_fitting_type/NEO_df_type type
type: String that specifies the type of auxiliary basis used for the Resolution of the Identity (RI) method.

•

0/aims: Not using RI method. Instead, electrostatic potential are calculated with a similar approach as described in [36], and exchange are calculated with four center two electron integrals between selected proton pairs.
•

2/even_tempered/ET: Even-tempered spherical gaussian basis whose radial part is defined as $\psi_{i,l}(r)=C_{i}*r^{l}*e^{-\alpha*\beta^{i-1}*r^{2}}$ , where $i$ runs from 1 to $N$ (NEO_df_n), $l$ runs from 0 to $l_{max}$ (NEO_df_l_max), $C_{i}$ is the automatically generated normalize constant, $\alpha$ and $\beta$ is defined by NEO_df_alpha and NEO_df_beta keyword. (Default)

By default, RI are used for calculating both coulomb and exchange terms. One can optionally use different ways to calculate exchange terms. See NEO_RI_exchange_type keyword for more detail. When NEO_density_fitting_type is set to 0 or ’aims’, density fitting for proton are done as a similar manar as electrons on grid points to evaluate coulomb interactions and atomic forces. In this case, proton exchange are calculated by the exact exchange with a radius cutoff.

Tag: NEO_rho_proton_type(control.in)

Usage: NEO_rho_proton_type type
type: String that specifies how proton density on grid points are calculated.

•

0/eigenvector: Using proton wavefunctions.
•

1/density_matrix: Using proton density matrix. (Default)
•

11/density_matrix_fast: Using proton density matrix but ultilize a sparsed storage and algorithm for linear algebra.
•

2/density_fitting/df: Using the density fitting approach.

Tag: NEO_epc_type(control.in)

Usage: NEO_epc_type type
type: Integer that specifies the type of eletron proton correlation (EPC) functionals.

•

0: Not using EPC functional. (Default)
•

171: Using the epc17-1 functional [45].
•

17 (or 172): Using the epc17-2 functional [318].

It has been shown that EPC functionals have significant impact on the total energy of the calculations. Including EPC will result in a better estimation of the spacial distribution of quantum protons as well as the zero-point energy caused by the quantization of protons.

Tag: NEO_initial_guess_type(control.in)

Usage: NEO_initial_guess_type type
type: String that specifies how the initial guess for proton state are determined.

•

1/potential: Starting the simulation from solving the proton states where protons can no feel each other. This option works great with single quantized proton cases.
•

2/read: Starting the simulation from eigenvectors provided in file "neo_eigenvectors.in".
•

3/random: Starting the simulation from n random eigenvectors, where n is the number of quantized protons. (Default)

In some special cases, a bad initial guess can lead to the simulation end up with a local minimum, where protons are very delocalized. When that happened, one can either change the NEO_initial_guess_type option or reorder the proton basis provided in NEO_basis to get a different initial guess.

Tag: NEO_RI_exchange_type(control.in)

Usage: NEO_RI_exchange_type type
type: String that specifies the type of RI method used especially for calculating the exchange terms in the simulation.

•

1/V: Using the RI-V method. (Should be avoid in extended systems)
•

3/LVL : Using the RI-LVL method. (Default)
•

32/LVL_fast: Using the RI-LVL method but ultilize a sparsed storage and algorithm for linear algebra. This option will be slower than option 3 for small systems but faster for large systems.
•

5/noRI: Not using RI method at all. This should be only used when once use NEO_density_fitting_type 0.
•

51/noRI_V: Directly calculating exchange terms using the wavefunctions of protons, but use RI-V to calculate electrostatic potentials.
•

53/noRI_LVL: Directly calculating exchange terms using the wavefunctions of protons, but use RI-LVL to calculate electrostatic potentials.

One can compare the results using option 3 and option 53 to check the performance of the current choice of auxiliary basis for the RI method.

Tag: NEO_occ_proton_type(control.in)

Usage: NEO_occ_proton_type type
type: Integer that specifies how the proton occupation number is treated.

•

1: Ground state. One protron per state starting from the lowest energy. (Default)
•

2: Read from FILE neo_occ.in, where first line indicate the number of occupied states, and the second line gives the occupation numbers from low energy state to high energy state separated by space.

Proton self consistent field loop (SCF) settings

Tag: NEO_scf_proton_max_step(control.in)

Usage: NEO_scf_proton_max_step N
Purpose: Defines the maximium proton SCF steps performed before updating the proton contribution and going back to electron SCF loops, if the proton SCF never meets the convergence criteria.
N is a integer. Default: 500.

The default value is set to be large for safety. Generally, N around 20 should be enough for convergence.

Tag: NEO_scf_accuracy_rho(control.in)

Usage: NEO_scf_proton_max_rho value
Purpose: Convergence criterion of proton SCF, based on proton charge density.
value is a small positive real number, against which the volume-integrated root-mean-square of the proton charge density between the present and previous proton SCF is checked. Default: 1e-6.

This convergence criterion is the most important SCF convergence for proton SCF cycles, one many consider to set this depend on the total number of protons.

Tag: NEO_scf_accuracy_eev(control.in)

Usage: NEO_scf_proton_max_eev value
Purpose: Convergence criterion of proton SCF, based on the change of proton eigenvalues.
value is a small positive real number, against which the change of the sum of proton eigenvalue between the present and previous proton SCF is checked. Default: 1e-3.

Tag: NEO_scf_accuracy_dm(control.in)

Usage: NEO_scf_proton_max_step value
Purpose: Convergence criterion of proton SCF, based on the change of proton density matrix.
value is a small positive real number, against which the largest changed element () in proton density matrix between the present and previous proton SCF is checked. Default: 1e-3.

This convergence criterion can be redundant when proton basis has linear dependence.

Tag: NEO_scf_accuracy_err(control.in)

Usage: NEO_scf_proton_max_err value
Purpose: Convergence criterion of proton SCF, based on the change of the error matrix defined by the DIIS procedure.
value is a small positive real number, against which the largest changed element in the error matrix between the present and previous proton SCF is checked. Default: 1e-5.

Tag: NEO_scf_diis_type(control.in)

Usage: NEO_scf_diis_type type
type: String that specifies the type of DIIS method used to enhance the convergence of proton scf cycle.

•

0/none: Not using DIIS method. Use the calculated proton Hamiltonian without any mixing method.
•

1/DIIS_new_notho: Using a altered DIIS method to update proton Hamiltonian before solving the proton KS equation.
•

2/DIIS_new: Similar to option 1, but using a orthonormal residual. This option has better performance when there are linear dependence in proton basis.
•

3/DIIS_notho: Using the original DIIS method to update proton Hamiltonian before solving the proton KS equation.
•

4/DIIS: Similar to option 4, but using a orthonormal residual. This option has better performance when there are linear dependence in proton basis. (Default)
•

5/simple_mixing: Simple mixing method. Mix $p$ of the new proton Hamiltonian with $1-p$ of the previous one. $p$ is set by neo_scf_diis_simple_mixing_G0.

For the simlations with extremely small proton basis sets (e.g. only $1s$ orbital), the DIIS protpagtion may fail. In that case, one can run the simulation without DIIS method.

Tag: neo_scf_diis_simple_mixing_G0(control.in)

Usage: neo_scf_diis_simple_mixing_G0 value
Purpose: Mixing parameter used by the simple mixing procedure.
value is a small positive real number. Default: 0.25.

Tag: neo_scf_diis_simple_mixing_G1(control.in)

Usage: neo_scf_diis_simple_mixing_G1 value
Purpose: Mixing parameter used by the DIIS procedure for mixing the new Hamiltonian with the previous one.
value is a small positive real number. Default: 1.0.

Tag: NEO_scf_diis_max_step(control.in)

Usage: NEO_scf_diis_max_step N
Purpose: Defines how many proton SCF steps are stored for the DIIS procedure. All steps are stored for the first N steps, and then replace the oldest with the newer one.
N is a integer. Default: 8.

Generally, the default value is sufficient for DIIS. One may increase N if the DIIS procedure performs badly.

Tag: NEO_constrain_position(control.in)

Usage: NEO_constrain_position type
Purpose: Controls whether using CNEO-DFT method and defines how the protons positions are constrained during the SCF procedure.
type is a string that determines the constraint mode in CNEO-DFT calculation.

•

F/.false.: No position constraint. (Default)
•

T/.true./force: Calculate the CNEO-DFT by a determined force.
•

iter/iterative: Iteratively solve the force of constrained proton.
•

iter-wannier: Iteratively solve the force of constrained proton with wannier process.

The wannier process should be use when protons are close. Currently the CNEO-DFT have conflicts with DIIS process, so the DIIS option should turn off accordingly.

Numerical cutoff values for speeding up

Tag: NEO_proton_hartree_rcut(control.in)

Usage: NEO_proton_hartree_rcut value
Purpose: Cutoff radius, within which the electrostatic potential of a quantum proton is evaluated and updated for SCF loops. For extended systems with Ewald method, all real space evaluations of electrostatic potential are limited in this radius as well.
value is a real number in Å. Default: 30 Bohr.

NEO_proton_hartree_rcut is the most important parameter in the NEO-DFT work flow and must be set. Based on tests, 15 to 20 Åis a reasonable choice for the balence of accuracy and performance.

Tag: NEO_proton_hartree_rclassic(control.in)

Usage: NEO_proton_hartree_rclassic value
Purpose: Cutoff radius, within which the electrostatic potential of a quantum proton is calculated analytically. Beyond the radius, the electrostatic potential is approximated by classical multipoles located at the quantum proton center.
value is a real number in Å. Default: -1.0, the radius cutoff for the multipole approximation is set automatically to ensure that beyond this radius the proton charge density is less than NEO_proton_hartree_charge_threshold_rclassic for all quantum protons.

This parameter should be less than NEO_proton_hartree_rcut. With the default setting, FHI-aims will print out NEO_proton_hartree_rclassic used in the calculation in its output. One can choose a larger number to ensure better numerical results.

Tag: NEO_basis_rcut(control.in)

Usage: NEO_basis_rcut value
Purpose: Cutoff radius, beyond which the value of proton wavefunctions are set as zero.
value is a real number in Å. Default: -1.0, the radius cutoff for proton basis is set automatically to ensure that beyond this radius its value is less than NEO_proton_basis_wf_threshold_rcut for all basis functions.

With the default setting, FHI-aims will print out NEO_basis_rcut used in the calculation in its output. One can choose a larger number to ensure better numerical results. Due to the localized nature of quantum proton, this value can be less than 5 Å.

Tag: NEO_exchange_rcut(control.in)

Usage: NEO_exchange_rcut value
Purpose: Cutoff radius, beyond which the exchange effect between the quantum proton pair is ignored.
value is a real number. Default: 4.

As discussed in our paper[317], the exchange term is mainly for counteract the self-interaction of quantum protons. So one can limited the exchange effect within a very small radius. It has been tested that, value can be set to very close to 0 if no ghost/empty protons are included. Otherwize, it is better to increase value to cover all ghost centers for a quantum proton.

Tag: NEO_rho_cut(control.in)

Usage: NEO_rho_cut value
Purpose: Cutoff value, less than which proton charge density is set as zero.
value is a small real number. Default: 1e-30.

Tag: NEO_proton_hartree_charge_threshold_rclassic(control.in)

Usage: NEO_proton_hartree_charge_threshold_rclassic value
Purpose: Parameter used to determine NEO_proton_hartree_rclassic by default. This will only be referenced if NEO_proton_hartree_rclassic is not set.
value is a small positive real number. Default: 1e-30.

Tag: NEO_proton_basis_wf_threshold_rcut(control.in)

Usage: NEO_proton_basis_wf_threshold_rcut value
Purpose: Parameter used to determine NEO_basis_rcut by default. This will only be referenced if NEO_basis_rcut is not set.
value is a small positive real number. Default: 1e-30.

Ewald method settings

Tag: NEO_use_ewald(control.in)

Usage: NEO_use_ewald flag
Purpose: Determines whether the Ewald method are used to calculate electrostatic potentials.
flag is a logical string, either .false. or .true.. Default: .false. for isolated systems, .true. for extended systems.

By default, for molecular cases, FHI-aims only update the electrostatic potential of a given quantum proton within a limited radius defined by NEO_proton_hartree_rcut. For periodical cases, the electrostatic potential on all grid points are updated with the Ewald method.

Tag: NEO_ewald_l_max(control.in)

Usage: NEO_ewald_l_max $l_{max}$
$l_{max}$ : Integer that specifies the highest angular momentum used in the Ewald Summation. Higher angular momentums are treated locally within NEO_proton_hartree_rcut. Default: same value as NEO_df_l_max.

Proton basis and auxiliary basis for Resolution of the Identity

Tag: NEO_basis_n(control.in)

Usage: NEO_basis_n N
N: Integer that specifies the rank of the Even-tempered spherical gaussian basis. See NEO_basis_type keyword for more detail. Default: 8.

Tag: NEO_basis_l_max(control.in)

Usage: NEO_basis_l_max $l_{max}$
$l_{max}$ : Integer that specifies the highest angular momentum of the Even-tempered spherical gaussian basis. See NEO_basis_type keyword for more detail. Default: 2.

Tag: NEO_basis_alpha(control.in)

Usage: NEO_basis_alpha alpha
alpha: Real value that specifies the inital exponent of the Even-tempered spherical gaussian basis. See NEO_basis_type keyword for more detail. Default: $2*\sqrt{2}$ .

Tag: NEO_basis_beta(control.in)

Usage: NEO_basis_beta beta
beta: Real value that specifies the interval of exponents of the Even-tempered spherical gaussian basis. See NEO_basis_type keyword for more detail. Default: $\sqrt{2}$ .

Tag: NEO_df_n(control.in)

Usage: NEO_df_n N
N: Integer that specifies the rank of the Even-tempered spherical gaussian auxiliary basis. See NEO_df_type keyword for more detail. Default: 10.

Tag: NEO_df_l_max(control.in)

Usage: NEO_df_l_max $l_{max}$
$l_{max}$ : Integer that specifies the highest angular momentum involved. See NEO_df_type keyword for more detail. Default: 3.

This parameter is also used when NEO_df_type is set to 0/.

Tag: NEO_df_alpha(control.in)

Usage: NEO_df_alpha alpha
alpha: Real value that specifies the inital exponent of the Even-tempered spherical gaussian auxiliary basis. See NEO_df_type keyword for more detail. Default: $2*\sqrt{2}$ .

Tag: NEO_df_beta(control.in)

Usage: NEO_df_beta beta
beta: Real value that specifies the interval of exponents of the Even-tempered spherical gaussian auxiliary basis. See NEO_df_type keyword for more detail. Default: $\sqrt{2}$ .

Tag: NEO_basis(control.in)

Usage: NEO_basis
Purpose: Defines the basis of quantum protons in NEO-DFT simulations.

Tag: NEO_basis_end(control.in)

Usage: NEO_basis_end
Purpose: End the defination of quantum proton basis in NEO-DFT simulations.

NEO_basis is use to invoke the defination of proton basis for individual protons, one can add NEO_basis_end to end the defination of basis set for safety. Examples to use these keywords are provided in "neo.in_example".

NEO_basis sub-tag: H(control.in)

Usage: H index pos_x pos_y pos_z
              [ l_1 alpha_1 weight_1 ]
              [ l_2 alpha_2 weight_2 ]
              [ … ]
              [ l_N alpha_N weight_N ]
Purpose: Adds a (ghost) quantum proton center in NEO-DFT calculation
index is a integer for the index of hydrogen atom (proton) declared by atom in geometry.in. Optionally, negtive number represents for adding a ghost center that do not have charge distributions, which usually works together with empty in geometry.in.
pos_x, pos_y, pos_z are optional parameters that redefine the position of proton center. By default, the position of proton center are read from geomerty.in if optionpos_x, pos_y, pos_z are not provided. Providing new position for proton centers are still in beta phase and not guaranteed to work.
l_i is a character, specifying the angular momentum (s, p, d, f, …) of the $i$ th gaussian type orbital.
alpha_i is a real number [in bohr^-2], specifying the exponents $\alpha_{i}$ of the $i$ th gaussian type orbital.
weight_i is a real number, specifying the weight of the $i$ th gaussian type orbital. For primitive gaussian basis, weight should be 1.0.

RT-NEO-TDDFT

RT-NEO-TDDFT simulations are invoked by setting NEO_type to "rt-tddft" along with other conventional RT-TDDFT settings such as RT_TDDFT_propagation. RT-NEO-TDDFT simulations are mostly handled by conventional RT-TDDFT Tags without necessity for additional NEO inputs.

Tag: NEO_rt_propagate_n_step(control.in)

Usage: NEO_rt_propagate_n_step N
Purpose: Defines how often proton propagation will be performed RT-NEO-TDDFT simulations.
N: Integer. Proton wavefuntion will be propagated every N electron steps. Default: 1.

Generally, N can be set to greater than one due to the larger mass of protons as compared to electrons. Energy conservation should be checked if N is set to a big value.

Tag: NEO_rt_ext_e_type(control.in)

Usage: NEO_rt_ext_e_type type
type: String that specifies the type of the external potential used in RT-NEO-TDDFT calculations.

•

-1/none: Ingore all external potentials inputs (read from electron systems) for the proton system.
•

0/scalar/length: Using scalar potential for protons, inputs are read from the electron RT-TDDFT setting. See RT_TDDFT_td_field for more detail. (Default)

Tag: NEO_rt_propagator_type(control.in)

Usage: NEO_rt_propagator_type type
type: String that specifies the type of the external potential used in RT-NEO-TDDFT calculations.

•

1/EM_iter: The exponential midpoint method. And solve iteratively until reach convergence criteria set by NEO_rt_precor_tol.
•

11/EM/exponential_midpoint: The exponential midpoint method. (Default)
•

3/RK4/runge_kutta_4: The Standard explicit Runge-Kutta 4 scheme. Works only for very small time steps but is a good choice for doing accurate reference simulations because its properties are well-known.

See RT_TDDFT_propagator for more descriptions.

Tag: NEO_rt_precor_tol(control.in)

Usage: NEO_rt_precor_tol tol
Purpose: Defines the convergence criteria of iterative EM propagator in RT-NEO-TDDFT simulations.
tol: Small Real value that specifies the convergence criteria. Default: 1e-5.

Tag: NEO_rt_update_occ_proton(control.in)

Usage: NEO_rt_update_occ_proton flag
Purpose: Determines whether to update proton occupation from file "neo_occ.in" after NEO-DFT before running RT-NEO-TDDFT.
flag is a logical string, either .false. or .true.. Default: .false..

If set to .true.., one needs to provide the "neo_occ.in" file, whose first line contains a integer neo_n_occ_max indicating max number of the neo states, and second line contains neo_n_occ_max count of positive real numbers indicating occupation number for each states. For example,
neo_n_occ_max
occ_1 occ_2 … occ_neo_n_occ_max

Tag: NEO_rt_update_occ_electron(control.in)

Usage: NEO_rt_update_occ_electron flag
Purpose: Determines whether to update electron occupation from file after NEO-DFT before running RT-NEO-TDDFT.
flag is a logical string, either .false. or .true.. Default: .false..

If set to .true., one needs to provide the inputfile. In the case without electron spin, the file is "neo_occ_electron.in". Otherwise "neo_occ_electron.1.in" and "neo_occ_electron.2.in" should be specified for each spin. The file should contains two columns with multiple rows. The first column indicates the electron states to update. The second column indicates the updated occupation number. For example,
              i_state_1 occ_1
              i_state_2 occ_2
              …
              i_state_N occ_N

Tag: NEO_rt_output_dipole(control.in)

Usage: NEO_rt_output_dipole flag
Purpose: Determines whether to write dipole to file (with surffix ".rt-tddft.dipole.dat") during RT-NEO-TDDFT.
flag is a logical string, either .false. or .true.. Default: .true..

Tag: NEO_rt_output_proton_center(control.in)

Usage: NEO_rt_output_proton_center flag1, flag2
Purpose: Determines whether to write proton center to file (with surffix ".rt-tddft.proton-center.*") during RT-NEO-TDDFT.
flag1, flag2 are logical string, either .false. or .true.. flag1 controls whether to write in ".dat" format. flag2 controls whether to write in ".xyz" format. Default: .false. for both.

Tag: NEO_rt_output_energy(control.in)

Usage: NEO_rt_output_energy flag
Purpose: Determines whether to write proton related energy to file (with surffix ".rt-tddft.energy.dat") during RT-NEO-TDDFT.
flag is a logical string, either .false. or .true.. Default: .true..

Tag: NEO_rt_output_n_step(control.in)

Usage: NEO_rt_output_n_step N
Purpose: Defines how often output files are written during RT-NEO-TDDFT simulations.
N: Integer. The output file are written every N electron steps. It is based on electron steps but not proton steps.
Default: Same as electron. See RT_TDDFT_propagation for how to control RT-TDDFT output.

Save and restart a RT-NEO-TDDFT simulation

In order to write an restart file or restart a RT-NEO-TDDFT calculation from a restart file, no additional inputs are required other than the regular RT-TDDFT Tags. See RT_TDDFT_restart_write and RT_TDDFT_restart_read for more instruction.
If RT_TDDFT_restart_write are set, NEO will write a restart file for the proton system whenever restart file is created for electron systems. The written NEO restart file will have the same name as the electronic one but with a prefix "neo."
If RT_TDDFT_restart_read are set, FHI-aims will run a NEO-SCF cycle, but the resulting eigenvectors, etc. will be overwritten from what is read in the restart file - based on this, the real-time simulation will be re-initialized. The NEO restart file should be named "neo.rt-tddft.save.in"

$\displaystyle L(t)$	$\displaystyle=\int\text{d}\bm{r}^{e}\int\text{d}\mathbf{k}\sum_{j}{\psi_{j% \mathbf{k}}^{e}}^{*}(\bm{r}^{e},t)\left[i\frac{\partial}{\partial t}+\frac{1}{% 2m^{e}}\nabla_{\bm{r}^{e}}^{2}\right]\psi_{j\mathbf{k}}^{e}(\bm{r}^{e},t)$	(3.124)
	$\displaystyle-\frac{1}{2}\int\text{d}\bm{r}^{e}\text{d}\bm{r}^{\prime e}\frac{% e^{2}}{\|\bm{r}^{e}-\bm{r}^{\prime e}\|}\rho^{e}(\bm{r}^{e},t)\rho^{e}(\bm{r}^{% \prime e},t)-E_{XC}^{e}[\rho^{e}]$
	$\displaystyle+\int\text{d}\bm{r}^{p}\sum_{j}{\psi_{j}^{p}}^{*}(\bm{r}^{p},t)% \left[i\frac{\partial}{\partial t}+\frac{1}{2M^{p}}\nabla_{\bm{r}^{p}}^{2}% \right]\psi_{j}^{p}(\bm{r}^{p},t)$
	$\displaystyle-\frac{1}{2}\int\text{d}\bm{r}^{p}\text{d}\bm{r}^{\prime p}\frac{% e^{2}}{\|\bm{r}^{p}-\bm{r}^{\prime p}\|}\rho^{p}(\bm{r}^{p},t)\rho^{p}(\bm{r}^{% \prime p},t)-E_{XC}^{p}[\rho^{p}]$
	$\displaystyle+\frac{1}{2}\int\text{d}\bm{r}^{e}\text{d}\bm{r}^{p}\frac{e^{2}}{% \|\bm{r}^{e}-\bm{r}^{p}\|}\rho^{e}(\bm{r}^{e},t)\rho^{p}(\bm{r}^{p},t)-E_{EPC}[% \rho^{e},\rho^{p}]$
	$\displaystyle+\sum_{I}^{N}\frac{1}{2}M_{I}\left[\frac{d}{dt}\bm{R}_{I}(t)% \right]^{2}-\sum_{I<J}\frac{Z_{I}Z_{J}e^{2}}{\|\bm{R}_{I}(t)-\bm{R}_{J}(t)\|}$
	$\displaystyle-\int\text{d}\bm{r}^{p}\rho^{p}(\bm{r}^{p},t)\sum_{I}^{N}\frac{Z_% {I}e^{2}}{\|\bm{r}^{p}-\bm{R}_{I}(t)\|}+\int\text{d}\bm{r}^{e}\rho^{e}(\bm{r}^{e% },t)\sum_{I}^{N}\frac{Z_{I}e^{2}}{\|\bm{r}^{e}-\bm{R}_{I}(t)\|}$