Force derivation with weighting function in orientation-dependent potentials

Hello everyone,

I’m working on a new pair potential that depends on particle orientation, and I would appreciate some feedback.

In my model, I’d like to modulate the interaction between particles using a weighting function—something similar to what Noguchi used in this paper (likely used in other publications as well, but this is the one I’m most familiar with).

The potential energy takes the form:

E = \sum_{i<j} U(r) w(r)

where ( w(r) ) is a weighting function that smoothly decays to zero at a cutoff distance ( r_\text{cut} ).

My main question is about how to correctly derive the force from this type of potential. I haven’t found many resources on this, even after searching the forum.

Specifically, should I differentiate the entire potential expression like this:

F_{\text{tot}} = - \nabla_{\vec{r}} \left[U(\hat{r})\, w(r)\right] = - w(r)\frac{\partial U(\hat{r})}{\partial \hat{r}} \frac{\partial \hat{r}}{\partial \vec{r}} - U(\hat{r})\frac{\partial w(r)}{\partial r} \frac{\partial r}{\partial \vec{r}} = w(r) F - U(\hat{r}) w^\prime(r) \hat{r}

or is it acceptable to scale the derived force directly:

F_{\text{tot}} = w(r)\, F

I came across this StackExchange post, which touches on a similar issue, but I’m not entirely convinced by the conclusion, since the full differentiation introduces an extra force term that shouldn’t be ignored.

I understand this may not be strictly LAMMPS-related, but I hope someone has encountered this kind of problem before, or perhaps there’s a LAMMPS pair style with a known implementation that could serve as a reference.

Best regards,
Pietro

This manual page should answer your questions: pair_style lj/gromacs command — LAMMPS documentation

You can also look at the source code: lammps/src/EXTRA-PAIR/pair_lj_gromacs.cpp at 53fec5563cf00bcdda454ccd22c460a39b368228 · lammps/lammps · GitHub

EDIT: Sorry, I ignored the fact that the potential is orientation-dependent.

1 Like