Basic paratec usage

Entering the crystal structure

The crystal structure is specified by the lattice vectors and the atomic coordinates. The lattice vectors are given as rows of coordinates with respect to a cartesian basis. The unit lenght is an atomic length unit, i.e. the Bohr radius.

In addition to the lattice vectors, the volume may be specified. If you provide a volume, the lattice vectors will be rescaled such that they span a unit cell with the given volume. If you leave out the volume statement, the lattice vectors are assumed to be given in atomic units.

The lattice vectors and optionally the volume must be enclosed within a begin latticevecs and end latticevecs statement.

Example:

begin latticevecs
coord       7.22362499411784  -5.04013632834720   0.00000000000000
coord       0.00000000000000   3.95068558900000   2.76187238600000
coord       0.00000000000000  -3.95068558900000   2.76187238600000
volume 157.3848
end latticevecs

The coordinates of the atoms must be enclosed by the statements begin coordinates and end coordinates. Each new element is indicated with a keyword newtype. All coord statements following a newtype statement refer to the new element, until the next newtype keyword is encountered.

The newtype keyword is followed by the atomic type, e.g. O,Si,C, by the spin polarization of the initial charge configuration for this atom, and the size of the muffin-tin sphere (used only for projected DOS, and, if unspecified, is set to 1.0 a.u by default). For instance

newtype Si 0.5 2.4

will define a Si atom. For setting up the initial charge density, a definition of spin polarization of $(\rho_{up} - \rho_{down})/(\rho_{up} +\rho_{down}) = 0.5$ will be used, with a default of 0.2. The muffin-tin sphere around the silicon atom will have a radius of 2.4 a.u. For non-spin-polarized calculations, set the spin polarization to 0.

By default, the coordinates following coord are understood to be relative to the basis vectors. To facilitate the input of molecular crystals, it is possible to give them in absolute, cartesian coordinates (a.u.) by specifying the keyword coordinates_absolute OUTSIDE the begin/end of the coordinates.

NOTE: This is a change from pre-5.0 versions of paratec! Previously one specified coordinates absolute, whereas now an underscore is required (i.e. coordinates_absolute)!

Example:

begin coordinates
newtype O 0.17
coord       0.15693515545994   0.10044534718500  -0.10044534718500
coord      -0.15693515545994  -0.10044534718500   0.10044534718500
newtype  Si 0.0
coord       0.5  0.5  0.5
end coordinates

The job variable

The job variable in the input file determines, what the code does:

job relax                   Does a structural relaxation job

relax_from_checkpoint       Continues relaxation from CHECKPOINT file.
                            This helps in case of a crash.

job relax_recycle_H         Relaxes, and uses Hessian matrix from
                            the CHECKPOINT file.  This improves efficiency
                            if you perform several relaxations on
                            similar systems.

job MD                      Does a molecular dynamics run

Default: job relax

If you don't want to relax, you simply set relax_max_iter to be 1.

NOTE: If you specify job relax_from_checkpoint and there is no CHECKPOINT available, paratec will abort without performing any calculation.

Deciding how much to relax

If the square modulus of the gradient (composed of forces and stress) falls below relax_accuracy, the relaxation stops.

Default: relax_accuracy 0.01

Don't choose it to be less then 0.001, since this is generally close to the accuracy with which the forces and stress is computed.

To avoid runaway jobs, the variable relax_max_iter determines the maximum number of relaxation steps to be performed. The statement

relax_max_iter 10

instructs the code to do at most 10 relaxation steps.

Default: relax_max_iter 0

Basic molecular dynamics variables

MD_time_step  default = 3.0        Time step for MD in femtoseconds

MD_temp       default = 500.0       Temperature for MD run in Kelvin

MD_max_iter   default = 100         Number of time steps in MD run

Plane wave code parameters

There are several parameters for the plane wave code which you will need to set correctly for each calculation.

Plane wave job

