Chemistry::Atom - Chemical atoms as objects in molecules

    use Chemistry::Atom;

    my $atom = new Chemistry::Atom(
        id => 'a1',
        coords => [$x, $y, $z],
        symbol => 'Br'
    );

    print $atom->print;

    $atom->symbol("C")->name("CA")->coords(1,2,3);

package Chemistry::Atom;

$VERSION = '0.37';
# $Id: Atom.pm,v 1.44 2009/05/10 19:37:58 itubert Exp $

# Considering to add the following attributes:
# mass_number (A)

use 5.006;
use strict;
use warnings;
use Scalar::Util 'weaken';
use Math::VectorReal qw(O vector);
use Math::Trig;
use Carp;
use base qw(Chemistry::Obj Exporter);
use List::Util qw(first);

our @EXPORT_OK = qw(distance angle dihedral angle_deg dihedral_deg);
our %EXPORT_TAGS = (
      all  => \@EXPORT_OK,
);



my $N = 0; # Number of atoms created so far, used to generate default IDs.

our @ELEMENTS = qw(
    n
    H                                                                   He
    Li  Be                                          B   C   N   O   F   Ne
    Na  Mg                                          Al  Si  P   S   Cl  Ar
    K   Ca  Sc  Ti  V   Cr  Mn  Fe  Co  Ni  Cu  Zn  Ga  Ge  As  Se  Br  Kr
    Rb  Sr  Y   Zr  Nb  Mo  Tc  Ru  Rh  Pd  Ag  Cd  In  Sn  Sb  Te  I   Xe
    Cs  Ba
        La  Ce  Pr  Nd  Pm  Sm  Eu  Gd  Tb  Dy  Ho  Er  Tm  Yb
            Lu  Hf  Ta  W   Re  Os  Ir  Pt  Au  Hg  Tl  Pb  Bi  Po  At  Rn
    Fr  Ra
        Ac  Th  Pa  U   Np  Pu  Am  Cm  Bk  Cf  Es  Fm  Md  No
            Lr  Rf  Db  Sg  Bh  Hs  Mt  Ds  Uuu Uub Uut Uuq Uup Uuh Uus Uuo
);

our %ELEMENTS;
for (my $i = 1; $i < @ELEMENTS; ++$i){
    $ELEMENTS{$ELEMENTS[$i]} = $i;
}
$ELEMENTS{D} = $ELEMENTS{T} = 1;

my %Atomic_masses = (
   "H" => 1.00794, "D" => 2.014101, "T" => 3.016049, "He" => 4.002602,
   "Li" => 6.941, "Be" => 9.012182, "B" => 10.811, "C" => 12.0107,
   "N" => 14.00674, "O" => 15.9994, "F" => 18.9984032, "Ne" => 20.1797,
   "Na" => 22.989770, "Mg" => 24.3050, "Al" => 26.981538, "Si" => 28.0855,
   "P" => 30.973761, "S" => 32.066, "Cl" => 35.4527, "Ar" => 39.948,
   "K" => 39.0983, "Ca" => 40.078, "Sc" => 44.955910, "Ti" => 47.867,
   "V" => 50.9415, "Cr" => 51.9961, "Mn" => 54.938049, "Fe" => 55.845,
   "Co" => 58.933200, "Ni" => 58.6934, "Cu" => 63.546, "Zn" => 65.39,
   "Ga" => 69.723, "Ge" => 72.61, "As" => 74.92160, "Se" => 78.96,
   "Br" => 79.904, "Kr" => 83.80, "Rb" => 85.4678, "Sr" => 87.62,
   "Y" => 88.90585, "Zr" => 91.224, "Nb" => 92.90638, "Mo" => 95.94,
   "Tc" => 98, "Ru" => 101.07, "Rh" => 102.90550, "Pd" => 106.42,
   "Ag" => 107.8682, "Cd" => 112.411, "In" => 114.818, "Sn" => 118.710,
   "Sb" => 121.760, "Te" => 127.60, "I" => 126.90447, "Xe" => 131.29,
   "Cs" => 132.90545, "Ba" => 137.327, "La" => 138.9055, "Ce" => 140.116,
   "Pr" => 140.90765, "Nd" => 144.24, "Pm" => 145, "Sm" => 150.36,
   "Eu" => 151.964, "Gd" => 157.25, "Tb" => 158.92534, "Dy" => 162.50,
   "Ho" => 164.93032, "Er" => 167.26, "Tm" => 168.93421, "Yb" => 173.04,
   "Lu" => 174.967, "Hf" => 178.49, "Ta" => 180.9479, "W" => 183.84,
   "Re" => 186.207, "Os" => 190.23, "Ir" => 192.217, "Pt" => 195.078,
   "Au" => 196.96655, "Hg" => 200.59, "Tl" => 204.3833, "Pb" => 207.2,
   "Bi" => 208.98038, "Po" => 209, "At" => 210, "Rn" => 222,
   "Fr" => 223, "Ra" => 226, "Ac" => 227, "Th" => 232.038,
   "Pa" => 231.03588, "U" => 238.0289, "Np" => 237, "Pu" => 244,
   "Am" => 243, "Cm" => 247, "Bk" => 247, "Cf" => 251,
   "Es" => 252, "Fm" => 257, "Md" => 258, "No" => 259,
   "Lr" => 262, "Rf" => 261, "Db" => 262, "Sg" => 266,
   "Bh" => 264, "Hs" => 269, "Mt" => 268, "Ds" => 271,
);

