# 分子二面角旋转 公式

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

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

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

######

BFS算法，确定2个平面

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

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

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

x=x1*sin（n/2）

y=y1*sin(n/2)

z=z1*sin(n/2)

w=cos(n/2)

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

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

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];
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;

}

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

https://blog.sciencenet.cn/blog-950202-779069.html

## 相关博文

GMT+8, 2021-12-8 12:49