Tutorial 2: Understanding electronic structure of molecules using NRLMOL

In this tutorial we will learn to optimize the geometries of diatomic molecules. We choose a case of H2 molecule. Again, as a reminder the CLUSTER file is the main input file of NRLMOL. It contains the minimal information to set up the calculation. See below the CLUSTER file for the hydrogen molecule.

GGA-PBE*GGA-PBE            # Exchange-correlation parametrization
GRP                        # 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

 Hydrogen molecule. Input auto generated by Perl script.

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 GRP. It refers to point group symmetry. The GRP means the symmetry point group operations will be read from the GRPMAT file. A few selected symmetry groups such as Ih, Td, Th, C3V etc. can also be provided here. In this case the GRPMAT corresponding to specified point group will be automatically generated by NRLMOL. Even though the H2 molecule posses symmetry we will perform calculation without symmetry. We will learn how to use symmetry to speed up calculations in the next tutorial. So keep “GRP” line as it is. The GRPMAT file will be created automatically and will contain the identity matrix only.

The third line specifies the number of inequivalent atoms in calculations. Since our calculation does not use symmetry, the number of inequivalent atoms. Since we are not making use of point group symmetry, the number of inequivalent atoms is same as total number of atoms, which is 2 in this case.

The lines following the third line contain information for each of the inequivalent atoms. The first three numbers are the cartesian positions in atomic units, its atomic number, and string ALL. In this example, we have placed one hydrogen atom at (0,0,Z) and the second one at (0,0,-Z). You can put the atoms anywhere. Atomic positions should be given in atomic units, Bohr. Followed by xyz coordinates is the atomic charge. 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 in example has two fields. The first field is 0.0 which means perform calculation for the neutral molecule. If it is 1 then the calculations will be performed for singly charged cation of N atoms. The next field which is 0 in this example corresponds to number of unpaired electrons in the system. There are no unpaired electrons. Lines after Charge and Moment line are ignored.

Now, copy the example input in the file called CLUSTER and run the calculation for H2 molecule using the following sequences of commands at the prompt.

$ rm -f GEOCNVRG SYMBOL ISYMGEN INPUT RUNS

$ ./cluster < print.SCF

Now, browse through the file print.SCF and look at the energies printed at each iteration. Also, look at the EVALUES file in which Kohn-Sham eigenvalues and occupation numbers are printed. Try to understand the SCF process. Now, let’s open the GEOCNVRG file:

  1.0000000000000000E-003
     CONVERGE FALSE
     ENERGY=   -1.113884538059298
     TOTAL GRADIENT=   0.4658517560326668
     LARGEST NUCLEAR GRADIENT=   0.3294070519678478

The number on the first line is the optimization threshold. If the RMS gradient falls below this value then optimization of structure is complete. The following line CONVERGENCE FALSE tells us that the structure is not optimized yet, that is, the atoms have still significant forces. The next line is the total energy in Hartree for the current geometry. The last line print out the maximum rms force. To get more information on the forces, look at the file FRCOUT.

                  -1.1138845381     1 Energy and Symbol NCALC
           0.00000000           0.00000000           0.50000000   1   1
  -0.593228680992D-07  -0.593228686397D-07   0.329406819469D+00
           0.00000000           0.00000000          -0.50000000   1   1
  -0.593227768567D-07  -0.593227770511D-07  -0.329407051968D+00  <--- MAX-FORCE
   0.557780120356D-07   0.557780123507D-07  -0.768552000929D-07
   0.000000000000D+00   0.000000000000D+00   0.000000000000D+00

The first line is the total energy in Hartree. It will also be printed in GEOCNVRG. Followed by these are atomic positions (x,y,z) components of atom 1, which is followed by the x,y,z component’s of forces. Similar, information is printed on subsequent lines. The second last line contains DIPOLE MOMENT while the last line prints out the external electric field. The H2 molecule does not have dipole moment and our calculations are in absence of electric field. Note that since we placed atoms on Z axis only the Z component of force is nonzero. The forces are significant and tell us that atoms should be moved away from each other. Since the forces are significant we need to repeat calculation until the forces vanish (fall below the threshold). So repeat following command until forces are zero which is easy to figure out by looking at GEOCNVRG.

$ ./cluster > print.SCF

You can also use following script to automate optimization. For the bash shell it is :

#!/bin/bash
  export niters=41

  export iter=1
 while [ $iter -le $niters ]; do

    #${mpirun}  -np ${NP}   ./bcluster > print.${iter}
    #${mpirun} -machinefile $PBS_NODEFILE -np ${NP} ./cluster  > print.${iter}
    ./cluster  > print.${iter}

    cat FRCOUT >> FRCOUT-TOT
    cp FRCOUT FRCOUT.${iter}
      if [ -e "EXIT" ]; then
            rm -f EXIT; exit
      fi
      if [ -e "GEOCNVRG" ]; then
        export x=`cat GEOCNVRG| grep CONV |awk '{print $2}'`
          if [ $x == 'TRUE' ]; then
             echo "Looks like geometry is converged at step ${iter}. \n Stopping..."
             exit
          fi
      fi
   let  iter=${iter}+1
 done # end_of_while_loop

Cut and paste the above script in a file called opt.bash. Then do the following to run it.

$ chmod +x opt.bash

$ ./opt.bash &

The first command will make the file opt.bash executable and the second one will run the program in background. The final GEOCNVRG file should look something like

 1.0000000000000000E-003
 CONVERGE TRUE
 ENERGY=   -1.166529242444089
 TOTAL GRADIENT=   4.6930900067498795E-004
 LARGEST NUCLEAR GRADIENT=   3.3185157688592808E-004

Note the second line says that CONVERGENCE is TRUE . This indicates that convergence was obtained for the optimization job, i.e. the structure is optimized. The final forces are in FRCOUT.

                  -1.1665292424     6 Energy and Symbol NCALC
          -0.00000200          -0.00000200           0.70838900   1   1
  -0.156998838711D-06  -0.156998838805D-06   0.331851502536D-03
          -0.00000200          -0.00000200          -0.70838900   1   1
  -0.156998838240D-06  -0.156998838613D-06  -0.331851502610D-03  <--- MAX-FORCE
   0.311218459833D-06   0.311218459219D-06  -0.930366894636D-13
   0.000000000000D+00   0.000000000000D+00   0.000000000000D+00

The final optimized coordinates are stored in “XMOL.DAT” (They are in Angstrom).

   2
Cluster output
    1    0.00000    0.00000    0.37486
    1    0.00000    0.00000   -0.37486

So the optimized bond length is 0.749 Angstrom.