diff --git a/src/ASPHERE/fix_npt_asphere.cpp b/src/ASPHERE/fix_npt_asphere.cpp index f09635211c..eb4b1f8901 100755 --- a/src/ASPHERE/fix_npt_asphere.cpp +++ b/src/ASPHERE/fix_npt_asphere.cpp @@ -119,6 +119,12 @@ void FixNPTAsphere::initial_integrate(int vflag) } factor_rotate = exp(-dthalf*eta_dot); + // update v of atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v + double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -140,7 +146,8 @@ void FixNPTAsphere::initial_integrate(int vflag) v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -157,7 +164,7 @@ void FixNPTAsphere::initial_integrate(int vflag) remap(0); - // x update by full step only for atoms in group + // update x by full step for atoms in group for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -171,7 +178,7 @@ void FixNPTAsphere::initial_integrate(int vflag) dtq = 0.5 * dtv; - // update angular momentum by 1/2 step for all particles + // update angular momentum by 1/2 step for atoms in group // update quaternion a full step via Richardson iteration // returns new normalized quaternion @@ -199,7 +206,11 @@ void FixNPTAsphere::final_integrate() int i; double dtfm; - // update v of only atoms in group + // update v,angmom of atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v double **v = atom->v; double **f = atom->f; @@ -223,7 +234,8 @@ void FixNPTAsphere::final_integrate() angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * factor_rotate; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); diff --git a/src/ASPHERE/fix_nvt_asphere.cpp b/src/ASPHERE/fix_nvt_asphere.cpp index ce69e97694..d30db19d9e 100755 --- a/src/ASPHERE/fix_nvt_asphere.cpp +++ b/src/ASPHERE/fix_nvt_asphere.cpp @@ -101,6 +101,12 @@ void FixNVTAsphere::initial_integrate(int vflag) eta += dtv*eta_dot; factor = exp(-dthalf*eta_dot); + // update v,x,angmom of atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v + double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -140,7 +146,8 @@ void FixNVTAsphere::initial_integrate(int vflag) } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -173,6 +180,12 @@ void FixNVTAsphere::final_integrate() { double dtfm; + // update v,angmom of atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v + double **v = atom->v; double **f = atom->f; double **angmom = atom->angmom; @@ -196,7 +209,8 @@ void FixNVTAsphere::final_integrate() } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp index d633b171f7..b042de4f7d 100644 --- a/src/fix_langevin.cpp +++ b/src/fix_langevin.cpp @@ -195,7 +195,12 @@ void FixLangevin::post_force_no_tally() double t_target = t_start + delta * (t_stop-t_start); double tsqrt = sqrt(t_target); - // apply damping and thermostat to appropriate atoms + // apply damping and thermostat to atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // test v = 0 since some computes mask non-participating atoms via v = 0 + // and added force has extra term not multiplied by v = 0 if (rmass) { double boltz = force->boltz; @@ -216,12 +221,8 @@ void FixLangevin::post_force_no_tally() } } - // invoke temperature since some computes require it to remove bias - // test v = 0 since some computes mask non-participating atoms via v = 0 - } else if (which == BIAS) { double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { gamma1 = -rmass[i] / t_period / ftm2v; @@ -252,12 +253,8 @@ void FixLangevin::post_force_no_tally() } } - // invoke temperature since some computes require it to remove bias - // test v = 0 since some computes mask non-participating atoms via v = 0 - } else if (which == BIAS) { double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { gamma1 = gfactor1[type[i]]; @@ -303,6 +300,11 @@ void FixLangevin::post_force_tally() double tsqrt = sqrt(t_target); // apply damping and thermostat to appropriate atoms + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // test v = 0 since some computes mask non-participating atoms via v = 0 + // and added force has extra term not multiplied by v = 0 if (rmass) { double boltz = force->boltz; @@ -326,12 +328,8 @@ void FixLangevin::post_force_tally() } } - // invoke temperature since some computes require it to remove bias - // test v = 0 since some computes mask non-participating atoms via v = 0 - } else if (which == BIAS) { double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { gamma1 = -rmass[i] / t_period / ftm2v; @@ -368,12 +366,8 @@ void FixLangevin::post_force_tally() } } - // invoke temperature since some computes require it to remove bias - // test v = 0 since some computes mask non-participating atoms via v = 0 - } else if (which == BIAS) { double tmp = temperature->compute_scalar(); - for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { gamma1 = gfactor1[type[i]]; diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp index a8880e470f..2efc81854d 100644 --- a/src/fix_npt.cpp +++ b/src/fix_npt.cpp @@ -388,6 +388,12 @@ void FixNPT::initial_integrate(int vflag) dilation[i] = exp(dthalf*omega_dot[i]); } + // update v and x of atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v + double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -408,7 +414,8 @@ void FixNPT::initial_integrate(int vflag) v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -431,7 +438,8 @@ void FixNPT::initial_integrate(int vflag) v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -449,7 +457,7 @@ void FixNPT::initial_integrate(int vflag) remap(0); - // x update by full step only for atoms in group + // update x by full step for atoms in group for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -475,6 +483,12 @@ void FixNPT::final_integrate() int i; double dtfm; + // update v of atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v + double **v = atom->v; double **f = atom->f; double *rmass = atom->rmass; @@ -494,7 +508,8 @@ void FixNPT::final_integrate() v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -517,7 +532,8 @@ void FixNPT::final_integrate() v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -602,7 +618,7 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag) // outermost level - update eta_dot and omega_dot, apply to v, remap box // all other levels - NVE update of v - // x,v updates only performed for atoms in group + // x,v updates performed for atoms in group if (ilevel == nlevels_respa-1) { @@ -635,7 +651,7 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag) dilation[i] = exp(dthalf*omega_dot[i]); } - // v update only for atoms in group + // update v for atoms in group if (rmass) { for (int i = 0; i < nlocal; i++) { @@ -663,57 +679,31 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag) } else { - // v update only for atoms in group + // update v for atoms in group if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - } - } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[i]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); - } + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; } } } else { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - } - } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); - } + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; } } } } - // innermost level - also update x only for atoms in group + // innermost level - also update x for atoms in group if (ilevel == 0) { for (int i = 0; i < nlocal; i++) { @@ -740,12 +730,12 @@ void FixNPT::final_integrate_respa(int ilevel) // outermost level - update eta_dot and omega_dot, // apply to v via final_integrate() // all other levels - NVE update of v - // v update only performed for atoms in group + // update v for atoms in group if (ilevel == nlevels_respa-1) final_integrate(); else { - // v update only for atoms in group + // update v for atoms in group double **v = atom->v; double **f = atom->f; @@ -757,48 +747,22 @@ void FixNPT::final_integrate_respa(int ilevel) if (igroup == atom->firstgroup) nlocal = atom->nfirst; if (rmass) { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / rmass[i]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - } - } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / rmass[i]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); - } + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; } } } else { - if (which == NOBIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - } - } - } else if (which == BIAS) { - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); - dtfm = dtf / mass[type[i]]; - v[i][0] += dtfm*f[i][0]; - v[i][1] += dtfm*f[i][1]; - v[i][2] += dtfm*f[i][2]; - temperature->restore_bias(i,v[i]); - } + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / mass[type[i]]; + v[i][0] += dtfm*f[i][0]; + v[i][1] += dtfm*f[i][1]; + v[i][2] += dtfm*f[i][2]; } } } diff --git a/src/fix_npt_sphere.cpp b/src/fix_npt_sphere.cpp index ac0acc26ef..3ba3bf9861 100644 --- a/src/fix_npt_sphere.cpp +++ b/src/fix_npt_sphere.cpp @@ -122,6 +122,12 @@ void FixNPTSphere::initial_integrate(int vflag) } factor_rotate = exp(-dthalf*eta_dot); + // update v of atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v + double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -146,7 +152,8 @@ void FixNPTSphere::initial_integrate(int vflag) v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -169,7 +176,8 @@ void FixNPTSphere::initial_integrate(int vflag) v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -187,7 +195,7 @@ void FixNPTSphere::initial_integrate(int vflag) remap(0); - // x update by full step only for atoms in group + // update x by full step for atoms in group for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -282,9 +290,13 @@ void FixNPTSphere::final_integrate() double dtfrotate = dtf / INERTIA; - // update v,omega for all particles + // update v,omega of atoms in group // d_omega/dt = torque / inertia // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v if (radius) { if (rmass) { @@ -305,7 +317,8 @@ void FixNPTSphere::final_integrate() factor_rotate; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -345,8 +358,9 @@ void FixNPTSphere::final_integrate() factor_rotate; } } - } else if (which == BIAS) { + } else { for (i = 0; i < nlocal; i++) { + double tmp = temperature->compute_scalar(); if (mask[i] & groupbit) { itype = type[i]; temperature->remove_bias(i,v[i]); @@ -388,7 +402,8 @@ void FixNPTSphere::final_integrate() factor_rotate; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; @@ -430,7 +445,8 @@ void FixNPTSphere::final_integrate() factor_rotate; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; diff --git a/src/fix_nve.cpp b/src/fix_nve.cpp index f8789c276e..3eb8c39a62 100644 --- a/src/fix_nve.cpp +++ b/src/fix_nve.cpp @@ -64,6 +64,8 @@ void FixNVE::initial_integrate(int vflag) { double dtfm; + // update v and x of atoms in group + double **x = atom->x; double **v = atom->v; double **f = atom->f; @@ -108,6 +110,8 @@ void FixNVE::final_integrate() { double dtfm; + // update v of atoms in group + double **v = atom->v; double **f = atom->f; double *mass = atom->mass; diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp index 44531aad16..d6e2364ce0 100644 --- a/src/fix_nvt.cpp +++ b/src/fix_nvt.cpp @@ -200,7 +200,11 @@ void FixNVT::initial_integrate(int vflag) factor = exp(-dthalf*eta_dot); - // update v and x of only atoms in group + // update v and x of atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v double **x = atom->x; double **v = atom->v; @@ -223,7 +227,8 @@ void FixNVT::initial_integrate(int vflag) x[i][2] += dtv * v[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -252,7 +257,8 @@ void FixNVT::initial_integrate(int vflag) x[i][2] += dtv * v[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -284,7 +290,11 @@ void FixNVT::final_integrate() { double dtfm; - // update v of only atoms in group + // update v of atoms in group + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v double **v = atom->v; double **f = atom->f; @@ -307,7 +317,8 @@ void FixNVT::final_integrate() v[i][2] = v[i][2]*factor + dtfm*f[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -330,7 +341,8 @@ void FixNVT::final_integrate() v[i][2] = v[i][2]*factor + dtfm*f[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -369,19 +381,20 @@ void FixNVT::final_integrate() for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0] * factor; - v[i][1] = v[i][1] * factor; - v[i][2] = v[i][2] * factor; + v[i][0] *= factor; + v[i][1] *= factor; + v[i][2] *= factor; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0] * factor; - v[i][1] = v[i][1] * factor; - v[i][2] = v[i][2] * factor; + v[i][0] *= factor; + v[i][1] *= factor; + v[i][2] *= factor; temperature->restore_bias(i,v[i]); } } @@ -392,19 +405,20 @@ void FixNVT::final_integrate() for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0] * factor; - v[i][1] = v[i][1] * factor; - v[i][2] = v[i][2] * factor; + v[i][0] *= factor; + v[i][1] *= factor; + v[i][2] *= factor; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); dtfm = dtf / mass[type[i]]; - v[i][0] = v[i][0] * factor; - v[i][1] = v[i][1] * factor; - v[i][2] = v[i][2] * factor; + v[i][0] *= factor; + v[i][1] *= factor; + v[i][2] *= factor; temperature->restore_bias(i,v[i]); } } @@ -492,7 +506,8 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag) v[i][2] = v[i][2]*factor + dtfm*f[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -516,7 +531,8 @@ void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag) } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); diff --git a/src/fix_nvt_sllod.cpp b/src/fix_nvt_sllod.cpp index 5f98838d4e..c2e61d1242 100644 --- a/src/fix_nvt_sllod.cpp +++ b/src/fix_nvt_sllod.cpp @@ -35,7 +35,6 @@ using namespace LAMMPS_NS; enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp -enum{NOBIAS,BIAS}; /* ---------------------------------------------------------------------- */ @@ -51,13 +50,17 @@ void FixNVTSlodd::init() if (!temperature->tempbias) error->all("Temperature for fix nvt/sllod does not have a bias"); + nondeformbias = 0; + if (strcmp(temperature->style,"temp/deform") != 0) nondeformbias = 1; + // check fix deform remap settings int i; for (i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"deform") == 0) { if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP) - error->all("Using fix nvt/sllod with inconsistent fix deform remap option"); + error->all("Using fix nvt/sllod with inconsistent fix deform " + "remap option"); break; } if (i == modify->nfix) @@ -82,10 +85,15 @@ void FixNVTSlodd::initial_integrate(int vflag) eta += dtv*eta_dot; factor = exp(-dthalf*eta_dot); - // update v and x of only atoms in group + // update v and x of atoms in group // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo // thermostat thermal velocity only // vdelu = SLLOD correction = Hrate*Hinv*vthermal + // for non temp/deform BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + + if (nondeformbias) double tmp = temperature->compute_scalar(); double **x = atom->x; double **v = atom->v; @@ -124,10 +132,16 @@ void FixNVTSlodd::final_integrate() { double dtfm; - // update v of only atoms in group + // update v of atoms in group // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo // thermostat thermal velocity only // vdelu = SLLOD correction = Hrate*Hinv*vthermal + // for non temp/deform BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v + + if (nondeformbias) double tmp = temperature->compute_scalar(); double **v = atom->v; double **f = atom->f; @@ -207,7 +221,9 @@ void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag) factor = exp(-dthalf*eta_dot); } else factor = 1.0; - // update v of only atoms in group + // update v of atoms in group + + if (nondeformbias) double tmp = temperature->compute_scalar(); double h_two[6],vdelu[3]; MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two); diff --git a/src/fix_nvt_sllod.h b/src/fix_nvt_sllod.h index b0a5955c22..ba47f77056 100644 --- a/src/fix_nvt_sllod.h +++ b/src/fix_nvt_sllod.h @@ -25,6 +25,9 @@ class FixNVTSlodd : public FixNVT { void initial_integrate(int); void final_integrate(); void initial_integrate_respa(int, int, int); + + private: + int nondeformbias; }; } diff --git a/src/fix_nvt_sphere.cpp b/src/fix_nvt_sphere.cpp index 3a6d9f45a4..be4f7d7740 100644 --- a/src/fix_nvt_sphere.cpp +++ b/src/fix_nvt_sphere.cpp @@ -125,9 +125,13 @@ void FixNVTSphere::initial_integrate(int vflag) double dtfrotate = dtf / INERTIA; - // update v,x,omega for all particles + // update v,x,omega of atoms in group // d_omega/dt = torque / inertia // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v if (radius) { if (rmass) { @@ -148,7 +152,8 @@ void FixNVTSphere::initial_integrate(int vflag) omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -188,7 +193,8 @@ void FixNVTSphere::initial_integrate(int vflag) omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; @@ -231,7 +237,8 @@ void FixNVTSphere::initial_integrate(int vflag) omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; @@ -273,7 +280,8 @@ void FixNVTSphere::initial_integrate(int vflag) omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; } } - } else if (which == BIAS) { + } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; @@ -306,7 +314,13 @@ void FixNVTSphere::final_integrate() int itype; double dtfm,dtirotate; - // update v of only atoms in group + // update v,omega of atoms in group + // d_omega/dt = torque / inertia + // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias + // for BIAS: + // calculate temperature since some computes require temp + // computed on current nlocal atoms to remove bias + // OK to not test returned v = 0, since factor is multiplied by v double **v = atom->v; double **f = atom->f; @@ -325,10 +339,6 @@ void FixNVTSphere::final_integrate() double dtfrotate = dtf / INERTIA; - // update v,omega for all particles - // d_omega/dt = torque / inertia - // 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias - if (radius) { if (rmass) { if (which == NOBIAS) { @@ -346,6 +356,7 @@ void FixNVTSphere::final_integrate() } } } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); @@ -380,6 +391,7 @@ void FixNVTSphere::final_integrate() } } } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; @@ -417,6 +429,7 @@ void FixNVTSphere::final_integrate() } } } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; @@ -453,6 +466,7 @@ void FixNVTSphere::final_integrate() } } } else { + double tmp = temperature->compute_scalar(); for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { itype = type[i]; diff --git a/src/fix_temp_berendsen.cpp b/src/fix_temp_berendsen.cpp index 5c8a7cbe84..adf260d3be 100644 --- a/src/fix_temp_berendsen.cpp +++ b/src/fix_temp_berendsen.cpp @@ -109,6 +109,9 @@ void FixTempBerendsen::end_of_step() t_target = t_start + delta * (t_stop-t_start); // rescale velocities by lamda + // for BIAS: + // temperature is current, so do not need to re-compute + // OK to not test returned v = 0, since lamda is multiplied by v double lamda = sqrt(1.0 + update->dt/t_period*(t_target/t_current - 1.0)); @@ -124,7 +127,7 @@ void FixTempBerendsen::end_of_step() v[i][2] *= lamda; } } - } else if (which == BIAS) { + } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { temperature->remove_bias(i,v[i]); diff --git a/src/fix_temp_rescale.cpp b/src/fix_temp_rescale.cpp index 8bf8d79a04..e738692ebf 100644 --- a/src/fix_temp_rescale.cpp +++ b/src/fix_temp_rescale.cpp @@ -114,6 +114,9 @@ void FixTempRescale::end_of_step() double t_target = t_start + delta * (t_stop-t_start); // rescale velocity of appropriate atoms if outside window + // for BIAS: + // temperature is current, so do not need to re-compute + // OK to not test returned v = 0, since factor is multiplied by v if (fabs(t_current-t_target) > t_window) { t_target = t_current - fraction*(t_current-t_target); @@ -133,8 +136,7 @@ void FixTempRescale::end_of_step() v[i][2] *= factor; } } - - } else if (which == BIAS) { + } else { energy += (t_current-t_target) * efactor; for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { @@ -178,7 +180,6 @@ int FixTempRescale::modify_param(int narg, char **arg) return 0; } - /* ---------------------------------------------------------------------- */ void FixTempRescale::reset_target(double t_new) diff --git a/src/lammps.cpp b/src/lammps.cpp index 77babb81f0..8bb7cf9f70 100644 --- a/src/lammps.cpp +++ b/src/lammps.cpp @@ -41,6 +41,7 @@ using namespace LAMMPS_NS; LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) { + printf("AAA\n"); memory = new Memory(this); error = new Error(this); universe = new Universe(this,communicator);