http://qs321.pair.com?node_id=639921


This code is related to a question (Doing rigid body rotations in Perl) I asked some time ago:

# # example generate a perfect spiral (helix) of (10-Atom)-rods # use Math::Trig; my $axis = [ 0, 0, 1 ]; # z axis my @tmpl = ( [0,0,0,'HEAD'], [1,0,0,'SEG'], [2,0,0,'SEG'], [3,0,0,'SEG'], [4,0,0,'SEG'], [5,0,0,'SEG'], [6,0,0,'SEG'], [7,0,0,'SEG'], [8,0,0,'SEG'], [9,0,0,'HEAD'] ); my ($nMol, $nPeriods) = (360, 12); # 360 Molecules make 12 helix tur +ns my $increment = $nPeriods*360.*(1./$nMol); # put that into output # - - - - - - - - - - - - - - - - - - - - - - - - - - - - # my @ensemble; for my $id (0..$nMol-1) { my $w = deg2rad ($nPeriods * 360. * ($id/$nMol)); my $m = setrotmatrix ( $w, $axis ); my $p = deepcopy_body( \@tmpl ); center_body ( $p ); translate_body( $p => [0,0,$id] ); rotate_body ( $p => $m ); push @ensemble, $p; } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - # { # do another rotation of the whole thing my $mat = setrotmatrix( deg2rad(36), normalize([ 1, 1, 0 ]) ); for my $en (@ensemble) { rotate_body( $en, $mat ) } } my $fn = sprintf "ihx-b%02d-n%04d-p%d.kar", 0+@tmpl, $nMol, $nPeriods +; open my $fh, '>',$fn or die "can't open $!"; printf $fh "%d\t # %f/%f deg/rad increment\n", @tmpl * @ensemble, $increment, deg2rad($increment); my($molecule, $idmol, $boffs) = (0,1,1); # coordinates for $molecule (@ensemble) { for my $coord (@$molecule) { print $fh "@$coord $idmol 1\n" } $idmol++; } # bonds for $molecule (@ensemble) { for my $n (0..@$molecule-2) { print $fh $boffs+$n, ' ', $boffs+$n+ +1, "\n" } $boffs += @tmpl; } print $fh "0 0\n"; close $fh; # - - - - - - - - - - - - - - - - - - - - - - - - - - - - # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - # # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +- - - - # # sub deg2rad { pi * $_[0] / 180. } => already defined in Math::Trig # sub rad2deg { 180. * $_[0] / pi } => already defined in Math::Trig sub setrotmatrix { my ($theta, $axis) = @_; my ($x, $y, $z) = @$axis; my $c = cos $theta; my $s = sin $theta, my $t = 1. - $c; return [ [$t * $x * $x + $c, $t * $x * $y + $s * $z, $t * $x * $z - $s + * $y], [$t * $x * $y - $s * $z, $t * $y * $y + $c, $t * $y * $z + $s + * $x], [$t * $x * $z + $s * $y, $t * $y * $z - $s * $x, $t * $z * $z + $c + ] ] } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +- - - - # sub rotate { my ($v, $m) = @_; my ($x, $y, $z) = @$v; $v->[0] = $m->[0][0]*$x + $m->[0][1]*$y + $m->[0][2]*$z ; $v->[1] = $m->[1][0]*$x + $m->[1][1]*$y + $m->[1][2]*$z ; $v->[2] = $m->[2][0]*$x + $m->[2][1]*$y + $m->[2][2]*$z ; return $v; } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +- - - - # sub veclength2 { my ($v) = @_; my ($x, $y, $z) = @$v; return $x**2 + $y**2 + $z**2 } sub veclength { my ($v) = @_; return sqrt( veclength2($v) ) } sub normalize { my ($v) = @_; my $len = veclength($v); $len = ($len < 1e-10) ? 0. : 1./$len; $v->[0] *= $len, $v->[1] *= $len, $v->[2] *= $len; return $v; } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +- - - - # sub deepcopy_body { my ($body) = @_; my $copy = []; for my $p ( @$body ) { push @$copy, [ @$p ] } return $copy } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +- - - - # sub translate_body { my ($body, $disp) = @_; for my $p (@$body) { $p->[0] += $disp->[0], $p->[1] += $disp->[1], $p->[2] += $disp->[2 +]; } return $body } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +- - - - # sub center_body { my ($body) = @_; my $cnt = @$body; my ($pcoord, @com) = 0, (0,0,0); for $pcoord (@$body) { $com[0] += $pcoord->[0], $com[1] += $pcoord->[1], $com[2] += $pcoord->[2] } $com[0] /= $cnt, $com[1] /= $cnt, $com[2] /= $cnt; for $pcoord (@$body) { $pcoord->[0] -= $com[0], $pcoord->[1] -= $com[1], $pcoord->[2] -= $com[2] } return $body } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +- - - - # sub rotate_body { my ($body, $m) = @_; for my $v (@$body) { rotate($v, $m); } return $body } # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +- - - - # __END__