#!/usr/bin/perl

use Math::Trig;

#Script que le as coordenadas cartesianas de uma molecula e fornece sua topologia.
#autor: Evanildo Laceda Jr.
#Intituto de Fisica, USP, Brasil. 

###################################
#Ultima modificacao:  02/07/2012 #
##################################
@comand = @ARGV;

$help =grep(/-h/,@comand);


if ($help==1){

&intro();
exit;

}

chomp($mol1 = "$ARGV[0]");
chomp($mol2 = "$ARGV[1]");



@nom1 = split(/\./,$mol1);
@nom2 = split(/\./,$mol2);
$ext1=$nom1[1];


open(in1,"<$mol1") || die "$!: $mol1";
open(in2,"<$mol2") || die "$!: $mol2";

@linhas1 = <in1>;
@linhas2 = <in2>;

&intro();

print "mol1: $mol1\n";
print "mol2: $mol2\n\n";


$arq_xray = $nom1[0] . "xVSx" . $nom2[0]. ".xraycon";

print "output: $arq_xray\n\n";

open(out,">$arq_xray");

print out "mol1: $mol1\n";
print out "mol2: $mol2\n\n";

print out "output: $arq_xray\n\n";

$i=0;

while ($linhas1[$i] !~ /\[ bonds \]/){ $i++;}

$i+=2;

print out  "\n[ bonds ]\n\n"; 

printf out "%6s %6s  %7s %7s %9s %9s\n",ai,aj,r1,r2,dr,"    (r1 - dr) (nm)";

# print "teste: $linhas1[$i]\n";

$sum_r=0;
while (($linhas1[$i] !~ /\[ angles \]/) || ($linhas1[$i] != "")){


@par1 = split(/ +/,$linhas1[$i]);
@par2 = split(/ +/,$linhas2[$i]);

$r1=$par1[3]; $r2=$par2[3];

$dr=$r2-$r1;
$sum_r += abs($dr);

$r_ajust = $r1 - $dr; 

printf out "%6s %6s %7.3f %7.3f %9.3f %9.4f\n",@par1[1],@par1[2],$r1,$r2,$dr,$r_ajust/10;

$i++;

}

printf out "; soma de \|dr\| = %8.3f\n",$sum_r;
printf "; soma de \|dr\| = %8.3f\n",$sum_r;

$i+=2;

print  out "\n[ angles ]\n\n"; 
printf out "%6s %6s %6s %9s %9s %11s %11s\n",ai,aj,ak,theta1,theta2,dtheta,"   (theta1 - dtheta)";

$sum_th=0;

while (($linhas1[$i] !~ /\[ torsions \]/) || ($linhas1[$i] != / +/)){

@par1 = split(/ +/,$linhas1[$i]);
@par2 = split(/ +/,$linhas2[$i]);

$th1=$par1[4]; $th2=$par2[4];

$dth = $th2-$th1;

$sum_th += abs($dth);
$th_ajust = $th1 - $dth; 

printf out "%6s %6s %6s  %9.2f %9.2f %11.2f %11.2f\n",@par1[1],@par1[2],@par1[3],$th1,$th2,$dth,$th_ajust;


$i++;

}

printf out "; soma de \|dtheta\| = %9.2f\n",$sum_th;
printf  "; soma de \|dtheta\| = %9.2f\n",$sum_th;

$i+=2;


print out  "\n[ torsions ]\n\n"; 
printf out "%4s %4s %4s %4s %9s %9s %11s\n",ai,aj,ak,al,phi1,phi2,dphi;

$sum_phi=0;

while  ($linhas1[$i]){
  
@par1 = split(/ +/,$linhas1[$i]);
@par2 = split(/ +/,$linhas2[$i]);

$phi1=$par1[5]; $phi2=$par2[5];

$dphi = $phi2-$phi1;

if ($dphi > 180.000) {$dphi=360-$dphi;}
if ($dphi < -180.000) {$dphi+=360;}

$sum_phi+=abs($dphi);

printf out "%6s %6s %6s %6s  %9.2f %9.2f %11.2f\n",@par1[1],@par1[2],@par1[3],@par1[4], $phi1,$phi2,$dphi;
$i++

}

printf out "; soma de \|dphi\| = %9.2f\n",$sum_phi;
printf "; soma de \|dphi\| = %9.2f\n\n",$sum_phi;

sub intro {
 print "--------------------------------------------------------------\n";
 print "X-raycon - Contrast Between Molecules in differents geometries\n";
 print "created by Evanildo Lacerda Jr.\n";
 print "--------------------------------------------------------------\n\n";
 
 print "command -> \$ xraycomp [option] [molA1.xray] [molA2.xray] \n\n";
 print "[options]  \n";
 print "-h   \t help\n\n";
}
