Tip
All input and output files can be downloaded here
.
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 B3LYPD3BJ/def2tzvp energy for a water molecule:
# This first calculation.
basis
def2tzvp
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 water1.inp n 4 > water1.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:
water1.out water1.mwfn
water1.out
is the output file for this calculation. You can find energies, molecular orbital (MO) information, spin, population, electric multipole moments in it:
SCF Energies
============
Kinetic energy: 76.16685643 Hartree
Electronnuclear attraction energy: 199.22128800 Hartree
Pseudopotential energy: 0.00000000 Hartree
Exchangecorrelation 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
HOMOLUMO (56) gap: 8.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

water1.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 water1.mwfn
$ 0
Explanation of The First Example
Now we will explain the input water1.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 inqbics/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 selfconsistent field (SCF) calculations are done. By its name, you may understand that:
charge
The charge of the molecule, like0
,+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 usespin2p1 3
, then 2 unpaired electrons will occupy 2 alpha orbitals.
type
The value isR
if a restricted Hartree–Fock (HF) or Kohn–Sham (KS) is needed, orU
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 andU
for other cases
grimmedisp ... end
This block controls how Grimme dispersion correction is applied. Here we usebj
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, saywater.xyz
orwater.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:
# This first calculation.
basis
def2tzvp
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:
Molecular Orbitals
==================
k = Gamma
Alpha HOMOLUMO (67) gap: 3.247 eV
Beta HOMOLUMO (45) 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 pc2 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
pc2
for H and O, and download it using GAUSSIAN format.
Attention
Replace all D
to E
in the basis set definitions!
Now we have two ways:
The first way: Simply copy it to
basis ... end
block:
# This first calculation.
basis
# From: https://www.basissetexchange.org/basis/pc2/format/gaussian94/?version=0& elements=1,8
H 0
S 4 1.00
0.754230E+02 0.240650E02
0.113500E+02 0.184870E01
0.259930E+01 0.897420E01
0.735130E+00 0.281110E+00
S 1 1.00
0.231670E+00 0.100000E+01
S 1 1.00
0.741470E01 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.535190E03
0.221730E+04 0.413750E02
0.504740E+03 0.212450E01
0.142870E+03 0.824530E01
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.192750E05
0.504740E+03 0.579640E04
0.142870E+03 0.794940E03
0.463000E+02 0.731250E02
0.163370E+02 0.405740E01
0.598280E+01 0.915940E01
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.689490E02
0.139350E+02 0.490050E01
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,
pc2OH.txt
. Assume you put it in/home/you/calc/pc2OH.txt
, then write this file name inbasis ... end
block:
# This first calculation.
basis
/home/you/calc/pc2OH.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 pc2 (defined above) for oxygen but ccpvdz for hydrogen, you can use the following basis ... end
:
basis
element # This indicates that Qbics will assign basis set element by element.
O /home/you/calc/pc2OH.txt
H ccpvdz
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.
Attention
You must write BOTH valence basis set and pseudopotential in the input file. If pseudopotential is not written explicitly in the input file, Qbics will NOT use it.
This is because for some elements, there exist both allelectron basis set and valenceonly basis set with pseudopotential at the same time, so you should make it clear in the input file.
The valence basis sets and core pseudopotentials are stored in qbics/basis
and qbics/pseudopotential
, respectively. For an element, you must use them consistently:
Valence Basis Set 
Pseudopotential 

def2X 
def2ecp 
(aug)ccXpp 
ccecp 
lanlX 
lanlecp 
For example, you can use def2TZVPP and defecp for any element, but NEVER use def2TZVPP and lanlecp together!
Assume you want to calculate B3LYPD3BJ energy for \(AuCl_3\), you can use def2 series:
basis
def2tzvp
end
pseudopotential
def2ecp
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 ccpvdz for chlorine and ccpvdzpp for gold, then the input is:
basis
element
Cl ccpvdz
Au ccpvdzpp
end
pseudopotential
ccecp
end
This means that for the valence and core part of gold, ccpvdz 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_3\), we want to use SBKJC valence basis set and pseudopotentials, then
Go to basis set exchange website: https://www.basissetexchange.org/, choose
SBKJCVDZ
for Au and Cl, and download it using GAUSSIAN format.
Paste basis set and pseudotential definitions in
basis ... end
andpseudopotential ... end
block, respectively:
basis
# From: https://www.basissetexchange.org/basis/sbkjcvdz/format/gaussian94/?version=0&elements=17,79
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/sbkjcvdz/format/gaussian94/?version=0&elements=17,79
CL 0
CLECP 2 10
d potential
1
1 4.8748300 3.4073800
sd potential
2
0 17.0036700 6.5096600
2 4.1038000 42.2778500
pd potential
2
0 8.9002900 3.4286000
2 3.5264800 22.1525600
****
AU 0
AUECP 4 60
g potential
1
1 4.3876300 10.7235800
sg potential
3
0 1.5563600 6.3561200
2 3.7159300 364.4403900
2 4.0679200 428.1975300
pg potential
3
0 1.1879800 4.4151800
2 3.0155100 136.5755000
2 3.5958800 194.2053500
dg potential
2
0 35.2500000 8.8819800
2 5.0230700 86.7661200
fg 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/SBJKCECP.txt
, respectively, and write these file names in the input file:
basis
/home/you/calc/SBJKC.txt
end
pseudopotential
/home/you/calc/SBJKCECP.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_{3}NH_{3}^{+}OH^{}\), we do a B3LYPD3BJ/631g(d) calculation:
basis
631g(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 ch3nh3oh1.mwfn
. We can use this as an initial guess for a more expensive calculation, like B3LYPD3BJ/6311g(2df,2pd) calculation:
basis
6311g(2df,2pd)
end
scf
charge 0
spin2p1 1
end
scfguess
type mwfn
file ch3nh3oh1.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:
Option 
Meaning 


Diagonalization of core Hamiltonian. NOT recommended. 

Superposition of atomic densities. Default. 

Read a wave function from a mwfn file given by 

Superposition of molecular fragment densities. See below. 
Attention
For an input file x.inp
, the program will rewrite x.mwfn
during the calculation, so if you want to keep the x.mwfn
, you can copy it to something like x1.mwfn
, and write
scfguess
type mwfn
file x1.mwfn
end
Then the program will read x1.mwfn
for initial guess but not overwrite x.mwfn
.
Guess From Fragment Superposition
Consider \(CH_{3}NH_{3}^{+}OH^{}\), it seems to be logical that one can use the superposition of \(CH_{3}NH_{3}^{+}\) and \(OH^{}\) as initial guess. This is very easy:
basis
6311g(2df,2pd)
end
scf
charge 0
spin2p1 1
end
scfguess
type fragden
frag +1 1 18
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 18
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 812 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.E6
density_cov 1.E6
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.E6
and 1.E6
, respectively. max_it
means the maximum number of SCF iterations. You can change these when you need.