Beefy Boxes and Bandwidth Generously Provided by pair Networks
good chemistry is complicated,
and a little bit messy -LW
 
PerlMonks  

mwah's scratchpad

by mwah (Hermit)
on Sep 19, 2007 at 16:10 UTC ( #639921=scratchpad: print w/ replies, xml ) Need Help??


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__



Log In?
Username:
Password:

What's my password?
Create A New User
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others romping around the Monastery: (15)
As of 2015-07-06 13:33 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    The top three priorities of my open tasks are (in descending order of likelihood to be worked on) ...









    Results (74 votes), past polls