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

Merge pull request #2367 from ndtrung81/rigid-langevin

Browse Source 
Fixed a bug in computing the langevin torques applied to rigid bodies
This commit is contained in:
Axel Kohlmeyer
2020-09-16 22:36:04 -04:00
committed by GitHub
parent d611b85d7d dfcdb61d87
commit c14fd3131a
2 changed files with 34 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
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,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]);
}
}

23
src/RIGID/fix_rigid_small.cpp
View File

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