This lesson aims at showing how to get the following physical properties :
This first lesson covers the sections 1, 3, 4 and 6 of the abinis_help file.
The very first step is a detailed tour of the input and output files : you
are like a tourist, and you discover a town in a coach. You will have a bit
more freedom after that first step ...
It is supposed that you have some good knowledge of UNIX/Linux.
This lesson should take about 2 hours to be done.
1.1.a In addition to the present windows, open the second windows. Go to the Tutorial directory (that we refer as ~ABINIT/Tutorial).
1.1.b You also need a working directory. So, you should create a subdirectory of this directory, whose name might be "Work" (so ~ABINIT/Tutorial/Work). Change the working directory of windows 2 to "Work":
You will do most of the actions of this tutorial in this working directory. Copy the file t1x.files in "Work" :
1.1.c Edit the t1x.files. It is not very long (only 6 lines). It gives the informations needed for the code to build other file names ... You will discover more about this file in the section 1.1 of the abinis_help file. Please, read it now (the third window shows up when you click on this link).
1.1.d Modify the first and second lines of t1x.files file, so that they read:
t11.in
t11.out
Later, you will again modify these lines, to treat more cases. Close the t1x.files file.
Then, copy the file
~ABINIT/Tutorial/t11.in in "Work" :
Also later, we will look at this file, and learn about its content. For now,
you will try to run the code. Its location is
../../abinis
So, in the Work directory, type :
You should get inoffensive differences, like differences in the name of input files or timing differences, like the following:
5c5 < Starting date : Tue 4 Jul 2000. --- > Starting date : Thu 22 Jun 2000. 7c7 < - input file -> t11.in --- > - input file -> ../t11.in 9,10c9,10 < - root for input files -> t1xi < - root for output files -> t1xo --- > - root for input files -> t11i > - root for output files -> t11o 214c214 < - Total cpu time (s,m,h): 4.7 0.08 0.001 --- > - Total cpu time (s,m,h): 4.6 0.08 0.001 221,229c221,228 (... and what comes after that is related only to timing ...). If you do not run on a PC under Linux, you might also have small numerical differences, on the order of 1.0d-10 at most.
If you get something else, you should ask for help !
Supposing everything went well, we will now detail the different steps that took place : how to run the code, what is in the "t11.in" input file, and, later, what is in the "t11.out" and "log" output files.
1.1.e Running the code is described in the section 1.2 of the abinis_help file. Please, read it now.
1.1.f It is now time to edit the t11.in file. You can have a first glance at it. It is not very long: about 40 lines, mostly comments. Do not try to understand everything immediately. After having gone through it, you should read general explanation about its content, and the format of such input files in the section 3.1 of the abinis_help file.
1.1.g You might now examine in more details some input variables. An alphabetically ordered index of all variables is provided, and their description is found in the following files:
It is now time to have a look at the two output files of the run.
1.1.h First, edit the "log" file. You can begin to read it. It is nasty. Jump to its end. You will find ther the number of WARNINGS and COMMENTS that were issued by the code during execution. You might try to find them in the file (localize the keywords "WARNING" or "COMMENT" in this file). Some of them are for the experienced user. For the present time, we will ignore them. You can find more information about messages in the log file in the section 6.1 of the abinis_help file.
1.1.i Then, edit the "t11.out" file. You find some general information about the output file in section 6.2 of the abinis_help file. You should also :
If the code does not stop there, the input parameters are consistent. At this stage, many default values have been provided, and the preprocessing is finished.
It is worth to come back to the echo of preprocessed input data. You should
first examine the "t11.in" file in more
details, and read the meaning of each of its variables in the
corresponding input variables
file, if it has not yet been done. Then, you should examine some variables
that were NOT defined in the input file, but that appear in the echo written
in "t11.out" :
- "nband" : its value is 2.
It is the number of electronic states that will be treated by the code. It has
been computed by counting the number of valence electrons in the unit cell (summing
the valence electrons brought by each pseudopotential) then occupying the lowest
states (look at the "occ" variable),
and adding some states (at least one, maybe more, depending on the size of the
system).
- "ngfft" : its value is 30
30 30 .
It is the number of points of the three-dimensional FFT grid. It has been derived
from "ecut" and the dimension
of the cell ("acell").
The maximal number of plane waves "mpw" is mentioned in the memory evaluation
section: it is 752.
Well, this is not completely right, as the code took advantage of the time-reversal
symmetry, valid for the k-point (0 0 0), to decrease the number of planewave
by about a factor of two.
The full set of plane waves is 1503 (see later in the "t11.out" file).
The code indicates the time-reversal symmetry by a value of istwfk=2
, instead of the usual istwfk=1
default.
- "nsym" : its value is 16.
It is the number of symmetries of the system. The 3x3 matrices symrel
define the symmetries operation. In this case, none of the symmetries is accompanied
by a translation, that would appear in the variable "tnons".
The code did an automatic analysis of symmetries.
They could alternatively be set by hand, or using the symmetry builder (to be
described later).
- "xangst" and "xred"
are alternative ways to "xcart"
to specify the positions of atoms within the primitive cell.
Now, you can start reading the description of the remaining of the t11.out file, in the section 6.3 of the abinis_help file. Look at the t11.out file at the same time.
1.1.j You have read completely an output file !
Could you answer the following questions ?
(answers are given at the end of the present file)
1.2.a Starting from now, everytime a new input variable is mentioned,
you should read the corresponding descriptive section in the ABINIT help.
We will now complete the description of the meaning of each term : there are still a few indications that you should be aware of, even if you will not use them in the tutorial. These might appear in the description of some input variables ... For this, you should read the section 3.2 of the abinis_help file.
1.2.b There are three methodologies to compute the optimal distance between the two Hydrogen atoms:
The interatomic distance in the t11.in file was 1.4 Bohr. Suppose you decide
to examine the interatomic distances from 1.0 Bohr to 2.0 Bohr, by steps of
0.05 Bohr. That is, 21 calculations.
If you are a UNIX guru, it will be easy for you to write a script that will
drive these 21 calculations, changing automatically the variable "xcart"
in the input file, and then gather all the data, in a convenient form to be
plotted.
Well, are you a UNIX guru ? If not, there is an easier path, all within ABINIT
!
This is the multi-dataset mode. Detailed explanations about it can be found in sections 3.3, 3.4, 3.5 and 3.6, of the abinis_help file.
1.2.c Now, can you write an input file that will do the computation described above (interatomic distances from 1.0 Bohr to 2.0 Bohr, by steps of 0.05 Bohr) ? You might start from t11.in . Try to define a series, and to use the "getwfk" input variable (the latter will make the computation much faster).
You should likely have a look at the section that describes the "irdwfk" and "getwfk" input variables: in particular, look at the meaning of getwfk -1
Also, define explicitely the number of states (or supercell "bands") to be one, using the input variables "nband". The input file ~ABINIT/Tutorial/t12.in is an example of file that will do the job, while ~ABINIT/Tutorial/Refs/t12.out is an example of output file. If you decide to use the ~ABINIT/Tutorial/t12.in file, do not forget to change the file names in the t1x.files file ...
So, you run the code with your input file (this might take one or two minutes), examine the output file quickly (there are many repetition of sections, for the different datasets), and get the output energies gathered in the final echo of variables :
etotal1 -1.0368223891E+00
etotal2 -1.0538645432E+00
etotal3 -1.0674504850E+00
etotal4 -1.0781904896E+00
etotal5 -1.0865814785E+00
etotal6 -1.0930286804E+00
etotal7 -1.0978628207E+00
etotal8 -1.1013539124E+00
etotal9 -1.1037224213E+00
etotal10 -1.1051483730E+00
etotal11 -1.1057788247E+00
etotal12 -1.1057340254E+00
etotal13 -1.1051125108E+00
etotal14 -1.1039953253E+00
etotal15 -1.1024495225E+00
etotal16 -1.1005310615E+00
etotal17 -1.0982871941E+00
etotal18 -1.0957584182E+00
etotal19 -1.0929800578E+00
etotal20 -1.0899835224E+00
etotal21 -1.0867972868E+00
You might try to plot these data. The minimum of energy in the above list is clearly between dataset 11 and 12, that is :
xcart11 -7.5000000000E-01 0.0000000000E+00 0.0000000000E+00
7.5000000000E-01 0.0000000000E+00 0.0000000000E+00
xcart12 -7.7500000000E-01 0.0000000000E+00 0.0000000000E+00
7.7500000000E-01 0.0000000000E+00 0.0000000000E+00
corresponding to a distance of H atoms between 1.5 Bohr and 1.55 Bohr. The forces vanish also between 1.5 Bohr and 1.55 Bohr :
fcart11 -5.4963645520E-03 0.0000000000E+00 0.0000000000E+00
5.4963645520E-03 0.0000000000E+00 0.0000000000E+00
fcart12 6.9585355532E-03 0.0000000000E+00 0.0000000000E+00
-6.9585355532E-03 0.0000000000E+00 0.0000000000E+00
From these two values, using a linear interpolation, one get the optimal value
of 1.522 Bohr .
Note that the number of SCF cycles drops from 7 to 5 when the wavefunctions
are read from the previous dataset.
1.3.a The other methodology is based on an automatic computation of
the minimum.
There are different algorithms to do that.
See the input variable "ionmov",
with values 2, 3 and 7. In the present case, with only one degree of freedom
to be optimized, the best choice
is ionmov 3 .
You have also to define the maximal number of timesteps for this optimization. Set the input variable "ntime" to 10, it will be largely enough. For the stopping criterion "tolmxf", use the reasonable value of 5.0d-4 Ha/Bohr. This defines the force treshhold to consider that the geometry is converged. The code will stop if the residual forces are below that value before reaching "ntime".
It is also worth to change the stopping criterion for the SCF cycle, in order to be sure that the forces generated for each trial interatomic distance are sufficiently converged. So, change "toldfe" in "toldff", and set the latter input variable to ten times smaller than "tolmxf". The input file ~ABINIT/Tutorial/t13.in is an example of file that will do the job, while ~ABINIT/Tutorial/Refs/t13.out is an example of output file. If you decide to use these files, do not forget to change the file names in the t1x.files" file ... So, you run the code with your input file (it should take less than one minute), examine quietly this file (which is much smaller than the t12.out file), and get some significant output data gathered in the final echo of variables :
etotal -1.1058360628E+00
fcart 1.8438010986E-04 0.0000000000E+00 0.0000000000E+00
-1.8438010986E-04 0.0000000000E+00 0.0000000000E+00
...
xcart -7.6091430410E-01 0.0000000000E+00 0.0000000000E+00
7.6091430410E-01 0.0000000000E+00 0.0000000000E+00
According to these data (see xcart), the optimal interatomic distance is about 1.520 Bohr, in good agreement with the estimation of t12.out . If you have time (this is to be done at home), you might try to change the stopping criteria, and redo the calculation, to see the level of convergence of the interatomic distance.
Note that the final value of fcart in your run might differ slightly from the one shown above (less than one percent change). Such a fluctuation is quite often observed for a value converging to zero (remember, we ask the code to determine the equilibrium geometry, that is, forces should be zero) when the same computation is done on different platforms.
1.4.a. We start from the optimized interatomic distance 1.522 Bohr, and make a run at fixed geometry. The input variable "prtden" must be set to 1. To understand correctly the content of the "prtden" description, it is worth to read a much more detailed description of the "files" file, in section 4 of the abinis_help file.
1.4.b The input file ~ABINIT/Tutorial/t14.in is an example of input file for a run that will print a density. If you decide to use this file, do not forget to change the file names in t1x.files. The run will take a few seconds.
The density will be output in the t1xo_DEN file. Try to edit it ... No luck ! This file is unformatted, not written using the ASCII code. Even if you cannot read it, its description is provided in the abinis_help. It contains first a header, then the density numbers. The description of the header is presented in section 6.4 of the abinis_help file, while the body of the _DEN file in presented in section 6.5. It is the appropriate time to read also the description of the potential files and wavefunctions files, as these files contain the same header as the density file, see sections 6.6 and 6.7
1.4.c Such a density file can be read by ABINIT, to restart a calculation
(see the input variable iscf,
when its value is -2), but more usually, by an utility called "cut3d". This
utility is available in the ABINIT package.
You might try to use it now, to generate two-dimensional
cuts in the density, and visualize the charge density contours.
Read the corresponding Cut3D help file.
Then, try to run cut3d to analyse t1xo_DEN.
You should first try to translate the unformatted density data to indexed formatted data,
by using option 6 in the adequate menu. Save the indexed formatted data to file t1xo_DEN_indexed.
Then, edit this file, to have an idea of the content of the _DEN files.
For further treatment, you might choose to select another option than 6. In particular,
if you have access to MATLAB, choose option 5. With minor modifications (set ngx=ngy=ngz to 30)
you will be able to use the file dim.m present in ~ABINIT/Tutorial to visualize the
3-Dimensional isosurfaces. Another option might be to use the XCrysDen software,
for which you need to use option 9.
1.5.a
The atomisation energy is the energy needed to separate a molecule
in its constituant atoms, each being neutral.
In the present case, one must compute first the total energy of an isolated
hydrogen atom. The atomisation energy will be the difference between the total
energy of H2 and twice the total energy of H .
There are some subtleties in the calculation of an isolated atom.
The input file ~ABINIT/Tutorial/t15.in is an example of file that will do the job, while ~ABINIT/Tutorial/Refs/t15.out is an example of output file. If you decide to use the t15.in file, do not forget to change the file names in the t1x.files file. The run lasts a few seconds.
You should read the output file, and note the tiny differences related with the spin-polarisation :
Eigenvalues (hartree) for nkpt= 1 k points, SPIN UP: kpt# 1, nband= 1, wtk= 1.00000, kpt= 0.0000 0.0000 0.0000 (reduced coord) -0.26422 Eigenvalues (hartree) for nkpt= 1 k points, SPIN DOWN: kpt# 1, nband= 1, wtk= 1.00000, kpt= 0.0000 0.0000 0.0000 (reduced coord) -0.11117
,Min spin pol zeta= 1.0000E+00 at reduced coord. 0.0000 0.0000 0.0000 , next min= 1.0000E+00 at reduced coord. 0.0333 0.0000 0.0000 ,Max spin pol zeta= 1.0000E+00 at reduced coord. 0.9667 0.9667 0.9667 , next max= 1.0000E+00 at reduced coord. 0.9333 0.9667 0.9667
The total energy is
etotal -4.7010531340E-01
while the total energy of the H2 molecule is (see test 13):
etotal -1.1058360629E+00
The atomisation energy is thus 0.1656 Ha.
At this stage, we can compare our results :
What is wrong ??
Well, are you sure that the input parameters that we did not discuss
are correct ? These are :
We will see in the next lesson how to address the choice of these parameters (except the pseudopotential).
Answers to the questions, section 1.1.j :
NOTE : there might be numerical differences, from platform to platform, in the quoted results !
Q1. 7 SCF cycles were needed :
iter Etot(hartree) deltaE(h) residm vres2 diffor maxfor
ETOT 1 -1.1013391225241 -1.101E+00 4.220E-04 8.396E+00 2.458E-02 2.458E-02
ETOT 2 -1.1034123727266 -2.073E-03 4.367E-09 1.668E+00 8.602E-03 3.318E-02
ETOT 3 -1.1037064870489 -2.941E-04 1.836E-05 3.207E-01 4.922E-03 3.810E-02
ETOT 4 -1.1037182046373 -1.172E-05 1.090E-07 8.675E-02 3.620E-04 3.774E-02
ETOT 5 -1.1037224013769 -4.197E-06 1.436E-07 1.829E-04 3.593E-04 3.738E-02
ETOT 6 -1.1037224209642 -1.959E-08 1.123E-09 1.445E-05 3.106E-05 3.741E-02
ETOT 7 -1.1037224213176 -3.534E-10 6.528E-12 8.113E-07 4.102E-06 3.741E-02
At SCF step 7, etot is converged :
for the second time, diff in etot= 3.534E-10 < toldfe= 1.000E-06
Q2. Yes, the energy is more converged than toldfe, since the stopping criterion asked for the difference between successive evaluations of the energy to be smaller than toldfe twice in a row, while the evolution of the energy is nice, and always decreasing by smaller and smaller amounts.
Q3. These values are :
cartesian forces (hartree/bohr) at end:
1 -0.03740515236097 0.00000000000000 0.00000000000000
2 0.03740515236097 0.00000000000000 0.00000000000000
frms,max,avg= 2.1595875E-02 3.7405152E-02 0.000E+00 0.000E+00 0.000E+00 h/b
On the first atom (located at -0.7 0 0 in cartesian coordinates, in Bohr), the
force vector is pointing in the minus x direction, and in the plus x direction
for the second atom located at +0.7 0 0 .
The H2 molecule would like to expand ...
Q4. The eigenvalues (in Hartree) are mentioned at the lines
Eigenvalues (hartree) for nkpt= 1 k points: kpt# 1, nband= 2, wtk= 1.00000, kpt= 0.0000 0.0000 0.0000 (reduced coord) -0.36526 -0.01379
As mentioned in the abinis_help file, the
absolute value of eigenenergies is not meaningful. Only differences of eigenenergies,
as well as differences with the potential.
The difference is 0.35147 Hartree, that is 9.564 eV .
Moreover, remember that Kohn-Sham eigenenergies are formally NOT connected to
experimental excitation energies !
(Well, more is to be said later about this ...).
Q5. The maximum electronic density in electron per Bohr cube is reached at the mid-point between the two H atoms :
,Max el dens= 2.6907E-01 el/bohr^3 at reduced coord. 0.0000 0.0000 0.0000