Ground State Simulation of trans-Polyacetylene
This tutorial introduces the basic concepts of NEO-DFT calculations for periodic system geometries using FHI-aims. In this section, we will carry out a ground state calculation of a trans-polyacetylene chain with both hydrogen atoms treated as quantum protons using the NEO-DFT method. We will also compare the electron density of states (DOS) from conventional DFT and NEO-DFT calculations.
More details about the periodic NEO-DFT calculations performed in this tutorial are provided in 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 trans-polyacetylene chain with periodic boundary conditions (PBC). - Prepare the
control.infile: Introduce additional parameters for the NEO-DFT calculation, including density fitting, initial guess, etc. A validcontrol.infor NEO-DFT is provided at the end. - Run the calculation: Run the calculations with DFT, NEO-DFT, with and without the EPC functional.
- Check the results: Go through the output files for energies, positions, DOS.
Prepare the geometry.in File
The syntax for the geometry.in file remains the same as in the previous section, with only a few additions:
lattice_vector 2.469 0.0 0.0
lattice_vector 0.0 50.0 0.0
lattice_vector 0.0 0.0 50.0
atom -0.58442 -0.32662 0.0 C
atom 0.58442 0.32662 0.0 C
atom 0.58442 1.40362 0.0 H
atom -0.58442 -1.40362 0.0 H
The first three lines starting with lattice_vector define the three unit cell vectors \(a_i\) (\(i\)=1,2,3).
In our example (a one-dimensional chain with PBC), we use the experimental cell parameter of the [C2H2]n along the chain direction.
In the other two directions, we have a sufficiently large vacuum to prevent interactions between periodic images.
The following lines define the positions and species of the atoms in this system.
Fractional Coordinates
atom and atom_frac
One can only use either the atom or the atom_frac keyword in a geometry.in file.
Mixing these keyword in the same geometry.in file is not supported.
The position of atoms can also be defined in so-called fractional coordinates (atom_frac) \(f_1\), \(f_2\), \(f_3\).
Fractional coordinates are especially convenient when working with high-symmetry crystals.
The Cartesian coordinates \(x\), \(y\), \(z\) of each atom follow by multiplying fractional coordinates and lattice vectors:
\((x,y,z)= f_1 * a_1 + f_2 * a_2 + f_3 * a_3\).
Practical Recommendation
Visual verification of the structure is usually key to avoid simple mistakes. This is even more important for periodic systems, as nearest neighbors and atom distances cannot be identified easily without visualization.
Prepare the control.in File
The control.in file for NEO-DFT in a periodic system is similar to the one used for a non-periodic system.
In particular, the NEO input can remain unchanged from the molecular case.
In the NEO section, we will explore additional options including density fitting, initial guess strategies, and the selection of quantum protons.
Electron Input
First, we specify the runtime choices, which should look like this:
xc pbe
relativistic none
k_grid 100 1 1
output dos -20 10 300 0.1
We are aiming for a ground state NEO-DFT calculation using the PBE XC functional.
Additionally, in case of periodic systems, we need to specify the grid of crystal momentum (\(k\)) vectors used to sample the reciprocal lattice of the system. f
In FHI-aims, this \(k\)-space integration grid is specified by the keyword k_grid, followed by three integer numbers.
These numbers specify the number of k-points along each reciprocal lattice vector.
Their order corresponds to the order of lattice vectors specified in the geometry.in file.
FHI-aims will create a three-dimensional, evenly spaced, \(\Gamma\)-centered \(k\)-grid using the number of points specified for each reciprocal space direction.
\(k\)-Grid Density
Alternatively, the keyword k_grid_density followed by a single floating point number can be used. FHI-aims will again create again an even-spaced, \(\Gamma\)-centered k-grid, with the same density along all reciprocal lattice vectors.
There is no absolute recommendation regarding whether to use k_grid or k_grid_density.
- In general, the keyword
k_grid_densityis a practical option if one employs a sufficiently large k-point density. To obtain well-converged calculations (with respect to the k-grid), it is often preferable to fulfill the criteria \(n_i * a_i > 50 Å\), for the lattice vector length \(a_i\) and the number of k-points \(n_i\) along the corresponding \(k\)-space direction \(i\). Usingk_grid_densityprevents the accidental use of over-converged settings, e.g., for larger supercell or slab calculations, where usually a much smaller number of k-points along the longer lattice vectors is sufficient. - To guarantee consistent k-grid numbers between different systems, e.g., when comparing total energies between different supercell sizes, it is preferable to use the
k_gridkeyword, which allows one to exactly control the number of k-points along each \(k\)-space direction. For example, if one doubles the unit cell size in each direction, you can divide the number of k-points in all directions by a factor of 2.
The necessary \(k\)-space grid can vary based on the target physical observable; for example, optical properties or densities of states generally require denser k_grid settings than just the total energy.
The k_grid 100 1 1 used in our current example, indicates a 100x1x1 k-grid, is chosen to achieve a converged and smooth DOS.
For printing the DOS of the system, the general syntax is:
output dos e_min e_max n_points broadening
e-min and e-max, float), the third parameter (n_points, an integer) specifies the number of points within this interval, and the last parameter (broadening, a float) specifies the broadening in eV.
In our example the broadening parameter is 0.1 eV.
If the goal is to resolve the individual peaks in the DOS, you can experiment with smaller values such as 0.05 eV.
Different Methods for DOS
Here, the \(k\)-space integration for calculating the DOS is performed using the Gaussian broadening method, as indicated by dos.
For smaller unit cells, we typically recommend the tetrahedron method using dos_tetrahedron, which can resolve the individual peaks of the DOS more accurately for higher \(k\)-point densities.
Generally, most Post-Processing tools that work with DFT should also work with NEO-DFT. In this tutorial, we will only cover the DOS of the trans-polyacetylene in this tutorial. You can also explore other options like:
Band Structure
The general syntax is:
output band k_xstart k_ystart k_zstart k_xend k_yend k_zend n_points symbol_start symbol_end
control.in) and the number of points for the segment.
Optionally, one can specify the name of the start and end point at the end of the corresponding line.
Mulliken-projected DOS and Band Structure
The general syntax is:
output species_proj_dos e_min e_max n_points broadening
output band_mulliken k_xstart k_ystart k_zstart k_xend k_yend k_zend n_points symbol_start symbol_end
NEO Input
In this section, we will explore some additional options that could be helpful when your calculation does not coverge or has numerical issues.
NEO_type 1
NEO_basis_type primitive
NEO_epc_type 17
# NEO_initial_guess_type potential
NEO_scf_proton_max_step 20
# NEO_df_type ET
NEO_df_n 12
# NEO_df_l_max 3
NEO_proton_hartree_rcut 20
# NEO_use_ewald F
# 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 2.5 0.0 0.0
NEO_plot_cube ax2 0.0 5.0 0.0
NEO_plot_cube ax3 0.0 0.0 1.0
# define cube nsteps xyz
NEO_plot_cube ns 50 100 20
# 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
NEO_basis
# PB4-D 4s3p2d
H 3
s 1.9570 1.000
s 8.7340 1.000
s 16.010 1.000
s 31.997 1.000
p 9.4380 1.000
p 13.795 1.000
p 24.028 1.000
d 10.524 1.000
d 19.016 1.000
H 4
s 1.9570 1.000
s 8.7340 1.000
s 16.010 1.000
s 31.997 1.000
p 9.4380 1.000
p 13.795 1.000
p 24.028 1.000
d 10.524 1.000
d 19.016 1.000
NEO_basis_end
NEO Essential
To initiate NEO-DFT calculations, we set NEO_type to 1.
Rather than utilizing a preset proton basis set, we use the primitive option as NEO_basis_type, which allow us to define our own choice of NEO basis configuration in the NEO_basis section.
The NEO_epc_type can be set to either 17 or 0 for assessing the impact of the EPC functional.
Adjusting Initial Guesses to Avoid Local Minima
In certain instances, NEO-DFT calculations may terminate prematurely without reaching the ground state. Notably, this can occur in the absence of EPC during this tutorial.
To address this issue, you can modify the initial guess for the NEO-SCF calculation by setting NEO_initial_guess_type to potential. This option configures the initial proton density guess through solving a single particle Hamiltonian.
By default, NEO_initial_guess_type is set to random, which employs the first basis function for the initial guess.
Therefore, one can also change the initial guess by changing the order of the proton basis functions within NEO_basis.
NEO density fitting
The NEO_df_type parameter controls the density fitting, also know as Resolution of idensity (RI).
By default, it is set to even_tempered or ET, which uses an even-tempered spherical Guassian auxiliary basis set.
Its radial part is defined as \(\psi_{i,l}=C_i\cdot r^l\cdot e^{-\alpha\cdot\beta^{i-1}r^2}\), where:
- \(i\) ranges from \(1\) to
NEO_df_n, - \(l\) ranges from \(0\) to
NEO_df_l_max, - \(C_i\) is a normalized constant,
- \(\alpha\) and \(\beta\) are defined by
NEO_df_alphaandNEO_df_beta.
In this tutorial, we increase NEO_df_n from 10 to 12 for better numerical accuracy.
The NEO_df_l_max (\(3\)), NEO_df_alpha (\(2\sqrt2\)), and NEO_df_beta (\(\sqrt2\)) are kept as default.
Density Fitting in NEO
To evaluating eletrostatic (ES) interactions on grid points in real space, NEO-DFT employs a density fitting method, where proton density is approximated using a chosen auxiliary basis.
The short range ES interactions (within NEO_proton_hartree_rclassic) are then calculated with the density fitting approach using analytical formulas.
By default, the proton density are fit in to a set of even tempered Gaussian auxiliary basis functions. This method also allows for the calculation of the Hartree-Fock exchange. For more details please refer to our NEO-DFT paper.
One can also chose to use the same density fitting approach, as used in FHI-aims for electrons by setting NEO_df_type to aims.
Currently, this is the only option to evaluate energy gradients.
We will cover more in the next section
For systems with periodic boudary conditions, the Ewald method are used to calculate long-range ES interactions on all grid points in real space.
You may choose to turn off the Ewald method by setting NEO_use_ewald to .false. and increasing NEO_proton_hartree_rcut for a brute force ES potential evaluation within the defined range.
NEO basis
The NEO basis set for individual protons can be define in the NEO basis section, specified with NEO_basis and NEO_basis_end.
Define NEO Basis
The general syntax is
NEO_basis
H index
[ l_1 alpha_1 weight_1 ]
[ l_2 alpha_2 weight_2 ]
[ ... ]
[ l_N alpha_N weight_N ]
...
NEO_basis_end
indexis a integer representing the index of hydrogen atom (proton) declared by theatomkeyword in geometry.in. Optionally, negtive number ofindexrepresents for adding a ghost center that do not have charge distributions, which usually works together withemptyin geometry.in. We will have more discussion on this in the RT-NEO-TDDFT section.l_iis a character specifying the angular momentum (s,p,d,f, ...) of the \(i\)th Gaussian-type orbital.alpha_iis a real number (in bohr\(^{−2}\)) specifying the exponent \(\alpha_i\) of the \(i\)th Gaussian-type orbital.weight_iis a real number specifying the weight of the \(i\)th gaussian type orbital. For primitive gaussian basis, weight should be 1.0.
In this tutorial, we continue using the PB4-D basis set for both hydrogen atoms in the system.
The above discussions only cover some of the optional input parameters. Full documentation on the input can be found in the manual of FHI-aims.
Before running the relaxation, do not forget to attach the species defaults for C and H. In this tutorial, we use a intermediate defaults NAO basis.
A Valid control.in for NEO-DFT
xc pbe
relativistic none
k_grid 100 1 1
output dos -20 10 300 0.1
NEO_type 1
NEO_basis_type primitive
NEO_epc_type 17
# NEO_initial_guess_type potential
# density fitting
NEO_df_n 12
NEO_df_l_max 3
# cut_off
NEO_proton_hartree_rcut 20
# scf
NEO_scf_proton_max_step 20
# 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 2.5 0.0 0.0
NEO_plot_cube ax2 0.0 5.0 0.0
NEO_plot_cube ax3 0.0 0.0 1.0
# define cube nsteps xyz
NEO_plot_cube ns 50 100 20
# 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
NEO_basis
# PB4-D 4s3p2d
# pos in angstrom
H 3
s 1.9570 1.000
s 8.7340 1.000
s 16.010 1.000
s 31.997 1.000
p 9.4380 1.000
p 13.795 1.000
p 24.028 1.000
d 10.524 1.000
d 19.016 1.000
H 4
s 1.9570 1.000
s 8.7340 1.000
s 16.010 1.000
s 31.997 1.000
p 9.4380 1.000
p 13.795 1.000
p 24.028 1.000
d 10.524 1.000
d 19.016 1.000
NEO_basis_end
################################################################################
#
# 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 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
################################################################################
#
# 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 C atom (to be pasted into control.in file)
#
################################################################################
species C
# global species definitions
nucleus 6
mass 12.0107
#
l_hartree 6
#
cut_pot 4.0 2.0 1.0
basis_dep_cutoff 1e-4
#
radial_base 34 7.0
radial_multiplier 2
angular_grids specified
division 0.2187 50
division 0.4416 110
division 0.6335 194
division 0.7727 302
division 0.8772 434
# division 0.9334 590
# division 0.9924 770
# division 1.0230 974
# division 1.5020 1202
# outer_grid 974
outer_grid 434
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 2.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 1.
################################################################################
#
# 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.25 A, 1.5 A, 2.0 A, 3.0 A
#
################################################################################
# "First tier" - improvements: -1214.57 meV to -155.61 meV
hydro 2 p 1.7
hydro 3 d 6
hydro 2 s 4.9
# "Second tier" - improvements: -67.75 meV to -5.23 meV
hydro 4 f 9.8
# hydro 3 p 5.2
# hydro 3 s 4.3
for_aux hydro 5 g 14.4
# hydro 3 d 6.2
# "Third tier" - improvements: -2.43 meV to -0.60 meV
# hydro 2 p 5.6
# hydro 2 s 1.4
# hydro 3 d 4.9
# hydro 4 f 11.2
# "Fourth tier" - improvements: -0.39 meV to -0.18 meV
# hydro 2 p 2.1
# hydro 5 g 16.4
# hydro 4 d 13.2
# hydro 3 s 13.6
# hydro 4 f 17.6
# Further basis functions - improvements: -0.08 meV and below
# hydro 3 s 2
# hydro 3 p 6
# hydro 4 d 20
Run the Calculation
Once the input files are ready, execute the calculation.
This time we will do three separate calculations, by modifing the control.in slightly:
-
NEO-DFT with epc17-2:
- Write output to
C2H2_EPC.out
- Write output to
-
NEO-DFT without EPC:
- Create a new folder
- Copy the input and change
NEO_epc_typeto0 - Write output to
C2H2_noEPC.out
-
Conventional DFT:
- Create a new folder
- Copy the input file and change
NEO_typeto0 - Write output to
C2H2_noNEO.out.
Avoid overwritting
One calculation - one directory.
On the computer intended to run FHI-aims, you will need to manage files using a command line interface and in different folders ("directories") intended to organize your data.
We strongly recommend to create a new directory (with its own control.in and geometry.in files) for every new FHI-aims calculation.
Rerunning and/or continuing FHI-aims calculations in the same directory that was used before will overwrite output and input files. In any simulation (including non-trivial failed attempts), keeping the data can be essential to later understanding and/or reconstructing what happened in a particular calculation.
Check the Results
After all FHI-aims calculations are complete, your folders should contain:
Output Files
| File name | Description |
|---|---|
C2H2_*.out |
FHI-aims output stream |
KS_DOS_total_raw.dat |
Raw total DOS data |
KS_DOS_total.dat |
DOS shifted by the chemical energy of the electrons (Fermi level) |
proton_density.cube |
Cube file with proton density |
Output Stream
Similar to the molecular case, you can find the final total energy from the s.c.f. cycle near the end of the standard output C2H2_*.out.
To estimate the zero-point energy of the Quantum proton, compare the total energy differences among the various runs:
- DFT(
-2104.419eV), - NEO-DFT epc17-2 (
-2103.530eV), - NEO-DFT without EPC (
-2090.953eV).
The absence of EPC leads to a drastic overestimation of the zero-point energy.
The expectation value of the position operator for the quantum protons and their spread are displayed just before the total energy. In the NEO-DFT with epc17-2 case, we obtained:
Writing Proton expectation value x,y,z and variance
Block separation
H 0.5857038149 1.4364634027 -0.0000000234 0.0198753978 1 3
H -0.5857038149 -1.4364634027 0.0000000234 0.0198753978 2 4
Eigenstate separation
H 0.5534112607 1.3572645463 -0.0000000221 0.2779203917 1 -27.2074400352eV, Occ: 1.00000
H -0.5534112607 -1.3572645463 0.0000000221 0.2779203917 2 -27.2074400351eV, Occ: 1.00000
3 and 4 are degenerate, each with eigenvalues of -27.207 eV.
In this case the position observables of two degenerated eigenstates are not well defined.
Density of states
One can plot the DOS from the three runs in the same figure using a plotting software of your choice.
As an example, we use PYTHON with the matplotlib and numpy package.
import numpy as np
import matplotlib.pyplot as plt
# load the DOS
dos_0 = np.loadtxt("[path_to_your_DOS]/KS_DOS_total_noNEO.dat")
dos_1 = np.loadtxt("[path_to_your_DOS]/KS_DOS_total_noEPC.dat")
dos_2 = np.loadtxt("[path_to_your_DOS]/KS_DOS_total_EPC.dat")
# Plot the DOS data
plt.xlim(-15,10)
plt.xlabel("KS Energy(eV)")
plt.ylabel("DOS(arb.)")
dos = dos_1
plt.plot(dos[:,0],dos[:,1],color="red",linestyle="--",label="w/o epc",alpha=0.8,lw=2)
dos = dos_2
plt.plot(dos[:,0],dos[:,1],color="blue",label="epc-17-2",alpha=0.8,lw=2)
dos = dos_0
plt.plot(dos[:,0],dos[:,1],color="black",linestyle=":",label="DFT",lw=2)
plt.legend(frameon=False)
Treating hydrogen at quantum level shows some impact on the DOS, such as:
- Peak shift at \(-10\) eV, \(4\) eV, \(8\) eV,
- Peak split at \(-5\) eV.
The EPC functional does not show a significant impact on the DOS in this case.
You are encouraged to further explore additional observables using NEO-DFT. This concludes the second part of our tutorial!
Solutions
You can find all the solutions to the exercises above by clicking on the button below: