Bandstructure code from 1972

Electronic bandstructure of cubic YZn calculated with an APW code from the early 70's.

Are there still source codes of other old bandstructure programs available?


See further details about how to compile and use the code on nowadays hardware below.

In 1972, Hoffstein, Ray, and Belakhovsky have published a symmetrized augmented planewave code http://dx.doi.org/10.1016/0010-4655(72)90098-7 (Comput. Phys. Commun. 4 361 (1972)), which I've used to calculate the data shown in the figure above. Neither source nor binary must be distributed. You need access to the Computer Physics Communications Program Library to get the source code (tip: search for the authors to find the archive).

The program was written in Fortran IV (≈ Fortran 66) and the file you download is a punch code conversion containing the command line for the compilation, the source code, what is fed to standard input...
To extract the Fortran code and the input data run:

zcat acmj_v1_0.gz | tail +10 | head -1224 >sapw.f
zcat acmj_v1_0.gz | tail +1308 | head -313 >>sapw.f
zcat acmj_v1_0.gz | tail +1622 | head -506 >inph

We've skipped a few lines of the source code, which back in the days served for the calculation of a float encoded in a four character string in a special (but simple) way. The replacement of this routine for other hardware was actually worth a second publication by de Hosson in 1975 (http://dx.doi.org/10.1016/0010-4655(75)90091-0) ...

We'll take advantage of something else that was developed in 1972 - C - to implement this simple conversion:

/* conv.c */

double conv_(char *x)
{
    double sign = 1., c = 0., d;
    int i;

    for(i=0, d=1000.; i<4; i++, d/=10.) {
       switch(x[i]) {
       case '+':
           sign = 1.;
           break;
       case '-':
           sign = -1.;
           break;
       case 'A':
           c += sign * d * 0.86602540378;
           break;
       case 'B':
           c += sign * d * 0.5;
           break;
       case 'C':
           c += sign * d * 1.5;
       case 32:
           break;
       default:
           c += ((double) (x[i]-48)) * sign * d;
       }
    }

    return c;
}

Phew, what a masterpiece ;)

Compile with:

gcc -O2 -c conv.c
gfortran -O2 -fdefault-real-8 -o sapw sapw.f conv.o

The following python script can be used to calculate the dispersion from Γ to R:

from os import system,popen

e = []
for i in range(9):
  e.append({})

for s in (-2.0,-1.8,-1.5,-1.2,-1.0,-0.7,-0.5,-0.2,0.4,0.6,0.8):
  for i in range(9):
    k = float(i)/10.+0.1

    system('cat inph >inp')
    f = open('inp','a')
    print >>f,'       %g       %g       %g' % (k,k,k)
    print >>f,' -4  3  1  1  1 R 15 BANDS\n      %g     0.05   50' %s
    f.close()

    p = popen("./sapw <inp | grep EIGENVA | awk '{print $3}'")
    a = p.readline().strip()
    p.close()
    if len(a)>0:
      e[i][a] = 1

for i in range(9):
  k = float(i)/10.+0.1
  q = e[i].keys()
  q.sort()
  for a in q:
    print k, float(a)*2.0