Ground State Simulation of an Isolated FHF− Molecule
In this section, we will simulate the ground state of an isolated FHF− molecule using the NEO-DFT method, treating the hydrogen atom as a quantum proton. We will also plot the proton density and proton wavefunctions in a cube file.
For a deeper understanding of the NEO-DFT calculations performed in this tutorial, please refer to the following paper:
J. Xu, et al., "Nuclear–electronic orbital approach to quantization of protons in periodic electronic structure calculations", J. Chem. Phys. 156, 224111 (2022). https://doi.org/10.1063/5.0088427
Workflow
We are going to cover:
- Prepare the
geometry.infile: Set up the structure of the FHF- molecule. - Prepare the
control.infile: Configure the important parameters of the NEO-DFT calculation, including functionals, convergence criteria, basis sets, etc. A validcontrol.infile for NEO-DFT is provided at the end. - Run the calculation: Execute the calculation using multiple processors.
- Check the results: Review the output files for energies, positions, and densities.
Prepare the geometry.in File
Quantum nuclei / Quantum protons
In current implementations, only hydrogen atoms can be treated with NEO. The terms quantum nuclei and quantum protons are used interchangeably.
As with all calculations, it's important to start from a reasonable structure. Rather than performing a geometry optimization here, we use a pre-optimized structure of the FHF− molecule from existing numerical studies. If you are not familiar with how to obtain a reasonable initial molecular structure, please refer to the Basics of Running AIMS tutorial for details on geometry optimization.
An example of the geometry input file geometry.in is:
atom 0.0 0.0 -1.16407 F
atom 0.0 0.0 1.16407 F
atom 0.0 0.0 0.0 H
Positions of Quantum Nuclei
In NEO-DFT calculations, selected nuclei are treated as quantum particles, which do not have well-defined positions.
The corresponding coordinates in the geometry.in file are only used as the basis function centers (of both electrons and protons).
Generally speaking, the optimized structure from conventional DFT is a reasonable choice for this purpose. One can also do geometry optimization with NEO-DFT, which will be covered in Section 3.
Using an ASCII text editor, such as vim or emacs, allows us to work directly with the input files (which are ASCII files).
The syntax of the geometry.in file is pretty much self-explanatory.
A more detailed definition of the format can be found here or in the FHI-aims manual.
Prepare the control.in File
Our next step is to define the parameters that control the various aspects of a simulation with FHI-aims.
These parameters are provided in a file called control.in and should be located in the same folder as geometry.in.
Our objective in this tutorial is to carry out a single-point NEO-DFT calculation of the FHF- molecule.
For NEO-DFT calculations, we need to define parameters for electrons and quantum nuclei separately.
Electron Input
As in a conventional DFT calculation with FHI-aims, the electron part contains two sets of parameters: runtime choices and numerical settings. We will provide a brief overview of the electron input settings. For detailed instructions on configuring parameters for electrons, please refer to the Basics of Running AIMS tutorial.
For the runtime choices, we use the PBE exchange-correlation functional for this negatively charged system (with one addtional electron):
xc pbe
charge -1
For the numerical settings, we use the cc-pVDZ Gaussian basis sets for the H and F atoms to compare with other Gaussian quantum chemistry packages.
One can also use the atom-centered orbital (NAO) basis sets and numerical definitions provided in the species_defaults directory.
Species Defaults: Choice of Numerical Settings, Including the Basis Set
FHI-aims uses numerically tabulated NAO basis sets to represent the quantum mechanical orbitals and wave functions. These basis functions are specific to each chemical element (`species``) used in the simulation.
FHI-aims also employs a variety of other numerical parameters that are specific to each atom and/or chemical element, such as angular-momentum expansion order of the electrostatic (Hartree) potential and atom-centered real-space integration grids, etc.
Predefined files that contain these numerical definitions for each chemical element are located in species_defaults directory of the FHI-aims code distribution.
At least four different levels of precision are provided for each chemical element: light, intermediate, tight, and really_tight.
A more accurate choice of defaults increases the computational cost of the simulation.
To set the numerical settings, the desired Species Defaults files for all atoms needs to be appended to control.in.
You can excute the following commands:
cat [FHI-aims-directory]/species_defaults/non-standard/gaussian_tight_770/cc-pVDZ/01_H_default >> control.in
cat [FHI-aims-directory]/species_defaults/non-standard/gaussian_tight_770/cc-pVDZ/09_F_default >> control.in
[FHI-aims-directory] should be replaced with the proper path to your FHI-aims distribution on the computer to be used.
NEO Input - Essential
We also need to define the NEO runtime choices in control.in after the electron runtime choices.
All NEO-related parameters have the prefix NEO_/
NEO_type 1
NEO_epc_type 17
NEO_basis_preset_type pb4d
# NEO_basis
# H 3
# NEO_basis_end
In the first line NEO_type is set to 1 which enables the NEO-DFT calculation in FHI-aims.
This is the only essential parameter required to run NEO-DFT calculations.
By default, the NEO_type is set to 0, meaning no NEO calculations are not invoked.
The NEO_epc_type parameter defines the electron-proton correlation (EPC) functional used in the NEO-DFT simulations.
Here we use a local density approximation (LDA) type of functional: epc17-2.
By default,NEO_epc_type is set to 0 meaning no EPC functional is used.
Electron-Proton Correlation Functional
Electron-proton correlation functional describes the correlation between quantum electrons and quantum protons as a functional of electron density and proton density. It is shown to be important for accurately describing the zero-point energy with NEO-DFT.
Several electron−proton correlation functionals based on a multicomponent extension of the Colle−Salvetti formalism have been developed:
Y. Yang, et al., "Development of a practical multicomponent density functional for electron-proton correlation to produce accurate proton densities", J. Chem. Phys. 147, 114113 (2017). https://doi.org/10.1063/1.4996038
The epc17 functional is defined as:
\(E^c_\text{epc17}[\rho^p,\rho^e]=-\int \text d\boldsymbol r \frac{\rho^p(\boldsymbol r)\rho^e(\boldsymbol r)}{a-b\rho^p(\boldsymbol r)^{\frac{1}{2}}\rho^e(\boldsymbol r)^{\frac{1}{2}}+c\rho^p(\boldsymbol r)\rho^e(\boldsymbol r)}\)
In the current implementation, NEO_epc_type supports 17, 17-2 and 0.
The NEO_basis_present_type parameter sets the basis function set for quantum protons.
Here we used the pb4d proton basis set, which includes four s orbitals, three p orbtials and 2d orbitals.
By default, NEO_basis_present_type is set to pb4f2.
Proton Basis Functions
The electronic wave functions are expressed as linear combinations of grid-based NAOs the FHI-aims. The accurate description of the sharp electron peak close to the nucleus is a distinct advantage of the NAO approach over other types of basis sets.
For quantum protons, the cusp at the center is less steep, and primitive Gaussian type orbital (GTO) basis sets are commonly chosen to describe the proton wave functions. Several series of GTO nuclear basis sets with varying levels of angular momentum, denoted PB4, PB5, and PB6, have been developed using machine-learning optimization procedures:
Q. Yu, et al., "Development of nuclear basis sets for multicomponent quantum chemistry methods", J. Chem. Phys. 152, 244123 (2020). https://doi.org/10.1063/5.0009233
WIth current implementation, NEO_basis_present_type support pb4d, pb4f2 and pb5d.
In this tutorial, the pb4d proton basis set,
Instead of using NEO_basis_present_type, one can also define custom GTO basis sets for each quantum proton with the NEO_basis and NEO_basis_end parameters.
More details about this option will be covered in the next section
The NEO_basis and NEO_basis_end parameters are used to specify which hydrogen atoms are treated as quantum protons.
The indices correspond to the order of atoms in geometry.in.
If not provided, all the hydrogen atoms in geometry.in are treated with NEO.
NEO Input - Extra
Other optional parameters in control.in for controlling the NEO-DFT calculations are listed below.
# numerical cut_off s
# NEO_proton_hartree_rcut 20
# NEO_proton_hartree_rclassic 5
# scf convergence threshold
NEO_scf_proton_max_step 20
# NEO_scf_accuracy_dm 1e-5
# NEO_scf_accuracy_eev 1e-5
# NEO_scf_accuracy_rho 1e-6
# for plotting
# define cube center ("c","cen","center") in Angstrom
NEO_plot_cube c 0.0 0.0 0.0
# define cube origin ("o","org","origin") in Angstrom
#NEO_plot_cube o 0.0 0.0 0.0
# define cube axis
NEO_plot_cube ax1 1.0 0.0 0.0
NEO_plot_cube ax2 0.0 1.0 0.0
NEO_plot_cube ax3 0.0 0.0 2.0
# define cube nsteps xyz
NEO_plot_cube ns 20 20 40
# include only quantum proton in cube "qp","quantum_proton"
#NEO_plot_cube qp
# print densiy (d,df,density)
#NEO_plot_cube df outfilename
NEO_plot_cube df
# print densiy (wf,wave,wavefunction)
# NEO_plot_cube wf i_begin i_end outfilename
NEO_plot_cube wf
In the runtime part, the first part specifies numerical cutoff values (in Å) for speeding up the calculations. Using the default values should suffice. One can increasing these values for more accurate results at the cost of increased computational effort, or decrease them for faster calculations with potentially reduced accuracy.
Electrostatic Potential in NEO
the ES potential of a given proton at basis fucntion center is calculated seperately for different ranges:
- short-range interactions (within
NEO_proton_hartree_rclassic) are calculated using analytical formulas with density fitting. - mid-range interactions are determined by multipole formulas, where the charge distribution are considered as set of multipoles.
- long-range interactions (beyond
NEO_proton_hartree_rcut) are treated as those from a point charge (for isolated molecule) or computed via Ewald summation (with periodic boundary conditions).
The second section contains the parameters for NEO-SCF convergence, which are similar to their electron conterparts.
The NEO_scf_proton_max_step parameter controls the maximum number of proton SCF iterations within a single electron SCF step.
By default, NEO_scf_proton_max_step is set to 500 for safety.
We change it to 20 because it is not necessary for the proton SCF to converge in each individual electron step.
The output section controls the printing of cube files.
Comments that begins with # explains the meaning of each parameter.
In this tutorial, we will print two cube files: one fro the proton density and one for the proton wave function.
Full documentation on the input parameters can be found in the FHI-aims manual.
A Valid control.in for NEO-DFT
Your final control.in file should look as follows:
# runtime choices for electrons
xc pbe
charge -1
# runtime choices for NEO
NEO_type 1
NEO_epc_type 17
NEO_basis_preset_type pb4d
NEO_scf_proton_max_step 20
# for plotting
# define cube center ("c","cen","center") in Angstrom
NEO_plot_cube c 0.0 0.0 0.0
# define cube axis
NEO_plot_cube ax1 1.0 0.0 0.0
NEO_plot_cube ax2 0.0 1.0 0.0
NEO_plot_cube ax3 0.0 0.0 2.0
# define cube nsteps xyz
NEO_plot_cube ns 20 20 40
# print densiy (d,df,density)
NEO_plot_cube df
# print densiy (wf,wave,wavefunction)
NEO_plot_cube wf
################################################################################
#
# FHI-aims code project
# Igor Ying Zhang and Xinguo Ren, Fritz Haber Institute Berlin, 2012
#
# Suggested "cc-pVDZ" defaults for H atom (to be pasted into control.in file)
#
################################################################################
species H
# global species definitions
nucleus 1
mass 1.00794
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 0e-0
#
radial_base 24 7.0
radial_multiplier 6
angular_grids specified
division 0.1930 50
division 0.3175 110
division 0.4293 194
division 0.5066 302
division 0.5626 434
division 0.5922 590
# division 0.6227 974
# division 0.6868 1202
outer_grid 770
# outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 1 s 1.
# ion occupancy
ion_occ 1 s 0.5
################################################################################
#
# For exact comparison with all GTO-based codes, one or more of
# the following flags are needed:
#
include_min_basis .false.
pure_gauss .true.
#
# H cc-pVDZ
gaussian 0 3
13.0100000 0.0196850
1.9620000 0.1379770
0.4446000 0.4781480
gaussian 0 1 0.1220000
gaussian 1 1 0.7270000
################################################################################
#
# FHI-aims code project
# Igor Ying Zhang and Xinguo Ren, Fritz Haber Institute Berlin, 2012
#
# Suggested "cc-pVDZ" defaults for F atom (to be pasted into control.in file)
#
################################################################################
species F
# global species definitions
nucleus 9
mass 18.9984032
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 0e-0
#
radial_base 37 7.0
radial_multiplier 6
angular_grids specified
division 0.4014 110
division 0.5291 194
division 0.6019 302
division 0.6814 434
division 0.7989 590
# division 0.8965 770
# division 1.3427 974
# outer_grid 974
outer_grid 770
# outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 5.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 4.
################################################################################
#
# For exact comparison with all GTO-based codes, one or more of
# the following flags are needed:
#
include_min_basis .false.
pure_gauss .true.
#
# F cc-pVDZ
gaussian 0 8
14710.0000000 0.0007210
2207.0000000 0.0055530
502.8000000 0.0282670
142.6000000 0.1064440
46.4700000 0.2868140
16.7000000 0.4486410
6.3560000 0.2647610
1.3160000 0.0153330
gaussian 0 8
14710.0000000 -0.0001650
2207.0000000 -0.0013080
502.8000000 -0.0064950
142.6000000 -0.0266910
46.4700000 -0.0736900
16.7000000 -0.1707760
6.3560000 -0.1123270
1.3160000 0.5628140
gaussian 0 1 0.3897000
gaussian 1 3
22.6700000 0.0448780
4.9770000 0.2357180
1.3470000 0.5085210
gaussian 1 1 0.3471000
gaussian 2 1 1.6400000
Run the Calculation
To run a parallel simulation with N processes, you can excuate the following command:
mpirun -n N aims.x > FHF.out 2>&1
mpirun -n N aims.x | tee FHF.out
In these commands, we write the output of the simulation to a file named FHF.out file, but you can choose any other name that you prefer.
As a comparison, we will also perform a conventional DFT calculation by changing NEO_type 1 to NEO_type 0 in control.in.
Save the output file of this run to a file named FHF-noNEO.out
Note on Running
The binary name aims.x should be replaced with the actual name of the FHI-aims executable you compiled, including the corresponding path( i.e., the directory where the executable is located).
For example, on the writer's laptop, the current binary name is ~/Software/fhi-aims/bin/aims.240507.scalapack.mpi.x.
The mpirun command facilitates the parallel execution of the code and may have different names that depend on the particular computer system and/or queueing system used.
The actual name and usage are usually documented by the computing center in question.
Check the Results
After both the NEO-DFT and DFT FHI-aims calculations are complete, your folder should contain:
Output Files
| File name | Description |
|---|---|
FHF.out |
FHI-aims output stream of NEO-DFT |
FHF-noNEO.out |
FHI-aims output stream of DFT |
proton_density.cube |
Cube file with proton density |
proton_wf001.cube |
Cube file with proton wavefunction |
Output Stream
The FHI-aims output stream of NEO-DFT begins with some important information regarding how the calculation was executed.
It includes complete copies of control.in and geometry.in, details all self-consistent field (SCF) iterations, as well as NEO-DFT SCF iterations.
The aims.out ends with the line:
Have a nice day.
which is used as FHI-aims' indicator that the calculation converged as expected and has completed successfully.
Near the end of the standard output FHF.out, you can find the final total energy from the s.c.f. cycle.
You can estimate the zero-point energy of the quantum proton by comparing the difference in total energy between DFT (-5452.960 eV) and NEO-DFT (-5452.818 eV) runs.
Just before that, the expectation value of the position operator of quantum protons and their spread is displayed.
In our case, we got:
Writing Proton expectation value x,y,z and variance
Block separation
H -0.0000000000 -0.0000000000 -0.0000000000 0.0244473093 1 3
Eigenstate separation
H -0.0000000000 -0.0000000000 -0.0000000000 0.0244473093 1 -33.6753716546eV, Occ: 1.00000
3 is centered exactly at the origin with a spread of 0.024 Å\(^2\) and an eigenvalue of -33.675 eV.
This result can provide more information when you have more than one proton in the system.
Cube files
The two cube files include the spacial distribution of proton density and proton wave functions.
You can open the .cube file using the visualization software of your choice - e.g VMD, VESTA, and Avogadro, to just name a few free software packages.
The proton_density.cube should looks like:
In this image, the grey spheres represent the fluoride ions, the pink sphere is the quantum proton basis center, and the yellow isosurface represents the quantum proton density.
This concludes the first part of the tutorial!!!
Solutions
You can find all the solutions to the above exercises by clicking on the button below: