Geometry Optimization of an Isolated Water Molecule
In this section, we take an isolated H2O molecule as an example for a 'geometry optimization' with the hydrogen atom treated as quantum protons using the NEO-DFT method.
Geometry Optimization with NEO
For classical nuclei, "geometry optimization" refers to a directed search for a local minimum on the Born-Oppenheimer potential energy surface.
In NEO-DFT calculations, selected protons 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 electron and proton).
In the context of basis function centers for quantum nuclei, "geometry optimization" involves searching for the optimal positions of the quantum proton centers on a quantum potential energy surface contributed by both electrons and nuclei.
The details of computing energy gradients in NEO-DFT calculations are provided in the following paper:
J. Xu, et al., "Lagrangian Formulation of Nuclear-Electronic Orbital Ehrenfest Dynamics with Real-time TDDFT for Extended Periodic Systems". https://doi.org/10.48550/arXiv.2407.18842
Workflow
We are going to cover:
- Prepare the
geometry.infile: Set up a rough initial structure of a water molecule. - Prepare the
control.infile: Configure the parameters of the NEO-DFT calculation that requires energy gradient. A validcontrol.infor NEO-DFT is provided at the end. - Run the calculation: Execute the calculation with isotope effects.
- Check the results: Go through the output files for relaxed strcutres.
Prepare the geometry.in File
Even though it's important to start a calculation from a reasonable structure, the water molecule is simple enough to start from an educated guess.
A possible geometry input file geometry.in for an initial guess of the water molecule (not the final geometry) might be:
atom 0.00000 0.00000 0.00000 O
atom 0.00000 -1.00000 0.20000 H
atom 0.00000 1.00000 0.20000 H
Prepare the control.in File
The control.in file for geometry optimization with NEO-DFT similar to the one used for ground state calculation, with only a few extra parameters in the electron runtime choice.
However, performing geometry optimization requires evaluating energy gradients, which is not supported by the default RI method in the current NEO-DFT implemetation.
In this tutorial, we will discuss an alternative approach to enable energy gradient calculations.
Electron Input
First, the runtime choices should look like this:
xc pbe
relativistic none
relax_geometry bfgs 1.e-2
We use the PBE exchange correlation functional.
For the geometry relaxation, a variant of the bfgs algorithm (BFGS: Broyden–Fletcher–Goldfarb–Shanno; more specifically, a trust radius based version) are used.
A local minimum of the total energy is considered to be reached when the magnitude of all residual forces (energy gradients) on the atomic nuclei (or basis function centers for quantum nuclei) is smaller than a predefined threshold, i.e., \(1\cdot 10^{-2}\) eV/Å.
or a smaller threshold, the simulation can take longer to converge.
On the other hand, a larger threshold may be insufficient to find the actual (local) minimum structure.
Notes about Structure Relaxation in FHI-aims
The terms "local structure optimization" and "structure relaxation" are used interchangeably in the following. They mean the same thing.
- FHI-aims can calculate gradients of the total energy with respect to atomic positions ("forces") and lattice parameters ("stresses") for local-density, generalized-gradient and meta-generalized gradient density functionals, as well as hybrid density functionals. The later (stresses) are not supported by NEO-DFT.
- A variant of the Broyden–Fletcher–Goldfarb–Shanno (
bfgs) algorithm, enhanced by a trust-region method (trm), is employed in FHI-aims to find a local minimum of the potential energy surface. Such a minimum is considered to be found if the moduli of all forces (energy gradients) should be smaller than a small, configurable threshold value, typically smaller than \(5\cdot10^{-3}\) eV/Å. Within FHI-aims, the termsbfgsandtrmare used to denote the same algorithm. - It is highly worthwhile to understand at least the basics of the optimization algorithms mentioned above - follow the links.
Importantly, this shows that the
bfgsalgorithm relies on a good approximation of the Hessian matrix, that is, the second derivatives of the total energy with respect to all combinations of atomic and/or lattice coordinates for a given structure. A good initial guess for the Hessian matrix at the outset of a structure optimization can be immensely helpful to speed up the calculation. We use a variant of the initial guess suggested by Lindh and coworkers but if this does not work, a simple diagonal initial guess for the Hessian can also be used. - Even more importantly, the
bfgsalgorithm updates and improves its own guess of the Hessian matrix as the structure optimization progresses. Structure optimization is a step-wise process. At each optimization step, FHI-aims writes the current atomic/lattice geometry into a filegeometry.in.next_step. It also writes a filehessian.aims, which contains the current guess of the Hessian matrix determined by thebfgsalgorithm. - To start a new simulation with the relaxed structure,
geometry.in.next_stepshould be used as thegeometry.infile of the new simulation, andhessian.aimsshould be present in the new simulation folder. For example, one can copyold-simulation/geometry.in.next_stepintonew-simulation/geometry.in, prepare appropriatecontrol.in, and then start the new simulation.
NEO Input
As mentioned earlier, the default density fitting method (NEO_df_type set to even_tempered) does not support energy gradient calculations.
In this section, we will explore some other available options in control.in for NEO-DFT calculations, such as a new NEO_type and the aims option for NEO_density_fitting_type to enable energy gradients evaluatiions.
NEO_type 3
NEO_epc_type 17
NEO_basis_preset_type pb4d
NEO_scf_proton_max_step 20
# NEO_nuclei_type D
NEO_exchange_type noRI
# density fitting
NEO_df_type aims
NEO_df_l_max 6
# cut_off
NEO_proton_hartree_rcut 20
NEO Essential
In the first line, we set NEO_type to 3 which enables a NEO-DFT calculation.
The difference between NEO_type 1 and NEO_type 3 is:
NEO_type 1: FHI-aims first performs a conventional DFT calculation. Then it reinitializes and restarts a new NEO-DFT using the conventional DFT results as the starting point.NEO_type 3: The electron part and NEO part are initialized and run simultanously, without performing an initial conventional DFT calcuation as a starting point.
The rest of the parameters are set the same as previous tutorials. We use the epc17-2 as EPC functional and use the PB4-D preset basis set.
To study the isotope effect with NEO-DFT, you can use the NEO_nuclei_type keyword.
It offers two options H (default) and D, for hydrogen and deuterium, respectively.
Isotope Effect with NEO
The NEO_nuclei_type option only changes the mass of quantum nuclei in the NEO part of the NEO-DFT calculations.
It does not change the isotope property in the electron parts of the FHI-aims simulations.
For molecular dynamics simulations or any calculations that use the mass of nuclei, please also adjust the isotope tag in geometry.in accordingly.
Note: Isotope effect is a recently implemented feature and has not been thoroughly tested. Use this feature at your own risk!
NEO density fitting
To calculate forces for geometry optimization, the density fitting type NEO_df_type needs to be set to aims.
This option evaluates the electrostatic potential of quantum protons in the same manner as that of electrons in FHI-aims, using an real-space grid basis set.
For more information please refer to this paper:
V. Blum, et al., "Ab initio molecular simulations with numeric atom-centered orbitals", Comput. Phys. Commun. 180, 11 (2009). https://doi.org/10.1016/j.cpc.2009.06.022
With this option, we recommend to increase the value of NEO_df_l_max to match the l_hartree parameter in the species defaults to achieve better numerical accuracy.
By default, NEO_df_l_max is set to 3.
Since the aims density fitting method is not an RI method, we cannot calculate the Hartree-Fock exchange using this approach.
Setting NEO_exchange_type to noRI ensures that the exchange interaction of quantum protons is evaluated by the excat exchange in a short range.
Before running the relaxation, do not forget to attach the species defaults for O and H. In this tutorial, we will keep using the intermediate defaults NAO basis set.
A Valid control.in for NEO-DFT
xc pbe
relativistic none
relax_geometry bfgs 1.e-2
NEO_type 3
NEO_epc_type 17
NEO_exchange_type noRI
NEO_basis_preset_type pb4d
# density fitting
NEO_df_type 0
NEO_df_l_max 6
# cut_off
NEO_proton_hartree_rcut 20
# scf
NEO_scf_proton_max_step 20
# NEO_scf_accuracy_dm 1e-3
# NEO_scf_accuracy_eev 1e-3
NEO_scf_accuracy_rho 1e-4
################################################################################
#
# FHI-aims code project
# Volker Blum, 2017
#
# Suggested "intermediate" 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 1e-4
#
radial_base 24 7.0
radial_multiplier 2
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
################################################################################
#
# Suggested additional basis functions. For production calculations,
# uncomment them one after another (the most important basis functions are
# listed first).
#
# Basis constructed for dimers: 0.5 A, 0.7 A, 1.0 A, 1.5 A, 2.5 A
#
################################################################################
# "First tier" - improvements: -1014.90 meV to -62.69 meV
hydro 2 s 2.1
hydro 2 p 3.5
# "Second tier" - improvements: -12.89 meV to -1.83 meV
hydro 1 s 0.85
# hydro 2 p 3.7
# hydro 2 s 1.2
for_aux hydro 3 d 7
# "Third tier" - improvements: -0.25 meV to -0.12 meV
# hydro 4 f 11.2
# hydro 3 p 4.8
# hydro 4 d 9
# hydro 3 s 3.2
################################################################################
#
# FHI-aims code project
# Volker Blum, 2017
#
# Suggested "intermediate" defaults for O atom (to be pasted into control.in file)
#
################################################################################
species O
# global species definitions
nucleus 8
mass 15.9994
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 1e-4
#
radial_base 36 7.0
radial_multiplier 2
angular_grids specified
division 0.1817 50
division 0.3417 110
division 0.4949 194
division 0.6251 302
division 0.8014 434
# division 0.8507 590
# division 0.8762 770
# division 0.9023 974
# division 1.2339 1202
# outer_grid 974
outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 4.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 3.
################################################################################
#
# Suggested additional basis functions. For production calculations,
# uncomment them one after another (the most important basis functions are
# listed first).
#
# Constructed for dimers: 1.0 A, 1.208 A, 1.5 A, 2.0 A, 3.0 A
#
################################################################################
# "First tier" - improvements: -699.05 meV to -159.38 meV
hydro 2 p 1.8
hydro 3 d 7.6
hydro 3 s 6.4
# "Second tier" - improvements: -49.91 meV to -5.39 meV
hydro 4 f 11.6
# hydro 3 p 6.2
# hydro 3 d 5.6
for_aux hydro 5 g 17.6
# hydro 1 s 0.75
# "Third tier" - improvements: -2.83 meV to -0.50 meV
# ionic 2 p auto
# hydro 4 f 10.8
# hydro 4 d 4.7
# hydro 2 s 6.8
# "Fourth tier" - improvements: -0.40 meV to -0.12 meV
# hydro 3 p 5
# hydro 3 s 3.3
# hydro 5 g 15.6
# hydro 4 f 17.6
# hydro 4 d 14
# Further basis functions - -0.08 meV and below
# hydro 3 s 2.1
# hydro 4 d 11.6
# hydro 3 p 16
# hydro 2 s 17.2
Run the Calculation
In this tutorial, we will perform a single geometry optimization of the H2O moelcule using the control.in file provided above.
Optionally, you can test the isotope effect by setting NEO_nuclei_type to D (deuterium) and compare the differences in zero-point energy and structure at the ground state.
Check the Results
After the FHI-aims calculation is complete, your folder should contain:
Output Files
| File name | Description |
|---|---|
H2O.out |
FHI-aims output stream |
geometry.in.next_step |
Geometry file after relaxation |
hessian.aims |
BFGS Hessian matrix |
The relaxed structure is written to geometry.in.next_step and looks like this:
#
# This is the geometry file that corresponds to the current relaxation step.
# If you do not want this file to be written, set the "write_restart_geometry" flag to .false.
# aims_uuid : EC03F690-7193-4760-A415-979445838DB1
#
atom 0.00000000 0.00000000 -0.25619363 O
atom 0.00000000 -0.78016475 0.32809682 H
atom 0.00000000 0.78016475 0.32809682 H
#
# What follows is the current estimated Hessian matrix constructed by the BFGS algorithm.
# This is NOT the true Hessian matrix of the system.
# If you do not want this information here, switch it off using the "hessian_to_restart_geometry" keyword.
#
trust_radius 0.2000000030
hessian_file
With the NEO-DFT relaxation, our water structure goes from
to
A file hessian.aims is also written.
The geometry.in.next_step and hessian.aims files can be used to continue the relaxation with tighter settings of species defaults to obtain better-converged results.
Restart a Geometry Optimization
The files geometry.in.next_step and hessian.aims are the two key pieces needed to restart a structure relaxation from the results of an earlier relaxation.
Such a restart can be desirable for several reasons, e.g, in order to:
- Continue an unfinished relaxation that has stopped (for example, due to a queue wall time limit or because the earlier FHI-aims relaxation run encountered a problem and stopped with an error message).
- Use the data created by a pre-relaxation run using computationally cheaper "light" settings and follow up with a post-relaxation using more accurate "intermediate" or "tight" settings.
- Use the data generated by pre-relaxing with a cheaper density functional (e.g., PBE) as a starting point to initialize a post-relaxation with a more expensive density functional (e.g., PBE0).
The method to restart a structure relaxation from an earlier one is simple:
- Run the pre-relaxation, which generates intermediate files
geometry.in.next_stepandhessian.aims - After the pre-relaxation is finished, create a new directory for the followup relaxation
- Copy the
geometry.in.next_stepandhessian.aimsfiles from the pre-relaxation directory to the followup-relaxation directory. - Change the name of the copy of
geometry.in.next_steptogeometry.in. - Create a new
control.infile in the followup-relaxation directory - Rerun FHI-aims in the follow-up directory.
We have completed the NEO-DFT tutorial!
Solutions
You can find all the solutions to the exercises above by clicking on the button below: