Langevin thermostat with anisotropic friction is broken

[christophlohrmann]

Brownian themorstat is fine, though

import numpy as np
import espressomd

use_brownian = False

system = espressomd.System(box_l=3 * [10])
system.time_step = 0.1
system.cell_system.skin = 0.4
if use_brownian:
    system.thermostat.set_brownian(kT=0.0, gamma=1.0, seed=42)
    system.integrator.set_brownian_dynamics()
else:
    system.thermostat.set_langevin(kT=0.0, gamma=1.0, seed=42)

ext_force = 1.5
gamma_parallel = 0.8
gamma_ortho = 0.2
mass = 0.05
part = system.part.add(pos=3 * [0], 
                       gamma=[gamma_ortho, gamma_ortho, gamma_parallel],
                       mass=mass)

# check particle orientation unchanged
part.director = [0, 0, 1]
part.ext_force = [ext_force, 0, ext_force]
system.integrator.run(1000)
np.testing.assert_allclose(
    part.v, [ext_force / gamma_ortho, 0, ext_force / gamma_parallel]
)

# check after rotation: switch parallel and perpendicular
part.director = [1, 0, 0]
part.ext_force = [ext_force, 0, ext_force]
system.integrator.run(10000)
np.testing.assert_allclose(
    part.v, 
    [ext_force / gamma_parallel, 0, ext_force / gamma_ortho], 
    atol=1e-10
)

# check noise
kT = 0.1234
if use_brownian:
    system.thermostat.set_brownian(kT=kT, gamma=1.0, seed=42)
    system.integrator.set_brownian_dynamics()
else:
    system.thermostat.set_langevin(kT=kT, gamma=1.0, seed=42)
    
part.ext_force = 3 * [0.]
part.rotation = 3*[False]

def check_temperature():
    n_samples = 10000
    steps_per_sample = 3

    vels = []
    for _ in range(n_samples):
        system.integrator.run(steps_per_sample)
        vels.append(np.copy(part.v))
    
    # https://en.wikipedia.org/wiki/Thermal_velocity
    # use 1D-Formula for each direction separately 
    thermal_vel = np.sqrt(np.mean(np.array(vels)**2, axis=0))
    kT_measured = mass * thermal_vel**2

    np.testing.assert_allclose(kT_measured,kT, rtol=3e-2)

part.director = [0,0,1]
check_temperature()
part.director = [0,1,0]
check_temperature()

I only checked friction.
This

espresso/src/core/thermostats/langevin_inline.hpp

Line 84 in 6d08b88

const Utils::Matrix<double, 3, 3> noise_op = boost::qvm::diag_mat(pref_noise);

also looks wrong to me. The friction force is converted from body to lab frame, but the noise is not.

EDIT:
Added temperature check that fails due to not-body-to-lab-converted noise prefactor

[RudolfWeeber]

this ticket is suitable for new/early contributors, as most likely the fix is limited to a very few places in the Langevin thermostat code, pertaining to conversions between lab and the particle's co-rotating body frame.

[fweik]

@christophlohrmann White noise is invariant under rotation, so there is no need to transform it.

[christophlohrmann]

even for anisotropic particles? It just feels weird that the friction matrix is "full" and the noise matrix always a diagonal matrix in lab frame no matter the particle orientation

[fweik]

Well my statement is correct, but not related :-) You are right: it's about the prefactors, not the noise, which of course needs to be transformed. Would be maybe a good idea to put the calculation for the anisotropic case into a function and add a unit test...

[christophlohrmann]

couldn't agree more