==== usage ====
  * Copy wf2state to bin directory
  * copy particular cel.lvs into same directory as your *.xyz file
  * run with *.xyz file name as 1st parameter and shift by lattice vectors da, db, dc, example:
    <code> modxyz scan.xyz 0.5 0.5 0.0 </code>
  * script generate ''scan.xyz.s.xyz'' file

==== how it works ====
  * compute inverse matrix of lattice
  * transform geometry by the inverse lattice (make it rectangular)
  * shift by da,db,dz in this rectangularized space
  * modulo by cel size a,b,c 
  * transform back by forward lattice matrix

==== modxyz: ====
<code>
#! /bin/bash

#input="scan-.xyz"
input=$1
output=$input".s.xyz"

da=$2
db=$3
dc=$4

p=`cat $input | head -1`
let p=$p+2

awk -v p=$p -v da=$da -v db=$db -v dc=$dc '
function deter(a){
det= a[1,1]*(a[3,3]*a[2,2]-a[3,2]*a[2,3])-a[2,1]*(a[3,3]*a[1,2]-a[3,2]*a[1,3])+a[3,1]*(a[2,3]*a[1,2]-a[2,2]*a[1,3])
return det
}

function invert(a,b){
det = deter(a)

b[1,1] =  (a[3,3]*a[2,2]-a[3,2]*a[2,3]) / det
b[1,2] = -(a[3,3]*a[1,2]-a[3,2]*a[1,3]) /det
b[1,3] =  (a[2,3]*a[1,2]-a[2,2]*a[1,3]) /det

b[2,1] = -(a[3,3]*a[2,1]-a[3,1]*a[2,3]) /det
b[2,2] =  (a[3,3]*a[1,1]-a[3,1]*a[1,3]) /det
b[2,3] = -(a[2,3]*a[1,1]-a[2,1]*a[1,3]) /det

b[3,1] =  (a[3,2]*a[2,1]-a[3,1]*a[2,2]) /det
b[3,2] = -(a[3,2]*a[1,1]-a[3,1]*a[1,2]) /det
b[3,3] =  (a[2,2]*a[1,1]-a[2,1]*a[1,2]) /det
}

function Mvec(M,v,out){
out[1]=v[1]*M[1,1]+v[2]*M[1,2]+v[3]*M[1,3]
out[2]=v[1]*M[2,1]+v[2]*M[2,2]+v[3]*M[2,3]
out[3]=v[1]*M[3,1]+v[2]*M[3,2]+v[3]*M[3,3]
}

function Mvec2(M,v,out){
out[1]=v[1]*M[1,1]+v[2]*M[2,1]+v[3]*M[3,1]
out[2]=v[1]*M[1,2]+v[2]*M[2,2]+v[3]*M[3,2]
out[3]=v[1]*M[1,3]+v[2]*M[2,3]+v[3]*M[3,3]
}



BEGIN{
getline < "cel.lvs" 
A[1,1]=$1;A[1,2]=$2;A[1,3]=$3
getline < "cel.lvs"
A[2,1]=$1;A[2,2]=$2;A[2,3]=$3
getline < "cel.lvs"
A[3,1]=$1;A[3,2]=$2;A[3,3]=$3

invert(A,B)
}

#(NR==1){print $0}
#(NR==2){print $0}

(NR>0){ if ($4!=""){
v[1]=$2;v[2]=$3;v[3]=$4;
Mvec2(B,v,v_)

v_[1]=(v_[1]+da)%1.0
v_[2]=(v_[2]+db)%1.0
v_[3]=(v_[3]+dc)%1.0

#print $1,"   ",v_[1],"    ",v_[2],"   ",v_[3]
Mvec2(A,v_,v__)
print $1,"   ",v__[1],"    ",v__[2],"   ",v__[3]
}else{print $0}
}

END{
}

' $input > $output 
</code>