Total energy difference for Wolf-Method and Ewald-Summation in a simulation of SPC/E water

Dear all,

thanks a lot for investing so much time in clarifying this problem. I am
impressed by how active this forum is.

@ Stan:

Just to be sure, eventually "special_bonds lj/coul 0 0 .5" did not fix the
problem, right? I still don't completely understand why special bonds would

right.

be important for my system. For example, when I wrote my own simulation of
SPC/E water, I used Lagrange multipliers + SHAKE even without a harmonic
potential between the bonds. I can see, that the potential does not affect
the solution, but probably it helps the SHAKE algorithm to converge faster.

more importantly, fully excluded pairs are completely removed from the
neighbor list and thus (if there are a lot of them) this will speed up
the calculation.

If I understood "special_bonds" in the context of SPC/E correctly, the first
0 tells LAMMPS to ignore O-H pairwise interaction on a given molecule, the
second 0 tells it to ignore H-H pairwise interactions on the same molecule.

However, does the third coefficient, 0.5, even have a meaning for a 3-site
model?

no. it has none.

@ Axel:

Did the error eventually only effect energies or energies an forces meaning
that previous results employing the Wolf method should not be trusted?

that is a bit of a trick question. it should not affect the forces in
your case, since you use shake to project them out. it would, however,
affect any calculations with bonded interactions that do not use shake
(or fix rigid for that matter).

axel.

Thanks for clarifying that, Axel!

Best,
Peter

Hello,

I’m writing a fix (not a pair_style) browsing all atoms with the usual “for(ii=0;ii<inum;ii++)” and all the neighbors of each “ii” atom with the usual “for(jj=0;jj<inum;jj++)”.

For example, if atom #512 has been browsed as “ii” and #513 as “jj”, i need also inverse : atom #513 browsed as “ii” and #512 as “jj” (meaning “newton off”).

As i want to check all atoms in a sphere around atom “ii” and to calculate the distance rij, i need neighbors even in neighboring domains/regions.

  1. Currently i wrote this about the init :

void FixXXX::init()

{

int irequest = neighbor->request(this);

neighhbor->requests[irequest]->pair = 0;

neighhbor->requests[irequest]->fix = 1;

neighhbor->requests[irequest]->half = 0;

neighhbor->requests[irequest]->full = 1;

}

What should i add or modify in this init?

  1. Do i need in my FixXXX::post_force(int vflag) the “j &= NEIGHMASK” i can found in some pair_style?

  2. I found also a function “setmask” in some fix_style. Do i need it to access all needed neighbors?

Example :

int FixXXX:setmask()

{

int mask = 0;

mask |= POST_FORCE;

}

Thanks for your answers.

Best regards,

Xavier Bidault

PhD student – Laboratoire de Photonique d’Angers – Université d’Angers - France

Hello,

I’m writing a fix (not a pair_style) browsing all atoms with the usual
“for(ii=0;ii<inum;ii++)” and all the neighbors of each “ii” atom with the
usual “for(jj=0;jj<inum;jj++)”.
For example, if atom #512 has been browsed as “ii” and #513 as “jj”, i need
also inverse : atom #513 browsed as “ii” and #512 as “jj” (meaning “newton
off”).
As i want to check all atoms in a sphere around atom “ii” and to calculate
the distance rij, i need neighbors even in neighboring domains/regions.

sorry, but that is not a convincing argument for requiring a full
neighbor list. the same argument applies to all pair styles with a
pairwise additive potential and there you get by just fine with
computing each pair only once followed by simply communicating the
contributions to ghost atoms back to the original local atoms. so
unless you are doing something that is more sophisticated than what
you say, you won't need a full neighborlist and can do the same with a
half neighbor list followed by a reverse communication and save a ton
of time in the process.

for the questions below, i would like to remind you of a process
called "testing", where you just try out what you want to do, see
whether it produces the result you expect or not and then report back
to us with a more specific question, when it does not work as
expected. ok?

axel.