HELP FOR INPUT FILE FOR BIOSYM/CORNING PLANEWAVE CODE
This document explains the i/o parameters and format needed
for version 940205 (05 Feb 1994) of the
Corning electronic structure code (plane_wave)
Copyright (C) Douglas C. Allan
========================================
This code, in conjunction with "insightII" (the Biosym interface)
has been commercialized. This file is intended to be helpful in
using the code as a stand-alone in which the user creates the input
file with any editor, not with the Biosym interface. Biosym provides
documentation to users who utilize the commercial version.
The main executable, now called plane_wave, uses dynamic allocation
to free the user from specifying array dimensions. An estimate of
the job size is given in the log file which is helpful when your
job slightly exceeds machine memory.
The executable "newsp" does not yet have dynamic allocation so
users may need to redimension arrays within the fortran source,
which is provided in the file "newsp_wrap.f".
Given an input file (parameters described below) and the required
pseudopotential files, the user must also create a "files file"
which lists the names of all the files the job will require, including
the input file. The files file could look like:
input
output
wf.input
wf.output
14si.psp
in which the input file is called "input", the output will be put
into the file called "output", the input wavefunctions (if any) are
in "wf.input", and the output wavefunctions will be written to
"wf.output". Finally the needed pseudopotential for this job is
"14si.psp". Given these files, the stand-alone version of plane_wave
is run with the command
plane_wave none input files.file none 1 >& log &
where again the input file is called "input", the files file is
called "files.file", and standard out and standard error are piped
to the log file called "log". The user can specify any names he/she
wishes for any of these files. Arguments 1 and 4 are ignored, and
argument 5 must be "1" when plane_wave is used as a stand-alone.
========================================
The methods employed in this computer code to solve the
electronic structure problem are described in part in
``solution of schrodinger's equation for large systems'',
Phys. Rev. B {40}, 12255-12263 (1989) by
Michael P. Teter, Michael C. Payne, and Douglas C. Allan.
A more complete description is available in a Reviews
of Modern Physics article, ``Iterative minimization techniques
for ab initio total-energy calculations: molecular dynamics
and conjugate gradients'', M. C. Payne, M. P. Teter,
D. C. Allan, T. A. Arias, and J. D. Joannopoulos,
Rev. Mod. Phys. 64, 1045-1097 (1992).
This document lists the input parameters and file names to
be provided, and describes the output of the code.
It also describes other aspects of running the code, such as
other parameters which may affect the quality of results and
the use of several programs to prepare the input file.
This code has been written with great care to avoid using
fortran features which are not supported by the American
National Standard Institute (ANSI) X3.9-1978, i.e. it is written
to conform with ANSI standard fortran (Fortran 77).
The fortran NAMELIST feature is not supported by ANSI and so is
avoided. In its place the code uses a set of parsing subroutines
which provide for a namelist-style input file containing all
control parameters.
The parameters are input to the code from a single input file.
Each parameter value is provided by giving the name of the parameter
and then placing the numerical value(s) beside the name, separated by
one or more spaces or tabs. The parameter may be an integer or a
floating point number, and may actually represent an array of values.
If it represents an array, the next set of numbers separated by spaces
are taken as the values for the array. More description is given
below.
The names and definitions of all the parameters are given below
in alphabetical order. In the actual input file, these parameters
may be given in any order desired and more than one may be given
per line. Spaces or tabs are used to separate values and additional
spaces or tabs are ignored.
An as example of input, the parameter for length scales is called
"acell" and is an array acell(3) for the lengths of the primitive
translations in bohr atomic units. To input a typical Si diamond
lattice one would have the line
acell 10.25311 10.25311 10.25311
in the input file. This may equivalently be written
acell 3*10.25311
and will still be parsed correctly.
Multiple spaces are ignored, as is any text which does not contain
the character strings which correspond to input parameters.
A given variable is identified by the parser by having at least
one blank or tab before it and after it (again, multiple blanks
or tabs are irrelevant). Tabs are simply converted to one blank
per tab before parsing and are therefore equivalent to blanks.
To include comments it is recommended that they be placed to the
right of the comment character #; anything to the right of a "#"
on any line is simply ignored by the parser. Comments would not
otherwise cause trouble unless they inadvertantly contained
character strings which were the same as variable names
(e.g. " acell ").
The character "#" can also be used to "store" old values of variables
or place anything else of convenience into the file in such a way as
to be ignored by the parser when the data is read.
Case is irrelevant as the entire input string is mapped to
upper case before parsing, to remove case sensitivity.
More than one parameter per line may be given.
If a given parameter name is given more than once in the input
file, the code reads the value which occurs first in the file
and ignores all the others.
Below are the parameters in alphabetical order with a
description of their meaning.
acell real array acell(3) giving length scales by which
dimensionless primitive translations (in "rprim") are
to be multiplied. Given in bohr atomic units
(1 bohr=0.529177249 Angstroms). See description of rprim.
amu real array amu(ntype) giving masses in amu for each kind
of atom in cell. These masses are used in conducting
molecular dynamical atomic motion if ionmov=1.
Note that one may set all masses to 1 for certain cases
in which merely structural relaxation is desired and not
actual molecular dynamics.
Using 1986 recommended values, 1 amu = 1.6605402e-27 kg.
In these units the mass of Carbon 12 in amu is exactly 12.
********
Input variable apwsph is presently not supported--do not use
********
apwsph real array apwsph(3) giving length scales "acell" for
"planewave sphere", i.e. giving the length scales
actually used in conjunction with "ecut" to determine
the basis of planewaves used. Default is "acell" if not
given, and should equal "acell" unless made not equal to
acell for special reasons. apwsph is introduced to allow
for independent control over the planewave basis and the
unit cell dimensions.
asr integer parameter controlling evaluation of the
acoustic sum rule in effective charge calculations
within linear response theory.
0 => no acoustic sum rule imposed
1 => acoustic sum rule imposed "democratically" with
extra charge evenly distributed among atoms
2 => acoustic sum rule imposed "aristocratically" with
extra charge given proportionally to those atoms with
the largest effective charge
Default is 0.
charge real number used to establish charge balance between
the number of electrons filling the bands and the
nominal charge associated with the atomic cores.
The code adds up the number of atoms of each type
times zion for each type (call this the core charge).
zion is obtained for each type of atom from the psp files.
Then the code adds up the band occupancies (given in
array occ) for all bands at each k point, then multiplies
by the k point weight wtk at each k point. Call this
sum nelect for the number of electrons. It is then
required that (core charge) - nelect = charge,
where charge is provided on input. To treat a neutral
system, which is desired in nearly all cases, one must
use charge=0. To treat a system missing one electron
per unit cell, set charge=+1. The default is charge=0.
Irrelevant if iscf=0.
date six digit integer parameter giving the date in the
form yymmdd where yy=year number (00-99)
where year number is assumed to refer to
1900-1999 (I will worry about 2000 in a future release),
mm=month number (1-12), and dd=day number (1-31 or less).
This parameter is provided in the input file and
echoed in the output file for the convenience of the
user in keeping track of different jobs. It has to
be reset by hand in the input file to avoid a machine
dependent "date" call. Default is the version date
of the input/output code.
dedlnn real parameter giving a value for derivative
d(Etotal)/d(log(Npw)) for given value of ecut.
This value must be computed independently by making
several runs at fixed geometry and variable ecut,
generally within +/- 3% of the desired ecut,
and using the Etotal(npw) data to compute the
derivative. See enefit.f for a convenient fitting
routine. dedlnn is used to compute the Pulay correction
to the stress tensor using correction=(1/ucvol)*dedlnn.
Default dedlnn is 0, i.e. no correction.
dtion real parameter for controlling ion time steps.
If ionmov is set to 1 then molecular dynamics is
used to update atomic positions in response to
forces. The parameter dtion is a time step in
atomic units of time. (One atomic time unit is
2.418884e-17 seconds, which is the value of
Planck's constant in hartree*sec.)
In this case the atomic masses, in amu (given in array "amu"),
are used in Newton's equation and the viscosity and number
of time steps are provided to the code using input variables
"vis" and "ntime". The code actually converts
from masses in amu to masses in atomic units (in units
of electron masses) but the user enters masses in amu.
(The conversion from amu to atomic units (electron
masses) is 1822.88851 electron masses/amu.)
A typical good value for dtion is about 100.
This is the default value. The user must try several values
for dtion in order to establish the stable and efficient
choice for the accompanying amu, atom types and positions,
and vis (viscosity).
eatpol integer array of 2 elements which controls the range
of polarizations computed in phonon calculations.
These values are only relevant to phonon response function
calculations using the response function code.
May take values from 1 to natom, with eatpol(1)<=eatpol(2).
Default values are 1,1.
(respfn only)
ecut real parameter for kinetic energy cutoff in Hartree atomic
units (1 Ha=27.2113961 eV) which controls number
of planewaves at given k point by
(1/2)[(2 Pi)*(k+Gmax)]^2=ecut for Gmax.
All planewaves inside this "basis sphere" centered
at k are included in the basis.
This is the single parameter which can have an enormous
effect on the quality of a calculation; basically the larger
ecut is, the better converged the calculation is. For fixed
geometry, the total energy MUST always decrease as ecut is
raised because of the variational nature of the problem.
Usually one runs at least several calculations at various ecut
to investigate the convergence needed for reliable results.
enunit integer parameter: 0=>print eigenvalues in hartree;
1=>print eigenvalues in eV; 2=>print eigenvalues in
both hartree and eV. Default is 0 (hartree only).
ftoler real parameter setting a root mean square force tolerance
(in hartree/bohr) below which molecular dynamics or
BFGS structural relaxation iterations will stop.
See the parameter ionmov.
This is to be used when trying to equilibrate a
structure to its lowest energy configuration. A value of
about 0.002 hartree/bohr or smaller is suggested (this would
correspond to 0.103 eV/Angstrom). Default is 0.
intxc integer parameter:
0=> do "usual" xc quadrature on fft grid
1=> do higher accuracy xc quadrature using fft grid
and additional points at the centers of each cube
(doubles number of grid points)--the high accuracy version
is only valid for boxcut>=2. If boxcut < 2, then
code sets intxc=0 and gives warning.
(Note: currently disabled, 24 Jun 1993.)
ionmov integer parameter:
0=> do not move ions;
1=> move atoms using molecular dynamics with
optional viscous damping (friction linearly proportional
to velocity). The viscous damping is controlled by the
parameter "vis".
If actual undamped molecular dynamics is desired,
set vis to 0.
2=> conduct structural optimization using the
Broyden-Fletcher-Goldfarb-Shanno minimization (BFGS).
This is much more efficient for structural optimization
than viscous damping.
3=> Vanderbilt-Louie structural optimization
(not well tested yet)
Default for ionmov is 0.
ireadc integer parameter:
0=>initialize wavefunctions with random numbers for
ground state code; initialize first order wavefunctions
to 0's for response function code
1=>read ground state wavefunctions from a disk file produced
in a previous run for ground state code; read first order
wavefunctions from a disk file for response function code
Default is 0.
iscf integer parameter controlling self-consistency.
0 =>
non-self-consistent calculation is to be done;
in this case a fixed potential on a real space grid
(produced in a previous calculation) is read from a
disk file. Name of file must be given as indicated
below. iscf=0 would be used for
band structure calculations, to permit computation of
the eigenvalues of occupied and unoccupied states at
arbitrary k points in the fixed self consistent potential
produced by some integration grid of k points.
At this time the algorithm coded does not necessarily
converge if unoccupied states are included within
self-consistent iterations. To compute the eigenvalues
(and wavefunctions) of unoccupied states, the user should
save the self-consistent V(r) computed using only occupied
states and then run iscf=0 for the unoccupied states.
1 =>
perform usual self-consistent calculation.
(Present implementation accepts any iscf > 0.)
This is the usual choice for doing the usual ground state
calculations or for structural relaxation.
Other options:
-1 =>
like iscf=0 (non-selfconsistent) but reads occ and wtk, as
needed in Levine's optical response code.
-2 =>
again non-selfconsistent but reads electron density rho(r)
in real space instead of potential V(r). The user may
work exclusively with density files instead of potential
files for conducting band structure calculations.
A density file can be produced from a wavefunction file
using the rhofwf routine, or may be created with
the parameter prtrho (see its description).
-3 =>
like -2 but reads occ and wtk (for Levine).
4 =>
in Gonze's response function code, do response calculations
without local field corrections.
Default is 1.
ixc integer parameter controlling choice of exchange and
correlation (xc).
0=> NO xc;
1=> Teter rational polynomial which reproduces Perdew-Wang
which reproduces Ceperley-Alder (with spin-pol if wanted).
Others are not supported now but could be put back in
fairly easily (were in earlier versions).
Note that the choice made here should agree with the choice
made in generating the original pseudopotential, except
for ixc=0 (usually only used for debugging).
Default is ixc=1 (Teter parameterization).
kpt real array kpt(3,nkpt) containing k points in terms
of reciprocal space primitive translations (NOT in
cartesian coordinates!). Thus the current kpt contains
dimensionless numbers in terms of which the cartesian
coordinates would be
k_cartesian = k1*G1+k2*G2+k3*G3
where (k1,k2,k3) represent the dimensionless "reduced"
coordinates and G1, G2, G3 are the cartesian coordinates
of the primitive translation vectors. G1,G2,G3 are related
to the choice of direct space primitive translation vectors
made in "rprim" below. Note that an overall norm for the k
points is supplied by the next argument. This allows
one to avoid supplying many digits for the k points to
represent such points as (1,1,1)/3.
Also note: one of the algorithms used to set up the sphere
of G vectors for the basis assumes that the input kpt
lies inside the first Brillouin zone. This will cause
unpredictable trouble if an input kpt lies outside the
first zone, so until further notice please confine values
of kpt to those lying inside the first zone. The code should
correctly handle components in the range [-1,1], so the
remapping is easily done by adding or subtracting 1 from
each component until it is in the range [-1,1]. That is,
given the k point normalization kptnrm described below,
each component must lie in [-kptnrm,kptnrm].
kptnrm real parameter established a normalizing denominator
for each k point. The k point coordinates as fractions
of reciprocal lattice translations are therefore
kpt(mu,ikpt)/kptnrm. kptnrm defaults to 1 and can
be ignored by the user. It is introduced to avoid
the need for many digits in representing numbers such as 1/3.
machin integer parameter telling machine on which the
code is being run. At present there are not machine
dependencies which rely on "machin" but they are
anticipated in the future. The present assignment is
machin machine
0 any (default)
1 ibm 3090 with CMS
2 IBM rs/6000
3 Cray YMP
natom integer parameter giving number of atoms in unit cell.
Default is 1 but you will obviously want to input this
value explicitly.
nband integer parameter giving number of bands, occupied plus
possibly unoccupied, under consideration. Gives total
number of bands for which wavefunctions are being computed
along with eigenvalues. Default is 1.
NOTE: if the parameter occopt (see below) is set to 2,
then nband must be an array nband(nkpt) giving the
number of bands explicitly for each k point. This
option is provided in order to allow the number of
bands treated to vary from k point to k point.
ng integer array ng(3) giving size of fast fourier transform
(fft) grid in three dimensions. Each number must be
composed of the factors 2, 3, and 5 with other possible
restrictions for use on the IBM with ESSL DCFT.
The auxiliary main routine "getng" allows convenient
computation of the recommended values for ng given
any "rprim", "acell" and "ecut". When ng is made smaller
than recommended values, the code runs faster and the
equations in effect are approximated by a low pass fourier
filter. The code reports to standard output (unit 06) a
parameter "boxcut" which is the smallest ratio of the fft
box side to the G vector basis sphere diameter. When
boxcut is less than 2 the fourier filter approximation is being
used. When boxcut gets less than about 1.5 the
approximation may be too severe for realistic results
and should be tested against larger values of ng.
When boxcut is larger than 2, ng could be reduced without
loss of accuracy.
nkpt integer parameter giving number of k points in k point set.
Default is 1.
nline integer parameter giving maximum number of line minimizations
allowed in preconditioned conjugate gradient minimization
for each band. Usually should be set around 4-8. Line
minimizations will be stopped anyway when improvement
gets small. Default is 4.
nqpt integer parameter giving number of q points in phonon
calculation. Presently this may be no larger than 1.
0 means phonons are not being computed. Default is 0.
(respfn only)
nsppol integer parameter which must take the value either
1 or 2. 1 means an unpolarized calculation, 2 means
spin-polarized. The default is 1.
nstep integer parameter giving maximum number of preconditioned
conjugate gradient passes through all the bands at each
k point before current job must end. Job will end sooner
if convergence criteria are met.
Full convergence from random numbers if usually achieved in
10-15 iterations. Each can take from minutes to hours.
In certain difficult cases, usually related to a small or
zero bandgap, convergence performance may be much worse.
Default is 1.
NOTE that a choice of nstep=0 is permitted; this will
either read wavefunctions from disk (with ireadc=1) and
compute the total energy and stop, or else will
initialize the wavefunctions and compute the resulting
total energy. This is provided for testing purposes.
The stress tensor still gets computed with nstep=0.
nsym integer parameter giving number of space group symmetries
to be applied in this problem. Symmetries will be input in
array "symrel" and nonsymmorphic translations will be input
in array "tnons". If there is no symmetry in the problem
then set nsym to 1, because the identity is still a symmetry.
Default is 1.
ntime integer number of molecular dynamics time steps or BFGS
structural optimization steps to be done if ionmov=1 or 2
respectively. Note that at the present
the molecular dynamics is initialized with four
Runge-Kutta steps which costs some overhead in the startup.
ntime is ignored if ionmov=0.
Default is 5.
ntype integer parameter giving number of types of atoms. E.g. for
a homopolar system (e.g. pure Si) ntype is 1. Default is 1.
occ real array occ(nband) giving occupation numbers for all
bands in the problem. Typical band occupancy is either
2 or 0, but can be 1 for half-occupied band or other
choices in special circumstances. If occopt is not 2,
then the occupancies must be the same for each k point.
If occopt=2, then the band occupancies must be
provided explicitly for each band AND EACH k POINT, in an
array which runs over all bands and k points.
The order of entries in the array would correspond to
all bands at the first k point, then all bands at the
second k point, etc. The total number of array elements
which must be provided is nband(1)+nband(2)+...+nband(nkpt).
occ is set to 0's and ignored when iscf=0.
occopt integer option parameter controlling how input
parameters nband, occ, and wtk are handled.
occopt=0:
All k points have the same number of bands
and the same occupancy of bands. nband is given as a
single number, and occ(nband) is an array of nband
elements. The k point weights in array wtk(nkpt) are
automatically normalized by the code to add to 1.
occopt=2:
k points may optionally have different numbers of
bands and different occupancies. nband(nkpt) is given
explicitly as an array of nkpt elements.
occ() is given explicitly for all bands at each k point--
the total number of elements is the sum of nband(ikpt)
over all k points. The k point weights wtk(nkpt) are
NOT automatically normalized under this option.
Other options, such as automatic band filling without
need of specifying occ, and Fermi function band filling
for metals, are envisioned for the future.
The default is occopt=0.
prtpot Provide output of local potential, sum of local pseudo-
potential, Hartree, and xc potentials if set >=1.
Present version only gives output for static (ionmov=0)
calculation. Default is 0.
prtrho Provide output of electron density in real space rhor
on every time step of molecular dynamics. Output is to
file rhor.n in units of electrons/bohr^3, where "n" starts
at the input value of prtrho and gets incremented by 1
on each step. No output is provided by prtrho <= 0.
The file structure of the unformatted output file is
described below, where output is described.
Default for prtrho is 0.
prtvol integer parameter controlling volume of printed output.
Present choices are 0 or >0; a setting of prtvol>0
causes the eigenvalues to be printed for every iteration.
Default is 0.
qpt real array of 3 elements containing the q vector of
a phonon mode. See qptnrm for extra normalization.
(respfn only)
qptnrm real parameter which provides normalization of qpt.
Default is 1. Actual q vector is qpt/qptnrm.
(respfn only)
rfdir integer array of 3 elements which gives the directions for
response function calculations. Default is 0 0 0.
(respfn only)
rfmeth integer parameter which selects method used in
response function calculations. Default is 1.
(respfn only)
rfelfd integer parameter to turn on electric field response
function calculations. Default is 0.
(respfn only)
rfphon integer parameter to run phonon response function
calculations. Default is 0.
(respfn only)
rfstrs integer parameter to run stress response function
calculations (elastic constants). Default is 0.
(respfn only)
rf3 integer parameter to control response function calculation
of third order response. Default is 0.
(respfn only)
rprim real array rprim(3,3) giving, in columnwise entry, the three
dimensionless primitive translations in real space.
For example rprim 1 2 3 4 5 6 7 8 9 corresponds to input
of the three primitive translations R1=(1,2,3),
R2=(4,5,6), and R3=(7,8,9).
Each column of rprim gets multiplied by the appropriate length
scale acell(1), acell(2), or acell(3) respectively to give the
dimensional primitive translations in real space in
cartesian coordinates. Note carefully that the first
three numbers input are the first column of rprim, the next
three are the second, and the final three are the third.
This corresponds with the usual Fortran order for arrays.
The matrix whose columns are the reciprocal space primitive
translations is the inverse transpose of the matrix whose
columns are the direct space primitive translations.
The dimensional cartesian coordinates of the crystal
primitive translations Rp are
R1p(i)=rprim(i,1)*acell(1) for i=1,2,3 (x,y,and z)
R2p(i)=rprim(i,2)*acell(2) for i=1,2,3
and R3p(i)=rprim(i,3)*acell(3) for i=1,2,3.
sciss value of "scissors operator" shift of conduction band
eigenvalues (hartree), used in response function
calculations. Default is 0.
(respfn only)
symrel integer array symrel(3,3,nsym) giving "nsym" 3x3 matrices
expressing space group symmetries in terms of their action
on the direct (or real) space primitive translations.
It turns out that these can always be expressed as integers.
The name is a mnemonic for "symmetry in real space".
Always give the identity matrix even if no other symmetries
hold, e.g.
symrel 1 0 0 0 1 0 0 0 1
Also note that for this array as for all others the array
elements are filled in a columnwise order as is usual for
Fortran.
The relation between the above symmetry matrices symrel,
expressed in the basis of primitive translations, and
the same symmetry matrices expressed in cartesian coordinates,
is as follows. Denote the matrix whose columns are the
primitive translations as R, and denote the cartesian
symmetry matrix as S. Then
symrel = R(inverse) * S * R
where matrix multiplication is implied.
tn real array tn(3,natom) giving atomic locations within
unit cell in coordinates relative to real space primitive
translations (NOT in cartesian coordinates). Thus these
are fractional numbers typically between 0 and 1 and
are dimensionless. The cartesian coordinates of atoms
are given by
t_cartesian = t1*r1*a1+t2*r2*a2+t3*r3*a3
where (t1,t2,t3) are the reduced coordinates given in
columns of "tn", (r1,r2,r3) are the columns of
dimensionless array "rprim", and (a1,a2,a3) are the
elements of the array "acell" giving length scales in bohr.
If you prefer to work only with cartesian coordinates, you
may work entirely with "xcart" (below) and ignore tn, in
which case tn must be absent from the input file.
tnons real array tnons(3,nsym) giving nonsymmorphic translations
associated with the symmetries expressed in "symrel".
These may all be 0, or may be fractional (nonprimitive)
translations expressed relative to the real space
primitive translations. If the space group
is symmorphic then these are all 0.
type integer array type(natom) giving an integer label to
every atom in the unit cell to denote its type.
This array is restricted to start at the value 1
and be nondecreasing, i.e. atoms of a given type
are grouped together for my convenience in writing
the program. The array "type" has to agree with the
actual locations of atoms given in "tn" and the input
of pseudopotentials has to be ordered to agree with the
atoms identified in "type". The atomic numbers of the
atoms, given by array "zatnum" below, also must agree with
the type of atoms designated in "type".
useri This is a user-definable integer which the user may
input and then utilize in subroutines of his/her own
design. It is not used in the Corning code except as
modified by users. Default value is 0 .
userr This is a user-definable real number with the same purpose
as useri above. Default value is 0.0 .
vel real array vel(3,natom) giving starting velocities
of atoms, in cartesian coordinates, in bohr/atomic time
units (atomic time units given where dtion is described).
Irrelevant unless ionmov > 0.
Default is 3*natom 0's.
vis real parameter giving viscosity (atomic units) for linear
frictional damping term applied to molecular dynamics
when ionmov=1. Used for relaxing structures.
An atomic unit of viscosity is measured in
hartrees*(atomic time units)/bohr^2. Units are not
critical as this is a ficitious damping used to relax
structures. A typical value for silicon is 400 with
dtion of 350 and atomic mass 28 amu. Critical
damping is most desirable and is found only by
optimizing vis for a given situation. Default is 100.
wftol real parameter giving a convergence tolerance for the
largest squared "residual" (defined below) for any
given band. The squared residual is
<nk|(H-E)^2|nk>, E=<nk|H|nk>,
which clearly is nonnegative and goes to 0 as
the iterations converge to an eigenstate.
With the squared residual expressed in
Hartrees^2 (Hartrees squared), the largest squared
residual (called residm) encountered over all bands
and k points must
be less than wftol for iterations to halt due to
successful convergence. This also imposes a restriction
on taking an ion step; ion steps are not permitted
unless the largest squared residual is less than
wftol, ensuring accurate forces. A typical setting
for fair convergence is 10^-10 to 10^-12, with good
convergence from about 10^-14 to 10^-22.
I have found about 10^-12 to work out well unless
I am demanding superconvergence in the wavefunctions
for some particular purpose. Default is 1.d-10.
Note that for approximate structural optimization
using vis>0, wftol may be set relatively larger, in
the range of 1.d-05. This should be tested in any
given system. Also note that machine precision of
about 2.22e-16 permits cutting the residual to about
6.e-31 but no lower. This might be slightly improvable
with a better matrix diagonalizer for subspace
diagonalization.
wtk real array wtk(nkpt) giving k point weights. The
k point weights will have their sum normalized to 1
(unless occopt=2; see description of occopt)
within the program and therefore may be input with any
arbitrary normalization. This feature helps avoid the
need for many digits in representing fractional weights
such as 1/3.
wtk is set to 0's and ignored if iscf=0.
xcart real array xcart(3,natom) giving cartesian coordinates
of atomic within unit cell,in bohr. This information is
redundant with that supplied by array tn and is provided
in order that xcart be available on output for general
information and for molecular dynamics runs.
If tn is ABSENT from the input file and xcart is
provided, then the values of tn will be computed from
the provided xcart (i.e. the user may use xcart instead
of tn to provide starting coordinates).
Default is computation of xcart from input values of tn
if tn is provided, or else computation of tn from xcart
if tn is not provided. If both are provided, tn wins.
zatnum real array zatnum(ntype) giving atomic number for each
type of atom. Ordering of types should agree with
sequence of integers given in "type" and with
reading of pseudopotential files for each atom.
Also obviously must agree with intended atom locations
given in "tn". If "zatnum" does not agree with atomic
number as given in pseudopotential file, program writes
error message and stops.
Additional keyword:
Placing the word "exit" on the top line will cause the job to
end smoothly on the very next iteration. This functions because
the program closes and reopens the input file on every iteration
and checks the top line for the keyword "exit". The word must be
placed with spaces (blanks) on both sides. Thus placing exit on the top
line of the input file WHILE THE JOB IS ALREADY RUNNING will force
the job to end smoothly on the very next iteration.
This feature does not work on the ibm 3090 running cms.
------------------------------------------------------------------------
File names to be provided
In addition to the input file, there are a set of filenames
which must be provided to the code when it runs. These are listed
below. They may be provided from unit 05 interactively during
the run but are more typically provided by piping from a file
in unix or from within a command file or exec in vms or cms.
example use of file
input filename of file containing the input data described
above
output filename of file in which formatted output will be placed.
Error messages and other diagnostics also are sent to
unit 06; the unit 06 output can be ignored unless something
goes wrong. The code repeats a lot of information to both
unit 06 (terminal or log file) and to the named output file.
The unit 06 output is intended to be discarded if the run
completes successfully, with the named output file keeping
the record of the run in a nicer looking format.
wfinput filename of file containing input wavefunction coefficients
created from an earlier run. Must be given as a dummy name
even if no earlier wavefunction file exists. Will only be
opened and read if ireadc is 1 in "input" file. The
wavefunction files are unformatted and can be very large.
wfoutput filename of file in which unformatted wavefunction
coefficients are output. Usually want this name to
resemble the name of "output" file.
psp1 filename of first pseudopotential input file. The
pseudopotential data files are formatted. There must be
as many filenames provided sequentially here as there are
types of atoms in the system, and the order in which the
names are given establishes the identity of the atoms
in the unit cell.
Pseudopotentials may be generated by running a modified
version of Don Hamann's code ATOMPP--the modified version
provides output readable by the corning code.
The modified atompp should be available wherever the
corning code is available.
(psp2, psp3, ...)
(rhorinp) This file name would have to be provided if iscf=-2,
in order to read rho(r) data produced from a previous
run.
------------------------------------------------------------------------
Explanation of the output from the code
Output from the code goes to several places listed below.
(1) Unit 06 (standard output): a log file which echoes the
values of the input parameters and describes various
steps of the calculation, typically in much more detail
than is desired as a permanent record of the run.
This log file is intended to be informative in case of
an error or for a fuller description of the run. For
a successful run the user will generally delete the log
file afterwards.
(2) The output file: a nicely formatted output file to be
kept as the permanent record of the run. It starts with
a heading, followed by the echo of the input data.
Then it reports the real and reciprocal space translation
vectors, then the volume of the unit cell, then the boxcut
ratio, a heading for each pseudopotential
which has been input, and a description of the wavefunction
initialization (random number initialization or input from
a disk file). This includes a report of the number of
planewaves (npw) in the basis at each k point.
The average number of planewaves averaged over all k points
is reported in two forms, arithmetic average and geometric
average. The most appropriate is the geometric average
for use in estimating the Pulay correction to stress or
for making basis set corrections to the total energy
as described in Francis and Payne, J. Phys. C 2, 4395-4404
(1990). See variable dedlnn.
Next the code reports the cpu and wall time for initialization,
and then information for each iteration: the iteration
information is the iteration number, the total energy
(Etot) in hartree, the change in Etot since last iteration
(deltaE), the maximum squared residual residm over all
bands and k points (the residual measures the quality of
convergence and is described with wftol above),
the maximum change in the gradients of Etot with respect
to fractional coordinates (grvar, in hartree), and the
rms value of the gradients of Etot with respect to
fractional coordinates (grsum, in hartree).
If atom motion is selected (ionmov not equal 0) then
in between the above Etot lines are reports of the fractional
coordinates and the gradients of Etot.
After the iterations are completed the squared residuals
for each band are reported, k point by k point.
Then the fractional or reduced coordinates are given,
followed by the energy gradients, followed by
the cartesian coordinates in Angstroms, followed by
the cartesian forces in hartree/bohr and eV/Angstrom.
Also are given the rms force (frms) and the maximum
absolute value of any force component (fmax).
Next are the length scales of the unit cell in bohr and
in Angstroms.
Next are the eigenvalues of each band for each k point,
in eV or hartree or both depending on the choice of enunit.
Next are the components of the total energy broken down into
kinetic, Hartree, exchange and correlation (xc), local
pseudopotential, nonlocal pseudopotential,
local pseudopotential "core correction", and Ewald energies.
Next are separate contributions to the energy gradients:
Ewald, nonlocal pseudopotential, and local pseudopotential.
Next is the stress tensor, (1/ucvol) d(Etot)/d(strain(a,b))
for Etot=total energy per unit cell and a,b are x,y, or z
components of strain. The stress tensor
is given in cartesian coordinates in hartree/bohr^3.
The reported value contains a "Pulay stress" which reflects
the dependence of the basis set on the volume of the unit
cell. This may be estimated by a term on the diagonal of
the stress tensor equal to -(1/ucvol)*d(Etot)/d(log(npw)),
where d(Etot)/d(log(npw)) is estimated from a numerical
derivative at fixed geometry (by varying ecut) and is input
as variable dedlnn. The correction
term is discussed in P. Gomes Dacosta, O. H. Nielsen,
and K. Kunc, J. Phys. C: Solid State Phys. 19, 3163-3172
(1986) and in the Francis-Payne paper cited above. The
correction can be quite significant, especially if the
corrected stress itself is small (as in the case in
which the unit cell is already in equilibrium). To
correct for the Pulay stress, one must SUBTRACT the term
described above, resulting in a correction of NEGATIVE
sign to the calculated stress. (The derivative d(Etot)/d(log(npw))
is always negative, so the Pulay stress is always positive,
so the correction to the analytically calculated stress is
always negative.) With an input value for dedlnn, the correction
is computed within the code and both uncorrected and corrected
stresses are output.
The basics of the stress tensor are described in
O. H. Nielsen and Richard M. Martin, ``Stresses in
semiconductors: ab initio calculations on Si, Ge, and
GaAs'', Phys. Rev. B 32, 3792-3805 (1985).
Next is an echo of the parameters listed in the input file,
using the latest values e.g. for tn, vel, and xcart.
This is followed finally by the timing output: both cpu time
and "wall clock" time as provided by machine-dependent
timing calls within the code.
The total cpu and wall clock times are reported first,
in seconds, minutes, and hours for convenient checking at
a glance.
Next are the cpu and wall times for the principal time consuming
subroutine calls, each of which is independent of the others.
The sum of these times accounts for about 80% of the run
time. These subroutines are
(1) fft: the subroutines fft3d and fft3dp which perform
all the fast fourier transforms;
(2) xcfast: computes the exchange-correlation energy
and potential and sometimes derivative of potential;
(3) nonlop: computes <G|Vnonlocal|C>, the matrix elements
of the nonlocal pseudopotential;
(4) harsum: sums up a Hartree energy;
(5) prjbfa: projects one vector from another, the main
operation performed during orthogonalization;
(6) symmap: calculates the irreducible Brillouin zone
(when symmetry is being used);
(7) ewald: calculates the Ewald energy and its derivatives.
Next are the cpu and wall times for major subroutines which are
NOT independent of each other; for example loopkp conducts
the loop over k points and calls practically everything
else. These subroutines are
(1) loopkp: conducts loop over k points;
(2) hartre: computes Hartree energy and potential
(harsum above does not compute Hartree potential);
(3) getghc: computes <G|H|C>, matrix elements of full
Hamiltonian;
(4) ssdiag: conducts subspace diagonalization;
(5) mkresi: computes residuals <C|H H|C>-<C|H|C>^2
(3) The wavefunction output file: this is an unformatted data
file containing the planewaves coefficients of all the
wavefunctions.
The wf file consists of 4+ntype header lines (described below)
and the following data:
bantot=0 <-- counts over all bands
do isppol=1,nsppol
do ikpt=1,nkpt
write(unit) npw,nband <-- for each k point
write(unit) (eigen(iband+bantot),iband=1,nband)
do iband=1,nband
write(unit) (cg(i+...),i=1,2*npw) <-- single band and k
enddo
bantot=bantot+nband
enddo
enddo
(4) The density output file: this is an unformatted data
file containing the electron density.
It consists of 4+ntype header lines and then
write(unit) (rhor(ir),ir=1,ng(1)*ng(2)*ng(3))
where rhor is the electron density in electrons/bohr^3.
If the calculation is spin-polarized, the total density
is given as above and the spin-up electron density is given
in an additional record of the same form as above.
NOTE: a very advanced user can override the header lines of
the density file by replacing the 4+ntype header lines with
the single line:
write(unit) codvsn, 12345
where "codvsn" is an arbitrary integer (will not be used) and
"12345" is an integer which specifies that no header follows.
To identify the points in real space which correspond
with the index "ir" above, consider the following.
The first array value (ir=1) corresponds with the first
grid point which is at the origin of the unit cell,
(x=0, y=0, z=0). The next grid point (ir=2) lies along the
first primitive translation at the next fft grid point,
which is (1/ng(1))*acell(1)*rprim(mu,1). This is
1/ng(1) of the way along the first primitive translation.
The rest of the values up to ir=ng(1) lie along this
vector, at (ir-1)/ng(1) of the way along the first
primitive translation. The point at ir=ng(1)+1 lies
at 1/ng(2) along the second primitive translation.
The next points up to ir=ng(1)+ng(1) are displaced in the
direction of the second primitive translation by 1/ng(2)
and in the first translation by (ir-ng(1)-1)/ng(1).
This pattern continues until ir=ng(1)*ng(2). The next
point after that is displaced along the third primitive
translation by 1/ng(3), and so forth until ir varies all
the way from 1 to ng(1)*ng(2)*ng(3). This last point is
in the corner diagonally opposite from the origin, or right
alongside the origin if the whole grid is viewed as being
periodically repeated.
(5) The local potential file: also an unformatted file consisting
of the 4+ntype header lines and
write(unit) (vloc(ir),ir=1,ng(1)*ng(2)*ng(3))
where vloc is the sum of the hartree+exchange-correlation+
local pseudopotential on the real space grid in hartree
energy units. If spin-polarized, the xc potential in the
first record is for spin-up electrons and a second record
is written in which the xc potential is spin-down.
The underlying grid is as described above.
The header lines which accompany wf files, density files,
and potential files:
There are 4+ntype unformatted records which make up the header.
These are:
write(unit=header) codvsn,fform
write(unit=header) bantot,date,intxc,ixc,natom,(ng(mu),mu=1,3),
$ nkpt,nsppol,nsym,ntype,(acell(mu),mu=1,3),ecut,(rprimd(mu),mu=1,9)
write(unit=header) (nband(mu),mu=1,nkpt*nsppol),
$ (npwarr(mu),mu=1,nkpt),
$ (((symrel(mu,nu,isym),mu=1,3),nu=1,3),isym=1,nsym),
$ (type(mu),mu=1,natom),
$ (kpt(mu),mu=1,3*nkpt),
$ (occ(mu),mu=1,bantot),
$ ((tnons(mu,isym),mu=1,3),isym=1,nsym),
$ (zatnum(mu),mu=1,ntype)
(ntype lines, 1 for each psp...)
write(unit=unit) title,zatom,zion,pspdat,pspcod,pspxc,lmax,
$ lloc,mmax,r2well,e990,e999,nproj,rcpsp,logder,ekb1,ekb2,epspsp,
$ rchrg,fchrg,qchrg
(final record: residm and coordinates)
write(unit=unit) residm,(tn(i),i=1,3*natom)
codvsn,fform=integers
bantot,date,intxc,ixc,natom,ng,nkpt,nsppol,nsym,ntype=integers
acell,ecut,rprimd=double precision
nband,npwarr,symrel,type=integers
kpt,occ,tnons,zatnum=double precision
title=character string of 80 characters
zatom,zion=double precision
pspdat,pspcod,pspxc,lmax,lloc,mmax=integers
r2well,e990,e999=double precision
nproj=integer
rcpsp,logder,ekb1,ekb2,epspsp,rchrg,fchrg,qchrg=double precision
residm,tn=double precision
--------------------------------------------------------------------------
The following section describes the file structure used for
the pseudopotential files. These are presently formatted data files
which give potentials and projector functions on a real space
radial grid.
Firstly, the radial grid runs from index 0 to 2000 (2001 points),
with grid points given by the following equation (given as fortran):
nmax=2000
do j=0,nmax
x=dble(j)/dble(nmax)
r(j)=100.d0*(x+.01d0)**5-1.d-8
enddo
The psp file consists of a number of header lines followed by the
data on the radial grid. The header section is as follows:
title (single 80 character line)
zatom, zion, pspdat
pspcod, pspxc, lmax, lloc, mmax, r2well
...then, for l=0 to lmax, the following 2 lines:
l,e99.0,e99.9,nproj,rcpsp
rms,ekb1,ekb2,epsatm
...finally one more line:
rchrg,fchrg,qchrg
The data may be located anywhere on the line as long as it is provided
in the order indicated (it is read with free format).
In the case of Si with lmax=2, the header may look like the following
10 lines:
Si Fri Oct 08 11:18:59 1993
14.00000 4.00000 930920 zatom, zion, pspdat
1 1 2 2 2001 .00050 pspcod,pspxc,lmax,lloc,mmax,r2wel
l
0 19.464 25.000 2 1.8971118 l,e99.0,e99.9,nproj,rcpsp
.00112760 6.1457108933 4.4765165955 29.74712295 rms,ekb1,ekb2,epsatm
1 21.459 28.812 2 1.8971118 l,e99.0,e99.9,nproj,rcpsp
.00119946 3.2090654032 2.0935248528 19.11150542 rms,ekb1,ekb2,epsatm
2 8.223 21.459 0 1.8971118 l,e99.0,e99.9,nproj,rcpsp
.00098688 .0000000000 .0000000000 -3.97301006 rms,ekb1,ekb2,epsatm
1.70000000000000 .22513330685109 .96523597101781 rchrg,fchrg,qchr
g
zatom is the atomic number of the atom (14 for Si)
zion is the number of valence electrons (4 for Si)
pspdat is a code revision date (930920 for this case)
pspcod is another index describing the code (1 for this case)
pspxc is an index showing the choice of exchange-correlation (1)
lmax is the highest angular momentum for which a pseudopotential
is defined, which is also used for the local potential (2)
lloc is the angular momentum used for the local potential (2)
mmax is the number of grid points (2001)
r2well is the prefactor of a harmonic well sometimes used to bind
electrons which would otherwise be unbound in lda (.00050)
l is the angular momentum (0, 1, or 2 for Si for this case)
e99.0 is the planewave cutoff needed to converge to within 99.0% of the
kinetic energy of the atom (various numbers for various l)
e99.9 is the planewave cutoff needed to converge to within 99.9% of the
kinetic energy of the atom (various numbers for various l)
nproj is the number of projection functions used for each l (2)
rcpsp is the pseudopotential core radius
rms is a measure of pseudopotential quality reflecting the value of the
penalty function in designing the potential
ekb1, ekb2 are the Kleinman-Bylander energies for each projection
function for each l
epsatm is the integral Int[0 to Inf] (4*Pi*r*(r*V(r)+Zion))
rchrg is the core charge radius for additional core charge used to
match the xc contribution to the hardness matrix
fchrg is the prefactor of the core charge expression
qchrg is the total (integrated) core charge
Following the header are, for l=0 to lmax, the pseudopotentials
in the form of a title line followed by 667 lines of data, each line
containing three numbers so that the radial grid values from 0 to 2000
are given. The title line gives the value of l first followed by some
text. For the case of Si, e.g., for l=0, this line is
0 =l for Teter pseudopotential
...followed by the 667 lines, 3 numbers each, giving the l=0
potential on the radial grid described above.
Following the pseudopotentials are the first projection functions,
again given for each l with a title line followed by 667 data lines.
Following the first projection functions for each l are the second
projection functions, if any (determined by nproj).
Following the second projection functions, if any, are several lines
of additional data which is not read by plane_wave but is provided
to describe more about the details of the construction of the
pseudopotential. Omission of these lines does not affect the running
of plane_wave (at this time).
If you do not wish to use core charges, simply set fchrg to 0 and
use rchrg=1, qchrg=0.
If you wish to make a local potential, use lmax=lloc=0, nproj=0, and
you need not provide any projection function(s) (at this time).
The values of rms, ekb1, ekb2, epsatm, e99.0, e99.9 are used only
for information (at this time) so may be set to 0 when creating
a pseudopotential file. rcpsp is still used as the definition of the
pseudopotential core radius so it must be provided.
To best understand this data structure, it is recommended to study
several pseudopotential files and compare their contents with the
description given here.
--------------------------------------------------------------------------
The following section describes various parameters which affect
convergence and the numerical quality of calculations.
The list of input parameters is
acell amu dtion ecut ionmov kpt kptnrm ng nkpt nstep ntime rprim
vis wftol wtk
The list of dimension parameters within the main routine is
mqgrid
The technical design of the pseudopotential also affects the
quality of the results.
(1) The first issue regarding convergence is the number of planewaves
in the basis for a given set of atoms. Some atoms (notably those
in the first row or first transition series row) have relatively
deep pseudopotentials which require many planewaves for convergence.
In contrast are atoms like Si for which fewer planewaves are needed.
A typical value of "ecut" for Si might be 10 hartree for quite good
convergence, while the value for O might be 30-40 hartree or more
depending on the convergence desired and the design of the pseudo-
potential. NOTE: It is necessary in every new problem to TEST the
convergence by RAISING ecut for a given calculation until the
results being computed are constant to within some tolerance. This
is up to the user and is very important. For a given acell (apwsph)
and rprim, ecut is the parameter which controls the number of
planewaves. Of course if rprim or acell is varied then the number
of planewaves will also change.
Let me reiterate that extremely careful pseudopotential design can
optimize the convergence of e.g. the total energy within some range
of planewave number or ecut. It is appropriate to attempt to
optimize this convergence, especially for difficult atoms like
oxygen or copper, as long as one does not significantly compromise
the quality or transferability of the pseudopotential. There are
many people working on new techniques for optimizing convergence.
At Corning we have recently (August 1991) made progress designing
pseudopotentials with optimal convergence which are also optimal
in their transferability, so-called extended norm conserving
pseudopotentials.
For information on extended norm conservation, see
E. L. Shirley, D. C. Allan, R. M. Martin, and J. D. Joannopoulos,
Phys. Rev. B 40, 3652 (1989).
For information on optimizing the convergence of pseudopotentials,
see A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos,
Phys. Rev. B 41, 1227 (1990).
When the Corning implementation of improved pseudopotentials
is finalized, these potentials will be made available along with
the Corning electronic structure code.
(2) In addition to achieving convergence in the number of planewaves
in the basis, one must ensure that the iterations which solve the
electronic structure for a given set of atomic coordinates are also
converged. This convergence is controlled by the parameter wftol
and the parameter nstep. wftol gives a tolerance, in hartree^2,
for which the iterations will stop when the matrix element
<nk|(H-<nk|H|nk>)(H-<nk|H|nk>)|nk> for every band and k point n,k
is smaller than wftol. A typical value for adequate convergence is
10^-10; fairly good convergence might require 10^-12 or smaller.
The wavefunctions will be converged to roughly the square root of
wftol as defined above. The total energy will be converged, with
respect to preconditioned conjugate gradient iterations, to roughly
wftol itself. The nstep variable also controls convergence in
preconditioned conjugate gradient iterations by forcing the calculation
to stop whenever the number of such iterations exceeds nstep.
Usually one wants nstep to be set larger than needed to reach a given
wftol, or else one wants to restart insufficiently converged calculations
until wftol is reached.
Note that, if the gap in the system closes (e.g. due to defect formation
or if the system is metallic in the first place), the presently coded
algorithm may fail to converge. Convergence trouble during iterations
usually signals closure of the gap.
(3) For self consistent calculations (iscf=1) it is important to
test the adequacy of the k point integration. If symmetry is used
then one usually tests a set of "special point" grids. Otherwise
one tests the addition of more and more k points, presumably on
uniform grids, to ensure that a sufficient number has been included
for good k point integration. The parameter nkpt indicates how
many k points are being used, and their coordinates are given
by kpt and kptnrm, described above. The weight given to each k point
is provided by input variable wtk.
Systematic tests of k point integration are much more difficult than
tests of the adequacy of the number of planewaves. The difficulty I
refer to is simply the lack of a very systematic method for generating
k point grids for tests.
(4) It is possible to run calculations for which the fft box
is not quite large enough to avoid aliasing error in fft
convolutions. An aliasing error, or a fourier filter approximation,
is occurring when the output variable "boxcut" is less than 2.
boxcut is the smallest ratio of the fft box side to the planewave
basis sphere diameter. If this ratio is 2 or larger then e.g.
the calculation of the Hartree potential from the charge density
is done without approximation. At ratios smaller than 2, certain
of the highest fourier components are corrupted in the convolution.
If the basis is nearly complete, this fourier filter can be an
excellent approximation. In this case values of boxcut can be
as small as about 1.5 without incurring significant error. For
a given ecut, acell, and rprim, one should run tests for which
ng is large enough to give boxcut >= 2, and then one may try
smaller values of ng if the results are not significantly altered.
See the descriptions of these variables above.
(5) mqgrid is an internal array dimension (dimensioned in a
parameter statement in the main routine) which can affect the
quality of the calculation. It sets the grid spacing in the
variable "q" for the cubic spline interpolation used for reproducing
the local and nonlocal pseudopotentials Bessel function transforms
at all |G| or |k+G| values. The grid is set to have mqgrid
equally spaced values from q=0 to 1.2*gsqcut, where gsqcut is
(cutrad^2)*(2 * ecut)/(2 * Pi)^2 and cutrad is min(2, boxcut).
Thus the grid goes to a value 20% larger than the largest
(k+G)^2 expected. Until August 1991 I used mqgrid=301, but
recent tests indicate that mqgrid=1201 is needed in cases where
convergence studies to very large ecut are being conducted and
one wants up to 13 or 14 digits of the total energy to be reliably
computed. The effect of mqgrid on the calulation should be
checked by the user--it has not been systematically checked yet
by me. The source for the main routine is available to any user
and the value of mqgrid is adjustable by you for conducting your
own tests on your own problem.
(6) If you are running calculations to relax or equilibrate
structures, i.e. with ionmov=1 and possibly vis>0, then
the quality of your molecular dynamics will be affected by
the parameters amu, dtion, vis, ntime. Clearly if you want a
relaxed structure you must either run long enough or make repeated
runs until the largest force in the problem (output as fmax) is
smaller than what you will tolerate. If dtion is too large for
the given values of masses (amu) and viscosity (vis) then the
molecular dynamics will be unstable. If dtion is too small, then
the molecular dynamics will move inefficiently slowly.
A consensus exists in the community that forces larger than about
0.1 eV/Angstrom are really too large to consider the relaxation
to be converged. It is best for the user to get experience with
this in his/her own application.
More recently the option ionmov=2 has been made available
(3 June 1992). This uses the BFGS scheme for structural optimization
and is much more efficient than viscous damping for structural
relaxation.
-----------------------------------------------------------------------
In addition to the main routine for running these calculations,
there are several other main routines provided which perform
various useful functions. The source for these routines is provided
so that the user may modify his/her input or output to suit his/her needs,
or may alter array dimensions as needed. It is recommended that any
code modifications (other than array dimensions) be confined to the
insertion of user-written subroutine, so that when new versions of
these routines are provided the user may alter the new copy to add
the desired modifications with no inconvenience.
Note: only npwng and newsp are being carefully maintained.
The functionality of rhofwf is replaced by using corning with nstep=0.
(24 Jun 1993).
The additional main routines are the following:
bondlist
npwng
newsp
rhofwf
These routines have not been as extensively used by many users and
as extensively debugged as the main routine. In addition to the
short descriptions given below, there are additional comments in
the code listings.
bondlist Reads coordinate data either from an input file or from
a typical output file and makes a listing of bond lengths
and bond angles. This can be very helpful in debugging an
input file by checking bond connectivity.
npwng Computes the number of planewaves npw for a given set of
parameters and gives suggested fft grid dimensions ng for input
ecut, acell, and rprim. This is mainly useful to get
an idea of npw without running a full job, in order to
set dimensions correctly. The dimension mpw in various
main routines is the largest npw may be in a given run.
newsp This is an extremely useful program which allows conversion
of a wavefunction file from one ecut to another, one
set of lattice constants to another, or one set of k points
to another. It is used mainly to get new starting values
(for a new acell, ecut, or kpt) from existing converged
values. The code prompts the user to provide the name of
the original input file which produced the known wavefunction
file, then the name of a new input file which is intended
to be used in reading the new wavefunction file, then
the name of the existing unformatted wavefunction file,
and finally the name of the about-to-be-created new
wavefunction file. The second input file, for example,
would contain the new ecut, acell, and rprim. The
prompts provided at the terminal should make the use
of this subroutine self-explanatory. Note that if
a new k point grid is being created, the code will do
the crudest possible interpolation: it searches for the
nearest k point of the known wavefunction file and simply
copies that data to represent the new k point.
rhofwf Reads wavefunctions from standard unformatted wavefunction
file and computes the electron number density. This routine
again prompts the user for the names of files which contain
the needed data and again should be fairly self-explanatory.
The source rhofwf.f is available (in fortran) should the
user wish to modify the output format.
-----------------------------------------------------------------------
The Corning electronic structure code is written by Douglas C. Allan
and owned by Corning Incorporated of Corning, New York.
The ideas in the code have many sources; their implementation
in this code is mostly through the work of Mike Teter of Corning
Incorporated and Cornell University and Doug Allan of Corning
Incorporated, Cornell University, and The Ohio State University.
Capabilities related to response functions are developed by
Xavier Gonze of Louvain-la-Neuve Belgium.
Debugging has been aided by the efforts of many users who have tried
many unanticipated and sometimes downright remarkable things.
Please send questions and constructive criticisms of the code or
this documentation to
Douglas C. Allan
SP-FR-05
Corning Incorporated
Corning, NY 14831
telephone 607 974 3498
fax 607 974 3675
email allan_dc@corning.com
Correspondence by email is usually most convenient.