File Coverage

blib/lib/Chemistry/Atom.pm
Criterion Covered Total %
statement 219 268 81.7
branch 57 90 63.3
condition 22 35 62.8
subroutine 45 50 90.0
pod 30 34 88.2
total 373 477 78.2


line stmt bran cond sub pod time code
1             package Chemistry::Atom;
2              
3             our $VERSION = '0.39'; # VERSION
4             # $Id$
5              
6             =head1 NAME
7              
8             Chemistry::Atom - Chemical atoms as objects in molecules
9              
10             =head1 SYNOPSIS
11              
12             use Chemistry::Atom;
13              
14             my $atom = new Chemistry::Atom(
15             id => 'a1',
16             coords => [$x, $y, $z],
17             symbol => 'Br'
18             );
19              
20             print $atom->print;
21              
22             =head1 DESCRIPTION
23              
24             This module includes objects to describe chemical atoms.
25             An atom is defined by its symbol and its coordinates, among other attributes.
26             Atomic coordinates are described by a Math::VectorReal
27             object, so that they can be easily used in vector operations.
28              
29             =head2 Atom Attributes
30              
31             In addition to common attributes such as id, name, and type,
32             atoms have the following attributes, which are accessed or
33             modified through methods defined below: bonds, coords, internal_coords,
34             Z, symbol, etc.
35              
36             In general, to get the value of a property, use $atom->method without
37             any parameters. To set the value, use $atom->method($new_value). When setting
38             an attribute, the accessor returns the atom object, so that accessors can be
39             chained:
40              
41             $atom->symbol("C")->name("CA")->coords(1,2,3);
42              
43             =cut
44              
45             # Considering to add the following attributes:
46             # mass_number (A)
47              
48 16     16   57659 use 5.006;
  16         61  
49 16     16   65 use strict;
  16         25  
  16         248  
50 16     16   54 use warnings;
  16         29  
  16         379  
51 16     16   67 use Scalar::Util 'weaken';
  16         28  
  16         970  
52 16     16   7225 use Math::VectorReal qw(O vector);
  16         55365  
  16         1049  
53 16     16   8085 use Math::Trig;
  16         186111  
  16         2493  
54 16     16   144 use Carp;
  16         29  
  16         764  
55 16     16   89 use base qw(Chemistry::Obj Exporter);
  16         26  
  16         4860  
56 16     16   97 use List::Util qw(first);
  16         30  
  16         17828  
