Add force support for `HarmonicTorsion` #363

mlund · 2021-05-07T06:23:25Z

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 · 2021-05-08T05:54:09Z

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

SHervoe-Hansen

It is overall consistent yes, however, I am fearing there is a sign error that cancels out since all tests pass.

SHervoe-Hansen · 2021-05-08T09:59:04Z

src/bonds.cpp

+        if (cosine_angle > 1.0) { // is this needed?
+            cosine_angle = 1.0;
+        } else if (cosine_angle < -1.0) {
+            cosine_angle = -1.0;


I don't understand the need of these lines.

SHervoe-Hansen · 2021-05-08T10:14:10Z

src/bonds.cpp

+            sine_angle = std::numeric_limits<double>::epsilon();
+        }
+        const auto angle = std::acos(cosine_angle);
+        const auto du_dangle = 2.0 * half_force_constant * (angle - equilibrium_angle);


minus sign missing.
F = - dU/dr

I think it's just the direction of the vector (01 vs 10). I checked flipping the sign, and see large differences compared to MC.

SHervoe-Hansen · 2021-05-08T10:18:49Z

src/bonds.cpp

+        }
+        const auto angle = std::acos(cosine_angle);
+        const auto du_dangle = 2.0 * half_force_constant * (angle - equilibrium_angle);
+        auto force0 = du_dangle / sine_angle * d01_inv * (r21 * d21_inv - r01 * d01_inv * cosine_angle);


In more Eigen friendly manner we could maybe write it as:

force0 = du_dangle / r01 * norm(10x(10x12))

Note the cross product of vector 10 and 12 is used in both expressions so maybe calculate that first and take the norm of the cross product between the newly calculated vector and the vector required.

This might be faster as it uses one less trigonometric function (sqrt) in exchange of multiplications. However, your suggestion gives NaN and ctest -R unittests -V fails. Would you be willing to double-check with the current code?

SHervoe-Hansen · 2021-05-08T10:20:58Z

src/bonds.cpp

+        const auto angle = std::acos(cosine_angle);
+        const auto du_dangle = 2.0 * half_force_constant * (angle - equilibrium_angle);
+        auto force0 = du_dangle / sine_angle * d01_inv * (r21 * d21_inv - r01 * d01_inv * cosine_angle);
+        auto force2 = du_dangle / sine_angle * d21_inv * (r01 * d01_inv - r21 * d21_inv * cosine_angle);


force2 = du_dangle / r12 * norm(21x(10x12))

SHervoe-Hansen · 2021-05-10T12:59:45Z

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

mlund · 2021-05-10T15:08:32Z

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 · 2021-05-11T07:04:04Z

@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 · 2021-05-11T10:12:49Z

@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 · 2021-05-11T10:22:32Z

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 · 2021-05-11T17:01:46Z

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 · 2021-05-11T17:15:28Z

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?

Direction of force vector was unclear this commit brings renaming to explicitly give directions of input arguments and calculated forces. Unit test for NewCoulomb is also added. It seems OK, but the documentation in the external CoulombGalore library may need updating

mlund · 2021-05-11T19:44:04Z

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).

…d-force

SHervoe-Hansen · 2021-05-12T15:58:54Z

The commit 4e34050 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 · 2021-05-12T21:11:18Z

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.

SHervoe-Hansen

The code and physics now match and passes all tests. Ready for merge.

mlund added 13 commits May 4, 2021 08:57

Add force for HarmonicBond

9a5e25b

Update

5b2a757

Update

27b3426

Update

e520575

Update

43ac507

Merge branch 'master' into bond-force

d3ce821

Add bend example for combined MC/MD propagation

25963bc

Update install

fd1570e

Update to remove drift

26119a7

Update

66ee279

Update

87813da

Revert Eigen upgrade

39b23d7

Update

9ba8cf4

mlund added the enhancement 💪 label May 7, 2021

mlund added this to the Version 2.5.0 milestone May 7, 2021

mlund added bug 🐛 Broken or wrong functionality code quality ⚙️ Code refactoring / restructuring labels May 7, 2021

This was linked to issues May 7, 2021

Langevin dynamics implementation #282

Closed

Error running example langevin_lj.yml #305

Closed

mlund removed the bug 🐛 Broken or wrong functionality label May 7, 2021

mlund requested a review from SHervoe-Hansen May 8, 2021 05:52

SHervoe-Hansen reviewed May 8, 2021

View reviewed changes

mlund and others added 2 commits May 9, 2021 20:07

