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 filesdo 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)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).
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.
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.
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.
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.0×30.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.
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.
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 GA in the usual way. The code then generates a regular grid in the parallelpiped defined by the reprical lattice vectors: GA/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.
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 axisNB: 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 vectorNB: 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.