Tutorial 1: Molecular Geometry Optimization at the DFT Level

This is a very simple tutorial to get acquainted with electronic structure calculations using FLOSIC. This tutorial will explain how to run FLOSIC for molecular geometry optimizations using DFT.

The CLUSTER file is the main input file of FLOSIC. It contains the minimal information to set up a calculation. For this tutorial, we will use a CH4 molecule, which uses a CLUSTER file like the one shown below:

GGA-PBE*GGA-PBE          (DF TYPE EXCHANGE*CORRELATION)
NONE                     (TD, OH, IH, X, Y, XY, ... OR GRP)
5
                         (NUMBER OF INEQUIV. ATOMS IN CH4)
 0.0000  0.0000  0.0000 6 ALL
 1.3000  1.3000  1.3000 1 ALL
 1.3000 -1.3000 -1.3000 1 ALL
-1.3000  1.3000 -1.3000 1 ALL
-1.3000 -1.3000  1.3000 1 ALL
0.0 0.0                  (NET CHARGE AND NET SPIN)
 --------------OR-------------------
@XMOL.DAT
 IF YOU WISH TO START FROM AN XYZ XMOL FILE

We will now describe the input structure of this file.

The first line is GGA-PBE*GGA-PBE. It means that the exchange-correlation interactions in the systems are modeled within the generalized gradient approximation (GGA) using the Perdew-Burke-Ernzerhof (PBE) parametrization. This is the default functional used in NRLMOL. A few other functionals are also available.

The second line is NONE. It refers to point group symmetry of the molecule. For the purposes of the tutorial, we will not enforce symmetry. If you would like to use symmetry, a symmetry (TD,OH, etc.) can be selected in place of NONE. In these cases, the code will create a GRPMAT file containing the appropriate symmetry operations (each represented by a 3x3 matrix). If you would like to use symmetry operations directly from an existing GRPMAT file, replace NONE with GRP.

The third line contains 5. It specifies the number of inequivalent atoms in the calculation. We’re running a CH4 calculation. So the number of atoms is 5 (1 C and 4 H).

The fourth line contains the Cartesian coordinates in atomic units of the Carbon atom. Following the xyz coordinates is the nuclear charge. The example listed is for Carbon whose atomic number is 6. The string ALL means include all (that is 6 in this case) electrons into the calculation. The next 4 lines are the hydrogen atoms, which follow the same format.

The ninth line in the example file has two fields, charge and moment. The first field is 0.0 which means to perform the calculation for the neutral molecule. If it is 1/-1 then the calculations will be performed for a cation/anion of CH4. The next field, which is also 0.0 in this example, corresponds to the number of unpaired electrons in the system. CH4 is a closed shell system, so it has no unpaired electrons. Lines after the Charge and Moment line are ignored.

Now, create an empty directory and execute the code inside of it. Multiple files should be created, including a CLUS- TER file. Copy the input from this example into the file called CLUSTER, replacing the default text. Change the CALCTYPE option in the NRLMOL_INPUT.DAT file to “LBFGS” and run the calculation for CH4.

$ PATH_TO_FLOSIC/nrlmol_exe

Open the GEOCNVRG file. If you have done everything correctly then the energy should be -40.448650 Hartree.

A new atomic geometry will be appended to the SYMBOL file. SYMBOL is created from the data in CLUSTER. The new geometry was created by a gradient optimization routine (here: LBFGS). The file FRCOUT.G0 contains the atomic forces for the previous atomic geometry. Running the code again will carry out a calculation at the updated molecular geometry and a new total energy and new atomic forces will be computed. A new update of the atomic coordinates will also be written into SYMBOL. Repeating this process several times will result in a local minimum energy geometry. This happens when the maximum force falls below the criterion set in GEOCNVRG.