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

[jasmine969]

Hi LAMMPS developers and users,

With Python lammps module, I’ve known that the atom-style variable values can be obtained by extract_variable method. The question is how to set the atom-style variable values, for example, by passing the numpy ndarray to lammps.

This question is a continuation of Problem 2 in Loop cannot be run normally in python lammps module - LAMMPS / LAMMPS General Discussion - Materials Science Community Discourse. In details, I am running a simulation which has to be coupled with a complex python function, say mypyfunc. For each timestep, mypyfunc receives the input from lammps simulated results, v_res, and computed the external force exerted on each atom, fext, and passes fext to lammps. For the time being, I don’t know how to set the per-atom fext via the Python lammps module. The expedient is that I save fext to a csv file, and let lammps read this file, assign fext to an atomfile-style variable as follows:

from mypyfunc import mypyfunc, save_to_csv
from lammps import *

lmp = lammps()
lmp.file('initialize.in')
for i in range(Nt):
    lmp.command('run 1')
    v_res = np.array([lmp.numpy.extract_fix(
        'res', LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY, nrow=nrow, ncol=2
    ) for nrow in range(50)])
    fext = mypyfunc(v_res)
    save_to_csv(fext, 'fext_current.csv')
    lmp.command('variable fext delete')
    lmp.command('variable fext atomfile fext_current.csv')

I think this is very awkward and time-consuming, in that I must delete fext, re-create it, read a file from the disk, parse it, and re-assign it to fext. As @akohlmey suggested in Loop cannot be run normally in python lammps module - LAMMPS / LAMMPS General Discussion - Materials Science Community Discourse

The only way I can think to speed this up would be to not read forces from files, but compute them on-the-fly.

In practice, they are indeed computed on the fly, merely via Python. LAMMPS has very good support for Python programming, in terms of whether the Python lammps/pylammps module or the embedded LAMMPS PYTHON package. So I believe there must a way to set this per-atom force on the fly, rather than save it to a file and read it. But I haven’t found this method yet. The methods to set variables are only set_string_variable and set_internal_variable, but they are meant for string-style and internal-style variables and cannot satisfy my needs.

Any commens or advice will be appreciated!

[akohlmey]

