Is it possible to use lammps to perform dynamics with forces provided by an external (python) code?
If so, could you please suggest a viable strategy to implement it?
I have a machine learning code that provides the cartesian components of the force acting on each atom.
I am not interested in defining pair or bond terms.
Most of the computer time will be used by the machine learning code, so efficiency for the lammps part is not an issue. In other words, I’d rather have a simple solution than an efficient one.
Here is what I tried:
A) according to the manual, python codes can set string variables in lammps, so for N atoms I declare 3N ‘string’ variables of the form:
variable sf1x string 0.0 # string to contain fx for atom 1
and 3N ‘equal’ variables:
variable f1x equal v_sf1x # value of fx for atom 1
then I create N groups, one for each atom:
group a1 id 1
and activate a fix to set the force components
fix set1 a1 setforce v_f1x v_f1y v_f1z
next I declare a ‘python’ function to compute the forces:
python forces input 1 SELF format p file set_forces.py
here, the code in set_forces.py extracts the coordinates from lammps, computes the forces with a neural network, and sets the string variables with lines of the form:
lmp.set_variable(“sf1x”,"g" F[0,0])
the dynamics is then executed computing the forces at every step:
run 100 every 1 “python forces invoke”
This approach works only up to 31 atoms due to the MAX_GROUP limit of 32.
I tried recompiling after setting MAX_GROUP to a higher value, but groups past the 32nd are not properly constructed and contain too many atoms.
Is there a way to increase the number of groups, or to apply fixes directly to atoms?
B) If the forces are in vector f, then the python command
lmp.scatter_atoms(“f”,1,3,f)
will change the force values directly, i.e., without passing through the ‘string’ and ‘equal’ variables, the groups and the fixes.
Forces, however, are zeroed and recomputed inside lammps prior to each step, so this route does not work, at least without tampering with the ‘run’ command, which is probably beyond my capabilities.
Is there a simple way to interface python directly with a fix->post_force() function?
Thank you for any suggestion.
Cecco