#include "Test.h"

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

#include "Mach_eps.h"
#include "Timer.h"
#include "Lock.h"
#include "Atom_impl.h"
#include "Model.h"
#include "Path.h"
#include "Create.h"
#include "Force_field.h"

#include "Conventional_interactions.h"
#include "Commonplace_interactions.h"

#include "Arbitrary_atom_group.h"
#include "Arbitrary_bond_group.h"
#include "Default_bond.h"

namespace MM
{

class Covalent_interactions_test : public Test
{
    Conventional_interactions *link;

    Force_field     *force_field_;

public:
    Covalent_interactions_test (const Text & suite_name) 
        : Test (suite_name) 
    {
    }

    void run()
    {
        force_field_ = new Force_field ("AMBER94");     //fix

        //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();

        delete force_field_;
    }

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

CT-CT  310.0    1.526       JCC,7,(1986),230; AA, SUGARS
*/
    //void CT_CT_bond_test ()
    //{
    //    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);

    //    Arbitrary_atom_group atom_group;
    //    atom_group.add (atom1);
    //    atom_group.add (atom2);

    //    Default_bond bond (atom1, atom2);

    //    Arbitrary_bond_group bond_group;
    //    bond_group.add (bond);

    //    //Atom_group_cache        nb (atom_group);
    //    //Lock <Atom_group_cache> lock(nb);
    //    //Common_interactions  nb (atom_group, *force_field_);
    //    //Lock <Common_interactions> lock(nb);
    //    //Common_interactions  interactions (bond_group, *force_field_, nb);
    //    Conventional_interactions interactions 
    //        (atom_group, bond_group, *force_field_);

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

    //    // Potential -----------------------------
    //    // 310.0 * (1.526 - 1)^2 = 85.76956
    //    double P = interactions.potential ();
    //    DOUBLES_EQUAL (85.76956, P, 0.000001);

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

    //    DOUBLES_EQUAL (-326.12, fx, 0.000001);
    //    DOUBLES_EQUAL ( 0.,     fy, 0.000001);
    //    DOUBLES_EQUAL ( 0.,     fz, 0.000001);

    //    fx = atom2.kit().potential_force().x();
    //    fy = atom2.kit().potential_force().y();
    //    fz = atom2.kit().potential_force().z();

    //    DOUBLES_EQUAL ( 326.12, fx, 0.000001);
    //    DOUBLES_EQUAL ( 0.,     fy, 0.000001);
    //    DOUBLES_EQUAL ( 0.,     fz, 0.000001);

    //    // finite difference Gradient ------------------------------
    //    atom1.set_x (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();
    //    double F_x = (P - interactions.potential ()) / d_mach_eps_2;
    //    atom1.set_x (0);
    //    DOUBLES_EQUAL (-326.12, F_x, 0.0001);

    //    atom2.set_x (1 + d_mach_eps_2);
    //    //nb.cache_atoms ();
    //    interactions.cache_atoms ();
    //    F_x = (P - interactions.potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL ( 326.12, F_x, 0.0001);
    //}

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

CT  CT  HC  50.0    109.50  changed based on NMA nmodes

     h1
     |
c1 - c2

*/
    //void CT_CT_HC_angle_test () 
    //{
    //    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);

    //    Arbitrary_atom_group atom_group;
    //    atom_group.add (c1);
    //    atom_group.add (c2);
    //    atom_group.add (h1);

    //    Default_bond bond1 (c2, c1);
    //    Default_bond bond2 (c2, h1);

    //    Arbitrary_bond_group bond_group;
    //    bond_group.add (bond1);
    //    bond_group.add (bond2);

    //    Conventional_interactions  interactions 
    //        (atom_group, bond_group, *force_field_);

    //    // Potential -----------------------------
    //    // 50*(Pi*109.5/180 - Pi/2)^2 = 5.79153869369479
    //    interactions.potential ();
    //    double P = interactions.last_angle_potential ();
    //    DOUBLES_EQUAL (5.79153869369479, P, 0.000001);


    //    // finite difference Gradient ------------------------------
    //    // ñ1 (0,0,0)
    //    c1.set_x (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    double F_x = (P - interactions.last_angle_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_x, 0.0001);        c1.set_x (0);

    //    c1.set_y (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    double F_y = (P - interactions.last_angle_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (-34.0339, F_y, 0.0001);  c1.set_y (0);

    //    c1.set_z (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    double F_z = (P - interactions.last_angle_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_z, 0.0001);        c1.set_z (0);

    //    // ñ2 (1,0,0)
    //    c2.set_x (1 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_x = (P - interactions.last_angle_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (-34.0339, F_x, 0.0001);  c2.set_x (1);

    //    c2.set_y (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_y = (P - interactions.last_angle_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL ( 34.0339, F_y, 0.0001);  c2.set_y (0);

    //    c2.set_z (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_z = (P - interactions.last_angle_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_z, 0.0001);        c2.set_z (0);

    //    // h1 (1,1,0)
    //    h1.set_x (1 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_x = (P - interactions.last_angle_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL ( 34.0339, F_x, 0.0001);  h1.set_x (1);

