Tutorial 3: FOD optimization of a water molecule#
This tutorial explains how to do a FODs optimization taking the example of the water molecule.
Setting Up Input Structure and Input Parameters#
It is recommended to run a DFT calculation first to get the DFT optimized density. To do this, we start from the usual input:
GGA-PBE*GGA-PBE (DF TYPE EXCHANGE*CORRELATION)
NONE (TD, OH, IH, X, Y, XY, ... OR GRP)
3 (NUMBER OF ATOMS)
0.000 0.000 0.0 8 ALL (R, Z, ALL-ELECTRON)
1.443 1.128 0.0 1 ALL (R, Z, ALL-ELECTRON)
-1.443 1.128 0.0 1 ALL (R, Z, ALL-ELECTRON)
0.000 0.000 (NET CHARGE AND NET SPIN)
In the nrlmol_input.dat file we change CALCTYPEV from LBFGS to SCF-ONLY to switch off the molecular optimization.
We also turn on the unrestricted calculation setting SPNPOLV = 'Y'. Then we run NRLMOL in a new directory with only CLUSTER and nrlmol_input.dat files:
${PATH_TO_FLOSIC}/nrlmol_exe > print.DFT
If the calculation finished correctly, all outputs should be generated, and we should have the following in the RUNS file
0 1 ITBEG, NCALC
4 4 START: 0=SCR.NUC, 1=HAM, 2=POT, 3=LSF, 4=WFUNC, 5=WFUNC_FRAG
0 START HAMILTONIAN IS INTERPOLATED: 0=NO, 1=YES
Important
Check that you have 4 4 in the second line. This tells FLOSIC to use the wave function stored in WFOUT to start the next calculation.
The second number in the first line may be different if a molecular optimization was carried out.
Creating FODs for a FLOSIC Calculation#
Now, we have to include the FOD positions. Here, we can use the Monte Carlo FOD generator (fodMC). For that, we copy XMOL.xyz into a file named system and edit the second and last lines as shown:
3
angstrom fix1s
O 0.00000 0.00000 0.00000
H 0.76360 0.59691 0.00000
H -0.76360 0.59691 0.00000
con_mat
It might have an empty line at the end. The fodMC code and documentation can be found here. When fodMC is executed 3 files are generated: CLUSTER, FRMORB and Nuc_FOD.xyz. Thus, it will overwrite our CLUSTER file if executed in the same directory. In our case, it does not matter because whenever the SYMBOL file is present CLUSTER is ignored. FRMORB contains the FODs in Bohrs for the unrestricted calculation (first the up, then the down). You should get something like this:
5 5
0.0000000000000000 0.0000000000000000 0.0000000000000000
1.2265460953713494 0.95879731507086485 0.0000000000000000
-1.2265460953713494 0.95879731507086485 0.0000000000000000
9.8690853618931581E-005 -0.84709659045951247 -1.0165159441533127
-1.5276582823942571E-005 -0.84709659045951247 1.0165159488293354
0.0000000000000000 0.0000000000000000 0.0000000000000000
9.8690853618931581E-005 -0.84709659045951247 -1.0165159441533127
-1.5276582823942571E-005 -0.84709659045951247 1.0165159488293354
1.2265461383759642 0.95879734868779065 0.0000000000000000
-1.2265461383759642 0.95879734868779065 0.0000000000000000
If you do not have fodMC, you can copy-paste the previous FRMORB. It is important to always check how the FODs are distributed. For that, you can visualize Nuc_FOD.xyz with the software of your choice. The spin up FODs have the label X and the down FODs have the symbol He by default, so be careful if the He atom is included in your system.
Running the FLOSIC code again in this directory will now cause a FLO-SIC-PBE calculation to be run. (The existence of FRMORB is the flag for running a FLO-SIC calculation).
$ PATH_TO_FLOSIC/nrlmol_exe > print.001
If everything goes well, we get now all the FLOSIC output files. Check the SUMMARY file. It should look like this (note that some columns of this file do not fit on the page):
IT TRACE ETOT EKIN+ENONLOC CHARGE EDFT+SIC LOWEST
1 -42.416153386 -76.326968147 76.469793374 9.999999301 -76.326968147 0.000000000
2 -40.529422883 -76.287479323 74.155137275 9.999999484 -76.287479323 -76.326968147
3 -41.037801666 -76.360861763 75.120435771 9.999999577 -76.360861763 -76.287479323
4 -41.547677518 -76.386605548 76.071092878 9.999999652 -76.386605548 -76.360861763
5 -41.579433885 -76.387012272 76.156784862 9.999999660 -76.387012272 -76.386605548
6 -41.537373241 -76.387091483 76.105263133 9.999999659 -76.387091483 -76.387012272
7 -41.535132409 -76.387091342 76.101118384 9.999999659 -76.387091342 -76.387091483
8 -41.541270575 -76.387093972 76.111112778 9.999999661 -76.387093972 -76.387091342
IT TRACE ETOT EKIN+ENONLOC CHARGE EDFT+SIC LOWEST
1 0.000000000 -76.387093972 76.111112778 9.999999661 -76.325884480 0.000000000
2 -47.383872787 -76.377766430 75.943917202 9.999999707 -76.336806819 -76.325884480
3 -47.434397301 -76.376079958 76.027004027 9.999999708 -76.337098863 -76.336806819
4 -47.441559320 -76.375691097 76.055957777 9.999999706 -76.337109578 -76.337098863
5 -47.324136604 -76.375845538 75.998150088 9.999999702 -76.337049496 -76.337109578
6 -47.312036177 -76.376132863 76.000962346 9.999999700 -76.337038707 -76.337049496
7 -47.315234163 -76.376185290 76.002987681 9.999999700 -76.337041096 -76.337038707
8 -47.332770384 -76.376285169 76.024925697 9.999999700 -76.337087551 -76.337041096
9 -47.303951079 -76.376194186 75.999998287 9.999999700 -76.337026990 -76.337087551
10 -47.374868090 -76.376350868 76.064567671 9.999999702 -76.337132117 -76.337026990
11 -47.381393225 -76.376294921 76.083303549 9.999999703 -76.337128855 -76.337132117
12 -47.365402178 -76.376255647 76.070377576 9.999999703 -76.337133358 -76.337128855
13 -47.358580387 -76.376246217 76.060102340 9.999999702 -76.337132539 -76.337133358
14 -47.369783782 -76.376261648 76.069168121 9.999999702 -76.337133952 -76.337132539
The first block is for the normal DFT calculation with the 3rd and 6th column exactly the same (no SIC). The second block is for the FLOSIC run. In the standard output file – renamed print.001 in this example (see the execution line above) – search for the word ITERATION.
ITERATION 1
=============
READING OLD WAVEFUNCTIONS FROM FILE WFOUT
MREP = 1
N_OCC, NBASF = 5 73 73
N_OCC, NBASF = 5 73 73
There it states that the wave function is read from WFOUT and gives the occupied and total number of orbitals. After giving the occupancies, it calculates the Lowdin overlap eigenvalues.
LOWDEN OVERLAP EIGENVALUES:
0.751388 0.791811 1.06742 1.07216 1.31722
BACK FROM LOWSIC
CALLING DIAGGE: 1000 5
LOWDEN OVERLAP EIGENVALUES:
0.751388 0.791811 1.06742 1.07216 1.31722
These values correspond to the eigenvalues of the overlap matrix formed from the Fermi orbitals. This is diagonalized in the Lowdin process. When one or more of these eigenvalues are smaller than 1E-08, it means that two or more of the Fermi orbitals are identical and the calculation stops with a message of bad FOD positions. At the end of the SCF cycle we get the following summary of energy contributions. Some of them are in the SUMMARY file too.
SUMMARY OF ENERGY CONTRIBUTIONS:
================================
TOTAL ENERGY: -76.337134
NUCLEAR REPULSION: 9.082196
LOCAL POTENTIAL: -198.922706
MEAN-FIELD COULOMB: 46.647189
NONLOCAL POTENTIAL: 0.000000
KINETIC: 76.069168
LOCAL EXCHANGE: -8.108217
LOCAL CORRELATION: -0.660695
NONLOCAL EXCHANGE: -0.816267
NONLOCAL CORRELATION: 0.333071
EXTERNAL ELECTRIC FIELD: 0.000000
Once self-consistency is reached, the FOD forces are calculated and FOD positions are updated in FRMIDT (and/or in FRMORB) using the chosen optimization method. The default is scaled LBFGS, but it sometimes gets stuck or fails to obtain the next step. In those cases, we may switch off this variable in NRLMOL_INPUT.DAT SCALEDLBFGSV = ‘N’ to use the conjuged gradient for the FOD optimization. As the FODs are already updated and we also have the wave function and the RUNS file, we just run FLOSIC in the same directory to do the SCF cycle with the new FODs.
$ PATH_TO_FLOSIC/nrlmol_exe > print.002
We can sucessively do this or use a script similar to that in the previous tutorial. After doing 5 steps of FOD optimization we can see in fande.out how the energy is decreasing and converging and the forces diminishing (not at every step).
1 -76.337133952182 0.533401361320E-02 0.376016287388E-02
2 -76.336140893884 0.316012354712E-01 0.222645230322E-01
3 -76.337242375732 0.165238505484E-02 0.109072276621E-02
4 -76.337256822409 0.184877969533E-02 0.129662662282E-02
5 -76.337254671865 0.203173772471E-03 0.997893243061E-04
All the used FODs, obtained forces, and corrected total energies can be seen in the records file.