vasp_tools/angle_calc

   1 #!/bin/bash
   2
   3 file=$1
   4 atom=$2
   5 btom=$3
   6 ctom=$4
   7
   8 ((offset=8))
   9
  10 scale=`sed -n 2p $file`
  11 X1=`sed -n 3p $file | awk '{ print $1 }'`
  12 X2=`sed -n 3p $file | awk '{ print $2 }'`
  13 X3=`sed -n 3p $file | awk '{ print $3 }'`
  14 Y1=`sed -n 4p $file | awk '{ print $1 }'`
  15 Y2=`sed -n 4p $file | awk '{ print $2 }'`
  16 Y3=`sed -n 4p $file | awk '{ print $3 }'`
  17 Z1=`sed -n 5p $file | awk '{ print $1 }'`
  18 Z2=`sed -n 5p $file | awk '{ print $2 }'`
  19 Z3=`sed -n 5p $file | awk '{ print $3 }'`
  20 export X1 X2 X3
  21 export Y1 Y2 Y3
  22 export Z1 Z2 Z3
  23
  24 ((line1=atom+offset))
  25 ((line2=btom+offset))
  26 ((line3=ctom+offset))
  27
  28 temp="`sed -n ${line1}p $1`"
  29 xa=`echo $temp | awk '{ print $1 }'`
  30 ya=`echo $temp | awk '{ print $2 }'`
  31 za=`echo $temp | awk '{ print $3 }'`
  32
  33 temp="`sed -n ${line2}p $1`"
  34 xb=`echo $temp | awk '{ print $1 }'`
  35 yb=`echo $temp | awk '{ print $2 }'`
  36 zb=`echo $temp | awk '{ print $3 }'`
  37
  38 temp="`sed -n ${line3}p $1`"
  39 xc=`echo $temp | awk '{ print $1 }'`
  40 yc=`echo $temp | awk '{ print $2 }'`
  41 zc=`echo $temp | awk '{ print $3 }'`
  42
  43 echo -en "angle: "
  44 foo=`echo "$xa $ya $za $xb $yb $zb $xc $yc $zc $scale" | \
  45         awk ' \
  46         BEGIN {
  47                 X1=ENVIRON["X1"]; X2=ENVIRON["X2"]; X3=ENVIRON["X3"]
  48                 Y1=ENVIRON["Y1"]; Y2=ENVIRON["Y2"]; Y3=ENVIRON["Y3"]
  49                 Z1=ENVIRON["Z1"]; Z2=ENVIRON["Z2"]; Z3=ENVIRON["Z3"]
  50                 pi=3.14159265
  51         }
  52         {
  53                 X=sqrt(X1^2+X2^2+X3^2)
  54                 Y=sqrt(Y1^2+Y2^2+Y3^2)
  55                 Z=sqrt(Z1^2+Z2^2+Z3^2)
  56                 dx=$1-$4
  57                 dy=$2-$5
  58                 dz=$3-$6
  59                 if(dx>1/2)
  60                         dx-=1
  61                 if(dx<-1/2)
  62                         dx+=1
  63                 if(dy>1/2)
  64                         dy-=1
  65                 if(dy<-1/2)
  66                         dy+=1
  67                 if(dz>1/2)
  68                         dz-=1
  69                 if(dz<-1/2)
  70                         dz+=1
  71                 Dx=$1-$7
  72                 Dy=$2-$8
  73                 Dz=$3-$9
  74                 if(Dx>1/2)
  75                         Dx-=1
  76                 if(Dx<-1/2)
  77                         Dx+=1
  78                 if(Dy>1/2)
  79                         Dy-=1
  80                 if(Dy<-1/2)
  81                         Dy+=1
  82                 if(Dz>1/2)
  83                         Dz-=1
  84                 if(Dz<-1/2)
  85                         Dz+=1
  86                 dxt=dx*X1+dy*Y1+dz*Z1
  87                 dyt=dx*X2+dy*Y2+dz*Z2
  88                 dzt=dx*X3+dy*Y3+dz*Z3
  89                 Dxt=Dx*X1+Dy*Y1+Dz*Z1
  90                 Dyt=Dx*X2+Dy*Y2+Dz*Z2
  91                 Dzt=Dx*X3+Dy*Y3+Dz*Z3
  92                 sp=dxt*Dxt+dyt*Dyt+dzt*Dzt
  93                 d=sqrt(dxt^2+dyt^2+dzt^2)
  94                 D=sqrt(Dxt^2+Dyt^2+Dzt^2)
  95                 print sp/(d*D)
  96         }'`
  97
  98 ./acos $foo
  99
 100 echo
 101