    //    h1.set_y (1 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_y = (P - interactions.last_angle_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_y, 0.0001);        h1.set_y (1);

    //    h1.set_z (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_z = (P - interactions.last_angle_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_z, 0.0001);        h1.set_z (0);

    //    //log() << "###   G_y " << G_y << "\n"; //###   G_y 34.0339

    //    // Gradient ------------------------------
    //    double fx, fy, fz;
    //    interactions.cache_atoms ();
    //    interactions.cache_bonds ();
    //    interactions.zero_forces ();
    //    interactions.add_angle_force (*c2.kit().cache___);
    //    //interactions.add_force   ();  // Does not consider bonds gradient
    //    interactions.flush_force ();

    //    fx = c1.kit().potential_force().x();
    //    fy = c1.kit().potential_force().y();
    //    fz = c1.kit().potential_force().z();
    //    DOUBLES_EQUAL ( 0.,         fx, 0.0001);
    //    DOUBLES_EQUAL (-34.0339,    fy, 0.0001);
    //    DOUBLES_EQUAL ( 0.,         fz, 0.0001);

    //    fx = c2.kit().potential_force().x();
    //    fy = c2.kit().potential_force().y();
    //    fz = c2.kit().potential_force().z();
    //    DOUBLES_EQUAL (-34.0339,    fx, 0.0001);
    //    DOUBLES_EQUAL ( 34.0339,    fy, 0.0001);
    //    DOUBLES_EQUAL ( 0.,         fz, 0.0001);

    //    fx = h1.kit().potential_force().x();
    //    fy = h1.kit().potential_force().y();
    //    fz = h1.kit().potential_force().z();
    //    DOUBLES_EQUAL ( 34.0339,    fx, 0.0001);
    //    DOUBLES_EQUAL ( 0.,         fy, 0.0001);
    //    DOUBLES_EQUAL ( 0.,         fz, 0.0001);
    //}

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

**  CT  CT  **      9   1.40   0.0   3      JCC,7,(1986),230

h1     h2
|     /
c1 - c2

*/
    //void HC_CT_CT_HC_torsion_test ()
    //{
    //    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);

    //    Arbitrary_atom_group atom_group;
    //    atom_group.add (c1);
    //    atom_group.add (c2);
    //    atom_group.add (h1);
    //    atom_group.add (h2);

    //    Default_bond bond0 (h1, c1);
    //    Default_bond bond1 (c1, c2);
    //    Default_bond bond2 (c2, h2);

    //    Arbitrary_bond_group bond_group;
    //    bond_group.add (bond0);
    //    bond_group.add (bond1);
    //    bond_group.add (bond2);

    //    Conventional_interactions  interactions 
    //        (atom_group, bond_group, *force_field_);
    //    
    //    // Potential -----------------------------
    //    // 1.40 * (1 + cos(3 * Pi/4 - 0)) / 9 = 0.0455611665950794
    //    interactions.potential ();
    //    double P = interactions.last_torsion_potential ();
    //    DOUBLES_EQUAL (0.0455611665950794, P, 0.000001);


    //    // finite difference Gradient ------------------------------
    //    // h1 (0,1,0)
    //    h1.set_x (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    double F_x = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_x, 0.0001);        h1.set_x (0);

    //    h1.set_y (1 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    double F_y = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_y, 0.0001);        h1.set_y (1);

    //    h1.set_z (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    double F_z = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (-0.329983, F_z, 0.0001); h1.set_z (0);

    //    // ñ1 (0,0,0)
    //    c1.set_x (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_x = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_x, 0.0001);        c1.set_x (0);

    //    c1.set_y (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_y = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_y, 0.0001);        c1.set_y (0);

    //    c1.set_z (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_z = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL ( 0.329983, F_z, 0.0001); c1.set_z (0);

    //    // ñ2 (1,0,0)
    //    c2.set_x (1 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_x = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0.,        F_x, 0.0001); c2.set_x (1);

    //    c2.set_y (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_y = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL ( 0.164992, F_y, 0.0001); c2.set_y (0);

    //    c2.set_z (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_z = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (-0.164992, F_z, 0.0001); c2.set_z (0);

    //    // h2 (1,1,1)
    //    h2.set_x (1 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_x = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (0., F_x, 0.0001);        h2.set_x (1);

    //    h2.set_y (1 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_y = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL (-0.164992, F_y, 0.0001); h2.set_y (1);

    //    h2.set_z (1 + d_mach_eps_2);
    //    interactions.cache_atoms ();    interactions.potential   ();
    //    F_z = (P - interactions.last_torsion_potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL ( 0.164992, F_z, 0.0001); h2.set_z (1);

    //    //log() << "###   G_x " << G_x << "\n"; //###   G_x  0        
    //    //log() << "###   G_y " << G_y << "\n"; //###   G_y  0.164992 
    //    //log() << "###   G_z " << G_z << "\n"; //###   G_z -0.164992

    //    // Gradient ------------------------------
    //    double fx, fy, fz;
    //    interactions.cache_atoms ();
    //    interactions.cache_bonds ();
    //    interactions.zero_forces ();
    //    interactions.add_torsion_force (*bond1.kit().cache___);
    //    // Does not consider bonds and angles gradients
    //    //interactions.add_force   ();  
    //    interactions.flush_force ();

