#!/usr/bin/perl

#
# This program calculates boxvectors compliant with Gromacs.
# Written by Tsjerk A. Wassenaar, University of Groningen (2002)
#

$usage_message="
NAME
        Perl script boxmtx.pl TAW021025
 
SYNOPSIS
        boxmtx.pl x1 y1 z1 x2 y2 z2 x3 y3 z3

        This script is written to perform a box transformation,
        through which a given triclinic box is made compliant with
        Gromacs rules.

        It outputs the new, Gromacs formatted, boxvectors and a 
        rotation matrix which can be used to rotate the atom 
        coordinates.

        The file BOXVECTORS contains the new vectors in Gromacs 
        format:

                x1 y2 z3 y1 z1 x2 z2 x3 y3

        The file RMATRIX contains the rotation matrix:

                a11
                a12
                a13
                a21
                a22
                a23
                a31
                a32
                a33

        Use these to perform the rotation on the coordinates:

                newx = x * a11 + y * a12 + z * a13
                newy = x * a21 + y * a22 + z * a23
                newz = x * a31 + y * a32 + z * a33

        Of course these things could have been done more elegantly, 
        but it is again a clear example of ad-hoc programming.
 
AUTHOR
        Tsjerk Wassenaar, University of Groningen.
\n"; 

###
### Command line arguments and initial info
###