sub new {
    my $class = shift;
    my %args = @_;

    my $self = bless {
        id => $class->nextID(),
        coords => vector(0, 0, 0),
        Z => 0,
        symbol => '',
        bonds => [],
    }, $class;

    $self->$_($args{$_}) for (keys %args);
    $self;
}

sub nextID {
    "a".++$N; 
}

sub reset_id {
    $N = 0; 
}

sub Z {
    my $self = shift;

    if(@_) {
        $self->{symbol} = $ELEMENTS[$_[0]];
        $self->{Z} = $_[0];
        return $self;
    } else {
        return $self->{Z};
    }
}

sub symbol {
    my $self = shift;

    if(@_) {
        my $symbol = shift;
        $symbol =~ s/ //g;
        $self->{Z} = $ELEMENTS{$symbol};
        $self->{symbol} = $symbol;
        return $self;
    } else {
        return $self->{symbol};
    }
}

sub mass {
    my $self = shift;
    if (@_) {
        $self->{mass} = shift;
        return $self;
    } else {
        if (defined $self->{mass}) {
            return $self->{mass};
        } elsif (defined $self->{mass_number}) {
            if (eval { require Chemistry::Isotope } and 
                my $m = Chemistry::Isotope::isotope_mass(
                    $self->{mass_number}, $self->{Z})
            ) {
                return $m;
            } else {
                return $self->{mass_number};
            }
        } else {
            return $Atomic_masses{$self->symbol};
        }
    }
}

Chemistry::Obj::accessor('mass_number');

sub coords {
    my $self = shift;

    if(@_) {
        if (UNIVERSAL::isa($_[0], "Math::VectorReal")) {
            $self->{coords} = $_[0];
        } elsif (ref $_[0] eq "ARRAY") {
            $self->{coords} = vector(@{$_[0]});
        } else {
            $self->{coords} = vector(@_);
        }
    } else {
        return $self->{coords};
    }
    $self;
}

sub internal_coords {
    my $self = shift;

    if(@_) {
        if (UNIVERSAL::isa($_[0], "Chemistry::InternalCoords")) {
            $self->{internal_coords} = $_[0];
        } elsif (ref $_[0] eq "ARRAY") {
            require Chemistry::InternalCoords;
            $self->{internal_coords} = 
                Chemistry::InternalCoords->new($self, @{$_[0]});
        } else {
            require Chemistry::InternalCoords;
            $self->{internal_coords} = 
                Chemistry::InternalCoords->new($self, @_);
        }
    } else {
        return $self->{internal_coords};
    }
    $self;
}

sub x3 { ($_[0]->coords->array)[0] }
sub y3 { ($_[0]->coords->array)[1] }
sub z3 { ($_[0]->coords->array)[2] }

Chemistry::Obj::accessor('formal_charge');

Chemistry::Obj::accessor('formal_radical');

sub implicit_hydrogens { shift->hydrogens(@_) }

Chemistry::Obj::accessor('hydrogens');

sub total_hydrogens {
    my ($self) = @_;
    no warnings 'uninitialized';
    $self->hydrogens + grep { $_->symbol eq 'H' } $self->neighbors;
}

sub sprout_hydrogens {
    my ($self) = @_;
    for (1 .. $self->implicit_hydrogens || 0) {
        $self->parent->new_bond(
            atoms => [$self, $self->parent->new_atom(symbol => 'H')]);
    }
    $self->implicit_hydrogens(0);
}

sub collapse_hydrogens {
    my ($self) = @_;
    no warnings 'uninitialized';
    my $implicit = $self->implicit_hydrogens;
    for my $nei ($self->neighbors) {
        $nei->delete, $implicit++ if $nei->symbol eq 'H';
    }
    $self->implicit_hydrogens($implicit);
}

