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

This commit is contained in:
sjplimp
2009-11-17 21:18:23 +00:00
parent 7fc817dc15
commit 4790c5402e
13 changed files with 218 additions and 160 deletions

View File

@ -119,6 +119,12 @@ void FixNPTAsphere::initial_integrate(int vflag)
} }
factor_rotate = exp(-dthalf*eta_dot); 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 **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; 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]; 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -157,7 +164,7 @@ void FixNPTAsphere::initial_integrate(int vflag)
remap(0); 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
@ -171,7 +178,7 @@ void FixNPTAsphere::initial_integrate(int vflag)
dtq = 0.5 * dtv; 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 // update quaternion a full step via Richardson iteration
// returns new normalized quaternion // returns new normalized quaternion
@ -199,7 +206,11 @@ void FixNPTAsphere::final_integrate()
int i; int i;
double dtfm; 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 **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -223,7 +234,8 @@ void FixNPTAsphere::final_integrate()
angmom[i][2] = (angmom[i][2] + dtf * torque[i][2]) * factor_rotate; 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);

View File

@ -101,6 +101,12 @@ void FixNVTAsphere::initial_integrate(int vflag)
eta += dtv*eta_dot; eta += dtv*eta_dot;
factor = exp(-dthalf*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 **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -173,6 +180,12 @@ void FixNVTAsphere::final_integrate()
{ {
double dtfm; 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 **v = atom->v;
double **f = atom->f; double **f = atom->f;
double **angmom = atom->angmom; 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);

View File

@ -195,7 +195,12 @@ void FixLangevin::post_force_no_tally()
double t_target = t_start + delta * (t_stop-t_start); double t_target = t_start + delta * (t_stop-t_start);
double tsqrt = sqrt(t_target); 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) { if (rmass) {
double boltz = force->boltz; 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) { } else if (which == BIAS) {
double tmp = temperature->compute_scalar(); double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
gamma1 = -rmass[i] / t_period / ftm2v; 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) { } else if (which == BIAS) {
double tmp = temperature->compute_scalar(); double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]]; gamma1 = gfactor1[type[i]];
@ -303,6 +300,11 @@ void FixLangevin::post_force_tally()
double tsqrt = sqrt(t_target); double tsqrt = sqrt(t_target);
// apply damping and thermostat to appropriate atoms // 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) { if (rmass) {
double boltz = force->boltz; 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) { } else if (which == BIAS) {
double tmp = temperature->compute_scalar(); double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
gamma1 = -rmass[i] / t_period / ftm2v; 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) { } else if (which == BIAS) {
double tmp = temperature->compute_scalar(); double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
gamma1 = gfactor1[type[i]]; gamma1 = gfactor1[type[i]];

View File

@ -388,6 +388,12 @@ void FixNPT::initial_integrate(int vflag)
dilation[i] = exp(dthalf*omega_dot[i]); 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 **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; 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]; 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); 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]; 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -449,7 +457,7 @@ void FixNPT::initial_integrate(int vflag)
remap(0); 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
@ -475,6 +483,12 @@ void FixNPT::final_integrate()
int i; int i;
double dtfm; 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 **v = atom->v;
double **f = atom->f; double **f = atom->f;
double *rmass = atom->rmass; double *rmass = atom->rmass;
@ -494,7 +508,8 @@ void FixNPT::final_integrate()
v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2]; 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); 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]; 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); 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 // outermost level - update eta_dot and omega_dot, apply to v, remap box
// all other levels - NVE update of v // 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) { 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]); dilation[i] = exp(dthalf*omega_dot[i]);
} }
// v update only for atoms in group // update v for atoms in group
if (rmass) { if (rmass) {
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
@ -663,57 +679,31 @@ void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
} else { } else {
// v update only for atoms in group // update v for atoms in group
if (rmass) { if (rmass) {
if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) {
if (mask[i] & groupbit) { dtfm = dtf / rmass[i];
dtfm = dtf / rmass[i]; v[i][0] += dtfm*f[i][0];
v[i][0] += dtfm*f[i][0]; v[i][1] += dtfm*f[i][1];
v[i][1] += dtfm*f[i][1]; v[i][2] += dtfm*f[i][2];
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]);
}
} }
} }
} else { } else {
if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) {
if (mask[i] & groupbit) { dtfm = dtf / mass[type[i]];
dtfm = dtf / mass[type[i]]; v[i][0] += dtfm*f[i][0];
v[i][0] += dtfm*f[i][0]; v[i][1] += dtfm*f[i][1];
v[i][1] += dtfm*f[i][1]; v[i][2] += dtfm*f[i][2];
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]);
}
} }
} }
} }
} }
// innermost level - also update x only for atoms in group // innermost level - also update x for atoms in group
if (ilevel == 0) { if (ilevel == 0) {
for (int i = 0; i < nlocal; i++) { 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, // outermost level - update eta_dot and omega_dot,
// apply to v via final_integrate() // apply to v via final_integrate()
// all other levels - NVE update of v // 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(); if (ilevel == nlevels_respa-1) final_integrate();
else { else {
// v update only for atoms in group // update v for atoms in group
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -757,48 +747,22 @@ void FixNPT::final_integrate_respa(int ilevel)
if (igroup == atom->firstgroup) nlocal = atom->nfirst; if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (rmass) { if (rmass) {
if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) {
if (mask[i] & groupbit) { dtfm = dtf / rmass[i];
dtfm = dtf / rmass[i]; v[i][0] += dtfm*f[i][0];
v[i][0] += dtfm*f[i][0]; v[i][1] += dtfm*f[i][1];
v[i][1] += dtfm*f[i][1]; v[i][2] += dtfm*f[i][2];
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]);
}
} }
} }
} else { } else {
if (which == NOBIAS) { for (int i = 0; i < nlocal; i++) {
for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) {
if (mask[i] & groupbit) { dtfm = dtf / mass[type[i]];
dtfm = dtf / mass[type[i]]; v[i][0] += dtfm*f[i][0];
v[i][0] += dtfm*f[i][0]; v[i][1] += dtfm*f[i][1];
v[i][1] += dtfm*f[i][1]; v[i][2] += dtfm*f[i][2];
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]);
}
} }
} }
} }

