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

#include "M_DynaMix_mol_file.h"

#include <math.h>
#include <fstream>
#include <sstream>
#include <string>
#include <iomanip>
#include <map>

//#include "Defs.h"
#include "AMBER_energy.h"
#include "Create.h"
#include "Atom.h"
#include "Atom_kit.h"
#include "Bond.h"
#include "Model.h"
#include "Model_kit.h"
#include "Molecule_impl.h"
#include "C_str.h"
#include "Force_field.h"
#include "Mach_eps.h"

#include "Log.h"


namespace MM
{
//M_DynaMix_mol_file M_DynaMix_mol_file::prototype_(Prototype());
M_DynaMix_mol_file M_DynaMix_mol_file::prototype_(it_is_a_prototype);


M_DynaMix_mol_file::
M_DynaMix_mol_file()
{
}

M_DynaMix_mol_file::
M_DynaMix_mol_file (Prototype const &)
:
    Model_file (Prototype())
{
    add_extension ("mmol");//fix mmol
}

M_DynaMix_mol_file::
M_DynaMix_mol_file (Model & model)
:
    Model_file (model)
{
}

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

bool M_DynaMix_mol_file::
next_line (std::ifstream & in, std::string & line)
{
    while (std::getline (in, line))
    {
        if (line[0] != '#')
            return true;
    }
    return false;
}

void M_DynaMix_mol_file::
add (Text const & file_name)
{
    using namespace std;

    Molecule & mol = new_molecule();

    //ifstream 		in (file_name.c_str());
    //string     		text;
    int     atom_count = 0;
    int     line_count = 0;
    string  element;
    double  x, y, z, mass, charge;
    int     broken_line = -1;

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

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

    do {getline (file, line); ++line_count;} while (line[0] == '#');
    istringstream in (line);
    in >> atom_count;

    for (int i=0;  i<atom_count;  ++i)
    {
        if (getline (file, line))
        {
            ++line_count;

            if (line[0] == '#')
            {
                --i;
                continue;
            }

            istringstream in (line);
            in >> element >> x >> y >> z >> mass >> charge;
            if (in.fail())
            {
                broken_line = line_count;
                break;
            }

            string current_element = element;
            if (element.length() > 2)
            {
                current_element  = element[0];
                current_element += element[1];
            }
            if (!Element::is (current_element.c_str()))
            {
                current_element = current_element[0];

                if (!Element::is (current_element.c_str()))
                {
                    broken_line = line_count;
                    break;
                }
            }

            Atom & current_atom = new_atom (Element (current_element.c_str()));

            current_atom.set_x (x);
            current_atom.set_y (y);
            current_atom.set_z (z);

            mol.add_atom (current_atom);
        }
        else
        {
            broken_line = i;
            break;
        }
    }//*/

    if (broken_line != -1)
        to_user().error (Text() + 
            "File " + file_name + " is broken at line " + 
            (broken_line /*+ 3*/) + ".");

    Text formula = mol.molecular_formula();
    model().kit().set_MDynaMix_file (formula, file_name);
}

void M_DynaMix_mol_file::
check_file_name ()
{
    std::string filename = name().c_str();
    if( filename.find( ".mmol" ) == std::string::npos )//fix extension()
    //fix to end
        FLAW (Text(name()) + " is a wrong file name.\n"
              "Check file name extension.");
}

void M_DynaMix_mol_file::
save (bool)
{
    save (model(), model());
}

//fix remove
Text & global_sigma_rule   ()
{
    static Text instance;
    return instance;
}

Text & global_epsilon_rule ()
{
    static Text instance;
    return instance;
}

void M_DynaMix_mol_file::
save (const Atom_group & atoms, const Bond_group & bonds)
{
    using namespace std;
    check_file_name();

    int     i=0;
    bool    warnings_about_partial_charges = true;
    Text    warning;

    //Force_field & force_field = Force_field::singleton(); //fix
    Force_field & force_field = model().kit().force_field(); //fix
//    Own <Force_field> janitor (new Force_field);
//    Force_field & force_field = janitor();
    double  lIEps = 1 / force_field.eps;

    //fix remove
    global_sigma_rule()   = force_field.sigma_rule_;
    global_epsilon_rule() = force_field.epsilon_rule_;

    ofstream out (name().c_str());
    //out.width (10);

    out << "# Atoms" << endl;
//    out << "\t" << model().atom_count() << std::endl;
    out << "\t" << atoms.atom_count() << std::endl;

    out << "#    X          Y          Z         M          Q      sigma   epsilon Type Number\n";
    out << "#              (A)                  (u)        (e)      (A)     (kJ/M)\n";

    double r_min2sigma = 2. / pow (2, 1./6.);
    //double r_min2sigma = 1.781797436280678609480452411181;

    //for (i=0;  i<model().atom_count();  ++i)
    for (i=0;  i<atoms.atom_count();  ++i)
    {
        //Atom & current_atom = model().atom(i);
        const Atom & current_atom = atoms.atom(i);
        current_atom.kit().set_number (i+1);
        Text atom_name;                                      

        /*if (current_atom.kit().has_hyperchem_attributes())
            name = current_atom.kit().hyperchem().name();

        else if (current_atom.kit().has_pdb_attributes())
            name = current_atom.kit().pdb().name();

        else//*/
            //name = Text(current_atom.element().c_str()) + (i+1);
            atom_name = Text(current_atom.element().c_str());

        int mm_type = current_atom.kit().mm_type();
        double r_min = force_field.R_eps [mm_type]. R_star; 
        double e_min = force_field.R_eps [mm_type]. Eps;

        out << setprecision(7);
        out << setw(2) << left                     << atom_name.c_str() << right << " ";
        out << setw(10)<< fixed                    << current_atom.x()      << " ";
        out << setw(10)<< fixed                    << current_atom.y()      << " ";
        out << setw(10)<< fixed                    << current_atom.z()      << " ";
        out << setw(9) << fixed << setprecision(6) << current_atom.mass()   << " ";

        double original_charge = current_atom.charge();
        double charge          = original_charge;
        
        Force_field::Partial_charge const& ff_charge 
            = force_field.partial_charge (mm_type);
        
        if (ff_charge.was_defined() && 
            original_charge != ff_charge.value())
        {
            warning = "Original Q=";
            warning += original_charge;
            warning += " !!!";

            charge = ff_charge.value();
            
            if (warnings_about_partial_charges)
            {
                Text message;
                message += "Atomic partial charges are different from Force Field charges.\n";
                message += "Actually Force Field charges in use.\n";
                message += "Look in file ";
                message += name();
                to_user(). warning (message);
                warnings_about_partial_charges = false;
            }
        }
        charge *= lIEps;
        out << setw(8) << fixed << setprecision(5) << charge                << " ";
        
        out << setw(6) << fixed                    << r_min * r_min2sigma   << " ";   //sigma     
        out << setw(8) << fixed << setprecision(6) << e_min * Joule_Cal     << " ";   //epsilon   
        out << setw(3) << current_atom.kit().mm_type_c_str()                << " ";
        out << setw(4) << (i+1)                                             << " ";
        out << warning.c_str() << endl;
    }

    out << "# References" << endl;
    out << "\t" << 1 << endl;
    out << "File created by M.DynaMix.shell." << endl;


    //AMBER_energy energy (model(), model(), force_field);
    AMBER_energy energy (atoms, bonds, force_field);

    out << "# Bonds" << endl;
    out << "\t" << energy.bond_count() << std::endl;
    out << "#      N1     N2  R-eqv    Force    Type" << std::endl;

    for (i=0;  i<energy.bond_count();  ++i)
    {
        out << "0" << "  ";
        out << setw(6) << energy.bond(i).first_atom ().kit().number()  << " ";
        out << setw(6) << energy.bond(i).second_atom().kit().number()  << " ";

        out << setw(7) << fixed << setprecision(4) << energy.bond(i).equilibrium_distance()        << " ";
        out << setw(10)<< fixed << setprecision(3) << energy.bond(i).force_constant() * Joule_Cal  << "  ";
        
        ostringstream tor;
        tor << energy.bond(i).first_atom ().kit().mm_type_c_str()
            << "-"
            << energy.bond(i).second_atom().kit().mm_type_c_str() << " ";
        out << setw(6) << tor.str();

        if (fabs(energy.bond(i).force_constant()) < d_mach_eps_2)
            out << "  Zero force constant !!!";

        out << endl;
    }

    out << "# Angles" << endl;
    out << "\t" << energy.angle_count() << std::endl;
    out << "#   N1     N2     N3  A-eqv  Force    Type\n";

    for (i=0;  i<energy.angle_count();  ++i)
    {
        out << setw(6) << energy.angle(i).first_atom ().kit().number()         << " ";
        out << setw(6) << energy.angle(i).second_atom().kit().number()         << " ";
        out << setw(6) << energy.angle(i).third_atom ().kit().number()         << " ";

        out << setw(6) << fixed << setprecision(2) << energy.angle(i).equilibrium_value() * 180 / __Pi__ << " ";
        out << setw(7) << fixed << setprecision(3) << energy.angle(i).force_constant() * Joule_Cal       << "  ";

        ostringstream tor;
        tor << energy.angle(i).first_atom ().kit().mm_type_c_str()
            << "-"
            << energy.angle(i).second_atom().kit().mm_type_c_str()
            << "-"
            << energy.angle(i).third_atom ().kit().mm_type_c_str() << " ";
        out << setw(9) << tor.str();

        if (fabs(energy.angle(i).force_constant()) < d_mach_eps_2)
            out << "  Zero force constant !!!";

        out << endl;
    }

    out << "# Dihedrals" << endl;

    int term_count = 0;
    for (i=0;  i<energy.torsion_count();  ++i)
    {
        if (energy.torsion(i).is_valid())
        {
            int terms = energy.torsion(i).fourier_terms();

            for (int k=0;  k<terms;  ++k)
            {
                //if (energy.torsion(i).force_constant(k) < d_mach_eps_2)
                //    continue;
                ++term_count;
            }
        }
        else 
            ++term_count;
    }

    //return;
    //fix - reimplement in AMBER_energy
    // Improper tortions.
    for (i=0;  i<atoms.atom_count();  ++i)
    {
        const Atom & current_atom = atoms.atom(i);

        if (current_atom.bound_count() == 3)
        {
            const int mm_type_i = current_atom.bound(0).kit().mm_type();
            const int mm_type_j = current_atom.bound(1).kit().mm_type();
            const int mm_type_k = current_atom.         kit().mm_type();
            const int mm_type_l = current_atom.bound(2).kit().mm_type();

            int i, j, l;
            const Force_field::Imp_Tor *improper_tortion = force_field.imp_tor
                (mm_type_i, mm_type_j, mm_type_k, mm_type_l,   i, j, l);

            if (improper_tortion != 0)
                ++term_count;
        }
    }

    out << "\t" << term_count << std::endl;               
    out << "#   N1     N2     N3     N4    A-eqv    Force n    Type\n";

    for (i=0;  i<energy.torsion_count();  ++i)
    {
        if (energy.torsion(i).is_valid())
        {
            int terms = energy.torsion(i).fourier_terms();

            for (int k=0;  k<terms;  ++k)
            {
                //if (energy.torsion(i).force_constant(k) < d_mach_eps_2)
                //    continue;

                out << setw(6) << energy.torsion(i).first_atom ().kit().number() << " ";
                out << setw(6) << energy.torsion(i).second_atom().kit().number() << " ";
                out << setw(6) << energy.torsion(i).third_atom ().kit().number() << " ";
                out << setw(6) << energy.torsion(i).fourth_atom().kit().number() << " ";

                out << setw(8) << fixed << setprecision(2) 
                    << energy.torsion(i).phase_angle(k) * 180 / __Pi__  << " ";
                out << setw(8) << fixed << setprecision(4) 
                    << Joule_Cal * energy.torsion(i).force_constant(k) 
                        / energy.torsion(i).multiplicity(k)             << " ";
                out << energy.torsion(i).n_fold_symmetry(k)             << " ";

                ostringstream tor;
                tor << energy.torsion(i).first_atom ().kit().mm_type_c_str()
                    << "-"
                    << energy.torsion(i).second_atom().kit().mm_type_c_str()
                    << "-"
                    << energy.torsion(i).third_atom ().kit().mm_type_c_str()
                    << "-"
                    << energy.torsion(i).fourth_atom().kit().mm_type_c_str() << " ";
                out << setw(12) << tor.str();

                if (energy.torsion(i).force_constant(k) < d_mach_eps_2)
                    out << "  Zero force constant !";

                out << endl;
            }
        }
        else
        {
            out << setw(6) << energy.torsion(i).first_atom ().kit().number()   << " ";
            out << setw(6) << energy.torsion(i).second_atom().kit().number()   << " ";
            out << setw(6) << energy.torsion(i).third_atom ().kit().number()   << " ";
            out << setw(6) << energy.torsion(i).fourth_atom().kit().number()   << " ";

            out << setw(8) << "0" << " ";
            out << setw(8) << "0" << " ";
            out << "1 ";

            ostringstream tor;
            tor << energy.torsion(i).first_atom ().kit().mm_type_c_str()
                << "-"
                << energy.torsion(i).second_atom().kit().mm_type_c_str()
                << "-"
                << energy.torsion(i).third_atom ().kit().mm_type_c_str()
                << "-"
                << energy.torsion(i).fourth_atom().kit().mm_type_c_str() 
                << " " << flush;
            out << setw(12) << tor.str();

            //if (fabs(energy.angle(i).force_constant()) < d_mach_eps_2)
            out << "  No such torsion !!!";
            out << endl;
        }
    }

    //return;
    //fix - reimplement in AMBER_energy
    // Improper torsions.
    for (i=0;  i<atoms.atom_count();  ++i)
    {
        const Atom & current_atom = atoms.atom(i);

        if (current_atom.bound_count() == 3)
        {
            const int mm_type_i = current_atom.bound(0).kit().mm_type();
            const int mm_type_j = current_atom.bound(1).kit().mm_type();
            const int mm_type_k = current_atom.         kit().mm_type();
            const int mm_type_l = current_atom.bound(2).kit().mm_type();

#ifndef NDEBUG
            Element el_i = current_atom.bound(0).element();
            Element el_j = current_atom.bound(1).element();
            Element el_k = current_atom.element();
            Element el_l = current_atom.bound(2).element();

            Text ijkl;
            ijkl += el_i.c_str();
            ijkl += el_j.c_str();
            ijkl += el_k.c_str();
            ijkl += el_l.c_str();
#endif
            int i, j, l;
            const Force_field::Imp_Tor *improper_tortion = force_field.imp_tor
                (mm_type_i, mm_type_j, mm_type_k, mm_type_l,   i, j, l);

            if (improper_tortion != 0) //fix ijl order
            {
                out << setw(6) << current_atom.bound(0).kit().number() << " ";
                out << setw(6) << current_atom.bound(1).kit().number() << " ";
                out << setw(6) << current_atom         .kit().number() << " ";
                out << setw(6) << current_atom.bound(2).kit().number() << " ";

                out << setw(8) << fixed << setprecision(2) << 180 << " ";
                out << setw(8) << fixed << setprecision(4)
                    << Joule_Cal * improper_tortion->half_V   << " ";
                out << improper_tortion->N                    << " ";

                ostringstream tor;
                tor << current_atom.bound(0).kit().mm_type_c_str()
                    << "-"
                    << current_atom.bound(1).kit().mm_type_c_str()
                    << "-"
                    << current_atom.         kit().mm_type_c_str()
                    << "-"
                    << current_atom.bound(2).kit().mm_type_c_str() << " ";
                out << setw(12) << tor.str();
                
                out << " - Improper torsion. ";

                out << current_atom.bound(0).element().c_str() << "-";
                out << current_atom.bound(1).element().c_str() << "-";
                out << current_atom.         element().c_str() << "-";
                out << current_atom.bound(2).element().c_str();

                out << endl;
            }
        }
    }

    out << setprecision(6);
    out << "#  Scaling of 1-4 electrostatic interactions"   << endl;
    out << "sel14  " << force_field.E_Static_Scale          << endl;
    out << "#  Scaling of 1-4 Lennard-Jones interactions"   << endl;
    out << "slj14  " << force_field.VdW_Scale               << endl;

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

}//MM
