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

#include "Test.h"

#include "Log.h"

#include "AMBER_energy.h"
#include "Element.h"
#include "Atom_impl.h"
#include "Arbitrary_atom_group.h"
#include "Arbitrary_bond_group.h"
#include "MOT.h"

#include "Angle.h"

#include "Bond_energy_CT_CT_amber.h"
#include "Angle_energy_CT_CT_HC_amber.h"
#include "Torsion_energy_HC_CT_CT_HC_amber.h"
#include "Van_der_Waals_energy_CT_CT_amber.h"

#include "AMBER_bond_energy.h"
#include "AMBER_angle_energy.h"
#include "AMBER_torsion_energy.h"
#include "AMBER_Van_der_Waals_energy.h"

#include "Array_of_ref.h"
#include "Array_of_own.h"
#include "Path.h"

#include "Default_bond.h"

#include "Dihedral_angle.h"

namespace MM
{

class Amber_energy_test : public Test
{
public:
    Amber_energy_test( const Text & suite_name ) : Test( suite_name ) { }

    void run ()
    {
        CT_CT_bond_test         ();     count_tests();
        CT_CT_HC_angle_test     ();     count_tests();
        HC_CT_CT_HC_torsion_test();     count_tests();
        Van_der_Waals_test      ();     count_tests();
        ethane_test             ();     count_tests();
    }

