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