#最近需要扫描分子的二面角

第一印象，肯定用内坐标解决了。

内坐标有个问题，一定要自己手动写gjf，否则直接改变二面角的度数会出现结构紊乱的现象。

所以最终还是决定采用直角坐标系

################

Dihe(8-7-4-3)=0 到360度扫描，步长假设是5度。

######

我在这里写下伪代码，等程序写完以后贴上所有的代码。

把74之间的可旋转单键断开，

利用

BFS算法，确定2个平面

http://baike.baidu.com/link?url=a0ZwNALg6phDz8UjeBFfHyFIMX8tgNKd9f8wE6jiTAhsSfWKhyUWm3FLXP6D4wIZ

####################

计算74的向量，并转换为单位向量

vectot(7-4)=(x1,y1,z1);#x1**2+y1**2+z1**2=1

###################

假设74单键旋转n弧度

构建一个4元组

Quad(x,y,z,w)

x=x1*sin（n/2）

y=y1*sin(n/2)

z=z1*sin(n/2)

w=cos(n/2)

#############################

构建一个3X3  的M矩阵

###################################################

M=

[

  [1-2y^2-2z^2,2xy+2wz,2xz-2wy],

  [2xy-2wz,1-2xx-2zz,2yz+2wx],

  [2xz+2wy ,2yz-2wx,1-2xx-2yy]

]

###############################################

[x,y,z] x M 就能得到旋转以后的坐标了

############简陋的测试代码###########################

#!/usr/bin/perl -w

use strict;

use PDL;


sub storecoor
{
   

   my @text=@_;
   
   my %hash;
   my $atomid=1;
   foreach my $line(@text)
   {
#    C    0.0000000000        0.0000000000        0.0000000000
       if($line=~/^\s+\w+\s+(\S+)\s+(\S+)\s+(\S+)/)
       {
           $hash{$atomid}{'x'}=$1;
           $hash{$atomid}{'y'}=$2;
           $hash{$atomid}{'z'}=$3;
       
           $atomid++;
       }
   }
   return %hash;
}



my @graph2=(7,8,9,10,11,17,18,19);
my @graph1=(1,2,3,4,5,6,12,13,14,15,16);
my @dihe=(8,7,4,3);

my %hash; #to store coor
#$hash{$atomid}{x}        $hash{$atomid}{y}        $hash{$atomid}{z}

open FH,"benz_furan.gjf";

my @text=<FH>;
close(FH);

%hash=&storecoor(@text);

my @rotateaxis=($hash{$dihe[1]}{'x'}-$hash{$dihe[2]}{'x'},$hash{$dihe[1]}{'y'}-$hash{$dihe[2]}{'y'},$hash{$dihe[1]}{'z'}-$hash{$dihe[2]}{'z'});
my $distance=sqrt($rotateaxis[0]**2+$rotateaxis[1]**2+$rotateaxis[2]**2);
$rotateaxis[0]/=$distance;
$rotateaxis[1]/=$distance;
$rotateaxis[2]/=$distance;
use constant   PI=>3.14159265358979323846264338327950288419716939937510582097494459;

#print cos(PI/3);

my %hashtran=&transform(60); #graph2  transform 60 degree



   my $atomid=1;
   foreach my $line(@text)
   {
   
#    C    0.0000000000        0.0000000000        0.0000000000
       if($line=~/^(\s+\w+)\s+\S+\s+\S+\s+\S+/)
       {
   

           my $x=sprintf("%15.6f",$hashtran{$atomid}{'x'});

           my $y=sprintf("%15.6f",$hashtran{$atomid}{'y'});
           my $z=sprintf("%15.6f",$hashtran{$atomid}{'z'});
       
         print $1.$x.$y.$z;
         print "\n";
           $atomid++;
           next;
       }
       print $line;
   }











