2823
2808
%endblock Zmatrix
2826
The sections that can be used within the Zmatrix block are as follows:
2828
Firstly, all atomic positions must be specified within either a
2829
``\texttt{molecule}'' block or a ``\texttt{cartesian}'' block. Any
2830
atoms subject to constraints more complicated than ``do not change
2831
this coordinate of this atom'' must be specified within a
2832
``\texttt{molecule}'' block.
2834
\item \texttt{molecule}:
2836
There must be one of these blocks for each independent set of
2837
constrained atoms within the simulation.
2839
This specifies the atoms that make up each molecule and their
2840
geometry. In addition, an option of ``\texttt{fractional}'' or
2841
``\texttt{scaled}'' may be passed, which indicates that distances are
2842
specified in scaled or fractional units. In the absence of such an
2843
option, the distance units are taken to be the value of
2844
``\texttt{ZM.UnitsLength}''.
2846
A line is needed for each atom in the molecule; the format of each
2849
\noindent\texttt{ Nspecies i j k r a t ifr ifa ift}
2851
Here the values \texttt{Nspecies}, \texttt{i}, \texttt{j}, \texttt{k},
2852
\texttt{ifr}, \texttt{ifa}, and \texttt{ift} are integers and
2853
\texttt{r}, \texttt{a}, and \texttt{t} are double precision reals.
2855
For most atoms, \texttt{Nspecies} is the species number of the atom,
2856
\texttt{r} is distance to atom number \texttt{i}, \texttt{a} is the
2857
angle made by the present atom with atoms \texttt{j} and \texttt{i},
2858
while \texttt{t} is the torsional angle made by the present atom with
2859
atoms \texttt{k}, \texttt{j}, and \texttt{i}. The values \texttt{ifr},
2860
\texttt{ifa} and \texttt{ift} are integer flags that indicate whether
2861
\texttt{r}, \texttt{a}, and \texttt{t}, respectively, should be
2862
varied; 0 for fixed, 1 for varying.
2865
The first three atoms in a molecule are a special case. Because there
2866
are insufficient atoms defined to specify a distance/angle/torsion,
2867
the values are set differently. For atom 1, \texttt{r}, \texttt{a},
2868
and \texttt{t}, are the Cartesian coordinates of the atom. For the
2869
second atom, \texttt{r}, \texttt{a}, and \texttt{t} are the
2870
coordinates in spherical form of the second atom relative to the
2871
first: first the radius, then the polar angle (angle between the
2872
$z$-axis and the displacement vector) and then the azimuthal angle
2873
(angle between the $x$-axis and the projection of the displacement
2874
vector on the $x$-$y$ plane). Finally, for the third atom, the numbers
2875
take their normal form, but the torsional angle is defined relative to
2876
a notional atom 1 unit in the z-direction above the atom \texttt{j}.
2878
Secondly. blocks of atoms all of which are subject to the simplest of
2879
constraints may be specified in one of the following three ways,
2880
according to the units used to specify their coordinates:
2883
\item \texttt{cartesian}: This section specifies a block of atoms
2884
whose coordinates are to be specified in Cartesian coordinates. Again,
2885
an option of ``\texttt{fractional}'' or ``\texttt{scaled}'' may be
2886
added, to specify the units used; and again, in their absence, the
2887
value of ``\texttt{ZM.UnitsLength}'' is taken.
2889
The format of each atom in the block will look like:
2891
\noindent\texttt{ Nspecies x y z ix iy iz}
2893
Here \texttt{Nspecies}, \texttt{ix}, \texttt{iy}, and \texttt{iz} are
2894
integers and \texttt{x}, \texttt{y}, \texttt{z} are
2895
reals. \texttt{Nspecies} is the species number of the atom being
2896
specified, while \texttt{x}, \texttt{y}, and \texttt{z} are the
2897
Cartesian coordinates of the atom in whichever units are being
2898
used. The values \texttt{ix}, \texttt{iy} and \texttt{iz} are integer
2899
flags that indicate whether the \texttt{x}, \texttt{y}, and \texttt{z}
2900
coordinates, respectively, should be varied or not. A value of 0
2901
implies that the coordinate is fixed, while 1 implies that it should
2902
be varied. \textbf{NOTE}: When performing ``variable cell''
2903
optimization while using a Zmatrix format for input, the algorithm
2904
will not work if some of the coordinates of an atom in a
2905
\texttt{cartesian} block are variables and others are not (i.e.,
2906
\texttt{ix iy iz} above must all be 0 or 1). This will be fixed in
2907
future versions of the program.
2909
A Zmatrix block may also contain the following, additional, sections, which
2910
are designed to make it easier to read.
2912
\item \texttt{constants}: Instead of specifying a numerical value, it
2913
is possible to specify a symbol within the above geometry
2914
definitions. This section allows the user to define the value of the
2915
symbol as a constant. The format is just a symbol followed by the
2918
\noindent\texttt{ HOH 104.4}
2920
\item \texttt{variables}: Instead of specifying a numerical value, it
2921
is possible to specify a symbol within the above geometry
2922
definitions. This section allows the user to define the value of the
2923
symbol as a variable. The format is just a symbol followed by the
2926
\noindent\texttt{ HO1 0.956997}
2928
Finally, constraints must be specified in a \texttt{constraints} block.
2930
\item \texttt{constraint} This sub-section allows the user to create
2931
constraints between symbols used in a Z-matrix:
2936
Here var1 and var2 are text symbols for two quantities in the Z-matrix
2937
definition, and A and B are real numbers. The variables are related by
2938
$var1 = A*var2 + B$.
2940
An example of a Z-matrix input for a benzene molecule over a metal surface is:
2811
The sections that can be used within the Zmatrix block are as
2814
Firstly, all atomic positions must be specified within either a
2815
``\texttt{molecule}'' block or a ``\texttt{cartesian}'' block. Any
2816
atoms subject to constraints more complicated than ``do not change
2817
this coordinate of this atom'' must be specified within a
2818
``\texttt{molecule}'' block.
2823
There must be one of these blocks for each independent set of
2824
constrained atoms within the simulation.
2826
This specifies the atoms that make up each molecule and their
2827
geometry. In addition, an option of ``\texttt{fractional}'' or
2828
``\texttt{scaled}'' may be passed, which indicates that distances are
2829
specified in scaled or fractional units. In the absence of such an
2830
option, the distance units are taken to be the value of
2831
``\texttt{ZM.UnitsLength}''.
2833
A line is needed for each atom in the molecule; the format of each
2836
Nspecies i j k r a t ifr ifa ift
2839
Here the values \texttt{Nspecies}, \texttt{i}, \texttt{j}, \texttt{k},
2840
\texttt{ifr}, \texttt{ifa}, and \texttt{ift} are integers and
2841
\texttt{r}, \texttt{a}, and \texttt{t} are double precision reals.
2843
For most atoms, \texttt{Nspecies} is the species number of the atom,
2844
\texttt{r} is distance to atom number \texttt{i}, \texttt{a} is the
2845
angle made by the present atom with atoms \texttt{j} and \texttt{i},
2846
while \texttt{t} is the torsional angle made by the present atom with
2847
atoms \texttt{k}, \texttt{j}, and \texttt{i}. The values \texttt{ifr},
2848
\texttt{ifa} and \texttt{ift} are integer flags that indicate whether
2849
\texttt{r}, \texttt{a}, and \texttt{t}, respectively, should be
2850
varied; 0 for fixed, 1 for varying.
2853
The first three atoms in a molecule are a special case. Because there
2854
are insufficient atoms defined to specify a distance/angle/torsion,
2855
the values are set differently. For atom 1, \texttt{r}, \texttt{a},
2856
and \texttt{t}, are the Cartesian coordinates of the atom. For the
2857
second atom, \texttt{r}, \texttt{a}, and \texttt{t} are the
2858
coordinates in spherical form of the second atom relative to the
2859
first: first the radius, then the polar angle (angle between the
2860
$z$-axis and the displacement vector) and then the azimuthal angle
2861
(angle between the $x$-axis and the projection of the displacement
2862
vector on the $x$-$y$ plane). Finally, for the third atom, the numbers
2863
take their normal form, but the torsional angle is defined relative to
2864
a notional atom 1 unit in the z-direction above the atom \texttt{j}.
2866
Secondly. blocks of atoms all of which are subject to the simplest of
2867
constraints may be specified in one of the following three ways,
2868
according to the units used to specify their coordinates:
2871
This section specifies a block of atoms
2872
whose coordinates are to be specified in Cartesian coordinates. Again,
2873
an option of ``\texttt{fractional}'' or ``\texttt{scaled}'' may be
2874
added, to specify the units used; and again, in their absence, the
2875
value of ``\texttt{ZM.UnitsLength}'' is taken.
2877
The format of each atom in the block will look like:
2879
Nspecies x y z ix iy iz
2882
Here \texttt{Nspecies}, \texttt{ix}, \texttt{iy}, and \texttt{iz} are
2883
integers and \texttt{x}, \texttt{y}, \texttt{z} are
2884
reals. \texttt{Nspecies} is the species number of the atom being
2885
specified, while \texttt{x}, \texttt{y}, and \texttt{z} are the
2886
Cartesian coordinates of the atom in whichever units are being
2887
used. The values \texttt{ix}, \texttt{iy} and \texttt{iz} are integer
2888
flags that indicate whether the \texttt{x}, \texttt{y}, and \texttt{z}
2889
coordinates, respectively, should be varied or not. A value of 0
2890
implies that the coordinate is fixed, while 1 implies that it should
2891
be varied. \textbf{NOTE}: When performing ``variable cell''
2892
optimization while using a Zmatrix format for input, the algorithm
2893
will not work if some of the coordinates of an atom in a
2894
\texttt{cartesian} block are variables and others are not (i.e.,
2895
\texttt{ix iy iz} above must all be 0 or 1). This will be fixed in
2896
future versions of the program.
2898
A Zmatrix block may also contain the following, additional, sections, which
2899
are designed to make it easier to read.
2903
Instead of specifying a numerical value, it is possible to specify
2904
a symbol within the above geometry definitions. This section
2905
allows the user to define the value of the symbol as a
2906
constant. The format is just a symbol followed by the value:
2913
Instead of specifying a numerical value, it is possible to specify
2914
a symbol within the above geometry definitions. This section
2915
allows the user to define the value of the symbol as a
2916
variable. The format is just a symbol followed by the value:
2921
Finally, constraints must be specified in a \fdf*{constraints} block.
2924
\option[constraint]%
2925
This sub-section allows the user to create
2926
constraints between symbols used in a Z-matrix:
2931
Here \fdf*{var1} and \fdf*{var2} are text symbols for two
2932
quantities in the Z-matrix definition, and $A and $B are real
2933
numbers. The variables are related by $\fdf*{var1} = A*\fdf*{var2}
2938
An example of a Z-matrix input for a benzene molecule over a metal surface is:
2944
2 0 0 0 xm1 ym1 zm1 0 0 0
2945
2 1 0 0 CC 90.0 60.0 0 0 0
2946
2 2 1 0 CC CCC 90.0 0 0 0
2947
2 3 2 1 CC CCC 0.0 0 0 0
2948
2 4 3 2 CC CCC 0.0 0 0 0
2949
2 5 4 3 CC CCC 0.0 0 0 0
2950
1 1 2 3 CH CCH 180.0 0 0 0
2951
1 2 1 7 CH CCH 0.0 0 0 0
2952
1 3 2 8 CH CCH 0.0 0 0 0
2953
1 4 3 9 CH CCH 0.0 0 0 0
2954
1 5 4 10 CH CCH 0.0 0 0 0
2955
1 6 5 11 CH CCH 0.0 0 0 0
2957
3 0.000000 0.000000 0.000000 0 0 0
2958
3 0.333333 0.000000 0.000000 0 0 0
2959
3 0.666666 0.000000 0.000000 0 0 0
2960
3 0.000000 0.500000 0.000000 0 0 0
2961
3 0.333333 0.500000 0.000000 0 0 0
2962
3 0.666666 0.500000 0.000000 0 0 0
2963
3 0.166667 0.250000 0.050000 0 0 0
2964
3 0.500000 0.250000 0.050000 0 0 0
2965
3 0.833333 0.250000 0.050000 0 0 0
2966
3 0.166667 0.750000 0.050000 0 0 0
2967
3 0.500000 0.750000 0.050000 0 0 0
2968
3 0.833333 0.750000 0.050000 0 0 0
2969
3 0.000000 0.000000 0.100000 0 0 0
2970
3 0.333333 0.000000 0.100000 0 0 0
2971
3 0.666666 0.000000 0.100000 0 0 0
2972
3 0.000000 0.500000 0.100000 0 0 0
2973
3 0.333333 0.500000 0.100000 0 0 0
2974
3 0.666666 0.500000 0.100000 0 0 0
2975
3 0.166667 0.250000 0.150000 0 0 0
2976
3 0.500000 0.250000 0.150000 0 0 0
2977
3 0.833333 0.250000 0.150000 0 0 0
2978
3 0.166667 0.750000 0.150000 0 0 0
2979
3 0.500000 0.750000 0.150000 0 0 0
2980
3 0.833333 0.750000 0.150000 0 0 0
2990
xm1 CC -1.0 3.903229
2994
Here the species 1, 2 and 3 represent H, C, and the metal of the
2995
surface, respectively.
2997
(Note: the above example shows the usefulness of symbolic names
2998
for the relevant coordinates, in particular for those which are
2999
allowed to vary. The current output options for Zmatrix information
3000
work best when this approach is taken. By using a ``fixed'' symbolic
3001
Zmatrix block and specifying the actual coordinates in a ``variables''
3002
section, one can monitor the progress of the optimization and
3003
easily reconstruct the coordinates of intermediate steps in the
3006
\textit{Use:} Specifies the geometry of the system according to a Z-matrix
3007
\textit{Default value:} Geometry is not specified using a Z-matrix
3009
\item[\textbf{ZM.UnitsLength}] (\textit{length}):
3010
\index{ZM.UnitsLength@\textbf{ZM.UnitsLength}}
3011
Parameter that specifies the units of length used during Z-matrix input.
3013
\textit{Use:} This option allows the user to chose between inputing distances
3014
in Bohr or Angstroms within the Z-matrix data block.
3016
\textit{Default value:} \texttt{Bohr}
3018
\item[\textbf{ZM.UnitsAngle}] (\textit{angle}):
3019
\index{ZM.UnitsAngle@\textbf{ZM.UnitsAngle}}
3020
Parameter that specifies the units of angles used during Z-matrix input.
3022
\textit{Use:} This option allows the user to chose between inputing angles
3023
in radians or degrees within the Z-matrix data block.
3025
\textit{Default value:} \texttt{rad}
2942
2 0 0 0 xm1 ym1 zm1 0 0 0
2943
2 1 0 0 CC 90.0 60.0 0 0 0
2944
2 2 1 0 CC CCC 90.0 0 0 0
2945
2 3 2 1 CC CCC 0.0 0 0 0
2946
2 4 3 2 CC CCC 0.0 0 0 0
2947
2 5 4 3 CC CCC 0.0 0 0 0
2948
1 1 2 3 CH CCH 180.0 0 0 0
2949
1 2 1 7 CH CCH 0.0 0 0 0
2950
1 3 2 8 CH CCH 0.0 0 0 0
2951
1 4 3 9 CH CCH 0.0 0 0 0
2952
1 5 4 10 CH CCH 0.0 0 0 0
2953
1 6 5 11 CH CCH 0.0 0 0 0
2955
3 0.000000 0.000000 0.000000 0 0 0
2956
3 0.333333 0.000000 0.000000 0 0 0
2957
3 0.666666 0.000000 0.000000 0 0 0
2958
3 0.000000 0.500000 0.000000 0 0 0
2959
3 0.333333 0.500000 0.000000 0 0 0
2960
3 0.666666 0.500000 0.000000 0 0 0
2961
3 0.166667 0.250000 0.050000 0 0 0
2962
3 0.500000 0.250000 0.050000 0 0 0
2963
3 0.833333 0.250000 0.050000 0 0 0
2964
3 0.166667 0.750000 0.050000 0 0 0
2965
3 0.500000 0.750000 0.050000 0 0 0
2966
3 0.833333 0.750000 0.050000 0 0 0
2967
3 0.000000 0.000000 0.100000 0 0 0
2968
3 0.333333 0.000000 0.100000 0 0 0
2969
3 0.666666 0.000000 0.100000 0 0 0
2970
3 0.000000 0.500000 0.100000 0 0 0
2971
3 0.333333 0.500000 0.100000 0 0 0
2972
3 0.666666 0.500000 0.100000 0 0 0
2973
3 0.166667 0.250000 0.150000 0 0 0
2974
3 0.500000 0.250000 0.150000 0 0 0
2975
3 0.833333 0.250000 0.150000 0 0 0
2976
3 0.166667 0.750000 0.150000 0 0 0
2977
3 0.500000 0.750000 0.150000 0 0 0
2978
3 0.833333 0.750000 0.150000 0 0 0
2988
xm1 CC -1.0 3.903229
2992
Here the species 1, 2 and 3 represent H, C, and the metal of the
2993
surface, respectively.
2995
(Note: the above example shows the usefulness of symbolic names
2996
for the relevant coordinates, in particular for those which are
2997
allowed to vary. The current output options for Zmatrix information
2998
work best when this approach is taken. By using a ``fixed'' symbolic
2999
Zmatrix block and specifying the actual coordinates in a ``variables''
3000
section, one can monitor the progress of the optimization and
3001
easily reconstruct the coordinates of intermediate steps in the
3006
\begin{fdfentry}{ZM!UnitsLength}[string]<Bohr>
3008
Parameter that specifies the units of length used during Z-matrix
3011
Specify \fdf*{Bohr} or \fdf*{Ang} for the corresponding unit of length.
3015
\begin{fdfentry}{ZM!UnitsAngle}[string]<rad>
3017
Parameter that specifies the units of angles used during Z-matrix input.
3019
Specify \fdf*{rad} or \fdf*{deg} for the corresponding unit of angle.
3031
3024
\subsubsection{Output of structural information}
3078
3068
noted also in fdf form.
3081
Note that the geometry reported is the last one for which forces and
3082
stresses were computed.
3084
\item\textbf{NEXT\_ITER.UCELL.ZMATRIX file:}
3085
\index{NEXT\_ITER.UCELL.ZMATRIX}@{\textit{NEXT\_ITER.UCELL.ZMATRIX}}
3087
A file with the same format as \texttt{OUT.UCELL.ZMATRIX} but with
3088
a possibly updated geometry.
3090
\item The coordinates can be also accumulated
3091
in the \textit{Systemlabel}.MD or \textit{Systemlabel}.MDX files
3092
depending on \textbf{WriteMDhistory}.
3094
\item Additionally, several optional formats are supported:
3100
\item[\textbf{WriteCoorXmol}] (\textit{logical}):
3101
\index{WriteCoorXmol@\textbf{WriteCoorXmol}}\index{XMol@\textsc{XMol}}
3102
\index{JMol@\textsc{JMol}}
3103
If \texttt{.true.} it originates the writing of an extra file
3104
named \textit{SystemLabel}\texttt{.xyz} containing the final atomic
3105
coordinates in a format directly readable by \textsc{XMol}.\footnote{XMol
3106
is under \copyright\ copyright of Research Equipment Inc., dba Minnesota
3107
Supercomputer Center Inc.} Coordinates come out in {\AA}ngstr\"om
3108
independently of what specified in \textbf{AtomicCoordinatesFormat} and
3109
in \textbf{AtomCoorFormatOut}. There is a present \textsc{Java} implementation
3110
of \textsc{XMol} called \textsc{JMol}.
3112
\textit{Default value:} \texttt{.false.}
3115
\item[\textbf{WriteCoorCerius}] (\textit{logical}):
3116
\index{WriteCoorCerius@\textbf{WriteCoorCerius}}\index{Cerius2@\textsc{Cerius2}}
3117
If \texttt{.true.} it originates the writing of an extra file
3118
named \textit{SystemLabel}\texttt{.xtl} containing the final atomic
3119
coordinates in a format directly readable by
3120
\method{Cerius}.\footnote{\method{Cerius} is under \copyright\
3121
copyright of Molecular Simulations Inc.}
3122
Coordinates come out in
3123
\texttt{Fractional} format (the same as \texttt{ScaledByLatticeVectors})
3124
independently of what specified in \textbf{AtomicCoordinatesFormat} and
3125
in \textbf{AtomCoorFormatOut}.
3126
If negative coordinates are to be avoided, it has to be
3127
done from the start by shifting all the coordinates rigidly
3128
to have them positive, by using \textbf{AtomicCoordinatesOrigin}.
3129
See the \textsc{Sies2arc}\index{Sies2arc@\textsc{Sies2arc}} utility in the Util/
3130
directory for generating .arc files for CERIUS animation.
3132
\textit{Default value:} \texttt{.false.}
3135
\item[\textbf{WriteMDXmol}] (\textit{logical}):
3136
\index{WriteMDXmol@\textbf{WriteMDXmol}}\index{XMol@\textsc{XMol}}
3137
If \texttt{.true.} it causes the writing of an extra file
3138
named \textit{SystemLabel}\texttt{.ANI} containing all the atomic
3139
coordinates of the simulation in a format directly readable by
3140
\textsc{XMol} for animation.\index{animation} Coordinates come out in
3141
{\AA}ngstr\"om independently of what is specified in
3142
\textbf{AtomicCoordinatesFormat} and in \textbf{AtomCoorFormatOut}.
3143
This file is accumulative even for different runs.
3145
There is an alternative for animation by generating a .arc file for
3146
CERIUS. It is through the
3147
\textsc{Sies2arc}\index{Sies2arc@\program{Sies2arc}} postprocessing
3148
utility in the Util/ directory, and it
3149
requires the coordinates to be accumulated in the output file, i.e.,
3150
WriteCoorStep = \texttt{.true.}
3152
\textit{Default value:} \texttt{.false.}
3154
Note change with respect to previous versions. This option is no
3155
longer coupled to \textbf{WriteCoorStep}.
3071
Note that the geometry reported is the last one for which forces and
3072
stresses were computed.
3074
\item\file{NEXT\_ITER.UCELL.ZMATRIX}:%
3075
A file with the same format as \file{OUT.UCELL.ZMATRIX} but with
3076
a possibly updated geometry.
3078
\item The coordinates can be also accumulated
3079
in the \sysfile{MD} or \sysfile{MDX} files
3080
depending on \fdf{WriteMDHistory}.
3082
\item Additionally, several optional formats are supported:
3083
\begin{fdflogicalF}{WriteCoorXmol}
3084
\index{JMol@\textsc{JMol}}
3085
\index{XMol@\textsc{XMol}}
3086
\index{Molden@\textsc{Molden}}
3088
If \fdftrue\ it originates the writing of an extra file named
3089
\sysfile{xyz} containing the final atomic coordinates in a format
3090
directly readable by \method{XMol}.\footnote{XMol is under
3091
\copyright\ copyright of Research Equipment Inc., dba
3092
Minnesota Supercomputer Center Inc.} Coordinates come out in
3093
\AA ngstr\"om independently of what specified in
3094
\fdf{AtomicCoordinatesFormat} and in
3095
\fdf{AtomCoorFormatOut}. There is a present \method{Java}
3096
implementation of \method{XMol} called \method{JMol}.
3100
\begin{fdflogicalF}{WriteCoorCerius}
3101
\index{Cerius2@\textsc{Cerius2}}
3103
If \fdftrue it originates the writing of an extra file named
3104
\sysfile{xtl} containing the final atomic coordinates in a format
3105
directly readable by \method{Cerius}.\footnote{\method{Cerius} is
3106
under \copyright\ copyright of Molecular Simulations Inc.}
3107
Coordinates come out in \fdf*{Fractional} format (the same as
3108
\fdf*{ScaledByLatticeVectors}) independently of what specified
3109
in \fdf{AtomicCoordinatesFormat} and in
3110
\fdf{AtomCoorFormatOut}. If negative coordinates are to be
3111
avoided, it has to be done from the start by shifting all the
3112
coordinates rigidly to have them positive, by using
3113
\fdf{AtomicCoordinatesOrigin}. See the
3114
\program{Sies2arc}\index{Sies2arc@\textsc{Sies2arc}} utility in the
3115
\program{Util/} directory for generating \sysfile*{.arc} files for CERIUS animation.
3119
\begin{fdflogicalF}{WriteMDXmol}
3120
\index{XMol@\textsc{XMol}}
3121
\index{Molden@\textsc{Molden}}
3123
If \fdftrue\ it causes the writing of an extra file
3124
named \sysfile{ANI} containing all the atomic
3125
coordinates of the simulation in a format directly readable by
3126
\method{XMol} for animation.\index{animation} Coordinates come out in
3127
\AA ngstr\"om independently of what is specified in
3128
\fdf{AtomicCoordinatesFormat} and in \fdf{AtomCoorFormatOut}.
3129
This file is accumulative even for different runs.
3131
There is an alternative for animation by generating a \sysfile*{arc} file for
3132
\method{CERIUS}. It is through the
3133
\method{Sies2arc}\index{Sies2arc@\program{Sies2arc}} postprocessing
3134
utility in the Util/ directory, and it
3135
requires the coordinates to be accumulated in the output file, i.e.,
3136
\fdf{WriteCoorStep} \fdftrue.
3161
3144
\subsubsection{Input of structural information from external files}
3163
3146
The structural information can be also read from external files. Note
3164
that the NumberOfAtoms, NumberOfSpecies, and ChemicalSpeciesLabel
3165
options are still mandatory in the fdf file.
3170
\index{MD.UseSaveXV@\textbf{MD.UseSaveXV}} \index{reading saved data!XV}
3171
The logical variable \textbf{MD.UseSaveXV} instructs \siesta\ to
3172
read the atomic positions and velocities stored in file \sysfile{XV} by a previous run.
3174
If the required file does not exist, a warning is printed but the
3175
program does not stop. Overrides \textbf{UseSaveData}, but can be
3176
implicitly set by it.
3178
\textit{Default value:} \texttt{.false.}
3180
\item A .STRUCT\_IN file. The logical fdf variable \textbf{UseStructFile}
3181
(for historical reasons, \textbf{MD.UseStructFile} is also accepted,
3183
controls whether the structural information is read from an external
3184
file of name \textit{SystemLabel}.STRUCT\_IN. If \texttt{.true.}, all
3185
other structural information in the fdf file will be
3186
ignored.\index{UseStructFile@\textbf{UseStructFile}}
3187
\index{Systemlabel.STRUCT\_IN@{\textit{Systemlabel}.STRUCT\_IN}}
3189
The format of the file is implied by the following code:
3147
that \fdf{ChemicalSpeciesLabel} is mandatory in the fdf file.
3149
\begin{fdflogicalF}{MD!UseSaveXV}
3150
\index{reading saved data!XV}
3152
Logical variable which instructs \siesta\ to read the atomic
3153
positions and velocities stored in file \sysfile{XV} by a previous
3156
If the file does not exist, a warning is printed but the
3157
program does not stop. Overrides \fdf{UseSaveData}, but can be
3158
implicitly set by it.
3162
\begin{fdflogicalF}{UseStructFile}
3164
Controls whether the structural information is read from an external
3165
file of name \sysfile{STRUCT\_IN}. If \fdftrue, all other
3166
structural information in the fdf file will be ignored.
3168
The format of the file is implied by the following code:
3191
3169
\begin{verbatim}
3192
3170
read(*,*) ((cell(ixyz,ivec),ixyz=1,3),ivec=1,3) ! Cell vectors, in Angstroms
3279
3255
other specialized grids, see the Macroscopic Polarization and Density
3280
3256
of States sections.
3287
\item[\textbf{kgrid\_cutoff}] (\textit{real length}):
3288
\index{kgrid\_cutoff@\textbf{kgrid\_cutoff}}
3289
Parameter which determines
3290
the fineness of the k-grid used for Brillouin zone sampling.
3291
It is half the length of the smallest lattice vector of the supercell
3292
required to obtain the same sampling precision with a single k point.
3293
Ref: Moreno and Soler, PRB 45, 13891 (1992).
3295
\textit{Use:} If it is zero, only the gamma point is used. The resulting
3296
k-grid is chosen in an optimal way, according to the method of Moreno
3297
and Soler (using an effective supercell which is as spherical as
3298
possible, thus minimizing the number of k-points for a given
3299
precision). The grid is displaced for even numbers of effective mesh
3300
divisions. This parameter is not used if \textbf{kgrid\_Monkhorst\_Pack}
3301
is specified. If the unit cell changes during the calculation (for
3302
example, in a cell-optimization run, the k-point
3303
grid will change accordingly (see \textbf{ChangeKgridInMD} for the case
3304
of variable-cell molecular-dynamics runs, such as Parrinello-Rahman).
3305
This is analogous to the changes in the
3306
real-space grid, whose fineness is specified by an energy cutoff. If
3307
sudden changes in the number of k-points are not desired, then the
3308
Monkhorst-Pack data block should be used instead. In this case there
3309
will be an implicit change in the quality of the sampling as the cell
3310
changes. Both methods should be equivalent for a well-converged
3313
\textit{Default value:} \texttt{0.0 Bohr}
3316
\item[\textbf{kgrid\_Monkhorst\_Pack}] (\textit{data block}):
3317
\index{kgrid\_Monkhorst\_Pack@\textbf{kgrid\_Monkhorst\_Pack}}
3318
Real-space supercell, whose reciprocal unit cell is that of the
3319
k-sampling grid, and grid displacement for each grid coordinate.
3320
Specified as an integer matrix and a real vector:
3323
%block kgrid_Monkhorst_Pack
3258
\begin{fdfentry}{kgrid!Cutoff}[length]<$0.\,\mathrm{Bohr}$>
3260
Parameter which determines the fineness of the k-grid used for
3261
Brillouin zone sampling. It is half the length of the smallest
3262
lattice vector of the supercell required to obtain the same sampling
3263
precision with a single k point. Ref: Moreno and Soler, PRB 45,
3266
\textit{Use:} If it is zero, only the gamma point is used. The resulting
3267
k-grid is chosen in an optimal way, according to the method of Moreno
3268
and Soler (using an effective supercell which is as spherical as
3269
possible, thus minimizing the number of k-points for a given
3270
precision). The grid is displaced for even numbers of effective mesh
3271
divisions. This parameter is not used if \fdf{kgrid!MonkhorstPack}
3272
is specified. If the unit cell changes during the calculation (for
3273
example, in a cell-optimization run, the k-point
3274
grid will change accordingly (see \fdf{ChangeKgridInMD} for the case
3275
of variable-cell molecular-dynamics runs, such as Parrinello-Rahman).
3276
This is analogous to the changes in the
3277
real-space grid, whose fineness is specified by an energy cutoff. If
3278
sudden changes in the number of k-points are not desired, then the
3279
Monkhorst-Pack data block should be used instead. In this case there
3280
will be an implicit change in the quality of the sampling as the cell
3281
changes. Both methods should be equivalent for a well-converged
3286
\begin{fdfentry}{kgrid!MonkhorstPack}[block]<$\Gamma$-point>
3288
Real-space supercell, whose reciprocal unit cell is that of the
3289
k-sampling grid, and grid displacement for each grid coordinate.
3290
Specified as an integer matrix and a real vector:
3293
%block kgrid.MonkhorstPack
3324
3294
Mk(1,1) Mk(2,1) Mk(3,1) dk(1)
3325
3295
Mk(1,2) Mk(2,2) Mk(3,2) dk(2)
3326
3296
Mk(1,3) Mk(2,3) Mk(3,3) dk(3)
3327
%endblock kgrid_Monkhorst_Pack
3330
where \texttt{Mk(j,i)} are integers and \texttt{dk(i)} are usually
3331
either 0.0 or 0.5 (the program will warn the user if the displacements
3332
chosen are not optimal).
3333
The k-grid supercell is defined from \texttt{Mk}
3334
as in block \textbf{SuperCell} above, i.e.:
3335
$KgridSuperCell(ix,i) = \sum_j CELL(ix,j)*Mk(j,i)$.
3336
Note again that the matrix indexes are inverted: each input line
3337
gives the decomposition of a supercell vector in terms of the unit
3341
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{diagon}.
3342
The k-grid supercell is compatible and unrelated
3343
(except for the default value, see below)
3344
with the \textbf{SuperCell} specifier. Both supercells are given in
3345
terms of the CELL specified by the \textbf{LatticeVectors} block.
3346
If \texttt{Mk} is the identity matrix and \texttt{dk}
3347
is zero, only the $\Gamma$ point of the \textbf{unit} cell is used.
3348
Overrides \textbf{kgrid\_cutoff}
3350
\textit{Default value:} $\Gamma$ point of the (super)cell.
3351
(Default used only when \textbf{kgrid\_cutoff} is not defined).
3353
\item[\textbf{ChangeKgridInMD}] (\textit{boolean}):
3354
\index{ChangeKgridInMD@\textbf{ChangeKgridInMD}}
3356
If \texttt{.true.}, the k-point grid is
3357
recomputed at every iteration during MD runs that potentially
3358
change the unit cell: Parrinello-Rahman, Nose-Parrinello-Rahman, and
3359
Anneal. Regardless of the setting of this flag, the k-point grid
3360
is always updated at every iteration of a variable-cell optimization
3361
and after each step in a ``siesta-as-server'' run.
3363
\textit{Default value:} \texttt{.false.} for historical reasons. The
3364
rationale was to avoid sudden jumps in some properties when the
3365
sampling changes, but if the calculation is well-converged there
3366
should be no problems if the update is enabled.
3368
\item[\textbf{TimeReversalSymmetryForKpoints}] (\textit{boolean}):
3369
\index{TimeReversalSymmetryForKpoints@\textbf{TimeReversalSymmetryForKpoints}}
3371
If \texttt{.true.}, the k-points in the BZ generated by the methods above
3372
are paired as (k,-k) and only one member of the pair is retained. This
3373
symmetry is valid in the absence of external magnetic fields or
3374
spin-orbit interaction.
3376
\textit{Default value:} \texttt{.true.} unless the option \textbf{SpinSpiral}
3377
is used. In this case time-reversal-symmetry is broken explicitly.
3300
where \texttt{Mk(j,i)} are integers and \texttt{dk(i)} are usually
3301
either 0.0 or 0.5 (the program will warn the user if the displacements
3302
chosen are not optimal).
3303
The k-grid supercell is defined from \texttt{Mk}
3304
as in block \fdf{SuperCell} above, i.e.:
3305
$KgridSuperCell(ix,i) = \sum_j CELL(ix,j)*Mk(j,i)$.
3306
Note again that the matrix indexes are inverted: each input line
3307
gives the decomposition of a supercell vector in terms of the unit
3311
\textit{Use:} Used only if \fdf{SolutionMethod} \fdf*{diagon}. The
3312
k-grid supercell is compatible and unrelated (except for the default
3313
value, see below) with the \fdf{SuperCell} specifier. Both
3314
supercells are given in terms of the CELL specified by the
3315
\fdf{LatticeVectors} block. If \texttt{Mk} is the identity matrix
3316
and \texttt{dk} is zero, only the $\Gamma$ point of the
3317
\textbf{unit} cell is used. Overrides \fdf{kgrid!Cutoff}
3321
\begin{fdflogicalF}{ChangeKgridInMD}
3323
If \fdftrue, the $k$-point grid is recomputed at every
3324
iteration during MD runs that potentially change the unit cell:
3325
Parrinello-Rahman, Nose-Parrinello-Rahman, and Anneal. Regardless of
3326
the setting of this flag, the k-point grid is always updated at
3327
every iteration of a variable-cell optimization and after each step
3328
in a ``siesta-as-server'' run.
3330
It is defaulted to \fdffalse\ for historical reasons. The rationale
3331
was to avoid sudden jumps in some properties when the sampling
3332
changes, but if the calculation is well-converged there should be no
3333
problems if the update is enabled.
3338
\begin{fdflogicalT}{TimeReversalSymmetryForKpoints}
3340
If \fdftrue, the $k$-points in the BZ generated by the methods above
3341
are paired as ($k$,$-k$) and only one member of the pair is
3342
retained. This symmetry is valid in the absence of external magnetic
3343
fields or spin-orbit interaction.
3345
Note this defaults to \fdffalse\ if \fdf{Spin} is
3346
\fdf*{non-collinear}, \fdf*{spin-orbit} or \fdf{Spin!Spiral} is used.
3381
3351
\subsubsection{Output of k-point information}
3383
3352
\index{output!grid $\vec k$ points}
3384
3354
The coordinates of the $\vec k$ points used in the sampling
3385
are always stored in the file SystemLabel.KP.
3390
\item[\textbf{WriteKpoints}] (\textit{logical}):
3391
\index{WriteKpoints@\textbf{WriteKpoints}}\index{output!grid $\vec k$ points}
3392
If \texttt{.true.} it writes the coordinates of the $\vec k$ vectors
3393
used in the grid for $k$-sampling, into the main output file.
3395
\textit{Default value:} \texttt{.false.} (see \textbf{LongOutput})
3355
are always stored in the file \sysfile{KP}.
3357
\begin{fdflogicalF}{WriteKpoints}
3358
\index{output!grid $\vec k$ points}
3360
If \fdftrue\ it writes the coordinates of the $\vec k$ vectors used
3361
in the grid for $k$-sampling, into the main output file.
3363
Default depends on \fdf{LongOutput}.
3400
3369
\subsection{Exchange-correlation functionals}
3407
\item[\textbf{XC.functional}] (\textit{string}):
3408
\index{XC.functional@\textbf{XC.functional}}
3409
Exchange-correlation functional type. May be \texttt{LDA}
3410
(local density approximation, equivalent to \texttt{LSD}),
3411
\texttt{GGA} (Generalized Gradient Approximation), or
3412
\texttt{VDW} (van der Waals).
3414
\textit{Use:} Spin polarization is defined by SpinPolarized label for
3415
both \texttt{LDA} and \texttt{GGA}. There is no difference between \texttt{LDA}
3418
\textit{Default value:} \texttt{LDA}
3421
\item[\textbf{XC.authors}] (\textit{string}):
3422
\index{XC.authors@\textbf{XC.authors}}
3423
Particular parametrization of the
3424
exchange-correlation functional. Options are:
3426
\item \texttt{CA} (equivalent to \texttt{PZ}): \index{CA} \index{PZ}
3427
(Spin) local density approximation (LDA/LSD). \index{LDA} \index{LSD}
3428
Quantum Monte Carlo calculation of the homogeneous electron gas
3429
by D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. \textbf{45},566 (1980),
3430
as parametrized by J. P. Perdew and A. Zunger, Phys. Rev B \textbf{23}, 5075 (1981)
3431
\item \texttt{PW92}: \index{PW92}
3432
LDA/LSD, as parametrized by
3433
J. P. Perdew and Y. Wang, Phys. Rev B, \textbf{45}, 13244 (1992)
3434
\item \texttt{PW91}: \index{PW91}
3435
Generalized gradients approximation (GGA) \index{GGA} of Perdew and Wang.
3436
Ref: P\&W, J. Chem. Phys., \textbf{100}, 1290 (1994)
3437
\item \texttt{PBE}: \index{PBE}
3438
GGA of J. P. Perdew, K. Burke and M. Ernzerhof,
3439
Phys. Rev. Lett. \textbf{77}, 3865 (1996)
3440
\item \texttt{revPBE}: \index{revPBE}
3441
Modified GGA-PBE functional of Y. Zhang and W. Yang,
3442
Phys. Rev. Lett. \textbf{80}, 890 (1998)
3443
\item \texttt{RPBE}: \index{RPBE}
3444
Modified GGA-PBE functional of
3445
B. Hammer, L. B. Hansen and J. K. Norskov Phys. Rev. B \textbf{59}, 7413 (1999)
3446
\item \texttt{WC} \index{WC}
3447
Modified GGA-PBE functional of
3448
Z. Wu and R. E. Cohen, Phys. Rev. B \textbf{73}, 235116 (2006)
3449
\item \texttt{AM05}: \index{AM05}
3450
Modified GGA-PBE functional of
3451
R. Armiento and A. E. Mattsson, Phys. Rev. B \textbf{72}, 085108 (2005)
3452
\item \texttt{PBEsol}: \index{PBEsol}
3453
Modified GGA-PBE functional of
3454
J. P. Perdew et al, Phys. Rev. Lett. \textbf{100}, 136406 (2008)
3455
\item \texttt{PBEJsJrLO}:
3456
GGA-PBE functional with parameters $\beta, \mu$, and $\kappa$ fixed by
3457
the jellium surface (Js), jellium response (Jr),
3458
and Lieb-Oxford bound (LO) criteria, respectively, as described by
3459
L. S. Pedroza, A. J. R. da Silva, and K. Capelle,
3460
Phys. Rev. B \textbf{79}, 201106(R) (2009), and by
3461
M. M. Odashima, K. Capelle, and S. B. Trickey,
3462
J. Chem. Theory Comput. \textbf{5}, 798 (2009)
3463
\item \texttt{PBEJsJrHEG}:
3464
Same as PBEJsJrLO, with parameter $\kappa$ fixed by the Lieb-Oxford bound
3465
for the low density limit of the homogeneous electron gas (HEG)
3466
\item \texttt{PBEGcGxLO}:
3467
Same as PBEJsJrLO, with parameters $\beta$ and $\mu$ fixed by the
3468
gradient expansion of correlation (Gc) and exchange (Gx), respectively
3469
\item \texttt{PBEGcGxHEG}:
3470
Same as previous ones, with parameters $\beta,\mu$, and $\kappa$ fixed by
3471
the Gc, Gx, and HEG criteria, respectively.
3472
\item \texttt{BLYP} (equivalent to \texttt{LYP}): \index{BLYP}
3473
GGA with Becke exchange (A. D. Becke, Phys. Rev. A \textbf{38}, 3098 (1988))
3474
and Lee-Yang-Parr correlation
3475
(C. Lee, W. Yang, R. G. Parr, Phys. Rev. B \textbf{37}, 785 (1988)),
3476
as modified by B. Miehlich, A. Savin, H. Stoll, and H. Preuss,
3477
Chem. Phys. Lett. \textbf{157}, 200 (1989).
3478
See also B. G. Johnson, P. M. W. Gill and J. A. Pople,
3479
J. Chem. Phys. \textbf{98}, 5612 (1993). (Some errors were detected in this
3480
last paper, so not all of their expressions correspond exactly to those
3481
implemented in \siesta)
3482
\item \texttt{DRSLL} (equivalent to \texttt{DF1}): \index{vdW-DF1} \index{DRSLL}
3483
van der Waals \index{vdW} density functional (vdW-DF) \index{vdW-DF}
3484
of M. Dion, H. Rydberg, E. Schr\"{o}der, D. C. Langreth, and B. I. Lundqvist,
3485
Phys. Rev. Lett. \textbf{92}, 246401 (2004), with the efficient implementation of
3486
G. Rom\'an-P\'erez and J. M. Soler, Phys. Rev. Lett. \textbf{103}, 096102 (2009)
3487
\item \texttt{LMKLL} (equivalent to \texttt{DF2}): \index{vdW-DF2} \index{LMKLL}
3488
vdW-DF functional of Dion \textit{et al} (same as DRSLL)
3489
reparametrized by K. Lee, E. Murray, L. Kong, B. I. Lundqvist and
3490
D. C. Langreth, Phys. Rev. B \textbf{82}, 081101 (2010)
3492
vdW-DF functional of Dion \textit{et al} (same as DRSLL)
3493
with exchange modified by J. Klimes, D. R. Bowler, and A. Michaelides,
3494
J. Phys.: Condens. Matter \textbf{22}, 022201 (2010) (optB88-vdW version)
3496
vdW-DF functional of Dion \textit{et al} (same as DRSLL)
3497
with exchange modified by V. R. Cooper, Phys. Rev. B \textbf{81}, 161104 (2010)
3499
vdW-DF functional of Dion \textit{et al} (same as DRSLL)
3500
with exchange modified by
3501
K. Berland and P. Hyldgaard, Phys. Rev. B 89, 035412 (2014)
3503
vdW-DF functional of O. A. Vydrov and T. Van Voorhis,
3504
J. Chem. Phys. \textbf{133}, 244103 (2010)
3508
\textit{Use:} \textbf{XC.functional} and \textbf{XC.authors} must be compatible.
3510
\textit{Default value:} \texttt{PZ}
3512
\item[\textbf{XC.hybrid}] (\textit{data block}):
3513
\index{XC.hybrid@\textbf{XC.hybrid}}
3514
This data block allows the user to create a ``cocktail'' functional by
3515
mixing the desired amounts of exchange and correlation from each of
3516
the functionals described under XC.authors. Note that these ``mixed''
3517
functionals do \emph{not} have the exact Hartree-Fock exchange which
3518
is a key ingredient of the true ``hybrid'' functionals. The use of
3519
the word ``hybrid'' in the label is unfortunate in this regard, and
3520
might be deprecated in a future version.
3522
The first line of the block must contain the number of functionals to
3523
be mixed. On the subsequent lines the values of XC.functl and
3524
XC.authors must be given and then the weights for the exchange and
3525
correlation, in that order. If only one number is given then the same
3526
weight is applied to both exchange and correlation.
3528
The following is an example in which a 75:25 mixture of Ceperley-Alder
3529
and PBE correlation is made, with an equal split of the exchange
3372
\begin{fdfentry}{XC!Functional}[string]<LDA>
3374
Exchange-correlation functional type. May be \fdf*{LDA} (local
3375
density approximation, equivalent to \fdf*{LSD}), \fdf*{GGA}
3376
(Generalized Gradient Approximation), or \fdf*{VDW} (van der
3381
\begin{fdfentry}{XC!Authors}[string]<PZ>
3382
\newcommand\xcidx[1]{\index{exchange-correlation!#1}}
3384
Particular parametrization of the exchange-correlation
3385
functional. Options are:
3388
\fdf*{CA} (equivalent to \fdf*{PZ}): \xcidx{CA} \xcidx{PZ}
3389
(Spin) local density approximation (LDA/LSD). \xcidx{LDA}
3390
\xcidx{LSD} Quantum Monte Carlo calculation of the homogeneous
3391
electron gas by D. M. Ceperley and B. J. Alder,
3392
Phys. Rev. Lett. \textbf{45},566 (1980), as parametrized by
3393
J. P. Perdew and A. Zunger, Phys. Rev B \textbf{23}, 5075 (1981)
3396
\fdf*{PW92}: \xcidx{PW92}
3397
LDA/LSD, as parametrized by
3398
J. P. Perdew and Y. Wang, Phys. Rev B, \textbf{45}, 13244 (1992)
3401
\fdf*{PW91}: \xcidx{PW91}%
3402
Generalized gradients approximation (GGA) \xcidx{GGA} of Perdew and Wang.
3403
Ref: P\&W, J. Chem. Phys., \textbf{100}, 1290 (1994)
3406
\fdf*{PBE}: \xcidx{PBE}%
3407
GGA of J. P. Perdew, K. Burke and M. Ernzerhof,
3408
Phys. Rev. Lett. \textbf{77}, 3865 (1996)
3411
\fdf*{revPBE}: \xcidx{revPBE}%
3412
Modified GGA-PBE functional of Y. Zhang and W. Yang,
3413
Phys. Rev. Lett. \textbf{80}, 890 (1998)
3416
\fdf*{RPBE}: \xcidx{RPBE}%
3417
Modified GGA-PBE functional of
3418
B. Hammer, L. B. Hansen and J. K. Norskov Phys. Rev. B \textbf{59}, 7413 (1999)
3421
\fdf*{WC}: \xcidx{WC}%
3422
Modified GGA-PBE functional of
3423
Z. Wu and R. E. Cohen, Phys. Rev. B \textbf{73}, 235116 (2006)
3426
\fdf*{AM05}: \xcidx{AM05}%
3427
Modified GGA-PBE functional of
3428
R. Armiento and A. E. Mattsson, Phys. Rev. B \textbf{72}, 085108 (2005)
3431
\fdf*{PBEsol}: \xcidx{PBEsol}%
3432
Modified GGA-PBE functional of
3433
J. P. Perdew et al, Phys. Rev. Lett. \textbf{100}, 136406 (2008)
3436
\fdf*{PBEJsJrLO}: \xcidx{PBEJsJrLO}%
3437
GGA-PBE functional with parameters $\beta, \mu$, and $\kappa$ fixed by
3438
the jellium surface (Js), jellium response (Jr),
3439
and Lieb-Oxford bound (LO) criteria, respectively, as described by
3440
L. S. Pedroza, A. J. R. da Silva, and K. Capelle,
3441
Phys. Rev. B \textbf{79}, 201106(R) (2009), and by
3442
M. M. Odashima, K. Capelle, and S. B. Trickey,
3443
J. Chem. Theory Comput. \textbf{5}, 798 (2009)
3446
\fdf*{PBEJsJrHEG}: \xcidx{PBEJsJrHEG}%
3447
Same as PBEJsJrLO, with parameter $\kappa$ fixed by the Lieb-Oxford bound
3448
for the low density limit of the homogeneous electron gas (HEG)
3451
\fdf*{PBEGcGxLO}: \xcidx{PBEGcGxLO}%
3452
Same as PBEJsJrLO, with parameters $\beta$ and $\mu$ fixed by the
3453
gradient expansion of correlation (Gc) and exchange (Gx), respectively
3456
\fdf*{PBEGcGxHEG}: \xcidx{PBEGcGxHEG}%
3457
Same as previous ones, with parameters $\beta,\mu$, and $\kappa$ fixed by
3458
the Gc, Gx, and HEG criteria, respectively.
3461
\fdf*{BLYP} (equivalent to \fdf*{LYP}): \xcidx{BLYP}%
3462
GGA with Becke exchange (A. D. Becke, Phys. Rev. A \textbf{38}, 3098 (1988))
3463
and Lee-Yang-Parr correlation
3464
(C. Lee, W. Yang, R. G. Parr, Phys. Rev. B \textbf{37}, 785 (1988)),
3465
as modified by B. Miehlich, A. Savin, H. Stoll, and H. Preuss,
3466
Chem. Phys. Lett. \textbf{157}, 200 (1989).
3467
See also B. G. Johnson, P. M. W. Gill and J. A. Pople,
3468
J. Chem. Phys. \textbf{98}, 5612 (1993). (Some errors were detected in this
3469
last paper, so not all of their expressions correspond exactly to those
3470
implemented in \siesta)
3473
\fdf*{DRSLL} (equivalent to \fdf*{DF1}): \xcidx{vdW-DF1}
3475
van der Waals \xcidx{vdW} density functional (vdW-DF) \xcidx{vdW-DF}
3476
of M. Dion, H. Rydberg, E. Schr\"{o}der, D. C. Langreth, and B. I. Lundqvist,
3477
Phys. Rev. Lett. \textbf{92}, 246401 (2004), with the efficient implementation of
3478
G. Rom\'an-P\'erez and J. M. Soler, Phys. Rev. Lett. \textbf{103}, 096102 (2009)
3481
\fdf*{LMKLL} (equivalent to \fdf*{DF2}): \xcidx{vdW-DF2} \xcidx{LMKLL}%
3482
vdW-DF functional of Dion \textit{et al} (same as DRSLL)
3483
reparametrized by K. Lee, E. Murray, L. Kong, B. I. Lundqvist and
3484
D. C. Langreth, Phys. Rev. B \textbf{82}, 081101 (2010)
3487
\fdf*{KBM}: \xcidx{KBM}%
3488
vdW-DF functional of Dion \textit{et al} (same as DRSLL)
3489
with exchange modified by J. Klimes, D. R. Bowler, and A. Michaelides,
3490
J. Phys.: Condens. Matter \textbf{22}, 022201 (2010) (optB88-vdW version)
3493
\fdf*{C09}: \xcidx{C09}%
3494
vdW-DF functional of Dion \textit{et al} (same as DRSLL)
3495
with exchange modified by V. R. Cooper, Phys. Rev. B \textbf{81}, 161104 (2010)
3498
\fdf*{BH}: \xcidx{BH}%
3499
vdW-DF functional of Dion \textit{et al} (same as DRSLL)
3500
with exchange modified by
3501
K. Berland and P. Hyldgaard, Phys. Rev. B 89, 035412 (2014)
3504
\fdf*{VV}: \xcidx{VV}%
3505
vdW-DF functional of O. A. Vydrov and T. Van Voorhis,
3506
J. Chem. Phys. \textbf{133}, 244103 (2010)
3513
\begin{fdfentry}{XC!Hybrid}[block]
3515
This data block allows the user to create a ``cocktail'' functional by
3516
mixing the desired amounts of exchange and correlation from each of
3517
the functionals described under XC.authors. Note that these ``mixed''
3518
functionals do \emph{not} have the exact Hartree-Fock exchange which
3519
is a key ingredient of the true ``hybrid'' functionals. The use of
3520
the word ``hybrid'' in the label is unfortunate in this regard, and
3521
might be deprecated in a future version.
3523
The first line of the block must contain the number of functionals to
3524
be mixed. On the subsequent lines the values of XC.functl and
3525
XC.authors must be given and then the weights for the exchange and
3526
correlation, in that order. If only one number is given then the same
3527
weight is applied to both exchange and correlation.
3529
The following is an example in which a 75:25 mixture of Ceperley-Alder
3530
and PBE correlation is made, with an equal split of the exchange
3533
3534
%block XC.hybrid
3535
3536
LDA CA 0.5 0.75
3536
3537
GGA PBE 0.5 0.25
3537
3538
%endblock XC.hybrid
3540
\textit{Default value:} \texttt{not hybrid}
3546
3545
\subsection{Spin polarization}
3552
\item[\textbf{Spin}] (\textit{string}):
3553
\index{Spin@\textbf{Spin}}\index{spin}
3554
Choose spin configuration. The following are available:
3559
\item[non-polarized] Spin unpolarized
3561
\item[polarized] Spin polarized
3563
\item[non-collinear] Non collinear spin
3565
Refs: T. Oda et al, PRL, \textbf{80}, 3622 (1998);
3566
V. M. Garc\'{\i}a-Su\'arez et al, Eur. Phys. Jour. B \textbf{40}, 371 (2004);
3567
V. M. Garc\'{\i}a-Su\'arez et al, Journal of
3568
Phys: Cond. Matt \textbf{16}, 5453 (2004).
3570
\item[spin-orbit] Spin-orbit coupling
3572
This flag requires that the generated pseudopotential be relativistic.
3574
\emph{Currently this is in experimental use.}
3577
This flag has precedence from \textbf{SpinOrbit}, \textbf{NonCollinearSpin} and
3578
\textbf{SpinPolarized} while they may still be used.
3580
Certain options may not be used together with specific parallelization routines.
3581
For instance only a spin-polarized calculation may use the \fdf{Diag!ParallelOverK}
3584
\textit{Default value:} \texttt{non-polarized}
3586
\item[\textbf{FixSpin}] (\textit{logical}):
3587
\index{FixSpin@\textbf{FixSpin}}\index{spin}
3588
\index{fixed spin state}\index{LSD}
3589
If \texttt{.true.}, the calculation is done with a fixed value of the
3590
spin of the system, defined by variable \textbf{TotalSpin}.
3591
This option can only be used for collinear spin polarized
3594
\textit{Default value:} \texttt{.false.}
3596
\item[\textbf{TotalSpin}] (\textit{real}):
3597
\index{TotalSpin@\textbf{TotalSpin}}
3599
\index{fixed spin state}\index{LSD}
3600
Value of the imposed total spin polarization of the system (in units of the
3601
electron spin, 1/2). It is only used
3602
if \textbf{FixSpin} = \texttt{.true.}
3604
\textit{Default value:} 0.0
3606
\item[\textbf{SingleExcitation}] (\textit{logical}):
3607
\index{SingleExcitation@\textbf{SingleExcitation}}
3608
If true, \siesta\ calculates a very rough approximation to
3609
the lowest excited state by swapping the populations of the HOMO
3610
and the LUMO. If there is no spin polarisation, it is half swap only.
3611
It is done for the first spin component (up) and first k vector.
3613
\textit{Default value:} \texttt{.false.}
3548
\begin{fdfentry}{Spin}[string]<non-polarized>
3550
Choose the spin-components in the simulation.
3552
\note this flag has precedence from \fdf*{SpinOrbit}, \fdf*{NonCollinearSpin} and
3553
\fdf*{SpinPolarized} while they may still be used.
3556
\option[non-polarized]%
3557
Perform a calculation with spin-degeneracy (only one component).
3560
Perform a calculation with collinear spin (two spin components).
3562
\option[non-collinear]%
3563
Perform a calculation with non-collinear spin (4 spin components),
3566
Refs: T. Oda et al, PRL, \textbf{80}, 3622 (1998);
3567
V. M. Garc\'{\i}a-Su\'arez et al, Eur. Phys. Jour. B \textbf{40}, 371 (2004);
3568
V. M. Garc\'{\i}a-Su\'arez et al, Journal of
3569
Phys: Cond. Matt \textbf{16}, 5453 (2004).
3571
\option[spin-orbit]%
3572
Perform a calculation with spin-orbit coupling. This requires the
3573
pseudopotentials to be relativistic.
3575
See Sect.~\ref{sec:spin-orbit}.
3579
Certain options may not be used together with specific
3580
parallelization routines. For instance only a spin-polarized
3581
calculation may use the \fdf{Diag!ParallelOverK} option.
3585
\begin{fdflogicalF}{Spin!Fix}
3586
\index{fixed spin state}\index{LSD}
3587
\fdfindex{FixSpin}[|Spin!Fix]
3589
If \fdftrue, the calculation is done with a fixed value of the spin
3590
of the system, defined by variable \fdf{Spin!Total}. This option can
3591
only be used for collinear spin polarized calculations.
3595
\begin{fdfentry}{Spin!Total}[real]<$0$>
3596
\index{fixed spin state}\index{LSD}
3598
\fdfindex{TotalSpin}[|Spin!Total]
3600
Value of the imposed total spin polarization of the system (in units
3601
of the electron spin, 1/2). It is only used if \fdf{Spin!Fix} \fdftrue.
3605
\begin{fdflogicalF}{SingleExcitation}
3607
If \fdftrue, \siesta\ calculates a very rough approximation to the
3608
lowest excited state by swapping the populations of the HOMO and the
3609
LUMO. If there is no spin polarisation, it is half swap only. It is
3610
done for the first spin component (up) and first $k$ vector.
3618
3615
\subsection{Spin--Orbit coupling}
3616
\label{sec:spin-orbit}
3620
3618
\siesta\ includes the posibility to perform fully relativistic
3621
3619
calculations by means of the inclusion in the total Hamiltonian not
3747
3746
step; if mixing H, dHmax refers to H(out)-H(in) in the previous(?)
3750
\item When achieving convergence, the loop might be exited without a further
3751
mixing of the DM, thus preserving DM(out) for further processing
3752
(including the calculation of forces and the analysis of the
3753
electronic structure) (see the \textbf{MixAfterConvergence} option).
3755
\index{Variational character of E\_KS}
3757
\item It remains to be seen whether the forces, being computed
3758
``right'' on the basis of DM(out), exhibit somehow better convergence
3759
as a function of the scf step. In order to gain some more data and
3760
heuristics on this we have implemented a force-monitoring option,
3761
activated by setting to \texttt{.true.} the variable \fdf{SCF.MonitorForces}. The program will then print the maximum
3762
absolute value of the change in forces from one step to the
3749
\item When achieving convergence, the loop might be exited without a
3750
further mixing of the DM, thus preserving DM(out) for further
3751
processing (including the calculation of forces and the analysis of
3752
the electronic structure) (see the \fdf{SCF.MixAfterConvergence}
3755
\index{Variational character of E\_KS}
3757
\item It remains to be seen whether the forces, being computed
3758
``right'' on the basis of DM(out), exhibit somehow better
3759
convergence as a function of the scf step. In order to gain some
3760
more data and heuristics on this we have implemented a
3761
force-monitoring option, activated by setting to \fdftrue\ the
3762
variable \fdf{SCF.MonitorForces}. The program will then print the
3763
maximum absolute value of the change in forces from one step to the
3763
3764
next. Other statistics could be implemented.
3765
\item While the (mixed) DM is saved at every SCF step, as was standard
3766
practice, the final DM(out) overwrites the .DM file at the end of
3767
the SCF cycle. Thus it is still possible to use a ``mixed'' DM for
3768
restarting an interrupted loop, but a ``good'' DM will be used for
3769
any other post-processing.
3766
\item While the (mixed) DM is saved at every SCF step, as was
3767
standard practice, the final DM(out) overwrites the .DM file at the
3768
end of the SCF cycle. Thus it is still possible to use a ``mixed''
3769
DM for restarting an interrupted loop, but a ``good'' DM will be
3770
used for any other post-processing.
3774
3775
\subsubsection{Harris functional and basic options}
3780
\item[\textbf{Harris\_functional}] (\textit{logical}):
3781
\index{Harris\_functional@\textbf{Harris\_functional}}
3782
Logical variable to choose between self-consistent Kohn-Sham functional or
3783
non self-consistent Harris functional to calculate energies and forces.
3785
\item \texttt{.false.} : Fully self-consistent Kohn-Sham functional.
3786
\item \texttt{.true.} : Non self consistent Harris functional. Cheap but
3787
pretty crude for some systems. The forces are computed within the
3788
Harris functional in the first SCF step. Only implemented for LDA in
3789
the Perdew-Zunger parametrization. It really only applies to starting
3790
densities which are superpositions of atomic charge densities.
3792
When this option is choosen, the values of DM.UseSaveDM,
3793
MaxSCFIterations, SCFMustConverge and \fdf{DM.MixSCF1} are automatically set
3794
up to False, 1, False and False respectively, no matter whatever other
3795
specification are in the INPUT file.
3798
\textit{Default value:} \texttt{.false.}
3800
\item[\textbf{MinSCFIterations}] (\textit{integer}):
3801
\index{MinSCFIterations@\textbf{MinSCFIterations}}
3802
Minimum number of SCF\index{SCF} iterations per time step. In MD simulations
3803
this can with benefit be set to 3.
3805
\textit{Default value:} \texttt{0}
3807
\item[\textbf{MaxSCFIterations}] (\textit{integer}):
3808
\index{MaxSCFIterations@\textbf{MaxSCFIterations}}
3809
Maximum number of SCF\index{SCF} iterations per time step.
3811
\textit{Default value:} \texttt{50}
3813
\item[\textbf{SCFMustConverge}] (\textit{logical}):
3814
\index{SCFMustConverge@\textbf{SCFMustConverge}}
3815
Defines the behaviour if convergence is not reached in the maximum
3816
number of SCF iterations. The default is to update the forces, perform
3817
an MD or geometry optimisation step and carry on. When set to true
3818
the calculation will stop on the first SCF convergence failure.
3820
\textit{Default value:} \texttt{.false.}
3777
\begin{fdflogicalF}{Harris!Functional}
3779
Logical variable to choose between self-consistent Kohn-Sham
3780
functional or non self-consistent Harris functional to calculate
3781
energies and forces.
3783
\item \fdffalse: Fully self-consistent Kohn-Sham functional.
3784
\item \fdftrue: Non self consistent Harris functional. Cheap but
3785
pretty crude for some systems. The forces are computed within the
3786
Harris functional in the first SCF step. Only implemented for LDA in
3787
the Perdew-Zunger parametrization. It really only applies to starting
3788
densities which are superpositions of atomic charge densities.
3790
When this option is choosen, the values of \fdf{DM!UseSaveDM},
3791
\fdf{MaxSCFIterations}, \fdf{SCF!MustConverge} and \fdf{SCF.Mix!First} are automatically set
3792
up to False, 1, False and False respectively, no matter whatever other
3793
specification are in the INPUT file.
3798
\begin{fdfentry}{MinSCFIterations}[integer]<0>
3800
Minimum number of SCF\index{SCF} iterations per time step. In MD
3801
simulations this can with benefit be set to 3.
3805
\begin{fdfentry}{MaxSCFIterations}[integer]<50>
3807
Maximum number of SCF\index{SCF} iterations per time step.
3811
\begin{fdflogicalF}{SCF!MustConverge}
3813
Defines the behaviour if convergence is not reached in the maximum
3814
number of SCF iterations. The default is to update the forces, perform
3815
an MD or geometry optimisation step and carry on. When set to true
3816
the calculation will stop on the first SCF convergence failure.
3824
3821
\subsubsection{Mixing options}
3825
3822
\index{SCF!mixing}
3827
Reaching self-consistency is heavily relying on the used mixing
3829
Any user, will be adviced to test mixing options for their
3830
calculations as choosing optimal parameters can mean a difference of
3831
several orders of consistency cycles.
3824
Whether a calculation reaches self-consistency in a moderate number of
3825
steps depends strongly on the mixing parameters used. The available
3826
mixing options should be carefully tested for a given calculation
3827
type. This search for optimal parameters can repay itself handsomely
3828
by potentially saving many self-consistency steps in production runs.
3834
3831
\begin{fdfentry}{SCF.Mix}[string]<Hamiltonian|density|charge>
4835
\item[\textbf{DM.UseSaveDM}] (\textit{logical}):
4836
\index{DM.UseSaveDM@\textbf{DM.UseSaveDM}}
4837
\index{reading saved data!density matrix}
4838
Instructs to read the density matrix stored in file
4839
\texttt{SystemLabel}.DM by a previous run.
4841
\textit{Use:} If the required file does not exist, a warning is
4842
printed but the program does not stop. Overrides \textbf{UseSaveData}.
4844
\textit{Default value:} \texttt{.false.}
4847
\item[\textbf{DM.FormattedFiles}] (\textit{logical}):
4848
\index{DM.FormattedFiles@\textbf{DM.FormattedFiles}}
4849
\index{reading saved data!density matrix}
4850
Instructs to use formatted files for reading and writing
4851
the density matrix. In this case, the files are labelled
4852
\texttt{SystemLabel}.DMF.
4854
\textit{Use:} This makes for much larger files, and slower i/o. However, the
4855
files are transferable between different computers, which is
4856
not the case normally.
4858
\textit{Default value:} \texttt{.false.}
4861
\item[\textbf{DM.FormattedInput}] (\textit{logical}):
4862
\index{DM.FormattedInput@\textbf{DM.FormattedInput}}
4863
\index{reading saved data!density matrix}
4864
Instructs to use formatted files for reading the density
4867
\textit{Use:} Overrides the value of \textbf{DM.FormattedFiles}.
4869
\textit{Default value:} \texttt{.false.}
4872
\item[\textbf{DM.FormattedOutput}] (\textit{logical}):
4873
\index{DM.FormattedInput@\textbf{DM.FormattedOutput}}
4874
\index{reading saved data!density matrix}
4875
Instructs to use formatted files for writing the density
4878
\textit{Use:} Overrides the value of \textbf{DM.FormattedFiles}.
4880
\textit{Default value:} \texttt{.false.}
4883
\item[\textbf{DM.InitSpinAF}] (\textit{logical}):\index{spin}
4884
\index{spin!initialization}\index{ferromagnetic initial DM}
4885
\index{antiferromagnetic initial DM}
4886
\index{DM.InitSpinAF@\textbf{DM.InitSpinAF}}
4887
It defines the initial spin density for a spin polarized calculation.
4888
The spin density is initially constructed with the maximum possible
4889
spin polarization for each atom in its atomic configuration.
4890
This variable defines the relative orientation of the atomic
4894
\item \texttt{.false.} gives ferromagnetic order (all spins up).
4895
\item \texttt{.true.} gives antiferromagnetic order. Up and down are
4896
assigned according to order in the block
4897
\textbf{AtomicCoordinatesAndAtomicSpecies}: up for the odd atoms, down for even.
4900
\textit{Default value:} \texttt{.false.}
4903
\item[\textbf{DM.InitSpin}] (\textit{data block}):
4904
\index{DM.InitSpin@\textbf{DM.InitSpin}} It defines the
4905
initial spin density for a spin polarized calculation atom by atom.
4906
In the block there is one line per atom to be spin-polarized,
4907
containing the atom index (integer, ordinal in the block
4908
\textbf{AtomicCoordinatesAndAtomicSpecies}) and the desired
4909
initial spin-polarization (real, positive for spin up, negative for
4910
spin down). A value larger than possible will be reduced
4911
to the maximum possible polarization, keeping its sign.
4912
Maximum polarization can also be given by introducing the
4913
symbol \texttt{+} or \texttt{-} instead of the polarization value.
4914
There is no need to include a line for every atom, only for
4915
those to be polarized. The atoms not contemplated in the block will
4916
be given non-polarized initialization.
4917
For non-collinear spin, the spin direction may be specified for
4918
each atom by the polar angles theta and phi, given as the last
4919
two arguments in degrees. If not specified, theta=0 is assumed.
4920
\textbf{NonCollinearSpin} must be \texttt{.true.} to use the spin direction.
4826
\begin{fdflogicalF}{DM!UseSaveDM}
4827
\index{reading saved data!density matrix}
4829
Instructs to read the density matrix stored in file
4830
\sysfile{.DM} by a previous run.
4832
\siesta\ will continue even if \sysfile*{DM} is not found.
4836
\begin{fdflogicalF}{DM!FormattedFiles}
4837
\index{reading saved data!density matrix}
4839
Setting this alters the default for \fdf{DM!FormattedInput} and
4840
\fdf{DM!FormattedOutput}.
4841
Instructs to use formatted files for reading and writing
4842
the density matrix. In this case, the files are labelled
4845
Only usable if one has problems transferring files from one computer
4850
\begin{fdflogicalF}{DM!FormattedInput}
4851
\index{reading saved data!density matrix}
4853
Instructs to use formatted files for reading the density
4858
\begin{fdflogicalF}{DM!FormattedOutput}
4859
\index{reading saved data!density matrix}
4861
Instructs to use formatted files for writing the density
4866
\begin{fdflogicalF}{DM!InitSpinAF}
4867
\index{spin!initialization}%
4868
\index{ferromagnetic initial DM}%
4869
\index{antiferromagnetic initial DM}
4871
It defines the initial spin density for a spin polarized calculation.
4872
The spin density is initially constructed with the maximum possible
4873
spin polarization for each atom in its atomic configuration.
4874
This variable defines the relative orientation of the atomic
4877
If \fdffalse\ the initial spin-configuration is a ferromagnetic
4878
order (all spins up).
4879
If \fdftrue\ all odd atoms are initialized to spin-up, all even
4880
atoms are initialized to spin-down.
4884
\begin{fdfentry}{DM!InitSpin}[block]
4886
It defines the initial spin density for a spin polarized calculation
4887
atom by atom. In the block there is one line per atom to be
4888
spin-polarized, containing the atom index (integer, ordinal in the
4889
block \fdf{AtomicCoordinatesAndAtomicSpecies}) and the desired
4890
initial spin-polarization (real, positive for spin up, negative for
4891
spin down). A value larger than possible will be reduced to the
4892
maximum possible polarization, keeping its sign. Maximum
4893
polarization can also be given by introducing the symbol \texttt{+}
4894
or \texttt{-} instead of the polarization value. There is no need
4895
to include a line for every atom, only for those to be
4896
polarized. The atoms not contemplated in the block will be given
4897
non-polarized initialization. For non-collinear spin, the spin
4898
direction may be specified for each atom by the polar angles theta
4899
and phi, given as the last two arguments in degrees. If not
4900
specified, theta=0 is assumed. \fdf{Spin} \fdf*{non-collinear} must be
4901
\fdftrue\ to use the spin direction.
4925
4905
%block DM.InitSpin
4926
4906
5 -1. 90. 0. # Atom index, spin, theta, phi (deg)
4929
4909
%endblock DM.InitSpin
4932
\textit{Default value:} If present but empty, all atoms are not polarized.
4933
If absent, \textbf{DM.InitSpinAF} defines the polarization.
4935
\item[\textbf{DM.AllowReuse}] (\textit{logical}):
4936
\index{DM.AllowReuse@\textbf{DM.AllowReuse}}
4937
Controls whether density matrix information from previous geometry
4938
iterations is re-used to start the new geometry's SCF cycle.
4940
\textit{Default value:} \texttt{.true.}
4943
\item[\textbf{DM.AllowExtrapolation}] (\textit{logical}):
4944
\index{DM.AllowExtrapolation@\textbf{DM.AllowExtrapolation}} Controls
4945
whether the density matrix information from several previous
4946
geometry iterations is extrapolated to start the new geometry's SCF
4947
cycle. This feature is useful for molecular dynamics simulations
4948
and possibly also for geometry relaxations. The number of geometry
4949
steps saved is controlled by the variable \textbf{DM.History.Depth}
4951
\textit{Default values:} \texttt{.true.} for molecular-dynamics simulations,
4952
but \texttt{.false.} for now for geometry-relaxations (pending further
4953
tests which users are kindly requested to perform).
4956
\item[\textbf{DM.History.Depth}] (\textit{integer}):
4957
\index{DM.History.Depth@\textbf{DM.Hitory.Depth}}
4958
Sets the number of geometry steps for which density-matrix information
4959
is saved for extrapolation.
4961
\textit{Default value:} 4
4912
If this block is defined, but empty, all atoms are not polarized.
4913
This block has precedence over \fdf{DM!InitSpinAF}.
4917
\begin{fdflogicalT}{DM!AllowReuse}
4919
Controls whether density matrix information from previous geometry
4920
iterations is re-used to start the new geometry's SCF cycle.
4924
\begin{fdflogicalT}{DM!AllowExtrapolation}
4926
Controls whether the density matrix information from several
4927
previous geometry iterations is extrapolated to start the new
4928
geometry's SCF cycle. This feature is useful for molecular dynamics
4929
simulations and possibly also for geometry relaxations. The number
4930
of geometry steps saved is controlled by the variable
4931
\fdf{DM!History.Depth}.
4933
This is default \fdftrue\ for molecular-dynamics simulations, but
4934
\fdffalse, for now, for geometry-relaxations (pending further tests
4935
which users are kindly requested to perform).
4940
\begin{fdfentry}{DM!History.Depth}[integer]<1>
4942
Sets the number of geometry steps for which density-matrix information
4943
is saved for extrapolation.
4965
4949
\subsubsection{Initialization of the SCF cycle with charge densities}
4968
\item[\textbf{SCF.Read.Charge.NetCDF}] (\textit{logical}):
4969
\index{SCF.Read.Charge.NetCDF@\textbf{SCF.Read.Charge.NetCDF}}
4970
\index{reading saved data!charge density} Instructs Siesta to read
4971
the charge density stored in the netCDF file
4972
\file{Rho.IN.grid.nc}. This feature allows the easier re-use of
4952
\begin{fdflogicalF}{SCF.Read.Charge.NetCDF}
4953
\index{reading saved data!charge density}
4955
Instructs \siesta\ to read the charge density stored in the netCDF
4956
file \file{Rho.IN.grid.nc}. This feature allows the easier re-use of
4973
4957
electronic-structure information from a previous run. It is not
4974
4958
necessary that the basis sets are ``similar'' (a requirement if
4975
4959
density-matrices are to be read in).
4977
\textit{Use:} This is an experimental feature. Until robust checks are
4978
implemented, care must be taken to make sure that the FFT grids in the
4979
.grid.nc file and in Siesta are the same.
4981
\textit{Default value:} \texttt{.false.}
4983
\item[\textbf{SCF.Read.Deformation.Charge.NetCDF}] (\textit{logical}):
4984
\index{SCF.Read.Deformation.Charge.NetCDF@\textbf{SCF.Read.Deformation.Charge.NetCDF}} \index{reading saved
4985
data!deformation charge density} Instructs Siesta to read the
4986
deformation charge density stored in the netCDF file
4987
\file{DeltaRho.IN.grid.nc}. This feature allows the easier re-use of
4988
electronic-structure information from a previous run. It is not
4989
necessary that the basis sets are ``similar'' (a requirement if
4990
density-matrices are to be read in). The deformation charge is
4991
particularly useful to give a good starting point for slightly
4992
different geometries.
4994
\textit{Use:} This is an experimental feature. Until robust checks are
4995
implemented, care must be taken to make sure that the FFT grids in the
4996
.grid.nc file and in Siesta are the same.
4998
\textit{Default value:} \texttt{.false.}
4961
\note this is an experimental feature. Until robust checks are
4962
implemented, care must be taken to make sure that the FFT grids in
4963
the \sysfile*{grid.nc} file and in \siesta\ are the same.
4968
\begin{fdflogicalF}{SCF.Read.Deformation.Charge.NetCDF}
4969
\index{reading saved data!deformation charge density}
4971
Instructs Siesta to read the deformation charge density stored in
4972
the netCDF file \file{DeltaRho.IN.grid.nc}. This feature allows the
4973
easier re-use of electronic-structure information from a previous
4974
run. It is not necessary that the basis sets are ``similar'' (a
4975
requirement if density-matrices are to be read in). The deformation
4976
charge is particularly useful to give a good starting point for
4977
slightly different geometries.
4979
\note this is an experimental feature. Until robust checks are
4980
implemented, care must be taken to make sure that the FFT grids in
4981
the \sysfile*{grid.nc} file and in Siesta are the same.
5003
4986
\subsubsection{Output of density matrix and Hamiltonian}
5409
5401
%endblock GridCellSampling
5412
gives again a cubic sampling with half the original side length.
5413
It is not trivial to choose a right set of displacements so as
5414
to maximize the new 'effective' cutoff. It depends on the
5415
kind of cell. It may be automatized in the future, but it
5416
is now left to the user, who introduces the displacements
5417
manually through this block.
5419
The quantities which are symmetrized are: ($i$) energy terms
5420
that depend on the grid, ($ii$) forces, ($iii$) stress
5421
tensor, and ($iv$) electric dipole.
5423
The symmetrization is performed at the end of every SCF cycle. The
5424
whole cycle is done for the (0,0,0) displacement, and, when the
5425
density matrix is converged, the same (now fixed)
5426
density matrix is used to obtain the desired quantities at the
5427
other displacements (the density matrix itself is \textit{not}
5428
symmetrized as it gives a much smaller egg-box effect).
5429
The CPU time needed for each displacement
5430
in the \textbf{GridCellSampling} block
5431
is of the order of one extra SCF iteration.
5433
\textit{Default value:} Empty.
5436
\item[\textbf{EggboxRemove}] (\textit{data block}):\index{egg-box effect}
5437
\index{EggboxRemove@\textbf{EggboxRemove}}\index{rippling}
5439
For recovering translational invariance in an approximate way.
5441
It works by substracting from Kohn-Sham's total energy (and forces) an
5442
approximation to the eggbox energy, sum of atomic contributions.
5443
Each atom has a predefined eggbox energy depending on where it sits on
5444
the cell. This atomic contribution is species dependent and is
5445
obviously invariant under grid-cell translations. Each species
5446
contribution is thus expanded in the appropriate Fourier series.
5447
It is important to have a smooth eggbox, for it to
5448
be represented by a few Fourier components. A jagged egg-box (unless
5449
very small, which is then unimportant) is often an
5450
indication of a problem with the pseudo.
5452
In the block there is one line per Fourier component. The first integer
5453
is for the atomic species it is associated with. The other three
5454
represent the reciprocal lattice vector of the grid cell (in units
5455
of the basis vectors of the reciprocal cell). The real number is
5456
the Fourier coefficient in units of the energy scale given in
5457
\textbf{EggboxScale} (see below), normally 1 eV.
5459
The number and choice of Fourier components is free, as well as their
5460
order in the block. One can choose to correct only some species and not
5461
others if, for instance, there is a substantial difference in hardness
5462
of the cores. The 0 0 0 components will add a species-dependent
5463
constant energy per atom. It is thus irrelevant except if comparing
5464
total energies of different calculations, in which case they
5465
have to be considered with care (for instance by putting them all to zero,
5466
i.e. by not introducing them in the list).
5467
The other components average to zero representing no bias in the
5468
total energy comparisons.
5471
energies of the free atoms are put as 0 0 0 coefficients (with
5472
spin polarisation if adequate etc.) the corrected total energy
5473
will be the cohesive energy of the system (per unit cell).
5475
\textit{Example:} For a two species system, this example would give a quite
5476
sufficent set in many instances (the actual values of the Fourier
5477
coefficients are not realistic).
5403
gives again a cubic sampling with half the original side length. It
5404
is not trivial to choose a right set of displacements so as to
5405
maximize the new 'effective' cutoff. It depends on the kind of
5406
cell. It may be automatized in the future, but it is now left to the
5407
user, who introduces the displacements manually through this block.
5409
The quantities which are symmetrized are: ($i$) energy terms that
5410
depend on the grid, ($ii$) forces, ($iii$) stress tensor, and ($iv$)
5413
The symmetrization is performed at the end of every SCF cycle. The
5414
whole cycle is done for the (0,0,0) displacement, and, when the
5415
density matrix is converged, the same (now fixed) density matrix is
5416
used to obtain the desired quantities at the other displacements
5417
(the density matrix itself is \emph{not} symmetrized as it gives a
5418
much smaller egg-box effect). The CPU time needed for each
5419
displacement in the \fdf{GridCellSampling} block is of the order
5420
of one extra SCF iteration.
5424
\begin{fdfentry}{EggboxRemove}[block]
5425
\index{egg-box effect}%
5428
For recovering translational invariance in an approximate way.
5430
It works by substracting from Kohn-Sham's total energy (and forces)
5431
an approximation to the eggbox energy, sum of atomic contributions.
5432
Each atom has a predefined eggbox energy depending on where it sits
5433
on the cell. This atomic contribution is species dependent and is
5434
obviously invariant under grid-cell translations. Each species
5435
contribution is thus expanded in the appropriate Fourier series. It
5436
is important to have a smooth eggbox, for it to be represented by a
5437
few Fourier components. A jagged egg-box (unless very small, which
5438
is then unimportant) is often an indication of a problem with the
5441
In the block there is one line per Fourier component. The first
5442
integer is for the atomic species it is associated with. The other
5443
three represent the reciprocal lattice vector of the grid cell (in
5444
units of the basis vectors of the reciprocal cell). The real number
5445
is the Fourier coefficient in units of the energy scale given in
5446
\fdf{EggboxScale} (see below), normally 1 eV.
5448
The number and choice of Fourier components is free, as well as
5449
their order in the block. One can choose to correct only some
5450
species and not others if, for instance, there is a substantial
5451
difference in hardness of the cores. The 0 0 0 components will add a
5452
species-dependent constant energy per atom. It is thus irrelevant
5453
except if comparing total energies of different calculations, in
5454
which case they have to be considered with care (for instance by
5455
putting them all to zero, i.e. by not introducing them in the list).
5456
The other components average to zero representing no bias in the
5457
total energy comparisons.
5459
If the total energies of the free atoms are put as 0 0 0
5460
coefficients (with spin polarisation if adequate etc.) the corrected
5461
total energy will be the cohesive energy of the system (per unit
5464
\emph{Example:} For a two species system, this example would give a
5465
quite sufficent set in many instances (the actual values of the
5466
Fourier coefficients are not realistic).
5480
5468
%block EggBoxRemove
5481
5469
1 0 0 0 -143.86904
5482
5470
1 0 0 1 0.00031
5803
5782
The complete set of wavefunctions obtained during the last
5804
5783
iteration of the SCF loop will be written to a NetCDF file
5805
\texttt{WFS.nc} if the \fdf{Diag!UseNewDiagk} option is in effect.
5784
\file{WFS.nc} if the \fdf{Diag!UseNewDiagk} option is in effect.
5807
5786
The complete set of wavefunctions obtained during the last
5808
iteration of the SCF loop will be written to SystemLabel.fullBZ.WFSX
5809
\index{WFSX file!fullBZ.WFSX@\texttt{fullBZ.WFSX}}
5810
if the \textbf{COOP.write} option is in effect.
5813
\item[\textbf{WriteEigenvalues}] (\textit{logical}):
5814
\index{WriteEigenvalues@\textbf{WriteEigenvalues}}\index{output!eigenvalues}
5815
If \texttt{.true.} it writes the Hamiltonian eigenvalues for the sampling
5816
$\vec k$ points, in the main output file.
5817
If \texttt{.false.}, it writes them in the file \textit{Systemlabel}.EIG,
5818
which can be used by the \textsc{Eig2DOS}\index{Eig2DOS@\textsc{Eig2DOS}} postprocessing
5819
utility (in the Util/Eig2DOS directory) for obtaining the density of
5820
states.\index{density of states}
5822
\textit{Use:} Only if \textbf{SolutionMethod} is \texttt{diagon}.
5824
\textit{Default value:} \texttt{.false.} (see \textbf{LongOutput})
5787
iteration of the SCF loop will be written to \sysfile{fullBZ.WFSX}
5788
if the \fdf{COOP.Write} option is in effect.
5791
\begin{fdflogicalF}{WriteEigenvalues}
5792
\index{output!eigenvalues}
5794
If \fdftrue\ it writes the Hamiltonian eigenvalues for the sampling
5795
$\vec k$ points, in the main output file. If \fdffalse, it
5796
writes them in the file \sysfile{EIG}, which can be used
5797
by the \program{Eig2DOS}\index{Eig2DOS@\textsc{Eig2DOS}}
5798
postprocessing utility (in the Util/Eig2DOS directory) for obtaining
5799
the density of states.\index{density of states}
5801
\note this option only works for \fdf{SolutionMethod} which
5802
calculates the eigenvalues.
5828
5807
\subsubsection{Occupation of electronic states and Fermi level}
5829
5808
\label{electronic-occupation}
5834
\item[\textbf{OccupationFunction}](\textit{string}):
5835
\index{OccupationFunction@\textbf{OccupationFunction}}
5836
\index{ElectronicTemperature@\textbf{ElectronicTemperature}}
5837
String variable to select the function that determines the occupation
5838
of the electronic states. Two options are available:
5840
\item \texttt{FD}: The usual Fermi-Dirac occupation function is used.
5841
\item \texttt{MP}: The occupation function proposed by Methfessel and
5842
Paxton (Phys. Rev. B, \textbf{40}, 3616 (1989)), is used.
5844
The smearing of the electronic occupations is done, in both cases,
5845
using an energy width defined by the \textbf{ElectronicTemperature}
5846
variable. Note that, while in the case of Fermi-Dirac, the occupations
5847
correspond to the physical ones if the electronic temperature is
5848
set to the physical temperature of the system, this is not the case
5849
in the Methfessel-Paxton function. In this case, the tempeature
5850
is just a mathematical artifact to obtain a more accurate
5851
integration of the physical quantities at a lower cost. In
5852
particular, the Methfessel-Paxton scheme has the advantage
5853
that, even for quite large smearing temperatures, the
5854
obtained energy is very close to the physical energy at T=0.
5855
Also, it allows a much faster convergence with respect to
5856
k-points, specially for metals. Finally, the convergence to
5857
selfconsistency is very much improved (allowing the use
5858
of larger mixing coefficients).
5860
For the Methfessel-Paxton case, one can use relatively large
5861
values for the \textbf{ElectronicTemperature} parameter. How large
5862
depends on the specific system. A guide can be found in the
5863
article by J. Kresse and J. Furthm\"uller, Comp. Mat. Sci.
5864
\textbf{6}, 15 (1996).
5866
If Methfessel-Paxton smearing is used, the order of
5867
the corresponding Hermite polynomial expansion must also be chosen
5868
(see description of variable \textbf{OccupationMPOrder}).
5870
We finally note that, in both cases (FD and MP), once a finite
5871
temperature has been chosen, the relevant energy is not the Kohn-Sham
5872
energy, but the Free energy. In particular, the
5873
atomic forces are derivatives of the Free energy, not the KS
5874
energy. See R. Wentzcovitch \textit{et al.}, Phys. Rev. B \textbf{45},
5875
11372 (1992); S. de Gironcoli, Phys. Rev. B \textbf{51}, 6773 (1995);
5876
J. Kresse and J. Furthm\"uller, Comp. Mat. Sci.
5877
\textbf{6}, 15 (1996),
5881
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{diagon}
5883
\textit{Default value:} \texttt{FD}
5885
\item[\textbf{OccupationMPOrder}](\textit{integer}):
5886
\index{OccupationMPOrder@\textbf{OccupationMPOrder}}
5887
\index{OccupationFunction@\textbf{OccupationFunction}}
5888
Order of the Hermite-Gauss polynomial expansion for the
5889
electronic occupation functions in the Methfessel-Paxton
5890
scheme (see Phys. Rev. B \textbf{40}, 3616 (1989)).
5891
Specially for metals, higher order expansions provide better convergence
5892
to the ground state result, even with larger smearing
5893
temperatures, and provide also better convergence with k-points.
5896
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{diagon}
5897
and \textbf{OccupationFunction} = \texttt{MP}
5899
\textit{Default value:} \texttt{1}
5903
\item[\textbf{ElectronicTemperature}] (\textit{real temperature or energy}):
5904
\index{ElectronicTemperature@\textbf{ElectronicTemperature}}
5905
\index{OccupationFunction@\textbf{OccupationFunction}}
5906
Temperature for Fermi-Dirac or Methfessel-Paxton
5907
distribution. Useful specially for
5908
metals, and to accelerate selfconsistency in some cases.
5910
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{diagon}
5912
\textit{Default value:} \texttt{300.0 K}
5917
\subsubsection{Orbital minimization method (OMM)} \label{SolverOMM}
5919
The OMM is an alternative cubic-scaling solver that uses a minimization
5920
algorithm instead of direct diagonalization to find the occupied subspace.
5921
The main advantage over diagonalization is the possibility of iteratively
5922
reusing the solution from each SCF/MD step as the starting guess of the
5923
following one, thus greatly reducing the time to solution. Typically,
5924
therefore, the first few SCF cycles of the first MD step of a simulation
5925
will be slower than diagonalization, but the rest will be faster. The main
5926
disadvantages are that individual Kohn-Sham eigenvalues are not computed,
5927
and that only a fixed, integer number of electrons at each k point/spin is
5928
allowed. Therefore, only spin-polarized calculations with \textbf{FixSpin}
5929
are allowed, and \textbf{TotalSpin} must be chosen appropriately. For
5930
non-$\Gamma$ point calculations, the number of electrons is set to be equal
5931
at all k points. \textbf{NonCollinearSpin} calculations are not supported at
5811
\begin{fdfentry}{OccupationFunction}[string]<FD>
5813
String variable to select the function that determines the
5814
occupation of the electronic states. Two options are available:
5817
The usual Fermi-Dirac occupation function is used.
5820
The occupation function proposed by Methfessel and
5821
Paxton (Phys. Rev. B, \textbf{40}, 3616 (1989)), is used.
5824
The smearing of the electronic occupations is done, in both cases,
5825
using an energy width defined by the \fdf{ElectronicTemperature}
5826
variable. Note that, while in the case of Fermi-Dirac, the
5827
occupations correspond to the physical ones if the electronic
5828
temperature is set to the physical temperature of the system, this
5829
is not the case in the Methfessel-Paxton function. In this case, the
5830
tempeature is just a mathematical artifact to obtain a more accurate
5831
integration of the physical quantities at a lower cost. In
5832
particular, the Methfessel-Paxton scheme has the advantage that,
5833
even for quite large smearing temperatures, the obtained energy is
5834
very close to the physical energy at $T=0$. Also, it allows a much
5835
faster convergence with respect to $k$-points, specially for
5836
metals. Finally, the convergence to selfconsistency is very much
5837
improved (allowing the use of larger mixing coefficients).
5839
For the Methfessel-Paxton case, one can use relatively large values
5840
for the \fdf{ElectronicTemperature} parameter. How large depends
5841
on the specific system. A guide can be found in the article by
5842
J. Kresse and J. Furthm\"uller, Comp. Mat. Sci. \textbf{6}, 15
5845
If Methfessel-Paxton smearing is used, the order of the
5846
corresponding Hermite polynomial expansion must also be chosen (see
5847
description of variable \fdf{OccupationMPOrder}).
5849
We finally note that, in both cases (FD and MP), once a finite
5850
temperature has been chosen, the relevant energy is not the
5851
Kohn-Sham energy, but the Free energy. In particular, the atomic
5852
forces are derivatives of the Free energy, not the KS energy. See
5853
R. Wentzcovitch \textit{et al.}, Phys. Rev. B \textbf{45}, 11372
5854
(1992); S. de Gironcoli, Phys. Rev. B \textbf{51}, 6773 (1995);
5855
J. Kresse and J. Furthm\"uller, Comp. Mat. Sci. \textbf{6}, 15
5856
(1996), for details.
5860
\begin{fdfentry}{OccupationMPOrder}[integer]<1>
5862
Order of the Hermite-Gauss polynomial expansion for the electronic
5863
occupation functions in the Methfessel-Paxton scheme (see
5864
Phys. Rev. B \textbf{40}, 3616 (1989)). Specially for metals,
5865
higher order expansions provide better convergence to the ground
5866
state result, even with larger smearing temperatures, and provide
5867
also better convergence with k-points.
5869
\note only used if \fdf{OccupationFunction} is \fdf*{MP}.
5874
\begin{fdfentry}{ElectronicTemperature}[temperature/energy]<$300\,\mathrm{K}$>
5876
Temperature for Fermi-Dirac or Methfessel-Paxton
5877
distribution. Useful specially for metals, and to accelerate
5878
selfconsistency in some cases.
5884
\subsubsection{Orbital minimization method (OMM)}
5887
The OMM is an alternative cubic-scaling solver that uses a
5888
minimization algorithm instead of direct diagonalization to find the
5889
occupied subspace. The main advantage over diagonalization is the
5890
possibility of iteratively reusing the solution from each SCF/MD step
5891
as the starting guess of the following one, thus greatly reducing the
5892
time to solution. Typically, therefore, the first few SCF cycles of
5893
the first MD step of a simulation will be slower than diagonalization,
5894
but the rest will be faster. The main disadvantages are that
5895
individual Kohn-Sham eigenvalues are not computed, and that only a
5896
fixed, integer number of electrons at each k point/spin is
5897
allowed. Therefore, only spin-polarized calculations with
5898
\fdf{Spin!Fix} are allowed, and \fdf{Spin!Total} must be chosen
5899
appropriately. For non-$\Gamma$ point calculations, the number of
5900
electrons is set to be equal at all k points. Non-collinear
5901
calculations (see \fdf{Spin}) are not supported at present.
5934
5903
It is important to note that the OMM requires all occupied Kohn-Sham
5935
5904
eigenvalues to be negative; this can be achieved by applying a shift
5936
to the eigenspectrum, controlled by \textbf{ON.eta} (in this case,
5937
\textbf{ON.eta} simply needs to be higher than the HOMO level). If the OMM
5938
exhibits a pathologically slow or unstable convergence, this is almost
5939
certainly due to the fact that the default value of \textbf{ON.eta}
5940
(\texttt{0.0 eV}) is too low, and should be raised by a few eV.
5946
\item[\textbf{OMM.UseCholesky}] (\textit{logical}):
5947
\index{OMM.UseCholesky@\textbf{OMM.UseCholesky}}
5948
Select whether to perform a Cholesky factorization of the generalized eigenvalue
5949
problem; this removes the overlap matrix from the problem but also destroys the
5950
sparsity of the Hamiltonian matrix.
5952
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
5954
\textit{Default value:} \texttt{.true.}
5956
\item[\textbf{OMM.Use2D}] (\textit{logical}):
5957
\index{OMM.Use2D@\textbf{OMM.Use2D}}
5958
Select whether to use a 2D data decomposition of the matrices for parallel
5959
calculations. This generally leads to superior scaling for large numbers of MPI
5962
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
5964
\textit{Default value:} \texttt{.true.}
5966
\item[\textbf{OMM.UseSparse}] (\textit{logical}):
5967
\index{OMM.UseSparse@\textbf{OMM.UseSparse}}
5968
Select whether to make use of the sparsity of the Hamiltonian and overlap
5969
matrices where possible when performing matrix-matrix multiplications (these
5970
operations are thus reduced from $O(N^3)$ to $O(N^2)$ without loss of
5973
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}; not compatible with
5974
\textbf{OMM.UseCholesky}, \textbf{OMM.Use2D}, or non-$\Gamma$ point calculations
5976
\textit{Default value:} \texttt{.false.}
5978
\item[\textbf{OMM.Precon}] (\textit{integer}):
5979
\index{OMM.Precon@\textbf{OMM.Precon}}
5980
Number of SCF steps for \em{all} MD steps for which to apply a preconditioning
5981
scheme based on the overlap and kinetic energy matrices; for negative values the
5982
preconditioning is always applied. Preconditioning is usually essential for fast
5983
and accurate convergence (note, however, that it is not needed if a Cholesky
5984
factorization is performed; in such cases this variable will have no effect on
5987
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}; not used with
5988
\textbf{OMM.UseCholesky}
5990
\textit{Default value:} \texttt{-1}
5992
\item[\textbf{OMM.PreconFirstStep}] (\textit{integer}):
5993
\index{OMM.PreconFirstStep@\textbf{OMM.PreconFirstStep}}
5994
Number of SCF steps in the \emph{first} MD step for which to apply the preconditioning
5995
scheme; if present, this will overwrite the value given in \textbf{OMM.Precon} for the
5998
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}; not used with
5999
\textbf{OMM.UseCholesky}
6001
\textit{Default value:} \texttt{OMM.Precon}
6003
\item[\textbf{OMM.Diagon}] (\textit{integer}):
6004
\index{OMM.Diagon@\textbf{OMM.Diagon}}
6005
Number of SCF steps for \emph{all} MD steps for which to use a standard diagonalization
6006
before switching to the OMM; for negative values diagonalization is always used, and
6007
so the calculation is effectively equivalent to \textbf{SolutionMethod} = \texttt{diagon}.
6008
In general, selecting the first few SCF steps can speed up the calculation by removing
6009
the costly initial minimization (at present this works best for $\Gamma$ point
6012
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
6014
\textit{Default value:} \texttt{0}
6016
\item[\textbf{OMM.DiagonFirstStep}] (\textit{integer}):
6017
\index{OMM.DiagonFirstStep@\textbf{OMM.DiagonFirstStep}}
6018
Number of SCF steps in the \emph{first} MD step for which to use a standard diagonalization
6019
before switching to the OMM; if present, this will overwrite the value given in
6020
\textbf{OMM.Diagon} for the first MD step only.
6022
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
6024
\textit{Default value:} \texttt{OMM.Diagon}
6026
\item[\textbf{OMM.BlockSize}] (\textit{integer}):
6027
\index{OMM.BlockSize@\textbf{OMM.BlockSize}}
6028
Blocksize used for distributing the elements of the matrix over MPI processes. Specifically,
6029
this variable controls the dimension relating to the trial orbitals used in the minimization
6030
(equal to the number of occupied states at each k point/spin); the equivalent variable for
6031
the dimension relating to the underlying basis orbitals is controlled by \textbf{BlockSize}.
6033
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
6035
\textit{Default value:} chosen by \siesta\
6037
\item[\textbf{OMM.TPreconScale}] (\textit{real energy}):
6038
\index{OMM.TPreconScale@\textbf{OMM.TPreconScale}}
6039
Scale of the kinetic energy preconditioning (see C.~K.~Gan \emph{et al.}, Comput. Phys. Commun.
6040
\textbf{134}, 33 (2001)). A smaller value indicates more aggressive kinetic energy preconditioning,
6041
while an infinite value indicates no kinetic energy preconditioning. In general, the kinetic
6042
energy preconditioning is much less important than the tensorial correction brought about by
6043
the overlap matrix, and so this value will have fairly little impact on the overall performace
6044
of the preconditioner; however, too aggressive kinetic energy preconditioning can have a
6045
detrimental effect on performance and accuracy.
6047
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
6049
\textit{Default value:} \texttt{10.0 Ry}
6051
\item[\textbf{OMM.RelTol}] (\textit{real}):
6052
\index{OMM.RelTol@\textbf{OMM.RelTol}}
6053
Relative tolerance in the conjugate gradients minimization of the Kohn-Sham band energy (see
6056
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
6058
\textit{Default value:} $10^{-9}$
6060
\item[\textbf{OMM.Eigenvalues}] (\textit{logical}):
6061
\index{OMM.Eigenvalues@\textbf{OMM.Eigenvalues}}
6062
Select whether to perform a diagonalization at the end of each MD step to obtain the Kohn-Sham
6065
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
6067
\textit{Default value:} \texttt{.false.}
6069
\item[\textbf{OMM.WriteCoeffs}] (\textit{logical}):
6070
\index{OMM.WriteCoeffs@\textbf{OMM.WriteCoeffs}}
6071
Select whether to write the coefficients of the solution orbitals to file at the end of each MD
6074
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
6076
\textit{Default value:} \texttt{.false.}
6078
\item[\textbf{OMM.ReadCoeffs}] (\textit{logical}):
6079
\index{OMM.ReadCoeffs@\textbf{OMM.ReadCoeffs}}
6080
Select whether to read the coefficients of the solution orbitals from file at the beginning of
6081
a new calculation. Useful for restarting an interrupted calculation, especially when used in
6082
conjuction with \textbf{DM.UseSaveDM}. Note that the same number of MPI processes and values of
6083
\textbf{OMM.Use2D}, \textbf{OMM.BlockSize}, and \textbf{BlockSize} must be used when restarting.
6085
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
6087
\textit{Default value:} \texttt{.false.}
6089
\item[\textbf{OMM.LongOutput}] (\textit{logical}):
6090
\index{OMM.LongOutput@\textbf{OMM.LongOutput}}
6091
Select whether to output detailed information of the conjugate gradients minimization for each
6094
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OMM}
6096
\textit{Default value:} \texttt{.false.}
6101
\subsubsection{Order(N) calculations} \label{SolverON}
5905
to the eigenspectrum, controlled by \fdf{ON.eta} (in this case,
5906
\fdf{ON.eta} simply needs to be higher than the HOMO level). If the
5907
OMM exhibits a pathologically slow or unstable convergence, this is
5908
almost certainly due to the fact that the default value of
5909
\fdf{ON.eta} (\fdf*{0.0 eV}) is too low, and should be raised by
5912
\begin{fdflogicalT}{OMM!UseCholesky}
5914
Select whether to perform a Cholesky factorization of the
5915
generalized eigenvalue problem; this removes the overlap matrix from
5916
the problem but also destroys the sparsity of the Hamiltonian
5921
\begin{fdflogicalT}{OMM!Use2D}
5923
Select whether to use a 2D data decomposition of the matrices for
5924
parallel calculations. This generally leads to superior scaling for
5925
large numbers of MPI processes.
5929
\begin{fdflogicalF}{OMM!UseSparse}
5931
Select whether to make use of the sparsity of the Hamiltonian and
5932
overlap matrices where possible when performing matrix-matrix
5933
multiplications (these operations are thus reduced from $O(N^3)$ to
5934
$O(N^2)$ without loss of accuracy).
5936
\note not compatible with \fdf{OMM!UseCholesky},
5937
\fdf{OMM!Use2D}, or non-$\Gamma$ point calculations
5941
\begin{fdfentry}{OMM!Precon}[integer]<-1>
5943
Number of SCF steps for \emph{all} MD steps for which to apply a
5944
preconditioning scheme based on the overlap and kinetic energy
5945
matrices; for negative values the preconditioning is always
5946
applied. Preconditioning is usually essential for fast and accurate
5947
convergence (note, however, that it is not needed if a Cholesky
5948
factorization is performed; in such cases this variable will have no
5949
effect on the calculation).
5951
\note cannot be used with \fdf{OMM!UseCholesky}.
5956
\begin{fdfentry}{OMM!PreconFirstStep}[integer]<\fdfvalue{OMM!Precon}>
5958
Number of SCF steps in the \emph{first} MD step for which to apply
5959
the preconditioning scheme; if present, this will overwrite the
5960
value given in \fdf{OMM!Precon} for the first MD step only.
5964
\begin{fdfentry}{OMM!Diagon}[integer]<0>
5966
Number of SCF steps for \emph{all} MD steps for which to use a
5967
standard diagonalization before switching to the OMM; for negative
5968
values diagonalization is always used, and so the calculation is
5969
effectively equivalent to \fdf{SolutionMethod} \fdf*{diagon}.
5970
In general, selecting the first few SCF steps can speed up the
5971
calculation by removing the costly initial minimization (at present
5972
this works best for $\Gamma$ point calculations).
5976
\begin{fdfentry}{OMM!DiagonFirstStep}[integer]<\fdfvalue{OMM!Diagon}>
5978
Number of SCF steps in the \emph{first} MD step for which to use a
5979
standard diagonalization before switching to the OMM; if present,
5980
this will overwrite the value given in \fdf{OMM!Diagon} for the
5985
\begin{fdfentry}{OMM!BlockSize}[integer]<\fdfvalue{BlockSize}>
5987
Blocksize used for distributing the elements of the matrix over MPI
5988
processes. Specifically, this variable controls the dimension
5989
relating to the trial orbitals used in the minimization (equal to
5990
the number of occupied states at each k point/spin); the equivalent
5991
variable for the dimension relating to the underlying basis orbitals
5992
is controlled by \fdf{BlockSize}.
5996
\begin{fdfentry}{OMM!TPreconScale}[energy]<$10\,\mathrm{Ry}$>
5998
Scale of the kinetic energy preconditioning (see C.~K.~Gan \emph{et
5999
al.}, Comput. Phys. Commun. \textbf{134}, 33 (2001)). A smaller
6000
value indicates more aggressive kinetic energy preconditioning,
6001
while an infinite value indicates no kinetic energy
6002
preconditioning. In general, the kinetic energy preconditioning is
6003
much less important than the tensorial correction brought about by
6004
the overlap matrix, and so this value will have fairly little impact
6005
on the overall performace of the preconditioner; however, too
6006
aggressive kinetic energy preconditioning can have a detrimental
6007
effect on performance and accuracy.
6011
\begin{fdfentry}{OMM!RelTol}[real]<$10^{-9}$>
6013
Relative tolerance in the conjugate gradients minimization of the
6014
Kohn-Sham band energy (see \fdf{ON!Etol}).
6018
\begin{fdflogicalF}{OMM!Eigenvalues}
6020
Select whether to perform a diagonalization at the end of each MD
6021
step to obtain the Kohn-Sham eigenvalues.
6025
\begin{fdflogicalF}{OMM!WriteCoeffs}
6027
Select whether to write the coefficients of the solution orbitals to
6028
file at the end of each MD step.
6032
\begin{fdflogicalF}{OMM!ReadCoeffs}
6034
Select whether to read the coefficients of the solution orbitals
6035
from file at the beginning of a new calculation. Useful for
6036
restarting an interrupted calculation, especially when used in
6037
conjuction with \fdf{DM.UseSaveDM}. Note that the same number of
6038
MPI processes and values of \fdf{OMM!Use2D},
6039
\fdf{OMM!BlockSize}, and \fdf{BlockSize} must be used when
6045
\begin{fdflogicalF}{OMM!LongOutput}
6047
Select whether to output detailed information of the conjugate
6048
gradients minimization for each SCF step.
6053
\subsubsection{Order(N) calculations}
6103
6056
The Ordern(N) subsystem is quite fragile and only works for systems
6104
6057
with clearly separated occupied and empty states. Note also that the
6110
6063
ongoing to clean and automate the O(N) process, to make the solver
6111
6064
more user-friendly and robust.
6116
\item[\textbf{ON.functional}] (\textit{string}):
6117
\index{ON.functional@\textbf{ON.functional}}
6118
Choice of order-N minimization functionals:
6120
\item \texttt{Kim}:\index{Kim@\texttt{Kim}}
6121
Functional of Kim, Mauri and Galli, PRB 52, 1640 (1995).
6122
\item \texttt{Ordejon-Mauri}:\index{Ordejon-Mauri@\texttt{Ordejon-Mauri}}
6123
Functional of Ordej\'on et al, or Mauri et al, see PRB 51, 1456 (1995).
6124
The number of localized wave functions (LWFs) used must coincide with
6125
$N_{el}/2$ (unless spin polarized).
6126
For the initial assignment of LWF centers to atoms, atoms
6127
with even number of electrons, n, get n/2 LWFs. Odd atoms
6128
get (n+1)/2 and (n-1)/2 in an alternating sequence, ir order
6129
of appearance (controlled by the input in the atomic coordinates block).
6131
\item \texttt{files}:\index{files (ON.functional)@\texttt{files} (ON.functional)}
6132
Reads localized-function information from a file and
6133
chooses automatically the functional to be used.
6136
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{ordern}
6138
\textit{Default value:} \texttt{Kim}
6140
\item[\textbf{ON.MaxNumIter}] (\textit{integer}):
6141
\index{ON.MaxNumIter@\textbf{ON.MaxNumIter}}
6142
Maximum number of iterations
6143
in the conjugate minimization of the electronic
6144
energy, in each SCF cycle.
6146
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}
6148
\textit{Default value:} \texttt{1000}
6150
\item[\textbf{ON.etol}] (\textit{real}):
6151
\index{ON.etol@\textbf{ON.etol}}
6152
Relative-energy tolerance in the conjugate minimization of the electronic
6153
energy. The minimization finishes if
6154
\hspace{0.2truecm} $2 (E_n - E_{n-1}) / (E_n + E_{n-1}) \leq $ ON.etol.
6156
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}
6158
\textit{Default value:} $10^{-8}$
6160
\item[\textbf{ON.eta}] (\textit{real energy}):
6161
\index{ON.eta@\textbf{ON.eta}}
6162
Fermi level parameter of Kim
6163
\textit{et al.}. This should be in the energy gap, and tuned to obtain
6164
the correct number of electrons. If the calculation is spin polarised,
6165
then separate Fermi levels for each spin can be specified.
6167
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}
6169
\textit{Default value:} \texttt{0.0 eV}
6171
\item[\textbf{ON.eta\_alpha}] (\textit{real energy}):
6172
\index{ON.eta\_alpha@\textbf{ON.eta\_alpha}}
6173
Fermi level parameter of Kim \textit{et al.} for alpha spin electrons.
6174
This should be in the energy gap, and tuned to obtain
6175
the correct number of electrons. Note that if the Fermi
6176
level is not specified individually for each spin then the
6177
same global eta will be used.
6179
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}
6181
\textit{Default value:} \texttt{0.0 eV}
6183
\item[\textbf{ON.eta\_beta}] (\textit{real energy}):
6184
\index{ON.eta\_beta@\textbf{ON.eta\_beta}}
6185
Fermi level parameter of Kim \textit{et al.} for beta spin electrons.
6186
This should be in the energy gap, and tuned to obtain
6187
the correct number of electrons. Note that if the Fermi
6188
level is not specified individually for each spin then the
6189
same global eta will be used.
6191
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}
6193
\textit{Default value:} \texttt{0.0 eV}
6195
\item[\textbf{ON.RcLWF}] (\textit{real legth}):
6196
\index{ON.RcLWF@\textbf{ON.RcLWF}}\index{Localized Wave Functions}
6197
Localization redius for the Localized Wave Functions (LWF's).
6199
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}
6201
\textit{Default value:} \texttt{9.5 Bohr}
6203
\item[\textbf{ON.ChemicalPotential}] (\textit{logical}):
6204
\index{ON.ChemicalPotential@\textbf{ON.ChemicalPotential}}
6205
\index{Chemical Potential}
6206
Specifies whether to calculate an order-\textit{N} estimate of the
6207
Chemical Potential, by the projection method
6208
(Goedecker and Teter, PRB \textbf{51}, 9455 (1995);
6209
Stephan, Drabold and Martin, PRB \textbf{58}, 13472
6211
done by expanding the Fermi function (or density matrix)
6212
at a given temperature, by means of Chebyshev
6213
polynomials\index{Chebyshev Polynomials}, and imposing a
6214
real space truncation on the density matrix.
6215
To obtain a realistic estimate, the temperature
6216
should be small enough (typically, smaller than
6217
the energy gap), the localization range large enough
6218
(of the order of the one you would use for the Localized Wannier
6219
Functions), and the order of the polynomial expansion
6220
sufficiently large (how large depends on the temperature;
6223
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}.
6225
\textit{Default value:} \texttt{.false.}
6227
\textbf{Note:} This option does not work in parallel. An alternative
6228
is to obtain the approximate value of the chemical potential using
6229
an initial diagonalization.
6232
\item[\textbf{ON.ChemicalPotentialUse}] (\textit{logical}):
6233
\index{ON.ChemicalPotentialUse@\textbf{ON.ChemicalPotentialUse}}
6234
\index{Chemical Potential}
6235
Specifies whether to use the calculated estimate of the
6236
Chemical Potential, instead of the parameter
6237
\textbf{ON.eta}\index{ON.eta@\textbf{ON.eta}}
6238
for the order-\textit{N} energy functional minimization.
6239
This is useful if you do not know the position
6240
of the Fermi level, typically in the beginning
6241
of an order-\emph{N} run.
6243
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}.
6244
Overrides the value of \textbf{ON.eta}.
6245
Overrides the value of \textbf{ON.ChemicalPotential}, setting
6246
it to \texttt{.true.}.
6248
\textit{Default value:} \texttt{.false.}
6250
\textbf{Note:} This option does not work in parallel. An alternative
6251
is to obtain the approximate value of the chemical potential using
6252
an initial diagonalization.
6254
\item[\textbf{ON.ChemicalPotentialRc}] (\textit{real length}):
6255
\index{ON.ChemicalPotentialRc@\textbf{ON.ChemicalPotentialRc}}
6256
\index{Chemical Potential}
6257
Defines the cutoff radius for the density matrix or Fermi
6258
operator in the calculation of the estimate of the
6261
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}
6262
and \textbf{ON.ChemicalPotential} or \textbf{ON.ChemicalPotentialUse}
6265
\textit{Default value:} \texttt{9.5 Bohr}.
6267
\item[\textbf{ON.ChemicalPotentialTemperature}] (\textit{real temperature
6269
\index{ON.ChemicalPotentialTemperature@\textbf{ON.ChemicalPotentialTemperature}}
6270
\index{Chemical Potential}
6271
Defines the temperature to be used in the Fermi function expansion
6272
in the calculation of the estimate of the Chemical Potential.
6273
To have an accurate results, this temperature should be smaller
6274
than the gap of the system.
6276
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN},
6277
and \textbf{ON.ChemicalPotential} or \textbf{ON.ChemicalPotentialUse} =
6280
\textit{Default value:} \texttt{0.05 Ry}.
6282
\item[\textbf{ON.ChemicalPotentialOrder}] (\textit{integer}):
6283
\index{ON.ChemicalPotentialOrder@\textbf{ON.ChemicalPotentialOrder}}
6284
\index{Chemical Potential}
6285
Order of the Chebishev expansion to calculate the estimate
6286
of the Chemical Potential.
6288
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN},
6289
and \textbf{ON.ChemicalPotential} or \textbf{ON.ChemicalPotentialUse} =
6292
\textit{Default value:} \texttt{100}
6294
\item[\textbf{ON.LowerMemory}] (\textit{logical}):
6295
\index{ON.LowerMemory@\textbf{ON.LowerMemory}}\index{Lower order N memory}
6296
If .true., then a slightly reduced memory algorithm is used in the
6297
3-point line search during the order N minimisation. Only affects
6300
\textit{Use:} Used only if \textbf{SolutionMethod} = \texttt{OrderN}
6302
\textit{Default value:} \texttt{.false.}
6307
\textbf{Output of localized wavefunctions}
6066
\begin{fdfentry}{ON.functional}[string]<Kim>
6068
Choice of order-N minimization functionals:
6071
Functional of Kim, Mauri and Galli, PRB 52, 1640 (1995).
6073
\option[Ordejon-Mauri]%
6074
Functional of Ordej\'on et al, or Mauri et al, see PRB 51, 1456
6075
(1995). The number of localized wave functions (LWFs) used must
6076
coincide with $N_{el}/2$ (unless spin polarized). For the initial
6077
assignment of LWF centers to atoms, atoms with even number of
6078
electrons, $n$, get $n/2$ LWFs. Odd atoms get $(n+1)/2$ and
6079
$(n-1)/2$ in an alternating sequence, ir order of appearance
6080
(controlled by the input in the atomic coordinates block).
6083
Reads localized-function information from a file and chooses
6084
automatically the functional to be used.
6089
\begin{fdfentry}{ON.MaxNumIter}[integer]<1000>
6091
Maximum number of iterations in the conjugate minimization of the
6092
electronic energy, in each SCF cycle.
6096
\begin{fdfentry}{ON.Etol}[real]<$10^{-8}$>
6098
Relative-energy tolerance in the conjugate minimization of the
6099
electronic energy. The minimization finishes if
6100
$2 (E_n - E_{n-1}) / (E_n + E_{n-1}) \leq $ \fdf{ON.Etol}.
6104
\begin{fdfentry}{ON.eta}[energy]<$0\,\mathrm{eV}$>
6106
Fermi level parameter of Kim \textit{et al.}. This should be in the
6107
energy gap, and tuned to obtain the correct number of electrons. If
6108
the calculation is spin polarised, then separate Fermi levels for
6109
each spin can be specified.
6113
\begin{fdfentry}{ON.eta.alpha}[energy]<$0\,\mathrm{eV}$>
6115
Fermi level parameter of Kim \textit{et al.} for alpha spin
6116
electrons. This should be in the energy gap, and tuned to obtain
6117
the correct number of electrons. Note that if the Fermi level is not
6118
specified individually for each spin then the same global eta will
6123
\begin{fdfentry}{ON.eta.beta}[energy]<$0\,\mathrm{eV}$>
6125
Fermi level parameter of Kim \textit{et al.} for beta spin
6126
electrons. This should be in the energy gap, and tuned to obtain
6127
the correct number of electrons. Note that if the Fermi level is not
6128
specified individually for each spin then the same global eta will
6133
\begin{fdfentry}{ON.RcLWF}[length]<$9.5\,\mathrm{Bohr}$>
6134
\index{Localized Wave Functions}
6136
Localization redius for the Localized Wave Functions (LWF's).
6140
\begin{fdflogicalF}{ON.ChemicalPotential}
6141
\index{Chemical Potential}
6143
Specifies whether to calculate an order-$N$ estimate of the Chemical
6144
Potential, by the projection method (Goedecker and Teter, PRB
6145
\textbf{51}, 9455 (1995); Stephan, Drabold and Martin, PRB
6146
\textbf{58}, 13472 (1998)). This is done by expanding the Fermi
6147
function (or density matrix) at a given temperature, by means of
6148
Chebyshev polynomials\index{Chebyshev Polynomials}, and imposing a
6149
real space truncation on the density matrix. To obtain a realistic
6150
estimate, the temperature should be small enough (typically, smaller
6151
than the energy gap), the localization range large enough (of the
6152
order of the one you would use for the Localized Wannier Functions),
6153
and the order of the polynomial expansion sufficiently large (how
6154
large depends on the temperature; typically, 50-100).
6157
\note this option does not work in parallel. An alternative is to
6158
obtain the approximate value of the chemical potential using an
6159
initial diagonalization.
6164
\begin{fdflogicalF}{ON.ChemicalPotential.Use}
6165
\index{Chemical Potential}
6167
Specifies whether to use the calculated estimate of the Chemical
6168
Potential, instead of the parameter
6169
\fdf{ON.eta} for the
6170
order-$N$ energy functional minimization. This is useful if
6171
you do not know the position of the Fermi level, typically in the
6172
beginning of an order-$N$ run.
6174
\note this overrides the value of \fdf{ON.eta} and \fdf{ON.ChemicalPotential}.
6175
Also, this option does not work in parallel. An alternative
6176
is to obtain the approximate value of the chemical potential using
6177
an initial diagonalization.
6181
\begin{fdfentry}{ON.ChemicalPotential.Rc}[length]<$9.5\,\mathrm{Bohr}$>
6182
\index{Chemical Potential}
6184
Defines the cutoff radius for the density matrix or Fermi operator
6185
in the calculation of the estimate of the Chemical Potential.
6189
\begin{fdfentry}{ON.ChemicalPotential.Temperature}[temperature/energy]<$0.05\,\mathrm{Ry}$>
6190
\index{Chemical Potential}
6192
Defines the temperature to be used in the Fermi function expansion
6193
in the calculation of the estimate of the Chemical Potential. To
6194
have an accurate results, this temperature should be smaller than
6195
the gap of the system.
6199
\begin{fdfentry}{ON.ChemicalPotential.Order}[integer]<$100$>
6200
\index{Chemical Potential}
6202
Order of the Chebishev expansion to calculate the estimate of the
6207
\begin{fdflogicalF}{ON.LowerMemory}
6208
\index{Lower order N memory}
6210
If \fdftrue, then a slightly reduced memory algorithm is used in the
6211
3-point line search during the order N minimisation. Only affects
6217
\paragraph{Output of localized wavefunctions}
6309
6218
\index{Localized Wave Functions}
6310
6220
At the end of each conjugate gradient minimization of the energy
6311
6221
functional, the LWF's are stored on disk. These can be used as an
6312
6222
input for the same system in a restart, or in case something goes
7171
7041
There are two options to obtain the partial density of states
7172
7042
\begin{itemize}
7173
7043
\item Using the options below
7174
\item Using the \texttt{Util/COOP/mprop} program for the off-line analysis of
7044
\item Using the \program{Util/COOP/mprop} program for the off-line analysis of
7175
7045
the electronic structure in PDOS mode. This method allows the
7176
7046
flexibility of specifying energy ranges, orbitals, and resolutions
7177
7047
at will, without re-running \siesta. See Sec.~\ref{sec:coop}.
7183
\item[\textbf{ProjectedDensityOfStates}] (\textit{block}):
7184
\index{ProjectedDensityOfStates@\textbf{ProjectedDensityOfStates}}
7185
\index{output!projected density of states}
7187
Instructs to write the Total Density Of States (Total DOS) and the
7188
Projected Density Of States (PDOS) on the basis orbitals,
7189
between two given energies,
7190
in files \texttt{SystemLabel}.DOS and
7191
\texttt{SystemLabel}.PDOS, respectively.
7192
The block must be a single line with the energies of the range for
7194
(relative to the program's zero, i.e. the same as the eigenvalues
7195
printed by the program), the peak width (an energy) for broadening
7196
the eigenvalues, the number of points in the energy window,
7197
and the energy units.
7050
\begin{fdfentry}{ProjectedDensityOfStates}[block]
7051
\index{output!projected density of states}
7053
Instructs to write the Total Density Of States (Total DOS) and the
7054
Projected Density Of States (PDOS) on the basis orbitals, between
7055
two given energies, in files \sysfile{DOS} and \sysfile{PDOS},
7056
respectively. The block must be a single line with the energies of
7057
the range for PDOS projection, (relative to the program's zero,
7058
i.e. the same as the eigenvalues printed by the program), the peak
7059
width (an energy) for broadening the eigenvalues, the number of
7060
points in the energy window, and the energy units. An example is:
7201
7062
%block ProjectedDensityOfStates
7202
7063
-20.00 10.00 0.200 500 eV
7203
7064
%endblock ProjectedDensityOfStates
7206
By default the projected density of states is generated for the same
7207
grid of points in reciprocal space as used for the SCF calculation.
7208
However, a separate set of K-points, usually on a finer grid, can
7209
be generated using one of the options \textbf{PDOS.kgrid\_cutoff} or
7210
\textbf{PDOS.kgrid\_Monkhorst\_Pack}. The format of these options is
7211
exactly the same as for \textbf{kgrid\_cutoff} and
7212
\textbf{kgrid\_Monkhorst\_Pack}, respectively. Note that if a gamma
7213
point calculation is being used in the SCF part, especially as part
7214
of a geometry optimisation, and this is then to
7215
be run with a grid of K-points for the PDOS calculation it is more
7216
efficient to run the SCF phase first and then restart to perform the
7217
PDOS evaluation using the density matrix saved from the SCF phase.
7219
\textit{Use:} The two energies of the range must be ordered, with lowest
7222
\textit{Output:} The Total DOS is dumped into a file
7223
called \texttt{SystemLabel}.DOS. The format of this file is:
7225
Energy value, Total DOS (spin up), Total DOS (spin down)
7227
The Projected Density Of States for all the orbitals in the unit cell
7228
is dumped sequentially into a file called
7229
\texttt{SystemLabel}.PDOS. This file is structured using spacing and
7230
xml tags. A machine-readable (but not very human readable) xml file
7231
\texttt{pdos.xml} is also produced. Both can be processed by the
7232
program in \program{Util/pdosxml}. The \sysfile{PDOS} file can be
7233
processed by utilites in \program{Util/Contrib/APostnikov}.
7235
In all cases, the units for the DOS are (number of states/eV), and the
7236
Total DOS, $g \left(\epsilon\right)$, is normalized as follows:
7239
\int_{-\infty}^{+\infty} g \left(\epsilon\right) d\epsilon =
7240
\text{number of basis orbitals in unit cell}
7244
\textit{Default value:} PDOS not calculated nor written.
7067
By default the projected density of states is generated for the same
7068
grid of points in reciprocal space as used for the SCF calculation.
7069
However, a separate set of K-points, usually on a finer grid, can be
7070
generated using one of the options \fdf{PDOS.kgrid.Cutoff} or
7071
\fdf{PDOS.kgrid.MonkhorstPack}. The format of these options is
7072
exactly the same as for \fdf{kgrid.Cutoff} and
7073
\fdf{kgrid.MonkhorstPack}, respectively. Note that if a gamma
7074
point calculation is being used in the SCF part, especially as part
7075
of a geometry optimisation, and this is then to be run with a grid
7076
of K-points for the PDOS calculation it is more efficient to run the
7077
SCF phase first and then restart to perform the PDOS evaluation
7078
using the density matrix saved from the SCF phase.
7080
\note the two energies of the range must be ordered, with lowest
7083
The total DOS is stored in a file called \sysfile{DOS}. The format
7085
\begin{shellexample}
7086
Energy value, Total DOS (spin up), Total DOS (spin down)
7089
The Projected Density Of States for all the orbitals in the unit
7090
cell is dumped sequentially into a file called
7091
\sysfile{PDOS}. This file is structured using spacing and
7092
xml tags. A machine-readable (but not very human readable) xml file
7093
\file{pdos.xml} is also produced. Both can be processed by the
7094
program in \program{Util/pdosxml}. The \sysfile{PDOS} file can be
7095
processed by utilites in \program{Util/Contrib/APostnikov}.
7097
In all cases, the units for the DOS are (number of states/eV), and the
7098
Total DOS, $g(\epsilon)$, is normalized as follows:
7100
\int_{-\infty}^\infty g (\epsilon) d\epsilon =
7101
\text{number of basis orbitals in unit cell}
7248
7108
\subsubsection{Local density of states}
7250
7110
The LDOS is formally the DOS weighted by the amplitude of the
7251
7111
corresponding wavefunctions at different points in space, and is then
7252
7112
a function of energy and position. \siesta\ can output the LDOS
7253
integrated over a range of energies. This information can be used to
7113
integrated over a range of energies. This information can be used to
7254
7114
obtain simple STM images in the Tersoff-Hamann approximation (See
7255
7115
\program{Util/STM/simple-stm}).
7258
\item[\textbf{LocalDensityOfStates}] (\textit{block}):
7259
\index{LocalDensityOfStates@\textbf{LocalDensityOfStates}}
7260
\index{output!local density of states}
7261
Instructs to write the LDOS, integrated between two given energies,
7262
at the mesh used by DHSCF,
7263
in file \texttt{SystemLabel}.LDOS. This file can be read by routine IORHO,
7264
which may be used by an application program in later versions.
7265
The block must be a single line with the energies of the range for
7267
(relative to the program's zero, i.e. the same as the eigenvalues
7268
printed by the program) and their units.
7117
\begin{fdfentry}{LocalDensityOfStates}[block]
7118
\index{output!local density of states}
7120
Instructs to write the LDOS, integrated between two given energies,
7121
at the mesh used by DHSCF, in file \sysfile{LDOS}. This file can be
7122
read by routine IORHO, which may be used by an application program
7123
in later versions. The block must be a single line with the
7124
energies of the range for LDOS integration (relative to the
7125
program's zero, i.e. the same as the eigenvalues printed by the
7126
program) and their units. An example is:
7272
7128
%block LocalDensityOfStates
7274
7130
%endblock LocalDensityOfStates
7277
\textit{Use:} The two energies of the range must be ordered,
7279
File \texttt{SystemLabel}.LDOS is only written, not read, by siesta.
7281
\textit{Default value:} LDOS not calculated nor written.
7283
Note that the LDOS can also be computed using the selected-inversion
7284
technology within the PEXSI solver (see Sec.~\ref{pexsi-ldos}).
7133
\note the two energies of the range must be ordered, with lowest
7288
7140
\subsection{Options for chemical analysis}
7290
7142
\subsubsection{Mulliken charges and overlap populations}
7295
\item[\textbf{WriteMullikenPop}] (\textit{integer}):
7296
\index{WriteMullikenPop@\textbf{WriteMullikenPop}}\index{Mulliken population
7297
analysis}\index{output!Mulliken analysis}
7298
It determines the level of Mulliken population analysis printed:
7301
\item 1 = atomic and orbital charges
7302
\item 2 = 1 + atomic overlap pop.
7303
\item 3 = 2 + orbital overlap pop.
7305
The order of the orbitals in the population lists is defined
7306
by the order of atoms. For each atom, populations for PAO orbitals and
7307
double-$z$, triple-$z$, etc... derived from them are displayed first for
7308
all the angular momenta. Then, populations for perturbative polarization
7309
orbitals are written.
7310
Within a $l$-shell be aware that the order is not
7311
conventional, being $y$, $z$, $x$ for $p$ orbitals, and
7312
$xy$, $yz$, $z^2$, $xz$, and $x^2-y^2$ for $d$ orbitals.
7314
\textit{Default value:} \texttt{0} (see \textbf{LongOutput})
7317
\item[\textbf{MullikenInSCF}] (\textit{logical}):
7318
\index{MullikenInSCF@\textbf{MullikenInSCF}}
7319
If true, the Mulliken populations will be written for every SCF step
7320
at the level of detail specified in \textbf{WriteMullikenPop}. Useful
7321
when dealing with SCF problems, otherwise too verbose.
7323
\textit{Default value:} \texttt{.false.}
7144
\begin{fdfentry}{WriteMullikenPop}[integer]<0>
7145
\index{Mulliken population analysis}%
7146
\index{output!Mulliken analysis}
7148
It determines the level of Mulliken population analysis printed:
7154
atomic and orbital charges
7157
atomic, orbital and atomic overlap populations
7160
atomic, orbital, atomic overlap and orbital overlap populations
7163
The order of the orbitals in the population lists is defined by the
7164
order of atoms. For each atom, populations for PAO orbitals and
7165
double-$z$, triple-$z$, etc... derived from them are displayed first
7166
for all the angular momenta. Then, populations for perturbative
7167
polarization orbitals are written. Within a $l$-shell be aware that
7168
the order is not conventional, being $y$, $z$, $x$ for $p$ orbitals,
7169
and $xy$, $yz$, $z^2$, $xz$, and $x^2-y^2$ for $d$ orbitals.
7174
\begin{fdflogicalF}{MullikenInSCF}
7176
If \fdftrue, the Mulliken populations will be written for every SCF
7177
step at the level of detail specified in
7178
\fdf{WriteMullikenPop}. Useful when dealing with SCF problems,
7179
otherwise too verbose.
7327
7185
\subsubsection{Voronoi and Hirshfeld atomic population analysis}
7332
\item[\textbf{WriteHirshfeldPop}] (\textit{logical}):
7333
\index{WriteHirshfeldPop@\textbf{WriteHirshfeldPop}}\index{Hirshfeld population
7334
analysis}\index{output!Hirshfeld analysis}
7336
If \texttt{.true.}, the program calculates and prints the Hirshfeld
7337
``net'' atomic populations on each atom in the system. For a
7338
definition of the Hirshfeld charges, see Hirshfeld, Theo Chem Acta
7339
\textbf{44}, 129 (1977) and Fonseca et al, J. Comp. Chem. \textbf{25},
7340
189 (2003). Hirshfeld charges are more reliable than Mulliken
7341
charges, specially for large basis sets. The number printed is the
7342
total net charge of the atom: the variation from the neutral charge,
7343
in units of $|e|$: positive (negative) values indicate deficiency
7344
(excess) of electrons in the atom.
7346
\textit{Default value:} \texttt{.false.}
7349
\item[\textbf{WriteVoronoiPop}] (\textit{logical}):
7350
\index{WriteVoronoiPop@\textbf{WriteVoronoiPop}} %
7351
\index{Voronoi population analysis}\index{output!Voronoi analysis}
7353
If \texttt{.true.}, the program calculates and prints the Voronoi
7354
``net'' atomic populations on each atom in the system. For a
7355
definition of the Voronoi charges, see Bickelhaupt et al,
7356
Organometallics \textbf{15}, 2923 (1996) and Fonseca et al,
7357
J. Comp. Chem. \textbf{25}, 189 (2003). Voronoi charges are more
7358
reliable than Mulliken charges, specially for large basis sets. The
7359
number printed is the total net charge of the atom: the variation from
7360
the neutral charge, in units of $|e|$: positive (negative) values
7361
indicate deficiency (excess) of electrons in the atom.
7363
\textit{Default value:} \texttt{.false.}
7188
\begin{fdflogicalF}{WriteHirshfeldPop}
7189
\index{Hirshfeld population analysis}%
7190
\index{output!Hirshfeld analysis}
7192
If \fdftrue, the program calculates and prints the Hirshfeld ``net''
7193
atomic populations on each atom in the system. For a definition of
7194
the Hirshfeld charges, see Hirshfeld, Theo Chem Acta \textbf{44},
7195
129 (1977) and Fonseca et al, J. Comp. Chem. \textbf{25}, 189
7196
(2003). Hirshfeld charges are more reliable than Mulliken charges,
7197
specially for large basis sets. The number printed is the total net
7198
charge of the atom: the variation from the neutral charge, in units
7199
of $|e|$: positive (negative) values indicate deficiency (excess) of
7200
electrons in the atom.
7204
\begin{fdflogicalF}{WriteVoronoiPop}
7205
\index{Voronoi population analysis}%
7206
\index{output!Voronoi analysis}
7208
If \fdftrue, the program calculates and prints the Voronoi ``net''
7209
atomic populations on each atom in the system. For a definition of
7210
the Voronoi charges, see Bickelhaupt et al, Organometallics
7211
\textbf{15}, 2923 (1996) and Fonseca et al,
7212
J. Comp. Chem. \textbf{25}, 189 (2003). Voronoi charges are more
7213
reliable than Mulliken charges, specially for large basis sets. The
7214
number printed is the total net charge of the atom: the variation
7215
from the neutral charge, in units of $|e|$: positive (negative)
7216
values indicate deficiency (excess) of electrons in the atom.
7367
7220
The Hirshfeld and Voronoi populations (partial charges) are computed
7368
7221
by default only at the end of the program (i.e., for the final
7369
geometry, after self-consistency). The following options allow more control:
7375
\item[\textbf{PartialChargesAtEveryGeometry}] (\textit{logical}):
7376
\index{PartialChargesAtEveryGeometry@\textbf{PartialChargesAtEveryGeometry}}%
7222
geometry, after self-consistency). The following options allow more
7225
\begin{fdflogicalF}{PartialChargesAtEveryGeometry}
7377
7226
\index{Voronoi population analysis}%
7378
7227
\index{output!Voronoi analysis} %
7379
7228
\index{Hirshfeld population analysis} %
7380
7229
\index{output!Hirshfeld analysis}
7382
The Hirshfeld and Voronoi populations are computed after
7383
self-consistency is achieved, for all the geometry steps.
7385
\textit{Default value:} \texttt{.false.}
7387
\item[\textbf{PartialChargesAtEveryScfStep}] (\textit{logical}):
7388
\index{PartialChargesAtEveryScfStep@\textbf{PartialChargesAtEveryScfStep}}%
7231
The Hirshfeld and Voronoi populations are computed after
7232
self-consistency is achieved, for all the geometry steps.
7236
\begin{fdflogicalF}{PartialChargesAtEverySCFStep}
7389
7237
\index{Voronoi population analysis}%
7390
7238
\index{output!Voronoi analysis}%
7391
7239
\index{Hirshfeld population analysis}%
7392
7240
\index{output!Hirshfeld analysis}
7394
The Hirshfeld and Voronoi populations are computed for every
7395
step of the self-consistency process.
7397
\textit{Default value:} \texttt{.false.}
7242
The Hirshfeld and Voronoi populations are computed for every step of
7243
the self-consistency process.
7401
7247
\textbf{Performance note:}
7402
7248
The default behavior (computing at the end of the program) involves
7403
7249
an extra calculation of the charge density.
7407
7253
\subsubsection{Crystal-Orbital overlap and hamilton populations (COOP/COHP)}
7409
7255
\index{COOP/COHP curves}
7411
7257
These curves are quite useful to analyze the electronic structure to
7412
get insight about bonding characteristics. See the \texttt{Util/COOP}
7413
directory for more details. The \textbf{COOP.Write} option must be
7258
get insight about bonding characteristics. See the \program{Util/COOP}
7259
directory for more details. The \fdf{COOP.Write} option must be
7414
7260
activated to get the information needed.
7418
7263
\begin{itemize}
7419
\item Original COOP reference:
7420
Hughbanks, T.; Hoffmann, R., J. Am. Chem. Soc., 1983, 105, 3528.
7421
\item Original COHP reference: Dronskowski, R.; Blöchl, P. E., J. Phys. Chem., 1993, 97, 8617.
7422
\item A tutorial introduction: Dronskowski, R. Computational Chemistry of Solid State
7423
Materials; Wiley-VCH: Weinheim, 2005.
7424
\item Online material maintained by R. Dronskowski's group: \url{http://www.cohp.de/}
7265
Original COOP reference:
7266
Hughbanks, T.; Hoffmann, R., J. Am. Chem. Soc., 1983, 105, 3528.
7269
Original COHP reference: Dronskowski, R.; Blöchl, P. E., J. Phys. Chem., 1993, 97, 8617.
7272
A tutorial introduction: Dronskowski, R. Computational Chemistry of Solid State
7273
Materials; Wiley-VCH: Weinheim, 2005.
7276
Online material maintained by R. Dronskowski's group: \url{http://www.cohp.de/}
7431
\item[\textbf{COOP.Write}] (\textit{logical}):
7432
\index{COOP.Write@\textbf{COOP.Write}}
7433
\index{output!Information for COOP/COHP curves}
7434
Instructs the program to generate \texttt{SystemLabel.fullBZ}.WFSX (packed
7435
wavefunction file) and \texttt{SystemLabel}.HSX (H, S and X\_~{ij} file),
7436
to be processed by \texttt{Util/COOP/mprop} to generate COOP/COHP curves,
7437
(projected) densities of states, etc.
7439
The WFSX file is in a more compact form than the usual WFS, and the
7440
wavefunctions are output in single precision. The \texttt{Util/wfsx2wfs}
7441
program can be used to convert to the old format.
7442
The HSX file is in a more compact form than the usual HS, and the
7443
Hamiltonian, overlap matrix, and relative-positions array (which is
7444
always output, even for gamma-point only calculations) are in
7447
\textit{Default value:} \texttt{.false.}
7449
The user can narrow the energy-range used (and save some file space)
7450
by using the \textbf{WFS.Energy.Min} and \textbf{WFS.Energy.Max} options (both take an
7451
energy (with units) as extra argument), and/or the \textbf{WFS.Band.Min} and
7452
\textbf{WFS.Band.Max} options. Care should be taken to make sure that the
7453
actual values of the options make sense.
7455
Note that the band range options could also affect the output of
7456
wave-functions associated to bands (see section~\ref{sec:wf-bands}),
7457
and that the energy range options could also affect the output of
7458
user-selected wave-functions with the \textbf{WaveFuncKpoints} block (see
7459
section~\ref{sec:wf-output-user}).
7461
\item[\textbf{WFS.Energy.Min}] (\textit{physical energy}):
7462
\index{WFS.Energy.Min@\textbf{WFS.Energy.Min}}
7463
Specifies the lowest value of the energy (eigenvalue) of the wave-functions to be written to
7464
the file \texttt{SystemLabel.fullBZ}.WFSX for each k-point (all k-points
7465
in the BZ sampling are affected).
7467
\textit{Default value:} $-\infty$
7469
\item[\textbf{WFS.Energy.Max}] (\textit{physical energy}):
7470
\index{WFS.Energy.Max@\textbf{WFS.Energy.Max}}
7471
Specifies the highest value of the energy (eigenvalue) of the wave-functions to be written to
7472
the file \texttt{SystemLabel.fullBZ}.WFSX for each k-point (all k-points
7473
in the BZ sampling are affected).
7475
\textit{Default value:} $+\infty$
7477
\item[\textbf{WFS.Band.Min}] (\textit{integer}):
7478
\index{WFS.Band.Min@\textbf{WFS.Band.Min}}
7479
Specifies the lowest band index of the wave-functions to be written to
7480
the file (in this context) \texttt{SystemLabel.fullBZ}.WFSX for each k-point (all k-points
7481
in the BZ sampling are affected).
7483
\textit{Default value:} \texttt{1}
7485
\item[\textbf{WFS.Band.Max}] (\textit{integer}):
7486
\index{WFS.Band.Max@\textbf{WFS.Band.Max}}
7487
Specifies the highest band index of the wave-functions to be written to
7488
the file (in this context) \texttt{SystemLabel.fullBZ}.WFSX for each k-point (all k-points
7489
in the BZ sampling are affected).
7491
\textit{Default value:} \texttt{(the number of orbitals)}
7280
\begin{fdflogicalF}{COOP.Write}
7281
\index{output!Information for COOP/COHP curves}
7283
Instructs the program to generate \sysfile{fullBZ.WFSX} (packed
7284
wavefunction file) and \sysfile{HSX} (H, S and X\_~{ij} file),
7285
to be processed by \program{Util/COOP/mprop} to generate COOP/COHP curves,
7286
(projected) densities of states, etc.
7288
The \sysfile*{WFSX} file is in a more compact form than the usual
7289
\sysfile*{WFS}, and the wavefunctions are output in single
7290
precision. The \program{Util/wfsx2wfs} program can be used to
7291
convert to the old format. The HSX file is in a more compact form
7292
than the usual HS, and the Hamiltonian, overlap matrix, and
7293
relative-positions array (which is always output, even for
7294
gamma-point only calculations) are in single precision.
7296
The user can narrow the energy-range used (and save some file space)
7297
by using the \fdf{WFS.Energy.Min} and \fdf{WFS.Energy.Max} options
7298
(both take an energy (with units) as extra argument), and/or the
7299
\fdf{WFS.Band.Min} and \fdf{WFS.Band.Max} options. Care should be
7300
taken to make sure that the actual values of the options make sense.
7302
Note that the band range options could also affect the output of
7303
wave-functions associated to bands (see section~\ref{sec:wf-bands}),
7304
and that the energy range options could also affect the output of
7305
user-selected wave-functions with the \fdf{WaveFuncKPoints} block
7306
(see section~\ref{sec:wf-output-user}).
7311
\begin{fdfentry}{WFS.Energy.Min}[energy]<$-\infty$>
7313
Specifies the lowest value of the energy (eigenvalue) of the
7314
wave-functions to be written to the file
7315
\sysfile{fullBZ.WFSX} for each $k$-point (all $k$-points in
7316
the BZ sampling are affected).
7320
\begin{fdfentry}{WFS.Energy.Max}[energy]<$\infty$>
7322
Specifies the highest value of the energy (eigenvalue) of the
7323
wave-functions to be written to the file \sysfile{fullBZ.WFSX} for
7324
each $k$-point (all $k$-points in the BZ sampling are affected).
7497
7331
\subsection{Optical properties}
7503
\item[\textbf{OpticalCalculation}] (\textit{logical}):\index{Dielectric function,
7505
If specified, the imaginary part of the dielectric function
7506
will be calculated and stored in a file called \textit{Systemlabel}.\textbf{EPSIMG}.
7507
The calculation is performed using the simplest approach based on the
7508
dipolar transition matrix elements between different eigenfunctions
7509
of the self-consistent Hamiltonian. For molecules the calculation
7510
is performed using the position operator matrix elements, while
7511
for solids the calculation is carried out in the momentum space
7513
Corrections due to the non-locality of the pseudopotentials
7514
are introduced in the usual way.
7516
\textit{Default value:} \texttt{false}
7518
\item[\textbf{Optical.EnergyMinimum}] (\textit{real energy}):
7519
\index{Optical.EnergyMinimum@\textbf{Optical.EnergyMinimum}}
7520
This specifies the minimum of the energy range in which
7521
the frequency spectrum will be calculated.
7523
\textit{Default value:} 0 Ry.
7525
\item[\textbf{Optical.EnergyMaximum}] (\textit{real energy}):
7526
\index{Optical.EnergyMaximum@\textbf{Optical.EnergyMaximum}}
7527
This specifies the maximum of the energy range in which
7528
the frequency spectrum will be calculated.
7530
\textit{Default value:} 10 Ry.
7532
\item[\textbf{Optical.Broaden}] (\textit{real energy}):
7533
\index{Optical.Broaden@\textbf{Optical.Broaden}}
7534
If this is value is set then a Gaussian broadening will be
7535
applied to the frequency values.
7537
\textit{Default value:} 0 Ry.
7539
\item[\textbf{Optical.Scissor}] (\textit{real energy}):
7540
\index{Optical.Scissor@\textbf{Optical.Scissor}}
7541
Because of the tendency of DFT calculations to under estimate
7542
the band gap, a rigid shift of the unoccupied states, known as
7543
the scissor operator, can be added to correct the gap and
7544
thereby improve the calculated results. This shift is only
7545
applied to the optical calculation and no where else within
7548
\textit{Default value:} 0 Ry.
7550
\item[\textbf{Optical.NumberOfBands}] (\textit{integer}):
7551
\index{Optical.NumberOfBands@\textbf{Optical.NumberOfBands}}
7552
This option controls the number of bands that are included in
7553
the optical property calculation. Clearly this number must be
7554
larger than the number of occupied bands and less than or
7555
equal to the number of basis functions (which determines the
7556
number of unoccupied bands available). Note, while including
7557
all the bands may be the most accurate choice this will also
7558
be the most expensive!
7560
\textit{Default value:} All bands.
7562
\item[\textbf{Optical.Mesh}] (\textit{data block}):
7563
\index{Optical.Mesh@\textbf{Optical.Mesh}}
7564
This block contains 3 numbers that determine the mesh size
7565
used for the integration across the Brillouin zone. For
7571
%endblock Optical.Mesh
7574
The three values represent the number of mesh points in
7575
the direction of each reciprocal lattice vector.
7577
\textit{Default value:} Empty in general. For atoms
7578
or molecules a k-sampling of only one point is assumed.
7580
\item[\textbf{Optical.OffsetMesh}] (\textit{logical}):
7581
\index{Optical.OffsetMesh@\textbf{Optical.OffsetMesh}}
7582
If set to true, then the mesh is offset away from the
7583
gamma point for odd numbers of points.
7585
\textit{Default value:} \texttt{false}
7587
\item[\textbf{Optical.PolarizationType}] (\textit{string}):
7588
\index{Optical.PolarizationType@\textbf{Optical.PolarizationType}}
7589
This option has three possible values that represent the
7590
type of polarization to be used in the calculation. The options
7591
are \textbf{polarized}, which implies the application of an electric
7592
field in a given direction, \textbf{unpolarized}, which implies the
7593
propagation of light in a given direction, and \textbf{polycrystal}.
7594
In the case of the first two options a direction in space must
7595
be specified for the electric field or propagation using the
7596
\textit{Optical.Vector} data block.
7598
\textit{Default value:} \texttt{polycrystal}
7600
\item[\textbf{Optical.Vector}] (\textit{data block}):
7601
\index{Optical.Vector@\textbf{Optical.Vector}}
7602
This block contains 3 numbers that specify the vector direction
7603
for either the electric field or light propagation, for a polarized
7604
or unpolarized calculation, respectively. A typical block might look
7608
%block Optical.Vector
7610
%endblock Optical.Vector
7613
\textit{Default value:} Empty.
7333
\begin{fdflogicalF}{OpticalCalculation}
7334
\index{Dielectric function,optical absorption}
7336
If specified, the imaginary part of the dielectric function will be
7337
calculated and stored in a file called \sysfile{EPSIMG}. The
7338
calculation is performed using the simplest approach based on the
7339
dipolar transition matrix elements between different eigenfunctions
7340
of the self-consistent Hamiltonian. For molecules the calculation is
7341
performed using the position operator matrix elements, while for
7342
solids the calculation is carried out in the momentum space
7343
formulation. Corrections due to the non-locality of the
7344
pseudopotentials are introduced in the usual way.
7348
\begin{fdfentry}{Optical.Energy.Minimum}[energy]<$0\,\mathrm{Ry}$>
7350
This specifies the minimum of the energy range in which the
7351
frequency spectrum will be calculated.
7355
\begin{fdfentry}{Optical.Energy.Maximum}[energy]<$10\,\mathrm{Ry}$>
7357
This specifies the maximum of the energy range in which the
7358
frequency spectrum will be calculated.
7362
\begin{fdfentry}{Optical.Broaden}[energy]<$0\,\mathrm{Ry}$>
7364
If this is value is set then a Gaussian broadening will be applied
7365
to the frequency values.
7369
\begin{fdfentry}{Optical.Scissor}[energy]<$0\,\mathrm{Ry}$>
7371
Because of the tendency of DFT calculations to under estimate the
7372
band gap, a rigid shift of the unoccupied states, known as the
7373
scissor operator, can be added to correct the gap and thereby
7374
improve the calculated results. This shift is only applied to the
7375
optical calculation and no where else within the calculation.
7379
\begin{fdfentry}{Optical.NumberOfBands}[integer]<all bands>
7381
This option controls the number of bands that are included in the
7382
optical property calculation. Clearly this number must be larger
7383
than the number of occupied bands and less than or equal to the
7384
number of basis functions (which determines the number of unoccupied
7385
bands available). Note, while including all the bands may be the
7386
most accurate choice this will also be the most expensive!
7390
\begin{fdfentry}{Optical.Mesh}[block]
7392
This block contains 3 numbers that determine the mesh size used for
7393
the integration across the Brillouin zone. For example:
7397
%endblock Optical.Mesh
7399
The three values represent the number of mesh points in the
7400
direction of each reciprocal lattice vector.
7405
\begin{fdflogicalF}{Optical.OffsetMesh}
7407
If set to true, then the mesh is offset away from the gamma point
7408
for odd numbers of points.
7412
\begin{fdfentry}{Optical.PolarizationType}[string]<polycrystal>
7414
This option has three possible values that represent the type of
7415
polarization to be used in the calculation. The options are
7418
implies the application of an electric field in a given direction
7420
\option[unpolarized]%
7421
implies the propagation of light in a given direction
7423
\option[polycrystal]%
7424
In the case of the first two options a direction in space must be
7425
specified for the electric field or propagation using the
7426
\fdf{Optical.Vector} data block.
7432
\begin{fdfentry}{Optical.Vector}[block]
7434
This block contains 3 numbers that specify the vector direction for
7435
either the electric field or light propagation, for a polarized or
7436
unpolarized calculation, respectively. A typical block might look
7439
%block Optical.Vector
7441
%endblock Optical.Vector
7618
7449
\subsection{Macroscopic polarization}
7622
\item[\textbf{PolarizationGrids}] (\textit{data block}):
7623
\index{PolarizationGrids@\textbf{PolarizationGrids}}
7624
\index{bulk polarization}\index{Berry phase}
7625
If specified, the macroscopic polarization will be calculated using the
7626
geometric Berry phase approach (R.D. King-Smith, and D. Vanderbilt,
7627
PRB \textbf{47}, 1651 (1993)). In this method the electronic
7628
contribution to the macroscopic polarization, along a given direction,
7630
a discretized version of the formula
7633
P_{e,\parallel}={ifq_e \over 8\pi^3} \int_A d\textbf{k}_\perp
7452
\begin{fdfentry}{PolarizationGrids}[block]
7453
\index{bulk polarization}%
7456
If specified, the macroscopic polarization will be calculated using
7457
the geometric Berry phase approach (R.D. King-Smith, and
7458
D. Vanderbilt, PRB \textbf{47}, 1651 (1993)). In this method the
7459
electronic contribution to the macroscopic polarization, along a
7460
given direction, is calculated using a discretized version of the
7464
P_{e,\parallel}={ifq_e \over 8\pi^3} \int_A d\mathbf{k}_\perp
7634
7465
\sum_{n=1}^M \int_0^{|G_\parallel|} dk_{\parallel}
7635
\langle u_{\textbf{k} n} |{\delta \over \delta k_{\parallel}} |
7636
u_{\textbf{k} n} \rangle
7638
where $f$ is the occupation (2 for a non-magnetic system),
7639
$q_e$ the electron charge, $M$ is the number of occupied bands (the
7640
system \textbf{must} be an insulator), and $u_{\textbf{k} n}$ are
7641
the periodic Bloch functions. $\textbf{G}_\parallel$ is the shortest
7642
reciprocal vector along the chosen direction.
7466
\langle u_{\mathbf{k} n} |{\delta \over \delta k_{\parallel}} |
7467
u_{\mathbf{k} n} \rangle
7469
where $f$ is the occupation (2 for a non-magnetic system), $q_e$ the
7470
electron charge, $M$ is the number of occupied bands (the system
7471
\textbf{must} be an insulator), and $u_{\mathbf{k} n}$ are the
7472
periodic Bloch functions. $\mathbf{G}_\parallel$ is the shortest
7473
reciprocal vector along the chosen direction.
7644
As it can be seen in formula (\ref{pol_formula}), to compute each
7645
component of the polarization we must perform a surface integration
7646
of the result of a 1-D integral in the selected direction.
7647
The grids for the calculation along the direction of each of the
7648
three lattice vectors are specified in the block
7649
\textbf{PolarizationGrids}.
7475
As it can be seen in formula \eqref{pol_formula}, to compute each
7476
component of the polarization we must perform a surface integration
7477
of the result of a 1-D integral in the selected direction. The
7478
grids for the calculation along the direction of each of the three
7479
lattice vectors are specified in the block \fdf{PolarizationGrids}.
7651
7481
%block PolarizationGrids
7655
7485
%endblock PolarizationGrids
7658
All three grids must be specified, therefore a 3$\times$3 matrix of
7659
integer numbers must be given: the first row specifies the grid that
7660
will be used to calculate the polarization along the direction of the
7661
first lattice vector, the second row will be used for the calculation
7662
along the the direction of the second lattice vector, and the third
7663
row for the third lattice vector. The numbers in the diagonal of the
7664
matrix specifie the number of points to be used in the one dimensional
7665
line integrals along the different directions. The other numbers
7666
specifie the mesh used in the surface integrals. The last column
7667
specifies if the bidimensional grids are going to be diplaced from the
7668
origin or not, as in the Monkhorst-Pack algorithm (PRB \textbf{13}, 5188
7669
(1976)). This last column is optional. If the number of points in
7670
one of the grids is zero, the calculation will not be performed for
7671
this particular direction.
7673
For example, in the given example, for the computation in the
7674
direction of the first lattice vector, 15 points will be used for the
7675
line integrals, while a 3$\times$4 mesh will be used for the surface
7676
integration. This last grid will be displaced from the origin, so
7677
$\Gamma$ will not be included in the bidimensional integral. For the
7678
directions of the second and third lattice vectors, the number of
7679
points will be 20 and 2$\times$2, and 15 and 4$\times$4, respectively.
7681
It has to be stressed that the macroscopic polarization can only be
7682
meaningfully calculated using this approach for insulators.
7683
Therefore, the presence of an energy gap is necessary, and no band can
7684
cross the Fermi level. The program performs a simple check of this
7685
condition, just by counting the electrons in the unit cell ( the
7686
number must be even for a non-magnetic system, and the total spin
7687
polarization must have an integer value for spin polarized systems),
7688
however is the responsability of the user to check that the system
7689
under study is actually an insulator (for both spin components if spin
7692
The total macroscopic polarization, given in the output of the
7693
program, is the sum of the electronic contribution (calculated as the
7694
Berry phase of the valence bands), and the ionic contribution, which
7695
is simply defined as the sum of the atomic positions within the unit
7696
cell multiply by the ionic charges ($\sum_i^{N_a} Z_i \textbf{r}_i$). In
7697
the case of the magnetic systems, the bulk polarization for each spin
7698
component has been defined as
7700
\textbf{P}^\sigma = \textbf{P}_e^\sigma +
7701
{1 \over 2} \sum_i^{N_a} Z_i \textbf{r}_i
7703
$N_a$ is the number of atoms in the unit cell, and $\textbf{r}_i$ and
7705
are the positions and charges of the ions.
7707
It is also worth noting, that the macroscopic polarization given by
7708
formula (\ref{pol_formula}) is only defined modulo a ``quantum" of
7709
polarization (the bulk polarization per unit cell is only well defined
7710
modulo $fq_e$\textbf{R}, being \textbf{R} an arbitrary lattice
7711
vector). However, the experimentally observable quantities are
7712
associated to changes in the polarization induced by changes on the
7713
atomic positions (dynamical charges), strains (piezoelectric tensor),
7714
etc... The calculation of those changes, between different
7715
configurations of the solid, will be well defined as long as they are
7716
smaller than the ``quantum", i.e. the perturbations are small enough
7717
to create small changes in the polarization.
7719
\textit{Use:} Only compatible with \textbf{SolutionMethod} = diagon.\\
7720
\textit{Default value:} Empty. No calculation performed.
7722
\item[\textbf{BornCharge}] (\textit{logical}):\index{Born effective charges}
7723
\index{BornCharge@\textbf{BornCharge}}
7724
If true, the Born effective charge tensor is calculated for each atom
7725
by finite differences, by calculating the change in electric polarization
7726
(see \textbf{PolarizationGrids}) induced by the small displacements generated
7727
for the force constants calculation (see \textbf{MD.TypeOfRun} = \texttt{FC}):
7728
\begin{equation}\label {eq:effective_charge}
7729
Z^*_{i,\alpha,\beta}=\frac{\Omega_0}{e} \left. {\frac{\partial{P_\alpha}}
7730
{\partial{u_{i,\beta}}}}\right|_{q=0}
7732
where e is the charge of an electron and $\Omega_0$ is the unit cell volume.
7734
To calculate the Born charges it is necessary to specify both the Born
7735
charge flag and the mesh used to calculate the polarization, for example:
7737
%block PolarizationGrids
7741
%endblock PolarizationGrids
7745
The Born effective charge matrix is then written to the file
7746
\textit{SystemLabel}.\texttt{BC}.
7748
The method by which the polarization is calculated may introduce an arbitrary
7749
phase (polarization quantum), which in general is far larger than the change
7750
in polarization which results from the atomic displacement. It is removed
7751
during the calculation of the Born effective charge tensor.
7753
The Born effective charges allow the calculation of LO-TO splittings and
7754
infrared activities. The version of the Vibra utility code in which these
7755
magnitudes are calculated is not yet distributed with \siesta, but can be
7756
obtained form Tom Archer (archert@tcd.ie).
7758
\textit{Use:} Only used if \textbf{MD.TypeOfRun} is \texttt{FC}.
7760
\textit{Default value:} \texttt{false}
7488
All three grids must be specified, therefore a $3\times3$ matrix of
7489
integer numbers must be given: the first row specifies the grid that
7490
will be used to calculate the polarization along the direction of
7491
the first lattice vector, the second row will be used for the
7492
calculation along the the direction of the second lattice vector,
7493
and the third row for the third lattice vector. The numbers in the
7494
diagonal of the matrix specifie the number of points to be used in
7495
the one dimensional line integrals along the different
7496
directions. The other numbers specifie the mesh used in the surface
7497
integrals. The last column specifies if the bidimensional grids are
7498
going to be diplaced from the origin or not, as in the
7499
Monkhorst-Pack algorithm (PRB \textbf{13}, 5188 (1976)). This last
7500
column is optional. If the number of points in one of the grids is
7501
zero, the calculation will not be performed for this particular
7504
For example, in the given example, for the computation in the
7505
direction of the first lattice vector, 15 points will be used for
7506
the line integrals, while a $3\times4$ mesh will be used for the
7507
surface integration. This last grid will be displaced from the
7508
origin, so $\Gamma$ will not be included in the bidimensional
7509
integral. For the directions of the second and third lattice
7510
vectors, the number of points will be $20$ and $2\times2$, and $15$
7511
and $4\times4$, respectively.
7513
It has to be stressed that the macroscopic polarization can only be
7514
meaningfully calculated using this approach for insulators.
7515
Therefore, the presence of an energy gap is necessary, and no band
7516
can cross the Fermi level. The program performs a simple check of
7517
this condition, just by counting the electrons in the unit cell (
7518
the number must be even for a non-magnetic system, and the total
7519
spin polarization must have an integer value for spin polarized
7520
systems), however is the responsability of the user to check that
7521
the system under study is actually an insulator (for both spin
7522
components if spin polarized).
7524
The total macroscopic polarization, given in the output of the
7525
program, is the sum of the electronic contribution (calculated as
7526
the Berry phase of the valence bands), and the ionic contribution,
7527
which is simply defined as the sum of the atomic positions within
7528
the unit cell multiply by the ionic charges
7529
($\sum_i^{N_a} Z_i \mathbf{r}_i$). In the case of the magnetic
7530
systems, the bulk polarization for each spin component has been
7533
\mathbf{P}^\sigma = \mathbf{P}_e^\sigma +
7534
{1 \over 2} \sum_i^{N_a} Z_i \mathbf{r}_i
7536
$N_a$ is the number of atoms in the unit cell, and $\mathbf{r}_i$
7537
and $Z_i$ are the positions and charges of the ions.
7539
It is also worth noting, that the macroscopic polarization given by
7540
formula \eqref{pol_formula} is only defined modulo a ``quantum'' of
7541
polarization (the bulk polarization per unit cell is only well
7542
defined modulo $fq_e\mathbf{R}$, being $\mathbf{R}$ an arbitrary
7543
lattice vector). However, the experimentally observable quantities
7544
are associated to changes in the polarization induced by changes on
7545
the atomic positions (dynamical charges), strains (piezoelectric
7546
tensor), etc... The calculation of those changes, between different
7547
configurations of the solid, will be well defined as long as they
7548
are smaller than the ``quantum'', i.e. the perturbations are small
7549
enough to create small changes in the polarization.
7553
\begin{fdflogicalF}{BornCharge}
7554
\index{Born effective charges}
7556
If true, the Born effective charge tensor is calculated for each
7557
atom by finite differences, by calculating the change in electric
7558
polarization (see \fdf{PolarizationGrids}) induced by the small
7559
displacements generated for the force constants calculation (see
7560
\fdf{MD.TypeOfRun} \fdf*{FC}):
7562
\label{eq:effective_charge}
7563
Z^*_{i,\alpha,\beta}=\frac{\Omega_0}{e} \left. {\frac{\partial{P_\alpha}}
7564
{\partial{u_{i,\beta}}}}\right|_{q=0}
7566
where e is the charge of an electron and $\Omega_0$ is the unit cell
7569
To calculate the Born charges it is necessary to specify both the
7570
Born charge flag and the mesh used to calculate the polarization,
7573
%block PolarizationGrids
7577
%endblock PolarizationGrids
7581
The Born effective charge matrix is then written to the file
7584
The method by which the polarization is calculated may introduce an
7585
arbitrary phase (polarization quantum), which in general is far
7586
larger than the change in polarization which results from the atomic
7587
displacement. It is removed during the calculation of the Born
7588
effective charge tensor.
7590
The Born effective charges allow the calculation of LO-TO splittings
7591
and infrared activities. The version of the Vibra utility code in
7592
which these magnitudes are calculated is not yet distributed with
7593
\siesta, but can be obtained form Tom Archer (archert@tcd.ie).
7764
7599
\subsection[Maximally Localized Wannier Functions]%
7765
7600
{Maximally Localized Wannier Functions. \\
7766
7601
Interface with the \textsc{wannier90} code}
7768
\textsc{wannier90} (http://www.wannier.org) is a code to generate
7603
\program{wannier90} (http://www.wannier.org) is a code to generate
7769
7604
maximally localized wannier functions according to the original
7770
7605
Marzari and Vanderbilt recipe.
7772
7607
It is strongly recommended to read the original papers on which this
7773
method is based and the documentation of the \textsc{wannier90} code.
7608
method is based and the documentation of \program{wannier90} code.
7774
7609
Here we shall focus only on those internal \siesta\ variables
7775
7610
required to produce the files that will be processed
7776
by \textsc{wannier90}.
7611
by \program{wannier90}.
7778
7613
A complete list of examples and tests (including molecules, metals,
7779
7614
semiconductors, insulators, magnetic systems, plotting of Fermi surfaces
7786
7621
used and other possibly non-reproducible details of the
7787
7622
calculation. In what follows it is essential to maintain
7788
7623
consistency in the handling of the overlap and Bloch-funcion
7789
files produced and fed to \textsc{wannier90}.
7792
\item[\textbf{Siesta2Wannier90.WriteMmn}] (\textit{logical}):
7793
\index{Siesta2Wannier90.WriteMmn@\textbf{Siesta2Wannier90.WriteMmn}}
7795
This flag determines whether the overlaps between the periodic part
7796
of the Bloch states at neighbour k-points are computed and dumped into a file
7797
in the format required by \textsc{wannier90}.
7798
These overlaps are defined in Eq. (27) in the paper by
7799
N. Marzari \textit{et al.}, Review of Modern Physics \textbf{84}, 1419 (2012),
7800
or Eq. (1.7) of the Wannier90 User Guide, Version 2.0.1.
7802
The k-points for which the overlaps will be computed are read from a
7803
\texttt{.nnkp} file produced by \textsc{wannier90}. It is strongly
7804
recommended for the user to read the corresponding user guide.
7806
The overlap matrices are written in a file with extension \texttt{.mmn}.
7808
\textit{Default value:} \texttt{.false.}
7810
\item[\textbf{Siesta2Wannier90.WriteAmn}] (\textit{logical}):
7811
\index{Siesta2Wannier90.WriteAmn@\textbf{Siesta2Wannier90.WriteAmn}}
7813
This flag determines whether the overlaps between Bloch states
7814
and trial localized orbitals are computed and dumped into a file in the
7815
format required by \textsc{wannier90}.
7816
These projections are defined in Eq. (16) in the paper by
7817
N. Marzari \textit{et al.}, Review of Modern Physics \textbf{84}, 1419 (2012),
7818
or Eq. (1.8) of the Wannier90 User Guide, Version 2.0.1.
7820
The localized trial functions to use are taken from the \texttt{.nnkp}
7821
file produced by \textsc{wannier90}. It is strongly recommended for the
7822
user to read the corresponding user guide.
7824
The overlap matrices are written in a file with extension \texttt{.amn}.
7826
\textit{Default value:} \texttt{.false.}
7828
\item[\textbf{Siesta2Wannier90.WriteEig}] (\textit{logical}):
7829
\index{Siesta2Wannier90.WriteEig@\textbf{Siesta2Wannier90.WriteEig}}
7831
Flag that determines whether the Kohn-Sham eigenvalues
7832
(in eV) at each point in the Monkhorst-Pack mesh required by
7833
\program{wannier90} are written to file.
7834
This file is mandatory in \program{wannier90} if any of disentanglement,
7835
plot\_bands, plot\_fermi\_surface or hr\_plot options are set to true in
7836
the \program{wannier90} input file.
7838
The eigenvalues are written in a file with extension \sysfile{eigW}.
7839
This extension is chosen to avoid name clashes
7840
with \siesta's standard eigenvalue file in case-insensitive filesystems.
7842
\textit{Default value:} \texttt{.false.}
7844
\item[\textbf{Siesta2Wannier90.WriteUnk}] (\textit{logical}):
7845
\index{Siesta2Wannier90.WriteUnk@\textbf{Siesta2Wannier90.WriteUnk}}
7847
Produces \texttt{UNKXXXXX.Y} files which contain the periodic part of a
7848
Bloch function in the unit cell on a grid given by global unk\_nx,
7849
unk\_ny, unk\_nz variables. The name of the output files is assumed
7850
to have the previous form, where the \texttt{XXXXXX} refer to the k-point
7851
index (from 00001 to the total number of k-points considered), and
7852
the \texttt{Y} refers to the spin component (1 or 2)
7854
The periodic part of the Bloch functions is defined by
7856
u_{n \vec{k}} (\vec{r}) =
7857
\sum_{\vec{R} \mu} c_{n \mu}(\vec{k})
7858
e^{i \vec{k} \cdot ( \vec{r}_{\mu} + \vec{R} - \vec{r} )}
7859
\phi_{\mu} (\vec{r} - \vec{r}_{\mu} - \vec{R} ) ,
7862
\noindent where $\phi_{\mu} (\vec{r} - \vec{r}_{\mu} - \vec{R} )$
7863
is a basis set atomic orbital centered on atom $\mu$ in
7864
the unit cell $\vec{R}$, and $c_{n \mu}(\vec{k})$ are the coefficients
7865
of the wave function. The latter must be identical to the ones used
7866
for wannierization in $M_{mn}$. (See the above comment about
7869
\textit{Default value:} \texttt{.false.}
7871
\item[\textbf{Siesta2Wannier90.UnkGrid1}] (\textit{integer}):
7872
\index{Siesta2Wannier90.UnkGrid1@\textbf{Siesta2Wannier90.UnkGrid1}}
7874
Number of points along the first lattice vector in the grid where
7875
the periodic part of the wave functions will be plotted.
7877
\textit{Default value:}{number of points in the \siesta\ mesh along this lattice vector}
7879
\item[\textbf{Siesta2Wannier90.UnkGrid2}] (\textit{integer}):
7880
\index{Siesta2Wannier90.UnkGrid2@\textbf{Siesta2Wannier90.UnkGrid2}}
7882
Number of points along the second lattice vector in the grid where
7883
the periodic part of the wave functions will be plotted.
7885
\textit{Default value:}{number of points in the \siesta\ mesh along this lattice vector}
7887
\item[\textbf{Siesta2Wannier90.UnkGrid3}] (\textit{integer}):
7888
\index{Siesta2Wannier90.UnkGrid3@\textbf{Siesta2Wannier90.UnkGrid3}}
7890
Number of points along the third lattice vector in the grid where
7891
the periodic part of the wave functions will be plotted.
7893
\textit{Default value:}{number of points in the \siesta\ mesh along this lattice vector}
7895
\item[\textbf{Siesta2Wannier90.UnkGridBinary}] (\textit{logical}):
7896
\index{Siesta2Wannier90.UnkGridBinary@\textbf{Siesta2Wannier90.UnkGridBinary}}
7898
Flag that determines whether the periodic part of the wave function in the
7899
real space grid is written in binary format (default) or in ASCII format.
7901
\textit{Default value:} \texttt{.true.}
7903
\item[\textbf{Siesta2Wannier90.NumberOfBands}] (\textit{integer}):
7904
\index{Siesta2Wannier90.NumberOfBands@\textbf{Siesta2Wannier90.NumberOfBands}}
7906
In spin unpolarized calculations, number of bands that will be
7907
initially considered by \siesta\ to generate the information
7908
required by \textsc{wannier90}. Note that it should be at least as large as the
7909
index of the highest-lying band in the \textsc{wannier90} post-processing. For
7910
example, if the wannierization is going to involve bands 3 to 5, the
7911
\siesta\ number of bands should be at least 5. Bands 1 and 2
7912
should appear in a ``excluded'' list.
7914
\textit{Default value:} {Number of occupied bands. To avoid problems
7915
of interpretation, it is strongly advised to avoid relying on the
7916
defaults and to specify the number of bands explicitly.}
7919
\item[\textbf{Siesta2Wannier90.NumberOfBandsUp}] (\textit{integer}):
7920
\index{Siesta2Wannier90.NumberOfBandsUp@\textbf{Siesta2Wannier90.NumberOfBandsUp}}
7922
In spin-polarized calculations, number of bands with spin up that will be
7923
initially considered by \siesta\ to generate the information
7924
required by \textsc{wannier90}. (See above for details.)
7926
\textit{Default value:} {Number of occupied bands with spin up. To
7927
avoid problems of interpretation, it is strongly advised to avoid
7928
relying on the defaults and to specify the number of bands
7931
\item[\textbf{Siesta2Wannier90.NumberOfBandsDown}] (\textit{integer}):
7932
\index{Siesta2Wannier90.NumberOfBandsDown@\textbf{Siesta2Wannier90.NumberOfBandsDown}}
7934
In spin-polarized calculations, number of bands with spin down that will be
7935
initially considered by \siesta\ to generate the information
7936
required by \textsc{wannier90}. (See above for details.)
7939
\textit{Default value:} {Number of occupied bands with spin down. To
7940
avoid problems of interpretation, it is strongly advised to avoid
7941
relying on the defaults and to specify the number of bands
7624
files produced and fed to \program{wannier90}.
7627
\begin{fdflogicalF}{Siesta2Wannier90.WriteMmn}
7629
This flag determines whether the overlaps between the periodic part
7630
of the Bloch states at neighbour k-points are computed and dumped
7631
into a file in the format required by \program{wannier90}. These
7632
overlaps are defined in Eq. (27) in the paper by N. Marzari
7633
\textit{et al.}, Review of Modern Physics \textbf{84}, 1419 (2012),
7634
or Eq. (1.7) of the Wannier90 User Guide, Version 2.0.1.
7636
The k-points for which the overlaps will be computed are read from a
7637
\sysfile*{nnkp} file produced by \program{wannier90}. It is strongly
7638
recommended for the user to read the corresponding user guide.
7640
The overlap matrices are written in a file with extension
7645
\begin{fdflogicalF}{Siesta2Wannier90.WriteAmn}
7647
This flag determines whether the overlaps between Bloch states and
7648
trial localized orbitals are computed and dumped into a file in the
7649
format required by \program{wannier90}. These projections are
7650
defined in Eq. (16) in the paper by N. Marzari \textit{et al.},
7651
Review of Modern Physics \textbf{84}, 1419 (2012), or Eq. (1.8) of
7652
the Wannier90 User Guide, Version 2.0.1.
7654
The localized trial functions to use are taken from the
7655
\sysfile*{nnkp} file produced by \program{wannier90}. It is strongly
7656
recommended for the user to read the corresponding user guide.
7658
The overlap matrices are written in a file with extension
7663
\begin{fdflogicalF}{Siesta2Wannier90.WriteEig}
7665
Flag that determines whether the Kohn-Sham eigenvalues (in eV) at
7666
each point in the Monkhorst-Pack mesh required by
7667
\program{wannier90} are written to file. This file is mandatory in
7668
\program{wannier90} if any of disentanglement, plot\_bands,
7669
plot\_fermi\_surface or hr\_plot options are set to true in the
7670
\program{wannier90} input file.
7672
The eigenvalues are written in a file with extension \sysfile*{eigW}.
7673
This extension is chosen to avoid name clashes with \siesta's
7674
standard eigenvalue file in case-insensitive filesystems.
7678
\begin{fdflogicalF}{Siesta2Wannier90.WriteUnk}
7680
Produces \file{UNKXXXXX.Y} files which contain the periodic part
7681
of a Bloch function in the unit cell on a grid given by global
7682
unk\_nx, unk\_ny, unk\_nz variables. The name of the output files
7683
is assumed to have the previous form, where the \texttt{XXXXXX}
7684
refer to the k-point index (from 00001 to the total number of
7685
k-points considered), and the \texttt{Y} refers to the spin
7688
The periodic part of the Bloch functions is defined by
7690
u_{n \vec{k}} (\vec{r}) =
7691
\sum_{\vec{R} \mu} c_{n \mu}(\vec{k})
7692
e^{i \vec{k} \cdot ( \vec{r}_{\mu} + \vec{R} - \vec{r} )}
7693
\phi_{\mu} (\vec{r} - \vec{r}_{\mu} - \vec{R} ) ,
7695
where $\phi_{\mu} (\vec{r} - \vec{r}_{\mu} - \vec{R} )$ is a basis
7696
set atomic orbital centered on atom $\mu$ in the unit cell
7697
$\vec{R}$, and $c_{n \mu}(\vec{k})$ are the coefficients of the wave
7698
function. The latter must be identical to the ones used for
7699
wannierization in $M_{mn}$. (See the above comment about arbitrary
7704
\begin{fdfentry}{Siesta2Wannier90.UnkGrid1}[integer]<\nonvalue{mesh
7707
Number of points along the first lattice vector in the grid where
7708
the periodic part of the wave functions will be plotted.
7712
\begin{fdfentry}{Siesta2Wannier90.UnkGrid2}[integer]<\nonvalue{mesh
7715
Number of points along the second lattice vector in the grid where
7716
the periodic part of the wave functions will be plotted.
7720
\begin{fdfentry}{Siesta2Wannier90.UnkGrid3}[integer]<\nonvalue{mesh
7723
Number of points along the third lattice vector in the grid where
7724
the periodic part of the wave functions will be plotted.
7729
\begin{fdflogicalT}{Siesta2Wannier90.UnkGridBinary}
7731
Flag that determines whether the periodic part of the wave function
7732
in the real space grid is written in binary format (default) or in
7737
\begin{fdfentry}{Siesta2Wannier90.NumberOfBands}[integer]<occupied bands>
7739
In spin unpolarized calculations, number of bands that will be
7740
initially considered by \siesta\ to generate the information
7741
required by \program{wannier90}. Note that it should be at least as
7742
large as the index of the highest-lying band in the
7743
\program{wannier90} post-processing. For example, if the
7744
wannierization is going to involve bands 3 to 5, the \siesta\ number
7745
of bands should be at least 5. Bands 1 and 2 should appear in a
7748
\note you are highly encouraged to explicitly specify the number of
7753
\begin{fdfentry}{Siesta2Wannier90.NumberOfBandsUp}[integer]<\fdfvalue{Siesta2Wannier90.NumberOfBands}>
7755
In spin-polarized calculations, number of bands with spin up that
7756
will be initially considered by \siesta\ to generate the information
7757
required by \program{wannier90}.
7761
\begin{fdfentry}{Siesta2Wannier90.NumberOfBandsDown}[integer]<\fdfvalue{Siesta2Wannier90.NumberOfBands}>
7763
In spin-polarized calculations, number of bands with spin down that
7764
will be initially considered by \siesta\ to generate the information
7765
required by \program{wannier90}.
7947
7771
\subsection{Systems with net charge or dipole, and electric fields}
7949
\begin{fdfentry}{NetCharge}[real]<$0.$>%
7773
\begin{fdfentry}{NetCharge}[real]<$0$>%
7950
7774
\index{Charge of the system}
7951
7775
\index{Doping}%
7952
7776
\index{SCF!Doping}%
7961
7785
is not advised to do charged systems other than atoms and molecules
7962
7786
in SC, FCC or BCC cells, unless you know what you are doing.
7964
\textit{Use:} For example, the F$^-$ ion would have \fdf*{NetCharge
7965
-1}, and the Na$^+$ ion would have \fdf*{NetCharge 1}.
7788
\textit{Use:} For example, the F$^-$ ion would have \fdf{NetCharge}
7789
\fdf*{-1} , and the Na$^+$ ion would have \fdf{NetCharge} \fdf*{1}.
7966
7790
Fractional charges can also be used.
7973
\item[\textbf{SimulateDoping}] (\textit{boolean}):
7974
\index{SimulateDoping@\textbf{SimulateDoping}}\index{Slabs with net charge}
7976
This option instructs the program to add a background charge density
7977
to simulate doping. The new ``doping'' routine calculates the net
7978
charge of the system, and adds a compensating background charge that
7979
makes the system neutral. This background charge is constant at points
7980
of the mesh near the atoms, and zero at points far from the atoms.
7981
This simulates situations like doped slabs, where the extra electrons
7982
(holes) are compensated by oposite charges at the material (the
7983
ionized dopant impurities), but not at the vacuum. This serves to
7984
simulate properly doped systems in which there are large portions of
7985
vacuum, such as doped slabs.
7987
(See \texttt{Tests/sic-slab})
7989
\textit{Default value:} \texttt{.false.}
7991
\item[\textbf{ExternalElectricField}] (\textit{data block}):
7992
\index{ExternalElectricField@\textbf{ExternalElectricField}}
7993
It specifies an external electric field for molecules, chains and slabs.
7994
The electric field should be orthogonal to `bulk directions', like
7995
those parallel to a slab (bulk electric fields, like in
7996
dielectrics or ferroelectrics, are not allowed). If it is not, an
7997
error message is issued and the components of the field in bulk
7998
directions are suppressed automatically. The input is a
7999
vector in Cartesian coordinates, in the specified units. Example:
7795
\begin{fdflogicalF}{SimulateDoping}
7796
\index{Slabs with net charge}
7798
This option instructs the program to add a background charge density
7799
to simulate doping. The new ``doping'' routine calculates the net
7800
charge of the system, and adds a compensating background charge that
7801
makes the system neutral. This background charge is constant at
7802
points of the mesh near the atoms, and zero at points far from the
7803
atoms. This simulates situations like doped slabs, where the extra
7804
electrons (holes) are compensated by oposite charges at the material
7805
(the ionized dopant impurities), but not at the vacuum. This serves
7806
to simulate properly doped systems in which there are large portions
7807
of vacuum, such as doped slabs.
7809
(See \texttt{Tests/sic-slab})
7813
\begin{fdfentry}{ExternalElectricField}[block]
7815
It specifies an external electric field for molecules, chains and
7816
slabs. The electric field should be orthogonal to `bulk
7817
directions', like those parallel to a slab (bulk electric fields,
7818
like in dielectrics or ferroelectrics, are not allowed). If it is
7819
not, an error message is issued and the components of the field in
7820
bulk directions are suppressed automatically. The input is a vector
7821
in Cartesian coordinates, in the specified units. Example:
8002
7823
%block ExternalElectricField
8003
7824
0.000 0.000 0.500 V/Ang
8004
7825
%endblock ExternalElectricField
8007
Starting with version 4.0, applying an electric field perpendicular to
8008
a slab will by default enable the slab dipole
8009
correction\index{SlabDipoleCorrection@\textbf{SlabDipoleCorrection}}
8010
option. To reproduce older calculations, set this correction option
8011
explicitly to \texttt{.false.} in the input file.
8013
\textit{Default value:} zero field
8015
\item[\textbf{SlabDipoleCorrection}] (\textit{boolean}):
8016
\index{SlabDipoleCorrection@\textbf{SlabDipoleCorrection}}
8018
If \texttt{true}, \siesta\ calculates the electric field required to
8019
compensate the dipole of the system at every iteration of the
8020
self-consistent cycle. The potential added to the grid corresponds to
8021
that of a dipole layer at the middle of the vacuum layer. For slabs,
8022
this exactly compensates the electric field at the vacuum created by
8023
the dipole moment of the system, thus allowing to treat asymmetric
8024
slabs (including systems with an adsorbate on one surface) and compute
8025
properties such as the work funcion of each of the surfaces.
8027
NOTE: If the program is fed a starting density matrix from an
8028
uncorrected calculation (i.e., with an exagerated dipole), the first
8029
iteration might use a compensating field that is too big, with the
8030
risk of taking the system out of the convergence basin. In that case,
8031
it is advisable to use the \texttt{DM.MixSCF1}\index{DM.MixSCF1@\textbf{DM.MixSCF1}}\index{Slab dipole correction} option to request a mix of the
8032
input and output density matrices after that first iteration.
8034
(See \texttt{Tests/sic-slab})
8036
This will default to \texttt{true} if an external field is applied to a slab
8037
calculation, otherwise it will default to \texttt{false}.
8039
\textit{Default value:} \texttt{false}
7828
Starting with version 4.0, applying an electric field perpendicular
7829
to a slab will by default enable the slab dipole correction, see
7830
\fdf{SlabDipoleCorrection}. To reproduce older calculations, set
7831
this correction option explicitly to \fdffalse\ in the input file.
7835
\begin{fdflogicalF}{SlabDipoleCorrection}
7837
If \fdftrue, \siesta\ calculates the electric field required to
7838
compensate the dipole of the system at every iteration of the
7839
self-consistent cycle. The potential added to the grid corresponds
7840
to that of a dipole layer at the middle of the vacuum layer. For
7841
slabs, this exactly compensates the electric field at the vacuum
7842
created by the dipole moment of the system, thus allowing to treat
7843
asymmetric slabs (including systems with an adsorbate on one
7844
surface) and compute properties such as the work funcion of each of
7847
NOTE: If the program is fed a starting density matrix from an
7848
uncorrected calculation (i.e., with an exagerated dipole), the first
7849
iteration might use a compensating field that is too big, with the
7850
risk of taking the system out of the convergence basin. In that
7851
case, it is advisable to use the
7852
\fdf{SCF.Mix!First}\index{Slab dipole correction} option to request a mix of the input and
7853
output density matrices after that first iteration.
7855
(See \texttt{Tests/sic-slab})
7857
This will default to \fdftrue\ if an external field is applied to a
7858
slab calculation, otherwise it will default to \fdffalse.
8044
7863
\begin{fdfentry}{Geometry!Hartree}[block]%
8213
8033
\program{Util/Contrib/APostnikov}. See also \program{Util/Denchar} for
8214
8034
an alternative way to plot the charge density (and wavefunctions).
8220
\item[\textbf{SaveRho}] (\textit{logical}): \index{SaveRho@\textbf{SaveRho}}\index{output!charge density} Instructs to write the
8221
valence pseudocharge density at the mesh used by DHSCF, in file \sysfile{RHO}.
8223
\textit{Use:} File \texttt{SystemLabel}.RHO is only written, not read, by siesta.
8224
This file can be read by routine IORHO, which may be used by other
8225
application programs.
8227
If netCDF support is compiled in, the file \texttt{Rho.grid.nc} is produced.
8229
\textit{Default value:} \texttt{.false.}
8232
\item[\textbf{SaveDeltaRho}] (\textit{logical}):
8233
\index{SaveDeltaRho@\textbf{SaveDeltaRho}}\index{output!$\delta \rho(\vec r)$}
8234
Instructs to write $\delta \rho(\vec r) = \rho(\vec r) - \rho_{atm}(\vec r)$,
8235
i.e., the valence pseudocharge density minus the sum of atomic valence
8236
pseudocharge densities. It is done for the mesh points used by DHSCF and it
8237
comes in file \texttt{SystemLabel}.DRHO. This file can be read by routine IORHO,
8238
which may be used by an application program in later versions.
8240
\textit{Use:} File \texttt{SystemLabel}.DRHO is only written, not read, by siesta.
8242
If netCDF support is compiled in, the file \texttt{DeltaRho.grid.nc} is produced.
8244
\textit{Default value:} \texttt{.false.}
8247
\item[\textbf{SaveRhoXC}] (\textit{logical}): \index{SaveRhoXC@\textbf{SaveRhoXC}}\index{output!charge density} Instructs to write the
8248
valence pseudocharge density at the mesh, including the nonlocal
8249
core corrections used to calculate the exchange-correlation energy,
8250
in file \sysfile{RHOXC}.
8252
\textit{Use:} File \texttt{SystemLabel}.RHOXC is only written, not read, by siesta.
8254
If netCDF support is compiled in, the file \texttt{RhoXC.grid.nc} is produced.
8256
\textit{Default value:} \texttt{.false.}
8259
\item[\textbf{SaveElectrostaticPotential}] (\textit{logical}):
8260
\index{SaveElectrostaticPotential@\textbf{SaveElectrostaticPotential}}
8261
\index{output!electrostatic potential}
8262
Instructs to write the total electrostatic potential, defined as the
8263
sum of the hartree potential plus the local pseudopotential, at the
8265
in file \texttt{SystemLabel}.VH. This file can be read by routine IORHO,
8266
which may be used by an application program in later versions.
8268
\textit{Use:} File \texttt{SystemLabel}.VH is only written, not read, by siesta.
8270
If netCDF support is compiled in, the file \file{ElectrostaticPotential.grid.nc} is produced.
8272
\textit{Default value:} \texttt{.false.}
8274
\item[\textbf{SaveNeutralAtomPotential}] (\textit{logical}):
8275
\index{SaveNeutralAtomPotential@\textbf{SaveNeutralAtomPotential}}
8276
\index{output!electrostatic potential}
8277
Instructs to write the neutral-atom potential, defined as the
8278
sum of the hartree potential of a ``pseudo atomic valence charge''
8279
plus the local pseudopotential, at the mesh used by DHSCF,
8280
in file \texttt{SystemLabel}.VNA. It is written at the start of the
8281
self-consistency cycle, as this potential does not change.
8283
\textit{Use:} File \texttt{SystemLabel}.VNA is only written, not read, by siesta.
8285
If netCDF support is compiled in, the file \texttt{Vna.grid.nc} is produced.
8287
\textit{Default value:} \texttt{.false.}
8290
\item[\textbf{SaveTotalPotential}] (\textit{logical}):
8291
\index{SaveTotalPotential@\textbf{SaveTotalPotential}}
8292
\index{output!total potential}
8293
Instructs to write the valence total effective local potential
8294
(local pseudopotential + Hartree + Vxc), at the
8296
in file \texttt{SystemLabel}.VT. This file can be read by routine IORHO,
8297
which may be used by an application program in later versions.
8299
\textit{Use:} File \texttt{SystemLabel}.VT is only written, not read, by siesta.
8301
If netCDF support is compiled in, the file \texttt{TotalPotential.grid.nc} is produced.
8303
\textit{Default value:} \texttt{.false.}
8305
\textit{Side effect:} the vacuum level, defined as the effective
8306
potential at grid points with zero density, is printed in the standard
8307
output whenever such points exist (molecules, slabs) and either
8308
SaveElectrostaticPotential or SaveTotalPotential are true.
8309
In a symetric (nonpolar) slab, the work function can be computed as
8310
the difference between the vacuum level and the Fermi energy.
8313
\item[\textbf{SaveIonicCharge}] (\textit{logical}):
8314
\index{SaveIonicCharge@\textbf{SaveIonicCharge}}
8315
\index{output!ionic charge}
8316
Instructs to write the soft diffuse ionic charge at the
8318
in file \texttt{SystemLabel}.IOCH. This file can be read by routine IORHO,
8319
which may be used by an application program in later versions.
8320
Remember that, within the \siesta\ sign convention, the electron charge
8321
density is positive and the ionic charge density is negative.
8324
\textit{Use:} File \texttt{SystemLabel}.IOCH is only written, not read, by siesta.
8326
If netCDF support is compiled in, the file \texttt{Chlocal.grid.nc} is produced.
8328
\textit{Default value:} \texttt{.false.}
8330
\item[\textbf{SaveTotalCharge}] (\textit{logical}):
8331
\index{SaveTotalCharge@\textbf{SaveTotalCharge}}
8332
\index{output!total charge}
8333
Instructs to write the total charge density (ionic+electronic) at the
8335
in file \texttt{SystemLabel}.TOCH. This file can be read by routine IORHO,
8336
which may be used by an application program in later versions.
8337
Remember that, within the \siesta\ sign convention, the electron charge
8338
density is positive and the ionic charge density is negative.
8340
\textit{Use:} File \texttt{SystemLabel}.TOCH is only written, not read, by siesta.
8342
\textit{Default value:} \texttt{.false.}
8344
\item[\textbf{SaveBaderCharge}] (\textit{logical}):
8345
\index{SaveBaderCharge@\textbf{SaveBaderCharge}} \index{output!Bader
8346
charge} Instructs the program to save the charge density for
8347
further post-processing by a Bader-analysis program. This ``Bader
8348
charge'' is the sum of the electronic valence charge density and a set
8349
of ``model core charges'' placed at the atomic sites. For a given
8350
atom, the model core charge has the same shape as the ``local
8351
pseudopotential charge'' used elsewhere by \siesta\ (a generalized
8352
Gaussian), but confined to a radius of 1.0 Bohr, and integrating to
8353
the total core charge ($Z$-$Z_{\mathrm{val}}$). These core charges are
8354
needed to provide local maxima for the charge density at the atomic
8355
sites, which are not guaranteed in a pseudopotential calculation. For
8356
hydrogen, an artificial core of 1 electron is added, with a
8357
confinement radius of 0.4 Bohr. The Bader charge is projected on the
8358
grid points of the mesh used by DHSCF, and saved in file
8359
\texttt{SystemLabel}.BADER. This file can be post-processed by the
8360
program \texttt{Util/grid2cube} to convert it to the ``cube'' format,
8361
accepted by several Bader-analysis programs (for example, see
8362
\url{http://theory.cm.utexas.edu/bader/}). Due to the need to
8363
represent a localized core charge, it is advisable to use a moderately
8364
high MeshCutoff when invoking this option (300-500 Ry). The size of
8365
the ``basin of attraction'' around each atom in the Bader analysis
8366
should be monitored to check that the model core charge is contained
8369
The suggested way to run the Bader analysis with the Univ. of Texas
8370
code is to use both the RHO and BADER files (both in ``cube''
8371
format), with the BADER file providing the ``reference'' and the RHO
8372
file the actual significant valence charge data which is important
8373
in bonding. (See the notes for pseudopotential codes in the above
8374
web page.) For example, for the h2o-pop example:
8376
\texttt{bader h2o-pop.RHO.cube -ref h2o-pop.BADER.cube}
8379
If netCDF support is compiled in, the file \texttt{BaderCharge.grid.nc}
8382
\textit{Default value:} \texttt{.false.}
8386
\item[\textbf{SaveInitialChargeDensity}] (\textit{logical}):
8387
\index{SaveInitialChargeDensity@\textbf{SaveInitialChargeDensity}}\index{output!charge density}
8389
If ``true'', the program generates a \texttt{SystemLabel}.RHOINIT file
8390
(and a \texttt{RhoInit.grid.nc} file if netCDF support is compiled in)
8391
containing the charge density used to start the first
8036
\begin{fdflogicalF}{SaveRho}
8037
\index{output!charge density}
8039
Instructs to write the valence pseudocharge density at the mesh used
8040
by DHSCF, in file \sysfile{RHO}.
8042
\note file \sysfile*{RHO} is only written, not read, by siesta.
8043
This file can be read by routine IORHO, which may be used by other
8044
application programs.
8046
If netCDF support is compiled in, the file \file{Rho.grid.nc} is
8051
\begin{fdflogicalF}{SaveDeltaRho}
8052
\index{output!$\delta \rho(\vec r)$}
8055
$\delta \rho(\vec r) = \rho(\vec r) - \rho_{atm}(\vec r)$, i.e., the
8056
valence pseudocharge density minus the sum of atomic valence
8057
pseudocharge densities. It is done for the mesh points used by DHSCF
8058
and it comes in file \sysfile{DRHO}. This file can be
8059
read by routine IORHO, which may be used by an application program
8062
\note file \sysfile*{DRHO} is only written, not read, by siesta.
8064
If netCDF support is compiled in, the file \file{DeltaRho.grid.nc}
8069
\begin{fdflogicalF}{SaveRhoXC}
8070
\index{output!charge density}
8072
Instructs to write the valence pseudocharge density at the mesh,
8073
including the nonlocal core corrections used to calculate the
8074
exchange-correlation energy, in file \sysfile{RHOXC}.
8076
\textit{Use:} File \sysfile*{RHOXC} is only written, not read, by
8079
If netCDF support is compiled in, the file \file{RhoXC.grid.nc} is produced.
8083
\begin{fdflogicalF}{SaveElectrostaticPotential}
8084
\index{output!electrostatic potential}
8086
Instructs to write the total electrostatic potential, defined as the
8087
sum of the hartree potential plus the local pseudopotential, at the
8088
mesh used by DHSCF, in file \sysfile{VH}. This file can be read by
8089
routine IORHO, which may be used by an application program in later
8092
\textit{Use:} File \sysfile*{VH} is only written, not read, by
8095
If netCDF support is compiled in, the file
8096
\file{ElectrostaticPotential.grid.nc} is produced.
8100
\begin{fdflogicalF}{SaveNeutralAtomPotential}
8101
\index{output!electrostatic potential}
8103
Instructs to write the neutral-atom potential, defined as the sum of
8104
the hartree potential of a ``pseudo atomic valence charge'' plus the
8105
local pseudopotential, at the mesh used by DHSCF, in file
8106
\sysfile{VNA}. It is written at the start of the
8107
self-consistency cycle, as this potential does not change.
8109
\textit{Use:} File \sysfile*{VNA} is only written, not read, by
8112
If netCDF support is compiled in, the file \file{Vna.grid.nc} is
8117
\begin{fdflogicalF}{SaveTotalPotential}
8118
\index{output!total potential}
8120
Instructs to write the valence total effective local potential
8121
(local pseudopotential + Hartree + Vxc), at the mesh used by DHSCF,
8122
in file \sysfile{VT}. This file can be read by routine
8123
IORHO, which may be used by an application program in later
8126
\textit{Use:} File \sysfile*{VT} is only written, not read, by
8129
If netCDF support is compiled in, the file
8130
\file{TotalPotential.grid.nc} is produced.
8132
\note a side effect; the vacuum level, defined as the effective
8133
potential at grid points with zero density, is printed in the
8134
standard output whenever such points exist (molecules, slabs) and
8135
either \fdf{SaveElectrostaticPotential} or \fdf{SaveTotalPotential}
8136
are \fdftrue. In a symetric (nonpolar) slab, the work function can
8137
be computed as the difference between the vacuum level and the Fermi
8142
\begin{fdflogicalF}{SaveIonicCharge}
8143
\index{output!ionic charge}
8145
Instructs to write the soft diffuse ionic charge at the mesh used by
8146
DHSCF, in file \sysfile{IOCH}. This file can be read by routine
8147
IORHO, which may be used by an application program in later
8148
versions. Remember that, within the \siesta\ sign convention, the
8149
electron charge density is positive and the ionic charge density is
8152
\textit{Use:} File \sysfile*{IOCH} is only written, not read, by siesta.
8154
If netCDF support is compiled in, the file \file{Chlocal.grid.nc} is produced.
8158
\begin{fdflogicalF}{SaveTotalCharge}
8159
\index{output!total charge}
8161
Instructs to write the total charge density (ionic+electronic) at
8162
the mesh used by DHSCF, in file \sysfile{TOCH}. This file
8163
can be read by routine IORHO, which may be used by an application
8164
program in later versions. Remember that, within the \siesta\ sign
8165
convention, the electron charge density is positive and the ionic
8166
charge density is negative.
8168
\textit{Use:} File \sysfile*{TOCH} is only written, not read, by
8173
\begin{fdflogicalF}{SaveBaderCharge}
8174
\index{output!Bader charge}
8176
Instructs the program to save the charge density for further
8177
post-processing by a Bader-analysis program. This ``Bader
8178
charge'' is the sum of the electronic valence charge density and a
8179
set of ``model core charges'' placed at the atomic sites. For a
8180
given atom, the model core charge has the same shape as the
8181
``local pseudopotential charge'' used elsewhere by \siesta\ (a
8182
generalized Gaussian), but confined to a radius of 1.0 Bohr, and
8183
integrating to the total core charge
8184
($Z$-$Z_{\mathrm{val}}$). These core charges are needed to provide
8185
local maxima for the charge density at the atomic sites, which are
8186
not guaranteed in a pseudopotential calculation. For hydrogen, an
8187
artificial core of 1 electron is added, with a confinement radius
8188
of 0.4 Bohr. The Bader charge is projected on the grid points of
8189
the mesh used by DHSCF, and saved in file
8190
\sysfile{BADER}. This file can be post-processed by the
8191
program \program{Util/grid2cube} to convert it to the ``cube''
8192
format, accepted by several Bader-analysis programs (for example,
8193
see \url{http://theory.cm.utexas.edu/bader/}). Due to the need to
8194
represent a localized core charge, it is advisable to use a
8195
moderately high MeshCutoff when invoking this option (300-500
8196
Ry). The size of the ``basin of attraction'' around each atom in
8197
the Bader analysis should be monitored to check that the model
8198
core charge is contained in it.
8200
The suggested way to run the Bader analysis with the Univ. of
8201
Texas code is to use both the RHO and BADER files (both in
8202
``cube'' format), with the BADER file providing the ``reference''
8203
and the RHO file the actual significant valence charge data which
8204
is important in bonding. (See the notes for pseudopotential codes
8205
in the above web page.) For example, for the h2o-pop example:
8207
\begin{shellexample}
8208
bader h2o-pop.RHO.cube -ref h2o-pop.BADER.cube}
8211
If netCDF support is compiled in, the file \file{BaderCharge.grid.nc}
8217
\begin{fdflogicalF}{SaveInitialChargeDensity}
8218
\index{output!charge density}
8220
If \fdftrue, the program generates a \sysfile{RHOINIT}
8221
file (and a \file{RhoInit.grid.nc} file if netCDF support is
8222
compiled in) containing the charge density used to start the first
8392
8223
self-consistency step, and it stops. Note that if an initial density
8393
8224
matrix (DM file) is used, it is not normalized. This is useful to
8394
8225
generate the charge density associated to ``partial'' DMs, as
8395
created by progras such as \texttt{dm\_creator} and \texttt{dm\_filter}.
8397
\textit{Default value:} \texttt{.false.}
8226
created by progras such as \program{dm\_creator} and
8227
\program{dm\_filter}.
8403
8233
\subsection{Auxiliary Force field}
8840
8658
Several options for MD and structural optimizations are
8841
8659
implemented, selected by
8842
\textbf{MD.TypeOfRun} (\textit{string}):
8843
\index{MD.TypeOfRun@\textbf{MD.TypeOfRun}}
8847
\item \texttt{CG} Coordinate optimization by conjugate
8848
gradients). Optionally (see variable MD.VariableCell below), the
8849
optimization can include the cell vectors.
8851
\item \texttt{Broyden} Coordinate optimization by a modified Broyden
8852
scheme). Optionally, (see variable MD.VariableCell below), the
8853
optimization can include the cell vectors.
8855
\item \texttt{FIRE} Coordinate optimization by Fast Inertial Relaxation
8856
Engine (FIRE) (E. Bitzek et al, PRL 97, 170201, (2006)).
8857
Optionally, (see variable MD.VariableCell below), the
8858
optimization can include the cell vectors.
8860
\item \texttt{Verlet} Standard Verlet algorithm MD
8862
\item \texttt{Nose} MD with temperature controlled by means of a Nos\'e
8865
\item \texttt{ParrinelloRahman} MD with pressure controlled by
8866
the Parrinello-Rahman method
8868
\item \texttt{NoseParrinelloRahman} MD with temperature controlled
8869
by means of a Nos\'e thermostat and pressure controlled by
8870
the Parrinello-Rahman method
8872
\item \texttt{Anneal} MD with annealing to a desired
8873
temperature and/or pressure (see variable MD.AnnealOption below)
8875
\item \texttt{FC} Compute force constants matrix\index{Force Constants
8876
Matrix} for phonon calculations.
8879
\item \texttt{Forces} (Receive coordinates from, and return forces to, an
8880
external driver program, using MPI, Unix pipes, or Inet sockets for
8882
The routines in module fsiesta allow the user's program to perform
8883
this communication transparently, as if siesta were a conventional
8884
force-field subroutine. See \texttt{Util/SiestaSubroutine/README} for
8885
details. WARNING: if this option is specified without a driver
8886
program sending data, siesta may hang without any notice).
8888
See directory Util/Scripting \index{Scripting} for other driving options.
8892
\textit{Default value:} \texttt{CG} (except if \textbf{compat-pre-v4-dynamics}
8893
is set. See the notes below.
8895
Note that some options specified in later variables
8896
(like quenching) modify the behavior of these MD options.
8898
Appart from being able to act as a force subroutine for a driver program that
8899
uses module fsiesta, \siesta\ is also prepared to communicate with the
8900
i-PI code (see \texttt{http://epfl-cosmo.github.io/gle4md/index.html?page=ipi}).
8901
To do this, \siesta\ must be started after i-PI (it acts as a client of
8902
i-PI, communicating with it through Inet or Unix sockets), and the following
8903
lines must be present in the .fdf data file:
8660
\begin{fdfentry}{MD.TypeOfRun}[string]<CG>
8665
Coordinate optimization by conjugate gradients). Optionally (see
8666
variable \fdf{MD.VariableCell} below), the optimization can include the
8670
Coordinate optimization by a modified Broyden scheme). Optionally,
8671
(see variable \fdf{MD.VariableCell} below), the optimization can
8672
include the cell vectors.
8675
Coordinate optimization by Fast Inertial Relaxation Engine (FIRE)
8676
(E. Bitzek et al, PRL 97, 170201, (2006)). Optionally, (see
8677
variable \fdf{MD.VariableCell} below), the optimization can
8678
include the cell vectors.
8681
Standard Verlet algorithm MD
8684
MD with temperature controlled by means of a Nos\'e
8687
\option[ParrinelloRahman]%
8688
MD with pressure controlled by the Parrinello-Rahman method
8690
\option[NoseParrinelloRahman]%
8691
MD with temperature controlled by means of a Nos\'e thermostat and
8692
pressure controlled by the Parrinello-Rahman method
8695
MD with annealing to a desired temperature and/or pressure (see
8696
variable \fdf{MD.AnnealOption} below)
8699
Compute force constants matrix\index{Force Constants Matrix} for
8700
phonon calculations.
8704
(Receive coordinates from, and return forces to, an external
8705
driver program, using MPI, Unix pipes, or Inet sockets for
8706
communication. The routines in module fsiesta allow the user's
8707
program to perform this communication transparently, as if siesta
8708
were a conventional force-field subroutine. See
8709
\shell{Util/SiestaSubroutine/README} for details. WARNING: if
8710
this option is specified without a driver program sending data,
8711
siesta may hang without any notice).
8713
See directory Util/Scripting \index{Scripting} for other driving options.
8716
\note if \fdf{Compat!Pre-v4-Dynamics} is \fdftrue\ this will default
8719
Note that some options specified in later variables (like quenching)
8720
modify the behavior of these MD options.
8722
Appart from being able to act as a force subroutine for a driver
8723
program that uses module fsiesta, \siesta\ is also prepared to
8724
communicate with the i-PI code (see
8725
\texttt{http://epfl-cosmo.github.io/gle4md/index.html?page=ipi}).
8726
To do this, \siesta\ must be started after i-PI (it acts as a client
8727
of i-PI, communicating with it through Inet or Unix sockets), and
8728
the following lines must be present in the .fdf data file:
8905
8730
MD.TypeOfRun Master # equivalent to 'Forces'
8906
8731
Master.code i-pi # ( fsiesta | i-pi )
8907
8732
Master.interface socket # ( pipes | socket | mpi )
8908
8733
Master.address localhost # or driver's IP, e.g. 150.242.7.140
8909
8734
Master.port 10001 # 10000+siesta_process_order
8910
8735
Master.socketType inet # ( inet | unix )
8913
8742
\subsection{Compatibility with pre-v4 versions}
8914
8743
\index{Backward compatibility}
8946
8784
Zmatrix-specific options take precedence. The 'MD' prefix is
8947
8785
misleading but kept for historical reasons.
8950
\item[\textbf{MD.VariableCell}] (\textit{logical}):
8951
\index{MD.VariableCell@\textbf{MD.VariableCell}} \index{cell
8952
relaxation} If true, the lattice is relaxed together with the
8953
atomic coordinates. It allows to target hydrostatic pressures or
8954
arbitrary stress tensors. See \textbf{MD.MaxStressTol}, \textbf{MD.TargetPressure}, \textbf{MD.TargetStress}, \textbf{MD.ConstantVolume}, and \textbf{MD.PreconditionVariableCell}.
8956
\textit{Use:} Used only if MD.TypeOfRun is \texttt{CG} or \texttt{Broyden} or
8959
\textit{Default value:} \texttt{.false.}
8961
\item[\textbf{MD.ConstantVolume}] (\textit{logical}):
8962
\index{MD.ConstantVolume@\textbf{MD.ConstantVolume}}
8963
\index{constant-volume cell relaxation} If true, the cell volume is
8964
kept constant in a variable-cell relaxation: only the cell shape and
8965
the atomic coordinates are allowed to change. Note that it does not
8966
make much sense to specify a target stress or pressure in this case,
8967
except for anisotropic (traceless) stresses. See \textbf{MD.VariableCell}, \textbf{MD.TargetStress}.
8969
\textit{Use:} Used only if MD.TypeOfRun is \texttt{CG} or \texttt{Broyden} or
8970
\texttt{FIRE}, and MD.VariableCell is \texttt{.true.}.
8972
\textit{Default value:} \texttt{.false.}
8974
\item[\textbf{MD.RelaxCellOnly}] (\textit{logical}):
8975
\index{MD.RelaxCellOnly@\textbf{MD.RelaxCellOnly}} \index{relaxation of
8976
cell parameters only}
8978
If true, only the cell parameters are relaxed (by the Broyden or FIRE
8979
method, not CG). The atomic coordinates are re-scaled to the new
8980
cell, keeping the fractional coordinates constant. For Zmatrix
8981
calculations, the fractional position of the first atom in each
8982
molecule is kept fixed, and no attempt is made to rescale the bond
8983
distances or angles.
8985
\textit{Use:} Used only if MD.TypeOfRun is \texttt{FIRE} or \texttt{Broyden} and
8986
MD.VariableCell is \texttt{.true.}.
8988
\textit{Default value:} \texttt{.false.}
8990
\item[\textbf{MD.MaxForceTol}] (\textit{real force}):
8991
\index{MD.MaxForceTol@\textbf{MD.MaxForceTol}}
8992
Force tolerance in coordinate optimization.
8993
Run stops if the maximum atomic force is
8994
smaller than \textbf{MD.MaxForceTol} (see \textbf{MD.MaxStressTol}
8997
\textit{Use:} Used only if MD.TypeOfRun is \texttt{CG} or \texttt{Broyden} or
9000
\textit{Default value:} \texttt{0.04 eV/Ang}
9003
\item[\textbf{MD.MaxStressTol}] (\textit{real pressure}):
9004
\index{MD.MaxStressTol@\textbf{MD.MaxStressTol}}
9005
Stress tolerance in variable-cell CG optimization. Run stops
9006
if the maximum atomic force is smaller than \textbf{MD.MaxForceTol}
9007
and the maximum stress component is smaller than \textbf{MD.MaxStressTol}.
9009
\textit{Use:} Used only if MD.TypeOfRun is \texttt{CG} or \texttt{Broyden} or
9010
\texttt{FIRE}, and Md.VariableCell is \texttt{.true.}
9012
Special consideration is needed if used with Sankey-type basis sets, since
9013
the combination of orbital kinks at the cutoff radii and the finite-grid
9014
integration originate discontinuities in the
9015
stress components, whose magnitude depends on the cutoff radii (or
9016
energy shift) and the mesh cutoff. The tolerance has to be larger
9017
than the discontinuities to avoid endless optimizations if the target
9018
stress happens to be in a discontinuity.
9020
\textit{Default value:} \texttt{1.0 GPa}
9023
\item[\textbf{MD.NumCGsteps}] (\textit{integer}):
9024
\index{MD.NumCGsteps@\textbf{MD.NumCGsteps}}
9025
Maximum number of conjugate gradient (or Broyden) minimization
9026
moves (the minimization will stop
9027
if tolerance is reached before; see MD.MaxForceTol below).
9029
\textit{Use:} Used only if MD.TypeOfRun is \texttt{CG} or \texttt{Broyden}
9031
\textit{Default value:} \texttt{0}
9033
\item[\textbf{MD.MaxCGDispl}] (\textit{real length}):
9034
\index{MD.MaxCGDispl@\textbf{MD.MaxCGDispl}}
9035
Maximum atomic displacements in an optimization move.
9037
\textit{Use:} Used only if MD.TypeOfRun is \texttt{CG} or \texttt{Broyden} or
9038
\texttt{FIRE} (despite its name). For the Broyden optimization method, it is also
9039
possible to limit indirectly the \textit{initial\/} atomic displacements
9040
using \textbf{MD.Broyden.Initial.Inverse.Jacobian}. For the FIRE method,
9041
the same result can be obtained by choosing a small time step.
9043
Note that there are Zmatrix-specific options that override this option.
9045
\textit{Default value:} \texttt{0.2 Bohr}
9048
\item[\textbf{MD.PreconditionVariableCell}] (\textit{real length}):
9049
\index{MD.PreconditionVariableCell@\textbf{MD.PreconditionVariableCell}} A length to multiply to the strain
9050
components in a variable-cell optimization. The strain components
9051
enter the minimization on the same footing as the coordinates. For
9052
good efficiency, this length should make the scale of energy
9053
variation with strain similar to the one due to atomic
9054
displacements. It is also used for the application of the \textbf{MD.MaxCGDispl} value to the strain components.
9056
\textit{Use:} Used only if MD.TypeOfRun is \texttt{CG} or \texttt{Broyden} or
9057
\texttt{FIRE} and MD.VariableCell is \texttt{.true.}
9059
\textit{Default value:} \texttt{5.0 Ang}
9061
\item[\textbf{ZM.ForceTolLength}] (\textit{real force}):
9062
\index{ZM.ForceTolLength@\textbf{ZM.ForceTolLength}} Parameter that
9063
controls the convergence with respect to forces on Z-matrix lengths
9065
\textit{Use:} This option sets the convergence criteria for the forces that
9066
act on Z-matrix components with units of length.
9068
\textit{Default value:} \texttt{$0.00155574$ Ry/Bohr}
9070
\item[\textbf{ZM.ForceTolAngle}] (\textit{torque}):
9071
\index{ZM.ForceTolAngle@\textbf{ZM.ForceTolAngle}} Parameter that
9072
controls the convergence with respect to forces on Z-matrix angles
9074
\textit{Use:} This option sets the convergence criteria for the forces that
9075
act on Z-matrix components with units of angle.
9077
\textit{Default value:} \texttt{$0.00356549$ Ry/rad}
9079
\item[\textbf{ZM.MaxDisplLength}] (\textit{real length}):
9080
\index{ZM.MaxDisplLength@\textbf{ZM.MaxDisplLength}} Parameter that
9081
controls the maximum change in a Z-matrix length during an
9084
\textit{Use:} This option sets the maximum displacement for a Z-matrix length
9086
\textit{Default value:} \texttt{0.2 Bohr}
9088
\item[\textbf{ZM.MaxDisplAngle}] (\textit{real angle}):
9089
\index{ZM.MaxDisplAngle@\textbf{ZM.MaxDisplAngle}} Parameter that
9090
controls the maximum change in a Z-matrix angle during an
9093
\textit{Use:} This option sets the maximum displacement for a Z-matrix angle
9095
\textit{Default value:} \texttt{0.003 rad }
8787
\begin{fdflogicalF}{MD.VariableCell}
8788
\index{cell relaxation}
8790
If \fdftrue, the lattice is relaxed together with the atomic
8791
coordinates. It allows to target hydrostatic pressures or arbitrary
8792
stress tensors. See \fdf{MD.MaxStressTol},
8793
\fdf{MD.TargetPressure}, \fdf{MD.TargetStress},
8794
\fdf{MD.ConstantVolume}, and
8795
\fdf{MD.PreconditionVariableCell}.
8797
\note only compatible with \fdf{MD.TypeOfRun} \fdf*{CG},
8798
\fdf*{Broyden} or \fdf*{fire}.
8803
\begin{fdflogicalF}{MD.ConstantVolume}
8804
\index{constant-volume cell relaxation}
8806
If \fdftrue, the cell volume is kept constant in a variable-cell
8807
relaxation: only the cell shape and the atomic coordinates are
8808
allowed to change. Note that it does not make much sense to specify
8809
a target stress or pressure in this case, except for anisotropic
8810
(traceless) stresses. See \fdf{MD.VariableCell},
8811
\fdf{MD.TargetStress}.
8813
\note only compatible with \fdf{MD.TypeOfRun} \fdf*{CG},
8814
\fdf*{Broyden} or \fdf*{fire}.
8818
\begin{fdflogicalF}{MD.RelaxCellOnly}
8819
\index{relaxation of cell parameters only}
8821
If \fdftrue, only the cell parameters are relaxed (by the Broyden or
8822
FIRE method, not CG). The atomic coordinates are re-scaled to the
8823
new cell, keeping the fractional coordinates constant. For
8824
\fdf{Zmatrix} calculations, the fractional position of the first
8825
atom in each molecule is kept fixed, and no attempt is made to
8826
rescale the bond distances or angles.
8828
\note only compatible with \fdf{MD.TypeOfRun} \fdf*{Broyden} or \fdf*{fire}.
8832
\begin{fdfentry}{MD.MaxForceTol}[force]<$0.04\,\mathrm{eV/Ang}$>
8834
Force tolerance in coordinate optimization.
8835
Run stops if the maximum atomic force is
8836
smaller than \fdf{MD.MaxForceTol} (see \fdf{MD.MaxStressTol}
8841
\begin{fdfentry}{MD.MaxStressTol}[pressure]<$1\,\mathrm{GPa}$>
8843
Stress tolerance in variable-cell CG optimization. Run stops if the
8844
maximum atomic force is smaller than \fdf{MD.MaxForceTol} and the
8845
maximum stress component is smaller than \fdf{MD.MaxStressTol}.
8847
Special consideration is needed if used with Sankey-type basis sets,
8848
since the combination of orbital kinks at the cutoff radii and the
8849
finite-grid integration originate discontinuities in the stress
8850
components, whose magnitude depends on the cutoff radii (or energy
8851
shift) and the mesh cutoff. The tolerance has to be larger than the
8852
discontinuities to avoid endless optimizations if the target stress
8853
happens to be in a discontinuity.
8857
\begin{fdfentry}{MD.NumCGsteps}[integer]<0>
8859
Maximum number of conjugate gradient (or Broyden) minimization moves
8860
(the minimization will stop if tolerance is reached before; see
8861
\fdf{MD.MaxForceTol} below).
8865
\begin{fdfentry}{MD.MaxCGDispl}[length]<$0.2\,\mathrm{Bohr}$>
8867
Maximum atomic displacements in an optimization move.
8869
In the Broyden optimization method, it is also possible to limit
8870
indirectly the \textit{initial\/} atomic displacements using
8871
\fdf{MD.Broyden.Initial.Inverse.Jacobian}. For the \fdf*{FIRE} method, the
8872
same result can be obtained by choosing a small time step.
8874
Note that there are Zmatrix-specific options that override this option.
8878
\begin{fdfentry}{MD.PreconditionVariableCell}[length]<$5\,\mathrm{Ang}$>
8880
A length to multiply to the strain components in a variable-cell
8881
optimization. The strain components enter the minimization on the
8882
same footing as the coordinates. For good efficiency, this length
8883
should make the scale of energy variation with strain similar to the
8884
one due to atomic displacements. It is also used for the application
8885
of the \fdf{MD.MaxCGDispl} value to the strain components.
8890
\begin{fdfentry}{ZM.ForceTolLength}[force]<$0.00155574\,\mathrm{Ry/Bohr}$>
8892
Parameter that controls the convergence with respect to forces on
8898
\begin{fdfentry}{ZM.ForceTolAngle}[torque]<$0.00356549\,\mathrm{Ry/rad}$>
8900
Parameter that controls the convergence with respect to forces on
8905
\begin{fdfentry}{ZM.MaxDisplLength}[length]<$0.2\,\mathrm{Bohr}$>
8907
Parameter that controls the maximum change in a Z-matrix length
8908
during an optimisation step.
8912
\begin{fdfentry}{ZM.MaxDisplAngle}[angle]<$0.003\,\mathrm{rad}$>
8914
Parameter that controls the maximum change in a Z-matrix angle
8915
during an optimisation step.
9099
8921
\subsubsection{Conjugate-gradients optimization}
9208
9017
the termination conditions of these methods. They run for the number
9209
9018
of MD steps requested (this is arguably a bug).
9212
\item[\textbf{MD.Quench}] (\textit{logical}):
9213
\index{MD.Quench@\textbf{MD.Quench}}
9214
Logical option to perform a power quench during the molecular dynamics.
9215
In the power quench, each velocity component is set to
9216
zero if it is opposite to the corresponding force
9217
of that component. This affects atomic velocities,
9218
or unit-cell velocities (for cell shape optimizations).
9220
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{Verlet} or
9221
\texttt{ParrinelloRahman}.
9222
It is incompatible with Nose thermostat options.
9223
The quench option allows structural relaxations of
9224
only atomic coordinates (with \textbf{MD.TypeOfRun} = \texttt{Verlet})
9225
or atomic coordinates AND cell shape
9226
(with \textbf{MD.TypeOfRun} = \texttt{ParrinelloRahman}).
9227
\textbf{MD.Quench} is superseded by \textbf{MD.FireQuench} (see below).
9229
\textit{Default value:} \texttt{.false.}
9231
\item[\textbf{MD.FireQuench}] (\textit{logical}) (Deprecated)
9232
\index{MD.FireQuench@\textbf{MD.FireQuench}}
9234
SEE the new option \texttt{FIRE} for \texttt{MD.TypeOfRun}
9236
Logical option to perform a FIRE quench during a Verlet molecular dynamics
9237
run, as described by Bitzek \textit{et al.} in Phys. Rev. Lett. \textbf{97},
9238
170201 (2006). It is a relaxation algorithm, and thus the dynamics
9239
are of no interest per se: the initial time-step can be played with
9240
(it uses \textbf{MD.LengthTimeStep} as initial $\Delta t$),
9241
as well as the initial temperature (recommended 0) and the atomic
9242
masses (recommended equal). Preliminary tests seem to indicate that
9243
the combination of $\Delta t = 5$ fs and a value of 20 for the atomic
9244
masses works reasonably. The dynamics stops when the force
9245
tolerance is reached (\textbf{MD.MaxForceTol}). The other
9246
parameters controlling the algorithm (initial damping,
9247
increase and decrease thereof etc.) are hardwired in the code,
9248
at the recommended values in the cited paper,
9249
including $\Delta t_{max} = 10$ fs.
9251
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{Verlet}.
9252
It is incompatible with Nose thermostat options. No variable
9253
cell option implemented for this at this stage.
9254
\textbf{MD.FireQuench} supersedes \textbf{MD.Quench}. This option is
9255
deprecated. The new option \texttt{FIRE} for \texttt{MD.TypeOfRun} should be
9258
\textit{Default value:} \texttt{.false.}
9020
\begin{fdflogicalF}{MD.Quench}
9022
Logical option to perform a power quench during the molecular
9023
dynamics. In the power quench, each velocity component is set to
9024
zero if it is opposite to the corresponding force of that
9025
component. This affects atomic velocities, or unit-cell velocities
9026
(for cell shape optimizations).
9028
\note only applicable for \fdf{MD.TypeOfRun} \fdf*{Verlet} or
9029
\fdf*{ParrinelloRahman}.
9031
It is incompatible with Nose thermostat options.
9033
\note \fdf{MD.Quench} is superseded by \fdf{MD.FireQuench} (see
9039
\begin{fdflogicalF}{MD.FireQuench}
9041
See the new option \fdf*{FIRE} for \fdf{MD.TypeOfRun}.
9043
Logical option to perform a FIRE quench during a Verlet molecular
9044
dynamics run, as described by Bitzek \textit{et al.} in
9045
Phys. Rev. Lett. \textbf{97}, 170201 (2006). It is a relaxation
9046
algorithm, and thus the dynamics are of no interest per se: the
9047
initial time-step can be played with (it uses
9048
\fdf{MD.LengthTimeStep} as initial $\Delta t$), as well as the
9049
initial temperature (recommended 0) and the atomic masses
9050
(recommended equal). Preliminary tests seem to indicate that the
9051
combination of $\Delta t = 5$ fs and a value of 20 for the atomic
9052
masses works reasonably. The dynamics stops when the force tolerance
9053
is reached (\fdf{MD.MaxForceTol}). The other parameters
9054
controlling the algorithm (initial damping, increase and decrease
9055
thereof etc.) are hardwired in the code, at the recommended values
9056
in the cited paper, including $\Delta t_{max} = 10$ fs.
9058
Only available for \fdf{MD.TypeOfRun} \fdf*{Verlet}
9059
It is incompatible with Nose thermostat options. No variable
9060
cell option implemented for this at this stage.
9061
\fdf{MD.FireQuench} supersedes \fdf{MD.Quench}. This option is
9062
deprecated. The new option \fdf*{FIRE} for \fdf{MD.TypeOfRun} should be
9262
9071
\subsection{Target stress options}
9264
9073
Useful for structural optimizations and constant-pressure molecular
9270
\item[\textbf{MD.TargetPressure}] (\textit{real pressure}):
9271
\index{MD.TargetPressure@\textbf{MD.TargetPressure}}
9272
Target pressure for Parrinello-Rahman method, variable cell optimizations,
9273
and annealing options.
9275
\textit{Use:} Used only if MD.TypeOfRun =
9276
\texttt{ParrinelloRahman}, \texttt{NoseParrinelloRahman},
9277
\texttt{CG, Broyden, or FIRE} (variable cell), or \texttt{Anneal}
9278
(if MD.AnnealOption = \texttt{Pressure} or \texttt{TemperatureandPressure})
9280
\textit{Default value:} \texttt{0.0 GPa}
9283
\item[\textbf{MD.TargetStress}] (\textit{data block}): External or
9284
target stress tensor for variable cell optimizations. Stress
9285
components are given in a line, in the order \texttt{xx, yy, zz, xy,
9286
xz, yz}. In units of \textbf{MD.TargetPressure}, but with the
9287
opposite sign. For example, a uniaxial compressive stress of 2 GPa
9288
along the 100 direction would be given by
9290
MD.TargetPressure 2. GPa
9291
%block MD.TargetStress
9292
-1.0 0.0 0.0 0.0 0.0 0.0
9293
%endblock MD.TargetStress
9296
\textit{Use:} Used only if MD.TypeOfRun is \texttt{CG, Broyden, or FIRE} and
9297
MD.VariableCell is \texttt{.true.}
9299
\textit{Default value:} Hydrostatic target pressure:
9300
\texttt{-1., -1., -1., 0., 0., 0.}
9302
\item[\textbf{MD.RemoveIntramolecularPressure}] (\textit{logical}):
9303
\index{MD.RemoveIntramolecularPressure@\textbf{MD.RemoveIntramolecularPressure}}
9304
\index{removal of intramolecular pressure}
9306
If \texttt{.true.}, the contribution to the stress coming from the
9307
internal degrees of freedom of the molecules will be subtracted from
9308
the stress tensor used in variable-cell optimization or variable-cell
9309
molecular-dynamics. This is done in an approximate manner, using the
9310
virial form of the stress, and assumming that the ``mean force'' over
9311
the coordinates of the molecule represents the ``inter-molecular''
9312
stress. The correction term was already computed in earlier versions
9313
of \siesta\ and used to report the ``molecule pressure''. The
9314
correction is now computed molecule-by-molecule if the Zmatrix format
9317
If the intra-molecular stress is removed, the corrected static and
9318
total stresses are printed in addition to the uncorrected items.
9319
The corrected Voigt form is also printed.
9321
\textit{Default value:} \texttt{.false.}
9076
\begin{fdfentry}{MD.TargetPressure}[pressure]<$0\,\mathrm{GPa}$>
9078
Target pressure for Parrinello-Rahman method, variable cell
9079
optimizations, and annealing options.
9081
\note this is only compatible with
9082
\fdf{MD.TypeOfRun} \fdf*{ParrinelloRahman}, \fdf*{NoseParrinelloRahman},
9083
\fdf*{CG}, \fdf*{Broyden} or \fdf*{FIRE} (variable cell), or \fdf*{Anneal}
9084
(if \fdf{MD.AnnealOption} \fdf*{Pressure} or \fdf*{TemperatureandPressure}).
9089
\begin{fdfentry}{MD.TargetStress}[block]<$-1 -1 -1 0 0 0$>
9091
External or target stress tensor for variable cell optimizations.
9092
Stress components are given in a line, in the order \texttt{xx, yy,
9093
zz, xy, xz, yz}. In units of \fdf{MD.TargetPressure}, but
9094
with the opposite sign. For example, a uniaxial compressive stress
9095
of 2 GPa along the 100 direction would be given by
9097
MD.TargetPressure 2. GPa
9098
%block MD.TargetStress
9099
-1.0 0.0 0.0 0.0 0.0 0.0
9100
%endblock MD.TargetStress
9103
Only used if \fdf{MD.TypeOfRun} is \fdf*{CG}, \fdf*{Broyden} or
9104
\fdf*{FIRE} and \fdf{MD.VariableCell} is \fdftrue.
9109
\begin{fdflogicalF}{MD.RemoveIntramolecularPressure}
9110
\index{removal of intramolecular pressure}
9112
If \fdftrue, the contribution to the stress coming from the internal
9113
degrees of freedom of the molecules will be subtracted from the
9114
stress tensor used in variable-cell optimization or variable-cell
9115
molecular-dynamics. This is done in an approximate manner, using
9116
the virial form of the stress, and assumming that the ``mean force''
9117
over the coordinates of the molecule represents the
9118
``inter-molecular'' stress. The correction term was already computed
9119
in earlier versions of \siesta\ and used to report the ``molecule
9120
pressure''. The correction is now computed molecule-by-molecule if
9121
the Zmatrix format is used.
9123
If the intra-molecular stress is removed, the corrected static and
9124
total stresses are printed in addition to the uncorrected items.
9125
The corrected Voigt form is also printed.
9325
9131
\subsection{Molecular dynamics}
9328
9134
the cell vectors) in response to the forces (and stresses), using the
9329
9135
classical equations of motion.
9331
Note that the Zmatrix input option (see Sec.~\ref{sec:Zmatrix}) is not
9137
Note that the \fdf{Zmatrix} input option (see Sec.~\ref{sec:Zmatrix}) is not
9332
9138
compatible with molecular dynamics. The initial geometry can be
9333
9139
specified using the Zmatrix format, but the Zmatrix generalized
9334
9140
coordinates will not be updated.
9337
\item[\textbf{MD.InitialTimeStep}] (\textit{integer}):
9338
\index{MD.InitialTimeStep@\textbf{MD.InitialTimeStep}}
9339
Initial time step of the MD simulation.
9340
In the current version of \siesta\ it must be 1.
9342
\textit{Use:} Used only if MD.TypeOfRun is not \texttt{CG} or \texttt{Broyden}
9344
\textit{Default value:} \texttt{1}
9346
\item[\textbf{MD.FinalTimeStep}] (\textit{integer}):
9347
\index{MD.FinalTimeStep@\textbf{MD.FinalTimeStep}}
9348
Final time step of the MD simulation.
9350
\textit{Default value:} \texttt{1}
9352
\item[\textbf{MD.LengthTimeStep}] (\textit{real time}):
9353
\index{MD.LengthTimeStep@\textbf{MD.LengthTimeStep}}
9354
Length of the time step of the MD simulation.
9356
\textit{Default value:} \texttt{1.0 fs}
9358
\item[\textbf{MD.InitialTemperature}] (\textit{real temperature or energy}):
9359
\index{MD.InitialTemperature@\textbf{MD.InitialTemperature}}
9360
Initial temperature for the MD run. The atoms are assigned random
9361
velocities drawn from the Maxwell-Bolzmann distribution with the
9362
corresponding temperature. The constraint of zero center of
9363
mass velocity is imposed.
9365
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{Verlet, Nose,
9366
ParrinelloRahman, NoseParrinelloRahman}
9369
\textit{Default value:} \texttt{0.0 K}
9372
\item[\textbf{MD.TargetTemperature}] (\textit{real temperature or energy}):
9373
\index{MD.TargetTemperature@\textbf{MD.TargetTemperature}}
9374
Target temperature for Nose thermostat and annealing options.
9376
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{Nose, NoseParrinelloRahman}
9377
or \texttt{Anneal} (if \textbf{MD.AnnealOption} = \texttt{Temperature} or
9378
\texttt{TemperatureandPressure})
9380
\textit{Default value:} \texttt{0.0 K}
9384
\item[\textbf{MD.NoseMass}] (\textit{real moment of inertia}):
9385
\index{MD.NoseMass@\textbf{MD.NoseMass}}
9386
Generalized mass of Nose variable.
9387
This determines the time scale of the Nose variable
9388
dynamics, and the coupling of the thermal bath to
9389
the physical system.
9391
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{Nose} or
9392
\texttt{NoseParrinelloRahman}
9394
\textit{Default value:} \texttt{100.0 Ry*fs**2}
9396
\item[\textbf{MD.ParrinelloRahmanMass}] (\textit{real moment of inertia}):
9397
\index{MD.ParrinelloRahmanMass@\textbf{MD.ParrinelloRahmanMass}}
9143
\begin{fdfentry}{MD.InitialTimeStep}[integer]<$1$>
9145
Initial time step of the MD simulation. In the current version of
9146
\siesta\ it must be 1.
9148
Used only if \fdf{MD.TypeOfRun} is not \fdf*{CG} or \fdf*{Broyden}.
9152
\begin{fdfentry}{MD.FinalTimeStep}[integer]<$1$>
9154
Final time step of the MD simulation.
9159
\begin{fdfentry}{MD.LengthTimeStep}[time]<$1\,\mathrm{fs}$>
9161
Length of the time step of the MD simulation.
9165
\begin{fdfentry}{MD.InitialTemperature}[temperature/energy]<$0\,\mathrm K$>
9167
Initial temperature for the MD run. The atoms are assigned random
9168
velocities drawn from the Maxwell-Bolzmann distribution with the
9169
corresponding temperature. The constraint of zero center of mass
9170
velocity is imposed.
9172
\note only used if \fdf{MD.TypeOfRun} \fdf*{Verlet}, \fdf*{Nose},
9173
\fdf*{ParrinelloRahman}, \fdf*{NoseParrinelloRahman} or
9178
\begin{fdfentry}{MD.TargetTemperature}[temperature/energy]<$0\,\mathrm K$>
9180
Target temperature for Nose thermostat and annealing options.
9182
\note only used if \fdf{MD.TypeOfRun} \fdf*{Nose},
9183
\fdf*{NoseParrinelloRahman} or
9184
\fdf*{Anneal} if \fdf{MD.AnnealOption} is \fdf*{Temperature} or
9185
\fdf*{TemperatureandPressure}.
9189
\begin{fdfentry}{MD.NoseMass}[moment of inertia]<$100\,\mathrm{Ry\,fs^2}$>
9191
Generalized mass of Nose variable. This determines the time scale
9192
of the Nose variable dynamics, and the coupling of the thermal bath
9193
to the physical system.
9195
Only used for Nose MD runs.
9199
\begin{fdfentry}{MD.ParrinelloRahmanMass}[moment of inertia]<$100\,\mathrm{Ry\,fs^2}$>
9398
9201
Generalized mass of Parrinello-Rahman variable. This determines the
9399
9202
time scale of the Parrinello-Rahman variable dynamics, and its
9400
9203
coupling to the physical system.
9402
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{ParrinelloRahman}
9403
or \texttt{NoseParrinelloRahman}
9405
\textit{Default value:} \texttt{100.0 Ry*fs**2}
9408
\textit{Default value:} Same as \textbf{NumberOfAtoms}
9410
\item[\textbf{MD.AnnealOption}] (\textit{string}):
9411
\index{MD.AnnealOption@\textbf{MD.AnnealOption}} Type of annealing MD
9412
to perform. The target temperature or pressure are achieved by
9413
velocity and unit cell rescaling, in a given time determined by the
9414
variable \textbf{MD.TauRelax} below.
9416
\item \texttt{Temperature} (Reach a target temperature by velocity
9418
\item \texttt{Pressure} (Reach a target pressure by scaling of the unit
9419
cell size and shape)
9420
\item \texttt{TemperatureandPressure} (Reach a target temperature and
9421
pressure by velocity rescaling and by scaling of the unit cell size
9425
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{Anneal}
9427
\textit{Default value:} \texttt{TemperatureAndPressure}
9429
\item[\textbf{MD.TauRelax}] (\textit{real time}): \index{MD.TauRelax@\textbf{MD.TauRelax}} Relaxation time to reach target temperature and/or
9430
pressure in annealing MD. Note that this is a ``relaxation time'',
9431
and as such it gives a rough estimate of the time needed to achieve
9432
the given targets. As a normal simulation also exhibits
9433
oscillations, the actual time needed to reach the \textit{averaged}
9434
targets will be significantly longer.
9436
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{Anneal}
9438
\textit{Default value:} \texttt{100.0 fs}
9440
\item[\textbf{MD.BulkModulus}] (\textit{real pressure}):
9441
\index{MD.BulkModulus@\textbf{MD.BulkModulus}} Estimate (may be rough)
9442
of the bulk modulus of the system. This is needed to set the rate
9443
of change of cell shape to reach target pressure in annealing MD.
9445
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{Anneal}, when
9446
\textbf{MD.AnnealOption} = \texttt{Pressure} or \texttt{TemperatureAndPressure}
9448
\textit{Default value:} \texttt{100.0 Ry/Bohr**3}
9205
Only used for Parrinello-Rahman MD runs.
9209
\begin{fdfentry}{MD.AnnealOption}[string]<TemperatureAndPressure>
9211
Type of annealing MD to perform. The target temperature or pressure
9212
are achieved by velocity and unit cell rescaling, in a given time
9213
determined by the variable \fdf{MD.TauRelax} below.
9215
\option[Temperature]%
9216
Reach a target temperature by velocity rescaling
9219
Reach a target pressure by scaling of the unit cell size and shape
9221
\option[TemperatureandPressure]%
9222
Reach a target temperature and pressure by velocity rescaling and
9223
by scaling of the unit cell size and shape
9226
Only applicable for \fdf{MD.TypeOfRun} \fdf*{Anneal}.
9230
\begin{fdfentry}{MD.TauRelax}[time]<$100\,\mathrm{fs}$>
9232
Relaxation time to reach target temperature and/or pressure in
9233
annealing MD. Note that this is a ``relaxation time'', and as such
9234
it gives a rough estimate of the time needed to achieve the given
9235
targets. As a normal simulation also exhibits oscillations, the
9236
actual time needed to reach the \emph{averaged} targets will be
9237
significantly longer.
9239
Only applicable for \fdf{MD.TypeOfRun} \fdf*{Anneal}.
9243
\begin{fdfentry}{MD.BulkModulus}[pressure]<$100\,\mathrm{Ry/Bohr^3}$>
9245
Estimate (may be rough) of the bulk modulus of the system. This is
9246
needed to set the rate of change of cell shape to reach target
9247
pressure in annealing MD.
9249
Only applicable for \fdf{MD.TypeOfRun} \fdf*{Anneal}, when
9250
\fdf{MD.AnnealOption} is \fdf*{Pressure} or \fdf*{TemperatureAndPressure}
9453
9256
\subsection{Output options for dynamics}
9455
9258
Every time the atoms move, either during coordinate relaxation or
9456
molecular dynamics, their positions \textbf{predicted for next step} and
9457
\textbf{current} velocities are stored in file SystemLabel.XV. The shape
9458
of the unit cell and its associated 'velocity' (in Parrinello-Rahman
9459
dynamics) are also stored in this file.
9461
Other options follow.
9464
\item[\textbf{WriteCoorInitial}] (\textit{logical}):
9465
\index{WriteCoorInitial@\textbf{WriteCoorInitial}}
9466
\index{output!atomic coordinates!initial}
9467
It determines whether the initial atomic coordinates of the simulation are
9468
dumped into the main output file. These coordinates correspond to the
9469
ones actually used in the first step (see the section on precedence
9470
issues in structural input) and are output in Cartesian coordinates in
9473
It is not affected by the setting of \textbf{LongOutput}.
9475
\textit{Default value:} \texttt{.true.}
9478
\item[\textbf{WriteCoorStep}] (\textit{logical}): \index{WriteCoorStep@\textbf{WriteCoorStep}} \index{output!atomic coordinates!in a dynamics
9479
step} If \texttt{.true.} it writes the atomic coordinates to standard
9259
molecular dynamics, their positions \textbf{predicted for next step}
9260
and \textbf{current} velocities are stored in file \sysfile{XV}. The
9261
shape of the unit cell and its associated 'velocity' (in
9262
Parrinello-Rahman dynamics) are also stored in this file.
9264
\begin{fdflogicalT}{WriteCoorInitial}
9265
\index{output!atomic coordinates!initial}
9267
It determines whether the initial atomic coordinates of the
9268
simulation are dumped into the main output file. These coordinates
9269
correspond to the ones actually used in the first step (see the
9270
section on precedence issues in structural input) and are output in
9271
Cartesian coordinates in Bohr units.
9273
It is not affected by the setting of \fdf{LongOutput}.
9278
\begin{fdflogicalF}{WriteCoorStep}
9279
\index{output!atomic coordinates!in a dynamics step}
9281
If \fdftrue, it writes the atomic coordinates to standard
9480
9282
output at every MD time step or relaxation step. The coordinates are
9481
always written in the \textit{Systemlabel}.XV file, but overriden at
9482
every step. They can be also accumulated in the \textit{Systemlabel}.MD
9483
or \textit{Systemlabel}.MDX files depending on \textbf{WriteMDhistory}.
9486
\textit{Default value:} \texttt{.false.} (see \textbf{LongOutput})
9489
\item[\textbf{WriteForces}] (\textit{logical}): \index{WriteForces@\textbf{WriteForces}}\index{output!forces} If \texttt{.true.} it writes the
9490
atomic forces to the output file at every MD time step or relaxation
9491
step. Note that the forces of the last step can be found in the
9492
file \textit{Systemlabel}.FA. If constraints are used, the file
9493
\textit{SystemLabel}.\texttt{FAC} is also written.
9496
\textit{Default value:} \texttt{.false.} (see \textbf{LongOutput})
9498
\item[\textbf{WriteMDhistory}] (\textit{logical}):
9499
\index{WriteMDhistory@\textbf{WriteMDhistory}} \index{output!molecular
9500
dynamics!history} If \texttt{.true.} \siesta\ accumulates the
9501
molecular dynamics trajectory in the following files:
9504
\textit{Systemlabel}.MD : atomic coordinates and velocities (and lattice
9505
vectors and their time derivatives, if the dynamics implies variable
9506
cell). The information is stored unformatted for postprocessing with
9507
utility programs to analyze the MD trajectory.
9509
\textit{Systemlabel}.MDE : shorter description of the run, with energy,
9510
temperature, etc., per time step.
9512
These files are accumulative even for different runs.
9514
\textit{Default value:} \texttt{.false.}
9516
\index{output!molecular dynamics!history}
9518
The trajectory of a molecular dynamics run (or a conjugate gradient
9519
minimization) can be accumulated in different files: SystemLabel.MD,
9520
SystemLabel.MDE, and SystemLabel.ANI. The first file keeps the whole
9521
trajectory information, meaning positions and velocities at every time
9522
step, including lattice vectors if the cell varies. NOTE that the
9523
positions (and maybe the cell vectors) stored at each time step are
9524
the \textbf{predicted} values for the next step. Care should be taken if
9525
joint position-velocity correlations need to be computed from this
9526
file. The second gives global information (energy, temperature, etc),
9527
and the third has the coordinates in a form suited for XMol animation.
9528
See the \textbf{WriteMDhistory} and \textbf{WriteMDXmol} data descriptors
9529
above for information. \siesta\ always appends new information on
9530
these files, making them accumulative even for different runs.
9532
The \texttt{iomd} subroutine can generate both an unformatted file
9533
SystemLabel.MD (default) or ASCII formatted files SystemLabel.MDX and
9534
SystemLabel.MDC containing the atomic and lattice trajectories,
9535
respectively. Edit the file to change the settings if desired.
9283
always written in the \sysfile{XV} file, but overriden at
9284
every step. They can be also accumulated in the \sysfile*{MD}
9285
or \sysfile{MDX} files depending on \fdf{WriteMDHistory}.
9289
\begin{fdflogicalF}{WriteForces}
9290
\index{output!forces}
9292
If \fdftrue, it writes the atomic forces to the output file at every
9293
MD time step or relaxation step. Note that the forces of the last
9294
step can be found in the file \sysfile{FA}. If constraints are used,
9295
the file \sysfile{FAC} is also written.
9299
\begin{fdflogicalF}{WriteMDHistory}
9300
\index{output!molecular dynamics!history}
9302
If \fdftrue, \siesta\ accumulates the molecular dynamics trajectory
9303
in the following files:
9306
\sysfile{MD} : atomic coordinates and velocities (and lattice
9307
vectors and their time derivatives, if the dynamics implies
9308
variable cell). The information is stored unformatted for
9309
postprocessing with utility programs to analyze the MD trajectory.
9312
\sysfile{MDE} : shorter description of the run, with energy,
9313
temperature, etc., per time step.
9316
These files are accumulative even for different runs.
9318
\index{output!molecular dynamics!history}
9320
The trajectory of a molecular dynamics run (or a conjugate gradient
9321
minimization) can be accumulated in different files: SystemLabel.MD,
9322
SystemLabel.MDE, and SystemLabel.ANI. The first file keeps the whole
9323
trajectory information, meaning positions and velocities at every time
9324
step, including lattice vectors if the cell varies. NOTE that the
9325
positions (and maybe the cell vectors) stored at each time step are
9326
the \textbf{predicted} values for the next step. Care should be taken if
9327
joint position-velocity correlations need to be computed from this
9328
file. The second gives global information (energy, temperature, etc),
9329
and the third has the coordinates in a form suited for XMol animation.
9330
See the \fdf{WriteMDHistory} and \fdf{WriteMDXmol} data descriptors
9331
above for information. \siesta\ always appends new information on
9332
these files, making them accumulative even for different runs.
9334
The \program{iomd} subroutine can generate both an unformatted file
9335
\sysfile*{MD} (default) or ASCII formatted files \sysfile*{MDX} and
9336
\sysfile*{MDC} containing the atomic and lattice trajectories,
9337
respectively. Edit the file to change the settings if desired.
9540
9343
\subsection{Restarting geometry optimizations and MD runs}
9834
9637
\subsection{Phonon calculations}
9836
If \textbf{MD.TypeOfRun} = \texttt{FC}, \siesta\ sets up a special
9837
outer geometry loop that displaces individual atoms along the
9838
coordinate directions to build the force-constant matrix.
9639
If \fdf{MD.TypeOfRun} is \fdf*{FC}, \siesta\ sets up a special outer
9640
geometry loop that displaces individual atoms along the coordinate
9641
directions to build the force-constant matrix.
9839
9642
\index{output!molecular dynamics!Force Constants Matrix}
9841
9644
The output (see below) can be analyzed to extract phonon frequencies
9842
and vectors with the VIBRA\index{VIBRA} package in the \texttt{Util/Vibra}
9645
and vectors with the VIBRA\index{VIBRA} package in the \program{Util/Vibra}
9843
9646
directory. For computing the Born effective charges together with the
9844
force constants, see \textbf{BornCharge} \index{BornCharge@\textbf{BornCharge}}).
9847
\item[\textbf{MD.FCDispl}] (\textit{real length}):
9848
\index{MD.FCDispl@\textbf{MD.FCDispl}}
9849
Displacement to use for the computation of the force constant
9850
matrix\index{Force Constants Matrix} for phonon calculations.
9852
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{FC}.
9854
\textit{Default value:} \texttt{0.04 Bohr}
9856
\item[\textbf{MD.FCfirst}] (\textit{integer}):
9857
\index{MD.FCfirst@\textbf{MD.FCfirst}}
9858
Index of first atom to displace for the computation of the force constant
9859
matrix\index{Force Constants Matrix} for phonon calculations.
9861
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{FC}.
9863
\textit{Default value:} \texttt{1}
9865
\item[\textbf{MD.FClast}] (\textit{integer}):
9866
\index{MD.FClast@\textbf{MD.FClast}}
9867
Index of last atom to displace for the computation of the force constant
9868
matrix\index{Force Constants Matrix} for phonon calculations.
9870
\textit{Use:} Used only if \textbf{MD.TypeOfRun} = \texttt{FC}.
9873
The force-constants matrix is written in file \textit{SystemLabel}.\texttt{FC}.
9874
The format is the following: for the displacement of
9875
each atom in each direction, the forces on each of the other
9876
atoms is writen (divided by the value of the displacement),
9877
in units of eV/\AA$^2$. Each line has the forces in the $x$, $y$
9878
and $z$ direction for one of the atoms.
9880
(If constraints are used, the file \textit{SystemLabel}.\texttt{FCC} is
9647
force constants, see \fdf{BornCharge}.
9649
\begin{fdfentry}{MD.FCDispl}[length]<$0.04\,\mathrm{Bohr}$>
9651
Displacement to use for the computation of the force constant
9652
matrix\index{Force Constants Matrix} for phonon calculations.
9656
\begin{fdfentry}{MD.FCFirst}[integer]<$1$>
9658
Index of first atom to displace for the computation of the force
9659
constant matrix\index{Force Constants Matrix} for phonon
9664
\begin{fdfentry}{MD.FCLast}[integer]<\fdfvalue{MD.FCFirst}>
9666
Index of last atom to displace for the computation of the force
9667
constant matrix\index{Force Constants Matrix} for phonon
9672
The force-constants matrix is written in file \sysfile{FC}. The
9673
format is the following: for the displacement of each atom in each
9674
direction, the forces on each of the other atoms is writen (divided by
9675
the value of the displacement), in units of eV/\AA$^2$. Each line has
9676
the forces in the $x$, $y$ and $z$ direction for one of the atoms.
9678
If constraints are used, the file \sysfile{FCC} is also written.
9883
9681
\section{LDA+U}
9884
Important note: Current implementation is based on the
9885
simplified rotationally invariant LDA+U formulation
9886
of Dudarev and collaborators
9684
Important note: Current implementation is based on the simplified
9685
rotationally invariant LDA+U formulation of Dudarev and collaborators
9887
9686
[see, Dudarev \textit{et al.}, Phys. Rev. B \textbf{57}, 1505 (1998)].
9888
Although the input allows to define
9889
independent values of the U and J parameters for
9891
in the actual calculation
9892
the two parameters are combined to produce
9893
an effective Coulomb repulsion U$_{\mathrm{eff}}$=U-J. U$_{\mathrm{eff}}$
9894
is the parameter actually used in the calculations
9687
Although the input allows to define independent values of the U and J
9688
parameters for each atomic shell, in the actual calculation the two
9689
parameters are combined to produce an effective Coulomb repulsion
9690
$U_{\mathrm{eff}}=U-J$. $U_{\mathrm{eff}}$ is the parameter actually
9691
used in the calculations for the time being.
9897
For large or intermediate values of U$_{\mathrm{eff}}$ the convergence
9693
For large or intermediate values of $U_{\mathrm{eff}}$ the convergence
9898
9694
is sometimes difficult. A step-by-step increase of the
9899
value of U$_{\mathrm{eff}}$ can be advisable in such cases.
9905
\item[\textbf{LDAU.ProjectorGenerationMethod}] (\textit{integer}):
9906
Generation method of the LDA+U projectors. The LDA+U projectors are the
9907
localized functions used
9908
to calculate the local populations used in a Hubbard-like term
9909
that modifies the LDA Hamiltonian and energy. It is important
9910
to recall that LDA+U projectors should be quite localized functions.
9911
Otherwise the calculated populations loose their atomic character
9912
and physical meaning. Even more importantly,
9913
the interaction range can increase so much that jeopardizes
9914
the efficiency of the calculation.
9916
Two methods are currently implemented (accepted values are 1 and 2):
9918
\item 1. Projectors are slightly-excited numerical atomic orbitals
9919
similar to those used as an automatic basis set by SIESTA.
9920
The radii of these orbitals are controlled using
9921
the parameter LDAU.EnergyShift and/or the data included in the
9922
block LDAU.proj (quite similar to the data block PAO.Basis used
9923
to specify the basis set, see below).
9924
\item 2. Projectors are exact solutions of the pseudoatomic
9925
problem (and, in principle, are not strictly localized) which are
9926
cut using a Fermi function $1/\{1+\exp[(r-r_c)\omega]\}$.
9927
The values of $r_c$ and $\omega$ are controlled using
9928
the parameter LDAU.CutoffNorm and/or the data included in the
9932
\textit{Default value:} \texttt{2}
9934
\item[\textbf{LDAU.EnergyShift}] (\textit{physical }):
9935
Energy increase used to define the localization radius
9936
of the LDA+U projectors (similar to the parameter
9939
\textit{Use:} Used only if LDAU.ProjectorGenerationMethod is \texttt{1 }
9941
\textit{Default value:} \texttt{0.05 Ry}
9943
\item[\textbf{LDAU.CutoffNorm}] (\textit{real}):
9944
Parameter used to define the value of $r_c$ used in the Fermi distribution to
9945
cut the LDA+U projectors generated according to generation
9946
method 2 (see above). LDAU.CutoffNorm is
9947
the norm of the original pseudoatomic orbital
9948
contained inside a sphere of radius equal to $r_c$.
9950
\textit{Use:} Used only if LDAU.ProjectorGenerationMethod is \texttt{2 }
9952
\textit{Default value:} \texttt{0.90 }
9954
\item[\textbf{LDAU.proj}] (\textit{data block}):
9955
Data block used to specify the LDA+U projectors.
9958
\item If LDAU.ProjectorGenerationMethod is \texttt{1 }, the
9959
syntax is as follows:
9962
%block LDAU.proj # Define LDA+U projectors
9963
Fe 2 # Label, l_shells
9964
n=3 2 E 50.0 2.5 # n (opt if not using semicore levels),l,Softconf(opt)
9965
5.00 0.35 # U(eV), J(eV) for this shell
9967
0.95 # scaleFactor (opt)
9969
1.00 0.05 # U(eV), J(eV) for this shell
9970
0.00 # rc(Bohr) (if 0, automatic r_c from LDAU.EnergyShift)
9976
\item If LDAU.ProjectorGenerationMethod is \texttt{2 }, the
9977
sintaxis is as follows:
9980
%block LDAU.proj # Define LDAU projectors
9981
Fe 2 # Label, l_shells
9982
n=3 2 E 50.0 2.5 # n (opt if not using semicore levels),l,Softconf(opt)
9983
5.00 0.35 # U(eV), J(eV) for this shell
9984
2.30 0.15 # rc (Bohr), \omega(Bohr) (Fermi cutoff function)
9985
0.95 # scaleFactor (opt)
9987
1.00 0.05 # U(eV), J(eV) for this shell
9988
0.00 0.00 # rc(Bohr), \omega(Bohr) (if 0 r_c from LDAU.CutoffNorm
9989
%endblock LDAU.proj # and \omega from default value)
9993
\textit{Default values of the different parameters:} \\
9994
U = \texttt{0.0 eV} \\
9995
J = \texttt{0.0 eV} \\
9996
$\omega$ = \texttt{0.05 Bohr} \\
9997
$r_c$ = see LDAU.EnergyShift or LDAU.CutoffNorm depending on generation method\\
9998
Scale factor = \texttt{1.0 }
10000
\textit{Default value:} None.
10002
\item[\textbf{LDAU.FirstIteration}] (\textit{logical})
10003
If .true. local populations are calculated and
10004
Hubbard-like term is switch on in the first iteration.
10005
Useful if restarting a calculation reading a converged
10006
or an almost converged density matrix from file.
10008
\textit{Default value:} \texttt{.false.}
10010
\item[\textbf{LDAU.ThresholdTol}] (\textit{real}):
10011
Local populations only calculated and/or updated
10012
if the change in the density matrix elements (dDmax)
10013
is lower than LDAU.ThresholdTol.
10015
\textit{Default value:} \texttt{1.0$^{-2}$}
10017
\item[\textbf{LDAU.PopTol}] (\textit{real}):
10018
Convergence criterium for the LDA+U local populations.
10019
In the current implementation the
10020
Hubbard-like term of the Hamiltonian is only updated (except
10021
for the last iteration)
10022
if the variations of the local populations are larger than
10025
\textit{Default value:} \texttt{1.0$^{-3}$}
10027
\item[\textbf{LDAU.PotentialShift}] (\textit{logical}):
10028
If set to .true. the value given to the $U$ parameter in the
10030
is interpreted as a local potential shift. Recording
10031
the change of the local populations as a function of
10032
this potential shift, we can calculate the appropriate
10033
value of U for the system under study following the methology
10034
proposed by Cococcioni and Gironcoli in Phys. Rev. B \textbf{71},
10036
\textit{Default value:} \texttt{.false.}
9695
value of $U_{\mathrm{eff}}$ can be advisable in such cases.
9697
\begin{fdfentry}{LDAU.ProjectorGenerationMethod}[integer]<2>
9699
Generation method of the LDA+U projectors. The LDA+U projectors are
9700
the localized functions used to calculate the local populations used
9701
in a Hubbard-like term that modifies the LDA Hamiltonian and
9702
energy. It is important to recall that LDA+U projectors should be
9703
quite localized functions. Otherwise the calculated populations
9704
loose their atomic character and physical meaning. Even more
9705
importantly, the interaction range can increase so much that
9706
jeopardizes the efficiency of the calculation.
9708
Two methods are currently implemented:
9711
Projectors are slightly-excited numerical atomic orbitals similar
9712
to those used as an automatic basis set by \siesta. The radii of
9713
these orbitals are controlled using the parameter
9714
\fdf{LDAU.EnergyShift} and/or the data included in the block
9715
\fdf{LDAU.Proj} (quite similar to the data block \fdf{PAO.Basis}
9716
used to specify the basis set, see below).
9719
Projectors are exact solutions of the pseudoatomic problem (and,
9720
in principle, are not strictly localized) which are cut using a
9721
Fermi function $1/\{1+\exp[(r-r_c)\omega]\}$. The values of $r_c$
9722
and $\omega$ are controlled using the parameter \fdf{LDAU.CutoffNorm}
9723
and/or the data included in the block \fdf{LDAU.Proj}.
9729
\begin{fdfentry}{LDAU.EnergyShift}[energy]<$0.05\,\mathrm{Ry}$>
9731
Energy increase used to define the localization radius of the LDA+U
9732
projectors (similar to the parameter \fdf{PAO.EnergyShift}).
9734
\note only used when \fdf{LDAU.ProjectorGenerationMethod} is \fdf*{1}.
9738
\begin{fdfentry}{LDAU.CutoffNorm}[real]<$0.9$>
9740
Parameter used to define the value of $r_c$ used in the Fermi
9741
distribution to cut the LDA+U projectors generated according to
9742
generation method 2 (see above). \fdf{LDAU.CutoffNorm} is the norm of the
9743
original pseudoatomic orbital contained inside a sphere of radius
9746
\note only used when \fdf{LDAU.ProjectorGenerationMethod} is \fdf*{2}.
9751
\begin{fdfentry}{LDAU.Proj}[block]
9753
Data block used to specify the LDA+U projectors.
9757
If \fdf{LDAU.ProjectorGenerationMethod} is \fdf*{1}, the
9758
syntax is as follows:
9760
%block LDAU.Proj # Define LDA+U projectors
9761
Fe 2 # Label, l_shells
9762
n=3 2 E 50.0 2.5 # n (opt if not using semicore levels),l,Softconf(opt)
9763
5.00 0.35 # U(eV), J(eV) for this shell
9765
0.95 # scaleFactor (opt)
9767
1.00 0.05 # U(eV), J(eV) for this shell
9768
0.00 # rc(Bohr) (if 0, automatic r_c from LDAU.EnergyShift)
9773
If \fdf{LDAU.ProjectorGenerationMethod} is \fdf*{2}, the
9774
syntax is as follows:
9776
%block LDAU.Proj # Define LDAU projectors
9777
Fe 2 # Label, l_shells
9778
n=3 2 E 50.0 2.5 # n (opt if not using semicore levels),l,Softconf(opt)
9779
5.00 0.35 # U(eV), J(eV) for this shell
9780
2.30 0.15 # rc (Bohr), \omega(Bohr) (Fermi cutoff function)
9781
0.95 # scaleFactor (opt)
9783
1.00 0.05 # U(eV), J(eV) for this shell
9784
0.00 0.00 # rc(Bohr), \omega(Bohr) (if 0 r_c from LDAU.CutoffNorm
9785
%endblock LDAU.Proj # and \omega from default value)
9789
Certain of the quantites have default values:
9792
U & \fdf*{0.0 eV} \\
9793
J & \fdf*{0.0 eV} \\
9794
$\omega$ & \fdf*{0.05 Bohr} \\
9795
Scale factor & \fdf*{1.0}
9798
$r_c$ depends on \fdf{LDAU.EnergyShift} or \fdf{LDAU.CutoffNorm}
9799
depending on the generation method.
9803
\begin{fdflogicalF}{LDAU.FirstIteration}
9805
If \fdftrue, local populations are calculated and Hubbard-like term
9806
is switch on in the first iteration. Useful if restarting a
9807
calculation reading a converged or an almost converged density
9812
\begin{fdfentry}{LDAU.ThresholdTol}[real]<$0.01$>
9814
Local populations only calculated and/or updated if the change in the
9815
density matrix elements (dDmax) is lower than \fdf{LDAU.ThresholdTol}.
9819
\begin{fdfentry}{LDAU.PopTol}[real]<$0.001$>
9821
Convergence criterium for the LDA+U local populations. In the
9822
current implementation the Hubbard-like term of the Hamiltonian is
9823
only updated (except for the last iteration) if the variations of
9824
the local populations are larger than this value.
9828
\begin{fdflogicalF}{LDAU.PotentialShift}
9830
If set to \fdftrue, the value given to the $U$ parameter in the
9831
input file is interpreted as a local potential shift. Recording the
9832
change of the local populations as a function of this potential
9833
shift, we can calculate the appropriate value of U for the system
9834
under study following the methology proposed by Cococcioni and
9835
Gironcoli in Phys. Rev. B \textbf{71}, 035105 (2005).
10040
9839
\section{RT-TDDFT}\index{RT-TDDFT}\index{TDDFT}
10041
9840
Now it is possible to perform Real-Time Time-Dependent Density