#include "MolMeccano_file.h"

#ifdef MM_WIN32
#pragma warning (disable : 4786)
#endif

#include "Create.h"
#include "User.h"
#include "Model.h"
#include "Model_kit.h"
#include "Atom.h"
#include "Atom_kit.h"
#include "Element.h"
#include "Atom_numerator.h"
#include "Molecule.h"
#include "Molecular_numerator.h"
#include "Molecular_detector.h"
#include "Residue_numerator.h"
#include "Joint.h"
#include "Bond.h"
#include "Bond_kit.h"
#include "Fine_type.h"
#include "Torsion_type.h"
#include "Angle_unit.h"
#include "Distance_restraint.h"
#include "Log.h"

#include <fstream>
#include <sstream>
#include <string>
#include <iomanip>
#include <set>

namespace MM
{
using namespace std;

bool MolMeccano_file::apply_torsions_ = true;

MolMeccano_file MolMeccano_file::prototype_(it_is_a_prototype);

MolMeccano_file::
MolMeccano_file (Prototype const &)
:
    Model_file (Prototype())
{
    add_extension ("mlm");
}

MolMeccano_file::
MolMeccano_file (Model & model)
:
    Model_file (model),
    line_counter_       (-1),
    local_line_counter_ (-1),
    end_reached_        (false),
    detect_molecules_   (false), 
    set_residues_       (false)
{
    add_extension ("mlm");
}

Model_file * MolMeccano_file::
clone (Model & model) const
{
    MolMeccano_file *result = new MolMeccano_file (model);
    
    //for (int i=0;  i<extension_.size();  ++i)
    //    result->extension_.push_back (extension_[i]);
    
    return result;
}

Text const & MolMeccano_file::
type_name () const {static Text t("MolMeccano"); return t;}

void MolMeccano_file::
save (bool to_comment)
{
    using namespace std;
    check_file_name();
    int i, j;

    //Atom_numerator      atom_numerator     (model());
    Residue_numerator   residue_numerator  (model());
    Molecular_numerator numerate_molecules (model());

    ofstream out (name().c_str());

    out << "# MolMeccano format (*.mlm)."                                                       << endl;
    out                                                                                         << endl;
    if (to_comment)
    {
    out << "# Comments (#) and empty lines are ignored."                                        << endl;
    out << "# Keyword begins with '@'."                                                         << endl;
    out << "# '@' must be the first character in line."                                         << endl;
    out                                                                                         << endl;
    }

    if (to_comment)
    {
    out << "# @Model [name] [multiple mode]"                                                    << endl;
    out << "#"                                                                                  << endl;
    out << "# Multiple modes:"                                                                  << endl;
    out << "# separate	    - models in different windows (default)."                           << endl;
    out << "# merge 	    - combine to previous one."                                         << endl;
    out << "# conformation	- other conformation of privious model."                            << endl;
    out << "#                 Many conformations consist of trajectory."                        << endl;
    out << "#"                                                                                  << endl;
    }
    if (name() != "")
    {
    out << "@Model        " << name().c_str()                                                   << endl;
    }

    if (to_comment)
    {
    out << "# @Comment_line <text>"                                                             << endl;
    }
    if (model().kit().comment_line() != "")
    {
    out << "@Comment_line " << model().kit().comment_line().c_str()                             << endl;
    }

    if (to_comment)
    {
    out << "# @Alias <word>"                                                                    << endl;
    }
    if (model().kit().alias() != "")
    {
    out << "@Alias        " << model().kit().alias().c_str()                                    << endl;
    }

    if (to_comment)
    {
    out << "# @Short_alias <word>"                                                              << endl;
    }
    if (model().kit().short_alias() != "")
    {
    out << "@Short_alias  " << model().kit().short_alias().c_str()                              << endl;
    }

    if (to_comment)
    {
    out << "# @Force_field <name>"                                                              << endl;
    out << "#"                                                                                  << endl;
    out << "# Use force field described in 'ForceField' folder (optional)."                     << endl;
    }
    if (model().kit().has_force_field())
    {
    out << "@Force_field  " << model().kit().force_field().name().c_str()                       << endl;
    }
    out                                                                                         << endl;

    if (model().kit().boundary_conditions().is_a_box())
    {
    out << setprecision(15)
        << "@Box  " << model().kit().boundary_conditions().box().a_mod()
        << "  "     << model().kit().boundary_conditions().box().b_mod()
        << "  "     << model().kit().boundary_conditions().box().c_mod()                        << endl;
    }
    out                                                                                         << endl;

    if (to_comment)
    {
    out << "# @Table <type> [options] [comments]"                                               << endl;
    out << "#"                                                                                  << endl;
    out << "# First  string define programming types of colunns."                               << endl;
    out << "# Second string define logical types of colunns."                                   << endl;
    out << "# Other  strings consist of table records."                                         << endl;
    out << "# "                                                                                 << endl;
    //out << "# Each cell should have a value (at least the missed value '-')."                   << endl;
    out                                                                                         << endl;
    out                                                                                         << endl;
    out << "# @Table Atoms [count]"                                                             << endl;
    out << "#"                                                                                  << endl;
    out << "# Count - number of atoms (optional, for memory reservation only)"                  << endl;
    out << "#"                                                                                  << endl;
    out << "# X Y Z Q are obligatory columns."                                                  << endl;
    out << "# ID and Element should have a value or '-'."                                       << endl;
    out << "# All others are optional but if some column item has a value then all preceding "  << endl;
    out << "# items should have a value or '-'."                                                << endl;
    out << "#"                                                                                  << endl;
    out << "# Column 'Flag':"                                                                   << endl;
    out << "# s	     - selected"                                                                << endl;
    out << "# b	     - backbone"                                                                << endl;
    out << "# [0-7][0-7] - mark (in octal form)"                                                << endl;
    out << "#"                                                                                  << endl;
    out << "# Column 'ID' conteins unique atom identifier. "                                    << endl;
    out << "# Usually that is an atom number but it can be any unique word."                    << endl;
    out << "#"                                                                                  << endl;
    out << "# Column 'Type' conteins \"fine atom type\". "                                      << endl;
    out << "# \"Fine type\" is a most specific atom type avalable. "                            << endl;
    out << "# Every atom has one \"fine type\" but can have many especial molecular mechanic "  << endl;
    out << "# types accordin to rules of differnt force fields. \"Fine type\" can be converted" << endl;
    out << "# to molecular mechanics type in a unique fashion but not in opposite direction."   << endl;
    out << "#"                                                                                  << endl;
    }
    out << "@Table Atoms  " << model().atom_count()                                             << endl;
    out << "str  str double     double     double     double str            str str str  double double   double str"       << endl;
    out << "ID Element X          Y          Z          Q    Type           Mol Res Flag R*     Epsilon  Mass   Comment"   << endl;
    
    for (i=0;  i<model().atom_count();  ++i)
    {
        Atom     const& atom = model().atom(i);
        Atom_kit const& kit  = atom.kit();

        kit.set_number (i);              //fix to 'Text id'
        atom.set_local_id (i);              

        out << setprecision(5);
        out << setw(4)          << left            << i                                    << " ";
        out << setw(2)          << left            << atom.element().c_str()               << " ";
        out << setw(10)<< fixed << right           << atom.x()                             << " ";
        out << setw(10)<< fixed                    << atom.y()                             << " ";
        out << setw(10)<< fixed                    << atom.z()                             << " ";
        out << setw(7) << fixed << setprecision(4) << atom.charge()                        << " ";
        out << setw(14)<< fixed << left            << kit.fine_type_name()                 << " ";
        out << setw(3) << fixed << right           << kit.in_molecule().local_id().c_str() << " ";
        out << setw(3) << fixed                    << kit.in_residue ().local_id().c_str() << " ";
        out << setw(4) << fixed                    << kit.flag().text().c_str()            << " ";

        //fix realize
        out << setw(6) << fixed                    << "-" << " ";
        out << setw(8) << fixed                    << "-" << " ";
        out << setw(6) << fixed                    << "-" << " ";

        if (atom.kit().comment_line().is_empty())
        {
            if (to_comment)
            {
                if (atom.kit().has_mm_type())
                    out << atom.kit().mm_type_c_str();
            }
        }
        else
            out << setw(6) << fixed << atom.kit().comment_line().c_str();
        
        //out << setw(9) << fixed << setprecision(6) << current_atom.mass()   << " ";
        out << endl;
    }
    out << endl;

    if (to_comment)
    {
    out << "# @Table Bonds [count]"                                                         << endl;
    out << "#"                                                                              << endl;
    out << "# If force constants defined then they overwrite default force field constants."<< endl;
    out << "#"                                                                              << endl;
    }
    out << "@Table Bonds  "   << model().bond_count()                                       << endl;
    out << "str  str  str   double   double   str"                                          << endl;
    out << "ID1  ID2 Order R-eqv    Force    Comment"                                       << endl;

    //fix Type: Covalent Ion Hydrogen Virtual Elastic

    for (i=0;  i<model().bond_count();  ++i)
    {
        Bond const& bond = model().bond(i);

        if (bond.atom_count() == 2)
        {
            out << setw(4) << bond.atom(0).kit().number()  << " "; //fix to 'Text id'
            out << setw(4) << bond.atom(1).kit().number()  << " "; //fix to 'Text id'
            out << setw(4) << bond.order().c_char()  << " "; 

                //fix realize
            out << setw(8) << fixed                    << "-" << " ";
            out << setw(8) << fixed                    << "-" << " ";
            
            if (bond.kit().comment_line().is_empty())
            {
                if (to_comment)
                {
                    out << setw(8) << fixed << bond.order().c_str() << " "
                        << bond.atom(0).element().c_str()  << "-" 
                        << bond.atom(1).element().c_str();
                }
            }
            else
                out << bond.kit().comment_line().c_str();

            out << endl;
        }
        //fix else multicenter bond
    }
    out << endl;

    if (model().residue_count() > 0)
    {
        if (to_comment)
        {
        out << "# @Table Residues [count]"                             << endl;
        out << "#"                                                     << endl;
        out << "# This is an optional table for residues types setup." << endl;
        out << "#"                                                     << endl;
        }

        out << "@Table Residues  "   << model().residue_count()        << endl;
        out << "str  str   str   str       str"                        << endl;
        out << "ID   Alias Short Name      Comment"                    << endl;

        for (i=0;  i<model().residue_count();  ++i)
        {
            Residue const & res = model().residue(i);
        
            out << setw(4) << res.local_id     ().c_str() << " "; 
            out << setw(5) << res.alias        ().c_str() << " "; 
            out << setw(5) << res.short_alias  ().c_str() << " "; 
            out << setw(10)<< res.name         ().c_str() << " "; 
            out <<            res.comment_line ().c_str();
            out << endl;
        }
        out << endl;
    }

    if (model().joint_count() > 0)
    {
        if (to_comment)
        {
        out << "# @Table Joints"                                                        << endl;
        out << "#"                                                                      << endl;
        out << "# Define three atoms for connecting during the polimer chain building." << endl;
        out << "# First  atoms (stub atoms) will be deleted. "                          << endl;
        out << "# Second atoms (terminal atoms) will form new bond."                    << endl;
        out << "# Third  atoms (directional atoms) defines torsional angle."            << endl;
        out << "#"                                                                      << endl;
        }
        out << "@Table Joints"                                                          << endl;
        out << "str  str  str  str   str"                                               << endl;
        out << "ID1  ID2  ID3  ResID Comment"                                           << endl;

        set <Joint const*>     saved_joint;

        for (i=0;  i<model().residue_count();  ++i)
        {
            Residue const & res = model().residue (i);
        
            for (j=0;  j<res.joint_count();  ++j)
            {
                Joint const & joint = res.joint (j);
                saved_joint.insert (&joint);

                out << setw(4) << joint.stub_atom       ().local_id().c_str() << " "; 
                out << setw(4) << joint.terminal_atom   ().local_id().c_str() << " "; 
                out << setw(4) << joint.directional_atom().local_id().c_str() << " "; 

                if (&joint.in_residue() == &Residue ::none())
                    out << setw(4) << "-"                                     << " ";
                else
                    out << setw(4) << joint.in_residue().local_id().c_str()   << " "; 
                
                out <<            joint.comment_line ().c_str();
                out << endl;
            }
        }
        for (i=0;  i<model().joint_count();  ++i)
        {
            Joint const & joint = model().joint (i);

            set <Joint const*>::iterator it = saved_joint.find (& joint);
            if (it == saved_joint.end())
            {
                out << setw(4) << joint.stub_atom       ().local_id().c_str() << " "; 
                out << setw(4) << joint.terminal_atom   ().local_id().c_str() << " "; 
                out << setw(4) << joint.directional_atom().local_id().c_str() << " ";
                out << endl;
            }
        }
        out << endl;
    }

    //if (to_comment)
    //{
    //out << "# @Table Restraints [count]"                                                    << endl;
    //out << "#"                                                                              << endl;
    //out << "#"                                                                              << endl;
    //}
    //out << "@Table Restraints  "   << model().bond_count()                                  << endl;
    //out << "str  str  str   double   double   str"                                          << endl;
    //out << "ID1  ID2  Type  R-eqv    Force    Comment"                                      << endl;


    if (to_comment)
    {
    out << ""                                                                       << endl;
    out << "#@Table Torsion_types"                                                  << endl;
    out << "#str  str   str   str    str  str"                                      << endl;
    out << "#Type Fine1 Fine2 Fine3 Fine4 Comment"                                  << endl;
    out << ""                                                                       << endl;
    out << "#@Table Torsions"                                                       << endl;
    out << "#str  double str str"                                                   << endl;
    out << "#Type Value  Res Comment"                                               << endl;
    out << ""                                                                       << endl;
    out << "#@Table Distances"                                                      << endl;
    out << ""                                                                       << endl;
    out << ""                                                                       << endl;
    out << "#@Table Angles"                                                         << endl;
    out << ""                                                                       << endl;
    out << ""                                                                       << endl;
    out << "#@Table Dihedrals"                                                      << endl;
    out << ""                                                                       << endl;
    out << ""                                                                       << endl;
    out << "#@Table Groups"                                                         << endl;
    out << ""                                                                       << endl;
    out << ""                                                                       << endl;
    out << "# All after '@End' are ignored (optional)."                             << endl;
    out << "#"                                                                      << endl;
    }

    //out << "@Table MM_types"                                              << endl;
    //out << "str       str str  "       << endl;
    //out << "Fine_type pdb AMBER94"   << endl;


    out << "@End"                                                                   << endl;

    out << endl;

    #ifndef NDEBUG
    for (i=0;  i<model().atom_count();  ++i)
        model().atom(i).kit().set_number (-1);
    #endif//*/
}

/*void MolMeccano_file::
clear ()
{
    Model_file::clear ();
    //torsion_list_.clear ();
}//*/

void MolMeccano_file::
add (Text const & file_name)
{
    torsion_list_.clear ();
    //fix clear local IDs

    line_counter_     = 0;
    detect_molecules_ = false;
    set_residues_     = false;
    end_reached_      = false;

    atom_map_.      clear ();
    molecular_map_. clear ();
    residue_map_.   clear ();

    Own <Accept_new_fine_types> lock (accept_new_fine_types ());

    ifstream   file (file_name.c_str());
    string     line;

    if (file.fail())
    {
        to_user(). error ("Cannot open file " + file_name);
        return;
    }

    try
    {
        while (next_line (file, line) && !end_reached_)
        {
            //log () << "->" << line.c_str() << "<-\n";

            if (line[0] == '@')             // Is a Record Type Indicator (RTI)
                read_title_line (file, line);
            else
                read_data_line (file, line);
        }

        //check_integrity (); //fix
    }
    catch (Text & error)
    {
        to_user().error (Text() + 
            "File " + file_name + "\n"
            "is broken at line " + line_counter_ + ":\n" +
            error);
    }

    if (!end_reached_)
        to_user(). warning ("@End was not reached.");

    if (detect_molecules_)
    {
        Molecular_detector detector(model());
        detector.execute ();
    }

    if (set_residues_)
        for (int i=0;  i<model().molecule_count();  ++i)
            model().molecule(i).reset_residues();

    if (model().atom_count() > 0 && apply_torsions_)
        torsion_list_.apply (model());

    if (FF_.is_not_empty())
        model().kit().set_force_field (FF_);
}

void MolMeccano_file::
check_file_name ()
{
    std::string filename = name().c_str();
    if (filename.find (".mlm" ) == std::string::npos)

        FLAW (Text(name()) + 
              " is a wrong file name.\n"
              "Check file name extension.");//*/


    //fix move to parent (for all classes)
    //fix remove
    //for (int i=0;  i<extension_.size();  ++i)
    //    result->extension_.push_back (extension_[i]);
}

bool MolMeccano_file::
next_line (ifstream & in, string & line)
{
    while (std::getline (in, line))
    {
        ++line_counter_;
        if (line[0] != '#' &&                                      //not a comment
            line.find_first_not_of (" \t\r\n") != std::string::npos) //not a blank line
            return true;
    }
    return false;
}

void MolMeccano_file::
read_title_line (ifstream & in, string & line)
{
    local_line_counter_ = 0;

    string          type;
    istringstream   rti (line); // Record Type Indicator (RTI)
    
    rti >> type;
    if (rti.fail())
        throw Text("'") + line.c_str() + 
                "' - is not a MolMeccano Record Type Indicator.";
    
    //log () << type.c_str() << "\n";
    if (type == "@Table")
    {
        int     count = -1;
        string  table_type;
        rti >> table_type >> count;

        next_line (in, first_table_line_);  
        next_line (in, second_table_line_);
        //next_line (in, third_table_line_);

        if (table_type == "Atoms")
        {
            read_data_line_function_ = &MolMeccano_file::read_Atoms;
        }
        else if (table_type == "Bonds")
        {
            read_data_line_function_ = &MolMeccano_file::read_Bonds;
        }
        else if (table_type == "Residues")
        {
            read_data_line_function_ = &MolMeccano_file::read_Residues;
        }
        else if (table_type == "Joints")
        {
            read_data_line_function_ = &MolMeccano_file::read_Joints;
        }
        else if (table_type == "Restraints")
        {
            read_data_line_function_ = &MolMeccano_file::read_Restraints;
        }
        else if (table_type == "Torsion_types")
        {
            read_data_line_function_ = &MolMeccano_file::read_Torsion_types;
        }
        else if (table_type == "Torsions")
        {
            read_data_line_function_ = &MolMeccano_file::read_Torsions;
        }
        else if (table_type == "Distances")
        {
            read_data_line_function_ = &MolMeccano_file::read_Distances;
        }
        else if (table_type == "Angles")
        {
            read_data_line_function_ = &MolMeccano_file::read_Angles;
        }
        else if (table_type == "Dihedrals")
        {
            read_data_line_function_ = &MolMeccano_file::read_Dihedrals;
        }
        else if (table_type == "Groups")
        {
            read_data_line_function_ = &MolMeccano_file::read_Groups;
        }
        else if (table_type == "MM_types")
        {
            read_MM_types_columns (line);
            read_data_line_function_ = &MolMeccano_file::read_MM_types;
        }
        else
            throw Text("'") + line.c_str() + 
                    "' - is not a MolMeccano Record Type Indicator.";
    }
    else if (type == "@Model")
        read_Model (rti);
    else if (type == "@Alias")
        read_Alias (rti);
    else if (type == "@Short_alias")
        read_Short_alias (rti);
    else if (type == "@Comment_line")
        read_Comment_line (rti);
    else if (type == "@Force_field")
        read_ForceField (rti);
    else if (type == "@Box")
        read_Box (rti);
    else if (type == "@End")
        read_End (rti);
    else 
        throw Text("'") + line.c_str() + 
                "' - is not a MolMeccano Record Type Indicator.";
}

void MolMeccano_file::
//read_data_line (ifstream & in, string & line)
read_data_line (ifstream & , string & line)
{
    (this->*read_data_line_function_) (line);
}

void MolMeccano_file::read_Model (istringstream & ){}

void MolMeccano_file::
read_Alias (istringstream & in)
{
    string text;
    in >> text;
    model().kit().set_alias (text.c_str());
}

void MolMeccano_file::
read_Short_alias (istringstream & in)
{
    string text;
    in >> text;
    model().kit().set_short_alias (text.c_str());
}

void MolMeccano_file::
read_Comment_line (istringstream & in)
{
    string text;
    getline (in, text);
    int start = text.find_first_not_of (" ");
    if (start != string::npos)
        model().kit().set_comment_line (text.substr(start).c_str());
}

void MolMeccano_file::
read_ForceField (istringstream & in)
{
    string text;
    in >> text;
    FF_ = text.c_str();
}

void MolMeccano_file::
read_Box (istringstream & in)
{
    double a_mod, b_mod, c_mod;
    in >> a_mod >> b_mod >> c_mod;
    Boundary_conditions & boundary_conditions = 
        model().kit().boundary_conditions ();
    boundary_conditions.set_type ("Box");
    boundary_conditions.box().set (a_mod, b_mod, c_mod);
}

void MolMeccano_file::
read_End (istringstream & )
{
    end_reached_ = true;

    if (model().kit().comment_line() != "")
        to_user ().status (model().kit().comment_line());
    else
        to_user ().status ("");
}

void MolMeccano_file::
read_Atoms (string & line)
{
    istringstream   stream (line);

    static Text format = 
            "Format of atom data line should be:\n"
            "[ID] [Element] <X> <Y> <Z> <Q> [Type] "
            "[Mol] [Res] [Flag] [R*] [Epsilon] [Mass] [Comment]";//*/

    string  atom_id     = "-";
    string  element     = "-";
    double  x, y, z, charge;
    string  atom_type   = "-";
    string  mol_id      = "-";
    string  res_id      = "-";
    string  flag        = "-";
    string  R0_str      = "-";
    string  Epsilon_str = "-";
    string  mass_str    = "-";
    string  comment;

    stream  >> atom_id >> element >> x >> y >> z >> charge; 
    if (stream.fail())
        throw Text (line.c_str()) + "\n" + format;
    
    stream  >> atom_type >> mol_id  >> res_id  >> flag 
            >> R0_str >> Epsilon_str >> mass_str;

    if (stream)
    {
        string text;
        getline (stream, text);
        string::size_type start = text.find_first_not_of (" \t\n");
        if (start != string::npos)
        {            
            comment = text.substr (start);
        }
    }

    // Element
    if (element == "-")
        element = "Unknown";

    if (!Element::is(element.c_str()))
        throw Text (line.c_str()) + "\n'" + 
            element.c_str() + "' is unknown element.\n" + format;

    Atom & current_atom = new_atom (element.c_str());
    
    // ID
    current_atom.set_local_id (atom_id.c_str());

    if (atom_map_.find (atom_id) == atom_map_.end())
        atom_map_[atom_id] = &current_atom;
    else
        throw Text (line.c_str()) + 
            "\nAtom ID '" + atom_id.c_str() + "' already exist.";

    // <X> <Y> <Z> <Q>
    current_atom.set_x      (x);
    current_atom.set_y      (y);
    current_atom.set_z      (z);
    current_atom.set_charge (charge);

    // Fine atom type
    if (Fine_type::exist(atom_type.c_str()))
        current_atom.kit().set_fine_type (atom_type.c_str());

    else 
    {
        if (accept_new_fine_types_ == Accept_new_fine_types::Ask)
        {
            Text message;
            message += "'";
            message += atom_type.c_str();
            message += "' is unknown fine atom type.\n"
                "Do you want to accept all unknown type from this file?\n"
                "\n"
                "The list of accepted types will be "
                "in the Information->Log panel.";

            bool accept = to_user(). question (message);

            if (accept)
                accept_new_fine_types_ = Accept_new_fine_types::Tmp;
            else
                accept_new_fine_types_ = Accept_new_fine_types::No;
        }

        if (accept_new_fine_types_ == Accept_new_fine_types::Yes)
        {
            Fine_type::add (atom_type.c_str());
            current_atom.kit().set_fine_type (atom_type.c_str());
        }
        else if (accept_new_fine_types_ == Accept_new_fine_types::Tmp)
        {
            Fine_type::add (atom_type.c_str());
            current_atom.kit().set_fine_type (atom_type.c_str());

            Text message;
            //message += "<p>";
            message += "Temporary fine atom type   ";
            message += atom_type.c_str();
            message += "\n";
            log() << message;
       }
    }

    if (current_atom.kit().has_fine_type () &&
        model().kit().has_force_field ())
    {
        const char * mm_type = 
            Fine_type::to (current_atom.kit().fine_type_name (), 
                           model().kit().force_field().name().c_str(), 0);
        if (mm_type != 0)
            current_atom.kit().set_mm_type (mm_type);
    }

    // Molecule
    Molecule *mol;

    if (mol_id != "-" && !mol_id.empty())
    {
        map<string, Molecule*>::iterator mol_it = molecular_map_.find (mol_id);

        if (mol_it == molecular_map_.end())
        {
            mol = &new_molecule();
			mol->set_local_id (mol_id.c_str());
            molecular_map_[mol_id] = mol;
        }
        else
            mol = mol_it->second;

        mol->add_atom (current_atom);
    }
    else
    {
        detect_molecules_ = true;
        //mol = 0;
    }

    // Residue
    Residue * res;

    if (res_id != "-" && !res_id.empty())
    {
        set_residues_ = true;

        map <string, Residue*>::iterator res_it = residue_map_.find (res_id);

        if (res_it == residue_map_.end())
        {
            res = &new_residue();
			res->set_local_id (res_id.c_str());
            residue_map_[res_id] = res;
        }
        else
            res = res_it->second;

        res->add_atom (current_atom);
    }

    // Comment
    if (!comment.empty())
        current_atom.kit().set_comment_line (comment.c_str());

}

void MolMeccano_file::
read_Bonds (string & line)
{
    istringstream   stream (line);

    static Text format = 
            "Format of bond data line should be:\n"
            "<ID1> <ID2> [Order] [R-eqv] [Force] [Comment]";

    string  atom_id_1   = "-";
    string  atom_id_2   = "-";
    string  order       = "-";
    string  R_eqv       = "-";
    string  force       = "-";
    string  comment     = "-";

    stream  >> atom_id_1 >> atom_id_2; 

    if (stream.fail())
        throw Text (line.c_str()) + "\n" + format;

    stream  >> order >> R_eqv >> force; 

    // ID1
    map <string, Atom*>::iterator it = atom_map_.find (atom_id_1);

    if (it == atom_map_.end())
	{
        //throw 
		to_user().warning(
			Text (line.c_str()) + 
            "\nAtom with ID1 '" + atom_id_1.c_str() + "' does not found.");
		return;
	}

    Atom & atom_1 = *it->second;

    // ID2
    it = atom_map_.find (atom_id_2);

    if (it == atom_map_.end())
	{
        //throw 
		to_user().warning(
			Text (line.c_str()) + 
            "\nAtom with ID2 '" + atom_id_2.c_str() + "' does not found.");
		return;
	}

    Atom & atom_2 = *it->second;

    // Order
    Bond & bond = new_bond (atom_1, atom_2, Order (order.c_str()));

    if (stream)
    {
        string text;
        getline (stream, text);
        string::size_type start = text.find_first_not_of (" \t\n");
        if (start != string::npos)
        {
            comment = text.substr (start);
            bond.kit().set_comment_line (comment.c_str());
        }
    }
}

void MolMeccano_file::
read_Residues (string & line)
{
    istringstream   stream (line);

    static Text format = 
            "Format of residue data line should be:\n"
            "<ID> [Alias] [Short] [Name] [Comment]";

    string  id;
    string  alias       = "-";
    string  short_alias = "-";
    string  name        = "-";
    string  comment     = "-";

    stream  >> id; 

    if (stream.fail())
        throw Text (line.c_str()) + "\n" + format;

    stream  >> alias >> short_alias >> name >> comment; 

    Residue * res;
    map <string, Residue*>::iterator res_it = residue_map_.find (id);

    if (res_it == residue_map_.end())
    {
        res = &new_residue();
        residue_map_[id] = res;
    }
    else
        res = res_it->second;

    if (alias != "-" && !alias.empty())
        res->set_alias (alias.c_str());

    if (short_alias != "-" && !short_alias.empty())
        res->set_short_alias (short_alias.c_str());

    if (name != "-" && !name.empty())
        res->set_name (name.c_str());

    if (comment != "-" && !comment.empty())
        res->set_comment_line (comment.c_str());
}

void MolMeccano_file::
read_Joints (string & line)
{
    istringstream   stream (line);

    static Text format = 
            "Format of joint data line should be:\n"
            "<ID1> <ID2> <ID3> [ResID] [Comment]";

    string  id1, id2, id3;
    string  resId   = "-";
    string  comment = "-";

    stream  >> id1 >> id2 >> id3; 

    if (stream.fail())
        throw Text (line.c_str()) + "\n" + format;

    stream  >> resId >> comment; 

    // ID1
    map <string, Atom*>::iterator it = atom_map_.find (id1);

    if (it == atom_map_.end())
        throw Text (line.c_str()) + 
            "\nAtom with ID1 '" + id1.c_str() + "' does not found.";

    Atom & atom_1 = *it->second;

    // ID2
    it = atom_map_.find (id2);

    if (it == atom_map_.end())
        throw Text (line.c_str()) + 
            "\nAtom with ID2 '" + id2.c_str() + "' does not found.";

    Atom & atom_2 = *it->second;

    // ID3
    it = atom_map_.find (id3);

    if (it == atom_map_.end())
        throw Text (line.c_str()) + 
            "\nAtom with ID3 '" + id3.c_str() + "' does not found.";

    Atom & atom_3 = *it->second;


    Joint & joint = new_joint (atom_1, atom_2, atom_3);

    if (comment != "-" && !comment.empty())
        joint.set_comment_line (comment.c_str());

    if (resId != "-" && !resId.empty())
    {
        map <string, Residue*>::iterator res_it = residue_map_.find (resId);

        if (res_it == residue_map_.end())
            throw Text (line.c_str()) + 
                "\nResidue with ID '" + resId.c_str() + "' does not found.";
        else
            res_it->second->add_joint (joint);
    }
}

void MolMeccano_file::
read_Restraints (string & line)
{
    //out << "@Table Restraints  "   << model().bond_count()                                  << endl;
    //out << "str  str  str   double   double   str"                                          << endl;
    //out << "ID1  ID2  Type  R-eqv    Force    Comment"                                      << endl;

    istringstream   stream (line);

    static Text format = 
            "Format of restraint line should be:\n"
            "<ID1> <ID2> <Type> <R-eqv> <Force> [Comment]";

    string  atom_id_1   = "-";
    string  atom_id_2   = "-";
    string  type        = "-";
    double  R_eqv       = 0;
    double  force       = 0;
    string  comment     = "-";

    stream  >> atom_id_1 >> atom_id_2 >> type >> R_eqv >> force;
    if (stream.fail())
        throw Text (line.c_str()) + "\n" + format;

    // ID1
    map <string, Atom*>::iterator it = atom_map_.find (atom_id_1);

    if (it == atom_map_.end())
        throw Text (line.c_str()) + 
            "\nAtom with ID1 '" + atom_id_1.c_str() + "' does not found.";

    Atom & atom_1 = *it->second;

    // ID2
    it = atom_map_.find (atom_id_2);

    if (it == atom_map_.end())
        throw Text (line.c_str()) + 
            "\nAtom with ID2 '" + atom_id_2.c_str() + "' does not found.";

    Atom & atom_2 = *it->second;

    //Distance_restraint & restraint = 
        new_distance_restraint (atom_1, atom_2, type.c_str(), R_eqv, force);

    //if (stream)
    //{
    //    string text;
    //    getline (stream, text);
    //    string::size_type start = text.find_first_not_of (" \t\n");
    //    if (start != string::npos)
    //    {
    //        comment = text.substr (start);
    //        bond.kit().set_comment_line (comment.c_str());
    //    }
    //}
}

void MolMeccano_file::
read_Torsion_types (string & line)
{
    // @Table Torsion_types       
    // str  str   str   str    str  str   
    // Type Fine1 Fine2 Fine3 Fine4 Comment

    istringstream   stream (line);

    static Text format = 
            "Format of Torsion_types data line should be:\n"
            "<Torsion_type> <Fine_type1> <Fine_type2> <Fine_type3> <Fine_type4> [Comment]";

    string  torsion_type, fine_type1, fine_type2, fine_type3, fine_type4;
    string  comment;

    stream  >> torsion_type >> fine_type1 >> fine_type2 >> fine_type3 >> fine_type4; 

    if (stream.fail())
        throw Text (line.c_str()) + "\n" + format;

    stream  >> comment; 

    Torsion_type::add (torsion_type.c_str(), 
        fine_type1.c_str(), fine_type2.c_str(), 
        fine_type3.c_str(), fine_type4.c_str());
}

void MolMeccano_file::
read_Torsions (string & line)
{
    // @Table Torsions        
    // str  double str str    
    // Type Value  Res Comment

    istringstream   stream (line);

    static Text format = 
            "Format of Torsions data line should be:\n"
            "<Torsion_type> <Value> [Residue_ID] [Comment]";

    string  torsion_type, residue_ID = "-";
    double  value;
    string  comment;

    stream  >> torsion_type >> value;

    if (stream.fail())
        throw Text (line.c_str()) + "\n" + format;

    stream  >> residue_ID >> comment; 

    Angle_unit angle;
    angle.set_deg (value);
    torsion_list_.add (torsion_type.c_str(), angle.rad(), residue_ID.c_str());
}

void MolMeccano_file::read_Distances (string & ){}
void MolMeccano_file::read_Angles (string & ){}
void MolMeccano_file::read_Dihedrals (string & ){}
void MolMeccano_file::read_Groups (string & ){}

void MolMeccano_file::
read_MM_types_columns (string & line)
{
    istringstream   stream (second_table_line_);
    mm_types_column_.clear ();

    string first, current;
    stream >> first;

    if (first != "Fine_type")
        throw Text (line.c_str()) + 
            "\nFirst column in 'MM_types' table should be 'Fine_type'.";

    while (stream >> current)
        mm_types_column_.push_back (current);
}

void MolMeccano_file::
read_MM_types (string & line)
{
    istringstream   stream (line);
    string fine_type;
    string mm_type;
    string charge_str;
    double charge;
    
    stream >> fine_type;

    if (!Fine_type::exist (fine_type.c_str()))
        Fine_type::add (fine_type.c_str());

    for (std::vector <std::string>::size_type i=0;  
         (stream >> mm_type >> charge_str) && i<mm_types_column_.size();  ++i)
    {
        if (charge_str == "-")
            charge = Fine_type::no_charge;
        else
        {
            istringstream   charge_stream (charge_str);
            charge_stream >> charge;

            if (!charge_stream)
                break;
        }

        Fine_type::map (fine_type.c_str(), 
                        mm_type.c_str(),
                        charge,
                        mm_types_column_[i].c_str());
    }
}

}//MM