    //    fx = h1.kit().potential_force().x();
    //    fy = h1.kit().potential_force().y();
    //    fz = h1.kit().potential_force().z();
    //    DOUBLES_EQUAL ( 0.,         fx, 0.0001);
    //    DOUBLES_EQUAL ( 0.,         fy, 0.0001);
    //    DOUBLES_EQUAL (-0.329983,   fz, 0.0001);

    //    fx = c1.kit().potential_force().x();
    //    fy = c1.kit().potential_force().y();
    //    fz = c1.kit().potential_force().z();
    //    DOUBLES_EQUAL ( 0.,         fx, 0.0001);
    //    DOUBLES_EQUAL ( 0.,         fy, 0.0001);
    //    DOUBLES_EQUAL ( 0.329983,   fz, 0.0001);

    //    fx = c2.kit().potential_force().x();
    //    fy = c2.kit().potential_force().y();
    //    fz = c2.kit().potential_force().z();
    //    DOUBLES_EQUAL ( 0.,         fx, 0.0001);
    //    DOUBLES_EQUAL ( 0.164992,   fy, 0.0001);
    //    DOUBLES_EQUAL (-0.164992,   fz, 0.0001);

    //    fx = h2.kit().potential_force().x();
    //    fy = h2.kit().potential_force().y();
    //    fz = h2.kit().potential_force().z();
    //    DOUBLES_EQUAL ( 0.,         fx, 0.0001);
    //    DOUBLES_EQUAL (-0.164992,   fy, 0.0001);
    //    DOUBLES_EQUAL ( 0.164992,   fz, 0.0001);
    //}

    //void Van_der_Waals_test ()
    //{
    //    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);

    //    Arbitrary_atom_group atom_group;
    //    atom_group.add (c1);
    //    atom_group.add (c2);

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

    //    //  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
    //    interactions.potential ();
    //    double P = interactions.last_Van_der_Waals_potential ();
    //    DOUBLES_EQUAL (1042404.61845581201, P, 0.000001);

    //     // Commonplace_interactions
    //    {
    //        Commonplace_interactions interactions (atom_group, *force_field_);

    //        interactions.potential ();
    //        double P = interactions.last_Van_der_Waals_potential ();
    //        DOUBLES_EQUAL (1042404.61845581201, P, 0.000001);
    //    }

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

    //    DOUBLES_EQUAL (-1.25129e+007,   fx, 100);
    //    DOUBLES_EQUAL ( 0.,             fy, 100);
    //    DOUBLES_EQUAL ( 0.,             fz, 100);

    //    fx = c2.kit().potential_force().x();
    //    fy = c2.kit().potential_force().y();
    //    fz = c2.kit().potential_force().z();

    //    DOUBLES_EQUAL ( 1.25129e+007,   fx, 100);
    //    DOUBLES_EQUAL ( 0.,             fy, 100);
    //    DOUBLES_EQUAL ( 0.,             fz, 100);

    //    // finite difference Gradient ------------------------------
    //    c1.set_x (0 + d_mach_eps_2);
    //    interactions.cache_atoms ();
    //    double F_x = (P - interactions.potential ()) / d_mach_eps_2;
    //    c1.set_x (0);
    //    DOUBLES_EQUAL (-1.25129e+007,   F_x, 100);

    //    c2.set_x (1 + d_mach_eps_2);
    //    interactions.cache_atoms ();
    //    F_x = (P - interactions.potential ()) / d_mach_eps_2;
    //    DOUBLES_EQUAL ( 1.25129e+007,   F_x, 100);

    //}

    //void ethane_test ()
    //{
    //    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 atom_group;
    //    atom_group.add (c1);    atom_group.add (c2);    
    //    atom_group.add (h1);    atom_group.add (h2);
    //    atom_group.add (h3);    atom_group.add (h4);    
    //    atom_group.add (h5);    atom_group.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 bond_group;
    //    bond_group.add (cc );   
    //    bond_group.add (ch1);   bond_group.add (ch2);   bond_group.add (ch3);
    //    bond_group.add (ch4);   bond_group.add (ch5);   bond_group.add (ch6);

    //    Conventional_interactions  interactions 
    //        (atom_group, bond_group, *force_field_);

    //    double P = interactions.potential ();
    //    DOUBLES_EQUAL (2.427278, P, 0.00001);
    //    DOUBLES_EQUAL (0.06076,     interactions.last_bond_potential(),          0.00001 );
    //    DOUBLES_EQUAL (0.000128842, interactions.last_angle_potential(),         0.0000001);
    //    DOUBLES_EQUAL (2.10002,     interactions.last_torsion_potential(),       0.00001 );
    //    DOUBLES_EQUAL (0.266372,    interactions.last_Van_der_Waals_potential(), 0.00001 );
    //    DOUBLES_EQUAL (0.,          interactions.last_electrostatic_potential(), 0.00001 );
    //}
};

Covalent_interactions_test test_covalent_interactions ("correctness");

}//MM
