SeqQuest input: Setup phase data
Table of Contents
- Overview
- Input options
- More detailed explanations
Overview
This page gives a description of the SeqQuest "setup phase" input section.
It sets up atom positions and types, cell vectors, symmetries, Brillouin
Zone sampling, etc., that define a QM calculation.
It attempts to handle the many interdependencies among the various input
options (e.g., one needs a k-point sample for a bulk calculation, but not
for molecule), and therefore is highly structured.
The input is driven by (left-justified) keywords in a very strict order,
though most of the data associated with those keywords is free-format.
There are optional keywords/data, but they must appear in specific places
within the input file. This section is mandatory for a Quest input file.
The default energy unit is Rydberg (Ry), and default distance unit is bohr.
Input options
The setup phase section is mandatory for SeqQuest.
The setup phase input section appears directly after the command options
section, begins with the keyword "setup data" and ends with the
keyword "end setup phase data".
The following are the common input keywords in this section, and must
appear in the exact order indicated, if they appear.
There are many mandatory keyword/data inputs.
Optional keywords are enclosed in square brackets, e.g.
[scale ].
- setup data - begin setup phase data input section
- [titleN] - not required, but highly recommended
- N title lines (up to 9 lines, as given by "N")
- [functional] - flavor of DFT functional;
default is LDA
- fcnal (e.g., LDA, LDA-SP, PBE, PBE-SP, PW91, BLYP, AM05 ...)
- [spin polarization] - required iff spin-polarized fcnal chosen
- spin_pol (real, units of excess electrons of majority spin)
- [vdw_potential] - type of C6-based van der Waals correction (either Grimme D2 or Kim/Goddard ULG)
- vdwtype (where vdwtype=[OFF|DFT-D2|D2|ULG12|ULG] )
- [efield] - external applied electric field vector (Ry/au)
- efield_x efield_y efield_z (Ry/bohr)
- [dielectric constant] - bulk material static dielectric constant
- dielec_constant
- dimension of system (cluster=0, ... ,3=bulk)
- ndim (int)
- [coordinate units] - units for atomic coordinate vector input
- coord_units (LATTICE or CARTESIAN)
- [scale ] - global scale factor for system
- scalep (real)
- [scaleu] - scale factor, periodic directions
- scaleu (real)
- [scalex] - scale factor, x-direction
- scalex (real)
- [scaley] - scale factor, y-direction
- scaley (real)
- [scalez] - scale factor, z-direction
(e.g. c/a ratio)
- scalez (real)
- [strain] - 3x3 deformation matrix (no identity) to apply to cell
- strain(1,1) strain(2,1) strain(3,1)
strain(1,2) strain(2,2) strain(3,2)
strain(1,3) strain(2,3) strain(3,3) (3*real)
- [[strfac]] - optional scale factor for strain (only if have strain)
- strfac (real)
- primitive lattice vectors - supercell vectors (three)
- A_1 A_2 A_3
B_1 B_2 B_3
C_1 C_2 C_3 (3*real)
- grid dimensions - # intervals along lattice vectors
- n_grid_A n_grid_B n_grid_C (3*int)
- atom types - number of atom types (atom files to be read in)
- n_types (int)
- do i_type = 1, n_types
- atom file - location of atom file
(optionally, an alias, too)
- atomfile | alias_type=atomfile (file/alias=char strings)
- enddo
- [vdw_data] - only available if vdw_potential is selected above
- ... modify/configure C6 vdW parameters
- [masses] - atomic masses for each type (atomic units, C=12.0)
- do i_type = 1, n_types
- atom_mass (real)
- enddo
- [energies] - reference energy for each atom type (Rydberg), in order of atom files
- do i_type = 1, n_types
- energy_ref (real)
- enddo
- [ionopt] expert only, see charge details
- ionopt (0=jellium, 1=iffy LMCC/sphere, 2=full LMCC/guassian, 3=dble-layer b.c.)
- [charge] - net charge on the system (in atomic units)
- elchrg (real)
- [location of charge] expert only, see
charge details
- x_charge y_charge z_charge (3*real)
- number of atoms
- n_atom
- [geomfile] - (2.66) read atom/type/coordinates from file "lcao.geom_in"
- atom, types, and positions? three formats currently supported:
- do i_atom = 1, n_atom
- i_atom i_type x_atom y_atom z_atom (int,int,3*real)
- enddo
- -- OR --
- do i_atom = 1, n_atom
- i_atom alias_type x_atom y_atom z_atom (int,char,3*real)
- enddo
- -- OR --
- do i_atom = 1, n_atom
- label_atm alias_type x_atom y_atom z_atom (char,char,3*real)
- enddo
- [origin offset] expert only - coordinate origin shift wrt supercell
- x_shift y_shift z_shift (3*real)
- [wigner seitz origin] expert only
- see charge details
- x_wigner y_wigner z_wigner (3*real)
- [symmetry center] expert only - location of high-symmetry point
- x_symctr y_symctr z_symctr (3*real)
- if( ndim > 0 )then - periodic system
- [kgrid | kgrid0 | kgrid1 | kgrid2 | kgrid3 | kgridh]
- k-point sampling
- n_kgrid_A n_kgrid_B n_kgrid_C (3*int)
- endif
- [symops ] - number of symmetry definitions
- n_symops (int)
- [definitions of symmetry operations] - iff symops invoked
- do i_symop = 1, n_symops
- i_symtype x_symaxis y_symaxis z_symaxis
a_off b_off c_off (int, 6*real)
- enddo
- end setup phase data - end of setup phase input section
- [run phase input data]
- begin run phase data (optional)
Return to Top
Details of the setup data input
Density functionals
| non-spin | spin-polarized |
CAPZ LDA | CAPZ LDA | CAPZSP LDA-SP |
PBE GGA | PBE GGA | PBE-SP GGA-SP |
PW91 | PW91 | PW91SP |
BLYP | BLYP | BLYPSP |
AM05 | AM05 | AM05SP |
SeqQuest has both LDA and GGA flavors available. All current flavors
of DFT are available as both spin-polarized or not, and the spin and
non-spin flavors of a functional are chosen with different keywords.
The spin-polarized variant of the functional is triggered by the "SP".
This then requires input of the net spin polarization immediately
following the functional specification.
The units of the spin polarization are net excess electrons
of majority spin. Note: the convention for specifying a spin polarized
calculation may change in a future version.
The default (only) LDA functional is the Perder/Zunger parameterization
J. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981)
of the Ceperley/Alder electron gas results
D.M. Ceperley and B.J. Alder, Phys. Rev. Lett. 45, 566 (1980).
The default GGA is PBE
J.P. Perdew, K. Burke, and M. Ernzerhof,
Phys. Rev. Lett. 77, 3865 (1996); 78, 1396(E) (1997)
Also available are PW91,
J. P. Perdew, in Electronic Structure of Solids '91,
Ed. P. Ziesche and H. Eschrig (Akademie Verlag, Berlin, 1991) 11.
J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson and
C. Fiolhais,
Phys. Rev. B 46, 6671 (1992); 48, 4978(E) (1993).
BLYP, using Becke exchange
A.D. Becke, Phys. Rev. A 38, 3098 (1988)
with Lee/Yang/Parr correlation
C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 37, 785 (1988)
and AM05:
R. Armiento and A. E. Mattsson, Phys. Rev. B 72, 085108 (2005).
Return to Top
Dimension of system
SeqQuest is capable of calculations of finite molecule/cluster (ndim=0),
as well as periodic 1D chain (ndim=1), 2D slab (ndim=2), and bulk 3D (ndim=3)
calculations. Note that a molecule/cluster (ndim=0) calculation of a
Na-Cl molecule, for example, is formally different than a bulk molecular
crystal (ndim=3) calculation of Na-Cl, even if one uses the same primitive
lattice vectors!
The choice of dimension has implications throughout the input and calculation:
- Important: The scaling inputs do not operate on the
non-periodic primitive lattice vectors. E.g., the C lattice vector
would be unaffected by scaling functions in a slab (ndim=2) calculation.
- The "scaleu" operates on atomic coordinates only in periodic directions.
- There is no k-point sampling for a molecule (it must be gamma-point),
but there is for a periodic calculation.
- Integrals and densities are restricted from crossing non-periodic bounds.
- Electrostatic boundary conditions are applied differently at boundaries
depending on dimension: periodic v. non-periodic.
Important: The code has very definite ideas about which
directions are periodic and which go off to vacuum in chain and slab
calculations (the point is moot in molecular or bulk calculations).
For a slab (ndim=2) calculation, the periodic directions must be in
the x-y plane and the third, non-periodic, lattice vector
must be along the z-axis.
For 1D chain (ndim=1) calculation, the periodic directions must be along
the x-axis, and the non-periodic lattice vectors must be
normal to the x-axis, i.e., must be in the y-z plane.
Return to Top
Scale factors
The purpose of the scaling factors is to easily modify the length
dimensions without having to edit every single coordinate (lattice vectors,
atomic position vectors, k-points) in the input file. The scaling functions
act to scale every common coordinate in the problem. The "scale" keyword
scales everything, the scalex/scaley/scalez keywords scale the individual
coordinates (which can aid in such things as altering the c/a ratio for
a crystal input), and the "scaleu" function scales the periodic directions
only (e.g. the x-y coordinates for an ndim=2 slab calculation).
Important: the scale functions do not scale lattice vectors in
non-periodic directions, only along periodic directions. The
atomic positions may be scaled in the z-direction for a slab
calculation using the "scale " or 'scalez" keyword, but the third
lattice vector "C", the non-periodic lattice vector along the
z-axis is not scaled.
Return to Top
Lattice vectors
The code is a supercell code, and the primitive lattice vectors define
the supercell. Hence, even when doing a ndim=0 molecular calculation,
all three primitive lattice vectors must be specified.
The lattice vectors along periodic directions will be determined by
the physics of the problem, but the lattice vectors in non-periodic
directions are determined by computational considerations.
First, note discussion above about dimension of system. For a ndim=2 slab
calculation, the slab normal must be along the z-axis,
and so must the third lattice vector "C".
For a ndim=1 chain, the periodic direction
(given by the first lattice vector "A") must be the x-direction,
and the last two lattice vectors, B and C, must be in the y-z plane.
The extent of the non-periodic lattice vectors must be large enough to
contain the charge density of the system. A good rule of thumb is
that charge density extends 10-12 bohr around an atom. Hence, lattice
vectors should be set up so that the no atom is closer than 10-12 bohr
to a non-periodic boundary.
For example, in a single-layer slab calculation (all atoms at z=0),
the supercell should extend 10-12 bohr in both +z and -z,
requiring a supercell of total length 20-24 bohr in the z-direction.
In a multi-layer slab calculation, with atoms ranging from z=-4 (bohr)
to z=+4, the supercell needs to stretch from [-14,14] or
a total of 28 bohr.
Important: the coordinate origin for the atomic input is,
by default, positioned at the center of the supercell (strictly speaking,
this is for non-periodic directions only). For example, in the slab
calculation mentioned above, if the atom positions ranged from z=0
to z=+10, one would need a supercell that ranged from z=-10
to z=+20 to contain the density, a width of 30 bohr.
However, with z=0 defined as the center of the supercell, one would
need a supercell to stretch from z=-20 to z=+20,
a width of 40 bohr.
The supercell origin would need to be shifted to center the slab in the
supercell to be able to take advantage of using the smaller supercell.
Return to Top
Grid dimensions
The regular space (fft) grid on which many of integrals in Quest are
performed is defined by primitive grid vectors: the primitive lattice
vectors divided by the respective grid dimensions. Hence the grid
directions are aligned with the cell directions.
Choice of the dimensions of the grid is guided by a compromise
of two considerations. Greater accuracy in the numerical
quadratures requires larger grid dimensions (more grid points),
but more grid points means a more expensive (in both memory and cpu)
calculation. Experience has indicated that most calculations are
converged (assuming orthogonal vectors) with a grid interval of
0.30-0.35 bohr along each of the lattice vectors. With a cubic
30.0x30.0x30.0 box, for example, grid dimensions of 100x100x100 would
normally be well converged.
For faster calculations where cruder accuracy is sufficient, a spacing
of 0.4-0.5 bohr can work.
For very fine accuracy, and for problems containing bad actors like
oxygen, a very "hard" atom, a finer grid might be necessary to improve
the energy calculation or improve the accuracy of a force calculation
(for a geometry calculation).
In that event, grids with spacings as small as 0.25 or even 0.20 can be
useful.
Return to Top
Charge state calculations and multipole moments in supercells
In a supercell code, special consideration must be given to the
treatment of supercell effects in systems that have a net charge,
or have a significant dipole or higher multipole moment.
SeqQuest corrects the Coulomb potential (not just the energy) for the
contamination of the local potential by the artificial periodic images
introduced in the supercell approximation, reconstructing the local
potential correct for an isolated charge or dipole.
SeqQuest uses the "Local Moment CounterCharge" (LMCC) method for
correcting the supercell electrostatics. For systems with a vacuum
gap (molecules, or slabs), SeqQuest automatically applies a "vacuum LMCC":
"Local electrostatic moments and periodic boundary conditions",
P.A. Schultz, Phys. Rev. B 60, 1551 (1999).
For treating a defect charge (or dipole) in a periodic system,
a "defect-LMCC":
"Charged local defects in extended systems",
P.A. Schultz, Phys. Rev. Lett. 84, 1942 (2000),
can be used. For an illustration of an application of this method
to computing defect levels in silicon, see the following:
"Theory of defect levels and the 'band gap problem' in silicon"
P.A. Schultz, Phys. Rev. Lett. 96, 246401 (2006).
In a periodic slab (ndim=2) with a planar (normal) dipole,
or a molecule/cluster (ndim=0) calculation with a net charge or dipole,
the LMCC is invoked automatically (without any input on the part of the
user), and supercell effects up to a dipole are treated exactly.
Supercell errors from quadrupole and higher moments are not removed,
but those errors scale as L^{-5} or smaller, where L is the
linear dimension of the supercell.
Use of SeqQuest on molecules with a net charge or dipole, or a polar
slab calculation (with a planar dipole) should cite the first paper above.
For a charged (local) defect in a periodic system (e.g.,
a +1 state for the nitrogen substitutional defect in bulk silicon),
the situation is more complicated because of the need to make
workable boundary conditions. The defect-LMCC method is then called
for, but its use is very involved, and should be done only by those
who are expert in the use of SeqQuest and only after detailed
consultations. This is a very new, and not completely researched
method for treating boundary conditions for defects. Anyone attempting
to use this method for a charged bulk defect should read and understand,
in detail, the defect-LMCC paper. The keyword "ionopt" tells the code
which method to use. If the defect-LMCC method is invoked (ionopt=-2)
for a charged defect, then a presumptive location for the charge must
be specified (using the "rchrg " keyword) so that the code knows where
to locate its LMCC. In addition, to extract useful total energies,
an energy reference (i.e., electron chemical potential) must be
established. There are multiple options a user has to try and do
this, and all are confusing and error-prone. Any calculation of
charged defects should be considered highly "researchy".
WARNINGS about the LMCC:
- (1) In a periodic bulk defect calculation, only use the fixed
location (ionopt=-2) version of the LMCC.
- Computation of the dipole necessary to dynamically locate the
compensation charge is faulty in a bulk calculation.
- (2) Restrict the LMCC location to the origin (0,0,0).
- The latest version of the code (2.59) contains a bug so that
the calculation will fail with a non-origin LMCC location.
This will be very obvious - the SCF cycle will fail to converge.
Return to Top
Brillouin Zone (k-point) sampling
Brillouin Zone (BZ), or k-point, sampling is a necessary burden of
periodic calculations.
For molecular calculations (ndim=0), BZ sampling is moot, as assumed
to be the gamma point.
For periodic calculations (ndim=1,2,3),
a BZ sample is required, and the kgridZ keywords need to be invoked.
SeqQuest has only a very primitive scheme for building a BZ sample.
In the following discussion, we assume 3D bulk (ndim=3), but the
1D and 2D cases are handled analogously. Taking {A,B,C} as
the primitive lattice (supercell) vectors, we get the reciprocal
lattice vectors G_{A} in the usual way.
The code then generates a regular grid in the parallelpiped defined by
the reprical lattice vectors: G_{A}/n_grid_A, etc.
SeqQuest, by default, shifts this BZ grid off of the gamma-point
to the body-centered location (i.e., the "kgrid " keyword).
This offset is usually a more efficient (convergent) k-sample for a given
number of points than centering a k-grid at the gamma point.
For an orthorhombic cell (e.g. cubic) this is true,
but for a hexagonal cell, the resulting k-grid does not have the needed
symmetries about the body-center point in the grid.
The offset k-sample in the hexagonal plane is anisotropic, introducing
an artificial symmetry-breaking to the system.
To make the sampling symmetric for a hexagonal system (hexagonal in
the x-y plane!), we must use either
"kgrid2" or "kgridh" so that the k-grid is not shifted off
gamma in the x-y plane, and is shifted only in the z-direction.
Use of "kgrid3" uses a gamma-centered grid, i.e., without any shift.
SeqQuest automatically uses any symmetry specified in the input
to reduce the generated BZ sample to k-points in the irreducible BZ.
To get a gamma point sample for a periodic calculation, enter
"0 0 0" as input to the "kgrid " keyword. The "0" input
causes a single interval along that reciprocal lattice vector,
with the k-point coordinate at gamma in that direction.
For a 2D (ndim=2) calculation, n_kgrid_C must be zero (i.e. gamma
point in the non-periodic z-direction).
For a 1D calculation, both n_kgrid_C and n_kgrid_B
must be entered as zero.
For bulk defect calculations (and closed-shell molecular calculations),
a "discrete defect occupation" (DDO) scheme is available to replace the
conventional fermi-dirac method for occupying states. DDO enforces
closed shell (T=0K) occupancies, and, in bulk, prevents metallic
occupations. See PRL 96, 246401 (2006). The DDO scheme is invoked
by using the closed option in the
Run phase Section of the input file.
Return to Top
Symmetry
Symmetry in Quest is used to:
- reduce the BZ sample to points in the irreducible Brillouin
Zone (reducing the cost of the calculation),
- symmetrize the resulting density,
- enforce symmetry constraints on atomic positions during relaxation,
- and enforce symmetry on cell vectors during cell optimizations.
The code does not deduce the symmetry of the system, nor does it
know anything about space/point groups; the desired symmetries
must be specified in the input file by the user.
The user need not specify all the symmetry operations for a given
symmetry group in the input file, just those that are sufficient
to define a group (i.e., a minimum basis that expands to a group).
The code will automatically expand the entered elements into the full group.
The symmetry operations have a general input: a type of symmetry
(expressed as a rotation with a possible inversion) around a
symmetry axis, with a possible translation (to allow for
non-symorphic groups):
- A symmetry operation is defined using the following:
- i_symtype x_symaxis y_symaxis z_symaxis a_off b_off c_off
- i_symtype - Defines symmetry type.
- Allowed values: 1,2,3,4,6, -1,-2,-3,-4,-6
|i_symtype | = order of rotation (2-fold,3-fold, etc)
i_symtype < 0 = rotation plus inversion
i_symtype = -1 is simple inversion
i_symtype = 1 or -2 is a reflection (2-fold+inversion=reflection)
- x_symaxis y_symaxis z_symaxis - Defines rotation axis
- NB: Input is Cartesian x,y,z (even if coordinate=LATTICE)
Axis need not be normalized, but must be non-zero (except for inversion)
- a_off b_off c_off - Defines translation/offset vector
- NB: Input is in Lattice units (even if coordinate=CARTESIAN)
This convention allows for a very natural (and very compact) scheme
for describing any general point group or space group. For example,
the bulk cubic symmetry, a group with 48 total symmetry operations,
can be described with four input symmetry definitions:
a 4-fold rotation about the z-axis,
a 3-fold rotation about the 111-axis,
a 2-fold rotation about the 110-axis,
and an inversion:
symops - for cubic symmetry
4
definitions of symmetry operations
4 0.0 0.0 1.0 0.0 0.0 0.0
3 1.0 1.0 1.0 0.0 0.0 0.0
2 1.0 1.0 0.0 0.0 0.0 0.0
-1 0.0 0.0 0.0 0.0 0.0 0.0
These four symmetry operations are sufficient to define the group; SeqQuest
will internally expand these four input elements to generate the full group
of 48 symmetry elements.
Note that the inversion has a zero length axis. This is a special case, as
axes for all other symmetry types must be non-zero. The cubic group needs
no translations (note translation vectors are all zero) in its definitions,
but other (non-symorphic) symmetry groups will have translation.
To see examples of how to define the symmetries for a variety of different
crystal symmetries,
see symmetry input examples.
Return to Top
or Input Manual
or User Guides
or SeqQuest Home
Send questions and comments to:
Peter Schultz
at
paschul@sandia.gov
Last updated: August 22, 2016