File Coverage

blib/lib/Chemistry/Atom.pm
Criterion Covered Total %
statement 219 269 81.4
branch 56 90 62.2
condition 21 35 60.0
subroutine 45 50 90.0
pod 30 34 88.2
total 371 478 77.6


line stmt bran cond sub pod time code
1             package Chemistry::Atom;
2              
3             $VERSION = '0.37';
4             # $Id: Atom.pm,v 1.44 2009/05/10 19:37:58 itubert Exp $
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 17     17   63915 use 5.006;
  17         65  
  17         730  
49 17     17   94 use strict;
  17         32  
  17         546  
50 17     17   95 use warnings;
  17         35  
  17         799  
51 17     17   104 use Scalar::Util 'weaken';
  17         38  
  17         3691  
52 17     17   54683 use Math::VectorReal qw(O vector);
  17         169206  
  17         1532  
53 17     17   24815 use Math::Trig;
  17         608188  
  17         7043  
54 17     17   255 use Carp;
  17         41  
  17         13383  
55 17     17   118 use base qw(Chemistry::Obj Exporter);
  17         37  
  17         12277  
56 17     17   214 use List::Util qw(first);
  17         43  
  17         34413  
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 4790 my $class = shift;
131 536         1563 my %args = @_;
132              
133 536         1438 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         13696 $self->$_($args{$_}) for (keys %args);
142 536         3576 $self;
143             }
144              
145             sub nextID {
146 538     538 0 2610 "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 15 my $self = shift;
163              
164 7 100       22 if(@_) {
165 4         135 $self->{symbol} = $ELEMENTS[$_[0]];
166 4         14 $self->{Z} = $_[0];
167 4         26 return $self;
168             } else {
169 3         14 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 1515 my $self = shift;
182              
183 1091 100       2354 if(@_) {
184 531         889 my $symbol = shift;
185 531         658 $symbol =~ s/ //g;
186 531         1970 $self->{Z} = $ELEMENTS{$symbol};
187 531         1305 $self->{symbol} = $symbol;
188 531         2319 return $self;
189             } else {
190 560         2149 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 25 my $self = shift;
212 13 100       33 if (@_) {
213 2         7 $self->{mass} = shift;
214 2         5 return $self;
215             } else {
216 11 100       61 if (defined $self->{mass}) {
    100          
217 2         9 return $self->{mass};
218             } elsif (defined $self->{mass_number}) {
219 2 50 33     3 if (eval { require Chemistry::Isotope } and
  2         872  
220             my $m = Chemistry::Isotope::isotope_mass(
221             $self->{mass_number}, $self->{Z})
222             ) {
223 0         0 return $m;
224             } else {
225 2         9 return $self->{mass_number};
226             }
227             } else {
228 7         18 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 562 my $self = shift;
257              
258 27 100       53 if(@_) {
259 6 100       40 if (UNIVERSAL::isa($_[0], "Math::VectorReal")) {
    100          
260 2         6 $self->{coords} = $_[0];
261             } elsif (ref $_[0] eq "ARRAY") {
262 3         7 $self->{coords} = vector(@{$_[0]});
  3         11  
263             } else {
264 1         4 $self->{coords} = vector(@_);
265             }
266             } else {
267 21         227 return $self->{coords};
268             }
269 6         80 $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 0         0 $self->{internal_coords} =
303 0         0 Chemistry::InternalCoords->new($self, @{$_[0]});
304             } else {
305 0         0 require Chemistry::InternalCoords;
306 0         0 $self->{internal_coords} =
307             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 7 sub x3 { ($_[0]->coords->array)[0] }
324 2     2 1 467 sub y3 { ($_[0]->coords->array)[1] }
325 2     2 1 461 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 75 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 13 my ($self) = @_;
369 17     17   144 no warnings 'uninitialized';
  17         43  
  17         5392  
370 5         17 $self->hydrogens + grep { $_->symbol eq 'H' } $self->neighbors;
  2         7  
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     13 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         27 $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 17     17   98 no warnings 'uninitialized';
  17         50  
  17         3323  
399 2         10 my $implicit = $self->implicit_hydrogens;
400 2         9 for my $nei ($self->neighbors) {
401 4 50       12 $nei->delete, $implicit++ if $nei->symbol eq 'H';
402             }
403 2         11 $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   10929 my ($self, $symbol, $valence, $charge, $radical) = @_;
414 17     17   109 no warnings 'uninitialized';
  17         42  
  17         40651  
415              
416 21         45 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     265 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         8 $h_count++;
424             } elsif ($symbol eq 'C' && $charge) {
425 2         5 $h_count--;
426             } elsif ($symbol eq 'B' && $charge == -1) {
427 1         2 $h_count++;
428             }
429              
430             # some common radical situations
431 21 100 100     108 if ($radical == 1 or $radical == 3) {
    100          
432             # carbene, singlet or triplet
433 3         6 $h_count -= 2;
434             } elsif ($radical == 2) {
435             # radical (doublet)
436 2         5 $h_count--;
437             }
438              
439 21 50       41 $h_count = 0 if $h_count < 0;
440 21         56 $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 5 my ($self) = @_;
453 3         7 $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         9 my $h_count = $self->calc_implicit_hydrogens;
473 3         9 $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 9 my ($self) = @_;
495 3         4 my $valence = 0;
496 3         10 $valence += $_->order for $self->bonds;
497 3   100     12 $valence += $self->hydrogens || 0;
498 3         13 $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 13 my ($self) = @_;
510 6         11 my $valence = 0;
511 6         18 $valence += $_->order for $self->bonds;
512 6         26 $valence;
513             }
514              
515             # this method is for internal use only; called by $mol->add_bond
516             sub add_bond {
517 20     20 0 38 my $self = shift;
518 20         25 my $bond = shift;
519 20         27 my %seen;
520             #return if first { $_ eq $bond } @{$self->{bonds}};
521              
522 20         24 for my $atom (@{$bond->{atoms}}){ #for each atom...
  20         86  
523 40 100       113 if ($atom ne $self) {
524 20         62 my $b = {to=>$atom, bond=>$bond};
525 20         61 weaken($b->{to});
526 20         50 weaken($b->{bond});
527 20         21 push @{$self->{bonds}}, $b;
  20         91  
528             }
529             }
530             }
531              
532             # make sure the atom doesn't cause circular references
533             sub _weaken {
534 112     112   142 my $self = shift;
535 112         125 for my $b (@{$self->{bonds}}) {
  112         330  
536 168         387 weaken($b->{to});
537 168         435 weaken($b->{bond});
538             }
539 112         450 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 38 my ($self, $bond) = @_;
547 28         35 $self->{bonds} = [ grep { $_->{bond} ne $bond } @{$self->{bonds}} ];
  51         136  
  28         56  
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         34 $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 617 my $self = shift;
570 561 100       974 if (@_) {
571 553         1041 ($self->{parent}) = @_;
572 553         1330 weaken($self->{parent});
573 553         2369 $self;
574             } else {
575 8         33 $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 98 my $self = shift;
589 17         24 my $from = shift;
590 17         79 my @ret = ();
591              
592 17         23 for my $b (@{$self->{bonds}}) {
  17         166  
593 24 50 33     89 push @ret, $b->{to} unless $from && $b->{to} eq $from;
594             }
595 17         69 @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 50 my $self = shift;
607 30         42 my $from = shift;
608 30         57 my @ret = ();
609              
610 30         37 for my $b (@{$self->{bonds}}) {
  30         80  
611 35 50 33     154 push @ret, $b->{bond} unless $from && $b->{to} eq $from;
612             }
613 30         164 @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 473 my $self = shift;
651 4         9 my $obj = shift;
652 4         4 my $min_length;
653 4         5 my $closest_atom = $obj;
654              
655 4 100       23 if ($obj->isa('Chemistry::Atom')) {
    50          
    0          
656 3         7 my $v = $self->coords - $obj->coords;
657 3         174 $min_length = $v->length;
658             } elsif ($obj->isa('Math::VectorReal')) {
659 1         4 my $v = $self->coords - $obj;
660 1         39 $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       146 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 1046 @_ == 3 or croak "Chemistry::Atom::angle requires three atoms!\n";
688 8         10 my @c;
689 8         15 for my $a (@_) { # extract coordinates
690 24 50       50 ref $a or croak "Chemistry::Atom::angle: $a is not an object";
691 24 50       166 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         88 my $v1 = $c[0] - $c[1];
696 8         486 my $v2 = $c[2] - $c[1];
697 8 100       279 my $l = ($v1->length * $v2->length) or return 0;
698 4         194 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 692 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 538 @_ == 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       72 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         13 my $v1 = $c[0] - $c[1];
729 3         116 my $v2 = $c[2] - $c[1];
730 3         243 my $v3 = $c[3] - $c[2];
731 3         136 my $x1 = $v1 x $v2;
732 3         193 my $x2 = $v3 x $v2;
733 3         168 my $abs_dih = angle($x1, O(), $x2);
734 3 50       138 $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 9 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 17     17   336 no warnings;
  17         110  
  17         8637  
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 0         0 my $coords = $self->{coords}->stringify(
762             '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 6 my ($atom, $format) = @_;
804 17     17   109 no warnings 'uninitialized'; # don't care if some properties are undefined
  17         36  
  17         16835  
805 3   50     9 $format ||= "%f";
806 3         7 $format =~ s/%%/\\%/g; # escape %% with a \
807 3 0       5 $format =~ s/(?formal_charge || 0/eg; # %q
  0         0  
808 3         9 $format =~ s/(?symbol/eg; # %s
  1         4  
809 3         8 $format =~ s/(?Z/eg; # %Z
  1         4  
810 3         4 $format =~ s/(?name/eg; # %n
  0         0  
811 3         7 $format =~ s/(?hydrogens/eg; # %h
  0         0  
812 3         3 $format =~ s/(?valence/eg; # %v
  0         0  
813 3         4 $format =~ s/(?
814 0 0       0 $1 ? sprintf "%$1f", $atom->mass : $atom->mass/eg; # %m
815 3         9 $format =~ s/(?
816 1 50       7 $1 ? sprintf "%$1f", $atom->x3 : $atom->x3/eg; # %x
817 3         14 $format =~ s/(?
818 1 50       6 $1 ? sprintf "%$1f", $atom->y3 : $atom->y3/eg; # %y
819 3         13 $format =~ s/(?
820 1 50       6 $1 ? sprintf "%$1f", $atom->z3 : $atom->z3/eg; # %z
821 3         13 $format =~ s/(?id/eg; # %i
  0         0  
822 3         5 $format =~ s/\\(.)/$1/g; # other \ escapes
823 3         13 $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 VERSION
843              
844             0.37
845              
846             =head1 SEE ALSO
847              
848             L, L,
849             L, L,
850             L
851              
852             The PerlMol website L
853              
854             =head1 AUTHOR
855              
856             Ivan Tubert-Brohman Eitub@cpan.orgE
857              
858             =head1 COPYRIGHT
859              
860             Copyright (c) 2005 Ivan Tubert-Brohman. All rights reserved. This program is
861             free software; you can redistribute it and/or modify it under the same terms as
862             Perl itself.
863              
864             =cut
865