I am studying the code for “compute born/matrix” and focusing on the improper contribution. It seems to me that there is no method to calculate this contribution when a command like
“compute ID group-ID born/matrix improper”
is specified in the input script file, since there is no “compute_impropers()” to be called (or another action taking place) when the impflag is set to True.
I would greatly appreciate if someone could comment/guide me on that.

as the writer of the compute I think I can help you with that.

I think you know the litterature regarding the derivation of \frac{\partial{}^2U}{\partial{}\varepsilon_{i}\varepsilon_{j}} for angles and dihedral angles potentials, papers such as Van Workums’ and yours were the one I used to write the compute command.

However, in the case of impropers, two issues appeared:

I didn’t find derivations specifically done for improper angles. Assuming that they were the same as for dihedral is a bit of a leap of faith. I tried to make the derivation myself on paper but didn’t find enough time to finish and unfortunately went on to work with more urgent matter.

I am unsure that the atom ordering is consistent between LAMMPS improper styles. See the documentation for harmonic improper, distance improper, fourier improper or class2 improper. Sometimes J seems to be the atom out of plane, sometimes it is I… I think this issue was also discussed a lot in the mailing list, @akohlmey or @jewettaij might have things to say on that but as far as I understand, this is due to several people making independent dihedral styles implementations each recomputing the atoms order in its own way.

For these reasons, I have indeed set an impflag to set up the stones for future developments but this needs a bit more of brain juice than I expected. This was even before I made the 1st PR and got in touch with @athomps to implement the numerical difference method. That said, this latter method allowed me to check for the correctness of the dihedral implementation and might help to verify that the impropers indeed follow the same derivation, which could tackle the 1st problem.

But the 2nd point might lead to large refactoring of dihedral styles which needs discussion if considered, and will anyway break backward reproducibility of previous simulations.

If you’re interested as I saw your PR passing, I would gladly help. I also began working again on the compute to some extent and wanted to submit some more work.

Thanks for your detailed response @Germain. Indeed, there is no single “consistent” way of defining improper angles. Even if the same equation for dihedral angles holds for improper angles, the variation in the ordering of the atoms in improper styles will make the final result inaccurate.

Before going into a large refactoring to cater for the improper contribution, there are a few spots where the analytical implementation can be improved. One of them is the inclusion of cross terms for angles and dihedrals (expressions similar to CLASS2 potentials). Is this something you have worked/currently working on?

The refactoring and backward compatibility issues could be avoided, if we add a flag to the dihedral and improper styles that indicate the ordering. This could be useful also for other applications where the relative order matters. If “symbolic” constants are used, then people reading the source code can easily see what the intended order for a given style/class is.

Yeah, for the longest time, all implementation details for contributions to LAMMPS were left to the implementer and correct operations were considered the responsibility of the contributor only. This was particularly true for any style that was in a USER-XXX package folder and most true for styles in USER-MISC.

This changed when we started various refactoring efforts and added systematic testing and static code analysis. This exposed a lot of inconsistencies that were subsequently eliminated. Adding the GitHub review process prompted to check for inconsistencies now before anything is merged and has massively reduced the spread in code quality and overall improved consistency and reliability. This is why we some time ago did away with the USER- designation of packages. However, some of the inherited implementation inconsistencies are still present and not easy to remove because the code has been like this for a long time.

Looking forward, we are trying to keep things moving and are quite willing to make (source code) incompatible changes for as long as they remain internal and don’t change behavior. For changes in behavior, there has to be a convincing reason or it has to be applied to either a new feature or a feature that is very unusual or exotic and thus the impact is very limited.

To give an example of a known (legacy) inconsistency, you can look at bond style fene. In its current form, it is only fully correct for “lj” units and epsilon = sigma = 1.0 (check out the function “computing” the equilibrium distance). But since this style is used regularly and has been like this for a very long time, we have not (yet) changed it.

As I said I got quite busy on other project and had to put this one aside in term fo coding. I’ve not done modifications to the code since submitting it and have not been heavily working on the c++ files.

I suppose you the contributions you have in mind are the E_{x} contributing to the total E of a given dihedral? For now the method for dihedral returns two values: \partial{}_{\chi}E and \partial{}^2_{\chi}E used in the compute. That’s it. I think a separate function could be used to make a specific difference for individual terms. I personally did not see any interesting use case for that or a way to implement it but if you have one in mind I am curious about it. As far as I know, class2 is the only dihedral style implemented in LAMMPS that splits its energy in clearly defined energy terms, the others being modes of Fourier like sums to give periodic energy forms. But if making such an extraction function for class2 values, or any other potential form, is of interest, fine by me.

@Germain I am interested in expanding the “analytical” part for systems that can be described by PCFF (which is a class 2 force field) / OPLS -AA and TraPPE-UA. Apart from the contribution of cross terms (even for angle_style class2 type), one has to address the inclusion of the long-range coulombic interactions to this formulation. Not sure how far I can go.

@akohlmey Could the “born_matrix” method implemented in this compute be useful for the “compute hma” ? HMA requires a “single_hessian” method for the calculation of heat capacity which is pretty much the same as the “born_matrix”. Perhaps I should contact Andrew Schultz and David Kofke directly.

For what it’s worth, the only paper I know which is concerned with analytical derivation for electromagnetic terms is this Phys. Rev. E from Van Workum et al.

But to be fair, its formalism is really hard to follow and I do not know how to implement it in a clean way.
I wonder if the use of numdiff methods specifically for kspace forces with a dedicated virial compute could give a good enough approximation. But that might be at the cost of an increased computation time.

Let me take a look at that P3M stuff. As long as we have the Fourier and inverse kernels (and are happy to launch a million FFTs), it should be doable.

Yes, currently compute hma cv is only implemented for a single pair style lj/smooth/linear. By hooking in the numdiff variant of compute born_matrix, this could probably be extended to all energy models. It would be well worth discussing this with Andrew Schultz (ajs42 at buffalo.edu) and David Kofke (kofke at buffalo.edu) .

Thanks a lot for sharing this paper. I was not aware of it. I had a quick glance and indeed it is not a bed time reading (unless someone wants to get nightmares )