diff --git a/doc/src/fix_external.rst b/doc/src/fix_external.rst index 44dd0929ee..3c7bfa315a 100644 --- a/doc/src/fix_external.rst +++ b/doc/src/fix_external.rst @@ -173,9 +173,17 @@ stress/atom ` commands. The former can be accessed by :doc:`thermodynamic output `. The default setting for this fix is :doc:`fix_modify virial yes `. -This fix computes a global scalar which can be accessed by various -:doc:`output commands `. The scalar is the potential -energy discussed above. The scalar stored by this fix is "extensive". +This fix computes a global scalar, a global vector, and a per-atom array +which can be accessed by various :doc:`output commands `. +The scalar is the potential energy discussed above. The scalar stored +by this fix is "extensive". The vector is a custom vector set by the +external program using the :cpp:func:`lammps_fix_external_set_vector() +` and +:cpp:func:`lammps_fix_external_set_vector_length() +` calls of the LAMMPS library +interface or the equivalent call of the Python or Fortran modules. The +per-atom array has 3 column for each atom and is the applied external +force. No parameter of this fix can be used with the *start/stop* keywords of the :doc:`run ` command. diff --git a/src/fix_external.cpp b/src/fix_external.cpp index 349558fe86..1580dd27e1 100644 --- a/src/fix_external.cpp +++ b/src/fix_external.cpp @@ -30,31 +30,33 @@ enum{PF_CALLBACK,PF_ARRAY}; /* ---------------------------------------------------------------------- */ FixExternal::FixExternal(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - fexternal(nullptr), caller_vector(nullptr) + Fix(lmp, narg, arg), fexternal(nullptr), caller_vector(nullptr) { if (narg < 4) error->all(FLERR,"Illegal fix external command"); scalar_flag = 1; global_freq = 1; + peratom_flag = 1; + peratom_freq = 1; + size_peratom_cols = 3; extscalar = 1; energy_global_flag = energy_peratom_flag = 1; virial_global_flag = virial_peratom_flag = 1; thermo_energy = thermo_virial = 1; if (strcmp(arg[3],"pf/callback") == 0) { - if (narg != 6) error->all(FLERR,"Illegal fix external command"); + if (narg != 6) error->all(FLERR,"Incorrect number of args for fix external pf/callback"); mode = PF_CALLBACK; ncall = utils::inumeric(FLERR,arg[4],false,lmp); napply = utils::inumeric(FLERR,arg[5],false,lmp); if (ncall <= 0 || napply <= 0) - error->all(FLERR,"Illegal fix external command"); + error->all(FLERR,"Illegal values for ncall or napply in fix external pf/callback"); } else if (strcmp(arg[3],"pf/array") == 0) { - if (narg != 5) error->all(FLERR,"Illegal fix external command"); + if (narg != 5) error->all(FLERR,"Incorrect number of args for fix external pf/array"); mode = PF_ARRAY; napply = utils::inumeric(FLERR,arg[4],false,lmp); - if (napply <= 0) error->all(FLERR,"Illegal fix external command"); - } else error->all(FLERR,"Illegal fix external command"); + if (napply <= 0) error->all(FLERR,"Illegal value for napply for fix external pf/array"); + } else error->all(FLERR,"Unknown fix external mode {}", arg[3]); callback = nullptr; @@ -323,6 +325,7 @@ void FixExternal::grow_arrays(int nmax) { memory->grow(fexternal,nmax,3,"external:fexternal"); memset(&fexternal[0][0], 0, sizeof(double)*3*nmax); + array_atom = fexternal; } /* ---------------------------------------------------------------------- diff --git a/unittest/c-library/test_library_external.cpp b/unittest/c-library/test_library_external.cpp index 606d53d38b..63ef4eee32 100644 --- a/unittest/c-library/test_library_external.cpp +++ b/unittest/c-library/test_library_external.cpp @@ -123,6 +123,18 @@ TEST(lammps_external, callback) EXPECT_DOUBLE_EQ(reduce[5], 1.4); EXPECT_DOUBLE_EQ(reduce[6], 1.4); + double **fext = + (double **)lammps_extract_fix(handle, "ext", LMP_STYLE_ATOM, LMP_TYPE_ARRAY, 0, 0); + EXPECT_DOUBLE_EQ(fext[0][0], 10.0); + EXPECT_DOUBLE_EQ(fext[0][1], 10.0); + EXPECT_DOUBLE_EQ(fext[0][2], 10.0); + EXPECT_DOUBLE_EQ(fext[3][0], 10.0); + EXPECT_DOUBLE_EQ(fext[3][1], 10.0); + EXPECT_DOUBLE_EQ(fext[3][2], 10.0); + EXPECT_DOUBLE_EQ(fext[7][0], 10.0); + EXPECT_DOUBLE_EQ(fext[7][1], 10.0); + EXPECT_DOUBLE_EQ(fext[7][2], 10.0); + ::testing::internal::CaptureStdout(); lammps_close(handle); output = ::testing::internal::GetCapturedStdout(); @@ -193,6 +205,18 @@ TEST(lammps_external, array) EXPECT_DOUBLE_EQ(pe, 1.0 / 8.0); EXPECT_DOUBLE_EQ(press, 0.15416666666666667); + double **fext = + (double **)lammps_extract_fix(handle, "ext", LMP_STYLE_ATOM, LMP_TYPE_ARRAY, 0, 0); + EXPECT_DOUBLE_EQ(fext[0][0], 6.0); + EXPECT_DOUBLE_EQ(fext[0][1], 6.0); + EXPECT_DOUBLE_EQ(fext[0][2], 6.0); + EXPECT_DOUBLE_EQ(fext[3][0], 6.0); + EXPECT_DOUBLE_EQ(fext[3][1], 6.0); + EXPECT_DOUBLE_EQ(fext[3][2], 6.0); + EXPECT_DOUBLE_EQ(fext[7][0], 6.0); + EXPECT_DOUBLE_EQ(fext[7][1], 6.0); + EXPECT_DOUBLE_EQ(fext[7][2], 6.0); + ::testing::internal::CaptureStdout(); lammps_close(handle); output = ::testing::internal::GetCapturedStdout(); diff --git a/unittest/python/python-fix-external.py b/unittest/python/python-fix-external.py index c061b9f466..55934b8e79 100644 --- a/unittest/python/python-fix-external.py +++ b/unittest/python/python-fix-external.py @@ -1,6 +1,6 @@ import sys,os,unittest from ctypes import * -from lammps import lammps, LMP_STYLE_GLOBAL, LMP_TYPE_VECTOR +from lammps import lammps, LMP_STYLE_GLOBAL, LMP_STYLE_ATOM, LMP_TYPE_VECTOR, LMP_TYPE_ARRAY try: import numpy @@ -92,6 +92,17 @@ class PythonExternal(unittest.TestCase): self.assertAlmostEqual(reduce[5],6.5,14) self.assertAlmostEqual(reduce[6],-0.6,14) + fext = lmp.extract_fix("ext", LMP_STYLE_ATOM, LMP_TYPE_ARRAY) + self.assertAlmostEqual(float(fext[0][0]), 0.0) + self.assertAlmostEqual(float(fext[0][1]), 0.0) + self.assertAlmostEqual(float(fext[0][2]), 0.0) + self.assertAlmostEqual(float(fext[1][0]), 0.0) + self.assertAlmostEqual(float(fext[1][1]), 0.0) + self.assertAlmostEqual(float(fext[1][2]), 0.0) + self.assertAlmostEqual(float(fext[2][0]), 0.0) + self.assertAlmostEqual(float(fext[2][1]), 0.0) + self.assertAlmostEqual(float(fext[2][2]), 0.0) + lmp.command("run 10 post no") self.assertAlmostEqual(lmp.get_thermo("temp"),1.0/30.0,14) self.assertAlmostEqual(lmp.get_thermo("pe"),1.0/8.0,14) @@ -110,6 +121,17 @@ class PythonExternal(unittest.TestCase): val += lmp.extract_fix("ext",LMP_STYLE_GLOBAL,LMP_TYPE_VECTOR,nrow=i) self.assertAlmostEqual(val,15.0,14) + fext = lmp.extract_fix("ext", LMP_STYLE_ATOM, LMP_TYPE_ARRAY) + self.assertAlmostEqual(float(fext[0][0]), 10.0) + self.assertAlmostEqual(float(fext[0][1]), 10.0) + self.assertAlmostEqual(float(fext[0][2]), 10.0) + self.assertAlmostEqual(float(fext[1][0]), 10.0) + self.assertAlmostEqual(float(fext[1][1]), 10.0) + self.assertAlmostEqual(float(fext[1][2]), 10.0) + self.assertAlmostEqual(float(fext[2][0]), 10.0) + self.assertAlmostEqual(float(fext[2][1]), 10.0) + self.assertAlmostEqual(float(fext[2][2]), 10.0) + def testExternalArray(self): """Test fix external from Python with pf/array""" @@ -142,6 +164,16 @@ class PythonExternal(unittest.TestCase): force[i][2] = 0.0 lmp.fix_external_set_energy_global("ext", 0.5) lmp.fix_external_set_virial_global("ext",[0.5, 0.5, 0.5, 0.0, 0.0, 0.0]) + fext = lmp.extract_fix("ext", LMP_STYLE_ATOM, LMP_TYPE_ARRAY) + self.assertAlmostEqual(float(fext[0][0]), 0.0) + self.assertAlmostEqual(float(fext[0][1]), 0.0) + self.assertAlmostEqual(float(fext[0][2]), 0.0) + self.assertAlmostEqual(float(fext[1][0]), 0.0) + self.assertAlmostEqual(float(fext[1][1]), 0.0) + self.assertAlmostEqual(float(fext[1][2]), 0.0) + self.assertAlmostEqual(float(fext[2][0]), 0.0) + self.assertAlmostEqual(float(fext[2][1]), 0.0) + self.assertAlmostEqual(float(fext[2][2]), 0.0) lmp.command("run 5 post no") self.assertAlmostEqual(lmp.get_thermo("temp"),4.0/525.0,14) @@ -170,6 +202,16 @@ class PythonExternal(unittest.TestCase): self.assertEqual(npforce[0][0],6.0) self.assertEqual(npforce[3][1],6.0) self.assertEqual(npforce[7][2],6.0) + fext = lmp.extract_fix("ext", LMP_STYLE_ATOM, LMP_TYPE_ARRAY) + self.assertAlmostEqual(float(fext[0][0]), 6.0) + self.assertAlmostEqual(float(fext[0][1]), 6.0) + self.assertAlmostEqual(float(fext[0][2]), 6.0) + self.assertAlmostEqual(float(fext[1][0]), 6.0) + self.assertAlmostEqual(float(fext[1][1]), 6.0) + self.assertAlmostEqual(float(fext[1][2]), 6.0) + self.assertAlmostEqual(float(fext[2][0]), 6.0) + self.assertAlmostEqual(float(fext[2][1]), 6.0) + self.assertAlmostEqual(float(fext[2][2]), 6.0) ############################## if __name__ == "__main__":