I’m trying to implement the multi-state empirical valence bond theory method using LAMMPS and its Python API. I understand that it seems to already have an implementation in the form of RAPTOR, but our research group wants to have an easily modifiable/usable version (hence the Python).
I have a working version of the code (with only a few bugs) that I can send upon request if someone would like a look. Briefly, my current implementation is as follows:
- Create one main lammps object that is used to calculate the main potential energy and integrate the equations of motions/thermostat with any necessary fixes.
- After each step, a once-off new lammps object is created per EVB state that is identified, with this lammps object reused in subsequent steps. These new lammps objects are created by reading a restart file dumped by the main lammps object (faster way I could see to copy across all the topological information from one object to another).
- The topology of the new lammps object is changed through “set” and “create_bonds” commands to match the identified EVB state, and then the “run 0” command is used to allow the computation of the potential energy. From what I can tell, calling run 0 is quite slow assuming because some neighbourlists are changed/updated.
- With each subsequent step, given no reaction has happened yet, the EVB-state lammps objects have their positions and box data changed from the “scatter_atoms” method and the “change_box” command. If a reaction occurs, the corresponding lammps object with the new identified state writes a restart file which is read by all other lammps objects (including the main).
With that explained, my question is if there is a more efficient way to do this. The code runs quite slowly (around 1 ns/day for a 25 Angstrom box of water with a hydronium ion in it). The bottleneck seems to be after calling run 0 for each EVB state, so I was wondering if there is some command I have missed that can calculate the potential energy once-off without rebuilding neighbourlists/preparing the object for a simulation. I have also attempted some parallelisation using pylammpsmpi on github (GitHub - pyiron/pylammpsmpi: Parallel Lammps Python interface - control a mpi4py parallel LAMMPS instance from a serial python process or an Jupyter notebook) with little success. It doesn’t seem to like having multiple lammps objects active at one time, and I can’t imagine parallelisation will help with the current bottleneck.
Any tips would be much appreciated.