Новосибирский институт органической химии им. Н.Н. Ворожцова СО РАН
Лаборатория изучения механизмов органических реакций
 |
 |
 |
 |
 |
 |
rotepi
#!/usr/bin/perl -ws #use Data::Dump 'pp'; #$Data::Dump::LINEWIDTH = 80; our ($h,$help,$rot,$epi,$angle,$reflect,$xtb); if ($h || $help) { (my $program = $0) =~ s/^.*[\/\\]//; print <<HELP; Usage: $program -rot=1,2 -angle=90 file.xyz Rotates a fragment of a molecule in xyz format around a given bond. Dependencies: Grimme's xtb (for -xtb option) Options: -rot=i,j,k atomic numbers i and j defined the rotatable bond. The bond must not be part of a cycle. k is optional and may be found automatically. -rot does not work when j and k are in the same cycle. -angle rotation angle in degrees (default 180). -reflect reflection of transformed part about ijk plane Set -angle=0 for only reflection -epi=i,j,k epimerization about i (rotation of a fragment including atoms j and k around the bissector of j-i-k angle. j and k are optional and may be found automatically. -angle and -reflect work for -epi too. -xtb=n Rough optimization with xtb (for linux). This option is useful if the transformation of the molecule results in very close atoms (there will be a warning about this). \'-xtb\' is \'-xtb=10\' (10 optimization steps). HELP exit; } # Ковалентные радиусы из Dalton Trans., 2008, 2832-2838 my %radius = qw( H 0.31 He 0.28 Li 1.28 Be 0.96 B 0.84 C 0.73 N 0.71 O 0.66 F 0.57 Ne 0.58 Na 1.66 Mg 1.41 Al 1.21 Si 1.11 P 1.07 S 1.05 Cl 1.02 Ar 1.06 K 2.03 Ca 1.76 Sc 1.70 Ti 1.60 V 1.53 Cr 1.39 Mn 1.50 Fe 1.44 Co 1.38 Ni 1.24 Cu 1.32 Zn 1.22 Ga 1.22 Ge 1.20 As 1.19 Se 1.20 Br 1.20 Kr 1.16 Rb 2.20 Sr 1.95 Y 1.90 Zr 1.75 Nb 1.64 Mo 1.54 Tc 1.47 Ru 1.46 Rh 1.42 Pd 1.39 Ag 1.45 Cd 1.44 In 1.42 Sn 1.39 Sb 1.39 Te 1.38 I 1.39 Xe 1.40 Cs 2.44 Ba 2.15 La 2.07 Ce 2.04 Pr 2.03 Nd 2.01 Pm 1.99 Sm 1.98 Eu 1.98 Gd 1.96 Tb 1.94 Dy 1.92 Ho 1.92 Er 1.89 Tm 1.90 Yb 1.87 Lu 1.87 Hf 1.75 Ta 1.70 W 1.62 Re 1.51 Os 1.44 Ir 1.41 Pt 1.36 Au 1.36 Hg 1.32 Tl 1.45 Pb 1.46 Bi 1.48 Po 1.40 At 1.50 Rn 1.50 Fr 2.60 Ra 2.21 Ac 2.15 Th 2.06 Pa 2.00 U 1.96 Np 1.90 Pu 1.87 Am 1.80 Cm 1.69 ); #my $ij = shift; die "You must give -rot=i,j or -epi=i,j,k option\n" unless $rot or $epi; die "Inconsistent options -rot and -epi\n" if $rot and $epi; $angle = 180 unless defined $angle; my @mols = read_molden(); my $xtb_x; if ($xtb) { $xtb_x = `which xtb`; chomp $xtb_x; if (-x $xtb_x) { $xtb = 10 if $xtb==1; open L, '>', "tmp_xtb.inp" or die "Can't write to tmp_xtb.inp :$!\n"; print L "\$opt\n maxcycle=$xtb\n"; close L; } else { warn "No xtb found\n"; undef $xtb; } } MOL: foreach my $mol (@mols) { do {warn "less then 4 atoms\n"; next} if @$mol<5; my ($topo,%bonds) = find_bonds($mol); #pp $topo; pp %bonds; exit; my %oldbonds = copy_hh(%bonds); #$topo =~ s/;/\n/g; print "$topo\n"; if ($xtb) { write_molden($mol, "tmp_old.xyz"); `$xtb_x topo tmp_old.xyz > /dev/null 2>&1`; } ####################### Rotation ########################### if ($rot) { # Rotatable bonds my ($i,$j,$k) = split m/\D+/, $rot; die "Invalid -rot=i,j option\n" unless $i && $j; my @rot_bonds; my %bonds_rot = copy_hh(%bonds); foreach (keys %bonds) { my @ar = keys %{$bonds{$_}}; if (@ar==1) { delete $bonds_rot{$_}{$ar[0]}; delete $bonds_rot{$ar[0]}{$_}; } } foreach my $i (sort {$a<=>$b} keys %bonds_rot) { foreach my $j (sort {$a<=>$b} keys %{$bonds_rot{$i}}) { if ($i > $j) { delete $bonds_rot{$i}{$j}; next; } if (part($i,$j,%bonds)==@$mol-1) { delete $bonds_rot{$i}{$j}; next; } #print "$i-$j\n"; push @rot_bonds, "$i-$j"; } } if (! exists $bonds{$i}{$j}) { warn "Can not rotate: $i-$j does not appear to be a bond.\n"; warn "Rotatable bonds: @rot_bonds\n"; next; } my @part_j = part($i,$j,%bonds); my @part_i = part($j,$i,%bonds); #print "@part_j\n@part_i\n"; if (@part_j==@$mol-1) { warn "Can not rotate: bond $i-$j seems to be part of a cycle.\n"; warn "Rotatable bonds: @rot_bonds\n"; next; } unless ($k) { foreach (sort {$a<=>$b} keys %{$bonds{$j}}) { next if $_==$i; $k = $_; last; } } #print "$i,$j,$k\n"; re_orientation($mol,$i,$j,$k); #write_molden($mol); my $mol_j; $mol_j->[0] = {}; foreach (@part_j) { push @$mol_j, $mol->[$_]; $mol->[$_][3] *= -1 if $reflect; } rot_x($mol_j,$angle); # Ссылки в $mol_j - из $mol, поэтому поворот и в $mol тоже } ####################### Epimerizaion ########################### elsif ($epi) { my ($i,$j,$k) = split m/[,-]/, $epi; do {warn "Invalid -epi=$epi option\n"; next} if $i && $j && !$k; my $N = @$mol-1; if ($i && !$j && !$k) { # Only one atom is given my @near = find_near($mol,$i); #print "@near\n"; exit; my @jk; foreach (@near) { if (part($i,$_,%bonds)<$N) { push @jk, $_; } } #print "@jk\n"; exit; if (@jk>=2) { ($j,$k) = @jk; } else { # Try spiro SPIRO: for (my $n=0; $n<@near; $n++) { my $jj = $near[$n]; for (my $m=$n+1; $m<@near; $m++) { my $kk = $near[$m]; my %bbb = copy_hh(%bonds); delete $bbb{$i}{$kk}; delete $bbb{$kk}{$i}; my @part_jk = part($i,$jj,%bbb); if (@part_jk<$N && grep {m/^$kk$/} @part_jk) { ($j,$k) = ($jj,$kk); last SPIRO; } else { warn "Can't find j and k for i=$i\n"; next MOL; } } } } } #warn "$i,$j,$k\n"; if (! exists $bonds{$i}{$j}) { warn "$i-$j does not appear to be a bond.\n"; next; } if (! exists $bonds{$i}{$k}) { warn "$i-$k does not appear to be a bond.\n"; next; } my @part_j = part($i,$j,%bonds); #pp @part_j; my @part_k = part($i,$k,%bonds); #print "@part_j\n@part_k\n"; my @part_jk; if (@part_j<$N && @part_k<$N) { # Both i-j and i-k start open-chains push @part_jk, @part_j, @part_k; } elsif (@part_j==$N && @part_k==$N) { # Both i-j and i-k are in the same cycle - try spiro my %bbb = copy_hh(%bonds); delete $bbb{$i}{$k}; delete $bbb{$k}{$i}; @part_jk = part($i,$j,%bbb); unless (@part_jk<$N && grep {m/^$k$/} @part_jk) { warn "Unsuitable topology for -epi=$epi\n"; next MOL; } } else { warn "Unsuitable topology for -epi=$epi\n"; next MOL; } #print "@part_jk\n"; # XX on the bisector of the angle j-i-k re_orientation($mol,$i,$j,$k); my $xj = $mol->[$j][1]; my $xk = $mol->[$k][1]; my $yk = $mol->[$k][2]; # unit segment on i-j my ($x1j,$y1j) = (1,0); # unit segment on i-k my $x1k = 1/sqrt(1+($yk/$xk)**2); $x1k *= -1 if $xk < 0; my $y1k = $x1k*$yk/$xk; my $x = ($x1j+$x1k)/2; my $y = ($y1j+$y1k)/2; my $XX = ['XX', $x, $y, 0 ]; # # XX at the middle of j,k # my $XX = ['XX', ($mol->[$j][1] + $mol->[$k][1])/2, # ($mol->[$j][2] + $mol->[$k][2])/2, # ($mol->[$j][3] + $mol->[$k][3])/2 ]; push @$mol, $XX; #write_molden($mol); exit; re_orientation($mol,$i,@$mol-1,$j); my $mol_jk; $mol_jk->[0] = {}; foreach (@part_jk) { push @$mol_jk, $mol->[$_]; $mol->[$_][3] *= -1 if $reflect; } rot_x($mol_jk,$angle); # Ссылки в $mol_j - из $mol, поэтому поворот и в $mol тоже pop @$mol; } ############################################################# foreach (keys %{$mol->[0]}) { if ($_ ne 'Charge' or $_ ne 'Mult') { delete $mol->[0]{$_}; } } write_molden($mol) unless $xtb; (my $new_topo,%new_bonds) = find_bonds($mol); my @overlaps; foreach my $i (keys %new_bonds) { foreach my $j (keys %{$new_bonds{$i}}) { next if $i > $j; next if exists $oldbonds{$i}{$j}; push @overlaps, sprintf("%s %.2f","$mol->[$i][0]$i-$mol->[$j][0]$j",$new_bonds{$i}{$j}); } } if (@overlaps) { warn "New short distances: ", join(', ',@overlaps),"\n"; } # if ($new_topo ne $topo) { # warn "Topology has changed!\nOld: $topo\nNew: $new_topo\n"; # } if ($xtb) { write_molden($mol, "tmp_new.xyz"); `$xtb_x tmp_new.xyz --gff --silent --opt -I tmp_xtb.inp`; (my $xtb_mol) = read_molden('xtbopt.xyz'); #pp $xtb_mol; $xtb_mol->[0] = {}; $xtb_mol->[0]{Charge} = $mol->[0]{Charge} if exists $mol->[0]{Charge}; $xtb_mol->[0]{Mult} = $mol->[0]{Mult} if exists $mol->[0]{Mult}; write_molden($xtb_mol); # my @args = ($xtb_x, "tmp_new.xyz", '--gff', '--silent', '--opt', '-I', 'tmp_xtb.inp'); # #print "Run `@args`\n"; # system(@args) == 0 or die "System `@args` failed: $?\n"; } } foreach (qw/gfnff_charges gfnff_topo tmp_new.xyz tmp_old.xyz tmp_xtb.inp xtbopt.log xtbopt.xyz .xtboptok/) { -e && unlink; } sub find_bonds { # External hash %radius (Ковалентные радиусы из Dalton Trans., 2008, 2832-2838) my $mol = shift; my %bonds; my $mult = 1.3; for (my $i=1; $i<@$mol; $i++) { for (my $j=$i+1; $j<@$mol; $j++) { my $dist = dist($mol,$i,$j); if ($dist < ($radius{$mol->[$i][0]}+$radius{$mol->[$j][0]})*$mult) { $bonds{$i}{$j} = $dist; $bonds{$j}{$i} = $dist; } } } #pp %bonds; my ($topo,%topo); foreach my $at (sort {$a<=>$b} keys %bonds) { foreach (sort {$a<=>$b} keys %{$bonds{$at}}) { next if exists $topo{$_}{$at}; $topo{$at}{$_} = 1; } } while (my ($k,$v) = each %topo) { delete $topo{$k} unless (keys %$v); } #pp %topo; foreach (sort {$a<=>$b} keys %topo) { $topo .= ("$_-" . join(',', sort {$a<=>$b} keys %{$topo{$_}}) . ";"); } $topo =~ s/;$//; return $topo, %bonds; } sub find_near { # External hash %radius (Ковалентные радиусы из Dalton Trans., 2008, 2832-2838) my $mol = shift; my $at = shift; my @near; my $mult = 1.3; for (my $i=1; $i<@$mol; $i++) { next if $i==$at; my $dist = dist($mol,$i,$at); if ($dist < ($radius{$mol->[$i][0]}+$radius{$mol->[$at][0]})*$mult) { push @near, $i; } } return @near; } sub part { # Задается связь i-j и хэша хэшей связей # Выдает массив номеров атомов от j и дальше my $i = shift; my $j = shift; my %bonds = copy_hh(@_); #pp %bonds; delete $bonds{$i}{$j} if exists $bonds{$i}{$j}; delete $bonds{$j}{$i} if exists $bonds{$j}{$i}; #pp %bonds; my @part = ($j); my @shell = ($j); while (@shell) { my @shell_; foreach my $at (@shell) { foreach (keys %{$bonds{$at}}) { push @shell_, $_; delete $bonds{$at}{$_}; delete $bonds{$_}{$at}; } } #print "@shell\n"; # @shell = sort {$a==$i ? 1 : -1} @shell_; # sort - for -epi # for -epi: i should be at the end of shell @shell = grep {$_ != $i} @shell_; push @shell, $i if @shell != @shell_; push @part, @shell; } my %h; @part = grep {! $h{$_}++} @part; #print "@part\n"; return @part; } #sub part2 { # # Задаются связи i-j, i-k (как i,j,k) и хэш хэшей связей # # Выдает массив номеров атомов от j,k и дальше # my $i = shift; # my $j = shift; # my $k = shift; # my %bonds = copy_hh(@_); # #pp %bonds; # # delete $bonds{$i}{$j}; # delete $bonds{$j}{$i}; # delete $bonds{$i}{$k}; # delete $bonds{$k}{$i}; # #pp %bonds; # my @part = ($j); # my @shell = ($j); # while (@shell) { # my @shell_; # foreach my $at (@shell) { # foreach (keys %{$bonds{$at}}) { # push @shell_, $_; # delete $bonds{$at}{$_}; # delete $bonds{$_}{$at}; # } # } # #print "@shell\n"; ## @shell = sort {$a==$i ? 1 : -1} @shell_; # sort - for -epi # # for -epi: i should be at the end of shell # @shell = grep {$_ != $i} @shell_; # push @shell, $i if @shell != @shell_; # push @part, @shell; # } # my %h; # @part = grep {! $h{$_}++} @part; # #print "@part\n"; # return @part; #} sub copy_hh { # Независимая копия хэша хэшей my %hh = @_; my %new_hh; while (my ($k,$v) = each %hh) { $new_hh{$k} = {%$v}; } return %new_hh; } # Читает xyz. Параметры - имена xyz-файлов. Если параметров нет, то <>. # Возвращает массив найденных молекул. sub read_molden { local @ARGV = @_ ? @_ : @ARGV; my $num = qr/-?\d+(?:\.\d+)?/; my @mols; my $line; LOOP: while ($line || defined($line = <>)) { #print $line; if ($line =~/^\s*(\d+)\s*$/) { my @mol; my $N = $1; last LOOP if eof(); next LOOP if eof(ARGV); $line = <>; if ($line =~ /^\s*(.*?)\s*$/) { $mol[0]{Title} = $1; } if ($line =~ /\s($num)\s/o) { $mol[0]{Energy} = $1; } if ($line =~ /Symmetry\s+(\w+)/) { $mol[0]{Symmetry} = $1; } #($mol[0]) = $line =~ /($num)/o; for (my $i=1; $i<=$N; $i++) { last LOOP if eof(); next LOOP if eof(ARGV); $line = <>; #print $line; if ($line =~ /^\s*([A-Za-z]{1,2})\s+($num)\s+($num)\s+($num)\s*(.*)/io) { $mol[$i] = [$1,$2,$3,$4,$5]; $mol[$i][0] = ucfirst(lc $mol[$i][0]); } else { next LOOP; } } push @mols, \@mol; last LOOP if eof(); } else { undef $line; } } return @mols; } # Печатает xyz. Параметры -- список молекул, в конце м.б. имя файла. # Если последний элемент списка - имя файла (не ссылка на массив), # то печать в этот файл, иначе - на stdout sub write_molden { # my $fh; my $file; if (ref($_[-1]) ne 'ARRAY') { $file = pop @_; } my @ar; foreach my $mol (@_) { my $N = $#{$mol}; push @ar, " $N\n"; push @ar, "$mol->[0]{Name} " if exists $mol->[0]{Name}; #print $fh " Symmetry $mol->[0]{Symmetry} " if $mol->[0]{Symmetry}; push @ar, "\n"; for (my $i=1; $i<=$N; $i++) { my ($atom,$x,$y,$z) = @{$mol->[$i]}; push @ar, sprintf " %-2s %12.8f %12.8f %12.8f", $atom, $x, $y, $z; #printf $fh uc($atom) eq 'H' ? " %10.3f" : " %9.2f" , $ppm if $ppm; push @ar, "\n"; } } if ($file) { open L, '>', $file or die "Can't write to $file: $!\n"; print L @ar; close L; } else { print @ar; } } #sub write_molden { ## my $fh = \*STDOUT; # my $fh = STDOUT; # my $to_file; # if (ref($_[-1]) ne 'ARRAY') { # $to_file = 1; # my $file = pop @_; # open $fh, '>', $file or die "Can't write to $file: $!\n"; # } # # foreach my $mol (@_) { # my $N = $#{$mol}; # print $fh " $N\n"; # print $fh "$mol->[0]{Name} " if exists $mol->[0]{Name}; # #print $fh " Symmetry $mol->[0]{Symmetry} " if $mol->[0]{Symmetry}; # print $fh "\n"; # for (my $i=1; $i<=$N; $i++) { # my ($atom,$x,$y,$z) = @{$mol->[$i]}; # printf $fh " %-2s %12.8f %12.8f %12.8f", $atom, $x, $y, $z; # #printf $fh uc($atom) eq 'H' ? " %10.3f" : " %9.2f" , $ppm if $ppm; # print $fh "\n"; # } # } # $fh = \*STDOUT; # #close $fh if $to_file; #} # my ($x,$y,$z) = get_xyz($mol); # Возвращает список ссылок на массивы 1..$N (@x,@y,@z) молекулы # (нулевые элементы пустые) sub get_xyz { my $mol = shift; my (@x,@y,@z); my $N = $#{$mol}; for (my $i=1; $i<=$N; $i++) { $x[$i] = $mol->[$i][1]; $y[$i] = $mol->[$i][2]; $z[$i] = $mol->[$i][3]; } return (\@x,\@y,\@z); } # put_xyz($mol,\@x,\@y,\@z); # Помещает координаты из массивов 1..$N (@x,@y,@z) в молекулу sub put_xyz { my ($mol,$x,$y,$z) = @_; my $N = $#{$mol}; return undef if $#{$x} != $N || $#{$y} != $N || $#{$z} != $N; for (my $i=1; $i<=$N; $i++) { $mol->[$i][1] = $x->[$i]; $mol->[$i][2] = $y->[$i]; $mol->[$i][3] = $z->[$i]; } 1 } # centre_to_atom($mol,$i) # Помещает центр координат молекулы на атом $i sub centre_to_atom { my ($mol,$i) = @_; my $N = $#{$mol}; return undef if $i > $N; my ($xi,$yi,$zi) = ($mol->[$i][1],$mol->[$i][2],$mol->[$i][3]); for (my $n=1; $n<=$N; $n++) { $mol->[$n][1] -= $xi; $mol->[$n][2] -= $yi; $mol->[$n][3] -= $zi; } 1 } # ($A,$B,$C) = plane_normal($mol,$i,$j,$k); # Возвращает вектор нормаль к плоскости, проходящей через атомы $i, $j, $k. # Если атомов больше трех, плоскость через них проводится по наименьшим квадратам. sub plane_normal { my ($mol,@nums) = @_; my ($x,$y,$z) = get_xyz($mol); my @x = @$x; my @y = @$y; my @z = @$z; my ($A,$B,$C); return undef if @nums<3; if (@nums == 3) { my ($i,$j,$k) = @nums; #print "@nums\n"; $A = ($y[$j]-$y[$i])*($z[$k]-$z[$i])-($z[$j]-$z[$i])*($y[$k]-$y[$i]); $B = ($z[$j]-$z[$i])*($x[$k]-$x[$i])-($x[$j]-$x[$i])*($z[$k]-$z[$i]); $C = ($x[$j]-$x[$i])*($y[$k]-$y[$i])-($y[$j]-$y[$i])*($x[$k]-$x[$i]); #print "@nums $A $B $C\n"; } else { my ($X,$Y,$Z,$XX,$YY,$ZZ,$XY,$YZ,$XZ); foreach (@nums) { $X += $x[$_]; $Y += $y[$_]; $Z += $z[$_]; $XX += $x[$_]*$x[$_]; $YY += $y[$_]*$y[$_]; $ZZ += $z[$_]*$z[$_]; $XY += $x[$_]*$y[$_]; $YZ += $y[$_]*$z[$_]; $XZ += $z[$_]*$x[$_]; } # A*XX + B*XY + C*XZ = -X XX XY XZ | -X -X XY XZ XX -X XZ XX XY -X # A*XY + B*YY + C*YZ = -Y XY YY YZ | -Y -Y YY YZ XY -Y YZ XY YY -Y # A*XZ + B*YZ + C*ZZ = -Z XZ YZ ZZ | -Z -Z YZ ZZ XZ -Z ZZ XZ YZ -Z $A = $X*($YZ*$YZ-$YY*$ZZ) + $Y*($XY*$ZZ-$YZ*$XZ) + $Z*($YY*$XZ-$XY*$YZ); $B = $X*($XY*$ZZ-$XZ*$YZ) + $Y*($XZ*$XZ-$XX*$ZZ) + $Z*($XX*$YZ-$XY*$XZ); $C = $X*($YY*$XZ-$YZ*$XY) + $Y*($YZ*$XX-$XY*$XZ) + $Z*($XY*$XY-$YY*$XX); } my $normal = sqrt($A**2 + $B**2 +$C**2); $A /= $normal; $B /= $normal; $C /= $normal; #print "@nums $A $B $C\n"; return ($A,$B,$C); } # re_orientation($mol,$i,$j,$k); # Переориентирует молекулу так чтобы у атомов $i,$j,$k были координаты # $i: 0.00 0.00 0.00 # $j: f.ff 0.00 0.00 # $k: f.ff f.ff 0.00 # Dependecies: plane_normal, colinearity, put_xyz sub re_orientation { my ($mol,$i,$j,$k) = @_; centre_to_atom($mol,$i); #&write_molden; my ($x,$y,$z) = get_xyz($mol); my @x = @$x; my @y = @$y; my @z = @$z; my $xj = sqrt($x[$j]**2 + $y[$j]**2 + $z[$j]**2); my $xk = ($x[$k]*$x[$j] + $y[$k]*$y[$j] + $z[$k]*$z[$j]) / $xj; my $yk = sqrt($x[$k]**2 + $y[$k]**2 + $z[$k]**2 - $xk**2); my ($A1,$B1,$C1) = plane_normal($mol,$i,$j,$k); for (my $n=1; $n<=$#{$mol}; $n++) { if ($n != $i && $n != $j && $n != $k) { my ($xn,$yn,$zn); if (colinearity($mol,$j,$k,$n)) { # Try colinearity my $slope; my $xkj = $x[$k]-$x[$j]; my $ykj = $y[$k]-$y[$j]; my $zkj = $z[$k]-$z[$j]; if (abs($xkj)>=abs($ykj) && abs($xkj)>=abs($zkj)) { $slope = ($x[$n]-$x[$j])/$xkj; } elsif (abs($ykj)>=abs($zkj)) { $slope = ($y[$n]-$y[$j])/$ykj; } else { $slope = ($z[$n]-$z[$j])/$zkj; } $zn = 0; $yn = $slope*$yk; $xn = $slope*($xk-$xj)+$xj; } else { $xn = ($x[$n]*$x[$j] + $y[$n]*$y[$j] + $z[$n]*$z[$j])/$xj; $yn = ($x[$n]*$x[$k] + $y[$n]*$y[$k] + $z[$n]*$z[$k] - $xn*$xk) / $yk; my $sqr = $x[$n]**2 + $y[$n]**2 + $z[$n]**2 - $xn**2 - $yn**2; $sqr = 0 if $sqr < 0; $zn = sqrt($sqr); my ($A2,$B2,$C2) = plane_normal($mol,$j,$k,$n); my $Vec = ($B1*$C2-$C1*$B2)*($x[$k]-$x[$j])+ ($C1*$A2-$A1*$C2)*($y[$k]-$y[$j])+ ($A1*$B2-$B1*$A2)*($z[$k]-$z[$j]); $zn = -$zn if $Vec < 0; #print "$Vec\n"; } ($x[$n],$y[$n],$z[$n]) = ($xn,$yn,$zn); } } ($x[$j],$y[$j],$z[$j],$x[$k],$y[$k],$z[$k]) = ($xj,0,0,$xk,$yk,0); put_xyz($mol,\@x,\@y,\@z); #&write_molden; } sub colinearity { my ($mol,$i,$j,$k) = @_; my $eps = 1e-6; my $rij = sqrt(($mol->[$i][1]-$mol->[$j][1])**2+ ($mol->[$i][2]-$mol->[$j][2])**2+ ($mol->[$i][3]-$mol->[$j][3])**2); my $rik = sqrt(($mol->[$i][1]-$mol->[$k][1])**2+ ($mol->[$i][2]-$mol->[$k][2])**2+ ($mol->[$i][3]-$mol->[$k][3])**2); my $rjk = sqrt(($mol->[$k][1]-$mol->[$j][1])**2+ ($mol->[$k][2]-$mol->[$j][2])**2+ ($mol->[$k][3]-$mol->[$j][3])**2); my ($r1,$r2,$r3) = sort {$a<=>$b} ($rij,$rik,$rjk); if ($r1<$eps) { warn "Two atoms from $i,$j,$k coincide\n"; return 1; } if ($r1+$r2-$r3<$eps) { return 1; } return undef; } sub dist { my ($mol,$i,$j) = @_; my $rij = sqrt(($mol->[$i][1]-$mol->[$j][1])**2 + ($mol->[$i][2]-$mol->[$j][2])**2 + ($mol->[$i][3]-$mol->[$j][3])**2 ); return $rij; } sub rot_x { my ($mol,$degrees) = @_; my $rad = $degrees/180*3.14159265358979; my $cos = cos($rad); my $sin = sin($rad); for (my $i=1; $i<@$mol; $i++) { my ($y,$z) = ($mol->[$i][2],$mol->[$i][3]); $mol->[$i][2] = $sin*$z+$cos*$y; $mol->[$i][3] = $cos*$z-$sin*$y; } } sub copy_mol { my $mol = shift; my $new_mol; my $N = $#{$mol}; $new_mol->[0] = $mol->[0] if $mol->[0]; for (my $i=1; $i<=$N; $i++) { $new_mol->[$i] = [@{$mol->[$i]}]; } return $new_mol; }