Add support for fix external

[vchuravy]

No description provided.

[codecov-commenter]

Codecov Report

Merging #28 (e47f567) into main (4ea8190) will increase coverage by 14.08%.
The diff coverage is 90.42%.

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

@@             Coverage Diff             @@
##             main      #28       +/-   ##
===========================================
+ Coverage   35.48%   49.56%   +14.08%     
===========================================
  Files           2        3        +1     
  Lines         372      464       +92     
===========================================
+ Hits          132      230       +98     
+ Misses        240      234        -6     

Files	Coverage Δ
src/api.jl	24.71% <100.00%> (+8.81%)	⬆️
src/external.jl	95.45% <95.45%> (ø)
src/LAMMPS.jl	55.90% <76.92%> (+2.84%)	⬆️

[vchuravy]

@akohlmey thank you for your comments and feedback on lammps/lammps#3688

If you have some time I would appreciate your input here. I am pretty sure I am still getting lots of details wrong, but the broad strokes are here.

At least I am now computing "something"

julia> forces = extract_atom(lmp, "f")
3×10 Matrix{Float64}:
 0.0   0.174073   0.0  0.0  0.0  0.0  0.0  0.0077713  0.0  -0.181845
 0.0  -0.138762   0.0  0.0  0.0  0.0  0.0  0.0496266  0.0   0.0891358
 0.0   0.0863738  0.0  0.0  0.0  0.0  0.0  0.069386   0.0  -0.15576

[vchuravy]

Is this necessary? At least it allows me to use += later.

[akohlmey]

Yes and no. The expectation for the design of fix external is that it is used to couple some other software that can compute forces to LAMMPS, feed the atom positions and (global) indexes to it and let it compute the forces. Those are then copied to the fexternal array. Typically the program would have its own input file(s) and read the topology and initial geometry from that and then the data from fix external is used to update the positions (hence the atom-IDs since the order in the local arrays can change all the time) and then compute the forces.

So if you collect your forces in a buffer of your own, you only need to copy them. If you want to collect them directly into the array passed from LAMMPS you need to zero it out first.

Another item to take care of are forces between pairs of atoms that straddle subdomain boundaries (or periodic domain boundaries in case you run in serial). You cannot return forces on "ghost" atoms so you either have to use a full neighbor list (have each pair listed twice) or use "newton off" so that pairs across domain boundaries are listed twice and then store forces only with local atoms (hence the passing of nlocal). The alternative would be to implement some communication.

At this point, it is probably a good idea to read through sections 4.1 to 4.6 here: https://docs.lammps.org/Developer.html
Fix external is called at the "post_force" step.

[vchuravy]

Sofar I have been following the pair_python code

https://github.com/lammps/lammps/blob/730e5d2e64106f3e5357fd739b44c7eec19c7d2a/src/PYTHON/pair_python.cpp#L189-L202

But I am unsure how to obtain eflag, evflag and special_lj as an external fix.

[akohlmey]

Sofar I have been following the pair_python code

Please note that in the flow of control, a post_force fix will be run after the reverse communication. So you must not store forces on ghost atoms. That information is lost. You either have to require newton_pair off or use a full neighbor list and store forces only on local atoms. Or you have to implement your own reverse communication (not recommended).

But I am unsure how to obtain eflag, evflag and special_lj as an external fix.

For starters, you can always compute energy and virial contribution and use the corresponding functions in the library interface to set them. At a later point we can look into an optimization to make that optional. It is really only a small performance optimization to not compute and tally those, if they are not required. You cannot use the "ev_tally" function, but there are library functions that you can call instead. https://docs.lammps.org/Library_utility.html#_CPPv437lammps_fix_external_set_energy_globalPvPKcd

Special_lj and special_coul are a different story. I would have to add access to those to the library interface so you can request those via lammps_extract_global() https://docs.lammps.org/Library_properties.html#_CPPv421lammps_extract_globalPvPKc

I will add that change to lammps/lammps#3687

[akohlmey]

lammps/lammps@fa9062a
adds access to the special factor arrays.