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__
|