diff --git a/source/lmp/pair_deepmd.cpp b/source/lmp/pair_deepmd.cpp index f0e0f23096..f285bed8ea 100644 --- a/source/lmp/pair_deepmd.cpp +++ b/source/lmp/pair_deepmd.cpp @@ -344,11 +344,18 @@ PairDeepMD::PairDeepMD(LAMMPS *lmp) if (lmp->citeme) { lmp->citeme->add(cite_user_deepmd_package); } - if (strcmp(update->unit_style, "metal") != 0) { - error->all( - FLERR, - "Pair deepmd requires metal unit, please set it by \"units metal\""); + int unit_convert; + if (strcmp(update->unit_style, "metal") == 0) { + unit_convert = utils::NOCONVERT; + } else if (strcmp(update->unit_style, "real") == 0) { + unit_convert = utils::METAL2REAL; + } else { + error->all(FLERR, + "Pair deepmd requires metal or real unit, please set it by " + "\"units metal\" or \"units real\""); } + ener_unit_cvt_factor = + utils::get_conversion_factor(utils::ENERGY, unit_convert); restartinfo = 1; #if LAMMPS_VERSION_NUMBER >= 20201130 centroidstressflag = @@ -361,6 +368,8 @@ PairDeepMD::PairDeepMD(LAMMPS *lmp) pppmflag = 1; respa_enable = 0; writedata = 0; + unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); + cutoff = 0.; numb_types = 0; numb_types_spin = 0; @@ -576,7 +585,7 @@ void PairDeepMD::compute(int eflag, int vflag) { } if (eflag_atom) { for (int ii = 0; ii < nlocal; ++ii) { - eatom[ii] += deatom[ii]; + eatom[ii] += scale[1][1] * deatom[ii]; } } // Added by Davide Tisi 2020 @@ -590,15 +599,15 @@ void PairDeepMD::compute(int eflag, int vflag) { // vatom[ii][3] += 1.0 * dvatom[9*ii+3]; // vatom[ii][4] += 1.0 * dvatom[9*ii+6]; // vatom[ii][5] += 1.0 * dvatom[9*ii+7]; - cvatom[ii][0] += 1.0 * dvatom[9 * ii + 0]; // xx - cvatom[ii][1] += 1.0 * dvatom[9 * ii + 4]; // yy - cvatom[ii][2] += 1.0 * dvatom[9 * ii + 8]; // zz - cvatom[ii][3] += 1.0 * dvatom[9 * ii + 3]; // xy - cvatom[ii][4] += 1.0 * dvatom[9 * ii + 6]; // xz - cvatom[ii][5] += 1.0 * dvatom[9 * ii + 7]; // yz - cvatom[ii][6] += 1.0 * dvatom[9 * ii + 1]; // yx - cvatom[ii][7] += 1.0 * dvatom[9 * ii + 2]; // zx - cvatom[ii][8] += 1.0 * dvatom[9 * ii + 5]; // zy + cvatom[ii][0] += scale[1][1] * dvatom[9 * ii + 0]; // xx + cvatom[ii][1] += scale[1][1] * dvatom[9 * ii + 4]; // yy + cvatom[ii][2] += scale[1][1] * dvatom[9 * ii + 8]; // zz + cvatom[ii][3] += scale[1][1] * dvatom[9 * ii + 3]; // xy + cvatom[ii][4] += scale[1][1] * dvatom[9 * ii + 6]; // xz + cvatom[ii][5] += scale[1][1] * dvatom[9 * ii + 7]; // yz + cvatom[ii][6] += scale[1][1] * dvatom[9 * ii + 1]; // yx + cvatom[ii][7] += scale[1][1] * dvatom[9 * ii + 2]; // zx + cvatom[ii][8] += scale[1][1] * dvatom[9 * ii + 5]; // zy } } } @@ -628,7 +637,7 @@ void PairDeepMD::compute(int eflag, int vflag) { dvatom = all_atom_virial[0]; if (eflag_atom) { for (int ii = 0; ii < nlocal; ++ii) { - eatom[ii] += deatom[ii]; + eatom[ii] += scale[1][1] * deatom[ii]; } } // Added by Davide Tisi 2020 @@ -642,15 +651,15 @@ void PairDeepMD::compute(int eflag, int vflag) { // vatom[ii][3] += 1.0 * dvatom[9*ii+3]; // vatom[ii][4] += 1.0 * dvatom[9*ii+6]; // vatom[ii][5] += 1.0 * dvatom[9*ii+7]; - cvatom[ii][0] += 1.0 * dvatom[9 * ii + 0]; // xx - cvatom[ii][1] += 1.0 * dvatom[9 * ii + 4]; // yy - cvatom[ii][2] += 1.0 * dvatom[9 * ii + 8]; // zz - cvatom[ii][3] += 1.0 * dvatom[9 * ii + 3]; // xy - cvatom[ii][4] += 1.0 * dvatom[9 * ii + 6]; // xz - cvatom[ii][5] += 1.0 * dvatom[9 * ii + 7]; // yz - cvatom[ii][6] += 1.0 * dvatom[9 * ii + 1]; // yx - cvatom[ii][7] += 1.0 * dvatom[9 * ii + 2]; // zx - cvatom[ii][8] += 1.0 * dvatom[9 * ii + 5]; // zy + cvatom[ii][0] += scale[1][1] * dvatom[9 * ii + 0]; // xx + cvatom[ii][1] += scale[1][1] * dvatom[9 * ii + 4]; // yy + cvatom[ii][2] += scale[1][1] * dvatom[9 * ii + 8]; // zz + cvatom[ii][3] += scale[1][1] * dvatom[9 * ii + 3]; // xy + cvatom[ii][4] += scale[1][1] * dvatom[9 * ii + 6]; // xz + cvatom[ii][5] += scale[1][1] * dvatom[9 * ii + 7]; // yz + cvatom[ii][6] += scale[1][1] * dvatom[9 * ii + 1]; // yx + cvatom[ii][7] += scale[1][1] * dvatom[9 * ii + 2]; // zx + cvatom[ii][8] += scale[1][1] * dvatom[9 * ii + 5]; // zy } } if (out_freq > 0 && update->ntimestep % out_freq == 0) { @@ -719,6 +728,12 @@ void PairDeepMD::compute(int eflag, int vflag) { all_v_avg = sqrt(all_v_avg / 9); } if (rank == 0) { + all_v_max *= scale[1][1]; + all_v_min *= scale[1][1]; + all_v_avg *= scale[1][1]; + all_f_max *= scale[1][1]; + all_f_min *= scale[1][1]; + all_f_avg *= scale[1][1]; fp << setw(12) << update->ntimestep << " " << setw(18) << all_v_max << " " << setw(18) << all_v_min << " " << setw(18) << all_v_avg << " " << setw(18) << all_f_max << " " << setw(18) << all_f_min @@ -744,7 +759,7 @@ void PairDeepMD::compute(int eflag, int vflag) { displacements, MPI_DOUBLE, 0, world); if (rank == 0) { for (int dd = 0; dd < all_nlocal; ++dd) { - std_f_all[tagrecv[dd] - 1] = stdfrecv[dd]; + std_f_all[tagrecv[dd] - 1] = stdfrecv[dd] * scale[1][1]; } for (int dd = 0; dd < all_nlocal; ++dd) { fp << " " << setw(18) << std_f_all[dd]; @@ -838,7 +853,7 @@ void PairDeepMD::allocate() { continue; } setflag[i][j] = 1; - scale[i][j] = 1; + scale[i][j] = 1.0 * ener_unit_cvt_factor; } } } @@ -998,11 +1013,11 @@ void PairDeepMD::settings(int narg, char **arg) { iarg += 1; } else if (string(arg[iarg]) == string("relative")) { out_rel = 1; - eps = atof(arg[iarg + 1]); + eps = atof(arg[iarg + 1]) / ener_unit_cvt_factor; iarg += 2; } else if (string(arg[iarg]) == string("relative_v")) { out_rel_v = 1; - eps_v = atof(arg[iarg + 1]); + eps_v = atof(arg[iarg + 1]) / ener_unit_cvt_factor; iarg += 2; } else if (string(arg[iarg]) == string("virtual_len")) { virtual_len.resize(numb_types_spin); @@ -1174,7 +1189,7 @@ void PairDeepMD::coeff(int narg, char **arg) { for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo, i); j <= jhi; j++) { setflag[i][j] = 1; - scale[i][j] = 1.0; + scale[i][j] = 1.0 * ener_unit_cvt_factor; if (i > numb_types || j > numb_types) { char warning_msg[1024]; sprintf(warning_msg, @@ -1221,7 +1236,7 @@ double PairDeepMD::init_one(int i, int j) { } if (setflag[i][j] == 0) { - scale[i][j] = 1.0; + scale[i][j] = 1.0 * ener_unit_cvt_factor; } scale[j][i] = scale[i][j]; diff --git a/source/lmp/pair_deepmd.h b/source/lmp/pair_deepmd.h index feff28b9a4..e811bc99b9 100644 --- a/source/lmp/pair_deepmd.h +++ b/source/lmp/pair_deepmd.h @@ -131,6 +131,7 @@ class PairDeepMD : public Pair { tagint *tagsend, *tagrecv; double *stdfsend, *stdfrecv; std::vector type_idx_map; + double ener_unit_cvt_factor; }; } // namespace LAMMPS_NS diff --git a/source/lmp/tests/test_lammps.py b/source/lmp/tests/test_lammps.py index 2920615b8e..78eaf7ea4e 100644 --- a/source/lmp/tests/test_lammps.py +++ b/source/lmp/tests/test_lammps.py @@ -217,6 +217,8 @@ # https://github.com/lammps/lammps/blob/1e1311cf401c5fc2614b5d6d0ff3230642b76597/src/update.cpp#L193 nktv2p = 1.6021765e6 +nktv2p_real = 68568.415 +metal2real = 23.060549 sp.check_output( "{} -m deepmd convert-from pbtxt -i {} -o {}".format( @@ -244,9 +246,9 @@ def teardown_module(): os.remove(data_type_map_file) -def _lammps(data_file) -> PyLammps: +def _lammps(data_file, units="metal") -> PyLammps: lammps = PyLammps() - lammps.units("metal") + lammps.units(units) lammps.boundary("p p p") lammps.atom_style("atomic") lammps.neighbor("2.0 bin") @@ -254,7 +256,12 @@ def _lammps(data_file) -> PyLammps: lammps.read_data(data_file.resolve()) lammps.mass("1 16") lammps.mass("2 2") - lammps.timestep(0.0005) + if units == "metal": + lammps.timestep(0.0005) + elif units == "real": + lammps.timestep(0.5) + else: + raise ValueError("units should be metal or real") lammps.fix("1 all nve") return lammps @@ -273,6 +280,13 @@ def lammps_type_map(): lmp.close() +@pytest.fixture +def lammps_real(): + lmp = _lammps(data_file=data_file, units="real") + yield lmp + lmp.close() + + def test_pair_deepmd(lammps): lammps.pair_style(f"deepmd {pb_file.resolve()}") lammps.pair_coeff("* *") @@ -452,3 +466,186 @@ def test_pair_deepmd_type_map(lammps_type_map): expected_f[lammps_type_map.atoms[ii].id - 1] ) lammps_type_map.run(1) + + +def test_pair_deepmd_real(lammps_real): + lammps_real.pair_style(f"deepmd {pb_file.resolve()}") + lammps_real.pair_coeff("* *") + lammps_real.run(0) + assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real) + for ii in range(6): + assert lammps_real.atoms[ii].force == pytest.approx( + expected_f[lammps_real.atoms[ii].id - 1] * metal2real + ) + lammps_real.run(1) + + +def test_pair_deepmd_virial_real(lammps_real): + lammps_real.pair_style(f"deepmd {pb_file.resolve()}") + lammps_real.pair_coeff("* *") + lammps_real.compute("virial all centroid/stress/atom NULL pair") + for ii in range(9): + jj = [0, 4, 8, 3, 6, 7, 1, 2, 5][ii] + lammps_real.variable(f"virial{jj} atom c_virial[{ii+1}]") + lammps_real.dump( + "1 all custom 1 dump id " + " ".join([f"v_virial{ii}" for ii in range(9)]) + ) + lammps_real.run(0) + assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real) + for ii in range(6): + assert lammps_real.atoms[ii].force == pytest.approx( + expected_f[lammps_real.atoms[ii].id - 1] * metal2real + ) + idx_map = lammps_real.lmp.numpy.extract_atom("id") - 1 + for ii in range(9): + assert np.array( + lammps_real.variables[f"virial{ii}"].value + ) / nktv2p_real == pytest.approx(expected_v[idx_map, ii] * metal2real) + + +def test_pair_deepmd_model_devi_real(lammps_real): + lammps_real.pair_style( + "deepmd {} {} out_file {} out_freq 1 atomic".format( + pb_file.resolve(), pb_file2.resolve(), md_file.resolve() + ) + ) + lammps_real.pair_coeff("* *") + lammps_real.run(0) + assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real) + for ii in range(6): + assert lammps_real.atoms[ii].force == pytest.approx( + expected_f[lammps_real.atoms[ii].id - 1] * metal2real + ) + # load model devi + md = np.loadtxt(md_file.resolve()) + expected_md_f = np.linalg.norm(np.std([expected_f, expected_f2], axis=0), axis=1) + assert md[7:] == pytest.approx(expected_md_f * metal2real) + assert md[4] == pytest.approx(np.max(expected_md_f) * metal2real) + assert md[5] == pytest.approx(np.min(expected_md_f) * metal2real) + assert md[6] == pytest.approx(np.mean(expected_md_f) * metal2real) + expected_md_v = ( + np.std([np.sum(expected_v, axis=0), np.sum(expected_v2, axis=0)], axis=0) / 6 + ) + assert md[1] == pytest.approx(np.max(expected_md_v) * metal2real) + assert md[2] == pytest.approx(np.min(expected_md_v) * metal2real) + assert md[3] == pytest.approx( + np.sqrt(np.mean(np.square(expected_md_v))) * metal2real + ) + + +def test_pair_deepmd_model_devi_virial_real(lammps_real): + lammps_real.pair_style( + "deepmd {} {} out_file {} out_freq 1 atomic".format( + pb_file.resolve(), pb_file2.resolve(), md_file.resolve() + ) + ) + lammps_real.pair_coeff("* *") + lammps_real.compute("virial all centroid/stress/atom NULL pair") + for ii in range(9): + jj = [0, 4, 8, 3, 6, 7, 1, 2, 5][ii] + lammps_real.variable(f"virial{jj} atom c_virial[{ii+1}]") + lammps_real.dump( + "1 all custom 1 dump id " + " ".join([f"v_virial{ii}" for ii in range(9)]) + ) + lammps_real.run(0) + assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real) + for ii in range(6): + assert lammps_real.atoms[ii].force == pytest.approx( + expected_f[lammps_real.atoms[ii].id - 1] * metal2real + ) + idx_map = lammps_real.lmp.numpy.extract_atom("id") - 1 + for ii in range(9): + assert np.array( + lammps_real.variables[f"virial{ii}"].value + ) / nktv2p_real == pytest.approx(expected_v[idx_map, ii] * metal2real) + # load model devi + md = np.loadtxt(md_file.resolve()) + expected_md_f = np.linalg.norm(np.std([expected_f, expected_f2], axis=0), axis=1) + assert md[7:] == pytest.approx(expected_md_f * metal2real) + assert md[4] == pytest.approx(np.max(expected_md_f) * metal2real) + assert md[5] == pytest.approx(np.min(expected_md_f) * metal2real) + assert md[6] == pytest.approx(np.mean(expected_md_f) * metal2real) + expected_md_v = ( + np.std([np.sum(expected_v, axis=0), np.sum(expected_v2, axis=0)], axis=0) / 6 + ) + assert md[1] == pytest.approx(np.max(expected_md_v) * metal2real) + assert md[2] == pytest.approx(np.min(expected_md_v) * metal2real) + assert md[3] == pytest.approx( + np.sqrt(np.mean(np.square(expected_md_v))) * metal2real + ) + + +def test_pair_deepmd_model_devi_atomic_relative_real(lammps_real): + relative = 1.0 + lammps_real.pair_style( + "deepmd {} {} out_file {} out_freq 1 atomic relative {}".format( + pb_file.resolve(), + pb_file2.resolve(), + md_file.resolve(), + relative * metal2real, + ) + ) + lammps_real.pair_coeff("* *") + lammps_real.run(0) + assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real) + for ii in range(6): + assert lammps_real.atoms[ii].force == pytest.approx( + expected_f[lammps_real.atoms[ii].id - 1] * metal2real + ) + # load model devi + md = np.loadtxt(md_file.resolve()) + norm = np.linalg.norm(np.mean([expected_f, expected_f2], axis=0), axis=1) + expected_md_f = np.linalg.norm(np.std([expected_f, expected_f2], axis=0), axis=1) + expected_md_f /= norm + relative + assert md[7:] == pytest.approx(expected_md_f * metal2real) + assert md[4] == pytest.approx(np.max(expected_md_f) * metal2real) + assert md[5] == pytest.approx(np.min(expected_md_f) * metal2real) + assert md[6] == pytest.approx(np.mean(expected_md_f) * metal2real) + expected_md_v = ( + np.std([np.sum(expected_v, axis=0), np.sum(expected_v2, axis=0)], axis=0) / 6 + ) + assert md[1] == pytest.approx(np.max(expected_md_v) * metal2real) + assert md[2] == pytest.approx(np.min(expected_md_v) * metal2real) + assert md[3] == pytest.approx( + np.sqrt(np.mean(np.square(expected_md_v))) * metal2real + ) + + +def test_pair_deepmd_model_devi_atomic_relative_v_real(lammps_real): + relative = 1.0 + lammps_real.pair_style( + "deepmd {} {} out_file {} out_freq 1 atomic relative_v {}".format( + pb_file.resolve(), + pb_file2.resolve(), + md_file.resolve(), + relative * metal2real, + ) + ) + lammps_real.pair_coeff("* *") + lammps_real.run(0) + assert lammps_real.eval("pe") == pytest.approx(expected_e * metal2real) + for ii in range(6): + assert lammps_real.atoms[ii].force == pytest.approx( + expected_f[lammps_real.atoms[ii].id - 1] * metal2real + ) + md = np.loadtxt(md_file.resolve()) + expected_md_f = np.linalg.norm(np.std([expected_f, expected_f2], axis=0), axis=1) + assert md[7:] == pytest.approx(expected_md_f * metal2real) + assert md[4] == pytest.approx(np.max(expected_md_f) * metal2real) + assert md[5] == pytest.approx(np.min(expected_md_f) * metal2real) + assert md[6] == pytest.approx(np.mean(expected_md_f) * metal2real) + expected_md_v = ( + np.std([np.sum(expected_v, axis=0), np.sum(expected_v2, axis=0)], axis=0) / 6 + ) + norm = ( + np.abs( + np.mean([np.sum(expected_v, axis=0), np.sum(expected_v2, axis=0)], axis=0) + ) + / 6 + ) + expected_md_v /= norm + relative + assert md[1] == pytest.approx(np.max(expected_md_v) * metal2real) + assert md[2] == pytest.approx(np.min(expected_md_v) * metal2real) + assert md[3] == pytest.approx( + np.sqrt(np.mean(np.square(expected_md_v))) * metal2real + )