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         else
 231                 ((offset=8))
 232                 echo "tsd: no, offset = $offset"
 233         fi
 234         while [ "1" ]; do
 235
 236         echo -en "$count "
 237
 238                 if [ $count -le $nsi ]; then
 239                         color="Yellow"
 240                 else
 241                         color="Gray"
 242                 fi
 243                 export color
 244
 245                 ((gl=offset+count))
 246                 line="`sed -n ${gl}p $file`"
 247                 x=`echo $line | awk '{ print $1 }'`
 248                 y=`echo $line | awk '{ print $2 }'`
 249                 z=`echo $line | awk '{ print $3 }'`
 250
 251                 echo $x $y $z | awk '\
 252                 BEGIN {
 253                         x0=ENVIRON["x0"]; y0=ENVIRON["y0"]; z0=ENVIRON["z0"];
 254                         x1=ENVIRON["x1"]; y1=ENVIRON["y1"]; z1=ENVIRON["z1"];
 255                         X1=ENVIRON["X1"]; X2=ENVIRON["X2"]; X3=ENVIRON["X3"];
 256                         Y1=ENVIRON["Y1"]; Y2=ENVIRON["Y2"]; Y3=ENVIRON["Y3"];
 257                         Z1=ENVIRON["Z1"]; Z2=ENVIRON["Z2"]; Z3=ENVIRON["Z3"];
 258                         radius=ENVIRON["radius"]; scale=ENVIRON["scale"];
 259                         color=ENVIRON["color"]
 260                         x=0; y=0; z=0;
 261                         xt=0; yt=0; zt=0;
 262                         i=0; j=0; k=0;
 263                 }
 264                 {
 265                         for(i=-1;i<=1;i++) {
 266                                 for(j=-1;j<=1;j++) {
 267                                         for(k=-1;k<=1;k++) {
 268
 269                         nx=$1+i; ny=$2+j; nz=$3+k;
 270
 271                         xt=nx*X1+ny*Y1+nz*Z1;
 272                         yt=nx*X2+ny*Y2+nz*Z2;
 273                         zt=nx*X3+ny*Y3+nz*Z3;
 274
 275                         xt*=scale;
 276                         yt*=scale;
 277                         zt*=scale;
 278
 279                         if((xt>=x0)&&(yt>=y0)&&(zt>=z0)&&\
 280                            (xt<=x1)&&(yt<=y1)&&(zt<=z1)) {
 281                                 print "sphere { <"xt","zt","yt">, "radius" ";
 282                                 print "texture { pigment { color " color " } ";
 283                                 print "finish { phong 1 metallic } } }";
 284                         }
 285
 286                                         }
 287                                 }
 288                         }
 289                 }' >> temp.pov
 290
 291                 # see, whether there are bonds to draw
 292                 if [ "$ab" = "0" ]; then
 293                         ((count+=1))
 294                         continue
 295                 fi
 296                 ((cnt=1))
 297                 dobond=no
 298                 while [ $cnt -le $ab ]; do
 299                         anr=${anr[$cnt]}
 300                         if [ $anr -eq $count ]; then
 301                                 dobond=yes
 302                                 break
 303                         fi
 304                         ((cnt+=1))
 305                 done
 306                 if [ "$ab" = "-1" ]; then
 307                         dobond=yes
 308                 fi
 309                 if [ "$dobond" = "no" ]; then
 310                         ((count+=1))
 311                         continue;
 312                 fi
 313
 314                 # draw bonds
 315                 ((ic=1))
 316                 while [ "1" ]; do
 317
 318                         [ $ic -gt $count ] && break
 319
 320                         if [ $ic -eq $count ]; then
 321                                 ((ic+=1))
 322                                 continue
 323                         fi
 324
 325                         ((il=ic+offset))
 326                         line="`sed -n ${il}p $file`"
 327                         xb=`echo $line | awk '{ print $1 }'`
 328                         yb=`echo $line | awk '{ print $2 }'`
 329                         zb=`echo $line | awk '{ print $3 }'`
 330
 331                         echo $x $y $z $xb $yb $zb $cutoff | awk '\
 332                         BEGIN {
 333                         x0=ENVIRON["x0"]; y0=ENVIRON["y0"]; z0=ENVIRON["z0"];
 334                         x1=ENVIRON["x1"]; y1=ENVIRON["y1"]; z1=ENVIRON["z1"];
 335                         X1=ENVIRON["X1"]; X2=ENVIRON["X2"]; X3=ENVIRON["X3"];
 336                         Y1=ENVIRON["Y1"]; Y2=ENVIRON["Y2"]; Y3=ENVIRON["Y3"];
 337                         Z1=ENVIRON["Z1"]; Z2=ENVIRON["Z2"]; Z3=ENVIRON["Z3"];
 338                         bcr=ENVIRON["bcr"]; scale=ENVIRON["scale"];
 339                         }
 340                         {
 341                                 for(i=-1;i<=1;i++) {
 342                                         for(j=-1;j<=1;j++) {
 343                                                 for(k=-1;k<=1;k++) {
 344
 345                                 for(ii=-1;ii<=1;ii++) {
 346                                         for(ij=-1;ij<=1;ij++) {
 347                                                 for(ik=-1;ik<=1;ik++) {
 348
 349                                 nx=$1+i; ny=$2+j; nz=$3+k;
 350                                 xt=nx*X1+ny*Y1+nz*Z1;
 351                                 yt=nx*X2+ny*Y2+nz*Z2;
 352                                 zt=nx*X3+ny*Y3+nz*Z3;
 353                                 xt*=scale;
 354                                 yt*=scale;
 355                                 zt*=scale;
 356
 357                                 inx=$4+ii; iny=$5+ij; inz=$6+ik;
 358                                 ixt=inx*X1+iny*Y1+inz*Z1;
 359                                 iyt=inx*X2+iny*Y2+inz*Z2;
 360                                 izt=inx*X3+iny*Y3+inz*Z3;
 361                                 ixt*=scale;
 362                                 iyt*=scale;
 363                                 izt*=scale;
 364
 365                                 dist=sqrt((xt-ixt)^2+(yt-iyt)^2+(zt-izt)^2);
 366
 367                                 if((xt>=x0)&&(yt>=y0)&&(zt>=z0)&&\
 368                                    (xt<=x1)&&(yt<=y1)&&(zt<=z1)) {
 369                                         if(dist<=$7) {
 370                                                 print "cylinder {";
 371                                                 print "<"xt","zt","yt">,";
 372                                                 print "<"ixt","izt","iyt">,0.1";
 373                                                 print "pigment { color Blue }";
 374                                                 print "}";
 375                                         }
 376                                 }
 377
 378                                                 }
 379                                         }
 380                                 }
 381                                                 }
 382                                         }
 383                                 }
 384                         }' >> temp.pov
 385
 386                         ((ic+=1))
 387                         [ $ic -gt $ntot ] && break;
 388
 389                 done
 390
 391                 ((count+=1))
 392                 [ $count -gt $ntot ] && break;
 393
 394         done
 395
 396         echo
 397
 398         # manually drawing the 3x4 boundaries specified by argv ...
 399         draw_cyl $bx0 $by0 $bz0 $bx1 $by0 $bz0
 400         draw_cyl $bx0 $by0 $bz0 $bx0 $by1 $bz0
 401         draw_cyl $bx1 $by1 $bz0 $bx1 $by0 $bz0
 402         draw_cyl $bx0 $by1 $bz0 $bx1 $by1 $bz0
 403
 404         draw_cyl $bx0 $by0 $bz1 $bx1 $by0 $bz1
 405         draw_cyl $bx0 $by0 $bz1 $bx0 $by1 $bz1
 406         draw_cyl $bx1 $by1 $bz1 $bx1 $by0 $bz1
 407         draw_cyl $bx0 $by1 $bz1 $bx1 $by1 $bz1
 408
 409         draw_cyl $bx0 $by0 $bz1 $bx0 $by0 $bz0
 410         draw_cyl $bx0 $by1 $bz1 $bx0 $by1 $bz0
 411         draw_cyl $bx1 $by0 $bz1 $bx1 $by0 $bz0
 412         draw_cyl $bx1 $by1 $bz1 $bx1 $by1 $bz0
 413
 414         # add camera and light source
 415         cat >> temp.pov <<-EOF
 416 camera {
 417 EOF
 418         if [ -n "$ortographic" ]; then  cat >> temp.pov <<-EOF
 419 orthographic
 420 EOF
 421         fi
 422         cat >> temp.pov <<-EOF
 423 location $camloc
 424 look_at <$clx,$clz,$clz>
 425 }
 426 light_source { <0,100,-100> color White shadowless }
 427 EOF
 428
 429         # mv png
 430         $POVRAY temp.pov > /dev/null 2>&1
 431         mv temp.png `dirname $file`/
 432
 433 done
 434
 435 echo "done"