My name is Pablo García-Risueño, and I am writing on behalf of a research project focused on developing novel methods for constrained molecular dynamics. We have recently introduced a new algorithm that, in most cases, outperforms the current state of the art in terms of efficiency and/or accuracy. Moreover, it allows significant increases in the time step of the simulations. The article where we present the algorithm is available at the following link (JCTC):
The article linked above is based on a parallel implementation within the GROMACS package that supports both shared and distributed memory parallelism for CPUs. We are confident that this parallel implementation can be naturally migrated to LAMMPS, making it possible to efficiently impose connected constraints on arbitrarily large sets of atoms (for example, C–C bonds in a protein backbone), overcoming a current limitation of LAMMPS (1).
We are writing to ask whether you might be interested in collaborating in the inclusion of this new algorithm in LAMMPS. We do not have available manpower at the momentfor that, but we would be very happy to support you with any related implementation efforts.
Thank you very much for your time and consideration. With best regards.
(1) In LAMMPS, only small clusters of atoms can be constrained. This is sothe constraint calculation for a cluster can be performed by a single processor, to enable good parallel performance. A cluster is defined as a central atom connected to others in the cluster by constrained bonds. LAMMPS allows for the following kinds of clusters to be constrained: one central atom bonded to 1 or 2 or 3 atoms, or one central atom bonded to 2 others andthe angle between the three atoms alsoconstrained. This means water molecules or CH2 or CH3 groups may be constrained, but not all the C-C backbone bonds of a long polymer chain.
@PabloGR I suggest you post this as a feature request issue on the LAMMPS GitHub project as this will get much more visibility among LAMMPS developers. The LAMMPS core development team is also very thinly stretched as far as manpower goes (plus we all are committed to projects that provide funding), so I cannot make any promises that it will be picked up. Overcoming the known limitations of fix shake/rattle would certainly be a very welcome addition to the codebase (anybody that is reading this and in search of a worthy programming project should take the hint…).
I am not sure if you have looked in detail at the LAMMPS code for fix_shake: it is not the same as GROMACS’s SHAKE code. In effect, it writes and stores the constraint Hessian matrix, which enables quasi-Newton (secant method) iterations with only one matrix inversion per timestep (not per secant iteration). The SHAKE constraints are also iterated over per cluster, rather than constraint-by-constraint.
In my experience this means LAMMPS’s SHAKE timing scales very well with requested precision (although I have not tried to push it to machine precision as you all have!)
This Hessian implementation is an important consideration for extending the SHAKE code to arbitrary clusters, but it can be programmatically generated from the rows of the Jacobian matrix and will be diagonally dominated in most useful cases.
FYI, I have started a Copilot session to create an implementation plan for this new constraint solver for LAMMPS. No idea how well this will work. We made some experiments recently and had moderate successes with porting parts of LAMMPS to the OPENMP and KOKKOS packages and also with augmenting existing code with more complex derivations. The common denominator is that things usually worked quite well, if there was a sufficient template within LAMMPS or the change was often discussed on the web (e.g. in materials related to classes in math or computer science). I don’t expect much, but there also is very little risk and effort involved. I hope that by pointing to the GROMACS implementation and fix shake there is sufficient context to have at least something minimal that can be tested and refined.