/* ----------------------------------------------------------------------
   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
   http://lammps.sandia.gov, Sandia National Laboratories
   Steve Plimpton, sjplimp@sandia.gov

   Copyright (2003) Sandia Corporation.  Under the terms of Contract
   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
   certain rights in this software.  This software is distributed under
   the GNU General Public License.

   See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */

#include "atom.h"
#include "force.h"
#include "info.h"
#include "input.h"
#include "lammps.h"
#include "pair.h"
#include "utils.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"

#include <cstdio>
#include <cstring>
#include <mpi.h>

// whether to print verbose output (i.e. not capturing LAMMPS screen output).
bool verbose = false;

namespace LAMMPS_NS {
using ::testing::Eq;

// eV to kcal/mol conversion constant (CODATA 2018)
const double ev_convert = utils::get_conversion_factor(utils::ENERGY, utils::METAL2REAL);
const double rel_error  = 1.0e-13;

class PairUnitConvertTest : public ::testing::Test {
protected:
    LAMMPS *lmp;
    Info *info;
    double fold[4][3];

    void SetUp() override
    {
        const char *args[] = {"PairUnitConvertTest", "-log", "none", "-echo", "screen", "-nocite"};
        char **argv        = (char **)args;
        int argc           = sizeof(args) / sizeof(char *);
        if (!verbose) ::testing::internal::CaptureStdout();
        lmp = new LAMMPS(argc, argv, MPI_COMM_WORLD);
        if (!verbose) ::testing::internal::GetCapturedStdout();
        ASSERT_NE(lmp, nullptr);
        if (!verbose) ::testing::internal::CaptureStdout();
        info = new Info(lmp);
        lmp->input->one("units metal");
        lmp->input->one("dimension 3");
        lmp->input->one("region box block -4 4 -4 4 -4 4");
        lmp->input->one("create_box 2 box");
        lmp->input->one("create_atoms 1 single -1.1  1.2  0.0 units box");
        lmp->input->one("create_atoms 1 single -1.2 -1.1  0.0 units box");
        lmp->input->one("create_atoms 2 single  0.9  1.0  0.0 units box");
        lmp->input->one("create_atoms 2 single  1.0 -0.9  0.0 units box");
        lmp->input->one("pair_style zero 4.0");
        lmp->input->one("pair_coeff * *");
        lmp->input->one("mass * 1.0");
        lmp->input->one("write_data test_pair_unit_convert.data nocoeff");
        lmp->input->one("clear");
        if (!verbose) ::testing::internal::GetCapturedStdout();
    }

    void TearDown() override
    {
        if (!verbose) ::testing::internal::CaptureStdout();
        delete info;
        delete lmp;
        if (!verbose) ::testing::internal::GetCapturedStdout();
        remove("test_pair_unit_convert.data");
    }
};

TEST_F(PairUnitConvertTest, zero)
{
    // check if the prerequisite pair style is available
    if (!info->has_style("pair", "zero")) GTEST_SKIP();

    if (!verbose) ::testing::internal::CaptureStdout();
    lmp->input->one("units metal");
    lmp->input->one("read_data test_pair_unit_convert.data");
    lmp->input->one("pair_style zero 6.0");
    lmp->input->one("pair_coeff * *");
    lmp->input->one("run 0 post no");
    if (!verbose) ::testing::internal::GetCapturedStdout();

    // copy energy and force from first step
    double eold = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
    double **f  = lmp->atom->f;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 3; ++j)
            fold[i][j] = f[i][j];

    if (!verbose) ::testing::internal::CaptureStdout();
    lmp->input->one("clear");
    lmp->input->one("units real");
    lmp->input->one("read_data test_pair_unit_convert.data");
    lmp->input->one("pair_style zero 6.0");
    lmp->input->one("pair_coeff * *");
    lmp->input->one("run 0 post no");
    if (!verbose) ::testing::internal::GetCapturedStdout();

    double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
    EXPECT_NEAR(ev_convert * eold, enew, fabs(enew * rel_error));

    f = lmp->atom->f;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 3; ++j)
            EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
}

TEST_F(PairUnitConvertTest, lj_cut)
{
    // check if the prerequisite pair style is available
    if (!info->has_style("pair", "lj/cut")) GTEST_SKIP();

    if (!verbose) ::testing::internal::CaptureStdout();
    lmp->input->one("units metal");
    lmp->input->one("read_data test_pair_unit_convert.data");
    lmp->input->one("pair_style lj/cut 6.0");
    lmp->input->one("pair_coeff * * 0.01014286346782117 2.0");
    lmp->input->one("run 0 post no");
    if (!verbose) ::testing::internal::GetCapturedStdout();

    // copy energy and force from first step
    double eold = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
    double **f  = lmp->atom->f;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 3; ++j)
            fold[i][j] = f[i][j];

    if (!verbose) ::testing::internal::CaptureStdout();
    lmp->input->one("clear");
    lmp->input->one("units real");
    lmp->input->one("read_data test_pair_unit_convert.data");
    lmp->input->one("pair_style lj/cut 6.0");
    lmp->input->one("pair_coeff * * 0.2339 2.0");
    lmp->input->one("run 0 post no");
    if (!verbose) ::testing::internal::GetCapturedStdout();

    double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
    EXPECT_NEAR(ev_convert * eold, enew, fabs(enew * rel_error));

    f = lmp->atom->f;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 3; ++j)
            EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
}

TEST_F(PairUnitConvertTest, tersoff)
{
    // check if the prerequisite pair style is available
    if (!info->has_style("pair", "tersoff")) GTEST_SKIP();

    if (!verbose) ::testing::internal::CaptureStdout();
    lmp->input->one("units metal");
    lmp->input->one("read_data test_pair_unit_convert.data");
    lmp->input->one("pair_style tersoff");
    lmp->input->one("pair_coeff * * SiC.tersoff Si C");
    lmp->input->one("run 0 post no");
    if (!verbose) ::testing::internal::GetCapturedStdout();

    // copy energy and force from first step
    double eold = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
    double **f  = lmp->atom->f;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 3; ++j)
            fold[i][j] = f[i][j];

    if (!verbose) ::testing::internal::CaptureStdout();
    lmp->input->one("clear");
    lmp->input->one("units real");
    lmp->input->one("read_data test_pair_unit_convert.data");
    lmp->input->one("pair_style tersoff");
    lmp->input->one("pair_coeff * * SiC.tersoff Si C");
    lmp->input->one("run 0 post no");
    if (!verbose) ::testing::internal::GetCapturedStdout();

    double enew = lmp->force->pair->eng_vdwl + lmp->force->pair->eng_coul;
    EXPECT_NEAR(ev_convert * eold, enew, fabs(enew * rel_error));

    f = lmp->atom->f;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 3; ++j)
            EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));
}

} // namespace LAMMPS_NS

int main(int argc, char **argv)
{
    MPI_Init(&argc, &argv);
    ::testing::InitGoogleMock(&argc, argv);
    if ((argc > 1) && (strcmp(argv[1], "-v") == 0)) verbose = true;

    int rv = RUN_ALL_TESTS();
    MPI_Finalize();
    return rv;
}