Help for Corning LDA Code

This help file contains 1390 lines!


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:
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

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.
	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.
	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
	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

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

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


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
	(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

    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

(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)

acell,ecut,rprimd=double precision
kpt,occ,tnons,zatnum=double precision
title=character string of 80 characters
zatom,zion=double precision
r2well,e990,e999=double precision
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):
do j=0,nmax

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:
...finally one more line:

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
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

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
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

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


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:

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
 Corning Incorporated
 Corning, NY  14831
 telephone 607 974 3498
 fax       607 974 3675

Correspondence by email is usually most convenient.

Your comments and suggestions are appreciated.
[Previous] [Wilkins Home Page]
[OSU Physics] [OSU Chemical Physics] [Center for Materials Research]
[College of Mathematical & Physical Sciences] [Ohio State University]
Edited by: [August1997]