.. tip:: All input files can be downloaded: :download:`Files `. .. tip:: For more information of this section, please refer to these pages: - :doc:`../keywords/basis` - :doc:`../keywords/pseudopotential` - :doc:`../keywords/scf` - :doc:`../keywords/scfguess` - :doc:`../keywords/grimmedisp` - :doc:`../keywords/mol` Density Functional Theory Calculations ========================================= .. contents:: :local: This tutorial will lead you step by step to do density functional theory (DFT) calculations using Qbics. Here we only do energy calculations. For geometry optimization and molecular dynamics, please refer to :doc:`opt` and :doc:`md`. Example: Single Point Energy of Water ----------------------------------------- We will now do the first calculation using Qbics. The following input file will calculate B3LYP-D3BJ/def2-tzvp energy for a water molecule: .. code-block:: :caption: water-1.inp :linenos: # 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: .. code-block:: bash $ 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: .. code-block:: 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: .. code-block:: :caption: water-1.out :linenos: 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 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 : 0.00000; S = 0.00000 Calculated : 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 in `Multiwfn `_ format. You can do all wave function visualization and analysis with Multiwfn. A faster way is to use `Qbics-MolStar `_, an online app, you can also download it on your own Windows machine. Just drag ``water-1.mwfn`` into explorer, and it will be automatically loaded. Right-click :guilabel:`water-1.mwfn` and select :guilabel:`View Molecular Orbitals`: .. figure:: figs/a1.png Now, select :guilabel:`5: -0.313564 (occ=2, RHF)`, which is the HOMO of water, then MO will be rendered: .. figure:: figs/a2.png You can also view other MOs or properties. 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 defines 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 :math:`n` unpaired electrons, this total spin is :math:`S=\frac{n}{2}`, so spin multiplicity is :math:`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: .. code-block:: :caption: water-2.inp :linenos: # Triplet state of water. 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: .. code-block:: :caption: water-2.out :linenos: 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 : 2.00000; S = 1.00000 Calculated : 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. .. image:: figs/a3.png .. 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: .. code-block:: :caption: water-3.inp :linenos: # 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: .. code-block:: :caption: water-4.inp :linenos: # 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``: .. code-block:: :caption: water-5.inp :linenos: 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. Example: Heavy Elements with 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 all-electron basis set and valence-only 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: .. list-table:: :widths: 8 8 :header-rows: 1 * - Valence Basis Set - Pseudopotential * - def2-X - def2-ecp * - (aug-)cc-X-pp - cc-ecp * - lanlX - lanl-ecp 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\ :sub:`3`, you can use def2- series: .. code-block:: :caption: aucl3-1.inp :linenos: 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: .. code-block:: :caption: aucl3-2.inp :linenos: 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\ :sub:`3`, 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. .. image:: figs/a4.png - Paste basis set and pseudotential definitions in ``basis ... end`` and ``pseudopotential ... end`` block, respectively: .. code-block:: :caption: aucl3-3.inp :linenos: basis # From: https://www.basissetexchange.org/basis/sbkjc-vdz/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/sbkjc-vdz/format/gaussian94/?version=0&elements=17,79 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: .. code-block:: :caption: aucl3-4.inp :linenos: basis /home/you/calc/SBJKC.txt end pseudopotential /home/you/calc/SBJKC-ECP.txt end Example: Set Special 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\ :sub:`3`\ NH\ :sub:`3`\ :sup:`+`\ OH\ :sup:`-`, we do a B3LYP-D3BJ/6-31g(d) calculation: .. code-block:: :caption: ch3nh3oh-1.inp :linenos: 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: .. code-block:: :caption: ch3nh3oh-2.inp :linenos: 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: .. list-table:: :widths: 2 10 :header-rows: 1 * - Option - Meaning * - ``HCore`` - Diagonalization of core Hamiltonian. **NOT** recommended. * - ``AtmDen`` - Superposition of atomic densities. **Default.** * - ``MWFN`` - Read a wave function from a mwfn file given by ``file`` as initial guess. * - ``FragDen`` - 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 ``x-1.mwfn``, and write .. code-block:: bash scfguess type mwfn file x-1.mwfn end Then the program will read ``x-1.mwfn`` for initial guess but not overwrite ``x.mwfn``. Guess From Fragment Superposition ++++++++++++++++++++++++++++++++++++ Consider CH\ :sub:`3`\ NH\ :sub:`3`\ :sup:`+`\ OH\ :sup:`-`, it seems to be logical that one can use the superposition of CH\ :sub:`3`\ NH\ :sub:`3`\ :sup:`+` and OH\ :sup:`-` as initial guess. This is very easy: .. code-block:: :caption: ch3nh3oh-3.inp :linenos: 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. To determine the atomic indices, it is very convienent to use `Qbics-MolStar `_. Just drag ``ch3nh3oh-3.inp`` into explore, right-click :guilabel:`All`, then click :guilabel:`Add Component`, then click :guilabel:`Label` in :guilabel:`Representation`, then click :guilabel:`+ Create Component`, you can see clearly the atomic indices: .. image:: figs/a5.png SCF Convergence Condition ---------------------------- There are 3 options in ``scf ... end`` block to control the convergence condition: .. code-block:: 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.