#!/bin/csh -f
# code to multiply input POSCAR or CONTCAR file in three dimensions

if (! -e POTCAR || $#argv < 4) then
  set script = $0
  set script = $script:t
  if (! -e POTCAR) echo "\nprovide a POTCAR file ########\n"
err:
  cat <<koniec

usage: $script nx ny nz zscale {nodisp or xbsa or rasmol} {rscale} {zvac} < {POSCAR or CONTCAR}  > {new POSCAR}

- nx, ny, nz determines how often the cell shall be repeated in each direction
- zscale scales the z dimension of the supercell
- nodisp or xbsa or rasmol, if set, determines which kind of display to use
- rscale, if set, determines the geometry scale factor (useful to
    control bonds displayed in rasmol or xbsa. Only bond lengths below the sum of
    covalent radii are displayed.)
- {POSCAR or CONTCAR}  is the input file name in POSCAR format
- {new POSCAR}  is the name of a POSCAR type file that will be crated and can be
    used in a vasp run. This is a convenient way to construct a big supercell from
    a small one. Using zscale > nz allows for adding vacuum. The new
    POSCAR file is numbered and the species of each atom identified.
- zvac, if set, determines a z coordinate inside the vacuum such that slabs will appear continous
- an appropriate POTCAR file with the number of species matching that in the input
    file has to exist in the current directory

info on 3D display:
- the atomic positions will be displayed in 3D using rasmol or xbsa. The rasmol file
   is vasp.pdb, the xyz or xmol file created is vasp.xyz and vasp.bs
- in rasmol, determine distance by clicking on atoms. Change to angles by "set picking angle".
    rasmol manual = www.umass.edu/microbio/rasmol/distrib/rasman.htm
-  xbsa:   cursor to rotate, N to display distances, + to zoom, h for help

Please report problems to Roland

koniec
  exit 1
endif

set rscale = 1.0
if ($#argv == 6) set rscale = $6
set zvac = 1.0
if ($#argv == 7) set zvac = $7

cat >awk1.tmp<<koniec
BEGIN{ispect=0;ispecr=0; print "BEGIN{"}
/TITEL/{++ispect; 
  name=\$4; name = substr(name,1,2);
  if (substr(name,2,1) == "_") {name = substr(name,1,1)};
  print "specs["ispect"] = \""name"\";"
}
/RWIGS/{++ispecr; print "rads["ispecr"] = \""\$6"\";"}
koniec
awk -f awk1.tmp <POTCAR>awk.tmp

cat <<koniec >>awk.tmp
  OFMT = "%14.8f"; sfound=0; dfound=0; lr=15;
  natomi=0; natomo=0; nspec=0; ispec=0; ni[0]=0;
  format =  "ATOM  %.5d  %-2.2s          1    %8.3f%8.3f%8.3f%6.2f%6.2f\n";
#           ATOM      1  N   ASP L   1       4.060   7.307   5.186  1.00 51.58   
#   format(a6,i5,1x,a2,a3,a3,1x,a1,i4,a1,3x,3f8.3,2f6.2,1x,i3)
  formatc = "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1\n";
  formata = "%14.8f %14.8f %14.8f %s %s %s   %s %i\n";
# REMARK 1234567890123456789012345678901234567890123456789012345678901234567890
#        CRYST1   11.732   13.000  19.799   90.00  90.00  90.00 P 1        0
}
{x = substr(\$1, 1, lr); y = substr(\$2, 1, lr); z = substr(\$3, 1, lr)} 
{if (NR==1) {comm = \$0; print "scale:", sx, sy, sz, szz, int(sz), comm}
  else if (NR==2) {alat = \$1 * $rscale; print}
  else if (NR==3) {x1 = x * sx;  x2 = y * sx; x3 = z * sx; 
    print x1, x2, x3, \$4, \$5, \$6}
  else if (NR==4) {y1 = x * sy;  y2 = y * sy; y3 = z * sy; 
    print y1, y2, y3, \$4, \$5, \$6}
  else if (NR==5) {z1 = x * szz;  z2 = y * szz; z3 = z * szz; 
    print z1, z2, z3, \$4, \$5, \$6;
    if (x2 == 0 && x3 == 0 && y1 == 0 && y3 == 0 && z1 == 0 && z2 == 0) {
      printf formatc, x1*alat/sx, y2*alat/sy, z3*alat/sz, 90, 90, 90 >"vasp.pdb"}}
  else if (NR==6) {out = ""; ntot = 0;
    for (i=0; ++i <= NF && \$i ~ "[1-9]";)
      {ni[i] = \$i; no[i] = \$i*int(sx+.5)*int(sy+.5)*int(sz+.5);
       ntot += no[i];  out = out " " no[i]; nspec++};
    print ntot >"vasp.xyz";
    print comm >"vasp.xyz";
    print out}
  else if (NR==7) {print;
    if (\$1 ~ "[sS]") {sfound=1; NR--}
    else if (\$1 ~ "[dD]") {dfound=1}
    else if (\$1 ~ "[cC]") {dfound=2}}
  else if (NR > 7) {
    ++natomi;
    if (natomi > ni[ispec]) {++ispec;
      if (ispec > nspec) exit;
      natomi=1;  spec=specs[ispec]; rad=rads[ispec]
    };
    for (iz=0; iz < int(sz+.5); ++iz) {
      for (ix=0; ix < int(sx+.5); ++ix) {
        for (iy=0; iy < int(sy+.5); ++iy) {
          natom++;
          if (dfound == 1) {
            xx = (x+ix)/int(sx+.5); yy = (y+iy)/int(sy+.5); 
            if (z > $zvac) {z -= 1}; 
            if ((sz-szz)^2 > .00000001) {zz = (z+iz)/szz}
            else {zz = (z+iz)/int(sz+.5)};
            xa = alat * (xx * x1 + yy * y1 + zz * z1);
            ya = alat * (xx * x2 + yy * y2 + zz * z2);
            za = alat * (xx * x3 + yy * y3 + zz * z3)}
          else {
            xx = x + ix*x1/sx + iy*y1/sy + iz*z1/sz;
            yy = y + ix*x2/sx + iy*y2/sy + iz*z2/sz;
            zz = z + ix*x3/sx + iy*y3/sy + iz*z3/sz;
            xa = alat * xx;
            ya = alat * yy;
            za = alat * zz};
          if (xa == 0) xa = 0.00000001;
          if (ya == 0) ya = 0.00000001;
          if (za == 0) za = 0.00000001;
          printf formata, xx, yy, zz, \$4, \$5, \$6, spec, natom;
          printf format, natom, spec, xa, ya, za, 1.0, rad >"vasp.pdb"
          print spec, xa, ya, za >"vasp.xyz"
    }}}}
  else {print}
}
koniec
if (-e vasp.pdb) rm -f vasp.pdb
if (-e vasp.xyz) rm -f vasp.xyz
awk -f awk.tmp sx=$1 sy=$2 sz=$3 szz=$4

if (! -e ~/.rasmolrc) then
cat>~/.rasmolrc<<koniec
spacefill temperature
color cpk
set specular on
set background gray
set picking distance
zoom 90
set unitcell on
connect true
koniec
endif

#xyz2bs vasp.xyz
#if ("$5" == "xbsa") then
#  xbsa vasp.bs >& /dev/null &
##  echo "Use h to get help in the xbsa window"
#else if ("$5" != "nodisp") then
##   (xterm -fn screen15 +ls -e rasmol -script ~/.rasmolrc vasp.pdb &) >& /dev/null
#  rasmol -script ~/.rasmolrc vasp.pdb >& /dev/null &
#endif