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.

Setting Up Input Structure#

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:

 1GGA-PBE*GGA-PBE          (DF TYPE EXCHANGE*CORRELATION)
 2NONE                     (TD, OH, IH, X, Y, XY, ... OR GRP)
 35                        (NUMBER OF INEQUIV. ATOMS IN CH4)
 4 0.0000  0.0000  0.0000 6 ALL
 5 1.1860  1.1860  1.1860 1 ALL
 6 1.1860 -1.1860 -1.1860 1 ALL
 7-1.1860  1.1860 -1.1860 1 ALL
 8-1.1860 -1.1860  1.1860 1 ALL
 90.0 0.0                  (NET CHARGE AND NET SPIN)
10 --------------OR-------------------
11@XMOL.DAT
12 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 following lines contain the cartesian positions of the atoms in atomic units (A.U.), followed by number of electrons of each atom (e.g. 6 for Carbon). The string ALL means include all (i.e. 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: net charge and net spin. The first field is 0.0 which means to perform the calculation for the neutral molecule. If the first field was 1, then the calculations would be performed for a cation of CH4. The next field which is also 0.0 in this example corresponds to the number of unpaired electrons in the system, or the net spin in multiples of \(\frac{\hbar}{2}\). CH4 is a closed shell system, so it has no unpaired electrons. Lines after the Charge and Spin fields are ignored.

Running the Calculation#

Now, create an empty directory and run the nrlmol_exe executable inside of it. Multiple default files will be created:

  • AVRGDAT

  • CLUSTER

  • GEOCNVRG

  • NRLMOL_INPUT.DAT

  • PARASAV

  • RUNS

  • SPNORB

Copy the input from this example into the CLUSTER, replacing the default input. Run the calculation for CH4 and redirect output into some log file:

${PATH_TO_FLOSIC}/nrlmol_exe &>> log &

Note

You can call your log file whatever you want. The final ampersand (&) is to run the program in the background.

Open the GEOCNVRG file, you should see the following:

   1.0000000000000000E-003
 CONVERGE FALSE
 ENERGY=   -40.466969
 TOTAL GRADIENT= 0.113937D-01
 LARGEST NUCLEAR GRADIENT= 0.569685D-02

A new atomic geometry will be appended to the SYMBOL file; this file is created from the data in CLUSTER. The new geometry is created by a gradient optimization routine (either LBFGS or CONJUGATE-GRADIENT). The file FRCOUT.G0 contains the atomic forces and coordiantes of 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 to be reached. This happens when the maximum force falls below the criterion set in GEOCNVRG.