Unverified Commit 85097df2 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

add unit conversion to pair style tersoff/mod and tersoff/mod/c

parent 00332d29

doc/src/pair_tersoff_mod.rst

+1 −1
Original line number Diff line number Diff line
@@ -74,7 +74,7 @@ formulation of the V_ij term, where it contains an additional c0 term. 


.. math::



   V_{ij}  & = f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) + c_0 \right]

   V_{ij}  = f_C(r_{ij}) \left[ f_R(r_{ij}) + b_{ij} f_A(r_{ij}) + c_0 \right]



The modified cutoff function :math:`f_C` proposed by :ref:`(Murty) <Murty>` and

having a continuous second-order differential is employed. The

src/MANYBODY/pair_tersoff_mod.cpp

+11 −1
Original line number Diff line number Diff line
@@ -53,9 +53,14 @@ void PairTersoffMOD::read_file(char *file) 
  // open file on proc 0



  if (comm->me == 0) {

    PotentialFileReader reader(lmp, file, "Tersoff");

    PotentialFileReader reader(lmp, file, "TersoffMod", unit_convert_flag);

    char * line;



    // transparently convert units for supported conversions



    int unit_convert = reader.get_unit_convert();

    double conversion_factor = utils::get_conversion_factor(utils::ENERGY,

                                                            unit_convert);

    while((line = reader.next_line(NPARAMS_PER_LINE))) {

      try {

        ValueTokenizer values(line);
@@ -109,6 +114,11 @@ void PairTersoffMOD::read_file(char *file) 
        params[nparams].c4         = values.next_double();

        params[nparams].c5         = values.next_double();

        params[nparams].powermint = int(params[nparams].powerm);



        if (unit_convert) {

          params[nparams].biga *= conversion_factor;

          params[nparams].bigb *= conversion_factor;

        }

      } catch (TokenizerException & e) {

        error->one(FLERR, e.what());

      }

src/MANYBODY/pair_tersoff_mod_c.cpp

+12 −1
Original line number Diff line number Diff line
@@ -44,9 +44,14 @@ void PairTersoffMODC::read_file(char *file) 
  // open file on proc 0



  if (comm->me == 0) {

    PotentialFileReader reader(lmp, file, "Tersoff");

    PotentialFileReader reader(lmp, file, "TersoffModC", unit_convert_flag);

    char * line;



    // transparently convert units for supported conversions



    int unit_convert = reader.get_unit_convert();

    double conversion_factor = utils::get_conversion_factor(utils::ENERGY,

                                                            unit_convert);

    while((line = reader.next_line(NPARAMS_PER_LINE))) {

      try {

        ValueTokenizer values(line);
@@ -101,6 +106,12 @@ void PairTersoffMODC::read_file(char *file) 
        params[nparams].c5         = values.next_double();

        params[nparams].c0         = values.next_double();

        params[nparams].powermint = int(params[nparams].powerm);



        if (unit_convert) {

          params[nparams].biga *= conversion_factor;

          params[nparams].bigb *= conversion_factor;

          params[nparams].c0   *= conversion_factor;

        }

      } catch (TokenizerException & e) {

        error->one(FLERR, e.what());

      }

unittest/formats/test_pair_unit_convert.cpp

+85 −0
Original line number Diff line number Diff line
@@ -214,6 +214,91 @@ TEST_F(PairUnitConvertTest, tersoff) 
            EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));

}



TEST_F(PairUnitConvertTest, tersoff_mod)

{

    // check if the prerequisite pair style is available

    if (!info->has_style("pair", "tersoff/mod")) 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/mod");

    lmp->input->one("pair_coeff * * Si.tersoff.mod Si Si");

    lmp->input->one("run 0 post no");

    if (!verbose) ::testing::internal::GetCapturedStdout();



    // copy pressure, energy, and force from first step

    double pold;

    lmp->output->thermo->evaluate_keyword("press", &pold);

    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/mod");

    lmp->input->one("pair_coeff * * Si.tersoff.mod Si Si");

    lmp->input->one("run 0 post no");

    if (!verbose) ::testing::internal::GetCapturedStdout();



    double pnew;

    lmp->output->thermo->evaluate_keyword("press", &pnew);

    EXPECT_NEAR(pold, p_convert * pnew, fabs(pnew * rel_error));

    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_mod_c)

{

    // check if the prerequisite pair style is available

    if (!info->has_style("pair", "tersoff/mod/c")) 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/mod/c");

    lmp->input->one("pair_coeff * * Si.tersoff.modc Si Si");

    lmp->input->one("run 0 post no");

    if (!verbose) ::testing::internal::GetCapturedStdout();



    // copy pressure, energy, and force from first step

    double pold;

    lmp->output->thermo->evaluate_keyword("press", &pold);

    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/mod/c");

    lmp->input->one("pair_coeff * * Si.tersoff.modc Si Si");

    lmp->input->one("run 0 post no");

    if (!verbose) ::testing::internal::GetCapturedStdout();



    double pnew;

    lmp->output->thermo->evaluate_keyword("press", &pnew);

    EXPECT_NEAR(pold, p_convert * pnew, fabs(pnew * rel_error));

    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_table)

{

    // check if the prerequisite pair style is available

CRA Git | Maintained and supported by SUSTech CRA and CCSE