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.

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: wilkins@mps.ohio-state.edu [August1997]