I am trying to write a compute temp/omega measuring the temperature of
a group after the center of mass velocity and the angular velocity
have been removed. The modified code seems to be working both in
serial and parallel when I only try to compute the temperature, and it
also works in serial when I try to bias my thermostat with compute
temp/omega. But in parallel the code hangs just after thermo output
for timestep 0 has been printed on the screen. I'm not at all an
expert in programming and here I don't have a clue on how to debug
this... Could someone help me on this?
Enclosed are the modified files with regard to current lammps
distribution, and the input file generating the problem...
I finally found out that my issue came from using group functions (for
instance group->vcm) inside remove_bias and restore_bias, where the
MPI_Allreduce command hangs. This could be solved by calling
remove_bias_all and restore_bias_all in the thermostat.
However I noticed that remove_bias_all and restore_bias_all, though
present in each compute_*.cpp code, is not called by any thermostat in
the current lammps version. Only remove_bias and restore_bias are
used.
So I had two questions:
* why did you chose at some point to only call remove_bias and
restore_bias, i.e. remove and restore bias atom by atom, instead of
calling remove_bias_all and restore_bias_all?
* If I slightly modify fix_nh.cpp to call remove_bias_all and
restore_bias_all (see below), would it work with all the other
computes (since the corresponding code exists in all others
compute_*.cpp)? If I understood correctly I believe it should provide
the exact same results, but I just wanted to check...
The bias_all() methods are called by other computes,
like compute temp/sphere. The thermostats use
the per-particle bias methods, so I wouldn't change
that. The idea is that the compute, say compute temp/rotate
in your case, is supposed to do any global operations it
needs to calculate the bias, when it's compute_scalar()
method is invoked to calculate the temperature. E.g. the
group operations you mention. This can involve creating
some data structure with the bias values if needed. Then
when remove/restore bias is called on a per-particle basis
it can do what is needed. All the thermostats call compute_scalar()
before they call remove_bias(), so this works. See
compute temp/com for an example.
Is there some reason that model doesn't work for what you
want to do?
Thanks a lot for this clarification! Doing as you wrote seems to work
like a charm. Just to check: in order to remove the angular velocity,
I have to store a velocity bias for all atoms, so I define: "double
**vbiasall;" and "int maxbias;" in "compute_temp_rotate.h".
Then I just have to add to ComputeTempRotate():
maxbias = 0;
and before I affect vbiasall[i][0,1,2] in compute_scalar() AND
compute_vector() (or only in compute_scalar?):