GB should return only one torque

[fweik]

Currently GB returns one force, but two independent torques. This is not needed, as the opposing torque can be calculated (as is the opposing force). It should instead return a ParticleForce.

[jngrad]

How would you deduce the opposing torque? GB torques aren't symmetric.

[jngrad]

Here is a MWE to illustrate the asymmetry of the GB torques, based on an idea from @christophlohrmann :

import numpy as np
import espressomd

system = espressomd.System(box_l=3*[10])
system.cell_system.skin = 0.
system.time_step = 0.01
system.non_bonded_inter[0, 0].gay_berne.set_params(
    k1=1.2, k2=2.4, mu=2, nu=5, sig=1, eps=1, cut=3)

# particles are scattered along the x-axis and oriented parallel to the z-axis
p1 = system.part.add(id=0, rotation=(1, 1, 1), type=0, pos=(0, 0, 0))
p2 = system.part.add(id=1, rotation=(1, 1, 1), type=0, pos=(0.5, 0, 0))
print(p1.director)
print(p2.director)
# orient particle 2 parallel to the y-axis
p2.rotate(axis=(1, 0, 0), angle=np.pi / 2.)
for i in range(13):
    # rotate particle 2 around the z-axis to overlap with particle 1
    system.integrator.run(recalc_forces=True, steps=0)
    t1 = p1.torque_lab
    t2 = p2.torque_lab
    angle = i * 360 / 12
    print(f'{angle:>3.0f} -> p1: [{t1[0]:+.1f},{t1[1]:+.1f},{t1[2]:+.1f}] p2: [{t2[0]:+.1f},{t2[1]:+.1f},{t2[2]:+.1f}]')
    p2.rotate(axis=(0, 0, 1), angle=2 * np.pi / 12)

output:

  0 -> p1: [-0.0,+0.0,+0.0] p2: [+0.0,+0.0,+0.0]
 30 -> p1: [-0.0,-0.0,+0.0] p2: [+0.0,-0.0,-44358.8]
 60 -> p1: [-0.0,-0.0,+0.0] p2: [+0.0,-0.0,-187597.8]
 90 -> p1: [+0.0,-0.0,+0.0] p2: [-0.0,-0.0,+0.0]
120 -> p1: [+0.0,-0.0,+0.0] p2: [+0.0,-0.0,+187597.8]
150 -> p1: [-0.0,+0.0,+0.0] p2: [-0.0,+0.0,+44358.8]
180 -> p1: [-0.0,+0.0,+0.0] p2: [+0.0,+0.0,+0.0]
210 -> p1: [-0.0,-0.0,+0.0] p2: [+0.0,-0.0,-44358.8]
240 -> p1: [-0.0,-0.0,+0.0] p2: [-0.0,-0.0,-187597.8]
270 -> p1: [-0.0,+0.0,+0.0] p2: [+0.0,-0.0,-0.0]
300 -> p1: [+0.0,+0.0,+0.0] p2: [+0.0,-0.0,+187597.8]
330 -> p1: [+0.0,+0.0,+0.0] p2: [-0.0,-0.0,+44358.8]
360 -> p1: [+0.0,-0.0,+0.0] p2: [-0.0,-0.0,+0.0]

The opposing torque cannot be deduced from a linear transformation nor from a vector product.

[RudolfWeeber]

So how is it with angular momentum conservation? dL/dt =[T_1 +T_2] +[ F_1 \times r_1 +F_2 \times R_2] =0 Why is/shouldn’t that be fulfilled? Or is it not valid to combine internal angular momentum and angular momentum from translation in that way? @christophlohrmann T_1,2 are torques on particles 1,2, F_1,2 and r_1,2 are forces and positions. For conservative forces, the second bracket should be convertible to [F_12 \times r_12] with the pair force and distance vector

[jngrad]

So how is it with angular momentum conservation?

You are right. Looking at eqn 6 in Allen et Germano 2006, Molecular Physics, 104:20-21 (doi:10.1080/00268970601075238):

We can indeed recover the second torque from a cross product with the force. Numerical precision is around 1e-12 in python and machine precision in the core. Returning a ParticleForce will help remove a few pointers.