Gay-Berne has low test coverage

[fweik]

The Gay-Bernen potential is not properly tested and documented.

[Kischy]

@fweik One quick question: How can I run only a specific python test or testfile from bash? I only found make -j$(nproc) check_python in the wiki. And nothing on running a specific test or testfile.

[RudolfWeeber]

From: Christian Haege notifications@github.com Sent: 19 June 2019 20:09 To: espressomd/espresso espresso@noreply.github.com Cc: Subscribed subscribed@noreply.github.com Subject: Re: [espressomd/espresso] Gay-Berne has low test coverage (#2922)

https://github.com/fweik @fweik One quick question: How can I run only a specific python test or testfile from bash? I only found make -j$(nproc) check_python in the wiki. And nothing on running a specific test or testfile.

You can just run them like normal Espresso scripts

pypresso testsuite/python/interactions_nonbonded.py

[Kischy]

Thank you.

[Kischy]

I would like to add some function to the tests for better expressiveness. Where do you put functions like that? They will only be used in my tests.

[KaiSzuttor]

you can have a look into the testscripts in testsuite/python. What do you mean by they are only used in your tests?

[Kischy]

I usally put everything, even one liners, into extra functions, so it is expressive. I guess it is not neccessary here, nervermind.

[fweik]

We use the testing framework from the python standard library, with that you typically have a class for each test module. So if you want to have functions that you only use internally you can just keep them in the test class.

[Kischy]

Okay, thank you.

[RudolfWeeber]

Just to be sure: In case you have not seen that yet, testsuite/python/interactions-nonbonded.py:test_gb() is where the existing test case for the energy lives. Maybe, you can add to that. What is missing are test for the forces and torques. The expression for those need to be derived from the expression for the potential or found in the literature. Unfortunately, there are quite a view variations of the potential aroudn.

[Kischy]

Yes, I found it thank you. @RudolfWeeber I have no idea how to calculate the torques yet, but I'll figure it out.

[Kischy]

I have looked at the problem yesterday evening and I am positively confused now. Maybe someone with more experience in md can help me.

The problem is the following: When calculating the energy in gb.hpp in function 'gb_pair_energy' the cutoff of the potential is enforced like this:

  if (!(dist < ia_params->GB_cut))
    return 0.0;

This if fine, but at the end of the function the returned potential includes a energy contribution from the cutoff:

return E(r_eff(dist)) - E(r_eff(ia_params->GB_cut));

The second part of the return value E(r_eff(ia_params->GB_cut)) does not make complete sense to me anymore. Why isn't just E(r_eff(dist)) returned? This is not mentioned in the documentation (http://espressomd.org/html/doc/inter_non-bonded.html#gay-berne-interaction) nor it is part of the original potential and therefore I find it a bit counter intuitive.

What do you think? Is this kind of cutoff contribution something normal and to be expected?

[RudolfWeeber]

On Sun, Jun 23, 2019 at 03:21:24AM -0700, Christian Haege wrote:

Which other non-bonded interaction is depended on the angle between the interaction centers? The dipolar interaction. E.g., the short range part of the dipolar P3M.

[RudolfWeeber]

On Mon, Jun 24, 2019 at 12:48:28AM -0700, Christian Haege wrote:

I have looked at the problem yesterday evening and I am positively confused now. Maybe someone with more experience in md can help me.

The problem is the following: When calculating the energy in gb.hpp in function 'gb_pair_energy' the cutoff of the potential is enforced like this:

  if (!(dist < ia_params->GB_cut))
    return 0.0;

This if fine, but in the end of the function the returned potential includes a energy contribution from the cutoff:

return E(r_eff(dist)) - E(r_eff(ia_params->GB_cut));

The second part of the return value E(r_eff(ia_params->GB_cut)) does not make complete sense to me anymore. Why isn't just E(r_eff(dist)) returned? This is not mentioned in the documentation (http://espressomd.org/html/doc/inter_non-bonded.html#gay-berne-interaction) nor it is part of the original potential and therefore I find it a bit counter intuitive.

What do you think? Probably, like the normaal Lennard Jones interaction, GB is shifted, such that the energy is zero at the cutoff distance. The pure GB expression does not take that into account, so when returning the energy, the appropriate shift has to be added. The energy shift is exactly the energy at the cutoff distance.

[Kischy]

Thank you. The documentation should be updated to include this.

[RudolfWeeber]

Any news here?

[Kischy]

No, I am on holiday now and before I didn't have enough time to finish it. I hope to have time in August or September.

I wrote tests for the GB force using gradient checking and they fail. So the next step is to find the error in the source code.

[jngrad]

@Kischy please note the source code for GB changed substantially last week, using Vector3d to improve readability. The numerical results however should be unchanged.

[Kischy]

Ok, I'll take a look at the changes when I'm home. Thank you.

[Kischy]

Part 1: Errors in the force calculation

I finally found the errors in file "gay_berne.hpp" in the force calculation. I will do a pull request in the next few day.

Both errors are related to the parameter gay_berne.sig. The error does not matter if one uses a value of 1.0 for gay_berne.sig, in all other cases the error is dire!

The first error is in this line in function inline void add_gb_pair_force:

auto Koef2 = int_pow<3>(Sigma) * 0.5;

should be:

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

The second error is in these lines in function inline void add_gb_pair_force:

auto const dU_dr = E *
                       (Koef1 * Brhi2 * (Brack - BrackCut) -
                        Koef2 * Brhi1 * (Bra12 - Bra12Cut) 
                        - Bra12 * dist) /
                       sqr(dist);

should be:

auto const dU_dr = E *
                       (Koef1 * Brhi2 * (Brack - BrackCut) -
                        Koef2 * Brhi1 * (Bra12 - Bra12Cut) 
                         - Bra12 * dist /  ia_params->gay_berne.sig) /
                       sqr(dist);

These errors were not easy to find and as I said I will do a pull request in the next few days.

---

Part 2: Two other problems

First problem: In the force calculation (function inline void add_gb_pair_force) in file "gay_berne.hpp" there are two cases: The cases are

if (X < 1.25) { /*THE FORCE CALCULATION IS DONE HERE*/}
else {   
    Koef1 = 100;
    force += Koef1 * d;
  }

This distinction is, in my opinion, arbitrary and not expected. In my opinion this distinction should be removed and the force should always be calculated correctly, independent of the distance between the particles. If someone wants to cap the forces, they can do so using the usual way. What do you think?

Second problem: I do not exactly know how to test the torques, because I do not know how to exactly calculate the torques. I thought that the torques would just be the cross product of the force between the particles and the unit vector of a particle, but to be honest, I am not sure. Can someone tell me how the torque is calculated in general? Does anyone have a good reference for this?

Thank you kindly

[jngrad]

Thank you for looking into this.

These errors were not easy to find

The force function is less readable than the energy function due to the use of meaningless variable names for intermediate results (a, b, Plus1, Minus1, etc.). Updating the nomenclature would probably help a lot when comparing the code to the formula in the literature. If you used a paper to work out the math, please cite it in the Doxygen comment block.

This distinction is, in my opinion, arbitrary and not expected

This magic value of 1.25 was introduced 16 years ago. It's most likely a relic from the days where LJ-like potentials were heavily monitored with tracing/debugging macros. Nowadays we have force caps and other mechanisms to avoid large force fluctuations during the minimization step of densely packed liquids. I would say you can remove this conditional and always calculate the force.

Can someone tell me how the torque is calculated in general? Does anyone have a good reference for this?

Not sure if this will help, but Swope and Ferguson 1992 explored two chain rule derivations for an angle-dependent potential between two particle pairs connected by a bond and went through the analytical solution for the torque in great details (equations 21 to 43).

[Kischy]

Thank you for your fast reply.

The force function is less readable than the energy function due to the use of meaningless variable names for intermediate results (a, b, Plus1, Minus1, etc.). Updating the nomenclature would probably help a lot when comparing the code to the formula in the literature. If you used a paper to work out the math, please cite it in the Doxygen comment block.

I did not use a paper to calculate the forces. I just calculated the (negative) derivatives of the potential using Sympy. Then I wrote my equations using the meaningless nomenclature used in code and compared them. Therefore I can't cite anything. Concerning the nomenclature: The equations are quite long and therefore the original author of the code probably substituted a lot of it. I don't know if it can be improved much without complicating it to much. I'll have a look, if I have time.

Not sure if this will help, but Swope and Ferguson 1992 explored two chain rule derivations for an angle-dependent potential between two particle pairs connected by a bond and went through the analytical solution for the torque in great details (equations 21 to 43).

I'll have a look, thank you. To be more specific: My problem is not within the code of espresso. At the moment I assume it is correct. My problem is the testing. I want to setup the test for the torque, before checking the code. When calculating the force I used gradient checking, because it is independent of the actual derivative. For my test case I have the correct force, the vector between the particles and the directors of the particles: Now I want to be able to calculate the torque from that. How can I do that in the most general way? It does not need to be fast!

[Kischy]

I'll have a look, thank you. To be more specific: My problem is not within the code of espresso. At the moment I assume it is correct. My problem is the testing. I want to setup the test for the torque, before checking the code. When calculating the force I used gradient checking, because it is independent of the actual derivative. For my test case I have the correct force, the vector between the particles and the directors of the particles: Now I want to be able to calculate the torque from that. How can I do that in the most general way? It does not need to be fast!

I found a way to test the torque. No help needed at this time. Thank you.

[Kischy]

Quick question: Am I allowed to make a new file or at least a new class for the python integration tests of the GB potential, force and torque? It would be much cleaner, but I will code to your preference.

[jngrad]

We have a few exceptions, such as testsuite/python/thermalized_bond.py, but if possible we would prefer you keep the GB test logic in testsuite/python/interactions_non-bonded.py and testsuite/python/tests_common.py. This is easier to maintain, and you won't have to tinker with our CMake logic.

[Kischy]

Ok, thank you. I will create some helper functions for clarity, if that is ok.

[jngrad]

Sure, we use helper functions all the time. We write them as class methods (without the test_ prefix) if they can be used by other test methods, otherwise we place them inside the test methods. See for example the 4 helper functions defined in ObservableTests::test_BondDihedrals().

[Kischy]

Ty

[Kischy]

I did a pull request (#3091).

I did not change the nomenclature, because I could not think of a better one yet. The tests should be sufficient now and the bugs are fixed. Will the GB potential stay in espresso now that it is properly tested? Do I need to do something else?

If there are problems with the pull request, please tell me on #3091.