# check arguments 
if ( $#ARGV != 8 )
{
  die "$usage_message"; 
}

@v1 = @ARGV[0..2];
@v2 = @ARGV[3..5];
@v3 = @ARGV[6..8];

#
# To start with the procedure it is wise to rearrange lattice vectors in
# descending order (see H. Bekker, J. Comp. Sci., 18(15) p. 1930 (1997):
#
$K = sqrt( $v1[0]*$v1[0] + $v1[1]*$v1[1] + $v1[2]*$v1[2] );
$L = sqrt( $v2[0]*$v2[0] + $v2[1]*$v2[1] + $v2[2]*$v2[2] );
$M = sqrt( $v3[0]*$v3[0] + $v3[1]*$v3[1] + $v3[2]*$v3[2] );
#
#if ($K > $L && $L > $M)
# nothing to be done...
#else
#
if ($K > $M && $M > $L)
{
    @tmp = @v2;
    @v2 = @v3;
    @v3 = @tmp;
}
elsif ($L > $K && $K > $M)
{
    @tmp = @v1;
    @v1 = @v2;
    @v2 = @tmp;
}
elsif ($M > $K && $K > $L)
{
    @tmp = @v1;
    @v1 = @v3;
    @v3 = @v2;
    @v2 = @tmp;
}
elsif ($L > $M && $M > $K)
{
    @tmp = @v1;
    @v1 = @v2;
    @v2 = @v3;
    @v3 = @tmp;
}
elsif ($M > $L && $L > $K)
{
    @tmp = @v1;
    @v1 = @v3;
    @v3 = @tmp;
}

$vpx = $v1[1] * $v2[2] - $v2[1] * $v1[2];
$vpy = $v2[0] * $v1[2] - $v1[0] * $v2[2];
$vpz = $v1[0] * $v2[1] - $v2[0] * $v1[1];

#
# First calculate the final vectors. After the operation the matrix 
# will look like:
#
#   / a11 a12 a13 \       / x1 x2 x* \       / x'1 x'2  0  \
#  /               \     /            \     /               \
# (   a21 a22 a23   ) * (   y1 y2 y*   ) = (    0  y'2  0    )
#  \               /     \            /     \               /
#   \ a31 a32 a33 /       \ z1 z2 z* /       \  0   0  z'3 /
#
#
#   ONLY THE FIRST TWO VECTORS ARE USED FOR CALCULATION
#   A THIRD VECTOR, BEING THE VECTOR PRODUCT IS ADDED,
#   BECAUSE OF CONVENIENCE (ONLY 3X3 MATRICES).
#

#
# x'1 equals the lenght of the vector (x1, y1, z1): 
# sqrt( x1x1 + y1y1 + z1z1 )
#
$xn1 = sqrt( $v1[0]*$v1[0] + $v1[1]*$v1[1] + $v1[2]*$v1[2] );
# $yn1 = 0;
# $zn1 = 0;

#
# x'2 equals the projection of the vector (x2, y2, z2) on (x1, y1, z1):
# sqrt( x2x2 + y2y2 + z2z2 ) * cos alpha
# where alpha is the angle between both vectors:
# cos alpha = (x1x2 + y1y2 + z1z2) / 
# ( sqrt( x1x1 + y1y1 + z1z1 ) + sqrt( x2x2 + y2y2 + z2z2 ) )
#
$xn2 = ( $v1[0]*$v2[0] + $v1[1]*$v2[1] + $v1[2]*$v2[2] ) / $xn1; #$xn1 = length of vector 1

#
# y'2 is equal to square root of ( x2² + y2² + z2² ) - x'2²
# which corresponds to the square root of the difference between the 
# length of vector 2 squared and the projection of that vector on vector
# one squared ( c² = a² - b² )
#
$yn2 = sqrt( $v2[0] * $v2[0] + $v2[1] * $v2[1] + $v2[2] * $v2[2] - $xn2 * $xn2 );
# $zn2 = 0;

# $vnx = 0;
# $vny = 0;
#
# z'3 is equal to the length of the vector product of vectors 1 and two,
# being equal to sqrt( x3² + y3² + z3² )
#
$vnz = sqrt( $vpx*$vpx + $vpy*$vpy + $vpz*$vpz );

print "\n  X' = AX = \n\n"; 
printf "       / %9.2e %9.2e %9.2e \\\n", ($xn1, $xn2, 0);
printf "      / %31s \\\n";
printf "     (   %9.2e %9.2e %9.2e   )\n", (0, $yn2, 0);
printf "      \\ %31s /\n";
printf "       \\ %9.2e %9.2e %9.2e /\n\n", (0, 0, $vnz);

#
# The rotation matrix is calculated as follows:
# ( see http://mathworld.wolfram.com/EulerAngles.html )
#
# X' = AX
#
# A = X' * Inv X
#
# Inv X = ( 1 / Det X ) * Adj X
#
# Adj X = C =
#
# ( c11 c12 c13 )
# (             )
# ( c21 c22 c23 ) =>
# (             )
# ( c31 c32 c33 )
#
$c11 = $v2[1] * $vpz - $v2[2] * $vpy;
$c12 = $v2[2] * $vpx - $v2[0] * $vpz;
$c13 = $v2[0] * $vpy - $vpx * $v2[1];
$c21 = $v1[2] * $vpy - $v1[1] * $vpz;
$c22 = $v1[0] * $vpz - $v1[2] * $vpx;
$c23 = $v1[1] * $vpx - $v1[0] * $vpy;
$c31 = $v1[1] * $v2[2] - $v1[2] * $v2[1];
$c32 = $v1[2] * $v2[0] - $v1[0] * $v2[2];
$c33 = $v1[0] * $v2[1] - $v1[1] * $v2[0];

#
# Determinant = x1*(y2*vpz-z2*vpy) - x2*(y1*vpz-z1*vpy) + vpx*(y1*z2-z1*y2) =
# 
$Det = $v1[0] * $c11 + $v2[0] * $c21 + $vpx * $c31;

#
# Inv (X) = (1 / Det X) * C = D = 
#
# ( d11 d12 d13 )
# (             )
# ( d21 d22 d23 ) =>
# (             )
# ( d31 d32 d33 )
#
$d11 = $c11 / $Det;
$d12 = $c12 / $Det;
$d13 = $c13 / $Det;
$d21 = $c21 / $Det;
$d22 = $c22 / $Det;
$d23 = $c23 / $Det;
$d31 = $c31 / $Det;
$d32 = $c32 / $Det;
$d33 = $c33 / $Det;

#
# A = X' * D = 
#
# ( a11 a12 a13 )
# (             )
# ( a21 a22 a23 ) =>
# (             )
# ( a31 a32 a33 )
#
$a11 = $xn1 * $d11 + $xn2 * $d21;  # = $xn1 * $d11 + $xn2 * $d21 + $xn3 * $d31;
$a12 = $xn1 * $d12 + $xn2 * $d22;  # = $xn1 * $d12 + $xn2 * $d22 + $xn3 * $d32;
$a13 = $xn1 * $d13 + $xn2 * $d23;  # = $xn1 * $d13 + $xn2 * $d23 + $xn3 * $d33;
$a21 = $yn2 * $d21;                # = $yn1 * $d11 + $yn2 * $d21 + $yn3 * $d31;
$a22 = $yn2 * $d22;                # = $yn1 * $d12 + $yn2 * $d22 + $yn3 * $d32;
$a23 = $yn2 * $d23;                # = $yn1 * $d13 + $yn2 * $d23 + $yn3 * $d33;
$a31 = $vnz * $d31;                # = $zn1 * $d11 + $zn2 * $d21 + $zn3 * $d31;
$a32 = $vnz * $d32;                # = $zn1 * $d12 + $zn2 * $d22 + $zn3 * $d32;
$a33 = $vnz * $d33;                # = $zn1 * $d13 + $zn2 * $d23 + $zn3 * $d33;

#
# That's it!
#
# The Matrix (...Sounds like a nice title for a Movie, doesn't it?...):
#
print "  Rotation matrix A = \n\n";
printf "       / %9.2e %9.2e %9.2e \\\n", ($a11, $a12, $a13);
printf "      /%33s\\\n";
printf "     (   %9.2e %9.2e %9.2e   )\n", ($a21, $a22, $a23);
printf "      \\%33s/\n";
printf "       \\ %9.2e %9.2e %9.2e /\n\n", ($a31, $a32, $a33);

#
# The new tensor for the triclinic box,
# where $xm3, $ym3 and $zm3 form the rotated third vector
# of the initial (given) box.
#
$xn3 = $a11*$v3[0] + $a12*$v3[1] + $a13*$v3[2];
$yn3 = $a21*$v3[0] + $a22*$v3[1] + $a23*$v3[2];
$zn3 = $a31*$v3[0] + $a32*$v3[1] + $a33*$v3[2];

#
# To meet the conditions of Gromacs, that zm3 > 0:
#
if ($zn3 < 0)
{
    $xn3 = -$xn3;
    $yn3 = -$yn3;
    $zn3 = -$zn3;
}

#
# To meet the conditions of Gromacs, that
#
# | $xn2 | =< 1/2 * $xn1
#
# | $xn3 | =< 1/2 * $xn1
#
# | $yn3 | =< 1/2 * $yn2
#
# ( Lattice Reduction )
#
$xn2 = $xn2 ? $xn2 - int( $xn2 / $xn1 + $xn2 / ( 2 * abs( $xn2 )) ) * $xn1 : $xn2;
$xn3 = $xn3 ? $xn3 - int( $xn3 / $xn1 + $xn3 / ( 2 * abs( $xn3 )) ) * $xn1 : $xn3;
$yn3 = $yn3 ? $yn3 - int( $yn3 / $yn2 + $yn3 / ( 2 * abs( $yn3 )) ) * $yn2 : $yn3;

print "  New boxvectors (from angstrom -> nm: * 0.1):\n\n";
printf "       / %10.5f %10.5f %10.5f \\\n", ($xn1/10, $xn2/10, $xn3/10);
printf "      /%36s\\\n";
printf "     (   %10.5f %10.5f %10.5f   )\n", (0, $yn2/10, $yn3/10);
printf "      \\%36s/\n";
printf "       \\ %10.5f %10.5f %10.5f /\n\n", (0, 0, $zn3/10);

# And the rotation matrix for rotpdb.pl:
open MATRIX, ">RMATRIX";
printf MATRIX "%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n%f\n", ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33);
close MATRIX;

