So I have this weird error which I think is a result of my not understanding the neighbor list correctly. I am computing gradients atom by atom using our new force field. There are two parts to our gradients a bond part which applied a force along the line of the bond (like a morse) and a vector part which applies a force perpendicular to the bond. I then apply this each of the two atoms involved in the bond. For example in the data file i have water molecule, with an oxygen at position 1 and hydrogens in positions 2 and 3. There are only two bonds in the current calculation between O and the Hs. (We could include a HH term but for now we are not).
So I calculate the forces and then apply them as follows:
f[i][0:2]+=[fx,fy,fz] (This is actually three lines of code, but the ideas is the same).
and
f[j][0:2]-=[fx,fy,fz]
This works fine for i=1 and j-2 & 3. However when I get to i=2 (the first hydrogen), and am computing its gradient force on the oxygen (not necessarily the same as before in fact it will be a little different). The force on j (the oxygen) has been reset to zero. It has somehow lost the first two gradient contributions.
How should I correct this?
matthew
So I have this weird error which I think is a result of my not
understanding the neighbor list correctly. I am computing gradients atom by
atom using our new force field. There are two parts to our gradients a bond
part which applied a force along the line of the bond (like a morse) and a
vector part which applies a force perpendicular to the bond. I then apply
this each of the two atoms involved in the bond. For example in the data
file i have water molecule, with an oxygen at position 1 and hydrogens in
positions 2 and 3. There are only two bonds in the current calculation
between O and the Hs. (We could include a HH term but for now we are not).
So I calculate the forces and then apply them as follows:
f[i][0:2]+=[fx,fy,fz] (This is actually three lines of code, but the ideas
is the same).
and
f[j][0:2]-=[fx,fy,fz]
This works fine for i=1 and j-2 & 3. However when I get to i=2 (the first
hydrogen), and am computing its gradient force on the oxygen (not
necessarily the same as before in fact it will be a little different). The
force on j (the oxygen) has been reset to zero. It has somehow lost the
first two gradient contributions.
How should I correct this?
there is not enough information to give a specific answer. most likely
your code is making incorrect assumptions.
for examples, atoms don't have to have a specific order in the way you
address them. also, with periodic boundary conditions, the same atom
(i.e. the same atom ID) may be present multiple times on the same
processors or distributed across multiple processors and those forces
will by accumulated through a reverse communication.
you can identify those atoms reliably only via their atom IDs (i.e.
via atom->tag). please have a look at the TIP4P styles that need to do
this in order to compute the location of the M pseudo particle.
axel.
Axel,
Right but if I delete the J block and only do the following three lines:f[i][0:2]+=2[fx,fy,fz]
will that give the same result after the center of mass shift?
or do I need to do a translation from the jlist to the ilist?
or is the correct approach something I am not considering? (I’ll take a look at the Tip4p.)
thanks
matthew
Axel,
Right but if I delete the J block and only do the following three lines:
f[i][0:2]+=2[fx,fy,fz]
will that give the same result after the center of mass shift?
or do I need to do a translation from the jlist to the ilist?
or is the correct approach something I am not considering? (I’ll take a look at the Tip4p.)
I don’t know your code and I don’t know your model and thus I cannot say.
Axel
So I think I can put my question in a way that doesn’t refer to my code:
in the pair_lj_cut and the could_cut styles. There is the following code:
f[i][0] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delzfpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delzfpair;
}
What I don’t understand and I think it is something to do with the newton_pair treatment is if newton_pair is off why isn’t it:
f[i][0] += 2delxfpair;
f[i][1] += 2delyfpair;
f[i][2] += 2delzfpair;
If you have a pair of bound atoms and the rule fij=-fji then when you remove the net force on all atoms in your system then half of fij should get subtracted from fji right?
Matthew
So I think I can put my question in a way that doesn't refer to my code:
in the pair_lj_cut and the could_cut styles. There is the following code:
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
What I don't understand and I think it is something to do with the
newton_pair treatment is if newton_pair is off why isn't it:
f[i][0] += 2*delx*fpair;
f[i][1] += 2*dely*fpair;
f[i][2] += 2*delz*fpair;
this doesn't make any sense at all. if newton_pair is off *and* atom j
is a ghost atom, the same atom will be a local atom on a different
processor and the bond will be listed there as atom i with an index <
nlocal. and the atom i will become a ghost with an index >= nlocal.
If you have a pair of bound atoms and the rule fij=-fji then when you remove
the net force on all atoms in your system then half of fij should get
subtracted from fji right?
please note, that lj/cut uses a *half* neighbor list, i.e. each fully
local pair is listed only once.
newton_pair now affects how those neighbor lists are created for pairs
that straddle sub-domain boundaries, those are pairs where one of the
two atoms is a ghost atom and the other a local atom. the index i is
always local, so only j may be a ghost.
with newton_pair on, those pairs are listed only once globally, i.e.
the force on the ghost needs to be accumulated and then collected onto
its original local atom via the reverse communication. with
newton_pair off, forces are only stored on local atoms and no
communication is needed (but each pair straddling subdomains will be
computed twice).
when you use a full neighborlist, everything depends on whether you
skip over double counted pairs and how (see e.g. the non-bonded
interactions in pair_sw.cpp or pair_tersoff.cpp). if you don't skip,
then - for pairwise additive forces - you *only* add to atom i, and
only tally half the energy, because all pairs are listed twice and
newton's 3rd law is not assumed at all.
axel.
[...]
If you have a pair of bound atoms and the rule fij=-fji then when you remove
the net force on all atoms in your system then half of fij should get
subtracted from fji right?
i don't understand this statement. the rule fij=-fji *is* newton's 3rd law.
but why would i "remove the net force on all atoms"?
what the code is doing here, is to use the symmetry fij = - fji to
avoid having to compute the distance r_ij = r_ji and the force a
second time.
axel.