View File

@ -122,6 +122,12 @@ void FixNPTSphere::initial_integrate(int vflag)
} }
factor_rotate = exp(-dthalf*eta_dot); 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 **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; 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]; 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); 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]; 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -187,7 +195,7 @@ void FixNPTSphere::initial_integrate(int vflag)
remap(0); 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++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
@ -282,9 +290,13 @@ void FixNPTSphere::final_integrate()
double dtfrotate = dtf / INERTIA; double dtfrotate = dtf / INERTIA;
// update v,omega for all particles // update v,omega of atoms in group
// d_omega/dt = torque / inertia // d_omega/dt = torque / inertia
// 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias // 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 (radius) {
if (rmass) { if (rmass) {
@ -305,7 +317,8 @@ void FixNPTSphere::final_integrate()
factor_rotate; factor_rotate;
} }
} }
} else if (which == BIAS) { } else {
double tmp = temperature->compute_scalar();
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -345,8 +358,9 @@ void FixNPTSphere::final_integrate()
factor_rotate; factor_rotate;
} }
} }
} else if (which == BIAS) { } else {
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
double tmp = temperature->compute_scalar();
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
itype = type[i]; itype = type[i];
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -388,7 +402,8 @@ void FixNPTSphere::final_integrate()
factor_rotate; factor_rotate;
} }
} }
} else if (which == BIAS) { } else {
double tmp = temperature->compute_scalar();
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
itype = type[i]; itype = type[i];
@ -430,7 +445,8 @@ void FixNPTSphere::final_integrate()
factor_rotate; factor_rotate;
} }
} }
} else if (which == BIAS) { } else {
double tmp = temperature->compute_scalar();
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
itype = type[i]; itype = type[i];

View File

@ -64,6 +64,8 @@ void FixNVE::initial_integrate(int vflag)
{ {
double dtfm; double dtfm;
// update v and x of atoms in group
double **x = atom->x; double **x = atom->x;
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -108,6 +110,8 @@ void FixNVE::final_integrate()
{ {
double dtfm; double dtfm;
// update v of atoms in group
double **v = atom->v; double **v = atom->v;
double **f = atom->f; double **f = atom->f;
double *mass = atom->mass; double *mass = atom->mass;

View File

@ -200,7 +200,11 @@ void FixNVT::initial_integrate(int vflag)
factor = exp(-dthalf*eta_dot); 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 **x = atom->x;
double **v = atom->v; double **v = atom->v;
@ -223,7 +227,8 @@ void FixNVT::initial_integrate(int vflag)
x[i][2] += dtv * v[i][2]; x[i][2] += dtv * v[i][2];
} }
} }
} else if (which == BIAS) { } else {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -252,7 +257,8 @@ void FixNVT::initial_integrate(int vflag)
x[i][2] += dtv * v[i][2]; x[i][2] += dtv * v[i][2];
} }
} }
} else if (which == BIAS) { } else {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -284,7 +290,11 @@ void FixNVT::final_integrate()
{ {
double dtfm; 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 **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -307,7 +317,8 @@ void FixNVT::final_integrate()
v[i][2] = v[i][2]*factor + dtfm*f[i][2]; 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); 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]; 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -369,19 +381,20 @@ void FixNVT::final_integrate()
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]]; dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0] * factor; v[i][0] *= factor;
v[i][1] = v[i][1] * factor; v[i][1] *= factor;
v[i][2] = v[i][2] * factor; v[i][2] *= factor;
} }
} }
} else if (which == BIAS) { } else {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]]; dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0] * factor; v[i][0] *= factor;
v[i][1] = v[i][1] * factor; v[i][1] *= factor;
v[i][2] = v[i][2] * factor; v[i][2] *= factor;
temperature->restore_bias(i,v[i]); temperature->restore_bias(i,v[i]);
} }
} }
@ -392,19 +405,20 @@ void FixNVT::final_integrate()
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]]; dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0] * factor; v[i][0] *= factor;
v[i][1] = v[i][1] * factor; v[i][1] *= factor;
v[i][2] = v[i][2] * factor; v[i][2] *= factor;
} }
} }
} else if (which == BIAS) { } else {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]]; dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0] * factor; v[i][0] *= factor;
v[i][1] = v[i][1] * factor; v[i][1] *= factor;
v[i][2] = v[i][2] * factor; v[i][2] *= factor;
temperature->restore_bias(i,v[i]); 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]; 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);

