Paratec

http://www.nersc.gov/projects/paratec5.1 Documentation

Bernd Pfrommerhttp://civet.berkeley.edu/pfrommer/index.html David Roundyhttp://civet.berkeley.edu/droundy/index.html Young-Gui Yoonhttp://civet.berkeley.edu/yyoon/index.html

David Raczkowski http://www.nersc.gov/ davidra

May 10, 2001

2

Contents

1 How to get and run paratec 5 1.1 Where to find paratec . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 How to run paratec . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.2.1 Input files . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Output files . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.3 How to compile and run on different machines . . . . . . 6 1.2.4 Some tips on porting paratec . . . . . . . . . . . . . . . . 6 1.2.5 Some performance numbers on various machines . . . . . 9 1.2.6 How to change paratec . . . . . . . . . . . . . . . . . . . . 9

1.3 Changes in Paratec 5.0 . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Basic paratec usage 11 2.1 Entering the crystal structure . . . . . . . . . . . . . . . . . . . . 11 2.2 Relaxing the crystal . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.2.1 The job variable . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.2 Deciding how much to relax . . . . . . . . . . . . . . . . . 12

2.3 Plane wave code parameters . . . . . . . . . . . . . . . . . . . . . 13 2.3.1 Plane wave job . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.2 Cutoff Energies . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.3 K-point grid and gaussian smearing . . . . . . . . . . . . 14 2.3.4 Number of bands . . . . . . . . . . . . . . . . . . . . . . . 14 2.3.5 Pseudopotential format . . . . . . . . . . . . . . . . . . . 14

2.4 Congratulations! . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 Advanced relaxation parameters 17

3.1 Specifying the crystal structure . . . . . . . . . . . . . . . . . . . 17 3.1.1 Inversion and symmetry . . . . . . . . . . . . . . . . . . . 17 3.1.2 Specifying an initial strain . . . . . . . . . . . . . . . . . . 17

3.2 Relaxing a structure . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2.1 Controlling the relaxation process . . . . . . . . . . . . . 18 3.2.2 Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.2.3 Constraining the relaxation . . . . . . . . . . . . . . . . . 19 3.2.4 Initialization of the inverse of the Hessian matrix . . . . . 20 4 Advanced plane wave code parameters 21

4.0.5 Plane wave code jobs . . . . . . . . . . . . . . . . . . . . . 21 4.0.6 Exchange correlation functional . . . . . . . . . . . . . . . 22 4.0.7 Diagonalization Methods . . . . . . . . . . . . . . . . . . 22 4.0.8 Diagonalization . . . . . . . . . . . . . . . . . . . . . . . . 23 4.0.9 Kpoint grid . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.0.10 Number of bands . . . . . . . . . . . . . . . . . . . . . . . 24 4.0.11 Occupation numbers and gaussian smearing . . . . . . . . 24 4.0.12 Self-consistent mixing methods . . . . . . . . . . . . . . . 25 4.0.13 Self-consistent iteration . . . . . . . . . . . . . . . . . . . 26 4.0.14 Options for improved stress . . . . . . . . . . . . . . . . . 27

5 Getting extra outputs and other interesting stuff 29

5.0.15 External potential . . . . . . . . . . . . . . . . . . . . . . 29 5.1 NMR shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

5.1.1 Additional Files . . . . . . . . . . . . . . . . . . . . . . . 29 5.1.2 NMR options . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.1.3 Running NMR jobs . . . . . . . . . . . . . . . . . . . . . 31 5.1.4 Calculating the induced current . . . . . . . . . . . . . . . 31

5.2 NMR shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.3 Band structure computation . . . . . . . . . . . . . . . . . . . . . 32 5.4 Density of states (DOS) . . . . . . . . . . . . . . . . . . . . . . . 33 5.5 Momentum density . . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.6 GW Interface and Usage . . . . . . . . . . . . . . . . . . . . . . . 34

5.6.1 GW Interface . . . . . . . . . . . . . . . . . . . . . . . . . 34 5.6.2 Example: Si . . . . . . . . . . . . . . . . . . . . . . . . . . 36

5.7 Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.7.1 Tensor Data . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.7.2 Line plots . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.7.3 Slice plots . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 5.7.4 Visualizing wave functions . . . . . . . . . . . . . . . . . . 44 5.7.5 Visualizing the charge density . . . . . . . . . . . . . . . . 44 5.7.6 Plotting the potential . . . . . . . . . . . . . . . . . . . . 45 5.7.7 Ball and stick model . . . . . . . . . . . . . . . . . . . . . 45

5.8 Spin polarized calculations and magnetic field . . . . . . . . . . . 45

5.8.1 Applying a magnetic field . . . . . . . . . . . . . . . . . . 45 5.9 Reading additional input . . . . . . . . . . . . . . . . . . . . . . . 46 5.10 Getting additional output . . . . . . . . . . . . . . . . . . . . . . 46

6 Funky extra flags 49 6.1 Optimization flags . . . . . . . . . . . . . . . . . . . . . . . . . . 49 7 Helpful tools and stuff 51

7.1 Paratec Perl Library . . . . . . . . . . . . . . . . . . . . . . . . . 51 7.2 tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

7.2.1 gaussdos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 7.2.2 bsfix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

8 Erratum etc. 53 8.1 Known difficulties . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6 CONTENTS

Chapter 1

How to get and run paratec

1.1 Where to find paratec

Paratec is located in

/usr/local/codes/paratec,

and subdirectories below. … This is no longer true. At the moment, there is no standard repository for paratec. ¨⌢

Paratec documentation can be found on the web at http://www.nersc.gov/projects/paratec http://www.nersc.gov/projects/paratec

1.2 How to run paratec

1.2.1 Input files

Required files:

paratec             the executable
input               the input file
Si_POT.DAT          the pseudopotential files
O_POT.DAT           ...

1.2.2 Output files

STRUC_CONVERGE    contains the progress of the structural relaxation
STRUC_LOG         a detailed log about the relaxation steps performed
ELASTO_PHONON     information from the H matrix: Bulk modulus, phonon
frequencies
PW_LOG            a log about all the PW calls done
OUT               output file of the current plane wave job
CHECKPOINT        contains all the information to continue a

                  structural relaxation. Also has the Hessian matrix
*.dx              files for visualization with the IBM Data Explorer
*.kh              files for visualization with Khoros

Note: As of paratec 4.5c, the output file OUT. 0 has been changed to OUT.

1.2.3 How to compile and run on different machines

Having obtained a gzipped tar file of the source e.g. paratec.xxx.tar.gz, unpack the source by typing

gzip -d -c paratec.xxx.tar.gz | tar -xf -

Then enter the paratec directory and run the configuration script by typing

cd paratec
./configure

For a summary of the options accepted by the global Makefile, see Table 1.1 or just type

make

Option to make Function
paratec  Build Paratec
tools  Build related tools
html  Build HTML documentation
microbench  Build microbenchmark for FFTs and MMMs
test  Build debugging and performance tests
clean  Remove object and test files
dist  Change version number and pack into a tar file
dist-pub  above plus removes NMR and TDDFT code

Table 1.1: Options for global Paratec Makefile

The machines currently supported are listed in Table 1.2 for reference. For each machine, a working directory in the source tree is created e.g. paratec/src/para/machine where the value of machine is given in the second column of Table 1.2. The binary files are placed in paratec/bin/machine. This facility allows the same source tree to be used by more than one machine simultaneously without any conflicts arising between object files.

1.2.4 Some tips on porting paratec

Due to heavy use of optimized library routines, paratec is not easy to port. Here are some steps on the way to port to a new architecture.

1. go into the config directory, and copy an established sysvars.machine and blurb.machine. Edit to meet your needs. The blurb file is displayed during configuration, and can be empty.

Machine  machine Description
Sun  sun  Sun
IBM RS6000  rs6k  Single-node IBM RS6000
IBM SP2  sp2  IBM SP2 at Cornell
IBM SP3 (NERSC)  sp3  Glen Seaborg at NERSC
IBM SP3 (NPACI)  sp3  Blue Horizon in San Diego
Cray T3E  t3e  Cray T3E at NERSC
SGI o200  sgio200 SGI Origin 200
SGI o2000  sgio2000 SGI Origin 2000 at NCSA
SGI PowerChallenge sgipc  SGI PowerChallenge at NCSA
HP/Convex Exemplar cvx
Linux x86  I386  Linux on x86
DEC Alpha  alpha  DEC Alpha
Hitachi SR2201  sr2201  Hitachi SR2201 at Cambridge
SGI o2000  hodgkin SGI Origin 2000 using new libraries

