diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index dd40b0fe45..64c7a16ffd 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,23 @@ 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 + + MathExtra::transpose_matvec(ex_space[i],ey_space[i],ez_space[i],omega[i],wbody); + + // 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 + + MathExtra::matvec(ex_space[i],ey_space[i],ez_space[i],tbody,&langextra[i][3]); } } diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 69f5b2b004..b6d68a5222 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,23 @@ 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 + + MathExtra::transpose_matvec(ex_space,ey_space,ez_space,omega,wbody); + + // 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 + + MathExtra::matvec(ex_space,ey_space,ez_space,tbody,&langextra[ibody][3]); } }