View File

@ -35,7 +35,6 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{NO_REMAP,X_REMAP,V_REMAP}; // same as fix_deform.cpp 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) if (!temperature->tempbias)
error->all("Temperature for fix nvt/sllod does not have a bias"); 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 // check fix deform remap settings
int i; int i;
for (i = 0; i < modify->nfix; i++) for (i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"deform") == 0) { if (strcmp(modify->fix[i]->style,"deform") == 0) {
if (((FixDeform *) modify->fix[i])->remapflag != V_REMAP) 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; break;
} }
if (i == modify->nfix) if (i == modify->nfix)
@ -82,10 +85,15 @@ void FixNVTSlodd::initial_integrate(int vflag)
eta += dtv*eta_dot; eta += dtv*eta_dot;
factor = exp(-dthalf*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 // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo
// thermostat thermal velocity only // thermostat thermal velocity only
// vdelu = SLLOD correction = Hrate*Hinv*vthermal // 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 **x = atom->x;
double **v = atom->v; double **v = atom->v;
@ -124,10 +132,16 @@ void FixNVTSlodd::final_integrate()
{ {
double dtfm; double dtfm;
// update v of only atoms in group // update v of atoms in group
// remove and restore bias = streaming velocity = Hrate*lamda + Hratelo // remove and restore bias = streaming velocity = Hrate*lamda + Hratelo
// thermostat thermal velocity only // thermostat thermal velocity only
// vdelu = SLLOD correction = Hrate*Hinv*vthermal // 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 **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -207,7 +221,9 @@ void FixNVTSlodd::initial_integrate_respa(int vflag, int ilevel, int flag)
factor = exp(-dthalf*eta_dot); factor = exp(-dthalf*eta_dot);
} else factor = 1.0; } 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]; double h_two[6],vdelu[3];
MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two); MathExtra::multiply_shape_shape(domain->h_rate,domain->h_inv,h_two);

