Bug in Compute temp/chunk

Dear Developers,

I believe there is a bug in the remove_bias() and restore_bias() functions in ComputeTempChunk. This issue will cause a seg fault when using a temp/chunk with keyword “com yes” as the temperature compute in fix nvt. I think it’s just an off by one error:

L692 and L729 of compute_temp_chunk.cpp (May 14 2016 version):
int index = cchunk->ichunk[i];

should be changed to:
int index = cchunk->ichunk[i]-1;

The same correction may also be needed for remove_bias_all() and restore_bias_all() on L714 and L752. Also, I’m a little confused about whether these two functions are correct. As far as I can tell, vbias is never set in temp/chunk. I may be mistaken but I think L716, for example, should be changed to:

v[i][0] -= vcmall[index][0];

Thanks,
-David

hi david,

thanks a lot for the report. for future reference, i have filed it as
issue #74 in the LAMMPS project issue tracker on github:
https://github.com/lammps/lammps/issues/74

i think you are right with the off-by-one, but i will take a closer
look and also check out the other issues.

please stay tuned,
     axel.

hi again,

it looks like all the issues you list need to be addressed, and there
is also a potential division by zero problem for empty chunks.
the following patch should correct it.

diff --git a/src/compute_temp_chunk.cpp b/src/compute_temp_chunk.cpp
index 5a58535..87d6e11 100644
--- a/src/compute_temp_chunk.cpp
+++ b/src/compute_temp_chunk.cpp
@@ -472,9 +472,13 @@ void ComputeTempChunk::vcm_compute()
   MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world);

   for (i = 0; i < nchunk; i++) {
- vcmall[i][0] /= masstotal[i];
- vcmall[i][1] /= masstotal[i];
- vcmall[i][2] /= masstotal[i];
+ if (masstotal[i] > 0.0) {
+ vcmall[i][0] /= masstotal[i];
+ vcmall[i][1] /= masstotal[i];
+ vcmall[i][2] /= masstotal[i];
+ } else {
+ vcmall[i][0] = vcmall[i][1] = vcmall[i][2] = 0.0;
+ }
   }
}

@@ -689,7 +693,7 @@ void ComputeTempChunk::internal(int icol)

void ComputeTempChunk::remove_bias(int i, double *v)
{
- int index = cchunk->ichunk[i];
+ int index = cchunk->ichunk[i]-1;
   if (index < 0) return;
   v[0] -= vcmall[index][0];
   v[1] -= vcmall[index][1];
@@ -711,11 +715,11 @@ void ComputeTempChunk::remove_bias_all()

   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
- index = ichunk[i];
+ index = ichunk[i]-1;
       if (index < 0) continue;
- v[i][0] -= vbias[0];
- v[i][1] -= vbias[1];
- v[i][2] -= vbias[2];
+ v[i][0] -= vcmall[index][0];
+ v[i][1] -= vcmall[index][1];
+ v[i][2] -= vcmall[index][2];
     }
}

@@ -726,7 +730,7 @@ void ComputeTempChunk::remove_bias_all()

void ComputeTempChunk::restore_bias(int i, double *v)
{
- int index = cchunk->ichunk[i];
+ int index = cchunk->ichunk[i]-1;
   if (index < 0) return;
   v[0] += vcmall[index][0];
   v[1] += vcmall[index][1];
@@ -749,11 +753,11 @@ void ComputeTempChunk::restore_bias_all()

   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
- index = ichunk[i];
+ index = ichunk[i]-1;
       if (index < 0) continue;
- v[i][0] += vbias[0];
- v[i][1] += vbias[1];
- v[i][2] += vbias[2];
+ v[i][0] += vcmall[index][0];
+ v[i][1] += vcmall[index][1];
+ v[i][2] += vcmall[index][2];
     }
}

axel.

Looks good, thank you Axel.

-David

thanks, should be in a patch for LAMMPS soon.