Density Functional Theory Calculations

This tutorial will lead you step by step to do quantum chemical calculations using Qbics.

The First Example: Water

We will now do the first calculation using Qbics. The following input file will do B3LYP-D3BJ/def2-tzvp energy for a water molecule:

water-1.inp

# This first calculation.

basis
    def2-tzvp
end

scf
    charge  0
    spin2p1 1
    type    R # You do not have to write it. The program will determine it by itself.
end

grimmedisp
    type bj
end

mol
   O  0.    0.00000000   -0.11081188
   H  0.    -0.78397589   0.44324751
   H  0.    0.78397589   0.44324751
end

task
    energy b3lyp
end

Then run it:

$ qbics water-1.inp -n 4 > water-1.out

Here, -n 4 means Qbics will use 4 CPU cores for parallization. This calculation is very fast. After calculation, you will find 2 new files:

water-1.out water-1.mwfn

water-1.out is the output file for this calculation. You can find energies, molecular orbital (MO) information, spin, population, electric multipole moments in it:

water-1.out

SCF Energies
============
Kinetic energy:                              76.16685643 Hartree
Electron-nuclear attraction energy:        -199.22128800 Hartree
Pseudopotential energy:                       0.00000000 Hartree
Exchange-correlation energy:                 -7.55662337 Hartree
Electron Coulomb energy:                     46.77766964 Hartree
Electron exchange energy:                    -1.78638759 Hartree
Nuclear repulsion energy:                     9.15711600 Hartree
Grimme dispersion energy:                    -0.00057358 Hartree
--------------------------------------------------------------------------
SCF energy (E):                             -76.46323045 Hartree
Virial quotien (V/T):                        -2.00389112

Molecular Orbitals
==================
k = Gamma
HOMO-LUMO (5-6) gap: 8.990 eV8.990 eV
    #   Occupancies     Energies/Hartree
    1   2.000     -19.12516259
    2   2.000     -1.01223986
    3   2.000     -0.54372159
    4   2.000     -0.38293836
...
Spin
====
Expected   <S^2>: 0.00000; S = 0.00000
Calculated   <S^2>: 0.00000; S = 0.00000

Mulliken Populations
====================
     #  Symbol          Charge            Spin
--------------------------------------------------------
     1       O     -0.64331377      0.00000000
     2       H      0.32165689      0.00000000
     3       H      0.32165689      0.00000000
--------------------------------------------------------
   Sum             -0.00000000      0.00000000
--------------------------------------------------------

Electric Multipole Moments
==========================
     #                Total           Electronic              Nuclear Unit
----------------------------------------------------------------------------------------------
Charge:
     0              -0.0000             -10.0000              10.0000 |e|
Dipole moment:
     X              -0.0000              -0.0000               0.0000 Debye
     Y               0.0000               0.0000               0.0000 Debye
     Z               1.9936               1.9936              -0.0000 Debye
 Totla               1.9936                                           Debye
Quadrupole moment:
    XX              -7.5861              -7.5861               0.0000 Debye*Angstrom
    XY              -0.0000              -0.0000               0.0000 Debye*Angstrom
    XZ              -0.0000              -0.0000               0.0000 Debye*Angstrom
    YY              -4.1639             -10.0682               5.9043 Debye*Angstrom
    YZ               0.0000               0.0000              -0.0000 Debye*Angstrom
    ZZ              -6.4570              -8.8162               2.3592 Debye*Angstrom
----------------------------------------------------------------------------------------------

water-1.mwfn is the wave function file, which can be read by Qbics and Multiwfn. You can do all wave function visualization and analysis with Multiwfn and this file. For example, the following code will show MOs:

$ Multiwfn water-1.mwfn 
$ 0  

Explanation of The First Example

Now we will explain the input water-1.inp of the last example.

• Anything after # is treated as comments. You can write anything anywhere in the input file.
• basis ... end This block define the basis set used the calculation. You can find all available basis sets in qbics/basis. Just input the file name. You can also define your own basis set, or set basis set for different elements. See below.
• scf ... end This block controls how the self-consistent field (SCF) calculations are done. By its name, you may understand that:

• charge The charge of the molecule, like 0, +3, -1, etc.
• spin2p1 The spin multiplicity. For a molecule with n unpaired electrons, this total spin is $S=\frac{n}{2}$, so spin multiplicity is 2S + 1 = n + 1. For example, for a triplet state of water, you will use spin2p1 3, then 2 unpaired electrons will occupy 2 alpha orbitals.
• type The value is R if a restricted Hartree--Fock (HF) or Kohn--Sham (KS) is needed, or U is an unrestricted one is needed. If you do not write this, the program will determine it automatically according to spin multiplicity: R for singlet and U for other cases

