Chemistry::File::MDLMol - MDL molfile reader/writer

    use Chemistry::File::MDLMol;

    # read a molecule
    my $mol = Chemistry::Mol->read('myfile.mol');

    # write a molecule
    $mol->write("myfile.mol");

    # use a molecule as a query for substructure matching
    use Chemistry::Pattern;
    use Chemistry::Ring;
    Chemistry::Ring::aromatize_mol($mol);

    my $patt = Chemistry::Pattern->read('query.mol');
    if ($patt->match($mol)) {
        print "it matches!\n";
    }

    $atom->attr('mdlmol/test_sub');
    $bond->attr('mdlmol/test_sub');

    $patt->aromatic ?  $bond->aromatic 
        : (!$bond->aromatic && $patt->order == $bond->order);

package Chemistry::File::MDLMol;
$VERSION = '0.21';
# $Id: MDLMol.pm,v 1.13 2009/05/10 19:51:50 itubert Exp $

use base "Chemistry::File";
use Chemistry::Mol;
use Carp;
use List::Util;
use strict;
use warnings;

our $DEBUG = 0;

Chemistry::Mol->register_format(mdl => __PACKAGE__);

# some constants, based on tables from the file format specification

my %OLD_CHARGE_MAP = (
    1 => 3,
    2 => 2,
    3 => 1,
    4 => 0,
    5 => -1,
    6 => -2,
    7 => -3,
);

my %BOND_TYPE_EXPR = (
    4 => '($bond->aromatic)',
    5 => '($bond->order == 1 or $bond->order == 2)',
    6 => '($bond->order == 1 or $bond->aromatic)',
    7 => '($bond->order == 2 or $bond->aromatic)',
    8 => '(1)',                                         # any bond
);

my %BOND_TOPOLOGY_EXPR = (
    1 => '@{$bond->attr("ring/rings")||[]}',
    2 => '! @{$bond->attr("ring/rings")||[]}',
);


sub read_mol {
    my ($self, $fh, %opts) = @_;
    return if $fh->eof;

    %opts = ( slurp => 1, %opts );
    my $mol_class  = $opts{mol_class}  || "Chemistry::Mol";
    my $atom_class = $opts{atom_class} || $mol_class->atom_class;
    my $bond_class = $opts{bond_class} || $mol_class->bond_class;
    local $_;

    my $mol = $mol_class->new();

    # header
    my $name    = <$fh>; chomp $name;
    my $line2   = <$fh>; chomp $line2;
    my $comment = <$fh>; chomp $comment;
    $mol->name($name);
    $mol->attr("mdlmol/line2", $line2);
    $mol->attr("mdlmol/comment", $comment);

    # counts line
    defined ($_ = <$fh>) or croak "unexpected end of file";
    my ($na, $nb) = unpack("A3A3", $_);

    my %old_charges;
    my %old_radicals;

    # atom block
    for my $i (1 .. $na) { # for each atom...
        defined (my $line = <$fh>) or croak "unexpected end of file";
        my ($x, $y, $z, $symbol, $mass, $charge)
            = unpack("A10A10A10xA3A2A3", $line);
        
        $old_charges{$i} = $OLD_CHARGE_MAP{$charge}
            if $OLD_CHARGE_MAP{$charge};
        $old_radicals{$i} = 2
            if $charge && $charge == 4;
        $mol->new_atom(
            symbol         => $symbol, 
            coords         => [$x*1, $y*1, $z*1],
        );
    }

    # bond block
    for my $i (1 .. $nb) { # for each bond...
        no warnings 'numeric';
        defined ($_ = <$fh>) or croak "unexpected end of file";
        my ($a1, $a2, $type, $stereo, $topology, $rxn) 
            = map {$_*1} unpack("A3A3A3A3x3A3A3", $_);
        my $order = $type =~ /^[123]$/ ? $type : 1;
        my $bond = $mol->new_bond(
            type => $type, 
            atoms => [$mol->atoms($a1,$a2)],
            order => $order,
        );
        if ($mol->isa('Chemistry::Pattern')) {
            $self->bond_expr($bond, $i, $type, $topology);
        }
    }

    # properties block
    while (<$fh>) {
        if (/^M  END/ or /^\$\$\$\$/) {
            last;
        } elsif (/^M  (...)/) { # generic extended property handler
            if ($1 eq 'CHG' or $1 eq 'RAD'){ # XXX
                # clear old-style info
                %old_radicals = (); 
                %old_charges  = ();
            }

            my $method = "M_$1";
            $self->$method($mol, $_) if $self->can($method);
        }
    }

    # add old-style charges and radicals if they still apply
    while (my ($key, $val) = each %old_charges) {
        $mol->atoms($key)->formal_charge($val);
    }
    while (my ($key, $val) = each %old_radicals) {
        $mol->atoms($key)->formal_radical($val);
    }

    # make sure we get to the end of the file
    if ($opts{slurp}) {
        1 while <$fh>;
    }

    $mol->add_implicit_hydrogens;

    if ($mol->isa('Chemistry::Pattern')) {
        require Chemistry::Ring;
        Chemistry::Ring::aromatize_mol($mol);
    }

    return $mol;
}

