ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Fixed a bug in computing the langevin torques applied to rigid bodies

Browse Source
This commit is contained in:
Trung Nguyen
2020-09-15 15:27:24 -05:00
parent de777ce994
commit ae68becf4a
2 changed files with 54 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

29
src/RIGID/fix_rigid.cpp
View File

@ -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];
}
}