#!/bin/ksh
#
# Constructs the SKIN file to be used in minimizing the total
# energy of the strained C2/m lattice for finding C44 of an
# hcp crystal.  See
# http://web.archive.org/web/20111231202250/http://cst-www.nrl.navy.mil/bind/static/example10/c44/c44r.html
# for details.
#
#                                         -- mjm  23 June 1998
#
# This software and any accompanying documentation are released "as
# is." The U.S. Government makes no warranty of any kind, expressed
# or implied, concerning this software and any accompanying
# documentation, including without limitation, any warranties of
# merchantability or fitness for a particular purpose. In no event
# will the U.S. Government be liable for any damages, including any
# lost profits, lost savings, or other incidental or consequential
# damages arising out of the use, or inability of use, of this
# software or any accompanying documentation, even if informed in
# advance of the possibility of such damages.
#
# Input is
#
# Argument 1:  Volume
# Argument 2:  c/a ratio
# Argument 3:  square of the strain parameter
# Argument 4:  search radius h
# Argument 5:  initial guess for u parameter
# Argument 6:  initial guess for v parameter

# Location of static code (edit for your machine)

# static=/usr/local/mehl/skfits/p6/static.106/static
static=/home/dumpy/mehl/skfits/many/static.106/static

# Check Argument list:

if (( $# < 6 ))
then
	echo Usage:  latcon vol covera strainsq h u0 v0
	exit
fi

# Construct the points needed for the calculation.  I'm using awk.
# A perl expert could certainly do better.

vol=$1
covera=$2
strainsq=$3

# Make sure we have the appropriate accumulation file:

touch skeng.$strainsq

h=`echo $4 | awk '{printf "%18.15f\n", $1}'`

u1=`echo $5 | awk '{printf "%17.15f\n", $1}'`
v1=`echo $6 | awk '{printf "%17.15f\n", $1}'`

u2=`echo $u1 $h | awk '{printf "%17.15f\n", $1+$2}'`
v2=$v1

u3=`echo $u1 $h | awk '{printf "%17.15f\n", $1+0.309016994374945*$2}'`
v3=`echo $v1 $h | awk '{printf "%17.15f\n", $1+0.951056516295154*$2}'`

u4=`echo $u1 $h | awk '{printf "%17.15f\n", $1-0.809016994374950*$2}'`
v4=`echo $v1 $h | awk '{printf "%17.15f\n", $1+0.587785252292469*$2}'`

u5=$u4
v5=`echo $v1 $h | awk '{printf "%17.15f\n", $1-0.587785252292469*$2}'`

u6=$u3
v6=`echo $v1 $h | awk '{printf "%17.15f\n", $1-0.951056516295154*$2}'`

# Labels have to fit in 20 characters:

label1=`echo 1 $u1 $v1 | awk '{printf "%2g%9.6f%9.6f\n", $1, $2, $3}'`
label2=`echo 2 $u2 $v2 | awk '{printf "%2g%9.6f%9.6f\n", $1, $2, $3}'`
label3=`echo 3 $u3 $v3 | awk '{printf "%2g%9.6f%9.6f\n", $1, $2, $3}'`
label4=`echo 4 $u4 $v4 | awk '{printf "%2g%9.6f%9.6f\n", $1, $2, $3}'`
label5=`echo 5 $u5 $v5 | awk '{printf "%2g%9.6f%9.6f\n", $1, $2, $3}'`
label6=`echo 6 $u6 $v6 | awk '{printf "%2g%9.6f%9.6f\n", $1, $2, $3}'`

echo $label1 $label2 $label3 $label4 $label5 $label6


# Now construct the SKIN file:

cat <<EOF > SKIN
Mode=3
0.00500   5.000000000000000
../ti_par
Titanium HCP (A3) -- Monoclinic Strain
$label1
0.00000
-3
$vol   $covera
6
$strainsq
2
4 4 4
F
1   $u1  -$u1   $v1 0 0 0
1  -$u1   $u1  -$v1 0 0 0
NEWSYM=T
LATTIC=4
spcgp1.12
ILAT=T
-1313
kmesh.12
Titanium HCP (A3) -- Monoclinic Strain
$label2
0.00000
-3
$vol   $covera
6
$strainsq
2
4 4 4
F
1   $u2  -$u2   $v2 0 0 0
1  -$u2   $u2  -$v2 0 0 0
NEWSYM=F
Titanium HCP (A3) -- Monoclinic Strain
$label3
0.00000
-3
$vol   $covera
6
$strainsq
2
4 4 4
F
1   $u3  -$u3   $v3 0 0 0
1  -$u3   $u3  -$v3 0 0 0
NEWSYM=F
Titanium HCP (A3) -- Monoclinic Strain
$label4
0.00000
-3
$vol   $covera
6
$strainsq
2
4 4 4
F
1   $u4  -$u4   $v4 0 0 0
1  -$u4   $u4  -$v4 0 0 0
NEWSYM=F
Titanium HCP (A3) -- Monoclinic Strain
$label5
0.00000
-3
$vol   $covera
6
$strainsq
2
4 4 4
F
1   $u5  -$u5   $v5 0 0 0
1  -$u5   $u5  -$v5 0 0 0
NEWSYM=F
Titanium HCP (A3) -- Monoclinic Strain
$label6
0.00000
-3
$vol   $covera
6
$strainsq
2
4 4 4
F
1   $u6  -$u6   $v6 0 0 0
1  -$u6   $u6  -$v6 0 0 0
NEWSYM=F
EOF

# Now run the code

$static > /dev/null

cat SKENG

# Save for further use:

cat SKENG >> skeng.$strainsq

# Pick off the energies:

eng1=`grep '^ 1' SKENG | awk '{print $6}'`
eng2=`grep '^ 2' SKENG | awk '{print $6}'`
eng3=`grep '^ 3' SKENG | awk '{print $6}'`
eng4=`grep '^ 4' SKENG | awk '{print $6}'`
eng5=`grep '^ 5' SKENG | awk '{print $6}'`
eng6=`grep '^ 6' SKENG | awk '{print $6}'`

# echo $eng1
# echo $eng2
# echo $eng3
# echo $eng4
# echo $eng5
# echo $eng6

# Auxillary quantities needed to fix the fit:

# Constants are 1+1/sqrt(5) and 1-1/sqrt(5)

q1=`echo $eng3 $eng6 $eng4 $eng5 $eng1 | awk '{printf "%25.15g\n", 1.447213595499958*($1+$2)+0.552786404500042*($3+$4) - 4*$5}'`

q2=`echo $eng1 $eng2 $eng3 $eng4 $eng5 $eng6 | awk '{printf "%25.15g\n", 2*(.2*($2+$3+$4+$5+$6)-$1)}'`

q3=`echo $eng2 $eng1 | awk '{printf "%25.15g\n", $1-$2}'`

# Constant is 2/sqrt[10+2sqrt(5)]
q4=`echo $eng3 $eng6 | awk '{printf "%25.15g\n", 0.525731112119134*($1-$2)}'`

# Constant is 2/sqrt[10-2sqrt(5)]
q5=`echo $eng4 $eng5 | awk '{printf "%25.15g\n", 0.850650808352041*($1-$2)}'`

# echo "q1 = " $q1
# echo "q2 = " $q2
# echo "q3 = " $q3
# echo "q4 = " $q4
# echo "q5 = " $q5

# Now for the actual fitting constants:

a00=$eng1

# Constant is 1/(sqrt((3-sqrt(5))/8) + (sqrt(5)+1)/4)

a11=`echo $q4 $q5 | awk '{printf "%25.15g\n", 0.894427190999917*($1-$2)}'`

a02=`echo $q1 $q2 | awk '{printf "%25.15g\n", 0.5*($1-$2)}'`

a20=`echo $q2 $a02 | awk '{printf "%25.15g\n", $1-$2}'`

a10=`echo $q3 $a20 | awk '{printf "%25.15g\n", $1-$2}'`

# Constant is (sqrt(5)+1)/4

a01=`echo $q5 $a11 | awk '{printf "%25.15g\n", $1+0.809016994374947*$2}'`

# echo $a00
# echo $a10
# echo $a01
# echo $a20
# echo $a02
# echo $a11

# Note that the fit is in terms of scaled constants,
# ($ui-$u1)/$h and ($vi-$v1)/$h 

# Now check fit if so desired:

# echo `echo $a00 $a10 $a01 $a20 $a02 $a11 0 0 | awk '{printf "%15.9f\n", $1 + $7*($2+$7*$4) + $8*($3+$8*$5) + $6*$7*$8}'` $eng1

# echo `echo $a00 $a10 $a01 $a20 $a02 $a11 1 0 | awk '{printf "%15.9f\n", $1 + $7*($2+$7*$4) + $8*($3+$8*$5) + $6*$7*$8}'` $eng2

# echo `echo $a00 $a10 $a01 $a20 $a02 $a11 0.309016994374945 0.951056516295154 | awk '{printf "%15.9f\n", $1 + $7*($2+$7*$4) + $8*($3+$8*$5) + $6*$7*$8}'` $eng3

# echo `echo $a00 $a10 $a01 $a20 $a02 $a11 -0.809016994374950 0.587785252292469 | awk '{printf "%15.9f\n", $1 + $7*($2+$7*$4) + $8*($3+$8*$5) + $6*$7*$8}'` $eng4

# echo `echo $a00 $a10 $a01 $a20 $a02 $a11 -0.809016994374950 -0.587785252292469 | awk '{printf "%15.9f\n", $1 + $7*($2+$7*$4) + $8*($3+$8*$5) + $6*$7*$8}'` $eng5

# echo `echo $a00 $a10 $a01 $a20 $a02 $a11 0.309016994374945 -0.951056516295154 | awk '{printf "%15.9f\n", $1 + $7*($2+$7*$4) + $8*($3+$8*$5) + $6*$7*$8}'` $eng6

# Now find the extremum of this fit, in the reduced coordinate set.
# Note that if we do not have a good start we may find a maximum
#  rather than a minimum.

xmin=`echo $a10 $a01 $a20 $a02 $a11 | awk '{printf "%20.15f\n", ($2*$5-2*$1*$4)/(4*$3*$4-$5*$5)}'`

ymin=`echo $a10 $a01 $a20 $a02 $a11 | awk '{printf "%20.15f\n", ($1*$5-2*$2*$3)/(4*$3*$4-$5*$5)}'`

engf=`echo $a00 $a10 $a01 $a20 $a02 $a11 $xmin $ymin | awk '{printf "%15.9f\n", $1 + $7*($2+$7*$4) + $8*($3+$8*$5) + $6*$7*$8}'`

# echo $xmin $ymin $engf

# Construct "actual" point:

umin=`echo $u1 $h $xmin | awk '{printf "%20.15f\n", $1+$2*$3}'`
vmin=`echo $v1 $h $ymin | awk '{printf "%20.15f\n", $1+$2*$3}'`

# And finally, print the output in a form that can be used for the
#  next iteration:

echo $vol $covera $strainsq $h $umin $vmin $engf

# Clean up

rm SKIN SKENG SKOUT