Whenever the plane wave code is called for a given structure, it can perform several jobs. For a simple relaxation, all you will need to do is an scf (self-consistent field) loop:

begin pw_jobs
pw_job scf
end pw_jobs

Later, you may want to add extra pw_jobs to do things like compute bandstructures and potential plots, at which time you should read Section 5.1.1.

Cutoff Energies

The following two variables for the energy cutoffs are mandatory, and have no default value. All energies are in Rydbergs.

energy_cutoff             Cutoff for the Hamiltonian matrix
submatrix_energy_cutoff   Cutoff for the submatrix to be used for the
                          generation of the starting guesses.

Example:

energy_cutoff           70
submatrix_energy_cutoff 10

Your energy_cutoff will depend on your pseudopotential, and generally you can simply use a submatrix_energy_cutoff of 5 to 10 Ry with no problem.

Advice: take the submatrix cutoff about 1/5 of the energy_cutoff, but usually not larger than about 10 Rydbergs. The generation of the starting guess can become very time consuming if the submatrix is too large. Also, the memory requirements increase rapidly with submatrix_energy_cutoff, because the full submatrix is kept in memory. The amount of time spent should be equal to the time of about 3 CG iterations on the wavefunctions.

Very small values can cause problems. In order to prevent high symmetry local minima, one can add randomization to Fourier coeffficents larger than the submatrix cutoff energy by

randomize_diag_start_guess 1.0  Default: 0

K-point grid and gaussian smearing

For each crystal you will need to specify what k-point grid you want to use to sample the Brillouin zone.

Example:

k_grid  4 4 4
k_grid_shift 0.5 0.5 0.5

This specifies a Monkhorst-Pack grid of 4x4x4 points shifted off $\Gamma$ by $(\frac18,\frac18,\frac18)$ in reciprocal lattice vector units. This example is probably sufficient for most semiconductors with a smallish unit cell. If you have a large unit cell, you need fewer k-points, and if you have a metal, you will need a lot more.

If you are calculating a metal, you should use gaussian smearing. This smears out the fermi level, so that more of your k-points will be effectively on the fermi surface. The gaussian_smearing is specified in eV. If you are calculating an insulator, you should set gaussian_smearing to be some small value.

Example:

gaussian_smearing 0.05

Number of bands

You need to specify the number of bands to be computed. Of course, this should be greater than or equal to the number of actually filled bands.

Example:

number_bands 4

For insulators one should be equal to the number of filled bands. For metals, this should be more and is dependent on the system and the smearing method used. For the Grassmann_metal minimization method, 4 bands are added to this value.

Pseudopotential (PP) format

Starting with Paratec 5.1.1 has three different formats for PP input. For all formats the file must reside in chemical_symbol_POT.DAT files. The format is chosen by the line

pp_format         1 - original format in Fourier space
                  2 - L-W Wang's version of Martins code
                  3 - FHI98PP code

Format 1 requires a binary file format (created with the kbascbin executable created by "make tools"). Formats 2 and 3 are in asci format. Format 1 is incompatible with the pulay_tf mixing option. Formats 2 and 3 are in real space and include the atomic pseudowavefunctions. Format 3 is the suggested format for new users. Information and source can be found at http://www.fhi-berlin.mpg.de/th/fhi98md/fhi98PP.

For format 3, this additional keyword information must be given:

begin pseudopotential
pp_data  (l of local PP) + 1, occupation of s shell, p shell, d shell
                  for every atom type in the order they appear
end pseudopotential

e.g. For GaAs (Ga appearing as the first newtype in the coordinates and the s potential as local),

begin pseudopotential
pp_data   1  2.0 1.0 0.0
pp_data   1  2.0 3.0 0.0
end pseudopotential

Congratulations!

You now know (unless I left something out) the everything you need to perform simple calculations with paratec.

You may want to read the sections on diagonalization methods 5.3.1 and mixing methods 5.3.4, to learn how you may be able to make your calculations run faster.