Table 1.2: Machines currently supported

2. Edit the ”configure” script to include the new machine.

3. Have a look at the files below...

Files in src/para

File  advice
flibcalls.m4h  Has some system dependencies on the use

of blas routines, i.e Cray versus the rest of the world. Will most likely need no changes.

fortranize.h  Depending on the compiler, c routines must

have a trailing underscore, or be in capitals in order for the linker to find them. If the linker doesn’t find the c subroutines, have a closer look here.

all to all.fp  May need some tweeking if there are hickups
data gather.fp  in the MPI implementation. Look for the
gspace gather.fp  machine-specific ifdefs there. Most likely will
isort gather.fp  need no changes.

setup packinfo.fp symmetrize tensor global.fp

zdiag scalapack.fp  Scalapack on the Cray uses Shmem, and that
rdiag scalapack.fp  requires some ugly workarounds.

file in src/macros/fft

fft macros.m4h

This one is quite a pain, especially since it uses m4. The standard quotation characters ‘’ have been replaced by []. The file fourier transform.fp gets preprocessed into fourier transform.p.f. To check correctness of the macros, look at fourier transform.p.f Look up the man on the fourier transforms (e.g man fft), add new macros accordingly to fft macros.m4h.

Don’t forget to provide for workspace in fft workspace.fp, and for a list of acceptable FFT lengths in adjust fft.fp. Once you have that done, you are through the worst. Have a look at fourier transform.fp. If the fourier transform interfaces are awkward, you might have to make some changes there as well. Make use of optimized routines to perform the convolution, rather than using the vanilla loop. Try if the preprocessor option FAST CONVOLUTE results in any speedup - it usually does for machines where no multiple-FFT routines are available.

At some point, the code was running also on cray vector machines. This is no longer the case, although still all the code is there, and it should not be too much work to get it going again.

David’s Short List of Files you Have to Edit

• fft macros.m4h

• flibcalls.fp

• libcalls.h

• inputread.c

• fft workspaces.fp

Tips on porting FFTW

FFTW (the Fastest Fourier Transform in the West) is a nice highly optimized portable FFT library. Paratec can now be compiled to use FFTW on platforms where it is available with only minimal changes needed. This should make it easier to port paratec.

To make FFTW work on a new platform, make sure the sysvars file includes the path of the FFTW library (for example, see sysvars.sgio200). Secondly, you must make sure that if the platform also supports a native FFT package that M4FFTFLAG is defined (again, see sgio200), and that adjustfft.fp uses that value.

I think that is all that should be needed. When you run the configure script, of course, you must say yes to fftw.

1.2.5 Some performance numbers on various machines

For performance numbers, please refer to Bernd Pfrommer’s Master’s thesishttp://civet.berkeley.edu/pfrommer/csthesis/csthesis.html.

1.2.6 How to change paratec

If you make a change to paratec, please make note of the change you make in the CHANGES file. This way you (and others) can tell what has changed in each version.

When you want to archive a version of the code, or give it to someone else, use the make dist target of the makefile. This saves only the necesary files, and not all the object files or binaries you may have generated. It then tars up the results and gzips them for your convenience.

