rotation - fix -interaction

src/core/thermostats/brownian_inline.hpp

@@ -265,7 +265,7 @@ bd_drag_rot(Thermostat::GammaType const &brownian_gamma_rotation, Particle &p,
 
   Utils::Vector3d dphi = {};
   for (int j = 0; j < 3; j++) {
-    if (!p.is_fixed_along(j)) {
+    if (p.can_rotate_around(j)) {
       // only a conservative part of the torque is used here
 #ifndef PARTICLE_ANISOTROPY
       dphi[j] = p.torque()[j] * dt / gamma;
@@ -304,12 +304,11 @@ bd_drag_vel_rot(Thermostat::GammaType const &brownian_gamma_rotation,
 
   Utils::Vector3d omega = {};
   for (int j = 0; j < 3; j++) {
-    if (!p.is_fixed_along(j)) {
-      // only conservative part of the force is used here
-#ifndef PARTICLE_ANISOTROPY
-      omega[j] = p.torque()[j] / gamma;
-#else
+    if (p.can_rotate_around(j)) {
+#ifdef PARTICLE_ANISOTROPY
       omega[j] = p.torque()[j] / gamma[j];
+#else
+      omega[j] = p.torque()[j] / gamma;
 #endif // PARTICLE_ANISOTROPY
     }
   }
@@ -343,7 +342,7 @@ bd_random_walk_rot(BrownianThermostat const &brownian, Particle const &p,
   auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_ROT_INC>(
       brownian.rng_counter(), brownian.rng_seed(), p.id());
   for (int j = 0; j < 3; j++) {
-    if (!p.is_fixed_along(j)) {
+    if (p.can_rotate_around(j)) {
 #ifndef PARTICLE_ANISOTROPY
       if (sigma_pos > 0.0) {
         dphi[j] = noise[j] * sigma_pos * sqrt(dt);
@@ -378,7 +377,7 @@ bd_random_walk_vel_rot(BrownianThermostat const &brownian, Particle const &p) {
   auto const noise = Random::noise_gaussian<RNGSalt::BROWNIAN_ROT_WALK>(
       brownian.rng_counter(), brownian.rng_seed(), p.id());
   for (int j = 0; j < 3; j++) {
-    if (!p.is_fixed_along(j)) {
+    if (p.can_rotate_around(j)) {
       domega[j] = sigma_vel * noise[j] / sqrt(p.rinertia()[j]);
     }
   }

testsuite/python/brownian_dynamics.py

@@ -152,6 +152,32 @@ def test_07__virtual(self):
         np.testing.assert_almost_equal(np.copy(virtual.v), [0, 0, 0])
         np.testing.assert_almost_equal(np.copy(physical.v), [0, 0, 0])
 
+    @utx.skipIfMissingFeatures(["ROTATION", "EXTERNAL_FORCES"])
+    def test_fix_rotation(self):
+        system = self.system
+        system.time_step = 0.01
+        part = system.part.add(
+            pos=3 * [0.],
+            fix=3 * [True],
+            rotation=3 * [True])
+
+        # torque only
+        part.ext_torque = [0, 0, 1.3]
+        system.thermostat.set_brownian(
+            kT=0, gamma=1, gamma_rotation=1.5, act_on_virtual=False, seed=41)
+        system.integrator.set_brownian_dynamics()
+        system.integrator.run(3)
+        np.testing.assert_allclose(
+            part.omega_lab, [
+                0, 0, 1.3 / 1.5], atol=1e-14)
+
+        # noise only
+        part.ext_torque = 3 * [0.]
+        system.thermostat.set_brownian(
+            kT=1, gamma=1, gamma_rotation=1.5, act_on_virtual=False, seed=41)
+        system.integrator.run(3)
+        self.assertGreater(np.linalg.norm(part.omega_lab), 0.)
+
 
 if __name__ == "__main__":
     ut.main()