#! /usr/bin/perl -w

# USAGE apex.pl $mol.xyz
# ACA is the Tolman cone angle
# For tempo-R, make sure N in tempo is atom1 

   #####################################################################
   # Author Leaf Lin cylin@rsc.anu.edu.au
   # Work out the steric parameter of a legend.
   # tolman cone angle (ACA).
   # Please see http://rsc.anu.edu.au/~cylin/cone_angle.pdf for the notations and equations for the abbreviations.
   # 
   # This works only for carbon center radicals (legand) if you want to do non-carbon centered ones, please email me.
   # Please make sure you define your cone center. In this file, it's the Cl atom. You can modify the code in the "cone element" section to adjust your cone center.
   #####################################################################

   use Data::Dumper;

$pi = atan2(1,1) * 4;
$piover180 = $pi/180;
$print = 1;

   $nat = 0;
   (@A,@X,@Y,@Z,@B) = ();


   while (<>) {
      next if (m/molecule/);
      last if (m/end/);
      @F = split(/\s+/,$_);
      next unless(m/\S/ && @F > 2);
      $nat++;
      @_ = split(/\s+/,$_);
      ($z,$y,$x,$a) = reverse(@_); 
      $a = ToSymbol($a) if ($a =~ m/\d/);
      #print "$a\n$x  $y  $z\n";
      push(@A,$a); push(@X,$x); push(@Y,$y);  push(@Z,$z);
   }

   for $i (0..$nat-1) {
       for $j (0..$nat-1) {
           $b = ($X[$i]-$X[$j])**2 + ($Y[$i]-$Y[$j])**2 + ($Z[$i]-$Z[$j])**2;
           $b = sprintf("%6.4f",sqrt($b));
           $B[$i][$j] = $b;
       }
   }


   #PrintMatrix($nat,\@A,\@B,0);
   PrintMatrix($nat,\@A,\@B,1);


   #----------------------------------------------------------------------


   sub PrintMatrix {
       my ($n,$A,$B,$bond) = @_;
       my @A = @$A;
       my @B = @$B;
       my $maxcol = $bond ? 80 : 8;
       my ($a1,$a2);
       my @group=() ;
       my @groups=() ;
       my @cxs=() ;
       my ($a,$alpha,$b,$d,$c);
       my @sorted=() ;
       my @longest=() ;
       my @longestx=() ;
       my @widest=() ;
       my @widestx=() ;
       my ($g,$gs);
       my ($theta);
       my @new=();
       my $apexsum;
       my $vsum;
       my $test;
       my ($h, $r);


       my ($ic,$ir,$imin,$imax,$nrows);


       $nrows = int($n/$maxcol);
       $nrows++ if ($nrows*$maxcol != $n); 

       for $ic (1..$nrows) {
           $imin = ($ic-1)*$maxcol + 1;
           $imax = min($imin+$maxcol-1,$n);


           if ($print) {
               $bond ? print "     " : print "    ";

               for $ir ($imin..$imax) {
                   $bond ? printf " %-3s",$A[$ir-1] : 
                   printf "   %-3s   ",$A[$ir-1];
               }
               print "\n";
           } 
           for $ir ($imin..$imax) {
               $count[$ir-1]=0;
           }
           for $ir (1..$n) {
               if ($print) {
                   $bond ? printf "%-3s",$A[$ir-1] : 
                   printf "%-3s",$A[$ir-1];
               }
               for $is ($imin..$imax) {
                   if ($bond) {
                       ($a1,$a2) = ($A[$ir-1],$A[$is-1]);
                       $a1 =~ s/([A-Za-z]+).*/$1/;
                       $a2 =~ s/([A-Za-z]+).*/$1/;
                       if ($B[$is-1][$ir-1] < Bond($a1,$a2) &&
                           $B[$is-1][$ir-1] > 0           ) {
                           if ($print) {
                               print "   X";
                           }
                           $count[$ir-1] ++;
                       }else {
                           if ($print) {
                               print "   .";
                           }
                       }
                   }else {
                       if ($print) {
                           print "  $B[$is-1][$ir-1] ";
                       }
                   }
               }
               if ($print) {
                   print "\n";
               }
           }

           for $ir ($imin..$imax) {
               $new[$ir-1]=1;
               #########################################################
               # 
               # Cone Element
               #
               # Here, change the element to cone element!
               # if ($count[$ir-1] == 2 && $A[$ir-1] eq "O" && $B[1-1][$ir-1] < Bond("n","o")){  # Tempo-R, here manually move N to atom 1
               # if ( $A[$ir-1] eq "Br"){  # Br-R
               # if ($count[$ir-1] == 2 && $A[$ir-1] eq "S" ){  # DTC(Me)-R
               #########################################################
                       if ( $A[$ir-1] eq "Cl"){  # Cl-R
                   $cc = $ir;
                   if ($print) {
                       print "this is the cone center $A[$ir-1] $ir\n";
                   }
               }
           }
           for $is ($imin..$imax) {
               ($a1,$a2) = ($A[$cc-1],$A[$is-1]);
               if ($B[$is-1][$cc-1] < Bond($a1,$a2) &&
                   $B[$is-1][$cc-1] > 0 && $count[$is-1] == 4 ) {
                   if ($print) {
                       print "this is the radical center $A[$is-1] $is \n";
                   }
                   $rc = $is;
               }
           }
           for $is ($imin..$imax) {
               ($a1,$a2) = ($A[$rc-1],$A[$is-1]);
               if ($B[$is-1][$rc-1] < Bond($a1,$a2) &&
                   $B[$is-1][$rc-1] > 0           ) {
                   if ($is != $cc){
                       #print "group $A[$is-1] [$is-1] \n";
                       push (@group,$is);
                       $new[$is-1]=0;
                   }
               }
           }

           for $g(0..$#group){
               for $is ($imin..$imax) {
                   ($a1,$a2) = ($A[$group[$g]-1],$A[$is-1]);
                   if ($B[$is-1][$group[$g]-1] < Bond($a1,$a2) &&
                       $B[$is-1][$group[$g]-1] > 0           ) {
                       if ($is != $cc and $is != $rc){
                           #            print "group $g has $A[$is-1] [$is-1] \n";
                           push (@{$groups[$g]},$is);

                           $new[$is-1]=0;
                       }
                   }
               }
               push (@{$groups[$g]},$group[$g]);
           }

           for $i(1..$n){
               if ($i != $cc && $i != $rc){
                   for $g(0..$#group){
                       for $gs(0..$#{$groups[$g]}){
                           if (connected($i-1,$groups[$g][$gs]-1) and $new[$i-1] != 0){
                               push (@{$groups[$g]},$i);
                               $new[$i-1] = 0;
                           }
                       }
                   }
               }
           }

           for $g(0..$#group){
               if ($print){
                   print "group $g has ";
               }
               for $gs(0..$#{$groups[$g]}){
                   if ($print){
                       print "$groups[$g][$gs]\t";
                   }
               }
               if ($print){
                   print "\n";
               }
               $test += $#{$groups[$g]}+1;
           }
           if ($test != $n-2){
               print "Warning! Total group number ", $test+2," is not the same as the legand number! $n \n";
           }
           for $g(0..$#group){
               $longest[$g] = 1.;
               $widest[$g] = 0.;
               $narrowest[$g] = 2.;
               for $gs(0..$#{$groups[$g]}){
                   $x = $groups[$g][$gs];
                   $d[$x-1] = $B[$x-1][$cc-1];
                   if ($print){
                       print "bond $A[$cc-1]$cc $A[$x-1]$x = $B[$x-1][$cc-1]\n";
                   }
                   if ($d[$x-1] > $longest[$g]){
                       $longest[$g] = $d[$x-1];
                       $longestx[$g] = $x-1;
                   }

                   if (dot($rc-1,$cc-1,$x-1)/( $B[$cc-1][$x-1]*$B[$cc-1][$rc-1] ) < 1.0){
                       $alpha[$x-1] = acos(dot($rc-1,$cc-1,$x-1)/( $B[$cc-1][$x-1]*$B[$cc-1][$rc-1] ) );
                   }
                   else{
                       $alpha[$x-1] = 0.0;
                   }
                   if ($print){
                       print "angle $rc $cc $x =", $alpha[$x-1]/$piover180,"\n";

                       print "vdw $A[$x-1]$x =", vdwradii($A[$x-1]),"\n";
                   }

                   $a[$x-1] = asin(vdwradii($A[$x-1])/$d[$x-1]);

                   if ($a[$x-1]+$alpha[$x-1] > $widest[$g]){
                       $widest[$g] = $a[$x-1]+$alpha[$x-1];
                       $widestx[$g] = $x-1;
                   }
                   if ($alpha[$x-1] < $narrowest[$g]){
                       $narrowest[$g] = $alpha[$x-1];
                       $narrowestx[$g] = $x-1;
                   }
                   $c[$x-1] = $a[$x-1]+$alpha[$x-1];
                       print "a+alpha ",$x," is $c[$x-1]\n";
               }

           }
           $widestL = 0.0;
           $longestL = 1.0;
           for $g(0..$#group){
               if ($widest[$g] > $widestL){
                   $widestL = $widest[$g];
               }
               if ($longest[$g] > $longestL){
                   $longestL = $longest[$g];
               }

               print "group ", $g+1, " has longest bond = $longest[$g] and apex angle = $widest[$g]\t";

               $sinapex += sin($widest[$g]);
               # height = (d_y+d_y') * cos(0.5*b_y)
               # radii  = height * tan(0.5*c_x)


               $b = abs(0.5*$c[$widestx[$g]]-$alpha[$longestx[$g]]);
               $dyp = vdwradii($A[$longestx[$g]])/cos($b);

               $h = ($longest[$g]+$dyp) * cos($b);
               $r = $h * tan( 0.5*$widest[$g]);
               print "gives volume = ",$pi*$h*$r**2/3.,"\n";

               $v[$g]=$pi*$h*$r**2/3.;
               $apexsum += $widest[$g];
               $vsum +=$v[$g] ;

#               $bL = $alpha[$longestx[$g]];
#              $dypL = vdwradii($A[$longestx[$g]])/cos($bL);
#              $hL = ($longest[$g]+$dypL) * cos($bL);
#              $rL = $hL * tan( $widest[$g]);
#              $lcv[$g]=$pi*$hL*$rL**2/3.;

               if ($narrowestx[$g] ne $widestx[$g]){
                   $axy = acos(dot($narrowestx[$g],$cc-1,$widestx[$g])/( $B[$cc-1][$narrowestx[$g]]*$B[$cc-1][$widestx[$g]] ) ) + $a[$narrowestx[$g]]+$a[$widestx[$g]];
               }
               else{
                   $axy = $a[$narrowestx[$g]]+$a[$widestx[$g]];
               }

               $sga += $axy;
               $bxy = abs(0.5*$axy-$alpha[$longestx[$g]]);
               $dxyp = vdwradii($A[$longestx[$g]])/cos($bxy);
               $hxy = ($longest[$g]+$dxyp) * cos($bxy);
               $rxy = $hxy * tan( 0.5*$axy);
               $sgv += $pi*$hxy*$rxy**2/3.;

           }

           # $mcv=1.;
           #for $g(0..$#group){

               #if ($lcv[$g] > $mcv){
               #    $mcv = $lcv[$g];
               #}
               #$acv+=$lcv[$g];
               #}
               #$acv = $acv/($#group+1);


           print "SSA = $sinapex\n";
           print "MCA = ",2.*$widestL,"\n";
           print "ACA = ",2*$apexsum/($#group+1),"\n";
           print "SCA = $apexsum\n"; 
           print "SGA = $sga\n";

#          print "MCV = $mcv\n";
#          print "ACV = $acv\n";
           print "SCV = $vsum\n";
           print "SGV = $sgv\n";


       }


   } 

   sub numerically { $a <=> $b; }

   sub vdwradii{
       my ($a) = @_;

        #van der Waals radii are taken from A. Bondi, J. Phys. Chem., 68,
        #441-452, 1964, except the value for H, which is taken from R. S.
        #Rowland & R. Taylor, J. Phys. Chem., 100, 7384 - 7391, 1996. Radii
        #that are not available in either of these publications have RvdW =
        #2.00 A.

    my %vdW = ( h  => 1.09, he => 1.40, li => 1.82, be => 2.00, b  => 2.00,
                c  => 1.70, n  => 1.55, o  => 1.52, f  => 1.47, ne => 1.54,
                na => 2.27, mg => 1.73, al => 2.00, si => 2.10, p  => 1.80,
                s  => 1.80, cl => 1.75, ar => 1.88, k  => 2.75, ca => 2.00,
                sc => 2.00, ti => 2.00, v  => 2.00, cr => 2.00, mn => 2.00,
                fe => 2.00, co => 2.00, ni => 1.63, cu => 1.40, zn => 1.39,
                ga => 1.87, ge => 2.00, as => 1.85, se => 1.90, br => 1.85,
                kr => 2.02, x  => 0.00, i  => 1.98 );
       return $vdW{lc($a)};
   }

   sub tan { sin($_[0]) / cos($_[0]) }
   sub acos { atan2( sqrt(1 - $_[0] * $_[0]), $_[0] ) }
   sub asin { atan2($_[0], sqrt(1 - $_[0] * $_[0])) }

   sub connected {
       my ($a,$b) = @_;
       if ( $B[$a][$b] < Bond($A[$a],$A[$b]) && $B[$a][$b] > 0 ){
           return 1;
       }
       else{
           return 0;
       }
   }
   sub dot {
       my ($a,$b,$c) = @_;
       return ($X[$a]-$X[$b])*($X[$c]-$X[$b])+($Y[$a]-$Y[$b])*($Y[$c]-$Y[$b])+($Z[$a]-$Z[$b])*($Z[$c]-$Z[$b]);
   }

   sub min {
       my ($a,$b) = @_;
       return $b < $a ? $b : $a;
   }
   sub max {
       my ($a,$b) = @_;
       return $b > $a ? $b : $a;
   }


   sub Bond {
       my ($a1,$a2) = @_;


       #CSD covalent radii (Rcov) are used to assign bonded connections in
       #crystal structures. A bonded connection is established if the distance
       #(d) between two elements A and B is in the range from Rcov(A) + Rcov(B)
       #- t to Rcov(A) + Rcov(B) + t, where the tolerance (t) is usually set to
       #0.4 A. The tolerance and Rcov values can be altered by CSD Editors to
       #force bonds in certain situations. The Rcov values have been refined
       #over the years to minimise the need for editorial changes for specific
       #structures. Elements not yet encountered in the CSD have Rcov = 1.50 A.

  %cov = ( h  => 0.23, he => 1.50, li => 0.68, be => 0.35, b  => 0.83,
           c  => 0.68, n  => 0.68, o  => 0.68, f  => 0.64, ne => 1.50,
           na => 0.97, mg => 1.10, al => 1.35, si => 1.20, p  => 1.05,
           s  => 1.02, cl => 0.99, ar => 1.51, k  => 1.33, ca => 0.99,
           sc => 1.44, ti => 1.47, v  => 1.33, cr => 1.35, mn => 1.35,
           fe => 1.34, co => 1.33, ni => 1.50, cu => 1.52, zn => 1.45,
           ga => 1.22, ge => 1.17, as => 1.21, se => 1.22, br => 1.21,
           kr => 1.50, x => 0.00 , i  => 1.40);

       return $cov{lc($a1)}+$cov{lc($a2)} + 0.4 ;

   }


   sub ToSymbol {
       my $n = shift;
       my @sym = qw/H He Li Be B C N O F Ne Na Mg Al Si P S Cl Ar 
       K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr/;
       return $sym[$n-1];
   }