Update

8a45768

Reorganized forces for harmonic torsion

62ce105

Conditional run either method

564c449

mlund added 2 commits May 11, 2021 18:54

Fix second force calc. method (stefan)

033d580

Update

4dbed85

Update

2a6635a

SHervoe-Hansen added 2 commits May 12, 2021 17:44

Update to angular harmonic force and unit test

4e34050

Merge branch 'bond-force' of https://github.com/mlund/faunus into bon…

bb248a3

…d-force

SHervoe-Hansen and others added 2 commits May 12, 2021 18:13

Approximate evaluation of unit tests for dipoles

469baee

Minor update to respect if constexpr

f3e5767

mlund added 3 commits May 13, 2021 10:50

Fix unittest force sign

fac13eb

Add average distance to bend test

18a02e1

Clarify sign

2b3d885

SHervoe-Hansen approved these changes May 13, 2021

View reviewed changes

update manual

65ef4a5

mlund merged commit 8ada37f into master May 13, 2021

mlund deleted the bond-force branch May 13, 2021 10:12

mlund linked an issue May 17, 2021 that may be closed by this pull request

Force calculation for bonded interactions #283

Closed

3 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add force support for `HarmonicTorsion` #363

Add force support for `HarmonicTorsion` #363

mlund commented May 7, 2021 •

edited

Loading

mlund commented May 8, 2021

SHervoe-Hansen left a comment •

edited

Loading

SHervoe-Hansen May 8, 2021

SHervoe-Hansen May 8, 2021

mlund May 9, 2021

SHervoe-Hansen May 8, 2021 •

edited

Loading

SHervoe-Hansen May 8, 2021 •

edited

Loading

mlund May 9, 2021

SHervoe-Hansen May 8, 2021 •

edited

Loading

SHervoe-Hansen commented May 10, 2021

mlund commented May 10, 2021 •

edited

Loading

mlund commented May 11, 2021 •

edited

Loading

SHervoe-Hansen commented May 11, 2021 •

edited

Loading

SHervoe-Hansen commented May 11, 2021

mlund commented May 11, 2021

SHervoe-Hansen commented May 11, 2021

mlund commented May 11, 2021

SHervoe-Hansen commented May 12, 2021 •

edited

Loading

mlund commented May 12, 2021 •

edited

Loading

SHervoe-Hansen left a comment

Add force support for HarmonicTorsion #363

Add force support for HarmonicTorsion #363

Conversation

mlund commented May 7, 2021 • edited Loading

mlund commented May 8, 2021

SHervoe-Hansen left a comment • edited Loading

Choose a reason for hiding this comment

SHervoe-Hansen May 8, 2021

Choose a reason for hiding this comment

SHervoe-Hansen May 8, 2021

Choose a reason for hiding this comment

mlund May 9, 2021

Choose a reason for hiding this comment

SHervoe-Hansen May 8, 2021 • edited Loading

Choose a reason for hiding this comment

SHervoe-Hansen May 8, 2021 • edited Loading

Choose a reason for hiding this comment

mlund May 9, 2021

Choose a reason for hiding this comment

SHervoe-Hansen May 8, 2021 • edited Loading

Choose a reason for hiding this comment

SHervoe-Hansen commented May 10, 2021

mlund commented May 10, 2021 • edited Loading

mlund commented May 11, 2021 • edited Loading

SHervoe-Hansen commented May 11, 2021 • edited Loading

SHervoe-Hansen commented May 11, 2021

mlund commented May 11, 2021

SHervoe-Hansen commented May 11, 2021

mlund commented May 11, 2021

SHervoe-Hansen commented May 12, 2021 • edited Loading

mlund commented May 12, 2021 • edited Loading

SHervoe-Hansen left a comment

Choose a reason for hiding this comment

Add force support for `HarmonicTorsion` #363

Add force support for `HarmonicTorsion` #363

mlund commented May 7, 2021 •

edited

Loading

SHervoe-Hansen left a comment •

edited

Loading

SHervoe-Hansen May 8, 2021 •

edited

Loading

SHervoe-Hansen May 8, 2021 •

edited

Loading

SHervoe-Hansen May 8, 2021 •

edited

Loading

mlund commented May 10, 2021 •

edited

Loading

mlund commented May 11, 2021 •

edited

Loading

SHervoe-Hansen commented May 11, 2021 •

edited

Loading

SHervoe-Hansen commented May 12, 2021 •

edited

Loading

mlund commented May 12, 2021 •

edited

Loading