• grimmedisp ... end This block controls how Grimme dispersion correction is applied. Here we use bj to indicate Becke--Johnson dampling D3 correction.
• mol ... end This block gives the molecule coordinates in XYZ format. You can also simply give a coordinate file name in XYZ or PDB format, say water.xyz or water.pdb.
• task ... end This tells Qbics what it should do. energy b3lyp means it will use B3LYP to calculate energy. You can put several tasks in this block.

Based on the explanations above, you can now try to calculate the triplet state of water:

water-2.inp

# This first calculation.

basis
    def2-tzvp
end

scf
    charge  0
    spin2p1 3
    type    U # If you do not write it. The program will determine it by itself.
end

grimmedisp
    type bj
end

mol
   O  0.    0.00000000   -0.11081188
   H  0.   -0.78397589    0.44324751
   H  0.    0.78397589    0.44324751
end

task
    energy b3lyp
end

You can find MO occupations and spin here:

water-2.out

Molecular Orbitals
==================
k = Gamma
Alpha HOMO-LUMO (6-7) gap: 3.247 eV
Beta HOMO-LUMO (4-5) gap:  5.963 eV
              Alpha                Alpha          Beta                 Beta
    #   Occupancies     Energies/Hartree   Occupancies     Energies/Hartree
    1         1.000         -19.35044069         1.000         -19.31689276
    2         1.000          -1.21472533         1.000          -1.12845421
    3         1.000          -0.69876887         1.000          -0.66807740
    4         1.000          -0.57036857         1.000          -0.52597568
    5         1.000          -0.56380445         0.000          -0.30683178
    6         1.000          -0.10311895         0.000          -0.01035402
...
Spin
====
Expected   <S^2>: 2.00000; S = 1.00000
Calculated <S^2>: 2.00188; S = 1.00063

The spin of this calculation is 1.00063, which is very close to the theoretical value 1, indicating a very small spin contamination.

More Basis Set Configurations

Qbics hase provided flexible basis set configuration. This section will show it.

Define Your Own Basis Set

Assume you want to use a basis set called pc-2 for water. However, this is not available in qbics/basis, you can then define it by your self.

• Go to basis set exchange website: https://www.basissetexchange.org/, choose pc-2 for H and O, and download it using GAUSSIAN format.

Now we have two ways:

• The first way: Simply copy it to basis ... end block:

water-3.inp

# This first calculation.

basis
    # From: https://www.basissetexchange.org/basis/pc-2/format/gaussian94/?version=0&    elements=1,8
    H     0
    S    4   1.00
          0.754230E+02           0.240650E-02
          0.113500E+02           0.184870E-01
          0.259930E+01           0.897420E-01
          0.735130E+00           0.281110E+00
    S    1   1.00
          0.231670E+00           0.100000E+01
    S    1   1.00
          0.741470E-01           0.100000E+01
    P    1   1.00
          0.160000E+01           0.100000E+01
    P    1   1.00
          0.450000E+00           0.100000E+01
    D    1   1.00
          0.125000E+01           0.100000E+01
    ****
    O     0
    S    7   1.00
          0.147820E+05           0.535190E-03
          0.221730E+04           0.413750E-02
          0.504740E+03           0.212450E-01
          0.142870E+03           0.824530E-01
          0.463000E+02           0.236710E+00
          0.163370E+02           0.440390E+00
          0.598280E+01           0.364650E+00
     S    7   1.00
          0.221730E+04          -0.192750E-05
          0.504740E+03          -0.579640E-04
          0.142870E+03          -0.794940E-03
          0.463000E+02          -0.731250E-02
          0.163370E+02          -0.405740E-01
          0.598280E+01          -0.915940E-01
          0.167180E+01           0.209400E+00
    S    1   1.00
          0.646620E+00           0.100000E+01
    S    1   1.00
          0.216690E+00           0.100000E+01
    P    4   1.00
          0.604240E+02           0.689490E-02
          0.139350E+02           0.490050E-01
          0.415310E+01           0.182550E+00
          0.141580E+01           0.376330E+00
    P    1   1.00
          0.475490E+00           0.100000E+01
    P    1   1.00
          0.145290E+00           0.100000E+01
    D    1   1.00
          0.220000E+01           0.100000E+01
    D    1   1.00
          0.650000E+00           0.100000E+01
    F    1   1.00
          0.110000E+01           0.100000E+01
    ****
end

scf
    charge  0
    spin2p1 1
    type    R # You need not to write it. The program will determine it by itself.
end

grimmedisp
    type bj
end

mol
   O  0.    0.00000000   -0.11081188
   H  0.   -0.78397589    0.44324751
   H  0.    0.78397589    0.44324751
end

task
    energy b3lyp
end

Then do the calculation as usual.

• The second way: Save the basis set in a file called, say, pc-2-OH.txt. Assume you put it in /home/you/calc/pc-2-OH.txt, then write this file name in basis ... end block:

water-4.inp

# This first calculation.

basis
    /home/you/calc/pc-2-OH.txt
end

scf
    charge  0
    spin2p1 1
    type    R # You need not to write it. The program will determine it by itself.
end

grimmedisp
    type bj
end

mol
   O  0.    0.00000000   -0.11081188
   H  0.   -0.78397589    0.44324751
   H  0.    0.78397589    0.44324751
end

task
    energy b3lyp
end

Both ways will give the same result.

Different Basis Sets for Different Elements

Assume you want to use pc-2 (defined above) for oxygen but cc-pvdz for hydrogen, you can use the following basis ... end:

water-5.inp

basis
    element # This indicates that Qbics will assign basis set element by element.
    O /home/you/calc/pc-2-OH.txt
    H cc-pvdz
end

The word element in the first line of basis ... end block means that Qbics will assign basis set element by element. Then simply write the basis set for every element.

Use Pseudopotential

General Rules

For heavy elements, one should use pseudopotential for reasonable calculations.

The valence basis sets and core pseudopotentials are stored in qbics/basis and qbics/pseudopotential, respectively. For an element, you must use them consistently:

For example, you can use def2-TZVPP and def-ecp for any element, but NEVER use def2-TZVPP and lanl-ecp together!

Assume you want to calculate B3LYP-D3BJ energy for AuCl₃, you can use def2-series:

aucl3-1.inp

basis
    def2-tzvp
end

pseudopotential
    def2-ecp
end

scf
    charge  0
    spin2p1 1
end

grimmedisp
    type bj
end

mol
   Au                 0.00000000    0.00000000   0.
   Cl                 0.00000000   -2.33000000   0.
   Cl                 2.01783919    1.16500000   0.
   Cl                -2.01783919    1.16500000   0.
end

task
    energy b3lyp
end

Or, you may want to use cc-pvdz for chlorine and cc-pvdz-pp for gold, then the input is:

aucl3-2.inp

basis
    element
    Cl  cc-pvdz
    Au  cc-pvdz-pp
end

pseudopotential
    cc-ecp
end

This means that for the valence and core part of gold, cc-pvdz basis set and pseudopotential is used, respectively.

Define Your Own Pseudopotential

You can also define your own pseudopotential like we have done for basis set. Again for AuCl₃, we want to use SBKJC valence basis set and pseudopotentials, then

• Go to basis set exchange website: https://www.basissetexchange.org/, choose SBKJC-VDZ for Au and Cl, and download it using GAUSSIAN format.

• Paste basis set and pseudotential definitions in basis ... end and pseudopotential ... end block, respectively:

aucl3-3.inp

basis
    # From: https://www.basissetexchange.org/basis/sbkjc-vdz/format/gaussian9
    Cl     0
    SP   3   1.00
          2.225000               -.330980               -.126040
          1.173000               0.115280               0.299520
          0.385100               0.847170               0.583570
    SP   1   1.00
          0.130100               0.265340               0.340970
    ****
    Au     0
    SP   1   1.00
          1.502000               1.000000               1.000000
    SP   4   1.00
          7.419000               0.222546               0.019924
          4.023000              -1.086045               -.299997
          1.698000               1.156039               0.748919
          0.627100               0.518061               0.504023
    SP   1   1.00
          0.151500               1.000000               1.000000
    SP   1   1.00
          0.049250               1.000000               1.000000
    D    3   1.00
          3.630000               -.087402
          1.912000               0.468634
          0.842300               0.654805
    D    1   1.00
          0.375600               1.000000
    D    1   1.00
          0.154400               1.000000
    ****
end

pseudopotential
    # From: https://www.basissetexchange.org/basis/sbkjc-vdz/format/gaussian94/?version=0&elem
    CL     0
    CL-ECP     2     10
    d potential
      1
    1      4.8748300             -3.4073800
    s-d potential
      2
    0     17.0036700              6.5096600
    2      4.1038000             42.2778500
    p-d potential
      2
    0      8.9002900              3.4286000
    2      3.5264800             22.1525600
    ****
    AU     0
    AU-ECP     4     60
    g potential
      1
    1      4.3876300            -10.7235800
    s-g potential
      3
    0      1.5563600              6.3561200
    2      3.7159300           -364.4403900
    2      4.0679200            428.1975300
    p-g potential
      3
    0      1.1879800              4.4151800
    2      3.0155100           -136.5755000
    2      3.5958800            194.2053500
    d-g potential
      2
    0     35.2500000              8.8819800
    2      5.0230700             86.7661200
    f-g potential
      1
    0      1.6888100              6.2160300
    ****
end

scf
    charge  0
    spin2p1 1
end

grimmedisp
    type bj
end