# Gromacs format:
# (x1 y2 z3 y1 z1 x2 z2 x3 y3)
open BOX, ">BOXVECTORS";
printf BOX "%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n", ($xn1/10, $yn2/10, $zn3/10, 0, 0, $xn2/10, 0, $xn3/10, $yn3/10);
close BOX;

#
# Final check:
#
# AX = X'
#
# $DELTA is the sum over the differences between the determined 
# coordinates and those obtained by matrix multiplication.
# 
# $DELTA shouldn't be larger than 0.000001
# Usually $DELTA will not exceed 1e-10...
# Extreme cases with one vector very small and the other very large
# should be checked ( len1 / len2 > 10^7 ).
#

#$DELTA = abs(($a11 * $v1[0] + $a12 * $v1[1] + $a13 * $v1[2]) - $xn1) + abs(($a11 * $v2[0] + $a12 * $v2[1] + $a13 * $v2[2]) - $xn2) + abs(($a11 * $v3[0] + $a12 * $v3[1] + $a13 * $v3[2]) - $xn3) + abs(($a21 * $v1[0] + $a22 * $v1[1] + $a23 * $v1[2]) - $yn1) + abs(($a21 * $v2[0] + $a22 * $v2[1] + $a23 * $v2[2]) - $yn2) + abs(($a21 * $v3[0] + $a22 * $v3[1] + $a23 * $v3[2]) - $yn3) + abs(($a31 * $v1[0] + $a32 * $v1[1] + $a33 * $v1[2]) - $zn1) + abs(($a31 * $v2[0] + $a32 * $v2[1] + $a33 * $v2[2]) - $zn2) + abs(($a31 * $v3[0] + $a32 * $v3[1] + $a33 * $v3[2]) - $zn3);

#printf "%20.6e\t\t%20.6e\n", ($xn1 / sqrt( $xn2*$xn2 + $yn2*$yn2 ), $DELTA);

#print "The sum of the differences between determined and calculated\n coordinates is $DELTA.\n\n";