57              
58             our @EXPORT_OK = qw(distance angle dihedral angle_deg dihedral_deg);
59             our %EXPORT_TAGS = (
60             all => \@EXPORT_OK,
61             );
62              
63              
64              
65             my $N = 0; # Number of atoms created so far, used to generate default IDs.
66              
67             our @ELEMENTS = qw(
68             n
69             H He
70             Li Be B C N O F Ne
71             Na Mg Al Si P S Cl Ar
72             K Ca Sc Ti V Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
73             Rb Sr Y Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I Xe
74             Cs Ba
75             La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb
76             Lu Hf Ta W Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn
77             Fr Ra
78             Ac Th Pa U Np Pu Am Cm Bk Cf Es Fm Md No
79             Lr Rf Db Sg Bh Hs Mt Ds Uuu Uub Uut Uuq Uup Uuh Uus Uuo
80             );
81              
82             our %ELEMENTS;
83             for (my $i = 1; $i < @ELEMENTS; ++$i){
84             $ELEMENTS{$ELEMENTS[$i]} = $i;
85             }
86             $ELEMENTS{D} = $ELEMENTS{T} = 1;
87              
88             my %Atomic_masses = (
89             "H" => 1.00794, "D" => 2.014101, "T" => 3.016049, "He" => 4.002602,
90             "Li" => 6.941, "Be" => 9.012182, "B" => 10.811, "C" => 12.0107,
91             "N" => 14.00674, "O" => 15.9994, "F" => 18.9984032, "Ne" => 20.1797,
92             "Na" => 22.989770, "Mg" => 24.3050, "Al" => 26.981538, "Si" => 28.0855,
93             "P" => 30.973761, "S" => 32.066, "Cl" => 35.4527, "Ar" => 39.948,
94             "K" => 39.0983, "Ca" => 40.078, "Sc" => 44.955910, "Ti" => 47.867,
95             "V" => 50.9415, "Cr" => 51.9961, "Mn" => 54.938049, "Fe" => 55.845,
96             "Co" => 58.933200, "Ni" => 58.6934, "Cu" => 63.546, "Zn" => 65.39,
97             "Ga" => 69.723, "Ge" => 72.61, "As" => 74.92160, "Se" => 78.96,
98             "Br" => 79.904, "Kr" => 83.80, "Rb" => 85.4678, "Sr" => 87.62,
99             "Y" => 88.90585, "Zr" => 91.224, "Nb" => 92.90638, "Mo" => 95.94,
100             "Tc" => 98, "Ru" => 101.07, "Rh" => 102.90550, "Pd" => 106.42,
101             "Ag" => 107.8682, "Cd" => 112.411, "In" => 114.818, "Sn" => 118.710,
102             "Sb" => 121.760, "Te" => 127.60, "I" => 126.90447, "Xe" => 131.29,
103             "Cs" => 132.90545, "Ba" => 137.327, "La" => 138.9055, "Ce" => 140.116,
104             "Pr" => 140.90765, "Nd" => 144.24, "Pm" => 145, "Sm" => 150.36,
105             "Eu" => 151.964, "Gd" => 157.25, "Tb" => 158.92534, "Dy" => 162.50,
106             "Ho" => 164.93032, "Er" => 167.26, "Tm" => 168.93421, "Yb" => 173.04,
107             "Lu" => 174.967, "Hf" => 178.49, "Ta" => 180.9479, "W" => 183.84,
108             "Re" => 186.207, "Os" => 190.23, "Ir" => 192.217, "Pt" => 195.078,
109             "Au" => 196.96655, "Hg" => 200.59, "Tl" => 204.3833, "Pb" => 207.2,
110             "Bi" => 208.98038, "Po" => 209, "At" => 210, "Rn" => 222,
111             "Fr" => 223, "Ra" => 226, "Ac" => 227, "Th" => 232.038,
112             "Pa" => 231.03588, "U" => 238.0289, "Np" => 237, "Pu" => 244,
113             "Am" => 243, "Cm" => 247, "Bk" => 247, "Cf" => 251,
114             "Es" => 252, "Fm" => 257, "Md" => 258, "No" => 259,
115             "Lr" => 262, "Rf" => 261, "Db" => 262, "Sg" => 266,
116             "Bh" => 264, "Hs" => 269, "Mt" => 268, "Ds" => 271,
117             );
118              
119             =head1 METHODS
120              
121             =over 4
122              
123             =item Chemistry::Atom->new(name => value, ...)
124              
125             Create a new Atom object with the specified attributes.
126              
127             =cut
128              
129             sub new {
130 536     536 1 3519 my $class = shift;
131 536         900 my %args = @_;
132              
133 536         778 my $self = bless {
134             id => $class->nextID(),
135             coords => vector(0, 0, 0),
136             Z => 0,
137             symbol => '',
138             bonds => [],
139             }, $class;
140              
141 536         6791 $self->$_($args{$_}) for (keys %args);
142 536         1135 $self;
143             }
144              
145             sub nextID {
146 538     538 0 1216 "a".++$N;
147             }
148              
149             sub reset_id {
150 0     0 0 0 $N = 0;
151             }
152              
153              
154             =item $atom->Z($new_Z)
155              
156             Sets and returns the atomic number (Z). If the symbol of the atom doesn't
157             correspond to a known element, Z = undef.
158              
159             =cut
160              
161             sub Z {
162 7     7 1 13 my $self = shift;
163              
164 7 100       20 if(@_) {
165 4         75 $self->{symbol} = $ELEMENTS[$_[0]];
166 4         7 $self->{Z} = $_[0];
167 4         13 return $self;
168             } else {
169 3         11 return $self->{Z};
170             }
171             }
172              
173              
174             =item $atom->symbol($new_symbol)
175              
176             Sets and returns the atomic symbol.
177              
178             =cut
179              
180             sub symbol {
181 1091     1091 1 1199 my $self = shift;
182              
183 1091 100       1478 if(@_) {
184 531         573 my $symbol = shift;
185 531         729 $symbol =~ s/ //g;
186 531         874 $self->{Z} = $ELEMENTS{$symbol};
187 531         630 $self->{symbol} = $symbol;
188 531         762 return $self;
189             } else {
190 560         974 return $self->{symbol};
191             }
192             }
193              
194             =item $atom->mass($new_mass)
195              
196             Sets and returns the atomic mass in atomic mass units. Any arbitrary mass may
197             be set explicitly by using this method. However, if no mass is set explicitly
198             and this method is called as an accessor, the return value is the following:
199              
200             1) If the mass number is undefined (see the mass_number method below), the
201             relative atomic mass from the 1995 IUPAC recommendation is used. (Table stolen
202             from the Chemistry::MolecularMass module by Maksim A. Khrapov).
203              
204             2) If the mass number is defined and the L module is
205             available and it knows the mass for the isotope, the exact mass of the isotope
206             is used; otherwise, the mass number is returned.
207              
208             =cut
209              
210             sub mass {
211 13     13 1 41 my $self = shift;
212 13 100       30 if (@_) {
213 2         5 $self->{mass} = shift;
214 2         4 return $self;
215             } else {
216 11 100       32 if (defined $self->{mass}) {
    100          
217 2         7 return $self->{mass};
218             } elsif (defined $self->{mass_number}) {
219 2 100 66     5 if (eval { require Chemistry::Isotope } and
  2         1049  
220             my $m = Chemistry::Isotope::isotope_mass(
221             $self->{mass_number}, $self->{Z})
222             ) {
223 1         16022 return $m;
224             } else {
225 1         13 return $self->{mass_number};
226             }
227             } else {
228 7         14 return $Atomic_masses{$self->symbol};
229             }
230             }
231             }
232              
233             =item $atom->mass_number($new_mass_number)
234              
235             Sets or gets the mass number. The mass number is undefined unless is
236             set explicitly (this module does not try to guess a default mass number based
237             on the natural occurring isotope distribution).
238              
239             =cut
240              
241             Chemistry::Obj::accessor('mass_number');
242              
243             =item $atom->coords
244              
245             my $vector = $atom->coords; # get a Math::VectorReal object
246             $atom->coords($vector); # set a Math::VectorReal object
247             $atom->coords([$x, $y, $z]); # also accepts array refs
248             $atom->coords($x, $y, $z); # also accepts lists
249              
250             Sets or gets the atom's coordinates. It can take as a parameter a
251             Math::VectorReal object, a reference to an array, or the list of coordinates.
252              
253             =cut
254              
255             sub coords {
256 27     27 1 485 my $self = shift;
257              
258 27 100       62 if(@_) {
259 6 100       29 if (UNIVERSAL::isa($_[0], "Math::VectorReal")) {
    100          
260 2         7 $self->{coords} = $_[0];
261             } elsif (ref $_[0] eq "ARRAY") {
262 3         5 $self->{coords} = vector(@{$_[0]});
  3         10  
263             } else {
264 1         6 $self->{coords} = vector(@_);
265             }
266             } else {
267 21         106 return $self->{coords};
268             }
269 6         93 $self;
270             }
271              
272             =item $atom->internal_coords
273              
274             # get a Chemistry::InternalCoords object
275             my $ic = $atom->internal_coords;
276              
277             # set a Chemistry::InternalCoords object
278             $atom->internal_coords($vic);
279              
280             # also accepts array refs
281             $atom->internal_coords([8, 1.54, 7, 109.47, 6, 120.0]);
282              
283             # also accepts lists
284             $atom->internal_coords(8, 1.54, 7, 109.47, 6, 120.0);
285              
286             Sets or gets the atom's internal coordinates. It can take as a parameter a
287             Chemistry::InternalCoords object, a reference to an array, or the list of
288             coordinates. In the last two cases, the list elements are the following: atom
289             number or reference for distance, distance, atom number or reference for angle,
290             angle in degrees, atom number or reference for dihedral, dihedral in degrees.
291              
292             =cut
293              
294             sub internal_coords {
295 0     0 1 0 my $self = shift;
296              
297 0 0       0 if(@_) {
298 0 0       0 if (UNIVERSAL::isa($_[0], "Chemistry::InternalCoords")) {
    0          
299 0         0 $self->{internal_coords} = $_[0];
300             } elsif (ref $_[0] eq "ARRAY") {
301 0         0 require Chemistry::InternalCoords;
302             $self->{internal_coords} =
303 0         0 Chemistry::InternalCoords->new($self, @{$_[0]});
  0         0  
304             } else {
305 0         0 require Chemistry::InternalCoords;
306             $self->{internal_coords} =
307 0         0 Chemistry::InternalCoords->new($self, @_);
308             }
309             } else {
310 0         0 return $self->{internal_coords};
311             }
312 0         0 $self;
313             }
314              
315             =item $atom->x3, $atom->y3, $atom->z3
316              
317             Get the x, y or z 3D coordinate of the atom. This methods are just accessors
318             that don't change the coordinates. $atom->x3 is short for
319             ($atom->coords->array)[0], and so on.
320              
321             =cut
322              
323 2     2 1 8 sub x3 { ($_[0]->coords->array)[0] }
324 2     2 1 589 sub y3 { ($_[0]->coords->array)[1] }
325 2     2 1 610 sub z3 { ($_[0]->coords->array)[2] }
326              
327             =item $atom->formal_charge($charge)
328              
329             Set or get the formal charge of the atom.
330              
331             =cut
332              
333             Chemistry::Obj::accessor('formal_charge');
334              
335             =item $atom->formal_radical($radical)
336              
337             Set or get the formal radical multiplicity for the atom.
338              
339             =cut
340              
341             Chemistry::Obj::accessor('formal_radical');
342              
343             =item $atom->implicit_hydrogens($h_count)
344              
345             Set or get the number of implicit hydrogen atoms bonded to the atom.
346              
347             =cut
348              
349 21     21 1 76 sub implicit_hydrogens { shift->hydrogens(@_) }
350              
351             =item $atom->hydrogens($h_count)
352              
353             Set or get the number of implicit hydrogen atoms bonded to the atom
354             (DEPRECATED: USE C INSTEAD).
355              
356             =cut
357              
358             Chemistry::Obj::accessor('hydrogens');
359              
360             =item $atom->total_hydrogens($h_count)
361              
362             Get the total number of hydrogen atoms bonded to the atom (implicit +
363             explicit).
364              
365             =cut
366              
367             sub total_hydrogens {
368 5     5 1 12 my ($self) = @_;
369 16     16   129 no warnings 'uninitialized';
  16         30  
  16         2731  
370 5         12 $self->hydrogens + grep { $_->symbol eq 'H' } $self->neighbors;
  2         4  
371             }
372              
373             =item $atom->sprout_hydrogens
374              
375             Convert all the implicit hydrogens for this atom to explicit hydrogens. Note:
376             it does B generate coordinates for the new atoms.
377              
378             =cut
379              
380             sub sprout_hydrogens {
381 2     2 1 8 my ($self) = @_;
382 2   50     5 for (1 .. $self->implicit_hydrogens || 0) {
383 4         10 $self->parent->new_bond(
384             atoms => [$self, $self->parent->new_atom(symbol => 'H')]);
385             }
386 2         6 $self->implicit_hydrogens(0);
387             }
388              
389             =item $atom->collapse_hydrogens
390              
391             Delete neighboring hydrogen atoms and add them as implicit hydrogens for this
392             atom.
393              
394             =cut
395              
396             sub collapse_hydrogens {
397 2     2 1 6 my ($self) = @_;
398 16     16   98 no warnings 'uninitialized';
  16         33  
  16         2335  
399 2         5 my $implicit = $self->implicit_hydrogens;
400 2         9 for my $nei ($self->neighbors) {
401 4 50       10 $nei->delete, $implicit++ if $nei->symbol eq 'H';
402             }
403 2         9 $self->implicit_hydrogens($implicit);
404             }
405              
406             my %VALENCE_TABLE = (
407             Br => 1, Cl => 1, B => 3, C => 4, N => 3, O => 2, P => 3, S => 2,
408             F => 1, I => 1,
409             );
410              
411             # to make it easier to test
412             sub _calc_implicit_hydrogens {
413 21     21   8384 my ($self, $symbol, $valence, $charge, $radical) = @_;
414 16     16   111 no warnings 'uninitialized';
  16         33  
  16         23010  
415              
416 21         34 my $h_count = $VALENCE_TABLE{$symbol} - $valence;
417             # should account for non-kekulized aromatic bonds
418              
419             # some common charge situations
420 21 100 100     174 if (($symbol =~ /^(?:[NOSFI]|Cl|Br)$/) && $charge == -1) {
    100 100        
    100 100        
    100 66        
421 4         6 $h_count--;
422             } elsif ($symbol =~ /^[NOSP]$/ && $charge == 1) {
423 4         7 $h_count++;
424             } elsif ($symbol eq 'C' && $charge) {
425 2         6 $h_count--;
426             } elsif ($symbol eq 'B' && $charge == -1) {
427 1         14 $h_count++;
428             }
429              
430             # some common radical situations
431 21 100 100     72 if ($radical == 1 or $radical == 3) {
    100          
432             # carbene, singlet or triplet
433 3         9 $h_count -= 2;
434             } elsif ($radical == 2) {
435             # radical (doublet)
436 2         2 $h_count--;
437             }
438              
439 21 50       36 $h_count = 0 if $h_count < 0;
440 21         34 $h_count;
441             }
442              
443             =item $atom->calc_implicit_hydrogens
444              
445             Use heuristics to figure out how many implicit hydrogens should the atom have
446             to satisfy its normal "organic" valence. Returns the number of hydrogens but
447             does not affect the atom object.
448              
449             =cut
450              
451             sub calc_implicit_hydrogens {
452 3     3 1 4 my ($self) = @_;
453 3         6 $self->_calc_implicit_hydrogens(
454             $self->symbol, $self->explicit_valence,
455             $self->formal_charge, $self->formal_radical,
456             );
457             }
458              
459             =item $atom->add_implicit_hydrogens
460              
461             Similar to calc_implicit_hydrogens, but it also sets the number of implicit
462             hydrogens in the atom to the new calculated value. Equivalent to
463              
464             $atom->implicit_hydrogens($atom->calc_implicit_hydrogens);
465              
466             It returns the atom object.
467              
468             =cut
469              
470             sub add_implicit_hydrogens {
471 3     3 1 5 my ($self) = @_;
472 3         7 my $h_count = $self->calc_implicit_hydrogens;
473 3         6 $self->implicit_hydrogens($h_count);
474             }
475              
476             =item $atom->aromatic($bool)
477              
478             Set or get whether the atom is considered to be aromatic. This property may be
479             set arbitrarily, it doesn't imply any kind of "intelligent aromaticity
480             detection"! (For that, look at the L module).
481              
482             =cut
483              
484             Chemistry::Obj::accessor('aromatic');
485              
486             =item $atom->valence
487              
488             Returns the sum of the bond orders of the bonds in which the atom participates,
489             including implicit hydrogens (which are assumed to have bond orders of one).
490              
491             =cut
492              
493             sub valence {
494 3     3 1 7 my ($self) = @_;
495 3         6 my $valence = 0;
496 3         6 $valence += $_->order for $self->bonds;
497 3   100     8 $valence += $self->hydrogens || 0;
498 3         12 $valence;
499             }
500              
501             =item $atom->explicit_valence
502              
503             Like C, but excluding implicit hydrogen atoms. To get the raw number
504             of bonds, without counting bond orders, call $atom->bonds in scalar context.
505              
506             =cut
507              
508             sub explicit_valence {
509 6     6 1 18 my ($self) = @_;
510 6         10 my $valence = 0;
511 6         14 $valence += $_->order for $self->bonds;
512 6         15 $valence;
513             }
514              
515             # this method is for internal use only; called by $mol->add_bond
516             sub add_bond {
517 20     20 0 32 my $self = shift;
518 20         23 my $bond = shift;
519 20         25 my %seen;
520             #return if first { $_ eq $bond } @{$self->{bonds}};
521              
522 20         22 for my $atom (@{$bond->{atoms}}){ #for each atom...
  20         47  
523 40 100       86 if ($atom ne $self) {
524 20         85 my $b = {to=>$atom, bond=>$bond};
525 20         65 weaken($b->{to});
526 20         41 weaken($b->{bond});
527 20         24 push @{$self->{bonds}}, $b;
  20         56  
528             }
529             }
530             }
531              
532             # make sure the atom doesn't cause circular references
533             sub _weaken {
534 112     112   131 my $self = shift;
535 112         122 for my $b (@{$self->{bonds}}) {
  112         209  
536 168         328 weaken($b->{to});
537 168         272 weaken($b->{bond});
538             }
539 112         238 weaken($self->{parent});
540             }
541              
542             # This method is private. Bonds should be deleted from the
543             # mol object. These methods should only be called by
544             # $bond->delete_atoms, which is called by $mol->delete_bond
545             sub delete_bond {
546 28     28 0 42 my ($self, $bond) = @_;
547 28         36 $self->{bonds} = [ grep { $_->{bond} ne $bond } @{$self->{bonds}} ];
  51         88  
  28         52  
548             }
549              
550             =item $atom->delete
551              
552             Calls $mol->delete_atom($atom) on the atom's parent molecule.
553              
554             =cut
555              
556             sub delete {
557 7     7 1 13 my ($self) = @_;
558 7         36 $self->{parent}->_delete_atom($self);
559             }
560              
561             =item $atom->parent
562              
563             Returns the atom's containing object (the molecule to which the atom belongs).
564             An atom can only have one parent.
565              
566             =cut
567              
568             sub parent {
569 561     561 1 630 my $self = shift;
570 561 100       722 if (@_) {
571 553         741 ($self->{parent}) = @_;
572 553         1068 weaken($self->{parent});
573 553         841 $self;
574             } else {
575 8         19 $self->{parent};
576             }
577             }
578              
579             =item $atom->neighbors($from)
580              
581             Return a list of neighbors. If an atom object $from is specified, it will be
582             excluded from the list (this is useful if an atom wants to know who its
583             neighbor's neighbors are, without counting itself).
584              
585             =cut
586              
587             sub neighbors {
588 17     17 1 29 my $self = shift;
589 17         21 my $from = shift;
590 17         24 my @ret = ();
591              
592 17         22 for my $b (@{$self->{bonds}}) {
  17         32  
593 24 50 33     52 push @ret, $b->{to} unless $from && $b->{to} eq $from;
594             }
595 17         50 @ret;
596             }
597              
598             =item $atom->bonds($from)
599              
600             Return a list of bonds. If an atom object $from is specified, bonds to
601             that atom will be excluded from the list.
602              
603             =cut
604              
605             sub bonds {
606 30     30 1 62 my $self = shift;
607 30         40 my $from = shift;
608 30         45 my @ret = ();
609              
610 30         44 for my $b (@{$self->{bonds}}) {
  30         69  
611 35 50 33     102 push @ret, $b->{bond} unless $from && $b->{to} eq $from;
612             }
613 30         97 @ret;
614             }
615              
616             =item $atom->bonds_neighbors($from)
617              
618             Return a list of hash references, representing the bonds and neighbors from the
619             atom. If an atom object $from is specified, it will be excluded from the list.
620             The elements of the hash are 'to', and atom reference, and 'bond', a bond
621             reference. For example,
622              
623             for my $bn ($atom->bonds_neighbors) {
624             print "bond $bn->{bond} point to atom $bn->{to}\n";
625             }
626              
627             =cut
628              
629             sub bonds_neighbors {
630 0     0 1 0 my $self = shift;
631 0         0 my $from = shift;
632 0         0 my @ret = ();
633              
634 0         0 for my $b (@{$self->{bonds}}) {
  0         0  
635 0 0 0     0 push @ret, {%$b} unless $from && $b->{to} eq $from;
636             }
637 0         0 @ret;
638             }
639              
640             =item ($distance, $closest_atom) = $atom->distance($obj)
641              
642             Returns the minimum distance to $obj, which can be an atom, a molecule, or a
643             vector. In scalar context it returns only the distance; in list context it
644             also returns the closest atom found. It can also be called as a function,
645             Chemistry::Atom::distance (which can be exported).
646              
647             =cut
648              
649             sub distance {
650 4     4 1 597 my $self = shift;
651 4         7 my $obj = shift;
652 4         7 my $min_length;
653 4         7 my $closest_atom = $obj;
654              
655 4 100       25 if ($obj->isa('Chemistry::Atom')) {
    50          
    0          
656 3         9 my $v = $self->coords - $obj->coords;
657 3         166 $min_length = $v->length;
658             } elsif ($obj->isa('Math::VectorReal')) {
659 1         4 my $v = $self->coords - $obj;
660 1         52 $min_length = $v->length;
661             } elsif ($obj->isa('Chemistry::Mol')) {
662 0         0 my @atoms = $obj->atoms;
663 0 0       0 my $a = shift @atoms or return undef; # ensure there's at least 1 atom
664 0         0 $min_length = $self->distance($a);
665 0         0 $closest_atom = $a;
666 0         0 for $a (@atoms) {
667 0         0 my $l = $self->distance($a);
668 0 0       0 $min_length = $l, $closest_atom = $a if $l < $min_length;
669             }
670             } else {
671 0         0 croak "atom->distance() undefined for objects of type '", ref $obj,"'";
672             }
673 4 50       156 wantarray ? ($min_length, $closest_atom) : $min_length;
674             }
675              
676             =item $atom->angle($atom2, $atom3)
677              
678             Returns the angle in radians between the atoms involved. $atom2 is the atom in
679             the middle. Can also be called as Chemistry::Atom::angle($atom1, $atom2,
680             $atom3). This function can be exported. Note: if you override this method,
681             you may also need to override angle_deg or strange things may happen.
682              
683             =cut
684              
685             # $a2 is the one in the center
686             sub angle {
687 8 50   8 1 1163 @_ == 3 or croak "Chemistry::Atom::angle requires three atoms!\n";
688 8         12 my @c;
689 8         13 for my $a (@_) { # extract coordinates
690 24 50       41 ref $a or croak "Chemistry::Atom::angle: $a is not an object";
691 24 50       106 push @c, $a->isa("Chemistry::Atom") ? $a->coords :
    100          
692             $a->isa("Math::VectorReal") ? $a :
693             croak "angle: $a is neither an atom nor a vector!\n";
694             }
695 8         65 my $v1 = $c[0] - $c[1];
696 8         302 my $v2 = $c[2] - $c[1];
697 8 100       246 my $l = ($v1->length * $v2->length) or return 0;
698 4         177 acos(($v1 . $v2) / $l);
699             }
700              
701             =item $atom->angle_deg($atom2, $atom3)
702              
703             Same as angle(), but returns the value in degrees. May be exported.
704              
705             =cut
706              
707             sub angle_deg {
708 2     2 1 636 rad2deg(angle(@_));
709             }
710              
711             =item $atom->dihedral($atom2, $atom3, $atom4)
712              
713             Returns the dihedral angle in radians between the atoms involved. Can also be
714             called as Chemistry::Atom::dihedral($atom1, $atom2, $atom3, $atom4). May be
715             exported. Note: if you override this method, you may also need to override
716             dihedral_deg and angle or strange things may happen.
717              
718             =cut
719              
720             sub dihedral {
721 3 50   3 1 394 @_ == 4 or croak "Chemistry::Atom::dihedral requires four atoms!\n";
722 3         5 my @c;
723 3         7 for my $a (@_) { # extract coordinates
724 12 50       42 push @c, $a->isa("Chemistry::Atom") ? $a->coords :
    100          
725             $a->isa("Math::VectorReal") ? $a :
726             croak "angle: $a is neither an atom nor a vector!\n";
727             }
728 3         10 my $v1 = $c[0] - $c[1];
729 3         111 my $v2 = $c[2] - $c[1];
730 3         145 my $v3 = $c[3] - $c[2];
731 3         96 my $x1 = $v1 x $v2;
732 3         115 my $x2 = $v3 x $v2;
733 3         94 my $abs_dih = angle($x1, O(), $x2);
734 3 50       114 $v1 . $x2 > 0 ? $abs_dih : -$abs_dih;
735             }
736              
737             =item $atom->dihedral_deg($atom2, $atom3, $atom4)
738              
739             Same as dihedral(), but returns the value in degrees. May be exported.
740              
741             =cut
742              
743             sub dihedral_deg {
744 1     1 1 10 rad2deg(dihedral(@_));
745             }
746              
747             =item $atom->print
748              
749             Convert the atom to a string representation (used for debugging).
750              
751             =cut
752              
753             sub print {
754 0     0 1 0 my $self = shift;
755 0         0 my ($indent) = @_;
756              
757 16     16   141 no warnings;
  16         53  
  16         3519  
758 0   0     0 $indent ||= 0;
759 0         0 my $bonds = join " ", map {$_->id} $self->bonds;
  0         0  
760 0         0 my $neighbors = join " ", map {$_->id} $self->neighbors;
  0         0  
761             my $coords = $self->{coords}->stringify(
762 0         0 'x3: %g
763             y3: %g
764             z3: %g'
765             );
766              
767 0         0 my $ret = <
768             $self->{id}:
769             symbol: $self->{symbol}
770             name : $self->{name}
771             $coords
772             formal_charge: $self->{formal_charge}
773             bonds: "$bonds"
774             neighbors: "$neighbors"
775             EOF
776 0         0 $ret .= " attr:\n";
777 0         0 $ret .= $self->print_attr($indent+2);
778 0         0 $ret =~ s/^/" "x$indent/gem;
  0         0  
779 0         0 $ret;
780             }
781              
782             =item my $info = $atom->sprintf($format)
783              
784             Format interesting atomic information in a concise way, as specified by
785             a printf-like format.
786              
787             %s - symbol
788             %Z - atomic number
789             %n - name
790             %q - formal charge
791             %h - implicit hydrogen count
792             %v - valence
793             %i - id
794             %8.3m - mass, formatted as %8.3f with core sprintf
795             %8.3x - x coordinate, formatted as %8.3f with core sprintf
796             %8.3y - y coordinate, formatted as %8.3f with core sprintf
797             %8.3z - z coordinate, formatted as %8.3f with core sprintf
798             %% - %
799              
800             =cut
801              
802             sub sprintf {
803 3     3 1 10 my ($atom, $format) = @_;
804 16     16   130 no warnings 'uninitialized'; # don't care if some properties are undefined
  16         30  
  16         7954  
805 3   50     10 $format ||= "%f";
806 3         11 $format =~ s/%%/\\%/g; # escape %% with a \
807 3 0       6 $format =~ s/(?formal_charge || 0/eg; # %q
  0         0  
808 3         12 $format =~ s/(?symbol/eg; # %s
  1         25  
809 3         13 $format =~ s/(?Z/eg; # %Z
  1         4  
810 3         6 $format =~ s/(?name/eg; # %n
  0         0  
811 3         7 $format =~ s/(?hydrogens/eg; # %h
  0         0  
812 3         6 $format =~ s/(?valence/eg; # %v
  0         0  
813 3         6 $format =~ s/(?
814 0 0       0 $1 ? sprintf "%$1f", $atom->mass : $atom->mass/eg; # %m
815 3         12 $format =~ s/(?
816 1 50       6 $1 ? sprintf "%$1f", $atom->x3 : $atom->x3/eg; # %x
817 3         19 $format =~ s/(?
818 1 50       6 $1 ? sprintf "%$1f", $atom->y3 : $atom->y3/eg; # %y
819 3         15 $format =~ s/(?
820 1 50       6 $1 ? sprintf "%$1f", $atom->z3 : $atom->z3/eg; # %z
821 3         11 $format =~ s/(?id/eg; # %i
  0         0  
822 3         6 $format =~ s/\\(.)/$1/g; # other \ escapes
823 3         16 $format;
824             }
825              
826             =item $atom->printf($format)
827              
828             Same as $atom->sprintf, but prints to standard output automatically. Used
829             for quick and dirty atomic information dumping.
830              
831             =cut
832              
833             sub printf {
834 0     0 1   my ($atom, $format) = @_;
835 0           print $atom->sprintf($format);
836             }
837              
838             1;
839              
840             =back
841              
842             =head1 SOURCE CODE REPOSITORY
843              
844             L
845              
846             =head1 SEE ALSO
847              
848             L, L,
849             L, L,
850             L
851              
852             =head1 AUTHOR
853              
854             Ivan Tubert-Brohman Eitub@cpan.orgE
855              
856             =head1 COPYRIGHT
857              
858             Copyright (c) 2005 Ivan Tubert-Brohman. All rights reserved. This program is
859             free software; you can redistribute it and/or modify it under the same terms as
860             Perl itself.
861              
862             =cut
863