my %VALENCE_TABLE = (
    Br => 1, Cl => 1, B => 3, C => 4, N => 3, O => 2, P => 3, S => 2, 
    F => 1, I => 1,
);

# to make it easier to test
sub _calc_implicit_hydrogens {
    my ($self, $symbol, $valence, $charge, $radical) = @_;
    no warnings 'uninitialized';

    my $h_count = $VALENCE_TABLE{$symbol} - $valence;
    # should account for non-kekulized aromatic bonds

    # some common charge situations
    if (($symbol =~ /^(?:[NOSFI]|Cl|Br)$/) && $charge == -1) {
        $h_count--;
    } elsif ($symbol =~ /^[NOSP]$/ && $charge == 1) {
        $h_count++;
    } elsif ($symbol eq 'C' && $charge) {
        $h_count--;
    } elsif ($symbol eq 'B' && $charge == -1) {
        $h_count++;
    }

    # some common radical situations
    if ($radical == 1 or $radical == 3) {
        # carbene, singlet or triplet
        $h_count -= 2;
    } elsif ($radical == 2) {
        # radical (doublet)
        $h_count--;
    }

    $h_count = 0 if $h_count < 0;
    $h_count;
}

sub calc_implicit_hydrogens {
    my ($self) = @_;
    $self->_calc_implicit_hydrogens(
        $self->symbol, $self->explicit_valence, 
        $self->formal_charge, $self->formal_radical,
    );
}

sub add_implicit_hydrogens {
    my ($self) = @_;
    my $h_count = $self->calc_implicit_hydrogens;
    $self->implicit_hydrogens($h_count);
}

Chemistry::Obj::accessor('aromatic');

sub valence {
    my ($self) = @_;
    my $valence = 0;
    $valence += $_->order for $self->bonds;
    $valence += $self->hydrogens || 0;
    $valence;
}

sub explicit_valence {
    my ($self) = @_;
    my $valence = 0;
    $valence += $_->order for $self->bonds;
    $valence;
}

# this method is for internal use only; called by $mol->add_bond
sub add_bond {
    my $self = shift;
    my $bond = shift;
    my %seen;
    #return if first { $_ eq $bond } @{$self->{bonds}}; 

    for my $atom (@{$bond->{atoms}}){ #for each atom...
        if ($atom ne $self) {
            my $b = {to=>$atom, bond=>$bond};
            weaken($b->{to});
            weaken($b->{bond});
            push @{$self->{bonds}}, $b;
        }
    }
}

# make sure the atom doesn't cause circular references
sub _weaken {
    my $self = shift;
    for my $b (@{$self->{bonds}}) {
        weaken($b->{to});
        weaken($b->{bond});
    }
    weaken($self->{parent});
}

# This method is private. Bonds should be deleted from the 
# mol object. These methods should only be called by 
# $bond->delete_atoms, which is called by $mol->delete_bond
sub delete_bond {
    my ($self, $bond) = @_;
    $self->{bonds} = [ grep { $_->{bond} ne $bond } @{$self->{bonds}} ];
}

sub delete {
    my ($self) = @_;
    $self->{parent}->_delete_atom($self);
}

sub parent {
    my $self = shift;
    if (@_) {
        ($self->{parent}) = @_;
        weaken($self->{parent});
        $self;
    } else {
        $self->{parent};
    }
}

sub neighbors {
    my $self = shift;
    my $from = shift;
    my @ret = ();

    for my $b (@{$self->{bonds}}) {
        push @ret, $b->{to} unless $from && $b->{to} eq $from;
    }
    @ret;
}

sub bonds {
    my $self = shift;
    my $from = shift;
    my @ret = ();

    for my $b (@{$self->{bonds}}) {
        push @ret, $b->{bond} unless $from && $b->{to} eq $from;
    }
    @ret;
}

sub bonds_neighbors {
    my $self = shift;
    my $from = shift;
    my @ret = ();

    for my $b (@{$self->{bonds}}) {
        push @ret, {%$b} unless $from && $b->{to} eq $from;
    }
    @ret;
}

