vasp_tools/visualize_contcar

   1 #!/bin/sh
   2
   3 #
   4 # visualization script
   5 # author: frank.zirkelbach@physik.uni-augsburg.de
   6 #
   7
   8 # help function
   9 draw_cyl() {
  10         cat >> temp.pov <<-EOF
  11 cylinder {
  12 <$1, $3, $2>, <$4, $6, $5>, 0.05
  13 pigment { color White }
  14 }
  15 EOF
  16 }
  17 draw_bond() {
  18         cat >> temp.pov <<-EOF
  19 cylinder {
  20 <$1, $3, $2>, <$4, $6, $5>, $7
  21 pigment {
  22         color rgbf <0, 0, 1.0, 0.8>
  23         finish { phong 1 metallic }
  24 }
  25 }
  26 EOF
  27 }
  28
  29 # defaults
  30
  31 lc=5.429
  32 directory="doesnt_exist____for_sure"
  33 width="640"
  34 height="480"
  35 radius="0.6"
  36 x0="-0.6"; y0="-0.6"; z0="-0.6";
  37 x1="0.6"; y1="0.6"; z1="0.6";
  38 cx=""; cy=""; cz="";
  39 lx="0"; ly="-100"; lz="100";
  40 ortographic=""
  41 bx0=""; by0=""; bz0="";
  42 bx1=""; by1=""; bz1="";
  43 bcr="0.1";
  44 clx="0"; cly="0"; clz="0";
  45 extra=0
  46 displace=""
  47 mirror=0
  48 mx1=0; mx2=0; mx3=0;
  49 my1=0; my2=0; my3=0;
  50 mz1=0; mz2=0; mz3=0;
  51 ab=0
  52
  53 # parse argv
  54
  55 while [ "$1" ]; do
  56         case "$1" in
  57                 -d)             directory=$2;           shift 2;;
  58                 -w)             width=$2;               shift 2;;
  59                 -h)             height=$2;              shift 2;;
  60                 -r)             radius=$2;              shift 2;;
  61                 -nll)           x0=$2; y0=$3; z0=$4;    shift 4;;
  62                 -fur)           x1=$2; y1=$3; z1=$4;    shift 4;;
  63                 -c)             cx=$2; cy=$3; cz=$4;    shift 4;;
  64                 -L)             clx=$2; cly=$3; clz=$4; shift 4;;
  65                 -l)             lx=$2; ly=$3; lz=$4;    shift 4;;
  66                 -o)             ortographic=1;          shift 1;;
  67                 -b)             bx0=$2; by0=$3; bz0=$4;
  68                                 bx1=$5; by1=$6; bz1=$7; shift 7;;
  69                 -B)             bcr=$2;                 shift 2;;
  70                 -C)             lc=$2;                  shift 2;;
  71                 -e)             extra=1;                shift 1;;
  72                 -D)             displace=$2;            shift 2;;
  73                 -m)             mx1=$2; mx2=$3; mx3=$4;
  74                                 my1=$5; my2=$6; my3=$7;
  75                                 mz1=$8; mz2=$9; mz3=${10};
  76                                 mirror=1;               shift 10;;
  77                 -A)             ab=$2;                  shift 2;
  78                                 ((cnt=1))
  79                                 while [ $cnt -le $ab ]; do
  80                                         anr[$cnt]=$1;   shift 1;
  81                                         ((cnt+=1))
  82                                 done
  83                                 cutoff=$1;              shift 1;;
  84                 *)
  85                                 echo "options:"
  86                                 echo "########"
  87                                 echo "directory to progress:"
  88                                 echo "  -d <directory> (mandatory)"
  89                                 echo "png dim:"
  90                                 echo "  -w <width>"
  91                                 echo "  -h <height>"
  92                                 echo "atom size:"
  93                                 echo "  -r <radius>"
  94                                 echo "  -B <bond cylinder radius>"
  95                                 echo "unit cell:"
  96                                 echo "  -C <lattice constant>"
  97                                 echo "  -m <dimX> <dimY> <dimZ> (mirror atoms)"
  98                                 echo "visualization volume:"
  99                                 echo "  -nll <x> <y> <z> (near lower left)"
 100                                 echo "  -fur <x> <y> <z> (far upper right)"
 101                                 echo "  -o (ortographic)"
 102                                 echo "bounding box:"
 103                                 echo "  -b <x0> <y0> <z0> <x1> <y1> <z1>"
 104                                 echo "povray:"
 105                                 echo "  -c <x> <y> <z> (camera position)"
 106                                 echo "  -L <x> <y> <z> (camera look)"
 107                                 echo "  -l <x> <y> <z> (light source)"
 108                                 echo "bonds:"
 109                                 echo "  -ab <atom number> <cutoff> (auto bonds)"
 110                                 exit 1;;
 111         esac
 112 done
 113
 114 # calculation from lattic eunits to angstroms
 115
 116 [ "$lc" = "sic" ] && lc=4.359
 117 [ "$lc" = "si" ] && lc=5.480
 118 [ "$lc" = "c" ] && lc=3.566
 119
 120 #offset=`echo 0.125 \* $lc | bc`
 121 offset=0.0
 122
 123 x0=`echo $x0 \* $lc + $offset | bc`
 124 y0=`echo $y0 \* $lc + $offset | bc`
 125 z0=`echo $z0 \* $lc + $offset | bc`
 126 x1=`echo $x1 \* $lc + $offset | bc`
 127 y1=`echo $y1 \* $lc + $offset | bc`
 128 z1=`echo $z1 \* $lc + $offset | bc`
 129
 130 mx1=`echo $mx1 \* $lc + $offset | bc`
 131 my1=`echo $my1 \* $lc + $offset | bc`
 132 mz1=`echo $mz1 \* $lc + $offset | bc`
 133 mx2=`echo $mx2 \* $lc + $offset | bc`
 134 my2=`echo $my2 \* $lc + $offset | bc`
 135 mz2=`echo $mz2 \* $lc + $offset | bc`
 136 mx3=`echo $mx3 \* $lc + $offset | bc`
 137 my3=`echo $my3 \* $lc + $offset | bc`
 138 mz3=`echo $mz3 \* $lc + $offset | bc`
 139
 140 clx=`echo $clx \* $lc + $offset | bc`
 141 cly=`echo $cly \* $lc + $offset | bc`
 142 clz=`echo $clz \* $lc + $offset | bc`
 143
 144 if [ -n "$cx" -a -n "$cy" -a -n "$cz" ]; then
 145         cx=`echo $cx \* $lc + $offset | bc`
 146         cy=`echo $cy \* $lc + $offset | bc`
 147         cz=`echo $cz \* $lc + $offset | bc`
 148 fi
 149
 150 if [ -n "$bx0" ]; then
 151         bx0=`echo $bx0 \* $lc + $offset | bc`
 152         by0=`echo $by0 \* $lc + $offset | bc`
 153         bz0=`echo $bz0 \* $lc + $offset | bc`
 154         bx1=`echo $bx1 \* $lc + $offset | bc`
 155         by1=`echo $by1 \* $lc + $offset | bc`
 156         bz1=`echo $bz1 \* $lc + $offset | bc`
 157 fi
 158
 159 # povray command
 160 POVRAY="povray -W${width} -H${height} -d"
 161
 162 # camera location
 163 [ -n "$cx" -a -n "$cy" -a -n "$cz" ] && camloc="<$cx,$cz,$cy>"
 164
 165 # convert options
 166 COPTS="-font /usr/share/fonts/truetype/ttf-dejavu/DejaVuSans.ttf"
 167 COPTS="$COPTS -depth 8 -fill white -stroke blue -pointsize 24"
 168
 169 # do it ...
 170
 171 if [ -d $directory ]; then
 172         filesource=$directory/CONTCAR*
 173 fi
 174
 175 if [ -f $directory ]; then
 176         filesource=$directory
 177 fi
 178
 179 for file in $filesource; do
 180
 181         echo "working on $file ..."
 182
 183         cat > temp.pov <<-EOF
 184 #include "colors.inc"
 185 #include "textures.inc"
 186 #include "shapes.inc"
 187 #include "glass.inc"
 188 #include "metals.inc"
 189 #include "woods.inc"
 190 #include "stones.inc"
 191 EOF
 192
 193         # meta info
 194         scale=`sed -n 2p $file | awk '{ print $1 }'`
 195         line="`sed -n 3p $file`"
 196         X1=`echo $line | awk '{ print $1 }'`
 197         X2=`echo $line | awk '{ print $2 }'`
 198         X3=`echo $line | awk '{ print $3 }'`
 199         line="`sed -n 4p $file`"
 200         Y1=`echo $line | awk '{ print $1 }'`
 201         Y2=`echo $line | awk '{ print $2 }'`
 202         Y3=`echo $line | awk '{ print $3 }'`
 203         line="`sed -n 5p $file`"
 204         Z1=`echo $line | awk '{ print $1 }'`
 205         Z2=`echo $line | awk '{ print $2 }'`
 206         Z3=`echo $line | awk '{ print $3 }'`
 207         line="`sed -n 6p $file`"
 208         nsi=`echo $line | awk '{ print $1 }'`
 209         nc=`echo $line | awk '{ print $2 }'`
 210         ((ntot=nsi+nc))
 211
 212         echo "Si : $nsi"
 213         echo "C  : $nc"
 214         echo "tot: $ntot"
 215         echo "cutoff: $cutoff"
 216
 217         export scale
 218         export X1 X2 X3
 219         export Y1 Y2 Y3
 220         export Z1 Z2 Z3
 221
 222         export x0 y0 z0 x1 y1 z1 radius
 223
 224         # atoms
 225         ((count=1))
 226         checktsd=`sed -n 7p $file | awk '{ print $1 }'`
 227         if [ "$checktsd" = "Transformed" ]; then
 228                 ((offset=9))
 229                 echo "tsd: yes, offset = $offset"
 230                 # but maybe not
 231                 checkmore=`sed -n 8p $file | awk '{ print $1 }'`
 232                 if [ "$checkmore" = "Direct" ]; then
 233                         ((offset=8))
 234                         echo "tsd: yes, but anyways using offset $offset"
 235                 fi
 236         else
 237                 ((offset=8))
 238                 echo "tsd: no, offset = $offset"
 239         fi
 240         while [ "1" ]; do
 241
 242                 [ $count -gt $ntot ] && break
 243
 244                 if [ $count -le $nsi ]; then
 245                         color="Yellow"
 246                 else
 247                         color="Gray"
 248                 fi
 249                 export color
 250
 251                 ((gl=offset+count))
 252                 line="`sed -n ${gl}p $file`"
 253                 x=`echo $line | awk '{ print $1 }'`
 254                 y=`echo $line | awk '{ print $2 }'`
 255                 z=`echo $line | awk '{ print $3 }'`
 256
 257                 echo $x $y $z $count | awk '\
 258                 BEGIN {
 259                         x0=ENVIRON["x0"]; y0=ENVIRON["y0"]; z0=ENVIRON["z0"];
 260                         x1=ENVIRON["x1"]; y1=ENVIRON["y1"]; z1=ENVIRON["z1"];
 261                         X1=ENVIRON["X1"]; X2=ENVIRON["X2"]; X3=ENVIRON["X3"];
 262                         Y1=ENVIRON["Y1"]; Y2=ENVIRON["Y2"]; Y3=ENVIRON["Y3"];
 263                         Z1=ENVIRON["Z1"]; Z2=ENVIRON["Z2"]; Z3=ENVIRON["Z3"];
 264                         radius=ENVIRON["radius"]; scale=ENVIRON["scale"];
 265                         color=ENVIRON["color"]
 266                         x=0; y=0; z=0;
 267                         xt=0; yt=0; zt=0;
 268                         i=0; j=0; k=0;
 269                 }
 270                 {
 271                         for(i=-1;i<=1;i++) {
 272                                 for(j=-1;j<=1;j++) {
 273                                         for(k=-1;k<=1;k++) {
 274
 275                         nx=$1+i; ny=$2+j; nz=$3+k;
 276
 277                         xt=nx*X1+ny*Y1+nz*Z1;
 278                         yt=nx*X2+ny*Y2+nz*Z2;
 279                         zt=nx*X3+ny*Y3+nz*Z3;
 280
 281                         xt*=scale;
 282                         yt*=scale;
 283                         zt*=scale;
 284
 285                         if((xt>=x0)&&(yt>=y0)&&(zt>=z0)&&\
 286                            (xt<=x1)&&(yt<=y1)&&(zt<=z1)) {
 287                                 print "// atom " $4;
 288                                 print "sphere { <"xt","zt","yt">, "radius" ";
 289                                 print "texture { pigment { color " color " } ";
 290                                 print "finish { phong 1 metallic } } }";
 291                         }
 292
 293                                         }
 294                                 }
 295                         }
 296                 }' >> temp.pov
 297
 298                 # see, whether there are bonds to draw
 299                 if [ "$ab" = "0" ]; then
 300                         echo -en "$count "
 301                         ((count+=1))
 302                         continue
 303                 fi
 304                 ((cnt=1))
 305                 dobond=no
 306                 while [ $cnt -le $ab ]; do
 307                         anr=${anr[$cnt]}
 308                         if [ $anr -eq $count ]; then
 309                                 dobond=yes
 310                                 break
 311                         fi
 312                         ((cnt+=1))
 313                 done
 314                 if [ "$ab" = "-1" ]; then
 315                         if [ -n "`grep \/\/\ atom\ $count temp.pov`" ]; then
 316                                 dobond=yes
 317                         fi
 318                 fi
 319                 if [ "$dobond" = "no" ]; then
 320                         echo -en "$count "
 321                         ((count+=1))
 322                         continue;
 323                 fi
 324
 325                 # draw bonds
 326                 ((ic=1))
 327                 while [ "1" ]; do
 328
 329                         [ $ic -gt $count ] && break
 330
 331                         if [ $ic -eq $count ]; then
 332                                 ((ic+=1))
 333                                 continue
 334                         fi
 335
 336                         ((il=ic+offset))
 337                         line="`sed -n ${il}p $file`"
 338                         xb=`echo $line | awk '{ print $1 }'`
 339                         yb=`echo $line | awk '{ print $2 }'`
 340                         zb=`echo $line | awk '{ print $3 }'`
 341
 342                         echo $x $y $z $xb $yb $zb $cutoff | awk '\
 343                         BEGIN {
 344                         x0=ENVIRON["x0"]; y0=ENVIRON["y0"]; z0=ENVIRON["z0"];
 345                         x1=ENVIRON["x1"]; y1=ENVIRON["y1"]; z1=ENVIRON["z1"];
 346                         X1=ENVIRON["X1"]; X2=ENVIRON["X2"]; X3=ENVIRON["X3"];
 347                         Y1=ENVIRON["Y1"]; Y2=ENVIRON["Y2"]; Y3=ENVIRON["Y3"];
 348                         Z1=ENVIRON["Z1"]; Z2=ENVIRON["Z2"]; Z3=ENVIRON["Z3"];
 349                         bcr=ENVIRON["bcr"]; scale=ENVIRON["scale"];
 350                         }
 351                         {
 352                                 for(i=-1;i<=1;i++) {
 353                                         for(j=-1;j<=1;j++) {
 354                                                 for(k=-1;k<=1;k++) {
 355
 356                                 for(ii=-1;ii<=1;ii++) {
 357                                         for(ij=-1;ij<=1;ij++) {
 358                                                 for(ik=-1;ik<=1;ik++) {
 359
 360                                 nx=$1+i; ny=$2+j; nz=$3+k;
 361                                 xt=nx*X1+ny*Y1+nz*Z1;
 362                                 yt=nx*X2+ny*Y2+nz*Z2;
 363                                 zt=nx*X3+ny*Y3+nz*Z3;
 364                                 xt*=scale;
 365                                 yt*=scale;
 366                                 zt*=scale;
 367
 368                                 if((xt>=x0)&&(yt>=y0)&&(zt>=z0)&&\
 369                                    (xt<=x1)&&(yt<=y1)&&(zt<=z1)) {
 370
 371                                 inx=$4+ii; iny=$5+ij; inz=$6+ik;
 372                                 ixt=inx*X1+iny*Y1+inz*Z1;
 373                                 iyt=inx*X2+iny*Y2+inz*Z2;
 374                                 izt=inx*X3+iny*Y3+inz*Z3;
 375                                 ixt*=scale;
 376                                 iyt*=scale;
 377                                 izt*=scale;
 378
 379                                 dist=sqrt((xt-ixt)^2+(yt-iyt)^2+(zt-izt)^2);
 380
 381                                 #if((xt>=x0)&&(yt>=y0)&&(zt>=z0)&&\
 382                                 #   (xt<=x1)&&(yt<=y1)&&(zt<=z1)) {
 383                                         if(dist<=$7) {
 384                                                 print "cylinder {";
 385                                                 print "<"xt","zt","yt">,";
 386                                                 print "<"ixt","izt","iyt">,0.1";
 387                                                 print "pigment { color Blue }";
 388                                                 print "}";
 389                                         }
 390                                 #}
 391
 392                                 }
 393
 394                                                 }
 395                                         }
 396                                 }
 397                                                 }
 398                                         }
 399                                 }
 400                         }' >> temp.pov
 401
 402                         ((ic+=1))
 403                         [ $ic -gt $ntot ] && break;
 404
 405                 done
 406
 407                 echo -en "[$count] "
 408                 ((count+=1))
 409                 [ $count -gt $ntot ] && break;
 410
 411         done
 412
 413         echo
 414
 415         # manually drawing the 3x4 boundaries specified by argv ...
 416         if [ -n "$bx0" ]; then
 417         draw_cyl $bx0 $by0 $bz0 $bx1 $by0 $bz0
 418         draw_cyl $bx0 $by0 $bz0 $bx0 $by1 $bz0
 419         draw_cyl $bx1 $by1 $bz0 $bx1 $by0 $bz0
 420         draw_cyl $bx0 $by1 $bz0 $bx1 $by1 $bz0
 421
 422         draw_cyl $bx0 $by0 $bz1 $bx1 $by0 $bz1
 423         draw_cyl $bx0 $by0 $bz1 $bx0 $by1 $bz1
 424         draw_cyl $bx1 $by1 $bz1 $bx1 $by0 $bz1
 425         draw_cyl $bx0 $by1 $bz1 $bx1 $by1 $bz1
 426
 427         draw_cyl $bx0 $by0 $bz1 $bx0 $by0 $bz0
 428         draw_cyl $bx0 $by1 $bz1 $bx0 $by1 $bz0
 429         draw_cyl $bx1 $by0 $bz1 $bx1 $by0 $bz0
 430         draw_cyl $bx1 $by1 $bz1 $bx1 $by1 $bz0
 431         fi
 432
 433         # add camera and light source
 434         cat >> temp.pov <<-EOF
 435 camera {
 436 EOF
 437         if [ -n "$ortographic" ]; then  cat >> temp.pov <<-EOF
 438 orthographic
 439 EOF
 440         fi
 441         cat >> temp.pov <<-EOF
 442 location $camloc
 443 look_at <$clx,$clz,$clz>
 444 }
 445 light_source { <0,100,-100> color White shadowless }
 446 EOF
 447
 448         # mv png
 449         $POVRAY temp.pov > /dev/null 2>&1
 450         mv temp.png `dirname $file`/
 451
 452 done
 453
 454 echo "done"