#include "Test.h"

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

#include <iostream>
#include "MolMeccano_file.h"
//#include "Log.h"
#include "User.h"

#include "Timer.h"
#include "Atom.h"
#include "Model.h"
#include "Model_kit.h"
#include "Path.h"
#include "Create.h"
#include "Force_field.h"
#include "Rough_hydrophobic_interactions.h"

#include "Mach_eps.h"

// Radizicka, A. and Wolfenden, R. 1988. Biochemistry, v.27. pp1664-1670

namespace MM
{

class Rough_hydrophobic_interactions_test : public Test
{
    Rough_hydrophobic_interactions *link;

    //Force_field     *force_field_;

public:
    Rough_hydrophobic_interactions_test (const Text & suite_name) 
        : Test (suite_name) 
    {
    }
    
    void run()
    {

        //force_field_ = new Force_field ("AMBER94");//fix
        test_2_atoms ();
        test_AMBER94 ();        
        //delete force_field_;

        //force_field_ = new Force_field ("OPLS");//fix
        test_2_atoms ();
        test_OPLS    ();        
        //delete force_field_;
    }

    void test_2_atoms ()
    {
        Own <Model>     model (prototype <Model> (). clone ());
        model().kit().set_force_field ("AMBER94");
        Rough_hydrophobic_interactions 
            interactions (model()/*, *force_field_*/);
        model(). load (Path::tests() + "HH.hin" );

        //----------------------------------------------------------
        Atom & atom = model().atom(0);
        double x = atom.x();
        //double charge = atom.charge();
        atom.set_charge (1.);
        model().atom(1).set_charge (1.);

        // Gradient ------------------------------
        double fx;
        interactions.add_force ();
        fx = atom.kit().potential_force().x();
        //VdW -0.072293126800337879
        //DOUBLES_EQUAL (0.61881365438620939,   fx, 0.000001);
        //DOUBLES_EQUAL (-27.376505282579391,   fx, 0.00001);

        // finite difference Gradient ------------------------------
        double E = interactions.potential ();
        atom.set_x (x + d_mach_eps_2);
        //interactions.cache_atoms ();
        double F_x = (E - interactions.potential ()) / d_mach_eps_2;
        atom.set_x (x);
        DOUBLES_EQUAL (F_x,   fx, 0.00001);

        //----------------------------------------------------------
    }

