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