Tutorial 2: a FLOSIC calculation for H2#
In this tutorial we will learn to do a simple FLOSIC calculation with an H2 molecule. As a reminder, the CLUSTER file is the main input file of FLOSIC. It contains the minimal information to set up the calculation. Below you will find the CLUSTER file for the hydrogen molecule.
Setting Up Input Structure#
GGA-PBE*GGA-PBE # Exchange-correlation parametrization
NONE # Point group of system
2 # No. of atoms
0.0 0.0 0.5 1 ALL # x,y,z coordinates, Atomic number, ALL means all electron
0.0 0.0 -0.5 1 ALL # x,y,z coordinates, Atomic number, ALL means all electron
0.0 0.000 # Charge and Moment
As you saw from Tutorial 1, the CLUSTER file will
have the same Exchange Correlation functional GGA-PBE*GGA-PBE, and point group symmetry NONE.
Since we are not making use of point group symmetry, the number of inequivalent atoms is the same as total number of atoms,
which is 2 in this case. You can put the atoms anywhere; in this example, we have placed one Hydrogen atom at (0,0,Z) and
the second one at (0,0,-Z). Atomic positions should be given in atomic units, Bohr.
Following the xyz coordinates is the atomic number. The example listed is for hydrogen
whose atomic number is 1. The string ALL signifies that the calculations are to be performed
at the all-electron level. It is also possible to use pseudopotentials. Only the BHS pseudopotentials
are hardwired into the code. It is also possible to use user-supplied (numerical) pseudopotentials also, but
requires more work and is not recommended for beginners.
The last line has two fields: The first field for a 0.0 net charge (i.e. a neutral molecule) and the second is 0.0 for zero spin momentum
FLOSIC Calculation: The FRMORB File#
A file called FRMORB is also required, which contains the FOD positions. An example FRMORB file for H2 is shown below:
1 1
0.00 0.00 0.00
0.00 0.00 0.00
The first line contains two fields:
The first line is the number of spin up FODs (N=1)
The second is the number of spin down FODs (M=1)
For H2, we only have one electron in each spin channel. The up FOD coordinates (x,y,z) are the following N lines, and the down FODs are the subsequent M lines.
To avoid optimizing atomic geometries, change the CALCTYPE field in NRLMOL_INPUT.DAT from LBFGS to SCF-ONLY.
This is recommended for FLOSIC calculations when the FODs are being optimized.
Running the FLOSIC Calculation#
Now, copy the example CLUSTER and FRMORB into an empty directory, and run the calculation for the H2 molecule using the following command at the prompt.
${PATH_TO_FLOSIC}/nrlmol_exe &>> log &
Understanding the SUMMARY file#
Browse through SUMMARY and look at the energies printed at each iteration of the SCF cycle. You should see that the minimum total energy is reached at self-consistency. The threshold for self-consistency is set in NRLMOL_INPUT.DAT, in the variable SCFTOLV. The default value is 1.0D-6.
Once a self-consistent calculation finishes, the FOD forces are used in a gradient optimization scheme to update the FOD positions.
You can see this for yourself, in SUMMARY, that the TRACE (i.e. the 1st column) is set to 0.0000, but the
EDFT+SIC energy changes per iteration.
25 -70.449947088 -114.945747472 115.012262729 21.999999620 -115.887906354 -115.887910090
26 -70.449832213 -114.945752918 115.012120854 21.999999620 -115.887901358 -115.887910090
27 -70.449833257 -114.945755178 115.012088824 21.999999620 -115.887899304 -115.887910090
28 -70.449765288 -114.945756086 115.012053732 21.999999620 -115.887898772 -115.887910090
29 -70.449689889 -114.945756716 115.012010580 21.999999620 -115.887898574 -115.887910090
30 0.000000000 -114.945756716 115.012010580 21.999999620 -115.887898574 -115.887910090
31 0.000000000 -114.945756716 115.012010580 21.999999620 -115.993753670 -115.993753670
32 0.000000000 -114.945756716 115.012010580 21.999999620 -116.376476586 -116.376476586
33 0.000000000 -114.945756716 115.012010580 21.999999620 -116.377998509 -116.377998509
34 0.000000000 -114.945756716 115.012010580 21.999999620 -116.382675125 -116.382675125
35 0.000000000 -114.945756716 115.012010580 21.999999620 -116.387459689 -116.387459689
Other Important Output files to keep in mind#
Look at the EVALUES file in which Kohn-Sham eigenvalues and occupation numbers are printed.
The FOD forces, along their respective FOD positions, are displayed in the records file.
The fande.out file contains the (1) iterations, (2) total DFT+SIC energies, (3) square root of the sum of the squares of the FOD forces, and (4) the max FOD forces. When optimizing FODs, this is a good file to check for the convergence of FOD forces. There is also fande.dat which holds the last iteration of fande.out.
How to Further Optimize FODs#
Re-executing the code will run another SCF calculation to be performed using the updated FOD positions. A new total energy and new FOD forces will be calculated, and the FOD positions will - again - be updated. Repeating this process will result in the optimization of the FOD positions. Convergence can be gauged by the size by the largest FOD force. When this drops below a chosen convergence criterion, the FODs are optimized.
Note
Currently, the FOD convergence criterion is hard-coded and requires recompilation. Look into electronic_geometry_lbfgs.f
A simple iterative loop can repeatedly run the code until the calculation is optimized to your criterion.
For the example of H2, there is only one FOD of each spin. Placing the FODs at any position in space for such a case will give the same energy and the force on the FOD will therefore be zero.