jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12864 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2015-01-05 17:08:33 +00:00
parent a6c70d8a68
commit d194479a74
2 changed files with 7 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/RIGID/fix_rigid.cpp
View File

@ -1141,7 +1141,9 @@ int FixRigid::dof(int tgroup)
for (int i = 0; i < nlocal; i++)
if (body[i] >= 0 && mask[i] & tgroupbit) {
if (extended && eflags[i]) mcount[body[i]]++;
// do not count point particles or point dipoles as extended particles
// a spheroid dipole will be counted as extended
if (extended && (eflags[i] & ~(POINT | DIPOLE))) mcount[body[i]]++;
else ncount[body[i]]++;
}

50
src/RIGID/fix_rigid_nh.cpp
View File

@ -955,53 +955,11 @@ void FixRigidNH::nhc_press_integrate()
double FixRigidNH::compute_scalar()
{
int i,k,ibody;
double kt = boltz * t_target;
double energy,ke_t,ke_q,tmp,Pkq[4];
const double kt = boltz * t_target;
double energy;
int i;
// compute the kinetic parts of H_NVE in Kameraj et al (JCP 2005, pp 224114)
// translational kinetic energy
ke_t = 0.0;
for (ibody = 0; ibody < nbody; ibody++)
ke_t += 0.5 * masstotal[ibody] * (vcm[ibody][0]*vcm[ibody][0] +
vcm[ibody][1]*vcm[ibody][1] +
vcm[ibody][2]*vcm[ibody][2]);
// rotational kinetic energy
ke_q = 0.0;
for (ibody = 0; ibody < nbody; ibody++) {
for (k = 1; k < 4; k++) {
if (k == 1) {
Pkq[0] = -quat[ibody][1];
Pkq[1] = quat[ibody][0];
Pkq[2] = quat[ibody][3];
Pkq[3] = -quat[ibody][2];
} else if (k == 2) {
Pkq[0] = -quat[ibody][2];
Pkq[1] = -quat[ibody][3];
Pkq[2] = quat[ibody][0];
Pkq[3] = quat[ibody][1];
} else if (k == 3) {
Pkq[0] = -quat[ibody][3];
Pkq[1] = quat[ibody][2];
Pkq[2] = -quat[ibody][1];
Pkq[3] = quat[ibody][0];
}
tmp = conjqm[ibody][0]*Pkq[0] + conjqm[ibody][1]*Pkq[1] +
conjqm[ibody][2]*Pkq[2] + conjqm[ibody][3]*Pkq[3];
tmp *= tmp;
if (fabs(inertia[ibody][k-1]) < 1e-6) tmp = 0.0;
else tmp /= (8.0 * inertia[ibody][k-1]);
ke_q += tmp;
}
}
energy = (ke_t + ke_q) * mvv2e;
energy = FixRigid::compute_scalar();
if (tstat_flag) {