    void CT_CT_bond_test();
    void CT_CT_HC_angle_test();
    void HC_CT_CT_HC_torsion_test();
    void Van_der_Waals_test();
    void ethane_test();
};

Amber_energy_test test_amber_energy("correctness");

/**
E_bond = Sum ( K_r(r-r0)^2 )

B1  B2  KR     REQ    REFERENCE
CT  CT  310.0  1.526  "JCC,7,(1986),230;" AA, SUGARS
*/
void Amber_energy_test::
CT_CT_bond_test ()
{
    //Atom_impl atom1 (Element ("C"));    atom1.set_mm_type (MOT.atoi ("CT"));
    //Atom_impl atom2 (Element ("C"));    atom2.set_mm_type (MOT.atoi ("CT"));
    Atom_impl atom1 (Element ("C"));    atom1.kit().set_mm_type ("CT");
    Atom_impl atom2 (Element ("C"));    atom2.kit().set_mm_type ("CT");

    atom1.set_x (0);   atom1.set_y (0);   atom1.set_z (0);
    atom2.set_x (1);   atom2.set_y (0);   atom2.set_z (0);

    // 310.0 * (1.526 - 1)^2 = 85.76956
    Bond_energy_CT_CT_amber bond_energy (atom1, atom2);
    DOUBLES_EQUAL (85.76956, bond_energy.value (), 0.000001);

    //Length length (atom1, atom2);
    //DOUBLES_EQUAL (1., length.value (), 0.000001);
    DOUBLES_EQUAL (1., atom1.position(). distance(atom2.position()), 0.000001);

    Bond_length bond_length (atom1, atom2);
    DOUBLES_EQUAL (1., bond_length.value(), 0.000001);

    AMBER_bond_energy energy (atom1, atom2, 310.0, 1.526);
    DOUBLES_EQUAL (85.76956, energy.value(), 0.000001);

    //AMBER_bond_energy energy2 (atom1, atom2);
    //DOUBLES_EQUAL (85.76956, energy2.value(), 0.000001);
}

/**
E_bond_angle = Sum ( K_theta (theta-theta0)^2 )

A1  A2  A3  KTHETA  THETAEQ REFERENCE
CT  CT  HC  50.0    109.50  changed based on NMA nmodes
*/
void Amber_energy_test::
CT_CT_HC_angle_test ()
{
    //Atom_impl c1 (Element ("C"));   c1.set_mm_type (MOT.atoi ("CT"));
    //Atom_impl c2 (Element ("C"));   c2.set_mm_type (MOT.atoi ("CT"));
    //Atom_impl h1 (Element ("H"));   h1.set_mm_type (MOT.atoi ("HC"));

    Atom_impl c1 (Element ("C"));   c1.kit().set_mm_type ("CT");
    Atom_impl c2 (Element ("C"));   c2.kit().set_mm_type ("CT");
    Atom_impl h1 (Element ("H"));   h1.kit().set_mm_type ("HC");

    c1.set_x (0);  c1.set_y (0);  c1.set_z (0);
    c2.set_x (1);  c2.set_y (0);  c2.set_z (0);
    h1.set_x (1);  h1.set_y (1);  h1.set_z (0);

    // 50*(Pi*109.5/180 - Pi/2)^2 = 5.79153869369479
    Angle_energy_CT_CT_HC_amber angle_energy (c1, c2, h1);
    DOUBLES_EQUAL (5.79153869369479, angle_energy.value(), 0.000001);

    Angle angle (c1.position(), c2.position(), h1.position());
    DOUBLES_EQUAL (1.57079632679489, angle.value(), 0.000001);

    AMBER_angle_energy energy (c1, c2, h1, 50., 109.5 * 3.14159265358/180.);
    DOUBLES_EQUAL (5.79153869369479, energy.value(), 0.000001);

    //AMBER_angle_energy energy2( c1, c2, h1 );
    //DOUBLES_EQUAL (5.79153869369479, energy2.value(), 0.000001);
}

/**
E_dihedral = Sum Sum(F terms) ( V_n/2 (1 + cos(n phi - phi0)))

T1      T2      T3      T4      MULT    HALF_VN GAMMA   N      REFERENCE	
**      CT      CT      **      9       1.40    0.0     3      JCC,7,(1986),230
*/
void Amber_energy_test::
HC_CT_CT_HC_torsion_test ()
{
    //Atom_impl h1 (Element ("H"));   h1.set_mm_type (MOT.atoi ("HC"));
    //Atom_impl c1 (Element ("C"));   c1.set_mm_type (MOT.atoi ("CT"));
    //Atom_impl c2 (Element ("C"));   c2.set_mm_type (MOT.atoi ("CT"));
    //Atom_impl h2 (Element ("H"));   h2.set_mm_type (MOT.atoi ("HC"));

    Atom_impl h1 (Element ("H"));   h1.kit().set_mm_type ("HC");
    Atom_impl c1 (Element ("C"));   c1.kit().set_mm_type ("CT");
    Atom_impl c2 (Element ("C"));   c2.kit().set_mm_type ("CT");
    Atom_impl h2 (Element ("H"));   h2.kit().set_mm_type ("HC");

    h1.set_x (0);  h1.set_y (1);  h1.set_z (0);
    c1.set_x (0);  c1.set_y (0);  c1.set_z (0);
    c2.set_x (1);  c2.set_y (0);  c2.set_z (0);
    h2.set_x (1);  h2.set_y (1);  h2.set_z (1);

    // 1.40 * (1 + cos(3 * Pi/4 - 0)) / 9 = 0.0455611665950794
    Torsion_energy_HC_CT_CT_HC_amber torsion_energy (h1, c1, c2, h2);
    DOUBLES_EQUAL (0.0455611665950794, torsion_energy.value(), 0.000001);

    Dihedral_angle dihedral (h1.position(), c1.position(), 
                             c2.position(), h2.position());
    DOUBLES_EQUAL (3.14159265358 / 4., dihedral.value(), 0.000001);

    //fix restore
    //AMBER_torsion_energy energy (h1, c1, c2, h2);
    //DOUBLES_EQUAL (0.0455611665950794, energy.value(), 0.000001);
}

void Amber_energy_test::
Van_der_Waals_test ()
{

    //Atom_impl c1 (Element ("C"));   c1.set_mm_type (MOT.atoi ("CT"));
    //Atom_impl c2 (Element ("C"));   c2.set_mm_type (MOT.atoi ("CT"));

    Atom_impl c1 (Element ("C"));   c1.kit().set_mm_type ("CT");
    Atom_impl c2 (Element ("C"));   c2.kit().set_mm_type ("CT");

    c1.set_x (0);    c1.set_y (0);    c1.set_z (0);
    c2.set_x (1);    c2.set_y (0);    c2.set_z (0);

    //  Type        R*      Epsilon            Reference
    //  CT          1.9080  0.1094             Spellmeyer
    //
    // R(ij)* = R(i)* + R(j)*                       = 1.9080 + 1.9080       = 3.816
    //
    //                                          1/2
    // Epsilon(ij) = ( Epsilon(i) * Epsilon(j) )                            = 0.1094
    //
    //                                 12
    // A(ij) = Epsilon(ij) * ( R(ij)* )             = 0.1094 * 3.816^12     = 1043080.2307033
    //                                     6
    // B(ij) = 2 * Epsilon(ij) * ( R(ij)* )         = 2 * 0.1094 * 3.816^6  = 675.61224748799
    //

    //  E = A/1 - B/1 = 1043080.2307033 - 675.61224748799 = 1042404.61845581201
    Van_der_Waals_energy_CT_CT_amber Van_der_Waals_energy (c1, c2);
    DOUBLES_EQUAL (1042404.61845581201, Van_der_Waals_energy.value(), 0.000001);

    AMBER_Van_der_Waals_energy energy (c1, c2, 1043080.2307033, 675.61224748799);
    DOUBLES_EQUAL (1042404.61845581201, energy.value(), 0.000001);

    //AMBER_Van_der_Waals_energy energy2 (c1, c2);
    //DOUBLES_EQUAL (1042404.61845581201, energy2.value(), 0.000001);
}

void Amber_energy_test::
ethane_test()
{
    //atom 1 - C CT - 0 -0.7810696 -0.6815369 -0.02689989 4 2 s 3 s 4 s 5 s
    //atom 2 - C CT - 0 -0.7810696 0.8584631 -0.02689989 4 1 s 6 s 7 s 8 s
    //atom 3 - H HC - 0 0.006169164 -1.044866 -0.6874654 1 1 s

    //atom 4 - H HC - 0 -0.602613 -1.044866 0.98515 1 1 s

    //atom 5 - H HC - 0 -1.746755 -1.044874 -0.3783805 1 1 s
    //atom 6 - H HC - 0 -1.808733 1.221792 -0.02689989 1 2 s
    //atom 7 - H HC - 0 -0.2672464 1.221792 0.8630875 1 2 s
    //atom 8 - H HC - 0 -0.2672352 1.2218 -0.9168776 1 2 s
    //ChemCon            E 2.42727779173

    //HyperChem 
    //Total Energy=2.427278 Gradient=2.624787.
    //Bond=0.06076   Angle=0.000128842   Dihedral=2.10002 Vdw=0.266372 

    /*Atom_impl c1 (Element ("C"));   c1.set_mm_type (MOT.atoi ("CT"));
    Atom_impl c2 (Element ("C"));   c2.set_mm_type (MOT.atoi ("CT"));
    Atom_impl h1 (Element ("H"));   h1.set_mm_type (MOT.atoi ("HC"));
    Atom_impl h2 (Element ("H"));   h2.set_mm_type (MOT.atoi ("HC"));
    Atom_impl h3 (Element ("H"));   h3.set_mm_type (MOT.atoi ("HC"));
    Atom_impl h4 (Element ("H"));   h4.set_mm_type (MOT.atoi ("HC"));
    Atom_impl h5 (Element ("H"));   h5.set_mm_type (MOT.atoi ("HC"));
    Atom_impl h6 (Element ("H"));   h6.set_mm_type (MOT.atoi ("HC"));//*/

    Atom_impl c1 (Element ("C"));   c1.kit().set_mm_type ("CT");
    Atom_impl c2 (Element ("C"));   c2.kit().set_mm_type ("CT");
    Atom_impl h1 (Element ("H"));   h1.kit().set_mm_type ("HC");
    Atom_impl h2 (Element ("H"));   h2.kit().set_mm_type ("HC");
    Atom_impl h3 (Element ("H"));   h3.kit().set_mm_type ("HC");
    Atom_impl h4 (Element ("H"));   h4.kit().set_mm_type ("HC");
    Atom_impl h5 (Element ("H"));   h5.kit().set_mm_type ("HC");
    Atom_impl h6 (Element ("H"));   h6.kit().set_mm_type ("HC");

    c1.set_x (-0.7810696);     c1.set_y (-0.6815369);  c1.set_z (-0.02689989);
    c2.set_x (-0.7810696);     c2.set_y (0.8584631);   c2.set_z (-0.02689989);

    h1.set_x (0.006169164);    h1.set_y (-1.044866);   h1.set_z (-0.6874654);
    h2.set_x (-0.602613);      h2.set_y (-1.044866);   h2.set_z (0.98515);
    h3.set_x (-1.746755);      h3.set_y (-1.044874);   h3.set_z (-0.3783805);
    h4.set_x (-1.808733);      h4.set_y (1.221792);    h4.set_z (-0.02689989);
    h5.set_x (-0.2672464);     h5.set_y (1.221792);    h5.set_z (0.8630875);
    h6.set_x (-0.2672352);     h6.set_y (1.2218);      h6.set_z (-0.9168776);

    Arbitrary_atom_group atoms;
    atoms.add (c1);    atoms.add (c2);    
    atoms.add (h1);    atoms.add (h2);
    atoms.add (h3);    atoms.add (h4);    
    atoms.add (h5);    atoms.add (h6);

    Default_bond cc  (c1, c2);
    Default_bond ch1 (c1, h1);
    Default_bond ch2 (c1, h2);
    Default_bond ch3 (c1, h3);
    Default_bond ch4 (c2, h4);
    Default_bond ch5 (c2, h5);
    Default_bond ch6 (c2, h6);
    Arbitrary_bond_group bonds;
    bonds.add (cc );
    bonds.add (ch1);
    bonds.add (ch2);
    bonds.add (ch3);
    bonds.add (ch4);
    bonds.add (ch5);
    bonds.add (ch6);

    Force_field *force_field = new Force_field ("AMBER94"); 

    Own <Force_field> janitor;

    AMBER_energy energy (atoms, bonds, *force_field);
    DOUBLES_EQUAL (0.06076,     energy.bond_energy(),          0.00001 );
    DOUBLES_EQUAL (0.000128842, energy.angle_energy(),         0.0000001);
    DOUBLES_EQUAL (2.10002,     energy.torsion_energy(),       0.00001 );
    DOUBLES_EQUAL (0.266372,    energy.van_der_Waals_energy(), 0.00001 );
    DOUBLES_EQUAL (2.427278,    energy.value(),                0.00001 );
}

}//MM

