C\ :sub:`60` Molecule
=====================
.. |c60| replace:: C\ :sub:`60`
This example demonstrates a number of Qbox features in calculations
of the electronic structure of a |c60| molecule
The example includes:
* Calculation of the ground state energy
* Calculation of Kohn-Sham eigenvalues and energy gap
* Maximally Localized Wannier Functions
Ground state energy calculation
-------------------------------
The Kohn-Sham energy of the |c60| molecule is computed using the Local
Density Approximation (LDA).
Pseudopotential
^^^^^^^^^^^^^^^
The pseudopotential used for this example is defined in the file
:download:`C_HSCV_LDA-1.0.xml`.
This file must be present in the current directory when running the first
electronic structure calculation of the |c60| molecule.
Sample definition file
^^^^^^^^^^^^^^^^^^^^^^
The sample is defined in the file
:download:`c60.i`
where the unit cell is first defined using the :ref:`set_cmd` command and the
:ref:`cell_var` variable, the species 'carbon' is defined using the
:ref:`species_cmd` command, and
each atom is defined using the :ref:`atom_cmd` Qbox command.
.. code-block:: none
:caption: c60.i
:name: c60.i
set cell 20 0 0 0 20 0 0 0 20
species carbon C_HSCV_LDA-1.0.xml
atom C1 carbon 0.0000 1.3127 6.5487
atom C2 carbon 2.2125 2.6800 5.7037
atom C3 carbon 1.3674 4.8925 4.3362
atom C4 carbon -1.3674 4.8925 4.3362
atom C5 carbon -2.2125 2.6801 5.7036
atom C6 carbon 0.0000 -1.3124 6.5488
atom C7 carbon -2.2125 -2.6798 5.7037
atom C8 carbon 2.2125 -2.6798 5.7037
atom C9 carbon -1.3674 -4.8925 4.3364
atom C10 carbon 1.3674 -4.8925 4.3364
atom C11 carbon 4.3363 1.3674 4.8925
atom C12 carbon 4.3363 -1.3674 4.8925
atom C13 carbon 5.7037 2.2126 2.6800
atom C14 carbon 5.7037 -2.2124 2.6800
atom C15 carbon 6.5488 0.0000 1.3126
atom C16 carbon 2.6800 5.7037 2.2125
atom C17 carbon 4.8925 4.3363 1.3674
atom C18 carbon 1.3126 6.5488 -0.0000
atom C19 carbon 4.8925 4.3362 -1.3674
atom C20 carbon 2.6800 5.7037 -2.2125
atom C21 carbon -2.6800 5.7037 2.2125
atom C22 carbon -1.3126 6.5488 -0.0000
atom C23 carbon -4.8925 4.3363 1.3673
atom C24 carbon -2.6800 5.7036 -2.2125
atom C25 carbon -4.8925 4.3362 -1.3675
atom C26 carbon -4.3363 1.3675 4.8924
atom C27 carbon -5.7037 2.2125 2.6799
atom C28 carbon -4.3363 -1.3673 4.8925
atom C29 carbon -6.5488 0.0000 1.3126
atom C30 carbon -5.7037 -2.2125 2.6800
atom C31 carbon 1.3674 -4.8926 -4.3362
atom C32 carbon 2.2125 -2.6801 -5.7036
atom C33 carbon 0.0000 -1.3127 -6.5487
atom C34 carbon -2.2125 -2.6801 -5.7036
atom C35 carbon -1.3674 -4.8926 -4.3362
atom C36 carbon 2.6800 -5.7037 -2.2125
atom C37 carbon 1.3126 -6.5488 0.0001
atom C38 carbon 4.8925 -4.3363 -1.3674
atom C39 carbon 2.6800 -5.7036 2.2126
atom C40 carbon 4.8925 -4.3362 1.3674
atom C41 carbon 4.3363 -1.3675 -4.8924
atom C42 carbon 5.7037 -2.2125 -2.6800
atom C43 carbon 4.3363 1.3673 -4.8925
atom C44 carbon 6.5488 0.0000 -1.3126
atom C45 carbon 5.7037 2.2125 -2.6800
atom C46 carbon 0.0000 1.3124 -6.5488
atom C47 carbon 2.2125 2.6798 -5.7037
atom C48 carbon -2.2125 2.6798 -5.7037
atom C49 carbon 1.3674 4.8923 -4.3364
atom C50 carbon -1.3674 4.8923 -4.3364
atom C51 carbon -4.3363 -1.3674 -4.8925
atom C52 carbon -4.3363 1.3674 -4.8925
atom C53 carbon -5.7037 -2.2125 -2.6800
atom C54 carbon -5.7037 2.2125 -2.6800
atom C55 carbon -6.5488 0.0000 -1.3126
atom C56 carbon -2.6800 -5.7037 -2.2125
atom C57 carbon -4.8925 -4.3363 -1.3674
atom C58 carbon -1.3126 -6.5488 0.0000
atom C59 carbon -4.8925 -4.3362 1.3674
atom C60 carbon -2.6800 -5.7036 2.2125
The calculation of the electronic ground state is done using the
Qbox script :download:`c60gs.i`
.. code-block:: none
:linenos:
:caption: c60gs.i
:name: c60gs.i
c60.i
set ecut 35
set wf_dyn PSDA
set scf_tol 1.e-8
randomize_wf 0.01
run -atomic_density 0 40 10
set wf_diag T
run 0
save c60gs.xml
* In line 1, the script ``c60.i`` is invoked and the contents of the
script are executed, defining the unit cell, the species ``carbon``
and the position of all atoms.
* Line 2 defines the plane wave energy cutoff :ref:`ecut_var` to be 35 Ry.
This resizes the plane wave basis set accordingly.
Note that the unit used for the plane wave energy cutoff (Rydberg) is
an exception to the convention used otherwise in Qbox to use
atomic units (Hartree and bohr).
* Line 3 defines the algorithm used to update wave functions by setting
the variable :ref:`wf_dyn_var` to ``PSDA`` (Preconditioned Steepest
Descent with Anderson acceleration).
* Line 4 defines the tolerance used to determine convergence in
self-consistent iterations (variable :ref:`scf_tol_var`).
Iterations are considered
converged when the energy changes by less than 10\ :sup:`-8` (hartree).
* Line 5 uses the :ref:`randomize_wf_cmd` command to add random numbers to the
electronic wave function coefficients before starting the calculation. This
is necessary to avoid spurious stationary points in the optimization for
systems of high symmetry (such as |c60|). The argument is the amplitude
of the random perturbation used.
* Line 6 invokes the :ref:`run_cmd` command to perform the calculation.
The option ``-atomic_density`` specifies that the initial charge density
should be approximated by a sum of atomic charge densities. This is
recommended in calculations of the electronic structure of molecules placed
in large, empty unit cells.
The first argument (0) refers to the number of ionic steps (during which the
atoms are moved). This is zero for the present calculation since we are only
optimizing the electronic wave functions with fixed atoms. The second argument
(40) is the maximum number of self-consistent iterations performed.
The actual number of iterations may be smaller if the convergence
criterion :ref:`scf_tol_var` is reached before 40 iterations.
The third argument (10) is the number of iterations used to update wave
functions between updates of the charge density.
* Line 7 sets the :ref:`wf_diag_var` variable to T (True) to request the
calculation of the Kohn-Sham eigenvalues and eigenvectors in the next
command.
* Line 8 invokes the :ref:`run_cmd` with a single parameter of ``0``. This
computes the energy (and in this case, Kohn-Sham eigenvalues and eigenvectors)
but does not update the electronic wave functions.
* Line 9 uses the :ref:`save_cmd` command to save a full description of the
system (atomic species, atomic positions and wave functions) on file
``c60gs.xml`` for later use in other calculations.
The calculation can be run by invoking Qbox as::
$ mpirun -np 2 qb c60gs.i > c60gs.r
The Kohn-Sham eigenvalues can be extracted from the output file using the
following ``xml_grep`` command::
$ xml_grep eigenset c60gs.r
This produces the following output
.. code-block:: xml
-25.40100 -24.88202 -24.88199 -24.88182 -23.89485
-23.89475 -23.89231 -23.89221 -23.89211 -22.85343
-22.85333 -22.85327 -22.00744 -22.00434 -22.00424
-22.00412 -20.70444 -20.70438 -20.69686 -20.69681
-20.69677 -20.23972 -20.22780 -20.22758 -20.22739
-18.70004 -18.69993 -18.69752 -18.69739 -18.69726
-18.49738 -18.49730 -18.49725 -17.44521 -17.44517
-17.44512 -16.39258 -16.39250 -16.38423 -16.38405
-16.38380 -16.22562 -16.22550 -16.22547 -15.86992
-14.67072 -14.66484 -14.66479 -14.66466 -13.90474
-13.90467 -13.90461 -13.89639 -13.89635 -13.60521
-13.60518 -13.60505 -13.06592 -12.98945 -12.98941
-12.98936 -12.38356 -12.38353 -12.38341 -12.28592
-12.28587 -12.28570 -12.27684 -11.73226 -11.73223
-11.72320 -11.72310 -11.72300 -11.37361 -11.37354
-11.36955 -11.35558 -11.35541 -11.35520 -11.26133
-11.26126 -11.26116 -11.08084 -11.08072 -11.08067
-10.52467 -10.52403 -10.52394 -10.52382 -9.95413
-9.95405 -9.95401 -9.52843 -9.52834 -9.52233
-9.52224 -9.52218 -9.39338 -9.36151 -9.36141
-9.36134 -9.19187 -9.19175 -9.19172 -9.18587
-9.18577 -7.88591 -7.88583 -7.86100 -7.77139
-7.77133 -7.77128 -7.65305 -7.65290 -7.65281
-6.44126 -6.44120 -6.43436 -6.43429 -6.43423
Kohn-Sham eigenvalues are printed in units of electron-volt (eV).
The results show
that the five topmost eigenvalues belonging to the
:math:`H_g` irreducible representation of the :math:`I_h` group are close
to the value 6.44 eV. A slight symmetry breaking due to the lower (cubic)
symmetry of the unit cell separates the five eigenvalues into two groups,
centered at approximately 6.441 and 6.434 eV.
Structure optimization
----------------------
The initial positions used in the file ``c60.i`` are close but not
exactly equal to the equilibrium positions.
A small residual force on all atoms can be seen in the output of the
ground state calculation. This can be shown using e.g. the
command::
$ xml_grep atom/force c60gs.r
If we are only interested in knowing the *largest* force acting on any atom,
we can use the ``qbox_maxforce.py`` python script (in the ``util``
directory) as follows::
$ qbox_maxforce.py c60gs.r
which produces the following output::
3.628e-03 C38 3.649e-03 C50 -3.632e-03 C41
-3.621e-03 C57 3.648e-03 C49 -3.636e-03 C41
This indicates that the largest force in the *x* direction is acting on
atom ``C57``, while the largest forces in the *y* and *z* directions
are acting on atoms ``C49`` and ``C41`` respectively. The residual force
is small (of the order of ``3.0e-3`` (hartree/bohr), but it can be
further reduced by running a structure optmization.
The structure optimization starts from the previously saved file ``c60gs.xml``
that contains the results from the above
electronic ground state calculation.
The optimization is described in the Qbox script
:download:`c60cg.i` .
.. code-block:: none
:linenos:
:caption: c60cg.i
:name: c60cg.i
load c60gs.xml
set wf_dyn PSDA
set scf_tol 1.e-8
set atoms_dyn CG
run 20 20
Note that the tolerance ``scf_tol`` has been reduced to :math:`10^{-8}`
to get more accuracy in the calculations of forces.
The structure optimization can be run using e.g. the command::
$ mpirun -np 2 qb c60cg.i > c60cg.r
At the end of the optimization, the residual forces are reduced::
$ qbox_maxforce.py c60cg.r
::
-3.621e-03 C57 3.648e-03 C49 -3.636e-03 C41
3.011e-03 C44 3.005e-03 C18 3.019e-03 C1
3.070e-03 C41 3.076e-03 C23 -3.087e-03 C50
-1.877e-03 C55 -1.877e-03 C37 -1.881e-03 C33
1.390e-03 C44 1.392e-03 C22 1.393e-03 C6
2.097e-03 C11 -2.118e-03 C40 -2.111e-03 C50
-2.670e-03 C41 2.684e-03 C40 2.708e-03 C50
-4.443e-04 C57 4.544e-04 C50 -4.475e-04 C41
1.939e-04 C44 1.946e-04 C22 -1.959e-04 C33
-2.542e-04 C22 -2.494e-04 C33 -2.366e-04 C44
-9.379e-05 C54 -9.332e-05 C36 -9.183e-05 C47
9.303e-05 C44 -9.345e-05 C58 -9.396e-05 C33
6.163e-05 C44 -6.229e-05 C58 -6.284e-05 C33
-4.909e-05 C57 5.004e-05 C49 4.981e-05 C12
2.976e-05 C44 3.006e-05 C18 -3.039e-05 C46
2.287e-05 C37 2.348e-05 C46 2.298e-05 C15
1.143e-05 C22 1.173e-05 C33 -1.151e-05 C15
8.290e-06 C18 -8.480e-06 C33 -8.220e-06 C44
2.831e-05 C44 -2.895e-05 C58 -2.912e-05 C46
1.131e-05 C18 -1.148e-05 C33 -1.122e-05 C44
The residual forces are now of the order of :math:`10^{-5}` (hartree/bohr).
Maximally Localized Wannier Functions (MLWF)
--------------------------------------------
Maximally Localized Wannier Functions (MLWF) can me computed using the
:download:`c60mlwf.i` script
.. code-block:: none
:linenos:
:caption: c60mlwf.i
:name: c60mlwf.i
load c60gs.xml
compute_mlwf
The position of each MLWF center and its spread is printed, together with
the total dipole moment. The total dipole moment is the sum of the
ionic dipole moment (defined by the sum of valence charges located at the
atomic positions) and the electronic dipole moment (approximated by point
charges located at the MLWF centers). Dipole moments are printed in
units of electron*bohr.
The total dipole moment is zero (to a good approximation), as expected for
the |c60| molecule::
-0.000000 -0.001600 0.000000
0.000036 -0.000505 0.000155
0.000530