View File

@ -25,6 +25,9 @@ class FixNVTSlodd : public FixNVT {
void initial_integrate(int); void initial_integrate(int);
void final_integrate(); void final_integrate();
void initial_integrate_respa(int, int, int); void initial_integrate_respa(int, int, int);
private:
int nondeformbias;
}; };
} }

View File

@ -125,9 +125,13 @@ void FixNVTSphere::initial_integrate(int vflag)
double dtfrotate = dtf / INERTIA; double dtfrotate = dtf / INERTIA;
// update v,x,omega for all particles // update v,x,omega of atoms in group
// d_omega/dt = torque / inertia // d_omega/dt = torque / inertia
// 8 cases depending on radius vs shape, rmass vs mass, bias vs nobias // 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 (radius) {
if (rmass) { if (rmass) {
@ -148,7 +152,8 @@ void FixNVTSphere::initial_integrate(int vflag)
omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); 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]; 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
itype = type[i]; itype = type[i];
@ -231,7 +237,8 @@ void FixNVTSphere::initial_integrate(int vflag)
omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
itype = type[i]; itype = type[i];
@ -273,7 +280,8 @@ void FixNVTSphere::initial_integrate(int vflag)
omega[i][2] = omega[i][2]*factor + dtirotate*torque[i][2]; 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++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
itype = type[i]; itype = type[i];
@ -306,7 +314,13 @@ void FixNVTSphere::final_integrate()
int itype; int itype;
double dtfm,dtirotate; 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 **v = atom->v;
double **f = atom->f; double **f = atom->f;
@ -325,10 +339,6 @@ void FixNVTSphere::final_integrate()
double dtfrotate = dtf / INERTIA; 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 (radius) {
if (rmass) { if (rmass) {
if (which == NOBIAS) { if (which == NOBIAS) {
@ -346,6 +356,7 @@ void FixNVTSphere::final_integrate()
} }
} }
} else { } else {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);
@ -380,6 +391,7 @@ void FixNVTSphere::final_integrate()
} }
} }
} else { } else {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
itype = type[i]; itype = type[i];
@ -417,6 +429,7 @@ void FixNVTSphere::final_integrate()
} }
} }
} else { } else {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
itype = type[i]; itype = type[i];
@ -453,6 +466,7 @@ void FixNVTSphere::final_integrate()
} }
} }
} else { } else {
double tmp = temperature->compute_scalar();
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
itype = type[i]; itype = type[i];

View File

@ -109,6 +109,9 @@ void FixTempBerendsen::end_of_step()
t_target = t_start + delta * (t_stop-t_start); t_target = t_start + delta * (t_stop-t_start);
// rescale velocities by lamda // 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)); 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; v[i][2] *= lamda;
} }
} }
} else if (which == BIAS) { } else {
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]); temperature->remove_bias(i,v[i]);

View File

@ -114,6 +114,9 @@ void FixTempRescale::end_of_step()
double t_target = t_start + delta * (t_stop-t_start); double t_target = t_start + delta * (t_stop-t_start);
// rescale velocity of appropriate atoms if outside window // 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) { if (fabs(t_current-t_target) > t_window) {
t_target = t_current - fraction*(t_current-t_target); t_target = t_current - fraction*(t_current-t_target);
@ -133,8 +136,7 @@ void FixTempRescale::end_of_step()
v[i][2] *= factor; v[i][2] *= factor;
} }
} }
} else {
} else if (which == BIAS) {
energy += (t_current-t_target) * efactor; energy += (t_current-t_target) * efactor;
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
@ -178,7 +180,6 @@ int FixTempRescale::modify_param(int narg, char **arg)
return 0; return 0;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void FixTempRescale::reset_target(double t_new) void FixTempRescale::reset_target(double t_new)

View File

@ -41,6 +41,7 @@ using namespace LAMMPS_NS;
LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator)
{ {
printf("AAA\n");
memory = new Memory(this); memory = new Memory(this);
error = new Error(this); error = new Error(this);
universe = new Universe(this,communicator); universe = new Universe(this,communicator);