sub transform
{
   my %outhash=%hash;
   my @quadd=($rotateaxis[0]*sin(PI/6),
   $rotateaxis[1]*sin(PI/6),
   $rotateaxis[2]*sin(PI/6),
   cos(PI/6));
 my $x=$quadd[0];my $y=$quadd[1];my $z=$quadd[2];
 my $w=$quadd[3];
   
   my $mat=pdl
   [
      [1-2*$y*$y-2*$z*$z,2*$x*$y+2*$w*$z,2*$x*$z-2*$w*$y],
      [2*$x*$y-2*$w*$z,1-2*$x*$x-2*$z*$z,2*$y*$z+2*$w*$x],
      [2*$x*$z+2*$w*$y,2*$y*$z-2*$w*$x,1-2*$x*$x-2*$y*$y]
   
   
   
   ];
   
   foreach my $graphatomid(@graph2)
   {
       my $atomvec=pdl
       [
         [$hash{$graphatomid}{'x'},$hash{$graphatomid}{'y'},$hash{$graphatomid}{'z'}]
       ];
       
       my $output=$atomvec x $mat;
       $output=~s/\n//g;
       if($output=~/\[\s+\[\s+(\S+)\s+(\S+)\s+(\S+)\]\]/)
       {
           #print $1,$2,$3,"\n";
           $outhash{$graphatomid}{'x'}=$1;
           $outhash{$graphatomid}{'y'}=$2;
           $outhash{$graphatomid}{'z'}=$3;
           
       }
   
   }
   
   return %outhash;
   
   
}
   


=pod





















my @vector1=&vector(@dihe[0..2]);

print join "_",@vector1;
print "\n";

my @vector2=&vector(@dihe[1..3]);

print join "_",@vector2;



sub vector
{
   my @vector=@_;
   my $x1=$hash{$vector[0]}{'x'};#A
   my $y1=$hash{$vector[0]}{'y'};
   my $z1=$hash{$vector[0]}{'z'};
   my $x2=$hash{$vector[1]}{'x'};  #B
   my $y2=$hash{$vector[1]}{'y'};
   my $z2=$hash{$vector[1]}{'z'};
       my $x3=$hash{$vector[2]}{'x'}; #C
   my $y3=$hash{$vector[2]}{'y'};
   my $z3=$hash{$vector[2]}{'z'};
   my @output;
   my $x=1;
   #BA
   my $a1=($x1-$x2);my $b1=($y1-$y2);my $c1=$z1-$z2;
   #BC
   my $a2=($x3-$x2);my $b2=$y3-$y2; my $c2=$z3-$z2;
   
   #AC
   my $a3=($x3-$x1);my $b3=$y3-$y1;my $c3=$z3-$z1;
   
   
   my $y=($a2*$c1-$a1*$c2)/($b1*$c2-$b2*$c1);
   my $z=($a3+$b3*$y)/(-$c3);
 
 push @output,$x,$y,$z;
 return @output;
}



=cut

###################################################################

#####find all the atom connect ##############

sub contact
{
   my %matrix=%{$_[0]};
   my $begin=$_[1];
   my $atomnum=$_[2];
   my @output;

   my @color;#0 白色 1灰色 2黑色
   my @parent;#0
   my @grey;  #to store grey nodes;
   #all atom is white default
   
   
   
   ############adj: which atoms near the atom?#########
#    my @adjatoms=&adj(\%matrix,$head,$atomnum);
sub  adj
{
my %matrix=%{$_[0]};
my $begin=$_[1];

my $atomnum=$_[2];
my @adjatoms;
for(my $i=1;$i<=$atomnum;$i++)
{
  if($matrix{$i}{$begin}==1 || $matrix{$begin}{$i}==1)
  {
      push @adjatoms,$i;
  }
}

return @adjatoms ;
}
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   for(1..$atomnum)
   {
       $color[$_]="white";
   }
   
   #the first atom
   $color[$begin]="grey";
   $parent[$begin]="null";
   
   push @grey,$begin;
   
   while(@grey)
   {
       my $head=shift @grey;
       my @adjatoms=&adj(\%matrix,$head,$atomnum);
       foreach my $atomid(@adjatoms)
       {
           if($color[$atomid]eq "white")
           {
               $color[$atomid]="grey";
               $parent[$atomid]=$head;
               $matrix{$atomid}{$head}=0;  #to avoid  backtrace
               push @grey,$atomid;
           }
       }        
           $color[$head]="black";
   }
   
       for(1..$atomnum)
   {
       if($color[$_]  eq "black")
       {
           push @output,$_;
       }
   }
   
   return @output;
   
}
               
   
#############################################################