ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

get access to the applied external force from the fix

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-11-24 22:22:13 -05:00
parent 6c333837e0
commit 3f78ee72c6
4 changed files with 88 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
doc/src/fix_external.rst
View File

@ -173,9 +173,17 @@ stress/atom <compute_stress_atom>` commands. The former can be
accessed by :doc:`thermodynamic output <thermo_style>`. The default
setting for this fix is :doc:`fix_modify virial yes <fix_modify>`.
This fix computes a global scalar which can be accessed by various
:doc:`output commands <Howto_output>`. 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 <Howto_output>`.
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()
<lammps_fix_external_set_vector>` and
:cpp:func:`lammps_fix_external_set_vector_length()
<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 <run>` command.

17
src/fix_external.cpp
View File

@ -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;
}
/* ----------------------------------------------------------------------

24
unittest/c-library/test_library_external.cpp
View File

@ -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();

44
unittest/python/python-fix-external.py
View File

@ -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__":