sub distance {
    my $self = shift;
    my $obj = shift;
    my $min_length;
    my $closest_atom = $obj;

    if ($obj->isa('Chemistry::Atom')) {
        my $v = $self->coords - $obj->coords;
        $min_length = $v->length;
    } elsif ($obj->isa('Math::VectorReal')) {
        my $v = $self->coords - $obj;
        $min_length = $v->length;
    } elsif ($obj->isa('Chemistry::Mol')) {
        my @atoms = $obj->atoms;
        my $a = shift @atoms or return undef; # ensure there's at least 1 atom
        $min_length = $self->distance($a);
        $closest_atom = $a;
        for $a (@atoms) {
            my $l = $self->distance($a);
            $min_length = $l, $closest_atom = $a if $l < $min_length;
        }
    } else {
        croak "atom->distance() undefined for objects of type '", ref $obj,"'";
    }
    wantarray ? ($min_length, $closest_atom) : $min_length;
}

# $a2 is the one in the center
sub angle {
    @_ == 3 or croak "Chemistry::Atom::angle requires three atoms!\n";
    my @c;
    for my $a (@_) { # extract coordinates
        ref $a or croak "Chemistry::Atom::angle: $a is not an object";
        push @c, $a->isa("Chemistry::Atom") ? $a->coords :
            $a->isa("Math::VectorReal") ? $a : 
                croak "angle: $a is neither an atom nor a vector!\n";
    }
    my $v1 = $c[0] - $c[1];
    my $v2 = $c[2] - $c[1];
    my $l = ($v1->length * $v2->length) or return 0;
    acos(($v1 . $v2) / $l);
}

sub angle_deg {
    rad2deg(angle(@_));
}

sub dihedral {
    @_ == 4 or croak "Chemistry::Atom::dihedral requires four atoms!\n";
    my @c;
    for my $a (@_) { # extract coordinates
        push @c, $a->isa("Chemistry::Atom") ? $a->coords :
            $a->isa("Math::VectorReal") ? $a : 
                croak "angle: $a is neither an atom nor a vector!\n";
    }
    my $v1 = $c[0] - $c[1];
    my $v2 = $c[2] - $c[1];
    my $v3 = $c[3] - $c[2];
    my $x1 = $v1 x $v2;
    my $x2 = $v3 x $v2;
    my $abs_dih = angle($x1, O(), $x2);
    $v1 . $x2 > 0 ? $abs_dih : -$abs_dih;
}

sub dihedral_deg {
    rad2deg(dihedral(@_));
}

sub print {
    my $self = shift;
    my ($indent) = @_;

    no warnings; 
    $indent ||= 0;
    my $bonds = join " ", map {$_->id} $self->bonds;
    my $neighbors = join " ", map {$_->id} $self->neighbors;
    my $coords = $self->{coords}->stringify(
    'x3: %g
        y3: %g
        z3: %g'
    );

    my $ret = <<EOF;
$self->{id}:
    symbol: $self->{symbol}
    name  : $self->{name}
    $coords
    formal_charge: $self->{formal_charge}
    bonds: "$bonds"
    neighbors: "$neighbors"
EOF
    $ret .= "    attr:\n";
    $ret .= $self->print_attr($indent+2);
    $ret =~ s/^/"    "x$indent/gem;
    $ret;
}

sub sprintf {
    my ($atom, $format) = @_;
    no warnings 'uninitialized'; # don't care if some properties are undefined
    $format ||= "%f";
    $format =~ s/%%/\\%/g;              # escape %% with a \
    $format =~ s/(?<!\\)%q/$atom->formal_charge || 0/eg;        # %q
    $format =~ s/(?<!\\)%s/$atom->symbol/eg;                    # %s
    $format =~ s/(?<!\\)%Z/$atom->Z/eg;                         # %Z
    $format =~ s/(?<!\\)%n/$atom->name/eg;                      # %n
    $format =~ s/(?<!\\)%h/$atom->hydrogens/eg;                 # %h
    $format =~ s/(?<!\\)%v/$atom->valence/eg;                   # %v
    $format =~ s/(?<!\\)%(\d*\.?\d*)m/
                $1 ? sprintf "%$1f", $atom->mass : $atom->mass/eg;      # %m
    $format =~ s/(?<!\\)%(\d*\.?\d*)x/
                $1 ? sprintf "%$1f", $atom->x3 : $atom->x3/eg;          # %x
    $format =~ s/(?<!\\)%(\d*\.?\d*)y/
                $1 ? sprintf "%$1f", $atom->y3 : $atom->y3/eg;          # %y
    $format =~ s/(?<!\\)%(\d*\.?\d*)z/
                $1 ? sprintf "%$1f", $atom->z3 : $atom->z3/eg;          # %z
    $format =~ s/(?<!\\)%i/$atom->id/eg;                        # %i
    $format =~ s/\\(.)/$1/g;                             # other \ escapes
    $format;
}

sub printf {
    my ($atom, $format) = @_;
    print $atom->sprintf($format);
}

1;