From ae68becf4ad50f194814f9500ecaceda1e1d6fea Mon Sep 17 00:00:00 2001
From: Trung Nguyen <ndactrung@gmail.com>
Date: Tue, 15 Sep 2020 15:27:24 -0500
Subject: [PATCH] Fixed a bug in computing the langevin torques applied to
 rigid bodies

---
 src/RIGID/fix_rigid.cpp       | 29 +++++++++++++++++++++++++----
 src/RIGID/fix_rigid_small.cpp | 33 +++++++++++++++++++++++++++++----
 2 files changed, 54 insertions(+), 8 deletions(-)

diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp
index f74cce8233..aa105d5344 100644
--- a/src/RIGID/fix_rigid.cpp
+++ b/src/RIGID/fix_rigid.cpp
@@ -970,7 +970,7 @@ void FixRigid::apply_langevin_thermostat()
 {
   if (me == 0) {
     double gamma1,gamma2;
-
+    double wbody[3],tbody[3];
     double delta = update->ntimestep - update->beginstep;
     if (delta != 0.0) delta /= update->endstep - update->beginstep;
     t_target = t_start + delta * (t_stop-t_start);
@@ -991,12 +991,33 @@ void FixRigid::apply_langevin_thermostat()
 
       gamma1 = -1.0 / t_period / ftm2v;
       gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
-      langextra[i][3] = inertia[i][0]*gamma1*omega[i][0] +
+
+      // convert omega from space frame to body frame
+
+      wbody[0] = omega[i][0]*ex_space[i][0] + omega[i][1]*ex_space[i][1] +
+        omega[i][2]*ex_space[i][2];
+      wbody[1] = omega[i][0]*ey_space[i][0] + omega[i][1]*ey_space[i][1] +
+        omega[i][2]*ey_space[i][2];
+      wbody[2] = omega[i][0]*ez_space[i][0] + omega[i][1]*ez_space[i][1] +
+        omega[i][2]*ez_space[i][2];
+
+      // compute langevin torques in the body frame
+
+      tbody[0] = inertia[i][0]*gamma1*wbody[0] +
         sqrt(inertia[i][0])*gamma2*(random->uniform()-0.5);
-      langextra[i][4] = inertia[i][1]*gamma1*omega[i][1] +
+      tbody[1] = inertia[i][1]*gamma1*wbody[1] +
         sqrt(inertia[i][1])*gamma2*(random->uniform()-0.5);
-      langextra[i][5] = inertia[i][2]*gamma1*omega[i][2] +
+      tbody[2] = inertia[i][2]*gamma1*wbody[2] +
         sqrt(inertia[i][2])*gamma2*(random->uniform()-0.5);
+
+      // convert langevin torques from body frame back to space frame
+
+      langextra[i][3] = tbody[0]*ex_space[i][0] + tbody[1]*ey_space[i][0] +
+        tbody[2]*ez_space[i][0];
+      langextra[i][4] = tbody[0]*ex_space[i][1] + tbody[1]*ey_space[i][1] +
+        tbody[2]*ez_space[i][1];
+      langextra[i][5] = tbody[0]*ex_space[i][2] + tbody[1]*ey_space[i][2] +
+        tbody[2]*ez_space[i][2];
     }
   }
 
diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp
index 85cadc826e..93e9c2050d 100644
--- a/src/RIGID/fix_rigid_small.cpp
+++ b/src/RIGID/fix_rigid_small.cpp
@@ -818,6 +818,7 @@ void FixRigidSmall::initial_integrate(int vflag)
 void FixRigidSmall::apply_langevin_thermostat()
 {
   double gamma1,gamma2;
+  double wbody[3],tbody[3];
 
   // grow langextra if needed
 
@@ -837,12 +838,15 @@ void FixRigidSmall::apply_langevin_thermostat()
   double mvv2e = force->mvv2e;
   double ftm2v = force->ftm2v;
 
-  double *vcm,*omega,*inertia;
+  double *vcm,*omega,*inertia,*ex_space,*ey_space,*ez_space;
 
   for (int ibody = 0; ibody < nlocal_body; ibody++) {
     vcm = body[ibody].vcm;
     omega = body[ibody].omega;
     inertia = body[ibody].inertia;
+    ex_space = body[ibody].ex_space;
+    ey_space = body[ibody].ey_space;
+    ez_space = body[ibody].ez_space;
 
     gamma1 = -body[ibody].mass / t_period / ftm2v;
     gamma2 = sqrt(body[ibody].mass) * tsqrt *
@@ -853,12 +857,33 @@ void FixRigidSmall::apply_langevin_thermostat()
 
     gamma1 = -1.0 / t_period / ftm2v;
     gamma2 = tsqrt * sqrt(24.0*boltz/t_period/dt/mvv2e) / ftm2v;
-    langextra[ibody][3] = inertia[0]*gamma1*omega[0] +
+
+    // convert omega from space frame to body frame
+
+    wbody[0] = omega[0]*ex_space[0] + omega[1]*ex_space[1] +
+      omega[2]*ex_space[2];
+    wbody[1] = omega[0]*ey_space[0] + omega[1]*ey_space[1] +
+      omega[2]*ey_space[2];
+    wbody[2] = omega[0]*ez_space[0] + omega[1]*ez_space[1] +
+      omega[2]*ez_space[2];
+
+    // compute langevin torques in the body frame
+
+    tbody[0] = inertia[0]*gamma1*wbody[0] +
       sqrt(inertia[0])*gamma2*(random->uniform()-0.5);
-    langextra[ibody][4] = inertia[1]*gamma1*omega[1] +
+    tbody[1] = inertia[1]*gamma1*wbody[1] +
       sqrt(inertia[1])*gamma2*(random->uniform()-0.5);
-    langextra[ibody][5] = inertia[2]*gamma1*omega[2] +
+    tbody[2] = inertia[2]*gamma1*wbody[2] +
       sqrt(inertia[2])*gamma2*(random->uniform()-0.5);
+
+    // convert langevin torques from body frame back to space frame
+
+    langextra[ibody][3] = tbody[0]*ex_space[0] + tbody[1]*ey_space[0] +
+      tbody[2]*ez_space[0];
+    langextra[ibody][4] = tbody[0]*ex_space[1] + tbody[1]*ey_space[1] +
+      tbody[2]*ez_space[1];
+    langextra[ibody][5] = tbody[0]*ex_space[2] + tbody[1]*ey_space[2] +
+      tbody[2]*ez_space[2];
   }
 }