mol
   Au                 0.00000000    0.00000000   0.
   Cl                 0.00000000   -2.33000000   0.
   Cl                 2.01783919    1.16500000   0.
   Cl                -2.01783919    1.16500000   0.
end

task
    energy b3lyp
end

Note, you MUST add **** between pseudopential definitions for each element.

• Or, you can save both definitions in /home/you/calc/SBJKC.txt and /home/you/calc/SBJKC-ECP.txt, respectively, and write these file names in the input file:

aucl3-4.inp

basis
    /home/you/calc/SBJKC.txt
end

pseudopotential
    /home/you/calc/SBJKC-ECP.txt
end

Good SCF Initial Guess

Qbics supports several ways to set SCF initial guess. The default one is the "superposition of atomic densities", which works quiet well in most cases. Sometimes, one may or must use other ways.

Guess From Another Calculations

Consider CH₃NH₃⁺OH⁻, we do a B3LYP-D3BJ/6-31g(d) calculation:

ch3nh3oh-1.inp

basis
   6-31g(d)
end

scf
    charge  0
    spin2p1 1
end

grimmedisp
    type bj
end

mol
C                 -0.02909313   -0.38166253   -1.63963825
H                  0.82173570   -1.03039025   -1.62769450
H                 -0.03691367    0.17958863   -2.55059158
H                 -0.92347728   -0.96516294   -1.57252437
N                  0.04484969    0.58404262   -0.44233044
H                  0.93923384    1.16754303   -0.50944432
H                 -0.80597913    1.23277035   -0.45427419
H                  0.05267023    0.02279146    0.46862289
O                  0.06914388   -0.26694332    2.08310089
H                 -0.31212699   -0.12725829    2.95299779
end

task
    energy b3lyp
end

After calculation, we will obtain a file ch3nh3oh-1.mwfn. We can use this as an initial guess for a more expensive calculation, like B3LYP-D3BJ/6-311g(2df,2pd) calculation:

ch3nh3oh-2.inp

basis
   6-311g(2df,2pd)
end

scf
    charge  0
    spin2p1 1
end

scfguess
    type mwfn
    file ch3nh3oh-1.mwfn
end

grimmedisp
    type bj
end

mol
C                 -0.02909313   -0.38166253   -1.63963825
H                  0.82173570   -1.03039025   -1.62769450
H                 -0.03691367    0.17958863   -2.55059158
H                 -0.92347728   -0.96516294   -1.57252437
N                  0.04484969    0.58404262   -0.44233044
H                  0.93923384    1.16754303   -0.50944432
H                 -0.80597913    1.23277035   -0.45427419
H                  0.05267023    0.02279146    0.46862289
O                  0.06914388   -0.26694332    2.08310089
H                 -0.31212699   -0.12725829    2.95299779
end

task
    energy b3lyp
end

Here, scfguess ... end block controls the initial guess with type, which can be:

 scfguess
      type mwfn
      file x-1.mwfn 
   end

Guess From Fragment Superposition

Consider CH₃NH₃⁺OH⁻, it seems to be logical that one can use the superposition of CH₃NH₃⁺ and OH⁻ as initial guess. This is very easy:

ch3nh3oh-3.inp

basis
   6-311g(2df,2pd)
end

scf
    charge  0
    spin2p1 1
end

scfguess
    type fragden
    frag +1 1 1-8
    frag -1 1 9 10
end

grimmedisp
    type bj
end

mol
C                 -0.02909313   -0.38166253   -1.63963825
H                  0.82173570   -1.03039025   -1.62769450
H                 -0.03691367    0.17958863   -2.55059158
H                 -0.92347728   -0.96516294   -1.57252437
N                  0.04484969    0.58404262   -0.44233044
H                  0.93923384    1.16754303   -0.50944432
H                 -0.80597913    1.23277035   -0.45427419
H                  0.05267023    0.02279146    0.46862289
O                  0.06914388   -0.26694332    2.08310089
H                 -0.31212699   -0.12725829    2.95299779
end

task
    energy b3lyp
end

Here, type fragden means we will use superposition of fragments as initial guess. The fragments can be defined with frag:

frag +1 1 1-8

The first and second number is the charge and spin multiplicity of the fragment. Then is the atom indices for this fargment. It can be discontinuous like:

frag +1 1 1 2 5 8-12 15

In this case, the atom 1 2 5 8 9 10 11 12 15 will be a fragment. You can define any numbers of fragments.

SCF Convergence Condition

There are 3 options in scf ... end block to control the convergence condition:

scf
    energy_cov  1.E-6
    density_cov 1.E-6  
    max_it      100   
end

Here, energy_cov and density_cov means that the SCF is thought to converge when the energy and density change is smaller than 1.E-6 and 1.E-6, respectively. max_it means the maximum number of SCF iterations. You can change these when you need.