Tests and Bugfixes for the GB potential #3091

Kischy · 2019-08-21T14:38:24Z

Fixes #2922

Description of changes:

Added tests for the gb potential for the force and torque
Fixed bugs in gay_berne.hpp
Removed condition in gay_berne.hpp

PR Checklist

Tests?
- Interface
- Core
Docs?

…emoved condition.

codecov · 2019-08-21T15:15:34Z

Codecov Report

Merging #3091 into python will increase coverage by <1%.
The diff coverage is 100%.

@@           Coverage Diff           @@
##           python   #3091    +/-   ##
=======================================
+ Coverage      83%     83%   +<1%     
=======================================
  Files         530     530            
  Lines       26100   26102     +2     
=======================================
+ Hits        21915   21919     +4     
+ Misses       4185    4183     -2

Impacted Files	Coverage Δ
src/core/nonbonded_interactions/gay_berne.hpp	`98% <100%> (+2%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update d231cd3...975bc21. Read the comment docs.

jngrad · 2019-08-21T18:37:24Z

So far we've always used analytical solutions to calculate reference forces. Using a numerical approximation of the gradient is an interesting approach. I'm actually surprised that forces calculated with an integral step of 2e-7 can pass a self.assertAlmostEqual test at precision 1e-7. I've started the tests on non-amd64 architectures (pipeline 8879) to see if there is a risk of numerical instability.

After re-organizing my version of the force function to match the nomenclature of the variables in the energy function, I still can't follow the force code. For example, in the calculation of epsilon, at line:

espresso/src/core/nonbonded_interactions/gay_berne.hpp

Lines 70 to 72 in d79e9c1

    
           auto const E = 4 * ia_params->gay_berne.eps * 
        
                          pow(E1, ia_params->gay_berne.nu) * 
        
                          pow(E2, ia_params->gay_berne.mu);

I would expect variable E2 to be raised to the power of mu - 1 in a subsequent variable to evaluate the derivative (d epsilon / d r), but I can't find it. It would be very helpful if we had the formula of the GB force, possibly with the intermediate steps to obtain it from the derivative of energy. We could write them with LaTeX in a Doxygen page, like we did for the derivation of the bond angle force.

Kischy · 2019-08-22T08:18:36Z

I've started the tests on non-amd64 architectures (pipeline 8879) to see if there is a risk of numerical instability.

I saw that two of them failed. This one failed because the reaction_ensemble test had a timeout. The second one failed because of different timeouts and fails, non of which was the interactions_non-bonded.py. Do you know why those failed?

So far we've always used analytical solutions to calculate reference forces. Using a numerical approximation of the gradient is an interesting approach. I'm actually surprised that forces calculated with an integral step of 2e-7 can pass a self.assertAlmostEqual test at precision 1e-7.

I think this way testing derivatives can be very powerful, if the analytical solution is long or complicated. Although it is much slower to calculate the derivative this way, it is not so prone to error. To advertise the approach a little bit more: If I had just implemented the same analytical solution used in the c++ code for testing I would not have found the error. I borrowed the idea from machine learning, see here and here.

I would expect variable E2 to be raised to the power of mu - 1 in a subsequent variable to evaluate the derivative (d epsilon / d r), but I can't find it.

It is a little bit confusing. The change in power comes with this line:

espresso/src/core/nonbonded_interactions/gay_berne.hpp

Line 77 in d79e9c1

auto Koef1 = ia_params->gay_berne.mu / E2;

And with this line:

espresso/src/core/nonbonded_interactions/gay_berne.hpp

Line 78 in d79e9c1

auto Koef2 = int_pow<3>(Sigma) * 0.5 / int_pow<3>(ia_params->gay_berne.sig);

And also with these lines:

espresso/src/core/nonbonded_interactions/gay_berne.hpp

Lines 89 to 90 in d79e9c1

    
           auto const Bra12 = 6 * X6 * X * (2 * X6 - 1); 
        
           auto const Bra12Cut = 6 * Xcut6 * Xcut * (2 * Xcut6 - 1);

It would be very helpful if we had the formula of the GB force, possibly with the intermediate steps to obtain it from the derivative of energy.

If you want I can send you my jupyter notebook (using python 3) or the python script, where I did the calculations of the derivative.

We could write them with LaTeX in a Doxygen page, like we did for the derivation of the bond angle force.

The problem is that I do not know what paper the orginal author used, so I can't cite anything. I will try to work out a new nomenclatur for the problem and then we can put it in latex code, if you want.

* add extra intermediate variables * use same variable names for identical or very similar intermediate quantities in the force/energy kernels * simplify code logic with vector/scalar utility functions

jngrad · 2019-08-22T10:11:13Z

The non-amd tests are known to sometimes false trigger, that's why we don't run them automatically. The gradient approach seems to work, though. Your UFLDL link states that an epsilon of 1e-4 will usually provide results at 4 digits precision, but in my experience a test is safer when the reference force is at least one order of magnitude more precise than the test threshold delta. Maybe @KaiSzuttor can comment on that.

The book or article used by the original author is unknown, all we know is this code is 16 years old. Given the notation in the energy function and the LaTeX formula in the docs, I would bet on Cleaver et al. 1996, equations 3-4 and 7-9.

Here is my modified version of the GB force function. Although it looks more verbose, when compiled with -O3 it generates 23 less lines of assembly code in gcc and 13 less in Clang (9% resp. 6% less). We could probably reduce the assembly code further with more refactoring based on the analytical expression of the force. I've looked at several sources on the GB potential, but none of them give the expression of its derivative in Cartesian space. If you already worked out the derivation, you can share it here by attaching a file.

Kischy · 2019-08-22T10:29:20Z

The non-amd tests are known to sometimes false trigger, that's why we don't run them automatically. The gradient approach seems to work, though. Your UFLDL link states that an epsilon of 1e-4 will usually provide results at 4 digits precision, but in my experience a test is safer when the reference force is at least one order of magnitude more precise than the test threshold delta. Maybe @KaiSzuttor can comment on that.

Okay, I will change the epsilon of my force calculation and update the pull request.

Here is my modified version of the GB force function. Although it looks more verbose, when compiled with -O3 it generates 23 less lines of assembly code in gcc and 13 less in Clang (9% resp. 6% less). We could probably reduce the assembly code further with more refactoring based on the analytical expression of the force. I've looked at several sources on the GB potential, but none of them give the expression of its derivative in Cartesian space. If you already worked out the derivation, you can share it here by attaching a file.

I worked on a better nomenclature just now too. If it is okay with you I would combine my version with yours and we see how it fits. What do you think?

Kischy · 2019-08-22T10:30:17Z

I can't attach the files here, because the file type is not supported by github. I will send them to you by mail.

Kischy · 2019-08-22T11:04:25Z

Here is my modified version of the GB force function.

I think you're version is the best we can do. Can you just add the changed version to my pull request or should I do so?

Kischy · 2019-08-22T11:07:48Z

Your UFLDL link states that an epsilon of 1e-4 will usually provide results at 4 digits precision, but in my experience a test is safer when the reference force is at least one order of magnitude more precise than the test threshold delta.

I added a commit that changes the delta value to 1x10E-6. For some reason it does not show here yet.

jngrad · 2019-08-22T11:13:36Z

You can use my version. The simplest way to do it is to merge my branch into your branch (this is a trivial fast-forward), then add the new commit that changes the precision. Doing it the other way around (adding the precision commit, then merging my branch) would require a complex merge commit.

Also, I would not change delta to a bigger value, instead I would change epsilon to 1x10E-8.

jngrad · 2019-08-22T11:20:48Z

You probably forgot to push the new commit to the GitHub repository. If it's ok with you, I can add my new version into your branch on GitHub. Then on you computer, you'll have to discard the precision commit (for example with git reset --soft HEAD^ to rewind the branch history to d79e9c1), then pull your GitHub branch into your local repository.

Kischy · 2019-08-22T11:22:11Z

No, it is already pushed. I'll merge your commit, change the delta and epsilon and push again.

Kischy · 2019-08-22T11:56:58Z

Also, I would not change delta to a bigger value, instead I would change epsilon to 1x10E-8.

The value of 1x10E-8 does not work, probably because of the numerical errors mentioned in the UFLDL link . I think 1x10E-7 is the lowest we can go here. Therefore I would recommend changing the delta value to 1x10E-6. Is that fine with you or not?

KaiSzuttor · 2019-08-22T12:18:48Z

testsuite/python/tests_common.py

@@ -596,6 +596,20 @@ def gay_berne_potential(r_ij, u_i, u_j, epsilon_0, sigma_0, mu, nu, k_1, k_2):
    return 4. * epsilon * (rr**-12 - rr**-6)


+def calc_derivative(func, x, eps=1.0e-7, index=None):


could you please use scipy.misc.derivative. If that does not suit your needs you should definitely add a docstring to this function. Actually we could even use sympy. But we probably don't want to depend on that.

I will add a docstring. scipy.misc.derivative does not allow for a multidimensional input.

https://docs.scipy.org/doc/numpy/reference/generated/numpy.vectorize.html

Does that help? I need to manipulate the element at a specific index while the rest stays the same.

I added a docstring. I could have thrown an error if someone gives a multidimensional input and index=None, but I think it would be a "overkill". In the docstring it is stated, that the index needs to be specified if x is an array.

jngrad · 2019-08-22T12:19:26Z

I've just discussed it with Kai, and we're both fine with a larger delta.

Kischy · 2019-08-22T12:21:01Z

I've just discussed it with Kai, and we're both fine with a larger delta.

Thanks for the update.

Added docstring to calc_derivative changed delta value in test_gb

testsuite/python/interactions_non-bonded.py

jngrad · 2019-08-23T08:21:44Z

bors r+

3091: Tests and Bugfixes for the GB potential r=jngrad a=Kischy Fixes #2922 Description of changes: - Added tests for the gb potential for the force and torque - Fixed bugs in gay_berne.hpp - Removed condition in gay_berne.hpp PR Checklist ------------ - [x] Tests? - [ ] Interface - [x] Core - [ ] Docs? Co-authored-by: Christian Kischy Haege <[email protected]> Co-authored-by: Jean-Noël Grad <[email protected]> Co-authored-by: Christian Haege <[email protected]>

…/espresso into Bugfix_Tests_Gay_Berne_Issue_2922

bors · 2019-08-23T08:35:24Z

Canceled

KaiSzuttor · 2019-08-23T08:39:51Z

I was not yet done...

KaiSzuttor · 2019-08-23T08:40:26Z

I think we should have all reviewers approve the PR before merging

Kischy · 2019-08-23T08:42:33Z

You changed the interface of the gradient function so that x0: (3,)? Does that matter that much now that it is specific to the gb_test?

Could you also update the Docstring?
Point in N-dimensional space where the derivatives are calculated
to
Point in 3-dimensional space where the derivatives are calculated

KaiSzuttor · 2019-08-23T08:49:35Z

the scope of the function is limited to test_gb. Have you ever calculated a gradient in more than three dimensions?

Kischy · 2019-08-23T08:51:31Z

the scope of the function is limited to test_gb.

That's why I thought that it doesn't matter.

Have you ever calculated a gradient in more than three dimensions?

No

KaiSzuttor · 2019-08-23T08:52:24Z

I thought it's just more readable

Kischy · 2019-08-23T08:52:41Z

Yeah, you are right.

Kischy · 2019-08-23T08:53:59Z

Can you exchange every "derivative" with "gradient" in the Docstring? I think than we are done. :)

KaiSzuttor · 2019-08-23T09:12:43Z

sorry for the long discussion, but I think we overall improved readability. Thanks for your contribution!

Kischy · 2019-08-23T09:15:22Z

@KaiSzuttor I like your implementation of the gradiant function, it is quite compact. However it might be slower, which does not matter in this case. Just keep that in mind, if you ever use it heavily somewhere else.

Kischy · 2019-08-23T09:21:23Z

sorry for the long discussion, but I think we overall improved readability. Thanks for your contribution!

Yes, it looks better now. Thank you. I'm glad I could help. And no worries about the discussion.

KaiSzuttor · 2019-08-23T09:38:58Z

However it might be slower

Why?

Kischy · 2019-08-23T09:46:01Z

However it might be slower

Why?

There will be one lambda for each call to the function. That will make it slower in comparison. I shouldn't have brought it up, because it does not matter.

KaiSzuttor · 2019-08-23T10:15:06Z

bors r+

3091: Tests and Bugfixes for the GB potential r=KaiSzuttor a=Kischy Fixes #2922 Description of changes: - Added tests for the gb potential for the force and torque - Fixed bugs in gay_berne.hpp - Removed condition in gay_berne.hpp PR Checklist ------------ - [x] Tests? - [ ] Interface - [x] Core - [ ] Docs? Co-authored-by: Christian Kischy Haege <[email protected]> Co-authored-by: Jean-Noël Grad <[email protected]> Co-authored-by: Christian Haege <[email protected]> Co-authored-by: Kai Szuttor <[email protected]>

bors · 2019-08-23T11:13:24Z

Build succeeded

ICP GitLab CI

Kischy added 4 commits August 21, 2019 15:49

Issue #2922: Refactoring energy tests gb potential

efe22bd

Issue #2922: Added tests for gb forces and corrected gay_berne.hpp. R…

d5bc3f1

…emoved condition.

Issue #2922: Added gb torque tests. Some refactoring

100362e

Issue #2922: Some refactoring

7631e45

Kischy changed the title ~~Tests and Bugfixes for the GB potential #2922~~ Tests and Bugfixes for the GB potential Aug 21, 2019

Kischy mentioned this pull request Aug 21, 2019

Gay-Berne has low test coverage #2922

Closed

jngrad self-requested a review August 21, 2019 14:48

Kischy added 2 commits August 21, 2019 16:56

Issue #2922: Updated style to meet code formatting rules

7986e5b

#2922: Changed misleading variable name

d79e9c1

Formatting

4fd0c68

Refactor Gay-Berne force/energy kernels

e421226

* add extra intermediate variables * use same variable names for identical or very similar intermediate quantities in the force/energy kernels * simplify code logic with vector/scalar utility functions

#2922: Changed delta value

d0f9a5c

Kischy added 2 commits August 22, 2019 13:47

Resolved small merge conflict

4ee19cd

Merge fixes of gay berne potential

20e904c

KaiSzuttor reviewed Aug 22, 2019

View reviewed changes

#2922:

26889e0

Added docstring to calc_derivative changed delta value in test_gb

jngrad approved these changes Aug 22, 2019

View reviewed changes

KaiSzuttor reviewed Aug 22, 2019

View reviewed changes

testsuite/python/interactions_non-bonded.py Outdated Show resolved Hide resolved

Rename gradient function

377b216

KaiSzuttor added 2 commits August 23, 2019 10:31

testsuite: improve readability of gradient calculation.

9cb0b7c

Merge branch 'Bugfix_Tests_Gay_Berne_Issue_2922' of github.com:Kischy…

27b283d

…/espresso into Bugfix_Tests_Gay_Berne_Issue_2922

numpy as np.

9a0afdd

jngrad assigned KaiSzuttor Aug 23, 2019

testsuite: docstring.

975bc21

bors bot merged commit 975bc21 into espressomd:python Aug 23, 2019

Kischy deleted the Bugfix_Tests_Gay_Berne_Issue_2922 branch August 23, 2019 11:15

Kischy mentioned this pull request May 20, 2020

Segmentation fault with Espresso 4.1.1 using p3m, Gay-Berne potential, Lennard-Jones potential and virtual sites #3713

Closed

		@@ -596,6 +596,20 @@ def gay_berne_potential(r_ij, u_i, u_j, epsilon_0, sigma_0, mu, nu, k_1, k_2):
		return 4. * epsilon * (rr-12 - rr-6)


		def calc_derivative(func, x, eps=1.0e-7, index=None):

Tests and Bugfixes for the GB potential #3091

Tests and Bugfixes for the GB potential #3091

Conversation

Kischy commented Aug 21, 2019

PR Checklist

codecov bot commented Aug 21, 2019 • edited Loading

Codecov Report

jngrad commented Aug 21, 2019

Kischy commented Aug 22, 2019

jngrad commented Aug 22, 2019

Kischy commented Aug 22, 2019

Kischy commented Aug 22, 2019

Kischy commented Aug 22, 2019

Kischy commented Aug 22, 2019 • edited Loading

jngrad commented Aug 22, 2019

jngrad commented Aug 22, 2019

Kischy commented Aug 22, 2019

Kischy commented Aug 22, 2019

KaiSzuttor Aug 22, 2019

Choose a reason for hiding this comment

Kischy Aug 22, 2019

Choose a reason for hiding this comment

KaiSzuttor Aug 22, 2019

Choose a reason for hiding this comment

Kischy Aug 22, 2019

Choose a reason for hiding this comment

Kischy Aug 22, 2019

Choose a reason for hiding this comment

jngrad commented Aug 22, 2019

Kischy commented Aug 22, 2019

jngrad commented Aug 23, 2019

bors bot commented Aug 23, 2019

Canceled

KaiSzuttor commented Aug 23, 2019

KaiSzuttor commented Aug 23, 2019

Kischy commented Aug 23, 2019 • edited Loading

KaiSzuttor commented Aug 23, 2019

Kischy commented Aug 23, 2019

KaiSzuttor commented Aug 23, 2019

Kischy commented Aug 23, 2019

Kischy commented Aug 23, 2019

KaiSzuttor commented Aug 23, 2019

Kischy commented Aug 23, 2019

Kischy commented Aug 23, 2019 • edited Loading

KaiSzuttor commented Aug 23, 2019

Kischy commented Aug 23, 2019 • edited Loading

KaiSzuttor commented Aug 23, 2019

bors bot commented Aug 23, 2019

Build succeeded

codecov bot commented Aug 21, 2019 •

edited

Loading

Kischy commented Aug 22, 2019 •

edited

Loading

Kischy commented Aug 23, 2019 •

edited

Loading

Kischy commented Aug 23, 2019 •

edited

Loading

Kischy commented Aug 23, 2019 •

edited

Loading