confuse about the code

[Lisa_Lan_5a]

Hi,

I am confused about the different meaning of these two methods for force calculations. I thought they will get the same results, but the test show that when I use the first method, the system will explode under NPT ensemble during relax and the second way is good.

for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
…
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
…
The first mehod is :
twobody(&params[ijparam],rsq,gfinal,fpair,eflag,evdwl,f2w);
f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delz*fpair;

if (evflag) ev_tally_full(i,evdwl,0.0,fpair,delx,dely,delz);
}
}

And the other way is :

twobody(&params[ijparam],rsq,gfinal,fpair,eflag,evdwl,f2w);
fpair=fpair/2.0;
evdwl=evdwl/2.0;

f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delzfpair;
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delzfpair;
if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);

Thanks.

Lisa

[akohlmey]

lisa,

Hi,

I am confused about the different meaning of these two methods for force
calculations. I thought they will get the same results, but the test show
that when I use the first method, the system will explode under NPT ensemble
during relax and the second way is good.

did you also ask for a "full" neighbor list in the first case?
if not, then the "explosion" is completely understandable.

in the first case, you would need the neighbor list to contain
_all_ pairs of particles that are likely to be within the cutoff.

in the second case, you only need _unique_ pairs, since
you are applying newton's 3rd law.

please note, that this is different from the "newton" keyword
that handles the case of whether to apply newton's 3rd law
for "ghost" particles, i.e. atoms that do not "belong" to the
specific mpi task.

axel.

p.s.: if you do need to work with a "full" neighbor list,
you will have to enforce "newton off" for pair interactions.

[Lisa_Lan_5a]

Axel,

For the both cases, I use the same setup which is the same with the the stilling-weber potential. The only difference is that I do not skip half of the j atom in two body calculations.
void PairSW::init_style()
{
if (atom->tag_enable == 0)
error->all(“Pair style Stillinger-Weber requires atom IDs”);
if (force->newton_pair == 0)
error->all(“Pair style Stillinger-Weber requires newton pair on”);
// need a full neighbor list
int irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
}

So there should not have a explosion in the first case.

2011/7/20 Axel Kohlmeyer <[email protected]>

[akohlmey]

Axel,

For the both cases, I use the same setup which is the same with the the
stilling-weber potential. The only difference is that I do not skip half of
the j atom in two body calculations.
void PairSW::init_style()
{
if (atom->tag_enable == 0)
error->all("Pair style Stillinger-Weber requires atom IDs");
if (force->newton_pair == 0)
error->all("Pair style Stillinger-Weber requires newton pair on");
// need a full neighbor list
int irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
}

So there should not have a explosion in the first case.

well,

this is really hard to tell from just the few pieces
of information that you are providing.
there are _so_ many ways in how you this can be
messed up, particularly since you don't seems to be
very familiar with the lammps internals and also the
intricacies of implementing non-bonded interactions
for an MD code.

if you really want to get down to it, would recommend
to implement a version of lj/cut that uses the full neighborlist
and thus no newton's 3rd law. this is such an easy case, that
it is easy to track down problems in the physics and errors
in the implementation. it is rather difficult to do this with a
newly implemented code, where you don't really know how
bug-free or bug-ridden it is.

cheers,
    axel.

[Lisa_Lan_5a]

Axel,
I have worked out the pair_lj_cut.cpp to be the full neighbor list. It seems to be the same as what I thought to be. I remove the
if (newton_pair || j < nlocal) {
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delz*fpair;
}
and change ev_tally function to be ev_tally_full. Finally, I add the
irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
in the init_style function.

Though it make the same result with the half neighbor list, I found that there are so many things I don’t understand,
the eflag and vflag are the flag of compute per atom value or the global value, right? But I don’t see any change of the eflag and vflag. How can I know if it is per atom or global?

Why does the pair_lj_cut need to has the single function and calculate the force and energy, while the tersoff and sw don’t need to?

Thank you.

Lisa
2011/7/20 Axel Kohlmeyer <[email protected]>

[akohlmey]

Axel,
I have worked out the pair_lj_cut.cpp to be the full neighbor list. It seems
to be the same as what I thought to be. I remove the
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
and change ev_tally function to be ev_tally_full. Finally, I add the
irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
in the init_style function.

right.

Though it make the same result with the half neighbor list, I found that
there are so many things I don't understand,

the eflag and vflag are the flag of compute per atom value or the global

eflag is >0 when energy contributions need to be computed.
not computing them, when they are not needed, makes lammps faster.

same goes for the vflag, which is for virial contributions.
this can be done in different ways, depending on whether
per atom virial contributions are needed or not.

value, right? But I don't see any change of the eflag and vflag. How can I
know if it is per atom or global?

eflag and vflag are bitmasks and the exact task is encoded.
for the Pair::compute() method it only matters if they are
zero (no work to be done) or not (work to be done).

Why does the pair_lj_cut need to has the single function and calculate the
force and energy, while the tersoff and sw don't need to?

the Pair::single() is not required, but needed for calculations
like compute group/group
however, this can only be done for pairwise additive potentials.
for many-body potentials there is no pairwise decomposition
of the interactions (or else they would not be manybody potentials).

cheers,
    axel.

[Lisa_Lan_5a]

Ok, thank you for your reply. I have done a few tests and found that if there are only two atoms in the sample, and I run a NVE ensemble ,the results of the two methods are exactly the same. But when I increase the atoms to be 1350. the two methods get the same Pe and Ke, but the different Pressure.
Step Temp Press PotEng KinEng TotEng Volume
I. 0 0 -123960.15 -7017.6677 0 -7017.6677 16710.449

II. 0 0 -211168.39 -7017.6677 0 -7017.6677 16710.44

So what cause the difference and how can I decide which one is right?

The first method:
twobody(&params[ijparam],rsq,gfinal,fpair,eflag,evdwl,f2w);
f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delz*fpair;

if (evflag) ev_tally_full(i,evdwl,0.0,fpair,delx,dely,delz);

And the second one is :

twobody(&params[ijparam],rsq,gfinal,fpair,eflag,evdwl,f2w);
fpair=fpair/2.0;
evdwl=evdwl/2.0;
f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delzfpair;
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delzfpair;
if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0,fpair,delx,dely,delz);

2011/7/21 Axel Kohlmeyer <[email protected]>