Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

lmp: support unit real #2775

Merged
merged 16 commits into from
Sep 2, 2023
80 changes: 45 additions & 35 deletions source/lmp/pair_deepmd.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -344,11 +344,17 @@ 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\"");
}
conversion_factor = utils::get_conversion_factor(utils::ENERGY, unit_convert);
wanghan-iapcm marked this conversation as resolved.
Show resolved Hide resolved
restartinfo = 1;
#if LAMMPS_VERSION_NUMBER >= 20201130
centroidstressflag =
Expand All @@ -361,6 +367,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;
Expand Down Expand Up @@ -576,7 +584,7 @@ void PairDeepMD::compute(int eflag, int vflag) {
}
if (eflag_atom) {
for (int ii = 0; ii < nlocal; ++ii) {
eatom[ii] += deatom[ii];
eatom[ii] += conversion_factor * deatom[ii];
}
}
// Added by Davide Tisi 2020
Expand All @@ -590,15 +598,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] += conversion_factor * dvatom[9 * ii + 0]; // xx
cvatom[ii][1] += conversion_factor * dvatom[9 * ii + 4]; // yy
cvatom[ii][2] += conversion_factor * dvatom[9 * ii + 8]; // zz
cvatom[ii][3] += conversion_factor * dvatom[9 * ii + 3]; // xy
cvatom[ii][4] += conversion_factor * dvatom[9 * ii + 6]; // xz
cvatom[ii][5] += conversion_factor * dvatom[9 * ii + 7]; // yz
cvatom[ii][6] += conversion_factor * dvatom[9 * ii + 1]; // yx
cvatom[ii][7] += conversion_factor * dvatom[9 * ii + 2]; // zx
cvatom[ii][8] += conversion_factor * dvatom[9 * ii + 5]; // zy
}
}
}
Expand Down Expand Up @@ -628,7 +636,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] += conversion_factor * deatom[ii];
}
}
// Added by Davide Tisi 2020
Expand All @@ -642,15 +650,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] += conversion_factor * dvatom[9 * ii + 0]; // xx
cvatom[ii][1] += conversion_factor * dvatom[9 * ii + 4]; // yy
cvatom[ii][2] += conversion_factor * dvatom[9 * ii + 8]; // zz
cvatom[ii][3] += conversion_factor * dvatom[9 * ii + 3]; // xy
cvatom[ii][4] += conversion_factor * dvatom[9 * ii + 6]; // xz
cvatom[ii][5] += conversion_factor * dvatom[9 * ii + 7]; // yz
cvatom[ii][6] += conversion_factor * dvatom[9 * ii + 1]; // yx
cvatom[ii][7] += conversion_factor * dvatom[9 * ii + 2]; // zx
cvatom[ii][8] += conversion_factor * dvatom[9 * ii + 5]; // zy
}
}
if (out_freq > 0 && update->ntimestep % out_freq == 0) {
Expand Down Expand Up @@ -774,7 +782,7 @@ void PairDeepMD::compute(int eflag, int vflag) {
if (!atom->sp_flag) {
for (int ii = 0; ii < nall; ++ii) {
for (int dd = 0; dd < 3; ++dd) {
f[ii][dd] += scale[1][1] * dforce[3 * ii + dd];
f[ii][dd] += conversion_factor * dforce[3 * ii + dd];
wanghan-iapcm marked this conversation as resolved.
Show resolved Hide resolved
}
}
} else {
Expand All @@ -783,12 +791,14 @@ void PairDeepMD::compute(int eflag, int vflag) {
for (int ii = 0; ii < nall; ++ii) {
for (int dd = 0; dd < 3; ++dd) {
int new_idx = new_idx_map[ii];
f[ii][dd] += scale[1][1] * dforce[3 * new_idx + dd];
f[ii][dd] += conversion_factor * dforce[3 * new_idx + dd];
if (dtype[ii] < numb_types_spin && ii < nlocal) {
fm[ii][dd] += scale[1][1] * dforce[3 * (new_idx + nlocal) + dd] /
fm[ii][dd] += conversion_factor *
dforce[3 * (new_idx + nlocal) + dd] /
(hbar / spin_norm[dtype[ii]]);
} else if (dtype[ii] < numb_types_spin) {
fm[ii][dd] += scale[1][1] * dforce[3 * (new_idx + nghost) + dd] /
fm[ii][dd] += conversion_factor *
dforce[3 * (new_idx + nghost) + dd] /
(hbar / spin_norm[dtype[ii]]);
}
}
Expand All @@ -803,15 +813,15 @@ void PairDeepMD::compute(int eflag, int vflag) {

// accumulate energy and virial
if (eflag) {
eng_vdwl += scale[1][1] * dener;
eng_vdwl += conversion_factor * dener;
}
if (vflag) {
virial[0] += 1.0 * dvirial[0] * scale[1][1];
virial[1] += 1.0 * dvirial[4] * scale[1][1];
virial[2] += 1.0 * dvirial[8] * scale[1][1];
virial[3] += 1.0 * dvirial[3] * scale[1][1];
virial[4] += 1.0 * dvirial[6] * scale[1][1];
virial[5] += 1.0 * dvirial[7] * scale[1][1];
virial[0] += 1.0 * dvirial[0] * conversion_factor;
virial[1] += 1.0 * dvirial[4] * conversion_factor;
virial[2] += 1.0 * dvirial[8] * conversion_factor;
virial[3] += 1.0 * dvirial[3] * conversion_factor;
virial[4] += 1.0 * dvirial[6] * conversion_factor;
virial[5] += 1.0 * dvirial[7] * conversion_factor;
}
}

Expand Down
1 change: 1 addition & 0 deletions source/lmp/pair_deepmd.h
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,7 @@ class PairDeepMD : public Pair {
tagint *tagsend, *tagrecv;
double *stdfsend, *stdfrecv;
std::vector<int> type_idx_map;
double conversion_factor;
};

} // namespace LAMMPS_NS
Expand Down
55 changes: 52 additions & 3 deletions source/lmp/tests/test_lammps.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -244,17 +246,22 @@ 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")
lammps.neigh_modify("every 10 delay 0 check no")
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

Expand All @@ -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("* *")
Expand Down Expand Up @@ -452,3 +466,38 @@ 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)