New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Fixes for PairAIREBO::bondorderLJ #376
Conversation
Since Etmp (representing sum_kijl omega_kijl * w_ik * w_jl) is not reset between the forward and reverse pass, the value used by later calculation will be twice the expected values. One could instead reset Etmp between these passes, but there really is no reason to calculate it twice.
The bonderorderLJ function operates on a facticious distance |rij|, i.e. everything gets calculated "as if" atoms i and j were a given distance alpha apart. Mathematically, bondorderLJ is a function of rij (a vector), that is (in terms of the real distance Rij) rij = alpha * Rij/|Rij|. When we calculate the forces in bondorderLJ, we have to make sure to chain in this derivative whenever we calculate derivatives w.r.t. rij. The right correction, as it turns our, is Fij = alpha / |Rij| * (Identity(3,3) - Rij * Rij^T / |Rij|^2) * fij. This commit only fixes this for the p_ij^sigma pi terms, which were modified to separate out the d/drij derivative in the cosine calculation. Now, derivatives are taken w.r.t. the connecting edges instead of the edge points.
The code tries to make this distinction between the real distance (r23) and the facticious one (rij), but does not do so very well. It is better if those two variables have the same value everywhere, and apply the correction where necessary. The current way to use the values is incorrrect. Remove those calculations that effectively are derivatives w.r.t. |rij| (the facticious distance), is constant and thus the chained derivative (d|rij|/dRij) is always zero. Apply the corrections due to drij/dRij in the sum omega term.
Since we compute dvdw as d vdw / d rij, we have to also compute dslw as d slw / d rij. Currently, we compute -1/r d slw/d rij, which leads to incorrect results when the two are later combined. Alternatively, one could also modify dvdw to be -1/r d vdw/d rij, which would be a more standard way to do LJ calculations, but this way seems more consistent.
Hi Markus, Based on the analysis that you emailed to @sjplimp on Feb 21, these are highly likely to be correct corrections to incorrect code. They should be merged ASAP. Aidan |
@athomps in order to signal steve to merge a pull request, we use the following convention: you set the Assignees on the right side to sjplimp.
based on your recommendation i have now assigned this PR to steve to be merged into the next patch. |
As alluded to in #59, I believe that there are a number of issues with
PairAIREBO::bondorderLJ
, that might partially cause the numerical behaviour that we all observe. I included some really unclean code in a ZIP-File attachment in that comment to illustrate the changes I am proposing, and this pull-request essentially contains the same stuff in a more tidy manner.It all boils down to the fact that the handling of the facticious distance (
rij = alpha*Rij/|Rij|
) in there seems incorrect: The derivatives ford/drij
need to be corrected withdrij/dRij
, andd|rij|/dRij=0
.