New constraint solver

Dear LAMMPS community

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):

https://urldefense.com/v3/__https://pubs.acs.org/doi/10.1021/acs.jctc.5c01376__;!!D9dNQwwGXtA!UGpcLyle_O3Sc7o5DaTvU3Pn-R4C56AwrNF8PsLAhTkxZVyxELd6SDpBnSrY1IeIACqY0yv1OVRLlXAFUF4d20hzu1JC$

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.

1 Like

@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…).

Hi @PabloGR

Nice to see you here!

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.

Happy to chat more at [email protected] .

1 Like