Unverified Commit 4824c460 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

add unit conversion checking support to Force::open_potential()

parent 3c9b40a3

src/MANYBODY/pair_eam_cd.cpp

+8 −6
Original line number Diff line number Diff line
@@ -29,6 +29,8 @@ 
#include "memory.h"

#include "error.h"

#include "tokenizer.h"

#include "utils.h"

#include "fmt/format.h"



using namespace LAMMPS_NS;
@@ -40,6 +42,7 @@ PairEAMCD::PairEAMCD(LAMMPS *lmp, int _cdeamVersion) 
{

  single_enable = 0;

  restartinfo = 0;

  unit_convert_flag = utils::get_supported_conversions(utils::ENERGY);



  rhoB = NULL;

  D_values = NULL;
@@ -500,12 +503,11 @@ void PairEAMCD::read_h_coeff(char *filename) 
    FILE *fptr;

    char line[MAXLINE];

    char nextline[MAXLINE];

    fptr = force->open_potential(filename);

    if (fptr == NULL) {

      char str[128];

      snprintf(str,128,"Cannot open EAM potential file %s", filename);

      error->one(FLERR,str);

    }

    int convert_flag = unit_convert_flag;

    fptr = force->open_potential(filename, &convert_flag);

    if (fptr == NULL)

      error->one(FLERR,fmt::format("Cannot open EAMCD potential file {}",

                                   filename));



    // h coefficients are stored at the end of the file.

    // Skip to last line of file.

src/force.cpp

+30 −4
Original line number Diff line number Diff line
@@ -1011,7 +1011,7 @@ tagint Force::tnumeric(const char *file, int line, char *str) 
   if fails, search in dir specified by env variable LAMMPS_POTENTIALS

------------------------------------------------------------------------- */



FILE *Force::open_potential(const char *name)

FILE *Force::open_potential(const char *name, int *auto_convert)

{

  std::string filepath = utils::get_potential_file_path(name);
@@ -1024,9 +1024,35 @@ FILE *Force::open_potential(const char *name) 
      utils::logmesg(lmp, fmt::format("Reading potential file {} "

                                      "with DATE: {}\n", name, date));

    }

    if (!units.empty() && (units != unit_style)) {



    if (auto_convert == nullptr) {

      if (units != unit_style) {

        error->one(FLERR, fmt::format("Potential file {} requires {} units "

                                      "but {} units are in use", name, units,

                                      unit_style));

        return nullptr;

      }

    } else {

      if (units == unit_style) {

        *auto_convert = utils::NOCONVERT;

      } else {

        if ((units == "metal") && (unit_style == "real")

            && (*auto_convert & utils::METAL2REAL)) {

          *auto_convert = utils::METAL2REAL;

        } else if ((units == "real") && (unit_style == "metal")

            && (*auto_convert & utils::REAL2METAL)) {

          *auto_convert = utils::REAL2METAL;

        } else {

          error->one(FLERR, fmt::format("Potential file {} requires {} units "

                                    "but {} units are in use", name, units, unit_style));

                                        "but {} units are in use", name,

                                        units, unit_style));

          return nullptr;

        }

      }

      if (*auto_convert != utils::NOCONVERT)

        lmp->error->warning(FLERR, fmt::format("Converting potential file in "

                                               "{} units to {} units",

                                               units, unit_style));

    }

    return fopen(filepath.c_str(), "r");

  }

src/force.h

+1 −1
Original line number Diff line number Diff line
@@ -135,7 +135,7 @@ class Force : protected Pointers { 
  bigint bnumeric(const char *, int, char *);

  tagint tnumeric(const char *, int, char *);



  FILE *open_potential(const char *);

  FILE *open_potential(const char *, int *auto_convert = nullptr);



  bigint memory_usage();

unittest/formats/test_pair_unit_convert.cpp

+129 −0
Original line number Diff line number Diff line
@@ -224,6 +224,135 @@ TEST_F(PairUnitConvertTest, eam) 
            EXPECT_NEAR(ev_convert * fold[i][j], f[i][j], fabs(f[i][j] * rel_error));

}



TEST_F(PairUnitConvertTest, eam_alloy)

{

    // check if the prerequisite pair style is available

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

    lmp->input->one("pair_coeff * * AlCu.eam.alloy Al Cu");

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

    lmp->input->one("pair_coeff * * AlCu.eam.alloy Al Cu");

    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, eam_fs)

{

    // check if the prerequisite pair style is available

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

    lmp->input->one("pair_coeff * * FeP_mm.eam.fs Fe P");

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

    lmp->input->one("pair_coeff * * FeP_mm.eam.fs Fe P");

    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, eam_cd)

{

    // check if the prerequisite pair style is available

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

    lmp->input->one("pair_coeff * * FeCr.cdeam Cr Fe");

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

    lmp->input->one("pair_coeff * * FeCr.cdeam Cr Fe");

    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, sw)

{

    // check if the prerequisite pair style is available

CRA Git | Maintained and supported by SUSTech CRA and CCSE