Introduction
The purpose of this tutorial is to introduce the user to the effective use of the code for three-dimensional bulk problems. This tutorial uses the example of hexagonal gallium nitride to illustrate the use of the code for bulk problems. The input is rather flexible, and can use a very basic input file, but has options that lead to easier manipulation of the calculation in increasingly more sophisticated input files. This tutorial is presented in the form of a sequence of input files that illustrate the use of various features available in the code: how to take advantage of symmetries in the system to reduce the size of the calculation, how to create Brillouin Zone sampling appropriate to the symmetry of the system, how to use scaling functions in the input to simplify the description of the atom coordinates, etc.
Units
The default units for energy are Rydberg (1 Rydberg = 0.5 Hartree = 13.605 eV), and for distances are bohr (1 bohr = 0.52918 Ang = 0.052918 nm). All input coordinates will be in these units, though one can use the scaling functions in the code to aid in the conversions.
Input files
The code uses one main input file that the user must construct. Atom files (potential+basis), which will need to be listed in the input file, can be obtained from the atom libraries. The input file is a text-rich keyword-driven document that, once constructed, is self-documenting. The input is a sequence of keyword lines with the indicated data on following lines. Only the first few (6 or 8) characters of the keyword line are significant. The remainder of a keyword line can be used for additional documentation. The data itself, in general, is free-format. The input in the setup data section is highly structured, and requires its input in a very specific order. Other sections take their data in any order. If you miss putting something in the right place in the input, not to worry. The code will very happily process the input file, and stop when it cannot find required input in the right place. It echoes the input file as it reads it, and will politely tell you what it was looking for when it runs into a problem. Interspersed among the required inputs, are a variety of usually invisible optional input sections. The code uses a variety of defaults that can be overridden in the input file. The use of many optional inputs will be illustrated in the tutorial.
Example: A basic input file
Illustrated in this first example is a very basic file for doing a bulk hexagonal GaN calculation.
Important: While a calculation with this input file will run to completion, this input file is flawed and would not be a viable production calculation because of an incorrect use of the Brillouin Zone sampling, which will be discussed in the next example.
Command options
The input file begins with a sequence of instructions to the code about the sort of calculation it will be doing. In this example, we will "do" a "setup" (the iteration independent integrals), and "do iters" (self-consistency), but "no force" and no relaxation. The "setup data" is the instruction to the code to read the data that describes the system. Please do not try to "do iters" with "no setup"!
Setup phase data
This section begins with a title line(s). The keyword "notes" signals the beginning of the title/notes section that describes the nature of the problem. The keyword "end_notes" marks the end of the title/notes section The notes keyword and title lines are optional, but highly encouraged. Since the entire input file will be echoed into the output listing file, this provides a useful output record of the nature of the calculation. In this case, we are setting up a bulk hexagonal GaN problem, with a=6.13, and a c/a ratio of 1.630714474.
The "dimension" statement is the first required data in the input. We are doing a 3-dimensional bulk problem. If we were doing a molecular problem, we would enter a "0", or if we were doing a surface problem using a 2-dimensional slab, we would have entered a "2" here.
The next required input is "primitive lattice" – supercell vectors. Here, we define a unit cell that is hexagonal in the x-y plane with a lattice parameter a=6.13 bohr (=3.24 A), and an orthogonal repeat vector in the z-direction with c=6.13*c/a.
The next required input is "grid dimensions". The code evaluates many of its integrals on a regular space (fft) grid. This regular grid is defined by taking the primitive lattice vectors, and dividing them into grid intervals. Here we use a 24×24 grid in the hexagonal plane, and put 36 grid intervals along the z-direction. The more grid intervals, the more accurately the integrals are evaluated, but, also, the more expensive the calculation. A general rule of thumb: you want spacing between points of 0.30-0.40 bohr for most systems, and 0.20-0.30 for finer accuracy in systems containing "hard" atoms. Hard atoms include late first row atoms like N, O, F, or late 3d atoms. We have nitrogen in our system, and this example has a spacing of slightly less than 0.3 bohr. This can be tuned to the problem. If I wanted to compute an energy, a coarser grid might suffice. If I wished to evaluate a force or stress, a finer grid might be required. For a hexagonal system, it is a good idea to use a grid dimension in the plane that is divisible by three, so that the various high-symmetry points in the surface (i.e.), atop, hcp and fcc hollows) are at the same position with respect to the grid.
The next required input is number of "atom types". Here, we have only two (Ga and N), and the code requires input to specify the "atom file" to be used to obtain the atomic potential and basis information for each atom type. Here, we will read from files "n.atm" and "ga.atm" for nitrogen and gallium, respectively.
Next, the input requires the number of atoms in the unit cell (4), and then a list of each "atom, type, and position". The arguments for an atom input are the atom index (the order of the atom in the list, #1-#4 here), which type (in the order in which the types were listed above), and the positions of the atoms in Cartesian coordinates (in bohr). Atom 1 is of the first atom type (i.e., from "n.atm"), and is at (3.539,0.0,0.0257), and then next atom is the second atom type (i.e., from "ga.atm"), etc.
The final required input in the setup section (for a periodic calculation) is definition of the Brillouin Zone sampling. For molecular calculations this is unnecessary – one has only real integrals, or the "gamma-point", to worry about. In a (periodic) bulk calculation, the Schroedinger equation generally needs to be solved for a sampling of points in reciprocal space. It is possible to input the k-point sampling as an explicit list of k-vectors and their weights (using a different keyword/data sequence than shown here), but this example takes advantage of an optional feature to automatically generate a simple k-point sample. The "kgrid " keyword tells the code to generate a regular mesh in reciprocal space, and here we specify a 4x4x2 grid in the BZ. A k-sample with just the gamma-point would be specified as 0x0x0. The kgrid command generates this grid, and offsets this grid from the reciprocal space origin (i.e., the gamma point) by (1/2,1/2,1/2) in the reciprocal space mesh, i.e., the body-centered position.
The final, required statement in the setup phase section is the "end setup phase" line.
Run phase data
The next section begins with the line "run phase" and ends with the line "end run phase", and allows you to modify parameters that affect the run (SCF) phase of the calculation. In this example, we modify the convergence criterion to be 0.0005. This criterion does NOT refer to the convergence of the energy in the SCF cycle, but rather refers to the maximum change in an Hamiltonian matrix element (integrals between basis functions over potentials). The value specified here will converge the total energy to roughly 1.d-6 Ry for this problem, but this will vary between different systems (e.g., a covalent system with a big energy gap will converge differently than a metallic system). Note that the convergence criterion is an optional parameter, and that the code will provide a default if you do not override it with this input.
do setup do iters no force no relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 end_notes dimension of system (0=cluster ... 3=bulk) 3 primitive lattice vectors 5.308735725 -3.0650 0.000000000 0.000000000 6.1300 0.000000000 0.000000000 0.0000 9.996279726 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 3.539157150 0.0000 0.025788046 2 2 1.769578575 3.0650 1.268818179 3 1 1.769578575 3.0650 5.023927909 4 2 3.539157150 0.0000 6.266958042 kgrid 4 4 2 end setup phase data run phase input data convergence criterion 0.000500 end run phase data
Hexagonal symmetries and BZ sampling
In the previous we mentioned that there was a problem with the input file. The problem lies in the Brillouin Zone sampling generated by the "kgrid " command and the fact that we have a hexagonal system. This command, by construction, offsets the grid of k-points to a have the gamma-point at the body-centered position, which, for a hexagonal system, is a reciprocal space grid that is not symmetric about the gamma point. For cubic and many other systems, a k-sample offset from gamma is more efficient. However, for hexagonal systems, the resulting BZ mesh is not symmetric about the gamma point. To preserve the symmetry, we must not offset the k-space grid from the gamma-point in the x-y plane, and changing the "kgrid" command to "kgrid2" accomplishes this.
In addition, note that we have expanded our input notes on this system and we have added a second line of notes. This method can be used to have up to 20 lines of notes.
do setup do iters no force no relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: hex cell needs special hex Brillouin Zone sample: "kgrid2" end_notes dimension of system (0=cluster ... 3=bulk) 3 primitive lattice vectors 5.308735725 -3.0650 0.000000000 0.000000000 6.1300 0.000000000 0.000000000 0.0000 9.996279726 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 3.539157150 0.0000 0.025788046 2 2 1.769578575 3.0650 1.268818179 3 1 1.769578575 3.0650 5.023927909 4 2 3.539157150 0.0000 6.266958042 kgrid2 4 4 2 end setup phase data run phase input data convergence criterion 0.000500 end run phase data
Picking a density functional (and spin)
By default, the code runs a spin-restricted LDA calculation, using the Perdew-Zunger parameterization of the Ceperley-Alder electron gas results (CAPZ). However, the code is capable of spin-polarized and GGA calculations, too. Add the "functional" keyword between the notes section and dimension statement, and select the desired functional. The available functionals include: CAPZ, LDA (which defaults to CAPZ), PBE, PW91, BLYP and GGA (which defaults to PBE). Note that atom files are functional-specific. However, the code does not check that atom files and calculations be the same flavor of DFT. It is possible to use PBE atoms in PW91 calculations, for example, and get good results. However, mixing functionals (atoms and calculations) is not advised, and mixing GGA and LDA is dangerous.
Spin polarization is not appropriate for this GaN example, but in the general case, to invoke spin, one specifies a spin-polarized functional. All current DFT functionals are spin-capable. The valid choices include: LDA-SP, CAPZSP, PBE-SP, PW91SP, etc. If one specifies a spin polarized functional, additional input is *required* immediately following the functional. The keyword "spin polarization" with the net spin polarization must be given. The spin polarization is in units of excess electrons of majority spin, must be non-negative (it can be zero), but need not be an integer. For example, an LDA hydrogen atom (or a molecule with a single dangling bond)
functional LDA-SP spin polarization 1.0000
will have a polarization of exactly one, i.e., one more electron of majority spin than minority spin.
do setup do iters no force no relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: change functional to PBE flavor of GGA end_notes functional PBE dimension of system (0=cluster ... 3=bulk) 3 primitive lattice vectors 5.308735725 -3.0650 0.000000000 0.000000000 6.1300 0.000000000 0.000000000 0.0000 9.996279726 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 3.539157150 0.0000 0.025788046 2 2 1.769578575 3.0650 1.268818179 3 1 1.769578575 3.0650 5.023927909 4 2 3.539157150 0.0000 6.266958042 kgrid2 4 4 2 end setup phase data run phase input data convergence criterion 0.000500 end run phase data
Symmetry I: Point groups
Our crystal has symmetries, and our calculation should take advantage of it. Around the origin defined in the atomic coordinates, we have a three-fold rotation (about the z-axis), and a number of reflection planes. The code does not (yet) find the symmetries automatically from the input coordinates, nor can you tell it the space group and have it automatically generate a crystal. You must tell it explicitly what symmetries to use. We build symmetry into our calculation using the "symops" keyword. This is the number of symmetries which we use to define the symmetry group. Then we provide the "definitions" of these symmetry elements. Symmetries are defined by seven parameters: an integer and two three-vectors. The integer defines the nature of the symmetry, the first three-vector the axis of that symmetry, and the final three-vector is a translation vector which we will ignore – and set to (0,0,0) – for now. The integer describes the order of a rotation, an inversion, or a reflection, and can take the following valid values:
2, 3, 4, 6 = an N-fold rotation -2,-3,-4,-6 = an N-fold rotation followed by inversion -1 = an inversion
Note that a reflection is equivalent to "-2", or a two-fold rotation plus inversion. The definitions in our example:
3 0.0 0.0 1.0 0.0 0.0 0.0 -2 0.0 1.0 0.0 0.0 0.0 0.0
specify a three-fold rotation around (0,0,1), i.e., the z-axis, and a rotation around the y-axis plus inversion, which is equivalent to a reflection through the x-z plane. The code will internally take these symmetry elements and generate a full symmetry group. Further examples of symmetry construction can be found on another page.
The code takes advantage of symmetry in two ways. First, it eliminates redundant symmetry-equivalent k-points in the Brillouin Zone sample, and reduces the sample to be in the irreducible BZ. This makes the size of the calculation smaller, and therefore faster. Second, if the calculation does forces and a relaxation, it will symmetrize the force calculation (in a grid-based calculation, numerical noise may cause the calculation to artificially break symmetry) and conserves symmetry in geometry updates during an energy minimization calculation/relaxation.
do setup do iters no force no relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: take advantage of trigonal symmetry about 001 axis (symops). end_notes dimension of system (0=cluster ... 3=bulk) 3 primitive lattice vectors 5.308735725 -3.0650 0.000000000 0.000000000 6.1300 0.000000000 0.000000000 0.0000 9.996279726 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 3.539157150 0.0000 0.025788046 2 2 1.769578575 3.0650 1.268818179 3 1 1.769578575 3.0650 5.023927909 4 2 3.539157150 0.0000 6.266958042 kgrid2 4 4 2 symops 2 definitions sym ops (3-fold rotate around 001, reflection plane normal 010) 3 0.0 0.0 1.0 0.0 0.0 0.0 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data convergence criterion 0.000500 end run phase data
Scaling functions I: global scaling
If we wish to modify cell parameters (e.g., to vary the lattice parameter), changing all the lattice vectors and all the atomic positions that are given in Cartesian coordinates and units of bohr is not convenient. You would need to modify every coordinate vector in the input file, a tedious process very prone to error. The code provides some simple tools to make this kind of manipulation easier. A useful one is the "scale " or "scalep" command. This command (like all the other scaling commands) must appear before dimension and after primitive lattice vectors. This command causes every coordinate in the input file to be scaled by the specified value. Here, we use the scale command to abstract the bulk lattice parameter of 6.13 from the primitive lattice vectors and all the atomic coordinates. Now, to modify the bulk lattice parameter, we need only modify the scaling value, and need not touch the various coordinate vectors in the input.
do setup do iters no force no relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: insert global scaling factor (a=6.13), and scale coordinates end_notes dimension of system (0=cluster ... 3=bulk) 3 scale 6.13 primitive lattice vectors 0.86602540378 -0.5000 0.000000000 0.00000000000 1.0000 0.000000000 0.00000000000 0.0000 1.630714474 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.57735026918 0.0000 0.004206859 2 2 0.28867513459 0.5000 0.206985021 3 1 0.28867513459 0.5000 0.819564096 4 2 0.57735026918 0.0000 1.022342258 kgrid2 4 4 2 symops 2 definitions sym ops (3-fold rotate around 001, reflection plane normal 010) 3 0.0 0.0 1.0 0.0 0.0 0.0 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data convergence criterion 0.000500 end run phase data
Scaling functions II: direction-specific
GaN is a hexagonal system, and the unit cell is defined by two parameters: the lattice parameter "a", and the c/a ratio. Using "scale " in the previous example, the lattice parameter can be modified easily. The code has a direction-specific scaling function "scalez" through which all z-coordinates can be scaled by a specified value. It comes immediately after the "scale " input, and acts in addition to (rather than exclusive of) the global scaling, i.e., z-coordinates are scaled by both the global scaling and the z-scaling factors. In this example, we abstract the c/a ratio, 1.630714474, from every z-coordinate. Now we can easily modify both the lattice parameter and c/a ratio by modifying the two scaling values.
do setup do iters no force no relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: insert global z-scaling factor (c/a), and scale z-coordinates end_notes dimension of system (0=cluster ... 3=bulk) 3 scale 6.13 scalez 1.630714474 primitive lattice vectors 0.86602540378 -0.5000 0.0000 0.00000000000 1.0000 0.0000 0.00000000000 0.0000 1.0000 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.57735026918 0.0000 0.002579764 2 2 0.28867513459 0.5000 0.126929039 3 1 0.28867513459 0.5000 0.502579764 4 2 0.57735026918 0.0000 0.626929039 kgrid2 4 4 2 symops 2 definitions sym ops (3-fold rotate around 001, reflection plane normal 010) 3 0.0 0.0 1.0 0.0 0.0 0.0 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data convergence criterion 0.000500 end run phase data
Shift coordinates
The specification of a bulk system is invariant to an overall translation. Here we shift all the atomic z-coordinates to put the first nitrogen at z=0. Further, we reorder the atoms to put the two nitrogen atoms to the front of the list and the galliums to the end of the list. Our system does not change with these modifications, but this reordering makes more apparent that there is more symmetry in this system than we have used, and this order will allow us to take better advantage of it in a later example.
do setup do iters no force no relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: translate atoms to make additional symmetry more apparent end_notes dimension of system (0=cluster ... 3=bulk) 3 scale 6.13 scalez 1.630714474 primitive lattice vectors 0.86602540378 -0.5000 0.0000 0.00000000000 1.0000 0.0000 0.00000000000 0.0000 1.0000 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.57735026918 0.0000 0.000000000 2 1 0.28867513459 0.5000 0.500000000 3 2 0.28867513459 0.5000 0.124349275 4 2 0.57735026918 0.0000 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (3-fold rotate around 001, reflection plane normal 010) 3 0.0 0.0 1.0 0.0 0.0 0.0 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data convergence criterion 0.000500 end run phase data
Symmetry II: non-symmorphic groups
We had already incorporated a three-fold axis and reflection planes in the calculation. However, this four-atom unit cell has yet more symmetry. The two nitrogen atoms are equivalent to each other as are the two gallium atoms. This is not a point group symmetry, but involves a point group symmetry operation and a translation in space. To build this symmetry into the input requires use of the second three-vector in the symmetry definition. This second vector, zero for point group operations, describes the translation necessary for symmetry specification in non-symmorphic groups. The crystal is preserved in a 60-degree rotation (i.e., a 6-fold axis) plus a translation of half the third lattice vector (z-direction). This maps the one nitrogen(gallium) atom onto the other nitrogen(gallium). Note that how the atoms were positioned in the cell in the x-y plane (remember that our bulk system is invariant to global translations) was critical to be able to specify this symmetry. The symmetry center is assumed to be at the coordinate origin. A different, translated coordinate system would have required a different symmetry specification. Further note that the original previous symmetry group is a subset of this new expanded non-symmorphic group. A 60-degree (6-fold) rotation with translation by half a z-lattice vector, applied twice, equals a 120-degree (3-fold) rotation with a full lattice vector translation, equivalent to a simple 3-fold rotation. Hence, we can simple replace the 3-fold rotation definition with the 6-fold plus translation definition, leaving only two symmetry element definitions. This doubles the size of the resulting symmetry group, reduces the size the irreducible Brillouin Zone, and enforces the equivalence of the atoms within the unit cell.
do setup do iters no force no relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: incorporate non-symmorphic (translation) symmetry: N1->2 Ga3->4 end_notes dimension of system (0=cluster ... 3=bulk) 3 scale 6.13 scalez 1.630714474 primitive lattice vectors 0.86602540378 -0.5000 0.0000 0.00000000000 1.0000 0.0000 0.00000000000 0.0000 1.0000 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.57735026918 0.0000 0.000000000 2 1 0.28867513459 0.5000 0.500000000 3 2 0.28867513459 0.5000 0.124349275 4 2 0.57735026918 0.0000 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (6-fold rotate+1/2Lz shift, reflection plane normal 010) 6 0.0 0.0 1.0 0.0 0.0 0.5 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data convergence criterion 0.000500 end run phase data
SCF parameters
In this example we show how to modify parameters that are important in the self-consistency calculations. These parameters are modified in the run phase section of the input file. The input can be in any order within this section (contrast this with the strict order of the input in the setup section). In the first example, the convergence criterion had already been specified. Here we also specify the maximum number of iterations for the self-consistency step, a "temperature" (in Ry) for level occupations (fermi-broadened, this is only really important for metallic systems, not systems with a significant gap. It aids in convergence in systems with many levels near the fermi level, i.e., metals), and an initial "blend ratio". The blend parameter can occasionally aid convergence. It sets the proportion of output potential to mix into the input potential after the first step of the self-consistency cycle (after the first cycle, a modified Broyden blending scheme updates the potential). It must take values between zero and one. There are other parameters which are also accessible from the input file, but these are the most commonly needed.
do setup do iters no force no relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: add various interesting scf options end_notes dimension of system (0=cluster ... 3=bulk) 3 scale 6.13 scalez 1.630714474 primitive lattice vectors 0.86602540378 -0.5000 0.0000 0.00000000000 1.0000 0.0000 0.00000000000 0.0000 1.0000 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.57735026918 0.0000 0.000000000 2 1 0.28867513459 0.5000 0.500000000 3 2 0.28867513459 0.5000 0.124349275 4 2 0.57735026918 0.0000 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (6-fold rotate+1/2Lz shift,reflection plane normal 010) 6 0.0 0.0 1.0 0.0 0.0 0.5 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data iterations - maximum number of scf iterations allowed 35 temperature for occupations (Ry) 0.0030 blend ratio - initial blend of input/output potential in scf 0.30 convergence criterion 0.000500 end run phase data
Atomic relaxation
Simple enough, to perform a geometry relaxation: turn on a force calculation with "do force" and use "do relax" to tell the code to do an energy minimization. The calculation will eliminate forces to find equilibrium positions of the atoms.
do setup do iters do force do relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: set commands to do a full relaxation end_notes dimension of system (0=cluster ... 3=bulk) 3 scale 6.13 scalez 1.630714474 primitive lattice vectors 0.86602540378 -0.5000 0.0000 0.00000000000 1.0000 0.0000 0.00000000000 0.0000 1.0000 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.57735026918 0.0000 0.000000000 2 1 0.28867513459 0.5000 0.500000000 3 2 0.28867513459 0.5000 0.124349275 4 2 0.57735026918 0.0000 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (6-fold rotate+1/2Lz shift,reflection plane normal 010) 6 0.0 0.0 1.0 0.0 0.0 0.5 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data iterations 35 temperature for occupations (Ry) 0.0030 blend ratio 0.30 convergence criterion 0.000500 end run phase data
Relaxation parameters
The code, by default, will relax all atoms in the problem until the maximum force is less than an internally specified force threshold (0.003 Ry/bohr). The user can modify how the relaxation is performed in the optional geometry section Of the input. The "geometry relaxation" statement begins, and the "end geometry" statement ends, the section that allows access to a variety of parameters that control the relaxation phase of the calculation. These inputs can be in any order in the geometry section. Here, "gsteps" limits the number of geometry steps to 5, "gconv" sets the force convergence criterion to 0.001 Ry/bohr, and we "only relax" the last two atoms. The code knows to relax the last two atoms in the list because a negative number "-2" was given. It would relax the first N atoms in the list if a positive integer N were given. There is a scheme, not shown here, for picking more general sets of atoms to relax or freeze in place. There is not, however, a method for for applying more sophisticated constraints on the relaxation (e.g., fixing the z-coordinate of an atom, but allowing other coordinates of that atom to relax). In this gallium nitride example, the only free parameter is the Ga-N distance, and, therefore, only the positions of a nitrogen atom with respect to a gallium atom is significant. Symmetry, and the fact that a bulk system is insensitive to a global translation allows us to fix either the positions of the nitrogen atoms or the positions of the gallium atoms without loss of generality.
I choose to relax tha galliums atoms rather than the nitrogen atoms for numerical reasons. Nitrogen is the "harder" atom, and therefore the numerical evaluations of the forces on the nitrogen atom will be less accurate. The code evaluates integrals on a regular space grid and the spatially more rapidly varying orbitals on nitrogen lead to more inaccurate grid integrals for the orbitals on the nitrogen atom. We therefore choose to fix the nitrogen and allow the atoms whose forces are more accurately evaluated to relax.
do setup do iters do force do relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: modify relaxation optimization parameters end_notes dimension of system (0=cluster ... 3=bulk) 3 scale 6.13 scalez 1.630714474 primitive lattice vectors 0.86602540378 -0.5000 0.0000 0.00000000000 1.0000 0.0000 0.00000000000 0.0000 1.0000 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.57735026918 0.0000 0.000000000 2 1 0.28867513459 0.5000 0.500000000 3 2 0.28867513459 0.5000 0.124349275 4 2 0.57735026918 0.0000 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (6-fold rotate+1/2Lz shift,reflection plane normal 010) 6 0.0 0.0 1.0 0.0 0.0 0.5 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data iterations 35 temperature for occupations (Ry) 0.0030 blend ratio 0.30 convergence criterion 0.000500 geometry relaxation gsteps - restrict the number of geometries to relax 5 gconv - max force convergence criterion 0.0010 only relax last - N are fixed by symmetries, move only Ga -2 end geometry relaxation end run phase data
Lattice coordinates
For some problems, it may be convenient to specify positions of the atoms in lattice units rather than Cartesian coordinates. This is particularly true if one desires to compute elastic constants via fitting to stress tensor results from calculations with small strains. The "coordinates" keyword, immediately after the dimension input, and before any scaling input, allows the user to switch the coordinate system used for the input of atoms. The choice of "LATTICE" specifies lattice units (as opposed to the default "CARTESIAN"). The atoms positions are constructed from the input three-vector (x,y,z) by: xA + yB + zC, where A, B, and C are the primitive lattice vectors.
do setup do iters do force do relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: transform coordinates into lattice units end_notes dimension of system (0=cluster ... 3=bulk) 3 coordinates LATTICE scale 6.13 scalez 1.630714474 primitive lattice vectors 0.86602540378 -0.5000 0.0000 0.00000000000 1.0000 0.0000 0.00000000000 0.0000 1.0000 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.66666666667 0.33333333333 0.000000000 2 1 0.33333333333 0.66666666667 0.500000000 3 2 0.33333333333 0.66666666667 0.124349275 4 2 0.66666666667 0.33333333333 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (6-fold rotate+1/2Lz shift,reflection plane normal 010) 6 0.0 0.0 1.0 0.0 0.0 0.5 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data iterations 35 temperature for occupations (Ry) 0.0030 blend ratio 0.30 convergence criterion 0.000500 geometry relaxation gsteps - restrict the number of geometries to relax 5 gconv - max force convergence criterion 0.0010 only relax last - N are fixed by symmetries, move only Ga -2 end geometry relaxation end run phase data
Calculations with strain
For use in calculation of elastic constants (from manipulation of stress-strain relations), or to examine other properties due to small perturbations for a given cell volume, there is a feature in the input that imposes a specified strain on a system. The keyword is "strain" with an optional multiplicative factor "strfac".
Notes: The "strain" keyword can only be invoked if the atomic coordinates are specified in LATTICE units. The strain matrix must be symmetric (of course), and is specified without the identity, i.e., the code interally adds the identity matrix to the specified input. The net result of the example given below is a +0.2% strain along the z-direction.
do setup do iters do force do relax setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: strain end_notes dimension of system (0=cluster ... 3=bulk) 3 coordinates LATTICE scale 6.13 scalez 1.630714474 strain 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00 strfac 0.002 primitive lattice vectors 0.86602540378 -0.5000 0.0000 0.00000000000 1.0000 0.0000 0.00000000000 0.0000 1.0000 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.66666666667 0.33333333333 0.000000000 2 1 0.33333333333 0.66666666667 0.500000000 3 2 0.33333333333 0.66666666667 0.124349275 4 2 0.66666666667 0.33333333333 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (6-fold rotate+1/2Lz shift,reflection plane normal 010) 6 0.0 0.0 1.0 0.0 0.0 0.5 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data iterations 35 temperature for occupations (Ry) 0.0030 blend ratio 0.30 convergence criterion 0.000500 geometry relaxation gsteps - restrict the number of geometries to relax 5 gconv - max force convergence criterion 0.0010 only relax last - N are fixed by symmetries, move only Ga -2 end geometry relaxation end run phase data
Cell optimization
The code as able to do cell shape optimization automatically for all LDA and GGA (and spin-polarized GGA stress is completed as of v2.57) calculations, for any dimension of periodicity. In fact, one can equilibrate at zero pressure, to a fixed non-zero pressure, to a uniaxial stress, or even a fully general stress tensor!
To invoke cell optimization, add "do cell" to the instructions at the beginning of the input. The optional cell section, beginning with "cell optimization" and ending with "end cell", allows some control over cell optimization. This is illustrated here by setting the unit cell convergence threshold (uc_conv) to a maximum stress of 0.01 GPa. To apply an external pressure, use the keyword "pressure" and give the desired pressure in units of GPa.
Important: Cell optimization requires:
- Coordinates be in LATTICE units
- no scaling functions be used (i.e. all scaling is 1.0)
(remember to install scaling directly into lattice vectors (bohr).
do setup do iters do force do relax do cell setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: turn on 3D cell optimization end_notes dimension of system (0=cluster ... 3=bulk) 3 coordinates LATTICE primitive lattice vectors 5.308735725 -3.0650 0.000000000 0.000000000 6.1300 0.000000000 0.000000000 0.0000 9.996279726 grid dimensions 24 24 36 atom types 2 atom file n.atm atom file ga.atm number of atoms in unit cell 4 atom, type, position vector 1 1 0.66666666667 0.33333333333 0.000000000 2 1 0.33333333333 0.66666666667 0.500000000 3 2 0.33333333333 0.66666666667 0.124349275 4 2 0.66666666667 0.33333333333 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (6-fold rotate+1/2Lz shift,reflection plane normal 010) 6 0.0 0.0 1.0 0.0 0.0 0.5 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data iterations 35 temperature for occupations (Ry) 0.0030 blend ratio 0.30 convergence criterion 0.000500 geometry relaxation gsteps - restrict the number of geometries to relax 5 gconv - max force convergence criterion 0.0010 only relax last - N are fixed by symmetries, move only Ga -2 end geometry relaxation cell optimization uc_conv - stress convergence (in GPa, default=0.002 ) 0.01 pressure - external pressure (GPa for 3D) 10.0 end cell end run phase data
Aliasing atom types
Sometimes keeping track of which atom is which type can be confusing, particularly in larger input files with involved configurations. The code offers the ability to use named atom types in the input of the configurations, and to alias atom files to those atom types. There are two methods to alias the input atom files. One is explicit, as with the the N atom below. You provide the alias name and set it equal to the atom file: N = n.atm. Alternatively, the code constructs a default from the atom file name, by using the file name prefix. From /library/lda/Ga.atm it abstracts the alias "Ga". Note that the code ignores the directory specification. To use the alias in the input, simply replace the type numbers with the type alias. Note: this is all or nothing. Use type names only, or use type numbers only.
do setup do iters do force do relax do cell setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: aliasing atom types end_notes dimension of system (0=cluster ... 3=bulk) 3 coordinates LATTICE primitive lattice vectors 5.308735725 -3.0650 0.000000000 0.000000000 6.1300 0.000000000 0.000000000 0.0000 9.996279726 grid dimensions 24 24 36 atom types 2 atom file N = n.atm atom file /library/lda/Ga.atm number of atoms in unit cell 4 atom, type, position vector 1 N 0.66666666667 0.33333333333 0.000000000 2 N 0.33333333333 0.66666666667 0.500000000 3 Ga 0.33333333333 0.66666666667 0.124349275 4 Ga 0.66666666667 0.33333333333 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (6-fold rotate+1/2Lz shift,reflection plane normal 010) 6 0.0 0.0 1.0 0.0 0.0 0.5 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data iterations 35 temperature for occupations (Ry) 0.0030 blend ratio 0.30 convergence criterion 0.000500 geometry relaxation gsteps - restrict the number of geometries to relax 5 gconv - max force convergence criterion 0.0010 only relax last - N are fixed by symmetries, move only Ga -2 end geometry relaxation cell optimization uc_conv - stress convergence (in GPa, default=0.002 ) 0.01 pressure - external pressure (GPa for 3D) 10.0 end cell end run phase data
Labeling atom names
You can also name the atoms themselves. There are two methods to label the input atoms. You can use a numerical sequence, as above, or you can use character strings as with the atom below. The atom names can be very useful for identifying which atoms are special, to follow the positions of atoms in the course of a relaxation with a simple grep, or to use in selecting atom atoms in specifying constraints for geometry updates. To use the labels in the input, simply replace the atom numbers with the atom labels. Note: this is all or nothing. Use atom names only, or use atom numbers only.
do setup do iters do force do relax do cell setup data notes GaN bulk wurzite: a=6.13 bohr, c/a=1.630714474, GaN(Z)=3.755109729 Example: aliasing atom types end_notes dimension of system (0=cluster ... 3=bulk) 3 coordinates LATTICE primitive lattice vectors 5.308735725 -3.0650 0.000000000 0.000000000 6.1300 0.000000000 0.000000000 0.0000 9.996279726 grid dimensions 24 24 36 atom types 2 atom file N = n.atm atom file /library/lda/Ga.atm number of atoms in unit cell 4 atom, type, position vector Atom01 N 0.66666666667 0.33333333333 0.000000000 Atom02 N 0.33333333333 0.66666666667 0.500000000 Atom03 Ga 0.33333333333 0.66666666667 0.124349275 Atom04 Ga 0.66666666667 0.33333333333 0.624349275 kgrid2 4 4 2 symops 2 definitions sym ops (6-fold rotate+1/2Lz shift,reflection plane normal 010) 6 0.0 0.0 1.0 0.0 0.0 0.5 -2 0.0 1.0 0.0 0.0 0.0 0.0 end setup phase data run phase input data iterations 35 temperature for occupations (Ry) 0.0030 blend ratio 0.30 convergence criterion 0.000500 geometry relaxation gsteps - restrict the number of geometries to relax 5 gconv - max force convergence criterion 0.0010 only relax last - N are fixed by symmetries, move only Ga -2 end geometry relaxation cell optimization uc_conv - stress convergence (in GPa, default=0.002 ) 0.01 pressure - external pressure (GPa for 3D) 10.0 end cell end run phase data