From cccd6c441c93a7fe1ddd79f2038c80d5ac1a9586 Mon Sep 17 00:00:00 2001 From: Charlles Abreu Date: Tue, 5 Sep 2017 23:27:26 -0300 Subject: [PATCH] Conditional late computation added to rigid_nh and rigid_nh_small --- src/RIGID/fix_rigid_nh.cpp | 78 ++------------------------------ src/RIGID/fix_rigid_nh_small.cpp | 76 ++----------------------------- 2 files changed, 8 insertions(+), 146 deletions(-) diff --git a/src/RIGID/fix_rigid_nh.cpp b/src/RIGID/fix_rigid_nh.cpp index e67f376135..a182481d06 100644 --- a/src/RIGID/fix_rigid_nh.cpp +++ b/src/RIGID/fix_rigid_nh.cpp @@ -620,87 +620,17 @@ void FixRigidNH::final_integrate() akin_t = akin_r = 0.0; } - // sum over atoms to get force and torque on rigid body + // late calculation of forces and torques (if requested) - double **x = atom->x; - double **f = atom->f; - int nlocal = atom->nlocal; - - double xprd = domain->xprd; - double yprd = domain->yprd; - double zprd = domain->zprd; - if (triclinic) { - xy = domain->xy; - xz = domain->xz; - yz = domain->yz; + if (!earlyflag) { + if (langflag) apply_langevin_thermostat(); + compute_forces_and_torques(); } - int xbox,ybox,zbox; - double xunwrap,yunwrap,zunwrap,dx,dy,dz; - for (ibody = 0; ibody < nbody; ibody++) - for (i = 0; i < 6; i++) sum[ibody][i] = 0.0; - - for (i = 0; i < nlocal; i++) { - if (body[i] < 0) continue; - ibody = body[i]; - - sum[ibody][0] += f[i][0]; - sum[ibody][1] += f[i][1]; - sum[ibody][2] += f[i][2]; - - xbox = (xcmimage[i] & IMGMASK) - IMGMAX; - ybox = (xcmimage[i] >> IMGBITS & IMGMASK) - IMGMAX; - zbox = (xcmimage[i] >> IMG2BITS) - IMGMAX; - - if (triclinic == 0) { - xunwrap = x[i][0] + xbox*xprd; - yunwrap = x[i][1] + ybox*yprd; - zunwrap = x[i][2] + zbox*zprd; - } else { - xunwrap = x[i][0] + xbox*xprd + ybox*xy + zbox*xz; - yunwrap = x[i][1] + ybox*yprd + zbox*yz; - zunwrap = x[i][2] + zbox*zprd; - } - - dx = xunwrap - xcm[ibody][0]; - dy = yunwrap - xcm[ibody][1]; - dz = zunwrap - xcm[ibody][2]; - - sum[ibody][3] += dy*f[i][2] - dz*f[i][1]; - sum[ibody][4] += dz*f[i][0] - dx*f[i][2]; - sum[ibody][5] += dx*f[i][1] - dy*f[i][0]; - } - - // extended particles add their torque to torque of body - - if (extended) { - double **torque_one = atom->torque; - - for (i = 0; i < nlocal; i++) { - if (body[i] < 0) continue; - ibody = body[i]; - - if (eflags[i] & TORQUE) { - sum[ibody][3] += torque_one[i][0]; - sum[ibody][4] += torque_one[i][1]; - sum[ibody][5] += torque_one[i][2]; - } - } - } - - MPI_Allreduce(sum[0],all[0],6*nbody,MPI_DOUBLE,MPI_SUM,world); - // update vcm and angmom - // include Langevin thermostat forces // fflag,tflag = 0 for some dimensions in 2d for (ibody = 0; ibody < nbody; ibody++) { - fcm[ibody][0] = all[ibody][0] + langextra[ibody][0]; - fcm[ibody][1] = all[ibody][1] + langextra[ibody][1]; - fcm[ibody][2] = all[ibody][2] + langextra[ibody][2]; - torque[ibody][0] = all[ibody][3] + langextra[ibody][3]; - torque[ibody][1] = all[ibody][4] + langextra[ibody][4]; - torque[ibody][2] = all[ibody][5] + langextra[ibody][5]; // update vcm by 1/2 step diff --git a/src/RIGID/fix_rigid_nh_small.cpp b/src/RIGID/fix_rigid_nh_small.cpp index 2b55c73ec6..c23e341e6b 100644 --- a/src/RIGID/fix_rigid_nh_small.cpp +++ b/src/RIGID/fix_rigid_nh_small.cpp @@ -643,79 +643,11 @@ void FixRigidNHSmall::final_integrate() scale_r *= exp(-dtq * (pdim * mtk_term2)); } - // sum over atoms to get force and torque on rigid body + // late calculation of forces and torques (if requested) - double **x = atom->x; - double **f = atom->f; - int nlocal = atom->nlocal; - - double dx,dy,dz; - double unwrap[3]; - double *xcm,*fcm,*tcm; - - for (ibody = 0; ibody < nlocal_body+nghost_body; ibody++) { - fcm = body[ibody].fcm; - fcm[0] = fcm[1] = fcm[2] = 0.0; - tcm = body[ibody].torque; - tcm[0] = tcm[1] = tcm[2] = 0.0; - } - - for (i = 0; i < nlocal; i++) { - if (atom2body[i] < 0) continue; - Body *b = &body[atom2body[i]]; - - fcm = b->fcm; - fcm[0] += f[i][0]; - fcm[1] += f[i][1]; - fcm[2] += f[i][2]; - - domain->unmap(x[i],xcmimage[i],unwrap); - xcm = b->xcm; - dx = unwrap[0] - xcm[0]; - dy = unwrap[1] - xcm[1]; - dz = unwrap[2] - xcm[2]; - - tcm = b->torque; - tcm[0] += dy*f[i][2] - dz*f[i][1]; - tcm[1] += dz*f[i][0] - dx*f[i][2]; - tcm[2] += dx*f[i][1] - dy*f[i][0]; - } - - // extended particles add their torque to torque of body - - if (extended) { - double **torque = atom->torque; - - for (i = 0; i < nlocal; i++) { - if (atom2body[i] < 0) continue; - - if (eflags[i] & TORQUE) { - tcm = body[atom2body[i]].torque; - tcm[0] += torque[i][0]; - tcm[1] += torque[i][1]; - tcm[2] += torque[i][2]; - } - } - } - - // reverse communicate fcm, torque of all bodies - - commflag = FORCE_TORQUE; - comm->reverse_comm_fix(this,6); - - // include Langevin thermostat forces and torques - - if (langflag) { - for (int ibody = 0; ibody < nlocal_body; ibody++) { - fcm = body[ibody].fcm; - fcm[0] += langextra[ibody][0]; - fcm[1] += langextra[ibody][1]; - fcm[2] += langextra[ibody][2]; - tcm = body[ibody].torque; - tcm[0] += langextra[ibody][3]; - tcm[1] += langextra[ibody][4]; - tcm[2] += langextra[ibody][5]; - } + if (!earlyflag) { + if (langflag) apply_langevin_thermostat(); + compute_forces_and_torques(); } // update vcm and angmom