/*
    ChemCon - molecular mechanics and molecular graphics
    Copyright (C) 1998-2002  Alexei Nikitin

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/
#ifdef WIN32
#pragma warning(disable : 4786)
#endif

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

//#include <iofswd>
#include <map>
//#include "Force_field.h"
#include <iostream>
//#include "Force_field.h"

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

#include "Force_field.h"

#include "User.h"
#include "Log.h"

#include "MOT.h"
#include "File_ptr.h"
#include "Path.h"
#include "Common_interactions.h"

#include <math.h>//*/

//#include "Defs.h"

#define TESDE 0

namespace MM
{

namespace
{
    char buf[BUF_SIZE];
}
/*
//------------ Smooth ------------------------------------------------------

Smooth::
Smooth()
  :  bond_  ( 1 ),
    angle_    ( 1 ),
    torsion_  ( 1 ),
    impTor_  ( 1 ),
    VdWaals_ ( 1 ),
    e_static_( 1 )
  {
  }
//*/
//==================== Force_field ========================================

void Force_field::
start_work (bool warnings) const                         
{
    problem_ = 0; 
    show_user_warning_ = warnings;
    //if (show_user_warning_)
    //    std::ofstream error_file ((Path::report()+"FF_error.txt").c_str());
}

/*
Force_field & Force_field::
singleton()
{
    static Force_field *instance = new Force_field;
    return *instance;
    //static Force_field instance;
    //return instance;
}//*/

Force_field::Force_field()
:
    show_user_warning_(true),
    problem_(0),
    success_(false),
    //error_file_((Path::report()+"FF_error.txt").c_str()),
    Angle_par( 4000 ), //fix
    Torsion_par( 4000 )
{
    using namespace std;

    clear_bonds ();
    clear_VdW   ();

    std::ofstream error_file ((Path::report()+"FF_error.txt").c_str());

    Text file_name (Path::force_field());
    file_name += "default.ff";
    ifstream in (file_name.c_str());
    if (!in)
        FLAW (Text() + "Can not open file " + file_name);

    string force_field_name;
    //in >> force_field_name;
    //init (force_field_name.c_str());

    in >> force_field_name;
    
    while (in)
    {
        init (force_field_name.c_str());
        in >> force_field_name;
    }

    //string line;
    //locale loc ( "German_Germany" );

    //while (getline (in, line))
    //{
    //    if (isalpha (line[0], loc))
    //    {
    //    }
    //}

}

Force_field::
Force_field (Text const & type)
:
    show_user_warning_(true),
    problem_(0),
    success_(false),
    //error_file_((Path::report()+"FF_error.txt").c_str()),
    Angle_par( 4000 ), //fix
    Torsion_par( 4000 )
{
    clear_bonds ();
    clear_VdW   ();

    init (type);
}

void Force_field::
init (Text const & force_field_name)
{
    using namespace std;
    name_ = force_field_name;

    I_ = 0.;
    I_cut_ = 0.2;
    //eps            = 1.;      // Epsilon
    //E_Static_Scale = 1./1.2;  // Electrostatic scale    //fix 1/S
    //VdW_Scale      = .5;      // VdW scale

    Text file_name (Path::force_field());
    file_name += force_field_name;
    file_name += ".ff";
    ifstream in (file_name.c_str());
    if (!in)
    {
        Text message ("Force field '");
        message += force_field_name;
        message += "' was not found. Only atoms types were associated.";
        to_user().warning (message);
        //FLAW (Text() + "Can not open file " + file_name);
        return;
    }

    string exceptions = "";
    string line;
    getline (in, line);
    getline (in, line);
    while (getline (in, line))
    {
        stringstream in_line (line);
        string type, file, remark;

        in_line >> type >> file >> remark;
        
        Text file_name (Path::force_field());
        file_name += file.c_str();

        if      (type == "atom_type")
        {
            //log () << "atom_type from " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            ::mot().load (file_name.c_str()); //fix remove c_str()
        }
        else if (type == "bond_stretching")
        {
            //log () << "bond_stretching from " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            read_bonds (file_name.c_str());
        }
        else if (type == "angle_bending")
        {
            //log () << "angle_bending from " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            read_angles (file_name.c_str());
        }
        else if (type == "torsional")
        {
            //log () << "torsional from " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            read_torsions (file_name.c_str());
        }
        else if (type == "improper_torsional")
        {
            //log () << "improper_torsional from " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            read_imp_tor (file_name.c_str());
        }
        else if (type == "Van_der_Waals")
        {
            //log () << "Van_der_Waals from " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            read_VdW (file_name.c_str());
        }
        else if (type == "VdW_exceptions")
        {
            exceptions = file_name.c_str();
        }
        else if (type == "dielectric")
        {
            //log () << "dielectric " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            stringstream in_line_2 (line);
            in_line_2 >> type >> eps;
        }
        else if (type == "14_electrostatic")
        {
            //log () << "14_electrostatic " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            stringstream in_line_2 (line);
            in_line_2 >> type >> E_Static_Scale;
        }
        else if (type == "14_Van_der_Waals")
        {
            //log () << "14_Van_der_Waals " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            stringstream in_line_2 (line);
            in_line_2 >> type >> VdW_Scale;
        }
        else if (type == "partial_charges")
        {
            //log () << "partial_charges " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            stringstream in_line_2 (line);
            load_partial_charges (file_name);
        }
        else if (type == "induction")
        {
            //log () << "partial_charges " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            stringstream in_line_2 (line);
            read_induction (file_name);
        }
        else if (type == "sigma_rule")
        {
            //log () << "sigma_rule " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            stringstream in_line_2 (line);
            string sigma_rule;
            in_line_2 >> type >> sigma_rule;
            sigma_rule_ = sigma_rule.c_str();
        }
        else if (type == "epsilon_rule")
        {
            //log () << "epsilon_rule " 
            //       << file.c_str() << ". " << remark.c_str() << "\n";
            stringstream in_line_2 (line);
            string epsilon_rule;
            in_line_2 >> type >> epsilon_rule;
            epsilon_rule_ = epsilon_rule.c_str();
        }
        else if (type == "Rough_hydrophobic")
        {
            stringstream in_line_2 (line);
            in_line_2 >> type >> Rm_ >> l_ >> Rm_hph_ >> K_rough_;
        }
        else if (type == "I")
        {
            stringstream in_line_2 (line);
            in_line_2 >> type >> I_ >> I_cut_;
        }
        else 
        {
            log () << "===> Unknown keyword '" << type.c_str() << "' in line: ";
            log () << line.c_str() << "\n";
        }
    }//*/

    if (exceptions != "")
        read_VdW_exceptions (exceptions.c_str());

    save_as ("last");
}

void Force_field::
save_as (Text const & force_field_name) const
{
    using namespace std;

    Text file_name (Path::force_field());
    file_name += force_field_name;
    file_name += ".ff";
    ofstream out (file_name.c_str());
    if (!out)
        FLAW (Text() + "Can not open file " + file_name);

    out << "string          \tstring                     \tstring"                          << endl;
    out << "Type            \tFile                       \tRemarks"                         << endl;
    out << "atom_type       \t" << save_atom_types        (force_field_name).c_str() << "\t"<< endl;
    out << "bond_stretching   \t" << save_bond              (force_field_name).c_str() << "\t"<< endl;
    out << "angle_bending     \t" << save_angle             (force_field_name).c_str() << "\t"<< endl;
    out << "torsional         \t" << save_torsional         (force_field_name).c_str() << "\t"<< endl;
    out << "improper_torsional\t" << save_improper_torsions (force_field_name).c_str() << "\t"<< endl;
    out << "Van_der_Waals     \t" << save_Van_der_Waals     (force_field_name).c_str() << "\t"<< endl;
    out << "partial_charges   \t" << save_partial_charges   (force_field_name).c_str() << "\t"<< endl;
    out << "dielectric        \t" << eps                                               << "\t"<< endl;
    out << "14_electrostatic  \t" << E_Static_Scale                                      << "\t"<< endl;
    out << "14_Van_der_Waals  \t" << VdW_Scale                                         << "\t"<< endl;
    out << "sigma_rule        \t" << sigma_rule_.c_str()                               << "\t"<< endl;
    out << "epsilon_rule      \t" << epsilon_rule_.c_str()                             << "\t"<< endl;

}

Text Force_field::
save_atom_types (Text const & force_field_name) const
{
    using namespace std;

    Text file_name;
    file_name += force_field_name;
    file_name += "_0_atom_types.xls";
    ofstream out ((Path::force_field() + file_name).c_str());
    if (!out)
        FLAW (Text() + "Can not open file " + file_name);


    out << "string \tdouble \tstring\t"    << endl;
    out << "Type   \tMass   \tReference\t" << endl;
    for (int i=0;  i<mot().get_counter();  ++i)
    {
        out << ::mot().itoa (i) << "\t" << "\t" << endl;
    }

    return file_name;
}

Text Force_field::
save_bond (Text const & force_field_name) const
{
    return force_field_name + "_1_bond.xls";
}

Text Force_field::
save_angle (Text const & force_field_name) const
{
    return force_field_name + "_2_angle.xls";
}

Text Force_field::
save_torsional (Text const & force_field_name) const
{
    return force_field_name + "_3_torsional.xls";
}

Text Force_field::
save_improper_torsions (Text const & force_field_name) const
{
    return force_field_name + "_4_improper_torsions.xls";
}

Text Force_field::
save_Van_der_Waals (Text const & force_field_name) const
{
    return force_field_name + "_5_Van_der_Waals.xls";
}

//-------------------- Bonds ----------------------------------------------
void  Force_field::
bond (int type1, int type2, double & Kr, double & Req) const
{
    if (type1 == Micro_object_type::none ||  //fix
        type1 == Micro_object_type::undefined_type    || 
        type1 == Micro_object_type::any_type          ||
        type2 == Micro_object_type::none ||
        type2 == Micro_object_type::undefined_type    || 
        type2 == Micro_object_type::any_type)
    {
        if (!Common_interactions::ad_hoc())
        {
            Text str = "No such bond parameter:  ";
            str += ::mot().itoa( type1 );   str += " ";
            str += ::mot().itoa( type2 );   
            
            if (show_user_warning_)
            {
                to_user().warning (str);
                show_user_warning_ = false;
            }
            //log () << "=== " << str << "\n";
            
            //error_file_ << str.c_str() << "\n";
            std::ofstream error_file ((Path::report()+"FF_error.txt").c_str(), 
                std::ios::out | std::ios::app);
            error_file << str.c_str() << std::endl;
        }

        ++problem_;

        Kr  = 0.;
        Req = 0.;
        success_ = false;
        return; 
    }

    Kr  = bond_par[type1][type2].Kr;
    Req = bond_par[type1][type2].R_eq;

    if (show_user_warning_ && Kr == 0.)
    {
        Text str = "Zero bond parameter:  ";
        str += ::mot().itoa( type1 );   str += " ";
        str += ::mot().itoa( type2 );   
        to_user().warning (str);
        show_user_warning_ = false;
    }
    
    /*log () << "bond()\n"
	   << "type1=" << type1
	   << "type2=" << type2
	   << "Kr="    << Kr
	   << "Req="   << Req
	   << "bond_par[type1][type2].Kr="   << bond_par[type1][type2].Kr
	   << "bond_par[type1][type2].R_eq=" << bond_par[type1][type2].R_eq
	   << "\n";//*/
#ifdef TESHI
    if (Kr == 0)
    {
        char buf[128];
        strcpy (buf, MOT.itoa (type1));
        strcat (buf, "-");
        strcat (buf, MOT.itoa (type2));
        strcat (buf, " - No such bond type");
        throw errorM (buf);
    }
#endif
    success_ = true;
}

void    Force_field::  
clear_bonds ()
{
    for (int i=0;  i<TYPE_N;  ++i)
    for (int j=0;  j<TYPE_N;  ++j)
    {
        bond_par[i][j].Kr   = 0.;
        bond_par[i][j].R_eq = 0.;
    }
}

void    Force_field::  
read_bonds (const char * bonds_file)
  {
    using namespace std;
    char     b1[5], b2[5];
    double     kr, req;

    File_ptr fp( bonds_file, "rt" );

    fgets( buf, BUF_SIZE, fp );
    fgets( buf, BUF_SIZE, fp );

#if TESDE
    cout << endl << "Bonds File " << bonds_file << endl;
#endif

    while( fgets( buf, BUF_SIZE, fp ) != NULL )
      {
        sscanf( buf, "%4s %4s %lg%lg", b1, b2, &kr, &req );
#if TESDE
        cout << b1 << "\t" << b2 << "\t" << kr << "\t" << req << endl;
        cout << ::mot().atoi(b1) << "\t" << ::mot().atoi(b2) << "\t" << kr << "\t" << req << endl;
#endif
        bond_par[::mot().atoi(b1)][::mot().atoi(b2)].Kr   = kr;
        bond_par[::mot().atoi(b1)][::mot().atoi(b2)].R_eq = req;
        bond_par[::mot().atoi(b2)][::mot().atoi(b1)].Kr   = kr;
        bond_par[::mot().atoi(b2)][::mot().atoi(b1)].R_eq = req;
    }
  }

//----------------- Angles ---------------------------------------------------

void    Force_field::
angle( int type1, int type2, int type3, double & K, double & angle_eq ) const
  {
    Array< int > key;
    /*key.add_by_default( 3 );
    key[0] = type1;
    key[1] = type2;
    key[2] = type3;//*/
    key.reserve( 3 );
    key.push_back( type1 );
    key.push_back( type2 );
    key.push_back( type3 );

    const Angles *result = Angle_par.get( key );
    if( !result)
    {
        if(!Common_interactions::ad_hoc())
        {
            Text str = "No such angle parameter:  ";
            str += ::mot().itoa( type1 );   str += " ";
            str += ::mot().itoa( type2 );   str += " ";
            str += ::mot().itoa( type3 );
            
            if (show_user_warning_)
            {
                to_user().warning (str);
                show_user_warning_ = false;
            }
            //log () << "=== " << str << "\n";

            //error_file_ << str.c_str() << "\n";
            std::ofstream error_file ((Path::report()+"FF_error.txt").c_str(), 
                std::ios::out | std::ios::app);
            error_file << str.c_str() << std::endl;
        }

        ++problem_;

        K = 0.0;
        angle_eq = 0.0;
        success_ = false;
        return;
    }
    K        = result->K;
    angle_eq = result->angle_eq;
    success_ = true;
  }

void    Force_field::
read_angles (const char * angles_file)
  {
    using namespace std;
    char     a1[5], a2[5], a3[5];
    double     k, angle_0;
    File_ptr fp( angles_file, "rt" );

    #if TESDE
    cout << endl << "Angles File " << angles_file << endl;
    #endif

    fgets( buf, BUF_SIZE, fp );
    fgets( buf, BUF_SIZE, fp );

    while( fgets( buf, BUF_SIZE, fp ) != NULL )
      {
        sscanf( buf, "%4s %4s %4s %lg%lg", a1, a2, a3, &k, &angle_0 );
    #if TESDE
    cout << a1 << "\t" << a2 << "\t" << a3 << "\t" << k << "\t" << angle_0 << endl;
    #endif
        angle_0 *= __Pi__ / 180;
        Array< int > key;
        /*key.add_by_default( 3 );

        key[0] = MOT.atoi( a1 );  key[1] = MOT.atoi( a2 );  key[2] = MOT.atoi( a3 );//*/
        key.reserve( 3 );
        key.push_back( ::mot().atoi( a1 ));
        key.push_back( ::mot().atoi( a2 ));
        key.push_back( ::mot().atoi( a3 ));
        Angle_par.set( key, Angles( k, angle_0 ) );

        if( key[0] == key[2] )  continue;

        key[0] = ::mot().atoi( a3 );  key[1] = ::mot().atoi( a2 );  key[2] = ::mot().atoi( a1 );
        Angle_par.set( key, Angles( k, angle_0 ) );
      }
  }
//----------------- Torsions -------------------------------------------------

const Array< Force_field::Torsions > *    Force_field::
torsion( int type1, int type2, int type3, int type4 ) const
  {
    Array< int > key;
    /*key.add_by_default( 4 );
    key[0] = type1;  key[1] = type2;  key[2] = type3;  key[3] = type4;//*/
    key.reserve( 4 );
    key.push_back( type1 );
    key.push_back( type2 );
    key.push_back( type3 );
    key.push_back( type4 );

    const Array< Torsions > *result = Torsion_par.get( key );
    if( result )  return result;
    else
      {
        key[0] = Micro_object_type::any_type;
        key[3] = Micro_object_type::any_type;
        result = Torsion_par.get( key );
        if( !result )
        {
            if(!Common_interactions::ad_hoc())
            {
                Text str = "No such torsion parameter:  ";
                str += ::mot().itoa( type1 );   str += " ";
                str += ::mot().itoa( type2 );   str += " ";
                str += ::mot().itoa( type3 );   str += " ";
                str += ::mot().itoa( type4 );
            
                if (show_user_warning_)
                {
                    to_user().warning( str );
                    show_user_warning_ = false;
                }
                //log () << "=== " << str << "\n";
                //error_file_ << str.c_str() << "\n";
                std::ofstream error_file ((Path::report()+"FF_error.txt").c_str(), 
                    std::ios::out | std::ios::app);
                error_file << str.c_str() << std::endl;
            }
            ++problem_;
            success_ = false;
            return 0;      
        }
      }
    success_ = true;
    return result;
  }

void    Force_field::
read_torsions (const char * torsions_file)
  {
    using namespace std;
    char     t1[5], t2[5], t3[5], t4[5];
    int      m, n;
    double     v, fi;
    File_ptr fp( torsions_file, "rt" );

    #if TESDE
    cout << endl << "Torsions File  " << torsions_file << endl;
    #endif

    fgets( buf, BUF_SIZE, fp );
    fgets( buf, BUF_SIZE, fp );

    while( fgets( buf, BUF_SIZE, fp ) != NULL )
      {
        sscanf( buf, "%4s %4s %4s %4s %i%lg%lg%i",
                      t1, t2, t3, t4, &m, &v, &fi, &n );
        #if TESDE
        cout << t1 << "\t" << t2 << "\t" << t3 << "\t"<< t4 << "\t"
             << m << "\t" << v << "\t" << fi << "\t" << n << endl;
        #endif

        if (fi != 0 && fi != 180)
        {
            Text message = "Warning in Force_field::read_torsions: "
                           "Fi0 != 0 && Fi0 != 180.";
        
            //fix to_user().status (message);
            //fix log() << "=== " << message << "\n";
        }
        fi *= __Pi__ / 180;

        Array< int >      key;
        //key.add_by_default( 4 );
        key.reserve( 4 );

        Array< Torsions > data;
        data.push_back( Torsions( m, v, fi, n ) );

        //key[0] = MOT.atoi( t1 );  key[1] = MOT.atoi( t2 );
        //key[2] = MOT.atoi( t3 );  key[3] = MOT.atoi( t4 );
        int id1 = ::mot().atoi (t1);  key.push_back (id1);
        int id2 = ::mot().atoi (t2);  key.push_back (id2);
        int id3 = ::mot().atoi (t3);  key.push_back (id3);
        int id4 = ::mot().atoi (t4);  key.push_back (id4);
        Torsion_par.add( key, data );

        if( key[0] == key[3] && key[1] == key[2] )  continue;

        key[0] = ::mot().atoi( t4 );  key[1] = ::mot().atoi( t3 );
        key[2] = ::mot().atoi( t2 );  key[3] = ::mot().atoi( t1 );
        Torsion_par.add( key, data );
      }
  }

//----------------- Improper Torsions --------------------------------------------

const Force_field::Imp_Tor *  Force_field::
imp_tor(int t1,  int t2, int t3, int t4,
        int& a1, int& a2,        int & a4 ) const
  {
    const Imp_Tor    *tmp = 0;

    for( int i=0;  i<Imp_tor_par.size();  i++ )
      {
        int T1 = Imp_tor_par[i].T1;
        int T2 = Imp_tor_par[i].T2;
        int T3 = Imp_tor_par[i].T3;
        int T4 = Imp_tor_par[i].T4;
        if( t3 == T3 )
          {
            if( ::mot().fit(t1,T1) && ::mot().fit(t2,T2) && ::mot().fit(t4,T4) )
              {
                a1 = 1;  a2 = 2; a4 = 4;
                tmp =  &(Imp_tor_par[i]);
                continue;
              }
            if(  ::mot().fit(t1,T1) && ::mot().fit(t2,T4) && ::mot().fit(t4,T2) )
              {
                a1 = 1;  a2 = 4; a4 = 2;
                tmp =  &(Imp_tor_par[i]);
                continue;
              }
            if(  ::mot().fit(t1,T2) && ::mot().fit(t2,T1) && ::mot().fit(t4,T4) )
              {
                a1 = 2;  a2 = 1; a4 = 4;
                tmp =  &(Imp_tor_par[i]);
                continue;
              }
            if(  ::mot().fit(t1,T2) && ::mot().fit(t2,T4) && ::mot().fit(t4,T1) )
              {
                a1 = 2;  a2 = 4; a4 = 1;
                tmp =  &(Imp_tor_par[i]);
                continue;
              }
            if(  ::mot().fit(t1,T4) && ::mot().fit(t2,T2) && ::mot().fit(t4,T1) )
              {
                a1 = 4;  a2 = 2; a4 = 1;
                tmp =  &(Imp_tor_par[i]);
                continue;
              }
            if(  ::mot().fit(t1,T4) && ::mot().fit(t2,T1) && ::mot().fit(t4,T2)  )
              {
                a1 = 4;  a2 = 1; a4 = 2;
                tmp =  &(Imp_tor_par[i]);
                continue;
              }
          }
      }
    //if( !tmp )  to_user().error( "Force_field::imp_tor" );
    return tmp;
  }

void Force_field::
read_imp_tor (const char * imp_tor_file)    //fix check all strings reading
  {
    using namespace std;
    char    t1[5], t2[5], t3[5], t4[5];
    int     n;
    double    v;
    int     size = 2;//fix ?

#if TESDE
    cout << endl << "Improper Torsions File  " << imp_tor_file << endl;
#endif

    {
      File_ptr fp( imp_tor_file, "rt" );
      while ( fgets( buf, BUF_SIZE, fp ) != NULL )  size++;
    }

    File_ptr fp( imp_tor_file, "rt" );
    //fgets( buf, BUF_SIZE, fp );
    //fix Imp_tor_par.add_by_default( size );
    for( int k=0;  k<size;  ++k ) Imp_tor_par.push_back( Imp_Tor() );

    fgets( buf, BUF_SIZE, fp );
    fgets( buf, BUF_SIZE, fp );

    for( int i=0;  i<size-1;  i++ )
      {
        fgets( buf, BUF_SIZE, fp );
        sscanf( buf, "%4s %4s %4s %4s %lg%*g%i",
                      t1, t2, t3, t4, &v, &n );
#if TESDE
        cout << t1 << "\t" << t2 << "\t" << t3 << "\t"<< t4 << "\t"
              << v << "\t" << n << endl;
#endif
        Imp_tor_par[i].T1      = ::mot().atoi( t1 );
        Imp_tor_par[i].T2      = ::mot().atoi( t2 );
        Imp_tor_par[i].T3      = ::mot().atoi( t3 );
        Imp_tor_par[i].T4      = ::mot().atoi( t4 );
        Imp_tor_par[i].N       = n;
        Imp_tor_par[i].half_V  = v;
      }
  }

//----------------- Van der Waals --------------------------------------------

double  Force_field::ad_hoc_R_star_ = 1.6;
double  Force_field::ad_hoc_Eps_ = 0.020001234;

void    Force_field::
clear_VdW ()
{
    for( int k=0;  k<TYPE_N;  ++k ) 
        R_eps.push_back( R_Eps() );

    //R_eps[0].R_star = 0.; 
    //R_eps[0].Eps    = 0.; 

    //for (int j=1;  j<TYPE_N;  j++)

    for (int j=0;  j<TYPE_N;  j++)
    {
        //R_eps[j].R_star = 0; 
        //R_eps[j].Eps = 0; 
        R_eps[j].R_star = ad_hoc_R_star_; 
        R_eps[j].Eps = ad_hoc_Eps_; 
    }
}

void    Force_field::
read_VdW (const char * VdW_file)
{
    using namespace std;
    int      i = 0, j;
    char     type[TYPE_N][5];
    double   r[TYPE_N], eps[TYPE_N];
    File_ptr fp( VdW_file, "rt" );

#if TESDE
    cout << endl << "Atom\tR_star\tEps  in file  " << VdW_file << endl;
#endif

    //fix R_eps.add_by_default( TYPE_N );
    //for( int k=0;  k<TYPE_N;  ++k ) R_eps.push_back( R_Eps() );

    fgets( buf, BUF_SIZE, fp );
    fgets( buf, BUF_SIZE, fp );

    while( fgets( buf, BUF_SIZE, fp ) != NULL )
    {
        sscanf( buf, "%4s %lg%lg", type[i], &r[i], &eps[i] );
#if TESDE
        cout << type[i] << "\t" << r[i] << "\t" << eps[i] << endl;
#endif
        i++;
    }
    //for( j=0;  j<TYPE_N;  j++ ) { R_eps[j].R_star = 0; R_eps[j].Eps = 0; }
    for( j=0;  j<i;  j++ )
      {
        R_eps[::mot().atoi(type[j])].R_star = r[j];
        R_eps[::mot().atoi(type[j])].Eps    = eps[j];
      }

    if (sigma_rule_ == "arithmetical" && epsilon_rule_ == "geometrical")
    {
        for(i=0;  i< ::mot().get_counter();  ++i)
        for(j=0;  j< ::mot().get_counter();  ++j)
        {
            cu_VdW[i][j].A = VdW_par[i][j].A =
                pow(  R_eps[i].R_star + R_eps[j].R_star, 12 )*
                sqrt( R_eps[i].Eps*R_eps[j].Eps );
            cu_VdW[i][j].B = VdW_par[i][j].B =
                2*pow(  R_eps[i].R_star + R_eps[j].R_star, 6 )*
                sqrt( R_eps[i].Eps*R_eps[j].Eps );
            cu_VdW[j][i].A = VdW_par[j][i].A =
                pow(  R_eps[i].R_star + R_eps[j].R_star, 12 )*
                sqrt( R_eps[i].Eps*R_eps[j].Eps );
            cu_VdW[j][i].B = VdW_par[j][i].B =
                2*pow(  R_eps[i].R_star + R_eps[j].R_star, 6 )*
                sqrt( R_eps[i].Eps*R_eps[j].Eps );
        }
    }
    else if (sigma_rule_ == "geometrical" && epsilon_rule_ == "geometrical")
    {
        for(i=0;  i< ::mot().get_counter();  ++i)
        for(j=0;  j< ::mot().get_counter();  ++j)
        {
            double Aii =     pow (2. * R_eps[i].R_star, 12) * R_eps[i].Eps;
            double Ajj =     pow (2. * R_eps[j].R_star, 12) * R_eps[j].Eps;
            double Bii = 2 * pow (2. * R_eps[i].R_star, 6 ) * R_eps[i].Eps;
            double Bjj = 2 * pow (2. * R_eps[j].R_star, 6 ) * R_eps[j].Eps;
            double Aij = sqrt(Aii * Ajj); 
            double Bij = sqrt(Bii * Bjj); 

            cu_VdW[i][j].A = VdW_par[i][j].A = Aij;
            cu_VdW[i][j].B = VdW_par[i][j].B = Bij;
            cu_VdW[j][i].A = VdW_par[j][i].A = Aij;
            cu_VdW[j][i].B = VdW_par[j][i].B = Bij;
        }
    }
    else if (sigma_rule_   == "Waldman-Hagler" && 
             epsilon_rule_ == "Waldman-Hagler")
    {
        double R2sig = 1. / pow (2, 1./6.);

        //int HW  = ::mot().atoi("HW" );
        //int OW  = ::mot().atoi("OW" );
        //int H4C = ::mot().atoi("H4C");
        //int H3C = ::mot().atoi("H3C");
        //int HC  = ::mot().atoi("HC" );

        for(i=0;  i< ::mot().get_counter();  ++i)
        for(j=0;  j< ::mot().get_counter();  ++j)
        {
     //   A6 = -4.d0*( 2.*(((S1**3)*(S2**3))/(S1**6+S2**6))*sqrt(E1*E2) )*
     //*       ( (0.5*(S1**6+S2**6)) )
     //   B12 = 4.d0*( 2.*(((S1**3)*(S2**3))/(S1**6+S2**6))*sqrt(E1*E2) )*
     //*       ( (0.5*(S1**6+S2**6)) )**2

            double sigma_ii = R2sig * R_eps[i].R_star * 2.;
            double sigma_jj = R2sig * R_eps[j].R_star * 2.;
            double eps_ii   = R_eps[i].Eps;
            double eps_jj   = R_eps[j].Eps;

            //==========================================================
            // fix
            //if (i == HW && j == OW)
            //{
            //    sigma_ii = R2sig * 0.2 * 2.;
            //    sigma_jj = R2sig * R_eps[j].R_star * 2.;
            //    eps_ii   = 0.01;
            //    eps_jj   = R_eps[j].Eps;
            //}
            //if (j == HW && i == OW)
            //{
            //    sigma_ii = R2sig * R_eps[i].R_star * 2.;
            //    sigma_jj = R2sig * 0.2 * 2.;
            //    eps_ii   = R_eps[i].Eps;
            //    eps_jj   = 0.01;
            //}
            ////OW	1.5	    0.1553	
            ////HW	1.20	0.01	
            //// crush
            //if (i == OW && (j == H4C || j == H3C || j == HC))
            //{
            //    sigma_ii = R2sig * 1.6 * 2.;
            //    sigma_jj = R2sig * R_eps[j].R_star * 2.;
            //    eps_ii   = 0.01;
            //    eps_jj   = R_eps[j].Eps;
            //}
            //if (j == OW && (i == H4C || i == H3C || i == HC))
            //{
            //    sigma_ii = R2sig * R_eps[i].R_star * 2.;
            //    sigma_jj = R2sig * 1.6 * 2.;
            //    eps_ii   = R_eps[i].Eps;
            //    eps_jj   = 0.01;
            //}
            //==========================================================

#ifndef NDEBUG
            double sigma_ij = 
                pow ((pow(sigma_ii, 6) + pow(sigma_jj, 6))/2., 1./6.);
#endif
            double sigma_ij_6 = (pow(sigma_ii, 6) + pow(sigma_jj, 6))/2.;
            double eps_ij   = 2. * ((pow(sigma_ii, 3) * pow(sigma_jj, 3))/
                (pow(sigma_ii, 6)+pow(sigma_jj, 6))) * sqrt (eps_ii*eps_jj);

            double Aij = 4. * eps_ij * sigma_ij_6*sigma_ij_6; 
            double Bij = 4. * eps_ij * sigma_ij_6; 

            cu_VdW[i][j].A = VdW_par[i][j].A = Aij;
            cu_VdW[i][j].B = VdW_par[i][j].B = Bij;
            cu_VdW[j][i].A = VdW_par[j][i].A = Aij;
            cu_VdW[j][i].B = VdW_par[j][i].B = Bij;
        }
    }
    else
        to_user().fatal_error ("\""+sigma_rule_ + " and " + "\""+epsilon_rule_+
            "\" are unknown sigma-epsilon rules.");
    //{    
    //    for(i=0;  i< ::mot().get_counter();  ++i)
    //    for(j=0;  j< ::mot().get_counter();  ++j)
    //    {
    //        cu_VdW[i][j].A = VdW_par[i][j].A =
    //            pow(  R_eps[i].R_star + R_eps[j].R_star, 12 )*
    //            sqrt( R_eps[i].Eps*R_eps[j].Eps );
    //        cu_VdW[i][j].B = VdW_par[i][j].B =
    //            2*pow(  R_eps[i].R_star + R_eps[j].R_star, 6 )*
    //            sqrt( R_eps[i].Eps*R_eps[j].Eps );
    //        cu_VdW[j][i].A = VdW_par[j][i].A =
    //            pow(  R_eps[i].R_star + R_eps[j].R_star, 12 )*
    //            sqrt( R_eps[i].Eps*R_eps[j].Eps );
    //        cu_VdW[j][i].B = VdW_par[j][i].B =
    //            2*pow(  R_eps[i].R_star + R_eps[j].R_star, 6 )*
    //            sqrt( R_eps[i].Eps*R_eps[j].Eps );
    //    }
    //}
}

void    Force_field::
read_VdW_exceptions (const char * VdW_file)
{
    using namespace std;
    char     type1[5], type2[5];
    double   A, B;
    File_ptr fp( VdW_file, "rt" );

    fgets (buf, BUF_SIZE, fp);
    fgets (buf, BUF_SIZE, fp);

    while (fgets( buf, BUF_SIZE, fp ) != NULL)
    {
        sscanf( buf, "%4s%4s %lg%lg", type1, type2, &A, &B );

        int i = ::mot().atoi(type1);
        int j = ::mot().atoi(type2);

        cu_VdW[i][j].A = VdW_par[i][j].A = A;
        cu_VdW[i][j].B = VdW_par[i][j].B = B;
        cu_VdW[j][i].A = VdW_par[j][i].A = A;
        cu_VdW[j][i].B = VdW_par[j][i].B = B;
    }
}

//---------------- Partial charges ------------------

void    Force_field::                    
load_partial_charges (Text const & file)
{
    using namespace std;

    ifstream in (file.c_str());
    if (in.fail())
        FLAW (Text("Cannot open file ") + file);

    string   line;
    getline (in, line);
    getline (in, line);

    while (getline (in, line))
    {
        if (line[0] = '#')
            continue;

        istringstream   in_line (line);
        string          type, reference;
        double          charge;
        in_line >> type >> charge >> reference;
        partial_charge_[::mot().atoi(type.c_str())] = Partial_charge (charge);
    }
}

void    Force_field::                    
read_induction  (Text const & file)
{
    using namespace std;

    ifstream in (file.c_str());
    if (in.fail())
        FLAW (Text("Cannot open file ") + file);

    string   line;
    getline (in, line);
    getline (in, line);

    while (getline (in, line))
    {
        if (line[0] = '#')
            continue;

        istringstream   in_line (line);
        string          type;
        double          iq;
        in_line >> type >> iq;
        induction_charge_[::mot().atoi(type.c_str())] = Induction_charge(iq);
    }
}

Text    Force_field::                    
save_partial_charges (Text const & force_field_name) const
{
    using namespace std;

    Text file_name;
    file_name += force_field_name;
    file_name += "_6_partial_charges.xls";
    ofstream out ((Path::force_field() + file_name).c_str());
    if (!out)
        FLAW (Text() + "Can not open file " + file_name);


    out << "string \tdouble \tstring   " << endl;
    out << "Type   \tCharge \tReference" << endl;

    for (std::map <int, Partial_charge>::const_iterator it=partial_charge_.begin();
         it != partial_charge_.end();  ++it)
    {
        int mm_type = it->first;
        Partial_charge const & charge = it->second;

        out << ::mot().itoa(mm_type) << "\t";
        if (charge.was_defined())
            out << charge.value();
            
        out  << "\t" << endl;
    }

    return file_name;
}

//---------------- Improper torsions ------------------

Improper_torsion::
Improper_torsion()
  : atom1( 0 ), atom2( 0 ), _atom3_( 0 ), atom4( 0 ),
    N( 0 ), half_V( 0 )   { }

Improper_torsion::
Improper_torsion(   int a1, int a2, int _a3_, int a4,
                    int n, double h )
  : atom1( a1 ), atom2( a2 ), _atom3_( _a3_ ), atom4( a4 ),
    N( n ), half_V( h ) { }

std::ostream &
operator<<( std::ostream & out, const Improper_torsion & iT )
  {
    using namespace std;
    return  out << "iTor "
                << setw( 4 ) << iT.atom1    << " "
                << setw( 4 ) << iT.atom2    << " "
                << setw( 4 ) << iT._atom3_  << " "
                << setw( 4 ) << iT.atom4    << " "
                << setw( 4 ) << iT.N        << " "
                << setw( 6 ) << iT.half_V;
  }

}//MM