make dist also prompts you for a new version number. It modifies the current version number (which is printed out in the OUT and PW LOG files, and shows up at the top of this documentation. It also updates the time stamp, which shows up in the OUT and PW LOG files. Then it names the tarball after the new version number. You can, of course, give it the same name as the old version, but if you change the number you will be able to keep track of which binary originated from which version of the code.

WARNING! make dist assumes your source code directory is called “paratec”. If you keep it in a directory with any other name, these may end up saving the wrong source code, or doing something else nasty.

If you add a new file to paratec, make sure to add it to the lists in the Makefile for the dist target.

1.3 Changes in Paratec 5.0

On the whole, paratec 5.0 uses the same input files as earlier versions of paratec. However, because the input code was completely rewritten in FORTRAN 90, several of the input structures have changed. This is where you look to find out what changes you will have to make to your input files. A more comprehensive list can be found on the home pagehttp://www.nersc.gov/projects/paratec

• The file formerly named OUT. 0 is now simply OUT. This is only an issue if you have scripts that assume it is called OUT. 0.

• All block inputs (i.e. inputs that start with begin blah now must close with end blah. Formerly some of these structures were closed with simply

     end.

• The statement coordinates absolute has acquired an underscore, so it

     is now coordinates absolute.

12 CHAPTER 1. HOW TO GET AND RUN PARATEC

Basic paratec usage

2.1 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 (ρ_{up} − ρ_{down})/(ρ_{up} + ρ_{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

2.2 Relaxing the crystal

Every paratec job involves relaxing the crystal. If you don’t want to relax, you simply set relax max iter to be 1.

2.2.1 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.
   Default: job relax

NOTE: If you specify job relax from checkpoint and there is no CHECK-POINT available, paratec will abort without performing any calculation.

2.2.2 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

2.3 Plane wave code parameters

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

2.3.1 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 (selfconsistent 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 4.0.5.

2.3.2 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 never larger than about 10 Rydbergs. The generation of the starting guess is not optimized, and 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.

2.3.3 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 Γ by (\frac{1}{8}, \frac{1}{8}, \frac{1}{8}) 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

2.3.4 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

2.3.5 Pseudopotential format

Paratec 5.1 has three different formats for pseudopotential 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

Formats 2 and 3 are in real space and include the atomic pseudowavefunctions. Format 1 is incompatible with the pulay tf mixing option. 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 http://www.fhi-berlin.mpg.de/th/fhi98md/fhi98PP.

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

begin pseudopotential
pp_data   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),

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

2.4 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 4.0.7 and mixing methods 4.0.12, to learn how you may be able to make your calculations run faster.

18 CHAPTER 2. BASIC PARATEC USAGE

Advanced relaxation parameters

3.1 Specifying the crystal structure

3.1.1 Inversion and symmetry

If there is a center of inversion, the code automatically tries to shift the atomic coordinates such that the crystal has inversion symmetry. If you want to inhibit that, put a line saying:

no restore_inversion

Finally, you can specify whether you want the symmetry-operations to be generated by the code, or if they should be read from the file SYMMETRYOPS. If number symmetry ops is followed by -1, then the symmetry operations are determined by the code itself. If number symmetry ops is > 0, the matrices for the symmetry operation are read in from the file. The file format is the same as the printout in e.g. OUT. For an example file, see SYMMETRYOP-Shttp://civet.berkeley.edu/pfrommer/SYMMETRYOPS

   Default: number symmetry ops -1

If number symmetry ops is followed by -2, then the symmetry operations are determined by the code itself, but nonzero fractional translations are ignored. This could be useful for supercells.

3.1.2 Specifying an initial strain

You may optionally specify an initial deformation tensor using the format…

begin deformation
0.00 0.00 0.02
0.00 0.00 0.00
0.00 0.00 0.00

end deformation

This is the deformation tensor applied to your structure. It is defined by the following equation:

a = (1 + D) · a_{o}

where a are the deformed lattice vectors, and a_{o} are the lattice vectors specified in your input file. This can be an easy way to apply a shear strain to a cubic crystal, for instance.

There are times when it is useful to ensure that a given direction is held fixed during the relaxation process. This does not involve a constraint on the relaxation, as there are three rotational degrees of freedom. One can also arrange to maintain fixed a plane which contains the fixed direction.

To keep the z fixed, and also keep the zx plane from rotating, you could specify…

mugwump_vector     0 0 1
desideratum_vector 1 0 0

3.2 Relaxing a structure

For the impact of the job variable on the relaxation process, see section 2.2.1.

3.2.1 Controlling the relaxation process

If job relax or job relax from checkpoint or job relax recycle H is cho-

sen, the variable relax what allows different ways to relax a structure:

   relax what force

or

   relax what stress

or

   relax what force and stress
   Default: relax what force and stress

Another variable allows to manipulate the way the relaxation is done:

Set

   relax how gradient

to simply follow along the forces and stress, without updating the inverse of the Hessian matrix.

Default: update the Hessian

With the variable relax how the line-minimization can be altered also:

relax_how slow        # do 3 steps per line minimization always
relax_how normal      # do up to 3 steps, but only if necessary
relax_how fast        # dont do line minimization at all

   Default: relax how normal

Note: This can be combined with the gradient option with a statement like

relax how normal gradient

The parameter lambda limit allows you to limit the size of relaxation steps to a certain number. Normally, lambda should be on the order of 1. Sometimes, especially in the beginning of the iteration, and if the H matrix is initialized to small, lambda can become to big, and the algorithm bombs out. It is then advisable to limit lambda with a line like that:

   lambda limit 4.0

Default: lambda limit 1e6 (no limit)

The statement starting distortion followed by nine floating point numbers allows to set an initial distortion of the lattice. This feature is only for debugging, and not recommended for use.

3.2.2 Pressure

An external pressure (units are GPa) can be applied with:

   relax pressure 10.0

If an entry like

   relax pressure adjust

is found, the pressure is chosen from the first iteration, taking it to be 1/3 of the trace of the stress.

   Default: relax pressure 0.0

3.2.3 Constraining the relaxation

Paratec offers the possibility to fix atoms to their input coordinates, meaning they will not move during the relaxation process. This is simply done by zeroing out the Hellman-Feynman forces. To fix atoms to their original positions, use the stay put structure. The following example fixes the positions of atoms number 1, 3, and 10.

begin stay_put
stay_put 1
stay_put 3
stay_put 10
end stay_put

Default: all atoms are relaxed.

As of paratec 4.5, paratec also offers the option of fixing certain strain coordinates while allowing the others to change during the relaxation. This is done by using the keyword stay put lattice within the stay put structure. The following example fixes the ϵ_{xx} and ϵ_{yy} strains, while allowing the z axis to expand or contract.

begin stay_put
stay_put_lattice xx

stay_put_lattice yy
end stay_put

Default: all strain coordinates are relaxed.

NOTE: The syntax of the stay put structure has changed in version 5.0. The end at the end is now an end stay put.

3.2.4 Initialization of the inverse of the Hessian matrix

The Broyden scheme used for the relaxation of the forces and the stress updates an inverse of the Hessian matrix, which essentially describes the enthalpy parabola around the minimum. There are some open parameters which determine the values to which the matrix is initialized. They basically determine the length of the trial steps during the first few iterations:

relax_eps_factor   0.1     initialization perfactor for the stress part
of H
relax_coord_factor 0.1     initialization perfactor for the force part
of H

Try with the values above for the first run. Adjust later such that the first few steps give a lambda of about 0.5 to 1.

Default:

relax_coord_factor 1.0
relax_eps_factor   assuming the bulk modulus is about 1 Mbar

Alternatively, you can specify the physical quantities to which the above two variables correspond. This is the recommended way:

estimated_bulk_modulus  1.0                 Bulk modulus in MBar
estimated_optical_phonon_frequency 12.0     TO phonon frequency f at
                                            Gamma in THz)

Diamond would have a bulk modulus of about 4.4 MBar and around 35 THz for the optical phonon frequency.

For molecular solids, the option

decouple_coordinates yes

causes a different initialization of the Hessian matrix, which decouples the coordinates of the molecule from the lattice vectors to first order. This means that if the lattice vectors change 10 percent, the cartesian coordinates will change by about 1 percent. IF THIS OPTION IS USED, THE SYMMETRY OF THE CRYSTAL IS NOT PRESERVED!.

Default:

decouple_coordinates no

Advanced plane wave code parameters

4.0.5 Plane wave code jobs

Whenever the plane wave code is called for a given structure, it can perform several jobs, depending on the job list:

begin pw_jobs
pw_job   first_job
pw_job   second_job
...
pw_job   last_job
end pw_jobs

If you do a structural relaxation, the list of jobs will be executed for every configuration of the relaxation, not just for the final one. Allowed jobs are:

• scf perform a self consistent field calculation, and compute force/stress.

• nmr shift compute the diamagnetic susceptibility of insulators for the nmr shift, see section 5.2.

• nmr plot plot the diamagnetic susceptibility computed in a previous calculation, and read from the file CHI. It will also plot the induced current, with the B-field aligned along the “ket” direction.

• band structure compute and plot the band structure.

• pot plot If also the output flag potplot, is set, this job plots the local part of the effective KS potential as a line plot, and produces a file for the IBM data explorer also. The G=0 component of the potential is set to 0. The printed out eigenvalues and the fermi level include VXC(G=0), so you must subtract this offset to find the alignment with respect to the potential. pot plot can also be used to plot the charge density, but not in conjunction with the energy window features. To plot the total charge, set the output flag cdplot. Also, the electric field along a line and at the atomic positions is plotted if the output flag efield is set.

4.0.6 Exchange correlation functional

Currently three different exchange-correlation functionals are implemented: LDA (Ceperley-Alder) and GGA (Perdew-Wang 91, Perdew Burke Ernzerhof). Spin polarized calculation is supported in all cases. You can set the exchangecorrelation functional with:

   exchange correlation name of functional.

By default:

   exchange correlation ceperley alder

The GGA can be switched on with:

   exchange correlation perdew wang 91 or
   exchange correlation perdew burke ernzerhof

The shortened forms pw91 and pbe are also recognised options to exchange correlation. Note that for spin-polarized calculation you should put the following line in

”input” file

   number of spins 2

For non-spin-polarized calculation you should put

number of spins 1 (This is default)

Note that direct minimization (optimize insulator option) is not supported for spin-polarized calculation.

4.0.7 Diagonalization Methods

The iterative method for the diagonalization to obtain the electronic wavefunctions is chosen by

diagonalization_method    Mauri, Grassmann, or Grassmann_metal

The default is Mauri. It is good for insulators and semiconductors at low temperature. This is the original algorithm in previous versions of Paratec. Grassmann gives similar results with little advantage to either method, comaprisons can be found at http://www.nersc.gov/projects/paratec/METHODS/comp TE min.htm http://www.nersc.gov/projects/paratec/METHODS/comp TE min.htm. For metallic surfaces and long skinny cells, both of these algorithms suffer problems. The Grassmann metal algorithm was developed to handle these systems. It uses the occupations of the electronic wavefunctions to facilitate the diagonalization. A write-up on this method can be found at http://www.nersc.gov/projects/paratec/METHODS/scf metal min.htm http://www.nersc.gov/projects/paratec/METHODS/scf metal min.htm.

4.0.8 Diagonalization

The accuracy of the diagonalization scheme is determined by the keyword accuracy diag. The diagonalization stops when the gradient (i.e. the residual) is below accuracy diag. The accuracy of the eigenvalues is normally at least as good, often an order of magnitude better.

Example:

   accuracy diag 1e-10

(requires high accuracy)

   Default: accuracy diag 1e-12

Advice: The energy is normally converged to 1e-6 only, so an accuracy of 1e-7 is often sufficient. However, to get the forces accurate to n digits, one should require accuracy diag to be less than 1e-(2n). For example if you set relax accuracy to 1e-4, you must set accuracy diag to 1e-9 or smaller.

The variable max iter diag determines when to stop iterating on the trial eigenvectors, if no convergence is achieved. This variable should be set such that during the first two or three self-consistent iterations, the diagonalization is stopped by max iter diag. There is no point in converging the eigenvectors for a Hamiltonian which is not nearly close to self-consistency.

Example:

   max iter diag 30
   Default: max iter diag 5

With the variable bandgap, you can set the expected bandgap of an insulator in eV. In this case, only the minimum number of bands are computed, and the diagonalization scheme uses this bandgap for performance optimization. If you screw up bandgap, the performance can get terrible, and the diagonalization routine can bomb out, so be careful.

In case of matrix diagonalization bombouts, try to increase: diagsafety

2.0
   Default: diagsafety 1.0

You can also try randomizing the starting guess by a line like:

   randomize diag start guess 1.0
   Default: randomize diag start guess 0

For the direct energy minimization, it is often useful to checkpoint the wave functions and charge density during the iteration. This can be done with

   checkpoint emin frequency

Example:

checkpoint emin 50 writes the wave functions to disk every 50 iterations. By default, no checkpointing is done.

4.0.9 Kpoint grid

Specify

   number kpoints 0

to have the kpoints generated, or use a positive integer number, and the corresponding number of kpoints will be read from the file KPOINTS, which must be in the current working directory.

   Default: number kpoints 0

As a special hack, you can specify the number of kpoints to be negative. This will also generate kpoints, but will NOT REDUCE THEM:

   Default: number kpoints -1

Implemented is a k-point generation scheme a la Monkhorst-Pack. The parameters are k grid, followed by 3 integers, which give the number of grid points along the 3 reciprocal lattice vectors for the unreduced k-points. An entry called k grid shift followed by three real numbers allows to shift the k-grid.

Example:

   k grid 4 4 6
k grid shift 0.5 0.5 0.5

Shifts the 4x4x6 grid by 1/2,1/2,1/2 (of one grid cell)

   Default: k grid 1 1 1
k grid shift 0 0 0

4.0.10 Number of bands

The entry number bands gives the number of bands that are computed. It should always be slightly more than the minimum number of bands required, unless a wide-gap insulator is studied. If the number is more than the minimum, the Grassmann metal diagonalization is suggested. For systems with a large number of states at the Fermi surface (such as metallic surfaces), as many as twice the minimum number of bands is optimum. The computational effort grows linearly with the number of bands, so do not pick the number unreasonably large either.

Example:

   number bands 12

Default: There is no default

4.0.11 Occupation numbers and gaussian smearing

After the energy levels at the various kpoints have been computed, they are occupied with electrons according to their energy.

With the variable occupy levels you can control the way the occupation numbers are computed.

occupy_levels   normal      occupy up to fermi level, which is common to
                            all kpoints or
occupy_levels   fixed       occupy at each kpoint separately, don't
                            level out
occupy_levels   antiferro   I (David Roundy) don't know what this does,
                            but it is in the code, and apparently does
                            antiferromagntetic filling
   Default: occupy levels normal

You must specify the width of the gaussian smearing IN ELECTRON VOLTS that is done to the energy levels when they are occupied in the flevel subroutine. For insulators, choose a very small number, say 0.001 eV.

Example:

gaussian_smearing 0.1        sets the width of the gaussian to 0.1 eV
   Default: gaussian smearing 0.05

4.0.12 Self-consistent mixing methods

The potential mixing method that obtains a new input potential for the next self-consistent step is chosen by

mix_method             broyden, pulay_kerker, or pulay_tf
Default = pulay_kerker

If broyden, the following additional parameter is used.

mixing_energy_cutoff      Sets the broyden mixing energy cutoff.
Default = 5 Ryd

If pulay tf, the following additional parameters are used

PTF_num_init_PK      gives the number of pulay_kerker mixing
                     steps before the pulay_tf mixing is used.
Default = 2          best performance for Al 20 layer surface
PTF_max_iter_cg      gives the number of conjugate gradient
                     iterations for solving the TFW equations
default = 40

The Broyden method should probably never be used unless someone has a sentimental attachment. For any size relatively homogeneous system, pulay kerker should be used. For large inhomogeneous systems, pulay tf should be used; especially for systems with a regions of vacuum such as sufaces, nanotubes, and molecules. Additional information can be obtained at http://www.nersc.gov/projects/paratec/analysis.htm http://www.nersc.gov/projects/paratec/analysis.htm.

4.0.13 Self-consistent iteration

The accuracy with which the self-consistent potential is computed can be set

with the variable potential convergence criterion.

Example:

   potential convergence criterion 1e-4

Requires that the largest difference between the current and the previous potential (in G-space) be 1e-4.

   Default: potential convergence criterion 1e-6

(this default is good enough for an energy accuracy of at least 1e-6, since the error enters the potential linearly, but the energy quadratically, so the energy converges faster than the potential) Nothing smaller than 1e-6 is suggested.

If no convergence is achieved within max iter scfloop, the self-consistent iteration terminates, and goes on to compute the stress and forces.

Example:

   max iter scfloop 50
   Default: max iter scfloop 40

Screening and mixing parameters

For the first iteration of each plane wave programm call, the ionic pseudo cores have to be screened with an initial guess for the charge density. With

   screening type atomic

the initial screening charge density is calculated from the valence charge density of the atoms.

Setting screening type previous screens with the charge density from the file CD if that is available, or with the atomic charge otherwise. This means that the charge density from the previous call to the plane wave code is reused. The mixed potential itself is stored in the file VMIX and if available, may be re-used

by setting screening type vmix.
   Default: screening type previous

Broyden mixing

The mixing parameter variable determines the initialization of the inverse of the Hessian matrix. A small mixing parameter will lead to smaller movements of the charge density during the first few iterationss. Later, when the inverse of the Hessian matrix has been constructed, the mixing parameter does not change the algorithm much.

mixing_parameter 0.33          initialize H to 0.33 of standard size
   Default: mixing parameter 0.77

For smaller wave vectors, a Broyden type mixing is used, for larger wave numbers a linear mixing is employed. The energy cutoff for this Broyden mixing may be set with the parameter:

mixing_energy_cutoff   5.0   Sets Broyden mixing cutoff to 5.0Ry

Often a small fraction of the total cutoff is sufficient, say 5 Ry. Default: 5 Ry

Linear mixing

In some cases, Broyden mixing is too aggressive. In this case, you can switch to linear mixing with a line

linear mixing a1 a2 a3

which mixes the potentials according to the formula:

a_{mix} = a_{1} + a_{2}e^{−a_{3}|G⃗|^{2}}  (4.1)
v_{new}(G⃗) = a_{mix}v_{in}(G⃗) + (1 − a_{mix})v_{out}(G⃗)  (4.2)

That way, the higher wavenumbers are mixed stronger, while the tricky small G components are mixed less.

Defaults: linear mixing 0 0.6 0.05

To switch linear mixing on, a1 has to be > 0. This switches off Broyden mixing. a1=0.2 seems to be a good value to start with.

4.0.14 Options for improved stress

The stress generally converges slower than the total energy, because the stress is normally evaluated assuming constant number of plane waves, which is a bad approximation.

The code has two options to improve that. The first trick is due to Michel Cote and goes as follows. After the selfconsistent calculation has been performed at normal cutoff, a second calculation with higher cutoff is done, but WITHOUT allowing the charge density to relax. This turns out to be still a very good approximation, but the second calculation is somewhat expensive. The cutoff energy (in Rydberg) for the second run is set with the keyword

polished energy cutoff.

Example:

polished energy cutoff 50

Default: switched off The second trick is to modify the kinetic energy according to Bernasconi et

al. You can set the height A, position E_{0}, and width σ (all in Rydberg) of the step function that gets added:

|k⃗ + G⃗|^{2} → |k⃗ + G⃗|^{2} + A(1 + erf((|k⃗ + G⃗|^{2} − E_{0})/σ))

with the keyword modify kinetic energy A E_{0} σ.

Example:

modify kinetic energy 20.0 10.0 1.0.

If you just do

modify kinetic energy on

a reasonable setting is assumed: A = 2E_{cut}, E_{0} = E_{cut}, σ = E_{cut}/10.

Default: modify kinetic energy off

30 CHAPTER 4. ADVANCED PLANE WAVE CODE PARAMETERS

Getting extra outputs and other interesting stuff

5.0.15 External potential

Paratec has an option to apply an external (electric) potential of the form:

V = V_{0} sin(G⃗ · r⃗).  (5.1)

This is done with an entry like:

   electric potential kvec 1 0 0 1

The last number switches the field on (=1) or off (=0), the first 3 INTEGER numbers specify the wave number in reciprocal lattice coordinates.

The amplitude (in Rydbergs) is entered as:

   electric potential strength 0.3

5.1 NMR shift

The theory behind paratec’s calculation of NMR chemical shifts can be found in the paper by: Chris Pickard and Francesco Mauri Phys. Rev. B 63, 245101 (2001 ) The following two papers provide further information Gregor, Mauri, Car J. Chem Phys (1999) (GMC) Mauri, Pfrommer, Louie PRL (1996) (MPL) Valid for insulators only, spin polarized and non-spin polarized

5.1.1 Additional Files

Required Input files:

Si_AEPS.DAT          ascii files containing data for the reconstruction
O_AEPS.DAT           ...

32CHAPTER 5. GETTING EXTRA OUTPUTS AND OTHER INTERESTING STUFF

Output files:

BLOP.DAT             contains the Blochl operators
CHI                  contains the magnetic susceptibility in g space
CHKPTCHI.00x         checkpoint data for the run. Can be written after
                     each kpoint and used for a re-start. One file for
                     each processor
JGS                  contains the directly calculated current in g space

5.1.2 NMR options

The following is a brief list of all the paratec options directly related to NMR calculations. Further details are provided below.

• pw jobs nmr shift Computes current due to a perturbing magnetic field, the local diamagnetic susceptibility and the chemical shieldings at each nuclear position

• pw jobs nmr plot Plots the induced current

• nmr q 0.001 sets the q (in cartesian coordinates), which determines the periodicity of the applied external field. Default: nmr q 0.01

• checkpoint nmr The code will checkpoint the NMR shift calculation after each kpoint. You can recover from a crash by giving an input flags of

     chkptchi.
   • output flags eg output flags crystal current

– original Use the original MPL approach for calculating shielding. Do not use this option, it is kept for historical reasons only!

– crystal Use the approach of PM for calculating the shielding. This is the default option.

– eqn3 For molecules only. Use the Molecular approach of PM. Refers to equation 3 of GMC.

– eqn8 For molecules only. Use the Molecular Sum Rule approach of PM. Refers to equation 8 of GMC.

– current For pw jobs nmr shift this writes a file JGS containing the current in g space. For pw jobs nmr plot this reads from JGS and does line plots or slice plots as described in section 5.10. The magnetic field is in the direction of ket

– nmrshift To get a plot of the susceptibility (section 5.10).

   • input flags eg input flags wavefn chkptchi

– chkptchi This causes a restart from a CHKPTCHI file – wavefn Use the wavefunctions from a previous scf run as a starting guess for the first diagonalization. Recommended. Essential for a spin polarized calculation.

   • optimize flags eg optimize flags memory

– memory optimizes for memory in the calculation of NMR shifts, at a small expense in additional CPU time.

– nocgsolveguess switches OFF the better starting guesses for the conjugate gradient solver in the NMR subroutine. Using this option will generally result in SUBSTANTIALLY SLOWER execution, but will reduce memory cost.

5.1.3 Running NMR jobs

Generating the input files.

A version of the atm pseudopotential code and a script genpot are available to generate the AEPS.DAT files. The atm code can also be used to find the diamagnetic shielding contribution of the core states.

Before an nmr run an scf run must be done to obtain a good charge density. A accuracy diag 1e-12 is probably sufficient. As the nmr routine only works for insulators at the moment use optimze flags insulator.

Crystal Method This calculates the shielding tensors for both crystals and molecules (in the supercell approximation) using the method of PM. The method uses the wavefunctions at points separated by q to obtain a numerical derivative wrt q. Consequently a high degree of accuracy is needed in the wavefunctions and you should set accuracy diag 1e-14. Results should be stable wrt changing q. nmr q 0.001 is the default and appears to give stable results.

Eqn3 and Eqn8 Methods These are formulated for isolated molecules in the unit cell and are several times faster than the crystal method. A accuracy diag 1e-12 should be ok. The molecule must be centred on the origin of the unit cell.

Comparisons of shieldings calculated with the three methods can be found in PM.

5.1.4 Calculating the induced current

Paratec computes the induced current and uses it to calculate the magnetic susceptibility CHI. It is possible to obtain the current directly or from CHI. For a complete basis set the two currents will be identical (see appendix of GMC for details) in practice the two differ very slightly (most likely more due to the n on-local pseudopotential than the basis set). For completness it is possible for paratec to plot out the current dervived from either method. I (jry) recommend obtaining the current from CHI, though you should probably check it agrees with the direct method. A the moment there is a problem with current calculated with the eqn3 and eqn8 methods. This will be fixed soon.

Obtaining the current from CHI During pw jobs nmr shift the magnetic susceptibility is written to a file CHI. During pw jobs nmr plot CHI is reread and the induced current calculated. For plotting perposes the magnetic field is 1Tesla in the direction of ket eg ket 0 0 1 . Plots along a line or plane can be defined see section 5.10. As the default the current will be plotted as a 3 component vector. To plot the projection of the current on the direction bra set output flags curproj and, for example

bra 1 0 0.

Obtaining the current directly By setting output flags current the file JGS is written during pw jobs nmr shift.

If output flags current is present during pw jobs nmr plot plots can be

made as in the previous section. (This is only implemented for the crystal method, only trivial modifications are needed to make it work in eqn3 and eqn8)

5.2 NMR shift

If a pw job nmr shift entry is encountered, paratec computes the local diamagnetic susceptibility. The variable

   nmr q 0.001

sets the q⃗ (in cartesian coordinates), which determines the periodicity of the applied external field.

   Default: nmr q 0.01

Set the output flag nmrshift to get a plot of the susceptibility (section 5.10).

If a line

   checkpoint nmr

is entered, the code will checkpoint the NMR shift calculation after each kpoint. You can recover from a crash by giving an input flags of chkptchi.

5.3 Band structure computation

An entry pw job band structure computes the band structure. You should set screening type to previous, and provide a converged CD file for the starting density. THERE WILL BE NO SELFCONSISTENT ITERATION, i.e. the bandstructure is plotted corresponding to the CD file, or to the atomic valence charge if no intact CD file is found. Furthermore, YOU HAVE TO USE THE FERMILEVEL OF A PREVIOUS SCF RUN. You do so by giving a line like:

   fermi level fermi energy in Ryd

   Default: fermi level 0

If you want to use a different number of bands from that used in an earlier SCF loop, you may specify this number using:

number_alt_bands <number of bands for the bandstructure>

You also have to determine the lines in the BZ along which you wish to plot the bandstructure. This is done by specifying a sequence of lines in the following way:

begin bandstructure
label left_middle_right
kpoint startx  starty  startz   endx  endy  endz   number_of_bins
label left_middle_right
kpoint startx  starty  startz   endx  endy  endz   number_of_bins
..
label left_middle_right
kpoint startx  starty  startz   endx  endy  endz   number_of_bins
end bandstructure

NOTE: The syntax of this structure has changed in version 5.0. Previously this ended with just end, and this has changed, so that now you need to end

with end bandstructure.

All coordinates are given with respect to the reciprocal lattice vectors. The labels are mandatory; they have to be triples, separated by an underscore. Notice that xmgr allows greek characters using escape sequences (see example below).

The kpoints will be generated into the file KPOINTS. The output of the calculation can be found in the file BANDSTRUC. Energies are in eV, and are plotted against an x-coordinate that corresponds to the physical distances in cartesian k-space, i.e. 1/a.u. The file can be displayed with the popular graphics program “xmgr”.

Example:

begin bandstructure
label \8G\4_\8\D\4_X
kpoint 0.0 0.0 0.0         0.0 0.0 0.5            10
label _V_W
kpoint 0.0 0.0 0.5         0.0 0.5 0.5            10
end bandstructure

5.4 Density of states (DOS)

To compute the density of states, a self-consistent calculation has to be done first to get the wave functions. To generate a DOS, add to the output flags either dos or angdos. This will produce a file “DOS” in a suitable format for xmgr (it is an ascii text file, and has comments in there). A tetrahedron interpolation scheme is used, where each of the Monkhorst-Pack grid cells is divided into six tetrahedra. If you compute the angular momentum resolved DOS, you must also specify the radius of the muffin-tin spheres (see section

The density of states is output in units of states/eV/unit cell. It does not count both spin states. 2.1).

5.5 Momentum density

To compute the momentum density, first compute the selfconsistent charge density, and then do a band structure calculation. There, you specify those lines in the BZ along which you want to know the momentum density. To get the momentum density, you must add the flag momdens to the list of output flags. Example:

begin pw_jobs
pw_job scf
pw_job band_structure
end pw_jobs
output_flags momdens

You also need to specify the Fermi level according to a previous scf run:

fermi level fermi energy in Ryd

If you don’t specify this correctly, you will get garbage!

5.6 GW Interface and Usage

GW-interfaced Version of Paratec with f-states was implemented by Eric Chang. Further development was done by Xavier Blase, Gian-Marco Riganese and others.

5.6.1 GW Interface

To calculate dielectric function, we need 2 wavefunctions from paratec. They are called CWFE and CWFEq. CWFE requires lots of conductance bands, but CWFEq only requires wavefunction of bands (valence bands + 2).

Procedure:

(1) Generate a converged charge density CD, save this CD to a file CD.conv. You must use the same charge density each time when you perform the following jobs.

(2) Compute CWFE. In the input file of paratec, use the options

number_bands 90 # this number means the number of conductance bands used in
outputflags gwr # the summation of dielectric function calculation

If the system does not have inversion symmetry, you need to run the complex version:

number_bands 90 # this is a large number compared to the number of valence bands
outputflags gwc

After run paratec, rename the file GWR (or GWC for complex version) to CWFE. You may need to use fewer or more bands in the option number bands.

(3) Compute CWFEq. In the input file of paratec, use the options

number_bands 5         # this number is (number of valence bands + 2)
outputflags gwr        # use gwc if you want to run complex version of the GW code
gwshift 0.005 0.005 0.005

You should use number bands slightly more than the number of valence bands in the system. gwshift gives a wavefuction at slightly shifted grid. Rename the file GWR to CWFEq after paratec computation.

GWC or GWR —¿ CWFE for the shifted grid (need many bands) GWC or GWR —¿ CWFEq for slightly shifted+q grid (need only valence bands,need to read k-points from a file)

(4) Compute VXC and CD95. VXC is the information of exchange-correlation functional; CD95 is the charge density. In the input file of paratec, use the option

outputflags gwscreening

(5) Compute wfn0. In the input file of paratec, use the option

k_shift 0 0 0
number_bands 90        # test the convergence of this number
output_flags gwr       # use gwc if you want to run complex version of the GW code

Rename GWR or GWC to wfn0.

(6) Now you are ready to use CWFE and CWFEq to compute eps0mat and epsmat. eps0mat, epsmat, wfn0, VXC, CD95 are required to do GW calculation.

5.6.2 Example: Si

If you want to run the code for Si, do the following:
(1) mkdir Si
(2) cd Si
(3) mkdir paratec; mkdir xi; mkdir sig
(4) if you type pwd you should have the subdirectories:
 paratec  sig      xi
paratec -> has the LDA/pseudopotential/planewave code
sig -> has self-energy program
xi -> has screening program
 (ignore subdirectory epsilon)
(2) in your paratec directory you should have:
the following input files:
input.cd
input.cwfe
input.cwfeq
input.wfn0
the pseudopotential:
Si_POT.DAT
and the executable:
paratec.mpi
--- to make the executable, paratec.mpi, do the following:
   cd ~
   mkdir source
   cp paratec.tar source/
   cd source
   tar -xvf paratec.tar (untar the tar file)
   cd paratec
   make paratec
   cd ~/Si/paratec
   ln -s ~/source/paratec/bin/paratec.mpi . ! paratec.mpi is the
                                              executable

----To convert the ascii pseudopotential to the binary form in which
it is read, do the following:
   cd ~/source/paratec ! same directory as above
   make tools
   cd ~/Si/paratec
   ln -s ~/source/paratec/bin/kbascbin .
   ./kbascbin Si_POT.ASC.DAT Si_POT.DAT

Now we are ready to run paratec for four times, one time for each
of the four input files:
input.cd ----> compute the ground state density, and screening.
               use output_flags gwscreening, and save CD.
input.cwfe -----> compute the wavefunctions on the
                  4x4x4 grid, shifted 0.5 0.5 0.5. In this
                  case we compute 5-10 times the number
                  of valence bands. Notice the option in the input:
                  output_flags gwr. This means, produce real
                  wavefunctions to be read by the program in directory
                  Xi. Also, make sure that you converge the wavefunctions
                  to a high accuracy.
                  By the way, for a system with inversion symmetry,
                 you need to say, output_flags gwc.
input.cwfeq -----> compute the wavefunction on the
                   4x4x4 grid, shifted 0.5 0.5 0.5 plus some small
                   q-vector, specified by the option in the input,
                   e.g. gw_shift 0.00 0.00 0.001. Compute only one more
                  than the number of valence bands in this case. Use
                  option output_flags gwr.
input.wfn0  ---> computes wavefunctions on unshifted grid, 5-10 number
                 of valence wavefunctions, to be read by

In summary, we have the following new features in paratec
which enable the interface between paratec and GW:
gw_shift 0.00 0.00 0.001
output_flags  gwr,gwc, or gwscreening
gwr -> produces real wavefunctions
gwc -> produces complex wavefunctions
gwscreening -> produces CD95 and VXC (for the sigma code)

So do the following:

cd ~/Si/paratec  ! go to the right directory
cp input.cd input  ! compute ground state charge density
mpprun -n 8 ./paratec.mpi ! run paratec
cp CD CD.save             ! save ground state charge density
mv CD95 ../sig  ! CD95 is the charge density read by self-energy program
mv VXC ../sig   ! VXC is the exchange corr
e. read by self-energy program
cp input.cwfe input ! computes wavefunctions on 4x4x4 0.5 0.5 0.5 k grid
                    ! computes about 5-10 times number of valence bands
cp CD.save CD
mpprun -n 8 ./paratec.mpi ! run paratec
mv GWR ../xi/CWFE         ! move wavefunctions (GWR or GWC) to xi
cp input.cwfeq input ! computes wavefunctions on 4x4x4 0.5 0.5 0.5 +
                                          small q  k grid for xi
                    ! only need one more than number of valence bands
cp CD.save CD
mpprun -n 8 ./paratec.mpi ! run paratec
mv GWR ../xi/CWFEq
cp input.wfn0 input   ! computes wavefunctions for self-energy
                     ! 4x4x4 unshifted
                     ! 5-10 times number of valence bands
cp CD.save CD
mpprun -n 8 ./paratec.mpi
mv GWR ../sig/wfn0

(2) calculate static screening matrix:
  q is an element of 4x4x4 unshifted grid (8 q-points)
  k,k' are elements of 4x4x4 shifted grid (10 k-points)
  q is {k-k'| k,k' shifted}
  RPA screening (xi polarization)
   xi(G,G';q) =  Sum(cvk) [ M(cvkqG) M*(cvkqG') /(E_ck - E_vk+q) ]

     M(cvkqG) = <ck|exp(i(G+q).r)|vk+q>
  eps(G,G';q) = delta_GG' - [4pi/(q+G)^2 ]xi(G,G';q)
  if G=0  and q=0
  then we let
   eps(0,G';q) = delta_0G' - lim(q->0) [4pi/(q)^2] xi(0,G';q)

   program takes inverse of eps(G'G;q) for every q:
      -> epsmat (for q not equal to gamma)
      -> eps0mat (for q equal to gamma)

cd xi  ! you should have the following files:
           CWFE
           CWFEq
           xi0.inp input file
           xi0 (or cxi0 for complex version) executable
xi0.inp looks like
epsilon_cutoff          6.25  ! compute eps matrix up to G^2,G'^2 < 6.25 Ry)
number_bands          44  ! number of bands in summation (cond.+valence)
wavefunction_cutoff         12.0 ! same as in paratec
number_qpoints  8   ! number of q-points to compute eps(GG';q)
band_occupation 4*1 40*0 ! band occupation
begin qpoints
  0.00      0.00      0.001       1.0 1   ! q-point,divisor, 1 or 0
  0.00      0.00      0.25       1.0 0    !  copy from paratec
  0.00      0.00      0.50       1.0 0
  0.00      0.25      0.25       1.0 0
  0.00      0.25      0.50       1.0 0
  0.00      0.25      0.75       1.0 0
  0.00      0.50      0.50       1.0 0
  0.25      0.50      0.75       1.0 0
end
#Values for the 1st step
evs     0.0000
evdel   0.0000
ev0     0.0000
ecs     0.0000

ecdel   0.0000
ec0     0.0000
mpprun -n 8 ./xi0    ! run xi
mv eps0mat ../sig
mv epsmat ../sig
always check Xi(0,0;q->0) and eps^(-1)(0,0;q->0)
in xi0.log
1/eps^(-1)(0,0;q->) should be the macroscopic epsilon,
                    local field effects
1/(1-4*pi*Xi(0,0;q->0) should also be (roughly) the macroscopic epsilon,
                     non local field effects
(3)  Look at Hybertsen/Louie, PRB 34,5390.
formulas:
  (34a) we calculate Sigma_(SEX)
          has two terms:  1 + Omega/(DeltaE^2+wp^2)
              first term -> Gv (small v is bare coulomb interaction)
              second term -> is part of G(W-v)
  (34b) we calculate Sigma_(COH)
               other part of G(W-v)
We compute

     <nk|Sigma(E)|nk>, Sigma = i G_0 W
cd sig
in your directory you should have:
   eps0mat
   epsmat
   CD95
   VXC
   wfn0
   sig.inp
   sigma (or csigma)

looking at sig.inp:
number_kpoints       2  # of k-points in <nk|Sigma(E)|nk>
gmaxs     2.50000   # G-vector sum for screened interaction, G<gmaxs
gmaxb     5.00000   # G-vector sum for bare interaction G<gmaxb
ncband  44         # number of bands in summation (n1 in notation of paper)
nbnmin    1     # min band to calculate Sigma for (n in notation of paper)
nbnmax    8
dw        1.0   # don't change.
epscut    3.0   # don't change
ieps      0     # don't change
ecut     12.00  # wavefunction cutoff
number_offdiag  3 # number of offidagonal elements
begin kpoints
       0.0       0.0       0.0       2.0
       0.0       0.0       0.5       1.0
end
begin offdiag
    1  3               <1k|Sigma(E)|3k>
    1  2
    2  1
end
evs     0.0000               ! scissor shift
evdel   0.0000
ev0     0.0000
ecs     0.0000
ecdel   0.0000
ec0     0.0000
mpprun -n 8 ./sigma    ! run Sigma

output looks like this (sig.log) :
           k=   0.000   0.000   0.000
  i    elda      e0       x      sx      ch     sig      v0    efso    enew (ev) spin=1
  1   3.884   3.884 -17.107  13.287  -7.228 -11.048 -10.429   3.265   3.467
  2  15.775  15.775 -12.700   9.382  -8.410 -11.727 -11.229  15.277  15.383
  3  15.775  15.775 -12.700   9.382  -8.410 -11.727 -11.229  15.277  15.383
  4  15.775  15.775 -12.700   9.382  -8.410 -11.727 -11.229  15.277  15.383
  5  18.356  18.356  -5.672   3.732  -7.827  -9.767 -10.051  18.640  18.579
  6  18.356  18.356  -5.672   3.732  -7.827  -9.767 -10.051  18.640  18.579

44CHAPTER 5. GETTING EXTRA OUTPUTS AND OTHER INTERESTING STUFF

  7  18.356  18.356  -5.672   3.732  -7.827  -9.767 -10.051  18.640  18.579
  8  19.096  19.096  -5.770   4.011  -8.825 -10.584 -10.804  19.316  19.268
 1  1 --real-       -17.107  13.287  -7.228 -11.048 -10.429 del=  -0.620
 1  2 --real-        -0.075   0.122  -0.013   0.034  -0.000 del=   0.034
 2  1 --real-        -0.075   0.060   0.002  -0.012  -0.000 del=  -0.012
elda -> E_LDA
e0 -> E_LDA scissor shifted (linear interpolation)
x -> <nk|Gv|nk>, fock operator
sx +ch -> <nk|G(W-v)|nk>
sig = sx + ch
v0 -> <nk|v_xc|nk>
efs0 = e0 + sig - v0
enew = e0 + Z(sig - v0), where Z is close to 1
     Sigma depends on E
E^QP = E_0 + <nk|Sigma(E^QP) - vxc|nk>
 Assuming Sigma(E^QP) = Sigma(E_LDA) + Sigma'(E_LDA)(E^QP-E^LDA)
then we can show that
E^QP = E_0 + Z <nk|Sigma(E^LDA) - vxc|nk>
 where Z = 1/(1-Sigma'(E_LDA))
------------------------------------------------------------------
We calculate
Sigma(E) = G({Psi_LDA},{E_0}) * W({Psi_LDA,{E_0'})
{E_0} is given by E_LDA with scissor shift specified in sig.inp
{E_0'} is given by E_LDA with scissor shift specified in xi0.inp
{E_0} --> some approximation to E^QP.
How scissor shift is defined:

eval=eval+(evs+evdel*(eval-ev0))/
                        ((1.0d0-evdel)*ryd)
econd=econd+(ecs+ecdel*(econd-ec0))/
                        ((1.0d0-ecdel)*ryd)

Updating Xi : epsilon -> 1 + wp^2/E_av^2
doesn't matter.         E_av = average gap (which doesn't change much)

5.7 Visualization

By setting the appropriate output flags, various quantities can be visualized. Paratec can produce line plots (filename ??.line.number), 2D slice plots and files for the IBM Data Explorer (filename ??.dx). If you are using Khoroshttp://www.khoros.unm.edu/ instead of the Data Explorer, set the output flags khoros, which will produce a file with suffix .kh.

5.7.1 Tensor Data

Often, you want to plot only a certain component of a field or tensor field. You use a statement like:

tensor_bra 0 1 0
tensor_ket 0 0 1

to select the yz component (in cartesian coordinates) of a tensor. The “tensor” input is also used for specifying the alignment of the external B-field for plotting purposes. Default:

tensor_bra 1 0 0
tensor_ket 1 0 0

5.7.2 Line plots

The line plot structure determines along which lines there should be plotted:

Example:

begin line_plot
line 0.0 0.0 0.0  1.0 1.0 1.0   100    plot data along (111) on 100 mesh
points
end line_plot

The start and end coordinates of line plot are relative to the realspace lattice vectors.

5.7.3 Slice plots

The slice plot structure sets the planes on which data should be plotted. To define a plane specify the lower left, lower right and top left points. The number of grid points in the horizontal and vertical directions are also needed.

Example:

begin slice_plot
slice 0.0 0.0 0.0  1.0 0.0 0.0   1.0 0.0 1.0  50    100
end slice_plot
  (for a cubic lattice) plots data on the xz plane with
  50 grid points in the x direction and 100 in the z

The start and end coordinates of slice plot are relative to the realspace lattice vectors.

5.7.4 Visualizing wave functions

Paratec can produce Data Explorer and Khoroshttp://www.khoros.unm.edu/ files that visualize wavefunctions at particular kpoints, bands. The square of the wavefunctions is plotted in realspace. Wave function plotting is switched on

by setting pw job scf and output flags waveplot.

The wave functions to be plotted are specified as pairs of the irreducible kpoint number and band number.

Example:

begin plot_wave_function
wavefunction 1 12                  irreducible kpoint 1, band 12
wavefunction 2 10                  irreducible kpoint 2, band 10
end plot_wave_function

5.7.5 Visualizing the charge density

To produce a simple three-dimensional charge density plot, one simply adds the flag cdplot to the list of output flags, and runs an scf calculation. If an energy window is specified with the energy window structure:

begin energy_window
window (start_energy in eV) (end_energy in eV)
end energy_window

then the charge density for the plot is generated only by using wave functions with eigenvalues within the energy window (specified in eV). To find suitable values for start energy and end energy, look at the eigenvalues printed out in previous runs. A gaussian smearing is applied, and the occupation numbers are ignored. This feature only works if the wave functions have been computed, i.e. it does not work with the pot plot plane wave job.

Warning: the energy window specification does not include an adjustment for the fermi level. Also, unoccupied bands will be included in the charge density if your energy window includes energies above the fermi level.

5.7.6 Plotting the potential

Both line plotting and three-dimensional plotting of the effective Kohn-Sham potential (without the nonlocal pseudopotential part) are available via the pot plot plane-wave job. For more see section 4.0.5.

5.7.7 Ball and stick model

A ball-and-stick file BALLS.dx or BALLS.kh can be generated and viewed with the IBM Data Explorer or Khoroshttp://www.khoros.unm.edu/. This is done by adding the output flag ballnstick to the output flags field (section 5.10). To use the data explorer, please refer to the README file in /usr/local/codes/paratec/dx. To use Khoroshttp://www.khoros.unm.edu/ instead, set output flags khoros.

5.8 Spin polarized calculations and magnetic field

The number of spins to be used can be set with the variable number of spins:

   number of spins 1

(for non-spin polarized LDA)

or

   number of spins 2

(for spin-polarized LSDA)

   Default: number of spins 1

To get a spin-polarized starting configuration for the self-consistent charge, use the spin polarization factors where you input the crystal structure.

5.8.1 Applying a magnetic field

A magnetic field

H = H_{0} cos(G⃗ · r⃗)  (5.2)

can be applied to stabilize for example an antiferromagnetic phase. The magnetic field acts only on the spin, AND ALTERS THE ENERGY AND WAVEFUNCTIONS!! This is done with an entry like:

   magnetic field kvec 1 1 0 1

The last number switches the field on (=1) or off (=0), the first 3 INTEGER numbers specify the wave number in reciprocal lattice coordinates.

The amplitude (normalized such that μ_{b}H_{0} = 1 Ryd) is entered as:

   magnetic field strength 0.001

5.9 Reading additional input

The option input flags enables you to read additional input from various files. The flags for the various inputs you want to read are concatenated by underscore.

input flags flag1 flag2 flag3 (sets 3 input flags)

Implemented flags:

• wavefn reads wave functions from a previous calculation as a starting guess for the iterative diagonalization. It is generally safer not to use this option, because the eigenvectors can be really bad if the structure is different, causing the diagonalization routine to bomb out.

• chkptchi reads the intermediate results from a previous nmr shift calculation. The data must be in a file CHKPTCHI.x, where x is the number of the processor. Make sure you don’t read in garbage from the checkpoint file!

5.10 Getting additional output

The option output flags enables you to get additional numbers out of the code. The flags for the various outputs you want to obtain are concatenated by underscore.

output flags flag1 flag2 flag3 (sets 3 output flags)

Implemented flags:

• diagperf gives the performance of the diagonalization routine, resolved into the FFT part and the matrix-matrix multiplies.

• vqmc writes out a data file for the variational quantum monte carlo code.

• nmrshift plots the diamagnetic susceptibility χ in realspace. If requested, line plots are done also.

• timing prints timing information for various parts of the code. The time printed is system and user time together, per processor.

• memusage prints out memory consumption of large arrays in various subroutines. Not all subroutines print out their memory usage currently.

• cdplot produces charge density plots after a self-consistent field calculation.

• neighbors prints nearest distances in a.u. between the atoms in the unit cell.

• ballnsticks generates a file for DX or Khoroshttp://www.khoros.unm.edu/ which contains the positions of the balls.

• ballswrapped same as ballnsticks, but all balls are within the first unit cell. This is useful for volumetric plots.

• angles prints angles between atoms, e.g. bond angles.

• dos computes the density of states with the tetrahedron method. The output goes to the file “DOS”. The job variable has to be set to pw job scf. You can inhibit the selfconsistent iterations my setting max iter scfloop

1.

• angdos Same as dos, but also computes the local and angular momentum resolved density of states for each atom. The output goes to the file “DOS”. The job variable has to be set to pw job scf. You can inhibit the selfconsistent iterations my setting max iter scfloop 1.

• project Outputs the angular projections of the wave functions. Currently this calls two separate codes, one of which supports f-states, while the other supports projection onto multiple atoms. One outputs to ’project.dat’ while the other outputs to ’projections’.

• waveplot Writes the square modulus of particular wave functions into a .dx or .kh file.

• potplot Generates a graphical output file of the Kohn-Sham effective potential (without the nonlocal part) into a .dx or .kh file. In order for this to work, you must also specify a pw job pot plot. For more information, see Section 4.0.5.

• eigvals Writes the eigenvalues, the weights of the kpoints, and (if computed) the angular momentum decomposition for each kpoint and band into a file “EIG”

• efield Compute the electric field due to the ionic potential and the charge density. In order for this to work, you must also specify a pw job pot plot. For more information, see Section 4.0.5.

• wavefn tells the code to write a file with the eigenvectors at the end of a self-consistent field calculation (pw job scf).

WARNING: paratec names the wavefn files using only one digit for the k points, so if you plot the same band at multiple k points (which have an index greater than 9), the wavefunctions will overwrite all but the first. This is a bug, but you can work around it, and I (David Roundy) don’t feel like fixing it right now. Sorry.

• momdens write a file “MOMDENS” with the momentum density information. See section 5.5.

• khoros write all graphical output files for Khoroshttp://www.khoros.unm.edu/ (.kh) rather than DX (.dx).

• xyz Output structural data in xyz format, for use in the xmol visualization program.

• diagconv Gives the keyword rho then the residualˆ2, iteration number, and the trace of the Hamiltonian.

• fermisurf Outputs the eigenvalues in k space, so you can plot the fermi surface.

Funky extra flags

6.1 Optimization flags

optimize allows you to optimize for different resources, e.g. disk IO, memory usage, etc. The flags for the various optimizations you want to obtain are concatenated by underscores.

optimize flag1 flag2 flag3 (sets 3 optimization flags)

Implemented flags:

• insulator does a direct minimization of the energy functional rather than the usual scfloop. Only works for insulators! Also, the number of bands must be ONLY THE NUMBER OF OCCUPIED BANDS.

• memory optimizes for memory in the calculation of NMR shifts, at a small expense in additional CPU time.

• nocgsolveguess switches OFF the better starting guesses for the conjugate gradient solver in the NMR subroutine. Using this option will generally result in SUBSTANTIALLY SLOWER execution, but will reduce memory cost.

• nowfnreuse when relaxing with a fixed lattice, DOES NOT reuse the wave functions from the previous relaxation step.

• noprecond do not precondition in the direct energy minimization. Use this option only if you have convergence problems, since it will generally result in 2-3 times slower convergence.

• noadjustshift keeps the number of bands fixed during the calculation and should be used if you have trouble with complex charge densities.

The computational cost is generally dominated by the FFTs, and on parallel computers, a significant fraction of the FFT time is involved in communication between processors. Poor scaling results when latency dominates these communications, but the latency may be reduced (with an associated memory cost on some machines) by setting

number bands fft n

where n is the number of bands to FFT simultaneously (this option is disabled when using the SHMEM library).

Helpful tools and stuff

7.1 Paratec Perl Library

The paratec perl library is currently under construction. It is a set of perl modules designed to make it much easier to parse and modify paratec input and output files. Actually, it is designed to parse paratec input and output files, and only modify input files! ¨⌣

If you want to use these modules in your perl scripts, you just have to set the environment variable $PERL5LIB to be the path to the paratec/src/perllib directory (so that perl will know where to find the modules), and then follow the instructions below to actually use the code.

Here are the manuals for the paratec modules, which are automatically generated from Plain Old Documentation embedded in the module itself:

• input.pminput.html

• OUT.pmOUT.html

• PW LOG.pmPW_{L}OG.html The library also has a few utility modules to make it easier to deal with things like vectors and tensors. Here is the documentation to the utility files:

• vector.pmvector.html

• tensor.pmtensor.html

7.2 tools

7.2.1 gaussdos

This program uses the EIG file, which is produced when the output flag eigvals is set. It computes the density of states with a gaussian smearing. The input file has the following form:

0.2      gaussian smearing in eV
1000     number of mesh points for DOS

7.2.2 bsfix

In case you (or paratec!) screwed up with the fermi level for the band structure plot, you can use this awk script to shift the bandstructure plot around. Just look at the script, modify the fermi energy there, and run it with

bsfix >NEWBANDSTRUC

It will read the file BANDSTRUC, shift the bands according to the new fermi level, and write the output to NEWBANDSTRUC.

Erratum etc.

8.1 Known difficulties

1. Matrix diagonalization bombs out. Increase diagsafety (section 4.0.8).

2. The double unit cell problems: If you run a supercell of twice the size of the real cell, i.e. if there are NONPRIMITIVE translations within the unit cell left, you will encounter complex charge density problems. This is due to shortcomings in the routine that finds the symmetry operations. Workaround: Only use the primitive unit cell, or break the symmetry completely by slightly perturbing the atomic positions. Alternatively, use

the noadjustshift option to the optimize entry.

3. Loss of symmetry: After many structural relaxation steps (say 15 or so), roundoff errors break the symmetry of the system, and lead to strange results. Often, the forces are symmetrized incorrectly, leading to termination of the code. Workaround: Restart from last step after explicitly symmetrizing the structural parameters.

4. Bad charge density: If the program bombs out, it is often because a bad initial charge density has been used. There is no sanity checks being done on the CD file! Workaround: Remove old charge density (the CD file), and retry.

5. Structural relaxation scheme does not converge: There can be several reasons. Often, the system is about to undergo a phase transition, and is just intrinsically instable. Another possibility is that the accuracy of the diagonalization is not good enough (see the corresponding section), or the potential convergence criterion is too large.

6. Negative charge density. Should only appear when the initial charge is set up from the atomic charge density, or if a nonlinear core correction is used.

7. Complex charge density. Most of the times, this is because the coordinates of the atoms somewhat match the symmetry (say to 1e-8), but not quite. Fix: input lattice parameters and basis coordinates with very high accuracy, e.g. 1e-20. The gnu utility “bc” is very useful for this case, as it has arbitrary precision arithmetic.

8. NMR q problem: If the perturbation nmr q for the NMR calculation is set too small, you will get wrong results unless the accuracy for the wave function is also cranked up. This affects mostly the susceptibility (G=0). An accuracy of 1e-12 for the wave function and a nmr q of 0.01 is normally a good choice. In doubt, check if results are stable when nmr q is varied.