#include "Test.h"

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

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

#include "Mach_eps.h"
#include "Timer.h"
#include "Atom.h"
#include "Atom_kit.h"
#include "Model.h"
#include "Model_kit.h"
#include "Path.h"
#include "Create.h"
#include "Force_field.h"
#include "Common_interactions.h"
#include "Conventional_interactions.h"

//Arg+
//         1  ACE ARG NME

//        tinker.key  
//parameters      ../params/oplsua
//dielectric               1.0
//vdw-14-scale    8.0                       default 2.0 !!!
//SOLVATE STILL

//              **   57    56    **     1 2     0   3  X-C2-C2-X   
//    20        9    10    11    12       2.000 180/3


//Ser
//24	23	553	0.96	HO-OH   	SUG(OL)
//    11       10    11                     553.0000    0.9450 23-24

//Thr
//    12       10    12                     553.0000    0.9450
//0      10     12  0.9600   2313.752  23-24 


namespace MM
{

class OPLS_interactions_test : public Test
{
    Common_interactions *link;

    Force_field     *force_field_;

public:
    OPLS_interactions_test (const Text & suite_name) 
        : Test (suite_name) 
    {
    }
    
    void run()
    {
        force_field_ = new Force_field ("OPLS");//fix
        
        test_20_acids         ();       count_tests();
        test_HNCCO            ();       count_tests();

        delete force_field_;
    }

    void test_20_acids ()
    {
        MolMeccano_file::apply_torsions (false);

        Own <Model>     model (prototype <Model> (). clone ());

        //Own <Common_interactions>   janitor (Common_interactions::
        //    create(model(), model(), *force_field_));
        model().kit().set_force_field ("OPLS");
        Own <Common_interactions>   janitor 
            (new Conventional_interactions (model()/*, model(), *force_field_*/));

        Common_interactions & interactions = janitor();

        double E;
                                                        // Radzicka & Wolfenden
        model(). load (Path::aminoacids() + "Ala.mlm" );//  1.81 Alanine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.4274, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.5420, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.1707, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  0.0909, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-46.4867, 0.1);

        model().load (Path::aminoacids() + "Arg+.mlm");//-14.92 Arginine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.9476, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  2.1194, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        //fix DOUBLES_EQUAL (interactions.last_torsion_potential          (),  4.3439, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  2.3168, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-46.2368, 0.1);

        model().load (Path::aminoacids() + "Asn.mlm" );// -6.64 Asparagine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.6183, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  2.4240, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.8054, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  3.0923, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (), 69.3093, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-56.5737, 0.1);

        model().load (Path::aminoacids() + "Asp-.mlm");// -8.72 Aspartic acid 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.4731, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.9490, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3439, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (), 43.0555, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-61.1498, 0.1);

        model().load (Path::aminoacids() + "Cys.mlm" );//  1.28 Cysteine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();
        //model().save_as ("_.mmol");

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.4540, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.5454, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3439, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  0.4772, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-40.7842, 0.1);


        model().load (Path::aminoacids() + "Gln.mlm" );// -5.54 Glutamine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.6693, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  2.5853, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.8054, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  3.0923, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  5.1153, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-75.4332, 0.1);

        model().load (Path::aminoacids() + "Glu-.mlm");// -6.81 Glutamic acid 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.5240, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  2.1104, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3439, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  3.0243, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-52.5289, 0.1);

        model().load (Path::aminoacids() + "Gly.mlm" );//  0.94 Glycine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.3765, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.3310, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0000, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3414, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  0.2706, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-46.4867, 0.1);

        model().load (Path::aminoacids() + "His+.mlm");// -4.66 Histidine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        //DOUBLES_EQUAL (interactions.last_bond_potential             (),  2.0154, 0.0001);
        //DOUBLES_EQUAL (interactions.last_angle_potential            (),  3.9668, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0378, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3457, 0.0001);
        //DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (), -0.1204, 0.0001);
        //DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-20.1574, 0.1);

        model().load (Path::aminoacids() + "Ile.mlm" );//  4.92 Isoleucine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.5803, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.7748, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0379, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3460, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  6.8531, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-46.4867, 0.1);

        model().load (Path::aminoacids() + "Leu.mlm" );//  4.92 Leucine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.5803, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.9425, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0102, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3439, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),151.4232, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-46.4867, 0.1);

