EAM depends on precomputed data and that will work for compute group/group if used within a run and it may or may not work for compute pair/local but it will certainly not work for pair_write. Having this inconsistency and potential segmentation faults is why I called this a very bad idea™.
Effectively, you have a model that is not really compatible with the single function.
Hello Dr Axel,
I just want to inquire one facet about that you mentioned a few months back. I am very close to getting this to work and I want to calculate my energies to check if everything is a-ok… I have been able to implement intergranular potential energy calculation and frictional dissipation calculation but my final bottleneck is kinetic energy. As you mentioned in the past, I might have to rescale my density accordingly so that the mass is constant…
Where exactly must I instantiate this? I dumped the mass via the dump file and noticed that the radius changes but the mass is still constant (same as that of timestep t=0). I also checked fix_nve.cpp and noticed that they called mass and not density*spherical_volume so I’m guessing it won’t affect time integration.
Is your recommendation just for good measure because other functions might call mass as rho*V?.. and if so can I just include it at the end of the pair potential?
Thanks and hopefully this will be the last time I trouble you! 
Wherever you change the radius.
To explain:
When you use atom style sphere, the mass of each atom is no longer taken as the mass value associated with the atom type, but computed from density and radius and stored in the Atom::rmass array. In those cases Atom::rmass_flag is 1. The rmass value is set when the radius is set or changed either when reading a data file or when using the set command or when the radius is changed.
lammps]$ grep -n rmass\\[.*=.*radius src/*.cpp src/*/*.cpp
src/atom_vec_line.cpp:405: rmass[ilocal] = 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
src/atom_vec_line.cpp:430: rmass[ilocal] *= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
src/atom_vec_line.cpp:455: rmass[ilocal] /= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
src/atom_vec_sphere.cpp:128: if (radius_one > 0.0) rmass[ilocal] *= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
src/atom_vec_sphere.cpp:148: rmass[ilocal] = rmass_one / (4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one);
src/atom_vec_tri.cpp:622: rmass[ilocal] = 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
src/atom_vec_tri.cpp:647: rmass[ilocal] *= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
src/atom_vec_tri.cpp:675: rmass[ilocal] /= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
src/set.cpp:898: atom->rmass[i] = MY_PI*atom->radius[i]*atom->radius[i] * dvalue;
src/BPM/atom_vec_bpm_sphere.cpp:207: if (radius_one > 0.0) rmass[ilocal] *= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
src/BPM/atom_vec_bpm_sphere.cpp:237: rmass[ilocal] = rmass_one / (4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one);
src/FEP/fix_adapt_fep.cpp:532: rmass[i] = 4.0*MY_PI/3.0 * radius[i]*radius[i]*radius[i] * density;
src/KOKKOS/atom_vec_sphere_kokkos.cpp:2533: h_rmass[nlocal] = 4.0*MY_PI/3.0 * h_radius[nlocal]*h_radius[nlocal]*h_radius[nlocal];
The same, of course, applies if you would change the density.
The code in atom_vec_sphere.cpp is a bit confusing, you need to understand that in the data file you provide the diameter and the density. Those are stored in the radius and rmass fields. The AtomVecSphere::data_atom_post() function then converts those to radius and mass.
So if you change the radius without updating the rmass field you implicitly also change the density so that the mass remains constant. If that is what you want, you will be OK, but should document this behavior.
Thanks Dr Axel.
This is what I want so I will keep it as is, and probably explains why the mass is constant in my dump file without making any alterations to the density.
Thanks!