| [a0bcf1] | 1 | #! @PERL@ -w
|
|---|
| 2 | #
|
|---|
| 3 | # This code helps in creating nanotube structures by specifying a unit cell and some
|
|---|
| 4 | # multiplicative factors which are used to first generate a rectangular sheet, then
|
|---|
| 5 | # rolling it up for the tube and again for the torus. All intermediate steps may be
|
|---|
| 6 | # supplied by file and then only later geometries are created.
|
|---|
| 7 | #
|
|---|
| 8 | # The hows and whys are explained in the comments before each section.
|
|---|
| 9 | #
|
|---|
| 10 | # \author Frederik Heber
|
|---|
| 11 | # \date 15.06.07
|
|---|
| 12 |
|
|---|
| 13 | use POSIX;
|
|---|
| 14 |
|
|---|
| 15 | our $pi = 3.1415926535; # no pi in perl
|
|---|
| 16 | our $MYEPSILON = 1e-13; # machine precision
|
|---|
| 17 |
|
|---|
| 18 | our @vectorA1;
|
|---|
| 19 | our @vectorA2;
|
|---|
| 20 | our @vectorA3;
|
|---|
| 21 | our @vector = ( \@vectorA1, \@vectorA2, \@vectorA3);
|
|---|
| 22 | our @betrag;
|
|---|
| 23 | our @RecivectorA1;
|
|---|
| 24 | our @RecivectorA2;
|
|---|
| 25 | our @RecivectorA3;
|
|---|
| 26 | our @Recivector = ( \@RecivectorA1, \@RecivectorA2, \@RecivectorA3);
|
|---|
| 27 | our @TubevectorA1;
|
|---|
| 28 | our @TubevectorA2;
|
|---|
| 29 | our @TubevectorA3;
|
|---|
| 30 | our @Tubevector = ( \@TubevectorA1, \@TubevectorA2, \@TubevectorA3);
|
|---|
| 31 | our @Tubebetrag = (0,0,0);
|
|---|
| 32 | our @TubevectorInverseA1;
|
|---|
| 33 | our @TubevectorInverseA2;
|
|---|
| 34 | our @TubevectorInverseA3;
|
|---|
| 35 | our @TubevectorInverse = ( \@TubevectorInverseA1, \@TubevectorInverseA2, \@TubevectorInverseA3);
|
|---|
| 36 | our @TubebetragInverse = (0,0,0);
|
|---|
| 37 | our @sheetnr;
|
|---|
| 38 | our @buffer;
|
|---|
| 39 | our @buffer2;
|
|---|
| 40 | our @buffer3;
|
|---|
| 41 |
|
|---|
| 42 | # ++++++++ Functions +++++++++++++++++++++++++++
|
|---|
| 43 |
|
|---|
| 44 | # Convertes units using the unit programme
|
|---|
| 45 | # \param $unit unit
|
|---|
| 46 | # \param $SrcUnit source unit name
|
|---|
| 47 | # \param $DstUnit destination unit name
|
|---|
| 48 | # \return converted unit
|
|---|
| 49 | sub UnitConversion
|
|---|
| 50 | {
|
|---|
| 51 | local ($unit, $SrcUnit, $DstUnit) = ($_[0], $_[1], $_[2]); # make local variables from given parameters
|
|---|
| 52 | $return = `units \'$unit, $SrcUnit,\' \'$DstUnit\' | head -n 1 | awk -F"* " {'print $2'}`;
|
|---|
| 53 | return $return; # return value
|
|---|
| 54 | }
|
|---|
| 55 |
|
|---|
| 56 | # Adds commentary stuff (needed for further levels) to xyz files
|
|---|
| 57 | # \param $filename file name
|
|---|
| 58 | # \param $atomicnumber number of atoms in cell
|
|---|
| 59 | sub AddAtomicNumber
|
|---|
| 60 | {
|
|---|
| 61 | local ($filename, $atomicnumber) = ($_[0], $_[1]);
|
|---|
| 62 | our @vector;
|
|---|
| 63 | our @Recivector;
|
|---|
| 64 | our @sheetnr;
|
|---|
| 65 | my $date = `date`;
|
|---|
| 66 | chomp($date);
|
|---|
| 67 | # fill buffer
|
|---|
| 68 | open(file1,$filename);
|
|---|
| 69 | @buffer = <file1>;
|
|---|
| 70 | close(file1);
|
|---|
| 71 |
|
|---|
| 72 | # open for writing and prepend
|
|---|
| 73 | open(file2,">$filename");
|
|---|
| 74 | print file2 "$atomicnumber\n"; # atomic number
|
|---|
| 75 | print file2 "\tgenerated with Nanotube creator on $date\n"; # ... and comment
|
|---|
| 76 | print file2 @buffer; # append buffer
|
|---|
| 77 |
|
|---|
| 78 | # Add primitive vectors as comment
|
|---|
| 79 | print file2 "\n****************************************\n\n";
|
|---|
| 80 | print file2 "Primitive vectors\n";
|
|---|
| 81 | print file2 "a(1) = $vector[0][0]\t$vector[0][1]\t $vector[0][2]\n";
|
|---|
| 82 | print file2 "a(2) = $vector[1][0]\t$vector[1][1]\t $vector[1][2]\n";
|
|---|
| 83 | print file2 "a(3) = $vector[2][0]\t$vector[2][1]\t $vector[2][2]\n";
|
|---|
| 84 | print file2 "\nVolume = ".$volume."\n";
|
|---|
| 85 | print file2 "\nReciprocal Vectors\n";
|
|---|
| 86 | print file2 "b(1) = $Recivector[0][0]\t$Recivector[0][1]\t$Recivector[0][2]\t\n";
|
|---|
| 87 | print file2 "b(2) = $Recivector[1][0]\t$Recivector[1][1]\t$Recivector[1][2]\t\n";
|
|---|
| 88 | print file2 "b(3) = $Recivector[2][0]\t$Recivector[2][1]\t$Recivector[2][2]\t\n";
|
|---|
| 89 |
|
|---|
| 90 | close(file2); # close file
|
|---|
| 91 | }
|
|---|
| 92 |
|
|---|
| 93 | # Transposes a 3x3 matrix
|
|---|
| 94 | # \param **matrixref
|
|---|
| 95 | sub Transpose
|
|---|
| 96 | {
|
|---|
| 97 | local ($matrixref) = @_;
|
|---|
| 98 | for (my $i=0;$i<3;$i++) {
|
|---|
| 99 | for (my $j=0;$j<$i;$j++) {
|
|---|
| 100 | my $tmp = ${${$matrixref}[$i]}[$j];
|
|---|
| 101 | ${${$matrixref}[$i]}[$j] = ${${$matrixref}[$j]}[$i];
|
|---|
| 102 | ${${$matrixref}[$j]}[$i] = $tmp;
|
|---|
| 103 | }
|
|---|
| 104 | }
|
|---|
| 105 | }
|
|---|
| 106 |
|
|---|
| 107 | # Adds some more commentary stuff to xyz files
|
|---|
| 108 | sub AddSheetInfo
|
|---|
| 109 | {
|
|---|
| 110 | local ($filename, $majoraxis, $minoraxis, $noaxis, $n, $m, $radius, $length) = @_;
|
|---|
| 111 | # open for writing and prepend
|
|---|
| 112 | open(file2,">>$filename");
|
|---|
| 113 |
|
|---|
| 114 | # Add primitive vectors as comment
|
|---|
| 115 | print file2 "\n****************************************\n\n";
|
|---|
| 116 | print file2 "Axis $majoraxis\t$minoraxis\t$noaxis\n";
|
|---|
| 117 | print file2 "(n,m) $n $m\n";
|
|---|
| 118 | print file2 "factors $radius $length\n";
|
|---|
| 119 | close(file2);
|
|---|
| 120 | }
|
|---|
| 121 |
|
|---|
| 122 | # Reads either value from stdin or recognizes if old value shall be used again by simple return key.
|
|---|
| 123 | # \param $_[0] old value
|
|---|
| 124 | # \return either old value or newly entered one
|
|---|
| 125 | sub GetValue
|
|---|
| 126 | {
|
|---|
| 127 | $input = <INPUT>;
|
|---|
| 128 | chomp $input;
|
|---|
| 129 | if ($input) {
|
|---|
| 130 | return $input;
|
|---|
| 131 | } else {
|
|---|
| 132 | print $_[0]."\n";
|
|---|
| 133 | return $_[0];
|
|---|
| 134 | }
|
|---|
| 135 | }
|
|---|
| 136 |
|
|---|
| 137 | # approximatively finds the greatest common divisor of two real numbers.
|
|---|
| 138 | # \param @array 3x3 matrix with row unit cell vectors
|
|---|
| 139 | # \param $axis1 major axis number
|
|---|
| 140 | # \param $axis2 minor axis number
|
|---|
| 141 | # \return approximative greatest common divisor
|
|---|
| 142 | # taken from http://www.perlmonks.org/?node_id=56906
|
|---|
| 143 | sub gcf {
|
|---|
| 144 | my ($x, $y) = @_;
|
|---|
| 145 | ($x, $y) = ($y, $x % $y) while $y;
|
|---|
| 146 | return $x;
|
|---|
| 147 | }
|
|---|
| 148 |
|
|---|
| 149 | # Evaluates scalar product
|
|---|
| 150 | # \param vector array 1
|
|---|
| 151 | # \param vector array 2
|
|---|
| 152 | # \return skp
|
|---|
| 153 | sub skp
|
|---|
| 154 | {
|
|---|
| 155 | my ($vec1ref, $vec2ref) = (@_);
|
|---|
| 156 | my $result = 0;
|
|---|
| 157 | for (my $i=0;$i<3;$i++) {
|
|---|
| 158 | #print ${$vec1ref}[$i]." * ".${$vec2ref}[$i]."\n";
|
|---|
| 159 | $result += ${$vec1ref}[$i] * ${$vec2ref}[$i];
|
|---|
| 160 | }
|
|---|
| 161 | return $result;
|
|---|
| 162 | }
|
|---|
| 163 |
|
|---|
| 164 | # Calculate norm with help of skp()
|
|---|
| 165 | sub norm
|
|---|
| 166 | {
|
|---|
| 167 | return sqrt(skp(@_));
|
|---|
| 168 | }
|
|---|
| 169 |
|
|---|
| 170 | # Evaluates projection of one vector onto another
|
|---|
| 171 | # \param $vec1ref reference to 3 vector upon which is projected
|
|---|
| 172 | # \param $vec2ref reference to 3 vector which is projected
|
|---|
| 173 | # \return projection factor
|
|---|
| 174 | sub projection
|
|---|
| 175 | {
|
|---|
| 176 | my ($vec1ref, $vec2ref) = (@_);
|
|---|
| 177 | return skp(@_)/skp($vec1ref, $vec1ref);
|
|---|
| 178 | }
|
|---|
| 179 |
|
|---|
| 180 | # Matrix transformation
|
|---|
| 181 | # \param $matrixref reference to 3x3 matrix A
|
|---|
| 182 | # \param $vectorref reference to 3 vector b
|
|---|
| 183 | # \return reference to resulting 3 vector Ab = x
|
|---|
| 184 | sub MatrixTrafo
|
|---|
| 185 | {
|
|---|
| 186 | my ($matrixref, $vectorref) = @_;
|
|---|
| 187 | @return = (0,0,0);
|
|---|
| 188 | for (my $i=0;$i<3;$i++) {
|
|---|
| 189 | for (my $j=0;$j<3;$j++) {
|
|---|
| 190 | #print ${${$matrixref}[$i]}[$j]." * ".${$vectorref}[$i]."\n";
|
|---|
| 191 | $return[$j] += ${${$matrixref}[$i]}[$j] * ${$vectorref}[$i];
|
|---|
| 192 | }
|
|---|
| 193 | }
|
|---|
| 194 | return @return;
|
|---|
| 195 | }
|
|---|
| 196 |
|
|---|
| 197 | # Fixed GramSchmidt-Orthogonalization for 3 vectors
|
|---|
| 198 | # \param @orthvector reference to 3x3 matrix
|
|---|
| 199 | # \param @orthbetrag reference to 3 vector with vector magnitudes
|
|---|
| 200 | # \param @axis major-, minor- and noaxis for specific order for the GramSchmidt procedure
|
|---|
| 201 | sub Orthogonalize
|
|---|
| 202 | {
|
|---|
| 203 | local ($orthvector, $orthbetrag, $axis) = @_;
|
|---|
| 204 | local $betrag1;
|
|---|
| 205 | local $betrag2;
|
|---|
| 206 |
|
|---|
| 207 | # orthogonalize axis[0] vector
|
|---|
| 208 | #print "Orthvectors: ${${$orthvector}[0]}[0]\t${${$orthvector}[0]}[1]\t${${$orthvector}[0]}[2]\n";
|
|---|
| 209 | #print "Orthbetrag: ${$orthbetrag}[0]\n";
|
|---|
| 210 | #print "OrthAxis: ${$axis}[0]\t${$axis}[1]\t${$axis}[2]\n";
|
|---|
| 211 | $betrag1 = projection(${$orthvector}[${$axis}[1]],${$orthvector}[${$axis}[0]]);
|
|---|
| 212 | for (my $k=0;$k<3;$k++) {
|
|---|
| 213 | ${${$orthvector}[$axis[0]]}[$k] -= ${${$orthvector}[$axis[1]]}[$k]*$betrag1;
|
|---|
| 214 | }
|
|---|
| 215 | # orthogonalize axis[2] vector
|
|---|
| 216 | $betrag1 = projection(${$orthvector}[${$axis}[0]],${$orthvector}[${$axis}[2]]);
|
|---|
| 217 | $betrag2 = projection(${$orthvector}[${$axis}[1]],${$orthvector}[${$axis}[2]]);
|
|---|
| 218 | for (my $k=0;$k<3;$k++) {
|
|---|
| 219 | ${${$orthvector}[${$axis}[2]]}[$k] -= ${${$orthvector}[${$axis}[0]]}[$k]*$betrag1 + ${${$orthvector}[${$axis}[1]]}[$k]*$betrag2;
|
|---|
| 220 | }
|
|---|
| 221 | }
|
|---|
| 222 |
|
|---|
| 223 | # calculate scalar product of vectors and store.
|
|---|
| 224 | # \param $vector reference to 3x3 matrix with row vectors
|
|---|
| 225 | # \param $betrag reference to 3 vector for found magnitudes
|
|---|
| 226 | sub DetermineVectorLengths
|
|---|
| 227 | {
|
|---|
| 228 | local ($vector, $betrag) = @_;
|
|---|
| 229 | for (my $i=0; $i<3; $i++) {
|
|---|
| 230 | ${$betrag}[$i] = skp(${$vector}[$i],${$vector}[$i]);
|
|---|
| 231 | print "\t$i: ${$betrag}[$i]";
|
|---|
| 232 | }
|
|---|
| 233 | }
|
|---|
| 234 |
|
|---|
| 235 | # Add two vectors
|
|---|
| 236 | # \param $vec1ref reference to first vector
|
|---|
| 237 | # \param $vec2ref reference to second vector
|
|---|
| 238 | # \return reference to vector sum
|
|---|
| 239 | sub VectorAdd
|
|---|
| 240 | {
|
|---|
| 241 | my ($vec1ref, $vec2ref) = @_;
|
|---|
| 242 | @return = (0,0,0);
|
|---|
| 243 | for (my $k=0;$k<3;$k++) {
|
|---|
| 244 | $return[$k] = ${$vec1ref}[$k] + ${$vec2ref}[$k];
|
|---|
| 245 | }
|
|---|
| 246 | return @return;
|
|---|
| 247 | }
|
|---|
| 248 |
|
|---|
| 249 | # Determines the inverse of a 3x3 matrix
|
|---|
| 250 | # \param $matrix1ref matrix to be inverted
|
|---|
| 251 | # \param $matrix2ref inverted matrix on exit
|
|---|
| 252 | # formula taken from Wikipedia (Regulaere Matrix)
|
|---|
| 253 | sub MatrixInversion
|
|---|
| 254 | {
|
|---|
| 255 | local ($matrix1ref,$matrix2ref) = @_;
|
|---|
| 256 | my $det = Det($matrix1ref);
|
|---|
| 257 | #print "Determinant is $det\n";
|
|---|
| 258 | ${${$matrix2ref}[0]}[0] = (${${$matrix1ref}[1]}[1]*${${$matrix1ref}[2]}[2] - ${${$matrix1ref}[1]}[2]*${${$matrix1ref}[2]}[1])/$det;
|
|---|
| 259 | ${${$matrix2ref}[1]}[0] = (${${$matrix1ref}[0]}[2]*${${$matrix1ref}[2]}[1] - ${${$matrix1ref}[0]}[1]*${${$matrix1ref}[2]}[2])/$det;
|
|---|
| 260 | ${${$matrix2ref}[2]}[0] = (${${$matrix1ref}[0]}[1]*${${$matrix1ref}[1]}[2] - ${${$matrix1ref}[0]}[2]*${${$matrix1ref}[1]}[1])/$det;
|
|---|
| 261 | ${${$matrix2ref}[0]}[1] = (${${$matrix1ref}[1]}[2]*${${$matrix1ref}[2]}[0] - ${${$matrix1ref}[1]}[0]*${${$matrix1ref}[2]}[2])/$det;
|
|---|
| 262 | ${${$matrix2ref}[1]}[1] = (${${$matrix1ref}[0]}[0]*${${$matrix1ref}[2]}[2] - ${${$matrix1ref}[0]}[2]*${${$matrix1ref}[2]}[0])/$det;
|
|---|
| 263 | ${${$matrix2ref}[2]}[1] = (${${$matrix1ref}[0]}[2]*${${$matrix1ref}[1]}[0] - ${${$matrix1ref}[0]}[0]*${${$matrix1ref}[1]}[2])/$det;
|
|---|
| 264 | ${${$matrix2ref}[0]}[2] = (${${$matrix1ref}[1]}[0]*${${$matrix1ref}[2]}[1] - ${${$matrix1ref}[1]}[1]*${${$matrix1ref}[2]}[0])/$det;
|
|---|
| 265 | ${${$matrix2ref}[1]}[2] = (${${$matrix1ref}[0]}[1]*${${$matrix1ref}[2]}[0] - ${${$matrix1ref}[0]}[0]*${${$matrix1ref}[2]}[1])/$det;
|
|---|
| 266 | ${${$matrix2ref}[2]}[2] = (${${$matrix1ref}[0]}[0]*${${$matrix1ref}[1]}[1] - ${${$matrix1ref}[0]}[1]*${${$matrix1ref}[1]}[0])/$det;
|
|---|
| 267 |
|
|---|
| 268 | print "Checking inverse ... ";
|
|---|
| 269 | my $flag = 0;
|
|---|
| 270 | for (my $i=0;$i<3;$i++) {
|
|---|
| 271 | for (my $j=0;$j<3;$j++) {
|
|---|
| 272 | my $tmp = 0;
|
|---|
| 273 | for (my $k=0;$k<3;$k++) {
|
|---|
| 274 | $tmp += ${${$matrix1ref}[$i]}[$k]*${${$matrix2ref}[$j]}[$k];
|
|---|
| 275 | }
|
|---|
| 276 | if ($flag == 0) {
|
|---|
| 277 | if ($i == $j) {
|
|---|
| 278 | $flag = fabs(1 - $tmp) > $MYEPSILON ? 1 : 0;
|
|---|
| 279 | } else {
|
|---|
| 280 | $flag = fabs($tmp) > $MYEPSILON ? 1 : 0;
|
|---|
| 281 | }
|
|---|
| 282 | }
|
|---|
| 283 | #print "$tmp";
|
|---|
| 284 | #if ($j != 2) { print ","; }
|
|---|
| 285 | #else { print "\n"; }
|
|---|
| 286 | }
|
|---|
| 287 | #if ($i != 2) { print ")\t("; }
|
|---|
| 288 | #else { print ")\n"; }
|
|---|
| 289 | }
|
|---|
| 290 | if ($flag == 0) { print "ok\n"; } else { print "false\n"; }
|
|---|
| 291 | }
|
|---|
| 292 |
|
|---|
| 293 | # Evaluates determinant of 3x3 matrix
|
|---|
| 294 | # \param $matrixref matrix whose determinant is to be calculated
|
|---|
| 295 | sub Det
|
|---|
| 296 | {
|
|---|
| 297 | local ($matrixref) = @_;
|
|---|
| 298 | my $return = 0;
|
|---|
| 299 | $return = ${${$matrixref}[0]}[0] * (${${$matrixref}[1]}[1]*${${$matrixref}[2]}[2] - ${${$matrixref}[1]}[2]*${${$matrixref}[2]}[1])
|
|---|
| 300 | - ${${$matrixref}[1]}[1] * (${${$matrixref}[0]}[0]*${${$matrixref}[2]}[2] - ${${$matrixref}[0]}[2]*${${$matrixref}[2]}[0])
|
|---|
| 301 | + ${${$matrixref}[2]}[2] * (${${$matrixref}[0]}[0]*${${$matrixref}[1]}[1] - ${${$matrixref}[0]}[1]*${${$matrixref}[1]}[0]);
|
|---|
| 302 | return $return;
|
|---|
| 303 | }
|
|---|
| 304 |
|
|---|
| 305 | # Parse for a pattern and return line split up by spaces
|
|---|
| 306 | sub ParseFor
|
|---|
| 307 | {
|
|---|
| 308 | my ($pattern, $arrayref) = @_;
|
|---|
| 309 | my $line;
|
|---|
| 310 | foreach $line (@{$arrayref}) {
|
|---|
| 311 | if ($line =~ /$pattern/) {
|
|---|
| 312 | return (split(/\s+/,$line));
|
|---|
| 313 | }
|
|---|
| 314 | }
|
|---|
| 315 | }
|
|---|
| 316 |
|
|---|
| 317 | # Reads contents of a file into an line-by-line array
|
|---|
| 318 | # \param filename
|
|---|
| 319 | # \return reference to array
|
|---|
| 320 | sub ReadFile
|
|---|
| 321 | {
|
|---|
| 322 | open(FILE,$_[0]);
|
|---|
| 323 | @storage = <FILE>;
|
|---|
| 324 | close(FILE);
|
|---|
| 325 | return @storage;
|
|---|
| 326 | }
|
|---|
| 327 |
|
|---|
| 328 | # A sheet, specified by axis[0] and axis[1] with vectors in $matrixref, is diametered
|
|---|
| 329 | # \param $matrixref reference to 3x3 matrix with row vectors
|
|---|
| 330 | # \param $axis reference to 3 vector with permutation of axis indices [0,1,2]
|
|---|
| 331 | # \param $lengthfactor factor for axis[0]
|
|---|
| 332 | # \param $radiusfactor factor for axis[1]
|
|---|
| 333 | # \return biggest diameter of sheet
|
|---|
| 334 | sub DetermineBiggestDiameter
|
|---|
| 335 | {
|
|---|
| 336 | local ($matrixref, $axis, $lengthfactor, $radiusfactor) = @_;
|
|---|
| 337 | # get greatest diameter of new sheet
|
|---|
| 338 | my @diameter = (0,0);
|
|---|
| 339 | for(my $i=0;$i<3;$i++) {
|
|---|
| 340 | $diameter[0] += (${${$matrixref}[${$axis}[0]]}[$i]*${$lengthfactor}-${${$matrixref}[${$axis}[1]]}[$i]*${$radiusfactor})
|
|---|
| 341 | *(${${$matrixref}[${$axis}[0]]}[$i]*${$lengthfactor}-${${$matrixref}[${$axis}[1]]}[$i]*${$radiusfactor});
|
|---|
| 342 | $diameter[1] += (${${$matrixref}[${$axis}[0]]}[$i]*${$lengthfactor}+${${$matrixref}[${$axis}[1]]}[$i]*${$radiusfactor})
|
|---|
| 343 | *(${${$matrixref}[${$axis}[0]]}[$i]*${$lengthfactor}+${${$matrixref}[${$axis}[1]]}[$i]*${$radiusfactor});
|
|---|
| 344 | }
|
|---|
| 345 | if (($diameter[0] - $diameter[1]) > $MYEPSILON) {
|
|---|
| 346 | $biggest = 0;
|
|---|
| 347 | } else {
|
|---|
| 348 | $biggest = 1;
|
|---|
| 349 | }
|
|---|
| 350 | @diameter = (sqrt($diameter[0]),sqrt($diameter[1]));
|
|---|
| 351 | printf("\n\nMajor diameter of the sheet is %5.5f, minor diameter is %5.5f.\n",$diameter[$biggest],$diameter[($biggest+1)%2]);
|
|---|
| 352 |
|
|---|
| 353 | return $diameter[$biggest];
|
|---|
| 354 | }
|
|---|
| 355 |
|
|---|
| 356 |
|
|---|
| 357 | # ++++++++ Beginning +++++++++++++++++++++++++++
|
|---|
| 358 |
|
|---|
| 359 | print "Nanotube Creator.\n=================\n\n"; # Welcome msg
|
|---|
| 360 |
|
|---|
| 361 | my $file;
|
|---|
| 362 |
|
|---|
| 363 | # Read command line arguments
|
|---|
| 364 | if (!$ARGV[1]) {
|
|---|
| 365 | print "Nanotubes <file> <status>\n\t<file> is the basis used for the various xyz files that are created.\n\t<status> is either None, Cell, Sheet or Tube,\n\t\t whether the given file does not yet exist, specifies a unit cell, a sheet or a tube.\n";
|
|---|
| 366 | exit 255;
|
|---|
| 367 | } else {
|
|---|
| 368 | $file = $ARGV[0];
|
|---|
| 369 | $stage = $ARGV[1];
|
|---|
| 370 | $file =~ s/(.*)\..*/$1/;
|
|---|
| 371 | print "I will use \'$file' as base for the filenames.\n\n";
|
|---|
| 372 |
|
|---|
| 373 | }
|
|---|
| 374 |
|
|---|
| 375 | # for text input
|
|---|
| 376 | open(INPUT,'-');
|
|---|
| 377 |
|
|---|
| 378 | # Get primitive vectors, either ...
|
|---|
| 379 | if ($stage !~ /[Nn]on/) {
|
|---|
| 380 | # ... from end of cell file
|
|---|
| 381 | print "Opening ".$file.".Cell.xyz to read primitive vectors in lines beginning with a(i) = ...\n";
|
|---|
| 382 | open(CELL, $file.".Cell.xyz");
|
|---|
| 383 | @buffer = <CELL>;
|
|---|
| 384 | local $flag = 3; # Marks whether vectors were found or not
|
|---|
| 385 | for($i=0; $i<@buffer; $i++) {
|
|---|
| 386 | #print "SEEK: $buffer[$i]";
|
|---|
| 387 | if ($buffer[$i] =~ /a\(1\) /) {
|
|---|
| 388 | ($name, $equal, @vectorA1) = split(/\s+/,$buffer[$i]);
|
|---|
| 389 | print "$name was found: ($vector[0][0], $vector[0][1], $vector[0][2])\n";
|
|---|
| 390 | $flag--;
|
|---|
| 391 | } elsif ($buffer[$i] =~ /a\(2\) /) {
|
|---|
| 392 | ($name, $equal, @vectorA2) = split(/\s+/,$buffer[$i]);
|
|---|
| 393 | print "$name was found: ($vector[1][0], $vector[1][1], $vector[1][2])\n";
|
|---|
| 394 | $flag--;
|
|---|
| 395 | } elsif ($buffer[$i] =~ /a\(3\) /) {
|
|---|
| 396 | ($name, $equal, @vectorA3) = split(/\s+/,$buffer[$i]);
|
|---|
| 397 | print "$name was found: ($vector[2][0], $vector[2][1], $vector[2][2])\n";
|
|---|
| 398 | $flag--;
|
|---|
| 399 | }
|
|---|
| 400 | }
|
|---|
| 401 | # parse number of atems in unit cell
|
|---|
| 402 | $numbercell = $buffer[0]; # number of atoms is first number in first line
|
|---|
| 403 | chomp $numbercell;
|
|---|
| 404 | print "\nThere are $numbercell atoms in the unit cell.\n";
|
|---|
| 405 |
|
|---|
| 406 | if ($flag != 0) { $stage = "None"; } # Set to None if vectors could not be found
|
|---|
| 407 | }
|
|---|
| 408 |
|
|---|
| 409 | if ($stage =~ /[Nn]on/) {
|
|---|
| 410 | # ... or from INPUT
|
|---|
| 411 | # give some explanation to help imagination
|
|---|
| 412 | print "You will give the unit cell of the given substance.\nAfterwards, the programme will create a Sheet, a Tube and a Torus, each with their own xyz-file named accordingly.\n\n";
|
|---|
| 413 |
|
|---|
| 414 | print "Enter 1st primitive vector (x z y): ";
|
|---|
| 415 | $line = <INPUT>;
|
|---|
| 416 | chomp($line);
|
|---|
| 417 | @vectorA1 = split(/\s+/,$line);
|
|---|
| 418 | print "Enter 2nd primitive vector (x z y): ";
|
|---|
| 419 | $line = <INPUT>;
|
|---|
| 420 | chomp($line);
|
|---|
| 421 | @vectorA2 = split(/\s+/,$line);
|
|---|
| 422 | print "Enter 3rd primitive vector (x z y): ";
|
|---|
| 423 | $line = <INPUT>;
|
|---|
| 424 | chomp($line);
|
|---|
| 425 | @vectorA3 = split(/\s+/,$line);
|
|---|
| 426 | }
|
|---|
| 427 |
|
|---|
| 428 | print "Cell vectors: (@{$vector[0]})\t(@{$vector[1]})\t(@{$vector[2]})\n";
|
|---|
| 429 |
|
|---|
| 430 | # Calculate some remaining stuff which comes at the end of the xyz file
|
|---|
| 431 | # volume is determinant of primitive vectors seen as a matrix (column vectors)
|
|---|
| 432 | $volume = Det(\@vector);
|
|---|
| 433 | print "Volume is $volume\n";
|
|---|
| 434 | MatrixInversion(\@vector, \@Recivector);
|
|---|
| 435 | print "Reciprocal cell vectors: (@{$Recivector[0]})\t(@{$Recivector[1]})\t(@{$Recivector[2]})\n";
|
|---|
| 436 | $Recivolume = Det(\@Recivector);
|
|---|
| 437 |
|
|---|
| 438 |
|
|---|
| 439 | # ============ CELL ===========================
|
|---|
| 440 | # The cell is simply created by transforming relative coordinates within the cell
|
|---|
| 441 | # into cartesian ones using the unit cell vectors.
|
|---|
| 442 |
|
|---|
| 443 | if ($stage =~ /[Nn]on/) {
|
|---|
| 444 | # enter number of atoms in unit cell
|
|---|
| 445 | print "\nHow many atoms are in the unit cell: ";
|
|---|
| 446 | $numbercell = <INPUT>;
|
|---|
| 447 | chomp($numbercell);
|
|---|
| 448 |
|
|---|
| 449 | # open unit cell file
|
|---|
| 450 | open(XYZ,">$file".".Cell.xyz"); # for xyz output
|
|---|
| 451 |
|
|---|
| 452 | print "\nNext, you have to enter each atom in the cell as follows, e.g.\n";
|
|---|
| 453 | print "Enter \'ChemicalSymbol X Y Z\' (relative to primitive vectors): C 0.5 0.25 0.5\n\n";
|
|---|
| 454 | # Enter element and coordinates of each
|
|---|
| 455 | for ($nr = 0; $nr < $numbercell; $nr++) {
|
|---|
| 456 | do {
|
|---|
| 457 | print "Enter for ".($nr+1)." atom \'ChemicalSymbol X Y Z\' (relative to primitive vectors): ";
|
|---|
| 458 | ($name, @atom) = split(/\s+/,<INPUT>);
|
|---|
| 459 | } while (@atom < 3);
|
|---|
| 460 | my @x = MatrixTrafo(\@vector, \@atom);
|
|---|
| 461 | print XYZ $name."\t".($x[0])."\t".($x[1])."\t".($x[2])."\n";
|
|---|
| 462 | }
|
|---|
| 463 |
|
|---|
| 464 | close(XYZ);
|
|---|
| 465 | AddAtomicNumber($file.".Cell.xyz", $numbercell); # prepend atomic number and comment
|
|---|
| 466 | print "\nThere are $numbercell atoms in the created cell.\n";
|
|---|
| 467 |
|
|---|
| 468 | # Read cell file into buffer for later geometries
|
|---|
| 469 | @buffer = ReadFile("<$file".".Cell.xyz");
|
|---|
| 470 |
|
|---|
| 471 | $stage = "Cell"
|
|---|
| 472 | }
|
|---|
| 473 |
|
|---|
| 474 | # evaluate length of each unit cell vector
|
|---|
| 475 | print "\nVektorbetraege: ";
|
|---|
| 476 | DetermineVectorLengths(\@vector,\@betrag);
|
|---|
| 477 |
|
|---|
| 478 | # ============ SHEET ==========================
|
|---|
| 479 | # The sheet is a bit more complex. We read the cell in cartesian coordinates
|
|---|
| 480 | # from the file. Next, we have to rotate the unit cell vectors by the so called
|
|---|
| 481 | # chiral angle. This gives a slanted and elongated section upon the sheet of
|
|---|
| 482 | # periodically repeated original unit cells. It only matches up if the factors
|
|---|
| 483 | # were all integer! (That's why the rotation is discrete and the chiral angle
|
|---|
| 484 | # specified not as (cos alpha, sin alpha) but as (n,m)) Also, we want this
|
|---|
| 485 | # section to be rectangular, thus we orthogonalize the original unit vectors
|
|---|
| 486 | # to gain our (later) tube axis.
|
|---|
| 487 | # By looking at the biggest possible diameter we know whose original cells to
|
|---|
| 488 | # look at and check if their respective compounds (contained atoms) still reside
|
|---|
| 489 | # in the rotated, elongated section we need for the later tube.
|
|---|
| 490 | # Then in a for loop we go through every cell. By projecting the vector leading
|
|---|
| 491 | # from the origin to the specific atom down onto the major and minor axis we
|
|---|
| 492 | # know if it's still within the boundaries spanned by these rotated and elongated
|
|---|
| 493 | # (radius-, length factor) unit vectors or not. If yes, its coordinates are
|
|---|
| 494 | # written to file.
|
|---|
| 495 |
|
|---|
| 496 | if ($stage =~ /[Cc]ell/) {
|
|---|
| 497 | # get radial axis of tube
|
|---|
| 498 | print "\nSpecify the axis permutation that is going to be perpendicular to the sheet [tubeaxis, torusaxis, noaxis]: ";
|
|---|
| 499 | $line = <INPUT>;
|
|---|
| 500 | chomp($line);
|
|---|
| 501 | if (!$line) {
|
|---|
| 502 | $majoraxis = 0;
|
|---|
| 503 | $minoraxis = 1;
|
|---|
| 504 | $noaxis = 2;
|
|---|
| 505 | } else {
|
|---|
| 506 | ($majoraxis, $minoraxis, $noaxis) = split(/\s+/,$line);
|
|---|
| 507 | }
|
|---|
| 508 |
|
|---|
| 509 | # find GCD for the two "longest" vector magnitudes
|
|---|
| 510 | #if ($betrag[($noaxis+1)%3] < $betrag[($noaxis+2)%3]) {
|
|---|
| 511 | # $majoraxis = ($noaxis+2)%3;
|
|---|
| 512 | # $minoraxis = ($noaxis+1)%3;
|
|---|
| 513 | #} else {
|
|---|
| 514 | # $majoraxis = ($noaxis+1)%3;
|
|---|
| 515 | # $minoraxis = ($noaxis+2)%3;
|
|---|
| 516 | #}
|
|---|
| 517 | #$gcf = ceil(1000*$betrag[$majoraxis])/gcf(ceil(1000*$betrag[$majoraxis]), ceil(1000*$betrag[$minoraxis])); # factor to a desired precision and round as they may not always match
|
|---|
| 518 | #print "Major axis is $majoraxis, minor axis is $minoraxis and greatest common factor is $gcf: ".$betrag[$majoraxis]." == ".($betrag[$minoraxis]*$gcf)."\n";
|
|---|
| 519 |
|
|---|
| 520 | my $flag = 0;
|
|---|
| 521 | do { # the discrete angle
|
|---|
| 522 | print "\nNow specify the two natural numbers (n m) defining the chiral angle: ";
|
|---|
| 523 | ($n,$m) = split(/\s+/,<INPUT>);
|
|---|
| 524 | chomp($m);
|
|---|
| 525 |
|
|---|
| 526 | # Compute rotated cell vectors
|
|---|
| 527 | for ($i=0; $i<3; $i++) {
|
|---|
| 528 | $Tubevector[$majoraxis][$i] = $n*$vector[$majoraxis][$i] + $m*$vector[$minoraxis][$i];
|
|---|
| 529 | $Tubevector[$minoraxis][$i] = -$m*$vector[$majoraxis][$i] + $n*$vector[$minoraxis][$i];
|
|---|
| 530 | $Tubevector[$noaxis][$i] = $vector[$noaxis][$i];
|
|---|
| 531 | }
|
|---|
| 532 | print "\nTubevektorbetraege: ";
|
|---|
| 533 | DetermineVectorLengths(\@Tubevector,\@Tubebetrag);
|
|---|
| 534 |
|
|---|
| 535 | # calculate angle and radius
|
|---|
| 536 | $tubecircum = sqrt($Tubebetrag[$minoraxis]);
|
|---|
| 537 | $tubelength = sqrt($Tubebetrag[$majoraxis]);
|
|---|
| 538 | $skp = skp($vector[$majoraxis],$Tubevector[$majoraxis]);
|
|---|
| 539 | $chiralangle = acos($skp/sqrt($betrag[$majoraxis]*$Tubebetrag[$majoraxis]))/$pi*180;
|
|---|
| 540 |
|
|---|
| 541 | print "\nGive integer factors for radius and length of tube: ";
|
|---|
| 542 | $line = <INPUT>;
|
|---|
| 543 | chomp($line);
|
|---|
| 544 | ($radiusfactor, $lengthfactor) = split(/\s+/, $line);
|
|---|
| 545 |
|
|---|
| 546 | printf("\nThe chiral angle then is %5.5f degrees with tube radius %5.5f A and length %5.5f A, i.e. torus radius of %5.5f A.\n", $chiralangle, $tubecircum*$radiusfactor/$pi, $tubelength*$lengthfactor, $tubelength*$lengthfactor/$pi);
|
|---|
| 547 |
|
|---|
| 548 | print "Satisfied? ";
|
|---|
| 549 | $flag = <INPUT>;
|
|---|
| 550 | } while ($flag !~ /[yY]/);
|
|---|
| 551 | } else {
|
|---|
| 552 | # Read sheetnumbers from file
|
|---|
| 553 | @buffer2 = ReadFile($file.".Sheet.xyz");
|
|---|
| 554 | ($pattern, $majoraxis, $minoraxis, $noaxis) = ParseFor("Axis", \@buffer2);
|
|---|
| 555 | ($pattern, $n, $m) = ParseFor("(n,m)", \@buffer2);
|
|---|
| 556 | ($pattern, $radiusfactor, $lengthfactor) = ParseFor("factors", \@buffer2);
|
|---|
| 557 |
|
|---|
| 558 | # Compute rotated cell vectors
|
|---|
| 559 | for ($i=0; $i<3; $i++) {
|
|---|
| 560 | $Tubevector[$majoraxis][$i] = $n*$vector[$majoraxis][$i] + $m*$vector[$minoraxis][$i];
|
|---|
| 561 | $Tubevector[$minoraxis][$i] = -$m*$vector[$majoraxis][$i] + $n*$vector[$minoraxis][$i];
|
|---|
| 562 | $Tubevector[$noaxis][$i] = $vector[$noaxis][$i];
|
|---|
| 563 | }
|
|---|
| 564 | print "\nTubevektorbetraege: ";
|
|---|
| 565 | DetermineVectorLengths(\@Tubevector,\@Tubebetrag);
|
|---|
| 566 |
|
|---|
| 567 | # calculate angle and radius
|
|---|
| 568 | $tubecircum = sqrt($Tubebetrag[$minoraxis]);
|
|---|
| 569 | $tubelength = sqrt($Tubebetrag[$majoraxis]);
|
|---|
| 570 | $skp = skp($vector[$majoraxis],$Tubevector[$majoraxis]);
|
|---|
| 571 | $chiralangle = acos($skp/sqrt($betrag[$majoraxis]*$Tubebetrag[$majoraxis]))/$pi*180;
|
|---|
| 572 | printf("\nThe chiral angle then is %5.5f degrees with tube radius %5.5f A and length %5.5f A, i.e. torus radius of %5.5f A.\n", $chiralangle, $tubecircum*$radiusfactor/$pi, $tubelength*$lengthfactor, $tubelength*$lengthfactor/$pi);
|
|---|
| 573 | }
|
|---|
| 574 |
|
|---|
| 575 | # Now for the rotated sheet we need orthogonal major and minor axis, but with major unchanged
|
|---|
| 576 | @axis = ($majoraxis,$minoraxis,$noaxis);
|
|---|
| 577 | Orthogonalize(\@Tubevector,\@Tubebetrag,\@axis);
|
|---|
| 578 | # print new Tubebetraege
|
|---|
| 579 | print "\nTubevektorbetraege after Orthogonalization: ";
|
|---|
| 580 | DetermineVectorLengths(\@Tubevector,\@Tubebetrag);
|
|---|
| 581 |
|
|---|
| 582 | print "Tube vectors: (@{$Tubevector[0]})\t(@{$Tubevector[1]})\t(@{$Tubevector[2]})\n";
|
|---|
| 583 | MatrixInversion(\@Tubevector, \@TubevectorInverse);
|
|---|
| 584 | Transpose(\@TubevectorInverse);
|
|---|
| 585 | Transpose(\@Recivector);
|
|---|
| 586 | print "Inverted tube vectors: (@{$TubevectorInverse[0]})\t(@{$TubevectorInverse[1]})\t(@{$TubevectorInverse[2]})\n";
|
|---|
| 587 |
|
|---|
| 588 | # check touched original(!) cells in each direction (in a circle with biggest sheet diameter as radius)
|
|---|
| 589 | $biggestdiameter = DetermineBiggestDiameter(\@Tubevector, \@axis, \$lengthfactor, \$radiusfactor);
|
|---|
| 590 | for ($i=0;$i<3;$i++) {
|
|---|
| 591 | $sheetnr[$i] = 0;
|
|---|
| 592 | }
|
|---|
| 593 | for ($i=0;$i<3;$i++) {
|
|---|
| 594 | for ($j=0;$j<3;$j++) {
|
|---|
| 595 | if (fabs($vector[$i][$j]) > $MYEPSILON) {
|
|---|
| 596 | $tmp = ceil($biggestdiameter/fabs($vector[$i][$j]));
|
|---|
| 597 | } else {
|
|---|
| 598 | $tmp = 0;
|
|---|
| 599 | }
|
|---|
| 600 | $sheetnr[$j] = $sheetnr[$j] > $tmp ? $sheetnr[$j] : $tmp;
|
|---|
| 601 | }
|
|---|
| 602 | }
|
|---|
| 603 | print "Maximum indices to regard: @sheetnr\n";
|
|---|
| 604 | for ($i=0;$i<3;$i++) {
|
|---|
| 605 | print "For axis $i: (".($vector[$i][0]*$sheetnr[$i])."\t".($vector[$i][1]*$sheetnr[$i])."\t".($vector[$i][2]*$sheetnr[$i]).") with ".sqrt($betrag[$i])."\n";
|
|---|
| 606 | }
|
|---|
| 607 |
|
|---|
| 608 | if ($stage =~ /[Cc]ell/) {
|
|---|
| 609 | # open sheet file
|
|---|
| 610 | open(XYZ,">$file".".Sheet.xyz"); # for xyz output
|
|---|
| 611 |
|
|---|
| 612 | # Quicker: scan all atoms beforehand
|
|---|
| 613 | for (my $nr = 0; $nr < $numbercell; $nr++) {
|
|---|
| 614 | push (@names, "none"); # These two are necessary as with ($name, @atom) = split(..); the same variable would be used and thus overwritten!
|
|---|
| 615 | push (@atoms, [0,0,0]); # They sort of allocate the necessary memory
|
|---|
| 616 | ($names[$nr],$atoms[$nr][0],$atoms[$nr][1],$atoms[$nr][2]) = split(/\s+/, $buffer[2+$nr]); # take line by line
|
|---|
| 617 | print "Prereading $names[$nr]_$nr\t@{$atoms[$nr]}\n";
|
|---|
| 618 | }
|
|---|
| 619 |
|
|---|
| 620 | # Now create the sheet
|
|---|
| 621 | $numbersheet = 0;
|
|---|
| 622 | @index = (0,0,0);
|
|---|
| 623 | @projection = (0,0,0);
|
|---|
| 624 | for ($index[$majoraxis] = -$sheetnr[$majoraxis]; $index[$majoraxis] < $sheetnr[$majoraxis]; $index[$majoraxis]++) { # NOTE: minor axis may start from 0! Check on this later ...
|
|---|
| 625 | for ($index[$minoraxis] = 0; $index[$minoraxis] < $sheetnr[$minoraxis]; $index[$minoraxis]++) { # These are all the cells that need be checked on
|
|---|
| 626 | #for ($index[$minoraxis] = -$sheetnr[$minoraxis]; $index[$minoraxis] <= $sheetnr[$minoraxis]; $index[$minoraxis]++) { # These are all the cells that need be checked on
|
|---|
| 627 | # Calculate offset
|
|---|
| 628 | @offset = MatrixTrafo(\@vector, \@index);
|
|---|
| 629 |
|
|---|
| 630 | # print "Now dealing with $numbercell atoms in unit cell at R = (@offset): ";
|
|---|
| 631 | my @coord;
|
|---|
| 632 | for (my $nr = 0; $nr < $numbercell; $nr++) {
|
|---|
| 633 | # Create coordinates at atom site
|
|---|
| 634 | @coord = VectorAdd($atoms[$nr], \@offset);
|
|---|
| 635 |
|
|---|
| 636 | # project down on major and minor Tubevectors and check for length if out of sheet
|
|---|
| 637 | #@projection = MatrixTrafo(\@Recivector, \@coord);
|
|---|
| 638 | #print "length: $lengthfactor\tradius: $radiusfactor\tdepth: 1\taxis: $majoraxis $minoraxis $noaxis\tEps: $MYEPSILON ";
|
|---|
| 639 | #if ($numbersheet < 2) {
|
|---|
| 640 | #if ((($projection[$majoraxis] + $MYEPSILON) >= 0) && ((1 - $projection[$majoraxis]) > $MYEPSILON) &&
|
|---|
| 641 | # (($projection[$minoraxis] + $MYEPSILON) >= 0) && ((1 - $projection[$minoraxis]) > $MYEPSILON) &&
|
|---|
| 642 | # (($projection[$noaxis] + $MYEPSILON) >= 0) && ((1 - $projection[$noaxis]) > $MYEPSILON)) { # check if within rotated cell
|
|---|
| 643 | # $char = 'C';
|
|---|
| 644 | #} else {
|
|---|
| 645 | # $char = 'O';
|
|---|
| 646 | #}
|
|---|
| 647 | @projection = MatrixTrafo(\@TubevectorInverse, \@coord);
|
|---|
| 648 | #print "projection: $projection[$majoraxis] $projection[$minoraxis] $projection[$noaxis]\n";
|
|---|
| 649 | #print "Atom Nr. ".($numbersheet+1).": @projection\n";
|
|---|
| 650 | if ((($projection[$majoraxis] + $MYEPSILON) >= 0) && (($lengthfactor - $projection[$majoraxis]) > $MYEPSILON) &&
|
|---|
| 651 | (($projection[$minoraxis] + $MYEPSILON) >= 0) && (($radiusfactor - $projection[$minoraxis]) > $MYEPSILON) &&
|
|---|
| 652 | (($projection[$noaxis] + $MYEPSILON) >= 0) && ((1 - $projection[$noaxis]) > $MYEPSILON)) { # check if within rotated cell
|
|---|
| 653 | $numbersheet++;
|
|---|
| 654 | print XYZ "$names[$nr]\t".($coord[0])."\t".($coord[1])."\t".($coord[2])."\n";
|
|---|
| 655 | #print XYZ "$char\t".($coord[0])."\t".($coord[1])."\t".($coord[2])."\n";
|
|---|
| 656 | #print "$nr,";
|
|---|
| 657 | } #else {
|
|---|
| 658 | #$numbersheet++;
|
|---|
| 659 | #print XYZ "$char\t".($coord[0])."\t".($coord[1])."\t".($coord[2])."\n";
|
|---|
| 660 | #print "#$nr#, ";
|
|---|
| 661 | #}
|
|---|
| 662 | }
|
|---|
| 663 | #print "\n";
|
|---|
| 664 | }
|
|---|
| 665 | }
|
|---|
| 666 |
|
|---|
| 667 | close(XYZ);
|
|---|
| 668 | AddAtomicNumber($file.".Sheet.xyz",$numbersheet); # prepend atomic number and comment
|
|---|
| 669 | AddSheetInfo($file.".Sheet.xyz",$majoraxis,$minoraxis,$noaxis,$n,$m,$radiusfactor,$lengthfactor);
|
|---|
| 670 | print "\nThere are $numbersheet atoms in the created sheet.\n";
|
|---|
| 671 |
|
|---|
| 672 | $stage = "Sheet";
|
|---|
| 673 | }
|
|---|
| 674 |
|
|---|
| 675 | open (file3, $file.".Sheet.xyz");
|
|---|
| 676 | @buffer2 = <file3>;
|
|---|
| 677 | close(file3);
|
|---|
| 678 | $numbersheet = $buffer2[0];
|
|---|
| 679 | chomp($numbersheet);
|
|---|
| 680 |
|
|---|
| 681 | # ============ TUBE ===========================
|
|---|
| 682 | # The tube starts with the rectangular (due to the orthogonalization) sheet
|
|---|
| 683 | # just created (or read). Along the minor axis it is rolled up, i.e. projected
|
|---|
| 684 | # from a 2d surface onto a cylindrical surface (x,y,z <-> r,alpha,z). The only
|
|---|
| 685 | # thing that's a bit complex is that the sheet it not aligned along the cartesian
|
|---|
| 686 | # axis but along major and minor. That's why we have to transform the atomic
|
|---|
| 687 | # cartesian coordinates into the orthogonal tubevector base, do the rolling up
|
|---|
| 688 | # there (and regard that minor and major axis must not necessarily be of equal
|
|---|
| 689 | # length) and afterwards transform back again (where we need the $halfaxis due to
|
|---|
| 690 | # the above possible inequality).
|
|---|
| 691 |
|
|---|
| 692 | if ($stage =~ /[Ss]heet/) {
|
|---|
| 693 | # open tube file
|
|---|
| 694 | open(XYZ,">$file".".Tube.xyz"); # for xyz output
|
|---|
| 695 |
|
|---|
| 696 | # Prepend begin
|
|---|
| 697 | print XYZ "$buffer2[0]"; # write number of atoms ...
|
|---|
| 698 | print XYZ "$buffer2[1]"; # ... and comment
|
|---|
| 699 |
|
|---|
| 700 | # determine center of gravity
|
|---|
| 701 | # @cog = (0,0,0);
|
|---|
| 702 | # for ($nr=0;$nr<$buffer2[0];$nr++) { # for each atom
|
|---|
| 703 | # ($name, @atom) = split(/\s+/, $buffer2[2+$nr]);
|
|---|
| 704 | # if (@atom >= 3) {
|
|---|
| 705 | # for (my $i=0;$i<3;$i++) {
|
|---|
| 706 | # $cog[$i] += $atom[$i];
|
|---|
| 707 | # }
|
|---|
| 708 | # }
|
|---|
| 709 | # }
|
|---|
| 710 | # for (my $i=0;$i<3;$i++) {
|
|---|
| 711 | # $cog[$i] /= -$buffer2[0];
|
|---|
| 712 | # }
|
|---|
| 713 | # @cog_projected = MatrixTrafo(\@TubevectorInverse, \@cog);
|
|---|
| 714 |
|
|---|
| 715 | # print "\nCenter of Gravity at (@cog) and projected at (@cog_projected)\n";
|
|---|
| 716 |
|
|---|
| 717 | # determine half axis as tube vector not necessarily half same length
|
|---|
| 718 | for (my $i=0;$i<3;$i++) {
|
|---|
| 719 | $halfaxis[$i] = $radiusfactor*($tubecircum)/sqrt($Tubebetrag[$i]);
|
|---|
| 720 | }
|
|---|
| 721 | # Rolling up the sheet
|
|---|
| 722 | my @x = (0,0,0);
|
|---|
| 723 | print "\nRadiusfactor: $radiusfactor\t pi: $pi\taxis ($majoraxis,$minoraxis,$noaxis)\n";
|
|---|
| 724 | for ($nr=0;$nr<$buffer2[0];$nr++) { # for each atom
|
|---|
| 725 | ($name, @atom) = split(/\s+/, $buffer2[2+$nr]);
|
|---|
| 726 | if (@atom >= 3) {
|
|---|
| 727 | #printf("Atom $nr: (%5.2f,%5.2f,%5.2f)\t",$atom[0],$atom[1],$atom[2]);
|
|---|
| 728 | # transform atom coordinates in cartesian system to the axis eigensystem
|
|---|
| 729 | @x = MatrixTrafo(\@TubevectorInverse, \@atom);
|
|---|
| 730 | #@x = VectorAdd(\@y, \@cog_projected);
|
|---|
| 731 | #printf("projected: (%5.2f,%5.2f,%5.2f)\t",$x[0],$x[1],$x[2]);
|
|---|
| 732 |
|
|---|
| 733 | # roll up (project (x,y,z) on cylindrical coordinates (radius,arg,z))
|
|---|
| 734 | $arg = 2*$pi*$x[$minoraxis]/($radiusfactor) - $pi; # is angle
|
|---|
| 735 | $radius = 1/(2*$pi); # is length of sheet in units of axis vector, divide by pi to get radius (from circumference)
|
|---|
| 736 | #printf("arg: %5.2f (c%2.2f,s%2.2f)\t",$arg, cos($arg), sin($arg));
|
|---|
| 737 | $x[$minoraxis] = cos($arg)*$halfaxis[$minoraxis]*$radius; # due to the back-transformation from eigensystem to cartesian one
|
|---|
| 738 | $x[$noaxis] = sin($arg)*$halfaxis[$noaxis]*$radius; # as both vectors are not normalized additional betrag has to be taken into account!
|
|---|
| 739 | #printf("rotated: (%5.2f,%5.2f,%5.2f)\n",$x[0],$x[1],$x[2]);
|
|---|
| 740 | @atom = MatrixTrafo(\@Tubevector, \@x);
|
|---|
| 741 | print XYZ $name."\t".$atom[0]."\t".$atom[1]."\t".$atom[2]."\n";
|
|---|
| 742 | #print $name."\t".$atom[0]."\t".$atom[1]."\t".$atom[2]."\n";
|
|---|
| 743 | } else {
|
|---|
| 744 | print "Error reading atom Nr. $nr\n";
|
|---|
| 745 | }
|
|---|
| 746 | }
|
|---|
| 747 | # append rest
|
|---|
| 748 | for ($i=$nr+2;$i<@buffer2;$i++) {
|
|---|
| 749 | print XYZ "$buffer2[$i]";
|
|---|
| 750 | }
|
|---|
| 751 | print XYZ "tubecircum $tubecircum\n";
|
|---|
| 752 | close(XYZ);
|
|---|
| 753 | print "\nThere are $numbersheet atoms in the created tube.\n";
|
|---|
| 754 |
|
|---|
| 755 | $stage = "Tube";
|
|---|
| 756 | }
|
|---|
| 757 |
|
|---|
| 758 | # Read tubecircum from file
|
|---|
| 759 | open (file3, "<$file".".Tube.xyz");
|
|---|
| 760 | @buffer3 = <file3>;
|
|---|
| 761 | close(file3);
|
|---|
| 762 | ($pattern, $tubecircum) = ParseFor("tubecircum", \@buffer3);
|
|---|
| 763 |
|
|---|
| 764 |
|
|---|
| 765 | # ============ TORUS ==========================
|
|---|
| 766 | # The procedure for the torus is very much alike to the one used to make the
|
|---|
| 767 | # tube. Only the projection is not from 2d surface onto a cylindrical one but
|
|---|
| 768 | # from a cylindrial onto a torus surface
|
|---|
| 769 | # (x,y,z) <-> (cos(s)*(R+r*cos(t)), sin(s)*(R+rcos(t)), r*sin(t)).
|
|---|
| 770 | # Here t is the angle within the tube with radius r, s is the torus angle with
|
|---|
| 771 | # radius R. We get R from the tubelength (that's why we need lengthfactor to
|
|---|
| 772 | # make it long enough). And due to fact that we have it already upon a cylindrical
|
|---|
| 773 | # surface, r*cos(t) and r*sin(t) already reside in $minoraxis and $noaxis.
|
|---|
| 774 |
|
|---|
| 775 | if ($stage =~ /[Tt]ube/) {
|
|---|
| 776 |
|
|---|
| 777 | # open torus file
|
|---|
| 778 | open(XYZ,">$file".".Torus.xyz"); # for xyz output
|
|---|
| 779 |
|
|---|
| 780 | # Prepend begin
|
|---|
| 781 | print XYZ "$buffer3[0]"; # write number of atoms ...
|
|---|
| 782 | print XYZ "$buffer3[1]"; # ... and comment
|
|---|
| 783 |
|
|---|
| 784 | # determine center of gravity
|
|---|
| 785 | # @cog = (0,0,0);
|
|---|
| 786 | # for ($nr=0;$nr<$buffer3[0];$nr++) { # for each atom
|
|---|
| 787 | # ($name, @atom) = split(/\s+/, $buffer3[2+$nr]);
|
|---|
| 788 | # if (@atom >= 3) {
|
|---|
| 789 | # for (my $i=0;$i<3;$i++) {
|
|---|
| 790 | # $cog[$i] += $atom[$i];
|
|---|
| 791 | # }
|
|---|
| 792 | # }
|
|---|
| 793 | # }
|
|---|
| 794 | # for (my $i=0;$i<3;$i++) {
|
|---|
| 795 | # $cog[$i] /= -$buffer3[0];
|
|---|
| 796 | # }
|
|---|
| 797 | # @ocog_projected = MatrixTrafo(\@TubevectorInverse, \@cog);
|
|---|
| 798 |
|
|---|
| 799 | # print "\nCenter of Gravity at (@cog) and projected at (@cog_projected)\n";
|
|---|
| 800 |
|
|---|
| 801 | # Determine half axis
|
|---|
| 802 | for (my $i=0;$i<3;$i++) {
|
|---|
| 803 | $halfaxis[$i] = ($tubelength)/sqrt($Tubebetrag[$i]);
|
|---|
| 804 | }
|
|---|
| 805 | # Rolling up the sheet
|
|---|
| 806 | my @x = (0,0,0);
|
|---|
| 807 | #print "\nRadiusfactor: $radiusfactor\t pi: $pi\taxis ($majoraxis,$minoraxis,$noaxis)\n";
|
|---|
| 808 | for ($nr=0;$nr<$buffer3[0];$nr++) { # for each atom
|
|---|
| 809 | ($name, @atom) = split(/\s+/, $buffer3[2+$nr]);
|
|---|
| 810 | if (@atom >= 3) {
|
|---|
| 811 | #printf("Atom $nr: (%5.2f,%5.2f,%5.2f)\t",$atom[0],$atom[1],$atom[2]);
|
|---|
| 812 |
|
|---|
| 813 | # transform atom coordinates in cartesian system to the axis eigensystem and shift into center of gravity
|
|---|
| 814 | @x = MatrixTrafo(\@TubevectorInverse, \@atom);
|
|---|
| 815 | #@x = VectorAdd(\@y, \@cog_projected);
|
|---|
| 816 | #printf("projected: (%5.2f,%5.2f,%5.2f)\t",$x[0],$x[1],$x[2]);
|
|---|
| 817 |
|
|---|
| 818 | # roll up (project (x,z,y) on cylindrical coordinates (radius,arg,z))
|
|---|
| 819 | $arg = 2*$pi*$x[$majoraxis]/($lengthfactor)-$pi; # is angle
|
|---|
| 820 | $radius = ($lengthfactor/(2*$pi) + $x[$minoraxis]/$halfaxis[$minoraxis]); # is R + r*cos(t)
|
|---|
| 821 | #printf("arg: %5.2f (c%2.2f,s%2.2f)\t",$arg, cos($arg), sin($arg));
|
|---|
| 822 | $x[$majoraxis] = cos($arg)*$halfaxis[$majoraxis]*($radius); # as both vectors are not normalized additional betrag has to be taken into account!
|
|---|
| 823 | $x[$minoraxis] = sin($arg)*$halfaxis[$minoraxis]*($radius); # due to the back-transformation from eigensystem to cartesian one
|
|---|
| 824 | #printf("rotated: (%5.2f,%5.2f,%5.2f)\n",$x[0],$x[1],$x[2]);
|
|---|
| 825 | @atom = MatrixTrafo(\@Tubevector, \@x);
|
|---|
| 826 |
|
|---|
| 827 | print XYZ $name."\t".$atom[0]."\t".$atom[1]."\t".$atom[2]."\n";
|
|---|
| 828 | #print $name."\t".$atom[0]."\t".$atom[1]."\t".$atom[2]."\n";
|
|---|
| 829 | } else {
|
|---|
| 830 | print "Error reading atom Nr. $nr\n";
|
|---|
| 831 | }
|
|---|
| 832 | }
|
|---|
| 833 |
|
|---|
| 834 | # append rest
|
|---|
| 835 | for ($i=$nr+2;$i<@buffer3;$i++) {
|
|---|
| 836 | print XYZ "$buffer3[$i]";
|
|---|
| 837 | }
|
|---|
| 838 | print XYZ "tubecircum $tubecircum\n";
|
|---|
| 839 | close(XYZ);
|
|---|
| 840 |
|
|---|
| 841 | $stage = "Torus";
|
|---|
| 842 | }
|
|---|
| 843 |
|
|---|
| 844 | print "\nFinished with geometry $stage.\n";
|
|---|
| 845 |
|
|---|
| 846 | # close stdin and exit
|
|---|
| 847 | close(INPUT);
|
|---|
| 848 | exit 0;
|
|---|
| 849 |
|
|---|
| 850 | # ++++++++ This is the end +++++++++++++++++++++
|
|---|
