# short range interactions in Ewald sum

Dear Lammps users,

I have a question regarding the short and long range coulombic interactions. My question is why the short interactions (interactions inside the cuttoff; V=qi x qj x erfc(alpha rij)/rij ) is always calculated as positive value and long range interactions as negative values? I have wrote a simple code to calculate the coulombic interactions according to above formula and I am getting sometimes positive and sometimes negative numbers as the ensemble of pair charges (neg/neg , pos/pos, neg/pos) move in the system.

Many thanks

Dear Lammps users,

I have a question regarding the short and long range coulombic interactions. My question is why the short interactions (interactions inside the cuttoff; V=qi x qj x erfc(alpha rij)/rij ) is always calculated as positive value and long range interactions as negative values?

How do you conclude that?

I have wrote a simple code to calculate the coulombic interactions according to above formula and I am getting sometimes positive and sometimes negative numbers as the ensemble of pair charges (neg/neg , pos/pos, neg/pos) move in the system.

Most likely you made a mistake somewhere.

Axel

Deal Axel,
Thanks for the respond. I have written a simple routine for only the short range interactions (I copied the routine below). I am running this piece of code for a liquid (each molecule consists of 3 atoms of charges 0.25,-0.5, 0.25).
In my calculations the potential become positive for pair charges of 0.25/0.25 and -0.5/-0.5; and become negative for pair charges of 0.25/-0.5. These negative and positive numbers cancel out each other and I end up with very small number which fluctuate around zero. I ran the same system with lammps and the short range is a large positive number ( ~+1800 kcal/mol for all snapshots). I am really wondering what am I doing wrong here!!

uElect_shortRrange=0.0; //stores the electrostatic potential
for (int i=0; i< nmolec; i++) //zero the force
Force [i] = {0}; //zero the force

for (int ii=0; ii< nmolec -1; ii++) //nmolec is the number of molecules
{
for (int jj=ii+1; jj< nmolec ; jj++)
{
for (int iiSite=0 ; iiSite< natom; iiSite++) //natom is the number of atoms in the molecules (all molecules have the same number of atoms)
{
indexI=iiSite+ii* natom ; //find the index of atoms on molecule ii

for (int jjSite=0; jjSite< natom ; jjSite++)
{
indexJ=jjSite+jj* natom ; //find the index of atoms on molecule jj
dr.x=rx[indexI] - rx[indexJ];
dr.y=ry[indexI] - ry[indexJ];
dr.z=rz[indexI] - rz[indexJ];
dr.x-=fast_nearest(dr.x); //find the shortest distance in the periodic box
dr.y-=fast_nearest(dr.y); //find the shortest distance in the periodic box
dr.z-=fast_nearest(dr.z); //find the shortest distance in the periodic box
rijSqr = dr.xdr.x + dr.ydr.y + dr.z*drz;
rij = sqrt (rr);

if ( rij < rCut )
{

Ewald_force_mag = charge[indexI] * charge[indexJ] * { erfc (alpha * rij) / rijSqr /rij + 2. * alpha * exp (- alpha* alpha * rijSqr ) / sqrt(pi) / rijSqr };
Ewald_force.x= force_mag * dr.x; //force vector due to short range ewald
Ewald_force.x= force_mag * dr.y;
Ewald_force.x= force_mag * dr.z;

Force[ indexI ].x += Ewald_force.x;
Force[ indexI ].y += Ewald_force.x;
Force[ indexI ].z += Ewald_force.x;

Force[ indexJ ].x -= Ewald_force.x;
Force[ indexJ ].y -= Ewald_force.x;
Force[ indexJ ].z -= Ewald_force.x;

uElect_shortRrange += charge[indexI] * charge[indexJ] * erfc (alpha* rij)/ rij;

}
}
}
}
}

________________________________
Sorry, there was a little typo in the parameters
Arshia

From: Arshia Movahed <[email protected]...>
To: "[email protected]"
<[email protected]>
Sent: Monday, February 20, 2012 1:42 PM

Subject: Re: [lammps-users] short range interactions in Ewald sum

Deal Axel,
Thanks for the respond. I have written a simple routine for only the **short
range** interactions (I copied the routine below). I am running this piece
of code for a liquid (each molecule consists of 3 atoms of charges
0.25,-0.5, 0.25).
In my calculations the potential become positive for pair charges of
0.25/0.25 and -0.5/-0.5; and become negative for pair charges of 0.25/-0.5.
These negative and positive numbers cancel out each other and I end up with
very small number which fluctuate around zero. I ran the same system with
lammps and the short range is a large positive number ( ~+1800 kcal/mol for
all snapshots). I am really wondering what am I doing wrong here!!

where is the lammps input?

difficult to tell why there is a difference,
if i cannot see the difference.

axel.

Dear Axel,
I have attached my data file + input file. I am only comparing the short range coulombic part (the vdw part are not important for me).
Many Thanks
Arshia

in.class2 (1.59 KB)

test_dme.lammps05 (31.2 KB)

ok. thanks.

now that i see your lammps input, the explanation is very simple:

your test code doesn't seem to handle exclusions,
i.e. skip over bonded interactions when computing
the non-bonded contributions.

if you change the lammps input by replacing:

special_bonds lj/coul 0.0 0.0 1.0

with

special_bonds lj/coul 1.0 1.0 1.0

you'll get a negative real space coulomb energy.

of course, this will be a bad model for a real
simulation, but that is beside the point here.

axel.

Dear Axel,
Thank you very much for your respond.

But the way that for-loops have been set, they would not count any interaction on the same molecule. And since my molecules have only 3 sites, the bond and angle terms (bonded) are in fact excluded (and I do not have any torsion of courseâ€¦)

And yes you are right about changing the special bond in lammps; I did check that.

I have become very confused that why I can not get the same results are lammps!!! It seems that I am making a big mistake since the numbers I get are very very different (1800 vs ~ 0 kcal/mol!).

Many thanks
Aeshia

And another thing is that, I get neg and positive values (fluctuating around zero since numbers cancel out each other)!
Thanks again

And another thing is that, I get neg *and* positive values (fluctuating
around zero since numbers cancel out each other)!

i don't have the time to debug your code,
but my suggestion is to first disregard the
short range screening and compare against

lj/class2/coul/cut

with a suitable cutoff and make sure that
that part is working properly and then
check what happens when you re-introduce
the screening.

axel.

Dear Axel,
Are the Windows source code (to run on series) available (I know that the executable file is available on the website)? Something that one can easily debug?

(debugging the linux version is difficult for me since I do have the root privilege for our cluster and hence building the code in debug mode would make lots of difficultiesâ€¦)

Thanks

Dear Axel,
Are the Windows source code (to run on series) available (I know that the
executable file is available on the website)? Something that one can easily
debug?

the sources for windows are the same as for unix.

(debugging the linux version is difficult for me since I do have the
root privilege for our cluster and hence building the code in debug mode
would make lots of difficulties...)

this is *complete* nonsense.
you don't need root to compile in debug mode
(or compile *anything*).

my experience with compiling on different platforms
is rather that compiling on windows is a major PITA.

axel.