        model().load (Path::aminoacids() + "Lys+.mlm");// -5.55 Lysine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  1.0192, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  2.0979, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3439, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  1.2937, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-44.4699, 0.1);

        model().load (Path::aminoacids() + "Met.mlm" );//  2.35 Methionine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.5228, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  2.3269, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3439, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  0.9296, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-47.7307, 0.1);

        model().load (Path::aminoacids() + "Phe.mlm" );//  2.98 Phenylalanine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        //DOUBLES_EQUAL (interactions.last_bond_potential             (),  1.4745, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.9316, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3439, 0.0001);
        //DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  2.3958, 0.0001);
        //DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-42.6353, 0.1);

        model().load (Path::aminoacids() + "Pro.mlm" );//       Proline 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.7387, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  18.4162, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  1.8630, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  20.3014, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  14.4979, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-37.6723, 0.1);


        model().load (Path::aminoacids() + "Ser.mlm" );// -3.40 Serine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        //DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.5281, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.6833, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3439, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  0.2077, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-36.5271, 0.1);

        model().load (Path::aminoacids() + "Thr.mlm" );// -2.57 Threonine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        //DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.5791, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.6792, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.2421, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3499, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  4.2630, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-36.6868, 0.1);

        model().load (Path::aminoacids() + "Trp.mlm" );//  2.33 Tryptophan 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        //DOUBLES_EQUAL (interactions.last_bond_potential             (),  15.2620, 0.0001);
        //DOUBLES_EQUAL (interactions.last_angle_potential            (),  4.5062, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0378, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3457, 0.0001);
        //DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  2.2679, 0.0001);
        //DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-49.2252, 0.1);

        model().load (Path::aminoacids() + "Tyr.mlm" );// -0.14 Tyrosine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        //DOUBLES_EQUAL (interactions.last_bond_potential             (),  1.6806, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  2.1981, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0101, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3439, 0.0001);
        //DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  2.2434, 0.0001);
        //DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-48.1473, 0.1);

        model().load (Path::aminoacids() + "Val.mlm" );//  4.04 Valine 
        model().kit().set_force_field ("OPLS");
        model().kit().strip ();
        E = interactions.potential ();

        DOUBLES_EQUAL (interactions.last_bond_potential             (),  0.5294, 0.0001);
        DOUBLES_EQUAL (interactions.last_angle_potential            (),  1.7732, 0.0001);
        //DOUBLES_EQUAL (interactions.last_improper_torsion_potential (),  0.0102, 0.0001);
        DOUBLES_EQUAL (interactions.last_torsion_potential          (),  0.3437, 0.0001);
        DOUBLES_EQUAL (interactions.last_Van_der_Waals_potential    (),  4.7015, 0.0001);
        DOUBLES_EQUAL (interactions.last_electrostatic_potential    (),-46.4867, 0.1);

        MolMeccano_file::apply_torsions (true);
    }

    void test_HNCCO ()
    {
        Own <Model>     model (prototype <Model> (). clone ());

        Conventional_interactions interactions (model()/*, model(), *force_field_*/);

        model().load (Path::tests() + "HNCCO.hin");
        model().kit().set_force_field ("OPLS");

        // Gradient ------------------------------
        double fx, fy, fz;
        interactions.add_force ();
        //interactions.flush_force ();
        fx = model().atom(0).kit().potential_force().x();
        fy = model().atom(0).kit().potential_force().y();
        fz = model().atom(0).kit().potential_force().z();

        DOUBLES_EQUAL (-0.0055321397183183407,     fx, 0.000001);   //-1.7063697844441057
        DOUBLES_EQUAL (-0.0022478062046809022,     fy, 0.000001);
        DOUBLES_EQUAL (-0.0022342843810867485,     fz, 0.000001);

        // finite difference Gradient ------------------------------
        double P = interactions.potential ();
        
        double x = model().atom(0).x();
        model().atom(0).set_x (x + d_mach_eps_2);
        interactions.cache_atoms ();

        double F_x = (P - interactions.potential ()) / d_mach_eps_2;
        model().atom(0).set_x (x);
        DOUBLES_EQUAL (-0.0055321397183183407, F_x, 0.0001);
    }

};

OPLS_interactions_test  test_OPLS_interactions ("correctness");

}//MM
