How to set the value of per-atom variables via python?

[jasmine969]

After reviewing the time integration algorithm in fix sph (see the SPH document), here are some new insights.

The velocities used to compute forces are the extrapolated velocities (vest in the code and \tilde{\mathbf{v}} in the formula). Although the real velocities are updated in final_integrate, they are not to used to compute forces wherever. However, it should also be noted that the densities are updated in final_integrate, and they are used to compute the pressure in pair/sph/taitwater and pair/sph/taitwater/morris. So the pressure obtained in Verlet::setup will be newer/more correct than those obtained in Verlet::run. And newer pressure will lead to newer forces by

\boldsymbol{f}_i=m_i \frac{\mathrm{d} \boldsymbol{v}_i}{\mathrm{d} t}=-\sum_j m_i m_j\left(\frac{P_i}{\rho_i^2}+\frac{P_j}{\rho_j^2}+\Pi_{i j}\right) \nabla_i W(r_{ij},h)

To conclude, the difference may be due to newer densities and pressure instead of newer velocities.

However, if external force is added, e.g., via fix addforce, I have also found it will lead to incorrect results in Test 2, because fix addforce is carried out in Modify::post_force(), but it will be cleared in Verlet::setup by the force_clear() method. This will make external force vanish in initial_integrate. The remedy might be adding the Modify::post_force method in the Verlet::setup phase. But I don’t think it is too incorrect in Test 3, because I skip the force_clear() method in Verlet::setup, so that the external force is still valid in initial_integrate.

@akohlmey

In my suggestion I had already favored using fix external for its simplicity compared to your approach, but now I am realizing that this is your only option

Thanks for your reminder. I’ve just studied fix external, it seems applicable for my simulation. But fix external only computes a scalar representing the potential energy. I wonder how to output such per-atom force added by fix external to the dump file in an elegant way.