Add force support for `HarmonicTorsion`

[mlund]

This adds force support for HarmonicTorsion and includes some cleanup of BondData and force evaluation.

• Add bend example for MD/MC testing of harmonic bonds and angles
• BondData refactoring

[mlund]

@SHervoe-Hansen is the force calculation consistent with what you derived?

[SHervoe-Hansen]

It passes the tests in terms of RDF is the same for Langevin dynamics and Monte Carlo and does not return NaN.

[mlund]

It passes the tests in terms of RDF is the same for Langevin dynamics and Monte Carlo and does not return NaN.

Thanks, the cross-product code works in Debug mode only, but fails on both GCC and Clang if with RelWithDebMode. I tried updating Eigen from 3.3.7 --> 3.3.9 but to no avail. I kept the code, but enclosed it in a constexpr(false) which will be evaluated at compile-time for further testing.

[mlund]

@SHervoe-Hansen in the unittest,bonds.cpp:312, the distance function is defined as b - a, whereas in Cuboid::vdist() (and all other geometries, see energy.h:471) it's a - b. Is there a lucky cancellation of signs somewhere?

[SHervoe-Hansen]

@SHervoe-Hansen in the unittest,bonds.cpp:312, the distance function is defined as b - a, whereas in Cuboid::vdist() (and all other geometries, see energy.h:471) it's a - b. Is there a lucky cancellation of signs somewhere?

We are partly saved because of my poor memory in terms of vector notation: like we wanted to use the vector ab = b - a however I had named it ba so we might gain our sign cancelation from that.

I will go over the math and code ones more in the near future!

[SHervoe-Hansen]

I have been testing for more odd behaviour, it turns out the debug only code can run in release mode if you insert std::cout << force0; and std::cout << force2 before the calculation of force1, which I can only associate with uninitialised variable behaviour.

[mlund]

I have been testing for more odd behaviour, it turns out the debug only code can run in release mode if you insert std::cout << force0; and std::cout << force2 before the calculation of force1, which I can only associate with uninitialised variable behaviour.

I think we just ran into one of Eigen's common pitfalls :-) Fixed in latest commit. Still a bit worried about the force sign more generally in Faunus (bonds, pair potential).

[SHervoe-Hansen]

Still a bit worried about the force sign more generally in Faunus (bonds, pair potential).

Let me do a worked out example for unit testing and such that we can compare and I will also check the code for the angular potential. It seems the major problem is the correct usage of the calculateDistance function. We could maybe look into documentation for future purposes to avoid these mistakes?

[mlund]

Still a bit worried about the force sign more generally in Faunus (bonds, pair potential).

Let me do a worked out example for unit testing and such that we can compare and I will also check the code for the angular potential. It seems the major problem is the correct usage of the calculateDistance function. We could maybe look into documentation for future purposes to avoid these mistakes?

Yes, good idea with more unittests. I made a commit clarifying both vdist() and force() and renamed arguments to reflect directions. I think documentation for the external CoulombGalore library might need updating as it should state on which particle the resulting force acts upon (ping @bjornstenqvist).

[SHervoe-Hansen]

The commit https://github.com/mlund/faunus/pull/363/commits/4e34050d54e6552a20e700324328e32bfe49e56e now reflects the addition of a unit test checking the direction of the force vectors to be correct. The unit test system is illustrated here.

The angular harmonic potential now passes this test. Agreement between MC and MD however still remains to be tested in addition to the speed of the two implementations.

[mlund]

The angular harmonic potential now passes this test. Agreement between MC and MD however still remains to be tested in addition to the speed of the two implementations.

Now ctest -R bend fails and with the latest "cross product" variant we get the graph below (green), indicating a sign is flipped somewhere. You can switch version by simple changing the bool in if constexpr(); no need to uncomment.