You cannot set the value of atom style variables (same as you cannot set the value of equal style variables"). What you can set is a custom per-atom property created by an instance of fix property/atom. Into that storage you can store data with the “scatter()” method. Unfortunately, this method has still not been refactored to have embedded docstrings, its documentation is still a comment in the python code:

  # scatter vector of atom/compute/fix properties across procs
  # 2 variants to match src/library.cpp
  # name = atom property recognized by LAMMPS in atom->extract()
  # type = 0 for integer values, 1 for double values
  # count = number of per-atom valus, 1 for type or charge, 3 for x or f
  # assume data is of correct type and length, as created by gather_atoms()
  # NOTE: need to ensure are converting to/from correct Python type
  #   e.g. for Python list or NumPy or ctypes

  def scatter(self,name,dtype,count,data):
    if name: newname = name.encode()
    else: newname = None
    with ExceptionCheck(self):
      self.lib.lammps_scatter(self.lmp,newname,dtype,count,data)

Where details are in the C-library interface docs at: 1.1.6. Scatter/gather operations — LAMMPS documentation

This custom property can be accessed in a per-atom variable and then used as desired.

An alternate and probably better approach could be to use fix external with a suitable callback function to add custom forces to atoms directly. This avoids having to use variables and global operations altogether. The callback may be set from Python and be a Python function: 2.4. The lammps Python module — LAMMPS documentation
This callback can operate on NumPy arrays and thus it may be quite fast even though using Python, if programmed properly.

[Update:] If you need some inspiration, there is a python unittest program in unittest/python/python-fix-external.py which tries to use as many fix external features as possible.

[jasmine969]

Thanks for your detailed reply. I will study the scatter method and fix external carefully.

[srtee]

Additionally, you should consider if it’s feasible to:

1. Write a new fix in C++ that specifically implements this added force in LAMMPS, or
2. Migrate to a primarily Python-based simulation framework, such as OpenMM or ASE, where Numpy arrays can be used “natively”.

This is not just motivated by performance concerns. From a reproducibility standpoint, it is easier to reproduce a paper if the scientific software in it relies primarily on one programming environment – either C++, or Python. Relying on two environments means an interested researcher must set up both environments correctly, and their interactions, to replicate your results.

This is not purely an ideological concern either (although surely scientific reproducibility is a high priority!) – when a scientific paper is more easily reproduced it can inspire more follow-up studies and be more widely cited, which is good for individual scientific careers.

LAMMPS “mods” get more traction in the materials modelling community than either modifications of other software or new packages, and one reason is because LAMMPS’s core code deliberately makes it easy for other C++ code to “plug in” where needed, and then anybody who wants to replicate some esoteric technique just needs a C++ compiler (which you need for any serious scientific calculation anyway).

[jasmine969]

Thanks for your comments. Actually I’m using LAMMPS due to its SPH package (which is not owned by other MD software) together with some molecular dynamics toolkits. But perhaps it is better to write the mypyfunc with C++.

[jasmine969]

I tried the scatter_atoms (same as scatter for peratom properties) these days. To accelerate the computation, I use run 1 pre no instead of the original run 1. Differences were observed. This is interesting.

To test whether the differences were owing to the scatter_atoms, I switched to a simpler cavity_flow case in examples/PACKAGES/sph/cavity_flow. Three tests were carried out:
Test 1. run 10000
Test 2.

for _ in range(10000):
    lmp.command('run 1')

Test 3.

for _ in range(10000):
    lmp.command('run 1 pre no')

The results of Test 1 and Test 3 are exactly the same, but are different from those of Test 2.

I guess the reason may be that Test 2 would do one more force computation before running 1 timestep (in void Verlet::setup(int flag)). For many molecular dynamics problems, this will make no difference because the force is only the function of position, like the LJ potential. But in SPH, force is not only related to the position, but also to the velocity when computing the artificial viscous force. In final_integrate, the velocity will be updated. So the force computed in Verlet::setup after the last final_integrate can be different from that computed in Verlet::run before the last final_integrate.

So, can I say Test 2 should be more accurate because it carries out one more force computation after final_integrate? And if this is true, there seems a need to improve the algorithm for those simulations with forces correlated with velocities, so that the force can be updated once more.

[akohlmey]

@jasmine969

I guess the reason may be that

You should never guess in research. Rather you need to apply the scientific method of forming a hypothesis and then devising suitable experiments to confirm or disprove its validity and iterate until consistency.

@jasmine969

can I say Test 2 should be more accurate because it carries out one more force computation after final_integrate?

No. Please recall that LAMMPS uses the Velocity Verlet algorithm, where you first do a half step updating the velocities with the “old” forces pertaining to the current positions, then do a full step updating the positions, then recompute the forces, and then do another half step updating the velocities with the “new” forces for the updated positions.

In your Test 2 you use your externally computed force contribution, which is based on the “old” positions and “old” velocities, twice and it is not correct in both cases since in the first case the velocities are incorrect (that is a general problem for any pair style that uses forces dependent on velocities, i.e. also applies to DPD, but only is a “small” problem since it happens only during the integrate::setup phase and not during the regular time stepping) and in the second case both positions and velocities are not correct.
This makes me realize that your approach to computing external forces was flawed from the beginning (you didn’t mention that you were using SPH) and cannot be done correctly this way, regardless of whether you are using run with pre yes or pre no.

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, short of writing a few fix style or creating a modified pair style yourself, that will result in correctly computed and applied forces during a run. You must compute and apply your forces during the Force::compute() or Modify::post_force() step and in the latter case also make certain that there is no fix modifying velocities like fix temp/berendsen or fix temp/rescale (which are not recommended to use in the first place) that is defined before your instance of the force computation fix (i.e. external or custom fix).

[jasmine969]

Please note that the three tests were carried out using the cavity-flow case, exactly the same as that in examples/PACKAGES/sph/cavity_flow, where no external force is exerted.

from lammps import lammps

lmp = lammps()
lmp.file('examples/PACKAGES/sph/cavity_flow/cavity_flow.lmp')
lmp.command('run 10000')  # Test 1
# for _ in range(10000):
    # lmp.command('run 1')  # Test 2
    # lmp.command('run 1 pre no')  # Test 3

Just run the above code. No external force is computed and added

[akohlmey]

@jasmine969

where no external force is exerted.

That does not contradict my statements and conclusions.
Without external force, both Test 1 and Test 3 are correct.
But with external force neither variant will be correct.

[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.

[akohlmey]

@jasmine969

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.

Why would you want to have LAMMPS output that force, since you are computing it yourself?
You seem to be misunderstanding how fix external works. The callback function that you need to register, is being passed multiple arguments. In python this is (as shown in the unittest file I mentioned):

def callback_one(lmp, ntimestep, nlocal, tag, x, f):

• the lammps class instance
• the time step numnber
• the number of local atoms
• the positions of the local atoms as list of list indexed by atom index and x,y,z
• a list of lists for added forces with the same indexing

If you need additional properties to compute your forces, you can access them through the LAMMPS class instance. Fix external then calls the callback function in Modify::post_force to obtain the force to be added to the per-atom forces. You can also set contributions to energy and virial from the external force.

[srtee]

@jasmine969

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.

@akohlmey may know better, but I don’t believe this is correct based on the source code.

Verlet::setup (including minimal setup) calls modify->setup(vflag) in its second-last line.

This is after the force-clear (and “normal” force calculations via pair, bond, etc), so any changes at that point persist into the first initial_integrate.

Modify::setup(vflag) in turn calls all fix->setup(vflag) functions.

FixAddForce::setup(vflag) explicitly adds its customised force into the force arrays, by calling (FixAddForce’s) post_force(vflag) function.

Therefore, fix addforce’s added force is applied from the very initial forces in step 0 (albeit possibly having different calculated values due to any other information differences during setup vs post_force).

For such a command to only apply its force during step 1 and not step 0 would seem extremely counterintuitive to anyone doing fancy molecular dynamics and would indeed require immediate highlighting and massive debugging efforts (which is why I am confident that’s not how it works).

[jasmine969]

I want to output the external force to the dump file simply because I’m going to visualize it in OVITO during post-processing, so that I can do some analysis. For sure, the force has been computed by Python and I can obtain it. And there is no convenient way to add it to the list of attributes to be dumped through Python lammps module.

As you can see in the following three ways I have used, I can obtain this extenal force via Python in all three ways. But only when using fix external, I cannot write the force to the dump file by the dump command (for visualization in OVITO).

When I use the atomfile-style variable, delete and re-assign it via Python, I can dump it as follows:

// LAMMPS
variable FP atomfile force.csv
variable Fx atom -v_FP*v_cosQ
variable Fy atom -v_FP*v_sinQ

fix      1 wall addforce v_Fx v_Fy 0
dump     my_dump all custom 1 xxx.dump id type &
         x y z v_FP

# Python
v_res = lmp.extract_fix(...)
FP = mypyfunc(v_res)
save_to_csv(FP, 'force.csv')
lmp.command(f'variable FP delete')
lmp.command(f'variable FP atomfile force.csv')
lmp.command('run 1')

When I use the per-atom property added by fix property/atom, scatter it via Python, I can dump it as follows:

// LAMMPS
fix      FP wall property/atom d_FP
variable Fx atom -d_FP*v_cosQ
variable Fy atom -d_FP*v_sinQ

fix      1 wall addforce v_Fx v_Fy 0
dump     my_dump all custom 1 xxx.dump id type &
         x y z d_FP

# Python
v_res = lmp.extract_fix(...)
FP = mypyfunc(v_res)
d_FP = (n_partices * c_double)(*FP)
lmp.scatter_atoms('d_FP', 1, 1, d_FP)

But when I use fix external,

fix 1    all external pf/callback 1 1
// fix addforce is not needed, but how to dump the external force?
dump     my_dump all custom 1 xxx.dump id type &
         x y z ???

# Python
def callback(caller, ntimestep, nlocal, tag, x, fext):
    v_res = lmp.extract_fix(...)
    FP = mypyfunc(v_res)
    cosQ = lmp.extract_variable(...)
    sinQ = lmp.extract_variable(...)
    fext[:, 0] = FP * cosQ
    fext[:, 1] = FP * sinQ

Another question is whether the fext in the callback function is sorted by the global ID or the local ID. This matters since cosQ and sinQ I extract are variables sorted by the local ID, while FP is computed by v_res and is sorted by global ID.

[jasmine969]

That’s true. I check the source code once again and find modify->setup(...) after pair->compute(...) in Verlet::setup(...). Sorry for my carelessness.

[akohlmey]

@jasmine969

But only when using fix external, I cannot write the force to the dump file by the dump command (for visualization in OVITO).

First of all, I think OVITO can add data to a trajectory read from an external file via a Python modifier, so it should not be necessary. But then again, it should also be possible to use the previously mentioned fix property/atom to have a per-atom property added to which you can set/add the computed force from within your callback and then output it. Also, it would be rather straightforward to give access to the added force for fix external by adding a few lines of C++ code in the right places to fix_external.cpp.

@jasmine969

Another question is whether the fext in the callback function is sorted by the global ID or the local ID

Local index. That is implied by the fix handing nlocal to the callback function. Only global data read from files has a global ID, but those usually do not need to be sorted by ID.

[srtee]

@jasmine969

But when I use fix external,

You may be able to add two separate fix store/force, before and after the fix external, and then write per-atom variables that subtract the before from after to give you the added force as LAMMPS-native per-atom variables.

[akohlmey]

@akohlmey

Also, it would be rather straightforward to give access to the added force for fix external by adding a few lines of C++ code in the right places to fix_external.cpp.

I have done just that here: Collected small changes and fixes by akohlmey · Pull Request #4389 · lammps/lammps · GitHub

[jasmine969]

@srtee @akohlmey
Thanks for your kind suggestions and improvements. I haven’t tried Axel’s updated fix external yet. Now I have compared the computation efficiencies and the result values using the original fix external command with the aforementioned other ways. When I use the original fix external, I output the external force by adding a new peratom property by fix property/atom and scatter it in the callback function.

Test 1 used the atomfile-style variable and read data from a csv file. In Python, lmp.command('run 1') is executed repeatedly in a for-loop.

result1  		2	39.10
result2 		0	-10.88
result3 		2	2.00e-02
barycenter/m	-0.0075 (y) 0.1414 (z)
Time elapsed	00:02:22	LAMMPS for 00:00:33	mypyfunc for 00:01:48
ETA				00:57:15
Step (curr -- total): 200/5000 -- 200/5000

Test 2 used fix property/atom and scattered it before lmp.command('run 1')

result1 		2	39.10
result2     	0	-10.88
result3 		2	2.00e-02
barycenter/m	-0.0075 (y) 0.1414 (z)
Time elapsed	00:02:16	LAMMPS for 00:00:31	mypyfunc for 00:01:44
ETA				00:54:47
Step (curr -- total): 200/5000 -- 200/5000

As you can see, the results of Test 1 and Test 2 are exactly the same.

Test 3 was nearly the same as Test 2. The only modification is that here lmp.command('run 1 pre no') is executed in the for-loop, rather than the default pre yes.

result1 		2	39.25
result2 		0	-10.88
result3 		2	2.00e-02
barycenter/m	-0.0075 (y) 0.1414 (z)
Time elapsed	00:02:02	LAMMPS for 00:00:17	mypyfunc for 00:01:45
ETA				00:49:18
Step (curr -- total): 200/5000 -- 200/5000

Result 1 (39.25) is a little different from Test 1 & 2 (39.10), but the other results are exactly the same. Meanwhile, the time elasped is greatly reduced, from over 30 sec to merely 17 sec.

Test 4 used fix external in the LAMMPS script. In Python script, no for-loop was needed. Instead, I used lmp.command('run 5000') directly. This is a fluid-structure interaction simulation. I used mypyfunc to compute the forces exerted on the wall. The callback function is like

def callback(caller, step, nlocal, tag, x, fext):
    ave_stretch = np.array([caller.numpy.extract_fix(
    'ave_stretch', LMP_STYLE_GLOBAL, LMP_TYPE_ARRAY, nrow=nrow, ncol=2
    ) for nrow in range(n_sense)])
    FP = mypyfunc(ave_stretch)
    cosQ = caller.extract_variable('cosQ', group='wall', vartype=LMP_VAR_ATOM)
    sinQ = caller.extract_variable('sinQ', group='wall', vartype=LMP_VAR_ATOM)
    FP_per_atom = np.r_[-np.repeat(FP, n_per_ring), np.zeros(n_fluid)]
    caller.scatter_atoms('d_F_active_mag', 1, 1, (n_partices * c_double)(*FP_per_atom))
    FP_per_atom = FP_per_atom[tag - 1] # transform global FP to local FP
    fext[:, 0] = FP_per_atom * cosQ
    fext[:, 1] = FP_per_atom * sinQ

The simulation results were

result1 		5	40.80
result2 		2	-10.68
result3 		3	2.00e-02
barycenter/m	-0.0074 (y) 0.1414 (z)
Time elapsed	00:01:53	LAMMPS for 00:00:14	mypyfunc for 00:01:45
ETA				00:47:57
Step (curr -- total): 200/5000 -- 200/5000

Though it was much more computationally efficient, the results are quite different from the previous three tests. Compared with the other tests which require a for-loop, the sole change by introducing fix external I could come up with is that the external force is calculated during modify->post_force, rather than one timestep is finished. More investigation into the codes needs to be carried out to demonstrate which of the four approaches is more correct and reliable.