    void test_AMBER94 ()
    {
        //MolMeccano_file::apply_torsions (false);

        Own <Model>     model (prototype <Model> (). clone ());
        model().kit().set_force_field ("AMBER94");
        //Common_interactions 
        Rough_hydrophobic_interactions 
            //interactions (model(), model(), *force_field_);
            interactions (model()/*, *force_field_*/);
        
        double E;
        std::cout << "\n\n####### AMBER94 ##############\n";
        std::cout <<     "res chx-wat oct-wat still   rough\n";

                                                        // Radzicka & Wolfenden
        model(). load (Path::aminoacids() + "Ala.mlm" );//  1.81 Alanine 
        E = interactions.potential ();
        std::cout << "Ala   1.81   0.52  -12.3860 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Arg+.mlm");//-14.92 Arginine 
        E = interactions.potential ();
        std::cout << "Arg+ -14.92 -1.32  -75.5872 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Asn.mlm" );// -6.64 Asparagine 
        E = interactions.potential ();
        std::cout << "Asn  -6.64  -0.01  -27.3687 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Asp-.mlm");// -8.72 Aspartic acid 
        E = interactions.potential ();
        std::cout << "Asp- -8.72   ----  -79.7101 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Cys.mlm" );//  1.28 Cysteine 
        E = interactions.potential ();
        std::cout << "Cys   1.28   ----  -13.4418 " 
            << interactions.last_hydrophobic_potential() << "\n";


        model().load (Path::aminoacids() + "Gln.mlm" );// -5.54 Glutamine 
        E = interactions.potential ();
        std::cout << "Gln  -5.54  -0.07  -24.9896 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Glu-.mlm");// -6.81 Glutamic acid 
        E = interactions.potential ();
        std::cout << "Glu- -6.81  -0.79  -88.1826 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Gly.mlm" );//  0.94 Glycine 
        E = interactions.potential ();
        std::cout << "Gly   0.94   0.00  -13.2011 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "His+.mlm");// -4.66 Histidine 
        E = interactions.potential ();
        std::cout << "His+ -4.66   0.95  -74.4488 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Ile.mlm" );//  4.92 Isoleucine 
        E = interactions.potential ();
        std::cout << "Ile   4.92   2.04  -10.9896 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Leu.mlm" );//  4.92 Leucine 
        E = interactions.potential ();
        std::cout << "Leu   4.92   1.76  -10.7380 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Lys+.mlm");// -5.55 Lysine 
        E = interactions.potential ();
        std::cout << "Lys+ -5.55   0.08  -84.4491 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Met.mlm" );//  2.35 Methionine 
        E = interactions.potential ();
        std::cout << "Met   2.35   1.32  -14.1034 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Phe.mlm" );//  2.98 Phenylalanine 
        E = interactions.potential ();
        std::cout << "Phe   2.98   2.09  -11.4356 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Pro.mlm" );//       Proline 
        E = interactions.potential ();
        std::cout << "Pro   ---    ----  -11.7816 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Ser.mlm" );// -3.40 Serine 
        E = interactions.potential ();
        std::cout << "Ser  -3.40   0.04  -16.8811 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Thr.mlm" );// -2.57 Threonine 
        E = interactions.potential ();
        std::cout << "Thr  -2.57   0.27  -15.0150 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Trp.mlm" );//  2.33 Tryptophan 
        E = interactions.potential ();
        std::cout << "Trp   2.33   2.51  -14.8535 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Tyr.mlm" );// -0.14 Tyrosine 
        E = interactions.potential ();
        std::cout << "Tyr  -0.14   1.63  -16.2236 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Val.mlm" );//  4.04 Valine 
        E = interactions.potential ();
        std::cout << "Val   4.04   1.18  -11.7038 " 
            << interactions.last_hydrophobic_potential() << "\n";

        std::cout << "##############################\n";
        MolMeccano_file::apply_torsions (true);
    }
    void test_OPLS ()
    {
        //MolMeccano_file::apply_torsions (false);

        Own <Model>     model (prototype <Model> (). clone ());
        model().kit().set_force_field ("OPLS");

        //Common_interactions 
        Rough_hydrophobic_interactions 
            //interactions (model(), model(), *force_field_);
            interactions (model()/*, *force_field_*/);

        double E;
        std::cout << "\n\n########  OPLS  ##############\n";
        std::cout <<     "res chx-wat oct-wat still   rough\n";

                                                        // Radzicka & Wolfenden
        model(). load (Path::aminoacids() + "Ala.mlm" );//  1.81 Alanine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Ala   1.81   0.52  -13.0805 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Arg+.mlm");//-14.92 Arginine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Arg+ -14.92 -1.32  -80.5023 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Asn.mlm" );// -6.64 Asparagine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Asn  -6.64  -0.01  -28.4528 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Asp-.mlm");// -8.72 Aspartic acid 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Asp- -8.72   ----  -84.6869 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Cys.mlm" );//  1.28 Cysteine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Cys   1.28   ----  -15.2745 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Gln.mlm" );// -5.54 Glutamine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Gln  -5.54  -0.07  -24.3874 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Glu-.mlm");// -6.81 Glutamic acid 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Glu- -6.81  -0.79  -94.5898 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Gly.mlm" );//  0.94 Glycine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Gly   0.94   0.00  -13.7098 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "His+.mlm");// -4.66 Histidine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "His+ -4.66   0.95  -76.6452 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Ile.mlm" );//  4.92 Isoleucine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Ile   4.92   2.04  -12.0019 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Leu.mlm" );//  4.92 Leucine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Leu   4.92   1.76  -11.8656 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Lys+.mlm");// -5.55 Lysine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Lys+ -5.55   0.08  -95.8007 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Met.mlm" );//  2.35 Methionine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Met   2.35   1.32  -17.8222 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Phe.mlm" );//  2.98 Phenylalanine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Phe   2.98   2.09  -11.8131 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Pro.mlm" );//       Proline 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Pro   ---    ----  -12.7124 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Ser.mlm" );// -3.40 Serine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Ser  -3.40   0.04  -17.6971 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Thr.mlm" );// -2.57 Threonine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Thr  -2.57   0.27  -16.4848 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Trp.mlm" );//  2.33 Tryptophan 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Trp   2.33   2.51  -15.5688 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Tyr.mlm" );// -0.14 Tyrosine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Tyr  -0.14   1.63  -17.5330 " 
            << interactions.last_hydrophobic_potential() << "\n";

        model().load (Path::aminoacids() + "Val.mlm" );//  4.04 Valine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        std::cout << "Val   4.04   1.18  -12.5206 " 
            << interactions.last_hydrophobic_potential() << "\n";

        std::cout << "##############################\n";

        //----------------------------------------------------------
        Atom & atom = model().atom(3);
        double x = atom.x();
        //double charge = atom.charge();
        atom.set_charge (0.);

        // Gradient ------------------------------
        double fx/*, fy, fz*/;
        interactions.add_force ();
        fx = atom.kit().potential_force().x();
        //DOUBLES_EQUAL (-2.0851679991786898,   fx, 0.000001);

        // finite difference Gradient ------------------------------
        E = interactions.potential ();
        atom.set_x (x + d_mach_eps_2);
        //interactions.cache_atoms ();
        //double F_x = (E - interactions.potential ()) / d_mach_eps_2;
        atom.set_x (x);
        //DOUBLES_EQUAL (-2.0851679991786898,   F_x, 0.000001);

        //----------------------------------------------------------
        
        MolMeccano_file::apply_torsions (true);
    }
};

Rough_hydrophobic_interactions_test 
    test_rough_hydrophobic_interactions ("correctness");

}//MM
