From aa885a9d8df4a432fdc1141fc73934412ff0f0c6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 21 Jul 2021 17:06:21 -0400 Subject: [PATCH] make virial processing use the total global virial --- src/fix_external.cpp | 37 +++++++++++++++++++++++++++-------- src/library.cpp | 46 ++++++++++++++++++++++++++++++-------------- 2 files changed, 61 insertions(+), 22 deletions(-) diff --git a/src/fix_external.cpp b/src/fix_external.cpp index 99539e86cf..f95c30c28d 100644 --- a/src/fix_external.cpp +++ b/src/fix_external.cpp @@ -15,6 +15,7 @@ #include "fix_external.h" #include "atom.h" +#include "comm.h" #include "error.h" #include "memory.h" #include "update.h" @@ -185,10 +186,13 @@ void FixExternal::min_post_force(int vflag) /* ---------------------------------------------------------------------- caller invokes this method to set its contribution to global energy + this is the *total* energy across all MPI ranks of the external code + and must be set for all MPI ranks. unlike other energy/virial set methods: do not just return if eflag_global is not set b/c input script could access this quantity via compute_scalar() even if eflag is not set on a particular timestep + this function is compatible with CALLBACK and ARRAY mode ------------------------------------------------------------------------- */ void FixExternal::set_energy_global(double caller_energy) @@ -197,22 +201,32 @@ void FixExternal::set_energy_global(double caller_energy) } /* ---------------------------------------------------------------------- - caller invokes this method to set its contribution to global virial + caller invokes this method to set its contribution to the global virial + for all MPI ranks. the virial value is the *total* contribution across + all MPI ranks of the external code and thus we need to divide by the + number of MPI ranks since the tallying code expects per MPI rank contributions. + this function is compatible with PF_CALLBACK and PF_ARRAY mode ------------------------------------------------------------------------- */ void FixExternal::set_virial_global(double *caller_virial) { + const double npscale = 1.0/(double)comm->nprocs; for (int i = 0; i < 6; i++) - user_virial[i] = caller_virial[i]; + user_virial[i] = npscale * caller_virial[i]; } /* ---------------------------------------------------------------------- - caller invokes this method to set its contribution to peratom energy + caller invokes this method to set its contribution to peratom energy. + this is applied to the *local* atoms only. + this function is compatible with PF_CALLBACK mode only since it tallies + its energy contributions directly into the accumulator arrays. ------------------------------------------------------------------------- */ void FixExternal::set_energy_peratom(double *caller_energy) { if (!eflag_atom) return; + if ((mode == PF_ARRAY) && (comm->me == 0)) + error->warning(FLERR,"Can only set energy/atom for fix external in pf/callback mode"); int nlocal = atom->nlocal; for (int i = 0; i < nlocal; i++) @@ -221,6 +235,9 @@ void FixExternal::set_energy_peratom(double *caller_energy) /* ---------------------------------------------------------------------- caller invokes this method to set its contribution to peratom virial + this is applied to the *local* atoms only. + this function is compatible with PF_CALLBACK mode only since it tallies + its virial contributions directly into the accumulator arrays. ------------------------------------------------------------------------- */ void FixExternal::set_virial_peratom(double **caller_virial) @@ -228,6 +245,8 @@ void FixExternal::set_virial_peratom(double **caller_virial) int i,j; if (!vflag_atom) return; + if ((mode == PF_ARRAY) && (comm->me == 0)) + error->warning(FLERR,"Can only set virial/atom for fix external in pf/callback mode"); int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) @@ -236,8 +255,8 @@ void FixExternal::set_virial_peratom(double **caller_virial) } /* ---------------------------------------------------------------------- - caller invokes this method to set length of vector of values - assume all vector values are extensive, could make this an option + caller invokes this method to set length of global vector of values + assume all vector values are extensive. ------------------------------------------------------------------------- */ void FixExternal::set_vector_length(int n) @@ -252,14 +271,15 @@ void FixExternal::set_vector_length(int n) } /* ---------------------------------------------------------------------- - caller invokes this method to set Index value in vector - index ranges from 1 to N inclusive + caller invokes this method to set value for item at "index" in vector + index is 1-based, thus index ranges from 1 to N inclusively. + Must be called from all MPI ranks. ------------------------------------------------------------------------- */ void FixExternal::set_vector(int index, double value) { if (index > size_vector) - error->all(FLERR,"Invalid set_vector index in fix external"); + error->all(FLERR,"Invalid set_vector index ({} of {}) in fix external",index,size_vector); caller_vector[index-1] = value; } @@ -290,6 +310,7 @@ double FixExternal::compute_vector(int n) double FixExternal::memory_usage() { double bytes = 3*atom->nmax * sizeof(double); + bytes += 6*sizeof(double); return bytes; } diff --git a/src/library.cpp b/src/library.cpp index a6e4509c56..0c897538a1 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -4806,8 +4806,8 @@ void lammps_decode_image_flags(imageint image, int *flags) /* ---------------------------------------------------------------------- */ -/** Set the callback function for a fix external instance with the given ID. - Optionally also set the pointer to the calling object. +/** Set up the callback function for a fix external instance with the given ID. + \verbatim embed:rst Fix :doc:`external ` allows programs that are running LAMMPS through @@ -4821,7 +4821,15 @@ mode. The function has to have C language bindings with the prototype: void func(void *ptr, bigint timestep, int nlocal, tagint *ids, double **x, double **fexternal); -This is an alternative to the array mechanism set up by :cpp:func:`lammps_fix_external_get_force`. +The argument *ptr* to this function will be stored in fix external and +the passed as the first argument calling the callback function `func()`. +This would usually be a pointer to the LAMMPS instance, i.e. the same +pointer as the *handle* argument. This would be needed to call +functions that set the global or per-atom energy or virial +contributions. + +The callback mechanism is on of the two modes of fix external. The +alternative is the array mode set up by :cpp:func:`lammps_fix_external_get_force`. Please see the documentation for :doc:`fix external ` for more information about how to use the fix and how to couple it with an @@ -4865,16 +4873,24 @@ its library interface to add or modify certain LAMMPS properties on specific timesteps, similar to the way other fixes do. This function provides access to the per-atom force storage in the fix -to be added to the individual atoms when using the "pf/array" mode. The -*fexternal* array can be accessed similar to the "native" per-atom -*arrays accessible via the :cpp:func:`lammps_extract_atom` function. -Because the underlying data structures can change as atoms migrate -between MPI processes, this function should be always called immediately -before the forces are going to be set. +external instance to be added to the individual atoms when using the +"pf/array" mode. The *fexternal* array can be accessed similar to the +"native" per-atom *arrays accessible via the +:cpp:func:`lammps_extract_atom` function. Please note that the array +stores forces for *local* atoms, in the order determined by the neighbor +list build. Because the underlying data structures can change as well as +the order of atom as they migrate between MPI processes, this function +should be always called immediately before the forces are going to be +set to get an up-to-date pointer. You can use +e.g. :cpp:func:`lammps_get_natoms` to obtain the number of local atoms +and thus the dimensions of the returned force array (``double force[nlocal][3]``). -This is an alternative to the callback mechanism set up by -:cpp:func:`lammps_set_fix_external_callback` with out using the callback -mechanism to call out to the external program. +This is an alternative to the callback mechanism in fix external set up by +:cpp:func:`lammps_set_fix_external_callback`. The main difference is +that this mechanism can be used when forces can be pre-computed and the +control alternates between LAMMPS and the external command, while the +callback mechanism can call the external code to compute the force when +the fix is triggered and needs them. Please see the documentation for :doc:`fix external ` for more information about how to use the fix and how to couple it with an @@ -4914,9 +4930,11 @@ double **lammps_fix_external_get_force(void *handle, const char *id) This is a companion function to :cpp:func:`lammps_set_fix_external_callback` and :cpp:func:`lammps_fix_external_get_force` to also set the contribution -to the global energy from the external code. The value of the *eng* +to the global energy from the external code. The value of the *eng* argument will be stored in the fix and applied on the current and all -following timesteps until changed. +following timesteps until changed by another call to this function. +When running in parallel, the value is the per-MPI process contribution, +not the total energy. Please see the documentation for :doc:`fix external ` for more information about how to use the fix and how to couple it with an