Hello,
I need to rescale velocities and angular velocities with the Berendsen
thermostat, therfore I have modified "void
FixTempBerendsen::end_of_step()" and included to functions in the
"compute_temp_sphere" files:
double ComputeTempSphere::compute_scalar_rot_T()
double ComputeTempSphere::compute_scalar_trans_T()
which are both copies of
double ComputeTempSphere::compute_scalar(), but I have deleted the
rotational contribution to the temperature in "compute_scalar_trans_T" and
the translational in the other case. Of courese the dof are wrong so far,
but my problem is bigger:
The function void FixTempBerendsen::end_of_step() throws the
"if (t_current == 0.0)
error->all("Computed temperature for fix temp/berendsen cannot be
0.0")" event. I do not know why! I also included the
compute_scalar_rot_T()
compute_scalar_trans_T()
as virtual in the "compute.h" file just like the compute_scalar().
Here the modified files:
void FixTempBerendsen::end_of_step()
{
// double t_current = temperature->compute_scalar(); //modified
double t_current_rot = temperature->compute_scalar_rot_T(); //modified
double t_current = temperature->compute_scalar_trans_T(); //modified
if (t_current == 0.0)
error->all("Computed temperature for fix temp/berendsen cannot be 0.0
- why");
if (t_current_rot < 1.E-8) //modified
{t_current_rot = 1.E-8;} //modified
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
double 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));
double efactor = 0.5 * force->boltz * temperature->dof;
energy += t_current * (1.0-lamda*lamda) * efactor;
double **v = atom->v;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
v[i][0] *= lamda;
v[i][1] *= lamda;
v[i][2] *= lamda;
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
v[i][0] *= lamda;
v[i][1] *= lamda;
v[i][2] *= lamda;
temperature->restore_bias(i,v[i]);
}
}
}
/* Begin modified */
/*lamda = sqrt(1.0 + update->dt/t_period*(t_target/t_current_rot - 1.0));
double **omega = atom->omega;
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
omega[i][0] *= lamda;
omega[i][1] *= lamda;
omega[i][2] *= lamda;
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
//temperature->remove_bias(i,v[i]);
omega[i][0] *= lamda;
omega[i][1] *= lamda;
omega[i][2] *= lamda;
//temperature->restore_bias(i,v[i]);
}
}
}
/* End modified */
}
double ComputeTempSphere::compute_scalar_rot_T()
{
int i,itype;
invoked_scalar = update->ntimestep;
if (tempbias) {
if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar();
tbias->remove_bias_all();
}
double **v = atom->v;
double **omega = atom->omega;
double *radius = atom->radius;
double *rmass = atom->rmass;
double *mass = atom->mass;
double **shape = atom->shape;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// 4 cases depending on radius vs shape and rmass vs mass
// point particles will not contribute rotation due to radius or shape = 0
double t = 0.0;
if (radius) {
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) *
INERTIA*radius[i]*radius[i]*rmass[i];
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
itype = type[i];
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) *
INERTIA*radius[i]*radius[i]*mass[itype];
}
}
} else {
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
itype = type[i];
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) *
INERTIA*shape[itype][0]*shape[itype][0]*rmass[i];
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
itype = type[i];
t += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] +
omega[i][2]*omega[i][2]) *
INERTIA*shape[itype][0]*shape[itype][0]*mass[itype];
}
}
}
if (tempbias) tbias->restore_bias_all();
MPI_Allreduce(&t,&scalar_rot,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic || tempbias == 2) dof_compute();
scalar_rot *= 1. / ( 0.4 * ( 5. + (1. / tfactor) ) - 1. ); // tfactor =
... 1 / dof ... // Annahme: 5N - 5 , wollen 2N - 1
return scalar_rot;
}
double ComputeTempSphere::compute_scalar_trans_T()
{
int i,itype;
invoked_scalar = update->ntimestep;
if (tempbias) {
if (tbias->invoked_scalar != update->ntimestep) tbias->compute_scalar();
tbias->remove_bias_all();
}
double **v = atom->v;
double **omega = atom->omega;
double *radius = atom->radius;
double *rmass = atom->rmass;
double *mass = atom->mass;
double **shape = atom->shape;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// 4 cases depending on radius vs shape and rmass vs mass
// point particles will not contribute rotation due to radius or shape = 0
double t = 0.0;
if (radius) {
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
rmass[i];
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
itype = type[i];
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[itype];
}
}
} else {
if (rmass) {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
itype = type[i];
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
rmass[i];
}
} else {
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
itype = type[i];
t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) *
mass[itype];
}
}
}
if (tempbias) tbias->restore_bias_all();
MPI_Allreduce(&t,&scalar_trans,1,MPI_DOUBLE,MPI_SUM,world);
if (dynamic || tempbias == 2) dof_compute();
scalar_trans *= 1. / ( 0.6 * ( 5. + (1. / tfactor) ) - 1. ); // tfactor
= ... 1 / dof ... // Annahme: 5N - 5 // wollen 3N - 4
return scalar_trans;
}
scalar_trans and scalar_rot are defined in the header
Best
G. Rosenthal