Dear mailinglist,
Currently I'm implementing a new compute style, which should calculate
some averages in a local environment of an atom. More concrete: within a
sphere of radius force->pair->cutforce averaged values will be
calculated. One of them is a local stress:
(virial[i][X] + virial[i][Y] + virial[i][Z])
The sekeleton (run through neighbors etc) of the new compute is ripped of
the compute_centro_atom.cpp. The virial is accessed from the compute
compute_stress_atom.cpp.
The problem which occurs is the following: The serial version of the
code works smoothly, but the || code leads to artificial structures (one
can clearly see the edges of the decomposition). The stress for the
ghost atoms is too high. An example:
stress_avg_ghostsinneighborhood=88152.5
stress_avg_noghostsinneighborhood=60280 (real units)
From within the same new compute style local temperature or density can
be calculated without problems. Also if I visualize the pressure/atom
from compute_stress_atom.cpp itself, there are no artifacts. There seems
to be a problem in my code with accessing properties of ghost atoms
calculated by another compute style.
My question to the developers is now: Is there anything, which I have to
take into account, when I want to access the pressure tensor (or virial)
from a new compute style, which does a summation over neighbor lists? Do
the ghost atoms contain the current values, after I call
c_virial->compute_peratom()?
If there are more detailled questions, here are the relevant(?)
snipplets:
/************************/
void ComputeAvgAtom::init()
==> request neighbor lists as in compute_centro_atom.cpp
/************************/
double ComputeAvgAtom::ComputeStress (int i)
...
if (force->pair)
stress += force->pair->vatom[i][X] + force->pair->vatom[i][Y] +
force->pair->vatom[i][Z];
==> etc for all other force types
return stress
/************************/
void ComputeAvgAtom::compute_peratom ()
==> neighbor list stuff etc
...
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
==> init some vars
stress = ComputeStress (i);
...
if (mask[i] & groupbit) {
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (j >= nall) j %= nall;
rsq = SQ(x[i][X]-x[j][X]) + SQ(x[i][Y]-x[j][Y]) + SQ(x[i][Z]-x[j][Z]);
if (rsq < rcut_pot_sq) {
...
stress += ComputeStress (j);
...
vector_per_atom[i][0] = stress;
...
If there are even more question, I could pm a snapshot of my repository.
Best regards:
Gerolf