sub bond_expr {
    my ($self, $bond, $i, $type, $topology) = @_;
    my @bond_exprs;
    my $s = $BOND_TOPOLOGY_EXPR{$topology};
    push @bond_exprs, $s if $s;
    $s = $BOND_TYPE_EXPR{$type};
    push @bond_exprs, $s if $s;
    if (@bond_exprs) {
        my $expr = join " and ", @bond_exprs;
        my $sub_txt = <<SUB;
            sub {
                no warnings;
                my (\$patt, \$bond) = \@_;
                $expr;
            };
SUB
        print "MDLMol bond($i) sub: <<<<$sub_txt>>>\n" if $DEBUG;
        $bond->attr('mdlmol/test_sub' => $sub_txt);
        $bond->test_sub(eval $sub_txt);
    } else { # default bond sub
        $bond->test_sub(\&default_bond_test);
    }
}

sub default_bond_test {
    no warnings;
    my ($patt, $bond) = @_;
    $patt->aromatic ?  $bond->aromatic 
        : (!$bond->aromatic && $patt->order == $bond->order);
}


sub M_CHG {
    my ($self, $mol, $line) = @_;
    my ($m, $type, $n, %data) = split " ", $line;
    while (my ($key, $val) = each %data) {
        $mol->atoms($key)->formal_charge($val);
    }
}

sub M_RAD {
    my ($self, $mol, $line) = @_;
    my ($m, $type, $n, %data) = split " ", $line;
    while (my ($key, $val) = each %data) {
        $mol->atoms($key)->formal_radical($val);
    }
}

sub M_ALS {
    my ($self, $mol, $line) = @_;

    return unless $line =~ /^M  ALS (...)(...) (.)/;
    my ($n, $cnt, $exclude) = ($1, $2, $3);
    my @symbols;
    $exclude = $exclude =~ /^[Tt]$/ ? '!' : '';

    # parse the symbols out of the atom list line
    for my $i (0 .. $cnt-1) {
        my $s = substr $_, 16+$i*4, 4;
        $s =~ s/\s+//g;
        push @symbols, $s;
    }

    # save attr
    $mol->atoms($n)->attr('mdlmol/atom_list' => $_);

    # create test sub
    if ($mol->isa('Chemistry::Pattern')) {
        my $sub_txt = <<SUB;
            sub{
                my (\$patt, \$atom) = \@_;
                my \$sym = \$atom->symbol;
                $exclude (List::Util::first {\$sym eq \$_} \@symbols);
            };
SUB
        print "MDLMol atom($n) sub: <<<<$sub_txt>>>\n" if $DEBUG;
        print "MDLMol symbol list: (@symbols)\n" if $DEBUG;
        $mol->atoms($n)->attr('mdlmol/test_sub' => $sub_txt);
        $mol->atoms($n)->test_sub(eval $sub_txt);
    }
}

sub name_is {
    my ($self, $fname) = @_;
    $fname =~ /\.mol$/i;
}


sub file_is {
    my ($self, $fname) = @_;
    $fname =~ /\.mol$/i;
}


sub write_string {
    my ($self, $mol, %opts) = @_;

    no warnings 'uninitialized';
    my $s = sprintf "%s\n      perlmol   \n\n", $mol->name;
    $s .= sprintf "%3i%3i%3i%3i%3i%3i%3i%3i%3i%3i%3i%6s\n", 
        0+$mol->atoms, 0+$mol->bonds, 
        0, 0, 0, 0, 0, 0, 0, 0, 999, "V2000";   # "counts" line

    my $i = 1;
    my %idx_map;
    my @charged_atoms;
    my @radical_atoms;
    for my $atom ($mol->atoms) {
        my ($x, $y, $z) = $atom->coords->array;

        $s .= sprintf 
            "%10.4f%10.4f%10.4f %-3s%2i%3i%3i%3i%3i%3i%3i%3i%3i%3i%3i%3i\n",
            $x, $y, $z, $atom->symbol,
            0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
        push @charged_atoms, $i if $atom->formal_charge;
        push @radical_atoms, $i if $atom->formal_radical;
        $idx_map{$atom->id} = $i++;
    }

    for my $bond ($mol->bonds) {
        my ($a1, $a2) = map {$idx_map{$_->id}} $bond->atoms;
        $s .= sprintf "%3i%3i%3i%3i%3i%3i%3i\n", 
            $a1, $a2, $bond->order,
            0, 0, 0, 0;
    }
    
    while (@charged_atoms) {
        my $n = @charged_atoms > 8 ? 8 : @charged_atoms;
        $s .= "M  CHG  $n";
        for my $key (splice @charged_atoms, 0, $n) {
            $s .= sprintf "%4d%4d", $key, $mol->atoms($key)->formal_charge;
        }
        $s .= "\n";
    }
    while (@radical_atoms) {
        my $n = @radical_atoms > 8 ? 8 : @radical_atoms;
        $s .= "M  RAD  $n";
        for my $key (splice @radical_atoms, 0, $n) {
            $s .= sprintf "%4d%4d", $key, $mol->atoms($key)->formal_radical;
        }
        $s .= "\n";
    }

    $s .= "M  END\n";
    $s;
}



1;